Classical and quantum $f(R)$ cosmology: The big rip, the little rip and the little sibling of the big rip

The big rip, the little rip and the little sibling of the big rip are cosmological doomsdays predicted by some phantom dark energy models that could describe the future evolution of our Universe. When the universe evolves towards either of these future cosmic events all bounded structures and, ultimately, space-time itself are ripped apart. Nevertheless, it is commonly belief that quantum gravity effects may smooth or even avoid these classically predicted singularities. In this review, we discuss the classical and quantum occurrence of these riplike events in the scheme of metric $f(R)$ theories of gravity. The quantum analysis is performed in the framework of $f(R)$ quantum geometrodynamics. In this context, we analyse the fulfilment of the DeWitt criterion for the avoidance of these singular fates.


Introduction
The discovery that our universe is currently undergoing an accelerated expansion phase has supposed an inflexion point in our understanding of the physics of the cosmos. Even though this acceleration was first noticed more than 20 years ago [1,2], its mysteries are still not yet fully unraveled. On the contrary, understanding the mechanism involved in this phase represents one of the greatest milestones in modern cosmology. This behavior is so challenging because no form of matter we know from our ordinary experience can actually produce this phenomenon. In general relativity (GR), this phase is attributed to the existence of an exotic form of energy with a negative pressure that causes the repulsion of the matter content and, thus, pushes the universe into expansion. This exotic content is dubbed "dark energy" (DE). The simplest DE model, and still the one which best fits the observational data, is the standard ΛCDM theory, where DE plays the role of a cosmological constant. Nevertheless, this model contains various theoretical and observational unresolved puzzles. To mention some of them: the nature of dark matter (DM) and DE [3][4][5][6][7], the coincidence problem [8], the cosmological constant problem [9] and the new tensions in certain cosmological parameters. (For a review of the topic see, for instances, references [10,11] and references therein.) Therefore, one may expect the ΛCMD paradigm to be just a useful effective model. Hence, a great number of alternative models for DE describing the acceleration phase have been proposed. Some models are phantom DE [12,13], tachyonic matter [14,15], Chaplygin gas [16,17], holographic among other singularities such as the big bang, can be avoided as a result of quantum effects emerging as the universe approaches the classical singularity [66,72,73]. However, since the same classical background evolution 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. In this work, we gather and review the so far published results about the classical and quantum occurrence of these three riplike events when the description of gravity is that provided by f (R) metric theories of gravity. For that purpose, we shall focus on the quantum cosmology scheme given by the Wheeler-DeWitt (WDW) equation [74] being adapted for the case of f (R) gravity [75], namely f (R) quantum geometrodynamics. Consequently, we explore here the possibility of avoiding the BR singularity, and the LR and the LSBR abrupt events in f (R) quantum cosmology.
This manuscript is organized as follows. In Section 2, we review some phantom DE models predicting the classical appearance of future singularities and/or abrupt events. We also provide a new classification of the cosmic singularities and abrupt events found there. Moreover, we discuss the viability of those models from the point of view of cosmological observations. In Section 3, we consider that the classical background evolution found in the previous section can be equivalently described in the framework of metric f (R) theories of gravity, where in the latter case the expansion of the universe will have a purely geometrical origin. For that aim, we briefly discuss a reconstruction method for metric f (R) theories of gravity. Thereafter, in Sections 3.1-3.3, we apply this background reconstruction technique to the specific phantom DE models reviewed in Section 2.1. Thus, we present the group of metric f (R) theories of gravity predicting a classical fate à la BR, LR or LSBR. Section 4 is entirely devoted to the revision of the fate of classical singularities in f (R) quantum cosmology. Therefore, we summarize the f (R) quantum geometrodynamics approach in Section 4.1. Then, the resolution of the modified Wheeler-DeWitt equation corresponding to the BR, LR and LSBR doomsdays is addressed in Sections 4.2-4.4, respectively. We also discuss there the avoidance of these cosmological catastrophes by means of the DeWitt criterion. Finally, we check the validity of the approximations performed to solve the modified Wheeler-DeWitt equations in Appendices A and B. Table 1. Classification of the riplike cosmic events by means of the time of occurrence of the rip t rip , the scale factor a, the Hubble parameter H and its cosmic time derivativeḢ. Please note that the pseudorip has been included for the sake of completeness, as this corresponds to a mild event (before which all bounded structures are disintegrated) rather than to a curvature singularity.

A Phantom Dark-Energy Universe
Throughout this review we consider homogeneous and isotropic cosmological scenarios, which are described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric given by where t stands for the cosmic time, a represents the scale factor, ds 2 3 denotes the threedimensional spatial metric and we have used the geometric unit system 8πG = c = 1.
Assuming the content of the universe to be described by a perfect fluid with energy density ρ and pressure p, the Friedmann and Raychaudhuri equations are respectively, where the dot denotes the derivative with respect to the cosmic time t, H stands for the Hubble rate and k represents the spatial curvature (not fixed at this point).
Within the standard interpretation of the observational data, the cosmic fluid at present is constituted by three species: radiation, matter (baryonic and dark) and dark energy. The present concentrations of matter and dark energy are roughly Ω M,0 = 0.315 and Ω DE,0 = 0.685, respectively, whereas the contribution of radiation to the total content of the universe is negligible at the present time [37,38]. Therefore, DE is the dominant ingredient today. Moreover, because we aim to study the asymptotic future expansion of the universe, we assume that the DE density increases, remains constant or decreases more slowly than the non-relativistic matter energy density and the spatial curvature term k/a 2 . Thus, from a practical point of view, we can neglect the contribution of the other cosmic ingredients when studying the future asymptotic behavior of these models. Consequently, hereafter p and ρ denote the pressure and the energy density of the dark fluid. Then, the DE density evolves aṡ To evaluate the future fate of a DE dominated universe in the perfect fluid representation, an equation of state (EoS) for the dark fluid must be provided. We emphasize that presently observational data are compatible with the w = −1 value, though there is a significantly overlap with both possibilities w < −1 and w > −1. On the other hand, from a theoretical point of view, the value −1 in the w-axis represents a qualitative change in the properties of DE, as well as in the future evolution of the universe. Therefore, the cosmological constant case (w = −1) seems to play a special role from both theoretical and observational perspectives. Then, let us begin with a brief review of the phenomenology of the following EoS for the DE content, which is reminiscent of an expansion around a cosmological constant, being A a positive constant 1 . This equation of state was first introduced in references [64,68] and thoroughly discussed in terms of singularity occurrence in reference [52]. In this section, of the manuscript, we review the analysis performed in reference [52] and present a novel alternative classification of the singularities there found by means of the behavior of the scale factor a, the Hubble rate H and its cosmic time derivativeḢ. That is a metric classification such as the one presented in [65]. (See also the characterization given in [64] and the study on cosmic singularities presented in references [76,77].) We argue for the need for such a classification, instead of only addressing the evolution of the scale factor a and the DE density ρ, since different cosmic events such as the LR and LSBR cannot be differentiated in the latter picture.
From the continuity equation (4) it follows that the DE density evolves as for α = 1, and for the case of α = 1, where the subscript "0" denotes the current value of the corresponding quantity. In either case, the EoS parameter w reads Note that for a non-negative parameter A, the denominator in the r.h.s. in equation (8) is always positive 2 , whatever the value of α. Thus, the DE content modelled by the EoS (5) exhibits a phantom-like behavior when the parameter A is positive. Please note that this can be also seen directly from equation (4) for an expanding (H > 0) universe. On the other hand, note that α = 1 and α = 1 2 are special values on the α-line that delimited qualitatively alike cosmological behaviors [52]. This can be deduced from the time dependence of the scale factor obtained from equations (2) and (6). This is [52] ln for the case of α = 1 2 , and [52] ln a a 0 = 1 for α = 1 2 and α = 1, and for α = 1. B is a constant defined as where we have denoted by t some arbitrary (future) moment in the expansion history of the universe from which we can safely neglect the contribution of matter fields and, therefore, assume that DE is the only content of the cosmos. In addition, a represents the corresponding scale factor. Therefore, when α = 1 2 the scale factor asymptotically approaches a double exponential growth on the cosmic time. However, for α > 1 2 , it evolves as some function of t s − t, being t s the time of occurrence of the corresponding singularity. Moreover, the expression for t s depends on the value of α. Hence, for the sake of the discussion, we proceed to summarize the main results found in [52] for each case. Those are: α > 1, α = 1, 1 2 < α < 1, α = 1 2 and α < 1 2 . Additionally, to the conclusion presented in [52], we also compute the corresponding H andḢ variables in each case. This allows us to provide a complementary classification of the cosmic events found in [52] similar to that introduced in reference [65], see Table 2.
For the case of α > 1, the universe reaches a maximum size given by a max := a 0 exp 1 Furthermore, this size is reached in a finite time into the future At that moment, the Hubble rate and its comics time derivative diverge, since This asymptotic behavior corresponds to the occurrence of a BF singularity [58,59](see also reference [76]). That is a type III singularity in the notation of [64]. The case of α = 1 corresponds to a constant EoS parameter w = −1 − A. Accordingly, the size of the observable universe becomes infinite at a finite time from present epoch, namely t rip . This is see equation (11). Furthermore, given that the Hubble rate and its cosmic time derivative also diverge at t = t rip . Therefore, for this value of α, the universe evolves towards a classical BR singularity, alike to that first introduced in [12,13,50]. This corresponds to a type I singularity according to the notation in reference [64].
On the other hand, for 1 2 < α < 1, the scale factor diverge at finite cosmic time, see equation (10) where the ratio 2(α − 1)/(2α − 1) is now negative. This makes the ln a proportional to some power of 1/(t rip − t). The time at which the observable universe becomes infinite is The Hubble rate and its time derivative follows the same relations given in equations (15) and (16), respectively. Hence, these quantities also diverge along with the scale factor. Therefore, in a finite time into the future, the scale factor, the Hubble rate andḢ explode, whereas the EoS parameter w converges to −1 from below. Of course, this implies that the DE density and pressure blow up, as was found in [52]. This event is qualitatively equivalent to a BR singularity (see the singularities classified as type I in, for example, references [63][64][65]). Please note that this behavior was also found in reference [76], where it was dubbed "grand rip" (see type −1 singularities in reference [77]). Another possible choice for α corresponds to the interesting case of α = 1 2 . In this scenario, the scale factor asymptotically evolves as an exponential function of an exponential function on the cosmic time, i.e., a ≈ e e t , see equation (9). Accordingly, the Hubble rate and its cosmic time derivative are Thus, the scale factor, H, andḢ diverge at the infinite asymptotic future. This drives the universe towards a LR abrupt event, see classification in [63,65]. (See also reference [53].) In fact, the EoS (5) with α = 1 2 corresponds to the DE model for which the name "little rip" was first given [51], even though this cosmological behavior was already known from before [52] (see also [53] where the LR was found in brane cosmology and before that in some modified theories of gravity [78]).
Finally, for the case of α < 1 2 , the scale factor obeys the relation given in equation (10). However, since now the ratio 2(α − 1)/(2α − 1) is positive, then the equation for ln a reduces to a certain polynomial on the cosmic time. This makes the expansion of the observable universe to last indefinitely. Hence, no finite time singularities are present in this case. The Hubble rate and its cosmic time derivative read For 0 < α < 1 2 , both quantities tend to infinity at the infinite asymptotic future, thus, leading to the occurrence of a LR abrupt event. On the contrary,Ḣ remains constant for α = 0, or shrinks to zero when α < 0. This behavior corresponds to a final fate à la LSBR, see reference [55]. Note that since the DE density grows unbounded in both abrupt events, the LR and the LSBR, see equation (6), this distinction could not be done analyzing only the asymptotic behavior of a and ρ. Additionally note that given some entire number n, all the higher order derivatives of H up to the n-th order diverge when α = (n − 1)/2n. Therefore, we conclude that the phantom DE model described by the EoS (5) with a positive parameter A entails a great variety of cosmological singularities and abrupt events. We summarized the results from reference [52] and our new findings here in Table 2. We will present a more detailed analysis in another work. Table 2. Classification of the singularities and abrupt events found for the phantom DE model given by the general EoS (5), see reference [52], accordingly to the value of the parameter α, the time of occurrence of the singular behavior t s , the Hubble rate H and its cosmic time derivativeḢ. Please note that since the BR and grand rip [76] singularities have qualitatively the same behavior in terms of those quantities, we do not address the possible differences between both events in the following classification. Hence, we keep the term "big rip" for both of them. α t s a HḢ Event

Cosmological Constraints
Since this review is mainly devoted to the analysis of the classical and quantum fate of the BR, the LR and the LSBR events in metric f (R) theories of gravity, we should now restrict our attention to some particular phantom DE models on which to apply the reconstruction techniques in the next sections. Therefore, hereafter we shall consider only the following phantom DE models when addressing the occurrence of these cosmic events: • For the BR singularity we consider the phantom DE model with a constant EoS parameter w < −1 [50]. This model corresponds to the choice of α = 1 in the more general EoS (5) studied in the previous section. • For a DE model with a future LR abrupt event, we select the EoS for DE described in reference [51], which corresponds to the case α = 1 2 in (5). • For a universe doomed to evolve towards a LSBR abrupt event, we consider the DE content to be described by the EoS (5) with α = 0. This LSBR was first introduced in reference [55].
In the following, for the sake of concreteness, we shall consider the cosmological constraints obtained in reference [57] when binding those DE models with observational data. Thus, the results presented in the incoming sections are subjected to the observational constraints on the parameter A presented in reference [57]. These constraints are summarized in Table 3. Please note that the small values for A there obtained suggest that tiny deviations from the ΛCDM scenario are, in fact, the observationally preferred situation today [57]. Nevertheless, we recall that the asymptotic evolution of these DE models is not that of a de Sitter universe, since the corresponding DE pressure and energy density do not converge to a constant value but diverge. In the next section, we will obtain the group of metric f (R) theories of gravity able to reproduce the same asymptotic expansion history, which in GR corresponds to these particular phantom DE models. Table 3. Best fit found in reference [57] for the phantom DE models discussed in Section 2.1, where A is dimensionless for the case of the BR, and has units of inverse of meter and inverse of square meter for the LR and LSBR models, respectively.

Phantom Dark-Energy Models in f (R) Cosmology
For a given cosmological background evolution it is possible to find a family of alternative theories of gravity that leads to the same expansion history. The group of techniques used to perform such a background "reconstruction" task are dubbed as "reconstruction methods" (for a review of the topic see, for instance, reference [80] and references therein). In this part of the review, we shall focus our attention on reconstruction methods within the scheme of metric f (R) theories of gravity. Hence, we look for a metric f (R) theory of gravity able to reproduce the same superaccelerated expansion profile to that of the general relativistic model filled with phantom DE with the EoS (5), for the values of α selected in Section 2.1. It is worth noting that whereas in GR the accelerated phase is attributed to the existence of an exotic form of energy with a negative pressure (DE), in the setup of f (R) theories of gravity the same expansion has a rather geometrical origin. Furthermore, as the general relativistic model with the EoS (5) expands the cosmos towards some future doomsdays, then the resulting metric f (R) theory of gravity will lead to the same classical fate. For previous works on reconstruction techniques in f (R) gravity see, for instance, references [80][81][82][83][84][85][86]. See also references [87][88][89][90][91] for successful reconstruction of phantom DE-driven riplike events in metric f (R) theories of gravity.
Henceforth, we shall refer to two cosmological evolutions as equivalent at the background level if the corresponding geometrical variables H,Ḣ, R andṘ are identical [87]. In GR, the expansion of a homogeneous and isotropic universe is ruled by the Friedmann and Raychaudhuri equations (2) and (3), respectively. Accordingly, the scalar curvature reads where ρ and p denotes the total energy density and pressure, respectively, these are ρ = ρ DE + ρ M and p = p DE + p M . Later we will assume that they just correspond to the DE component when the cosmic matter is diluted. Additionally, from the continuity equation for the perfect fluid (4) and the Friedmann equation (2), it followṡ where we have assumed p = p(ρ). Consequently, the cosmic time derivative of the scalar curvature R readsṘ On the other hand, for the gravitational interaction being that provided by metric f (R) theories of gravity, the evolution of the universe is described by the action where S m stands for the minimally coupled matter fields. In this framework, the field equations no longer coincide with (2) and (3). In fact, the modified Friedmann equation reads being ρ m the energy density of the minimally coupled matter fields. Since we are interested in a metric f (R) theory able to reproduce the same background cosmological expansion as a certain general relativistic model, then the preceding expression can be considered to be a differential equation for some, a priori unknown, function f (R), where the coefficients are already fixed, i.e., when the geometrical quantities involved in equation (30) are set to be equal to those of the GR model that we want to reproduce. Thus, the background cosmological expansion of the resulting metric f (R) theory of gravity will be equivalent to that provided by the general relativistic model. Furthermore, as we are interested on the asymptotic future behavior of the universe, we can neglect the matter and spatial curvature contribution in equations (25) and (30), which will be (more) quickly redshifted with the superaccelerated expansion. Therefore, we drop again the subindex DE for the energy density and pressure. Hence, in the following sections we present the most general solution to this reconstruction procedure when considering the EoS (5) for the different values of α discussed in Section 2.1. Nevertheless, the obtained f (R) theories must be understood as useful asymptotic models for the study of the behavior of the universe near classical singularities, since a more realistic background reconstruction would imply the non-cancellation of the matter and spatial curvature contribution in equations (25) and (30).

BR Singularity in f (R) Gravity
For the BR singularity we considered the case of α = 1 in the DE general EoS (5). This is For that specific value of α, the cosmic time derivative of the Hubble rate can be re-expressed aṡ compare equations (18) and (19). Consequently, the scalar curvature and its cosmic time derivative reduces to Therefore, the modified Friedmann equation (30) reads where we have used the notation f R := d f /dR and f RR := d 2 f /dR 2 . The most general solution for the function f (R) was already obtained in reference [87]. That is being c + and c − arbitrary integration constants and In general, the parameter γ ± may take complex values. However, for A ≤ (10 − 4 √ 6)/3 ∪ A ≥ (10 + 4 √ 6)/3 , both branches are real valued. Please note that this is precisely the case for the values of A showed in Table 3. Hence, when the observational constraints [57] on the model are taken into account, both branches of the parameter γ ± are real valued.

The LR in f (R) Gravity
For the LR abrupt event, we consider here the EoS (5) for DE with α = 1 2 in (5). This is Note that for this choice of α, the Hubble rate in equation (21) is an exponential function on the cosmic time. Subsequently, its cosmic time derivative is proportional to itself. Thus, for the sake of simplicity, we denoted by β that proportionality constant. That iṡ where we have defined The curvature scalar and its cosmic time derivative in terms of the Hubble rate are As noted in reference [88] (se also [91]), the modified Friedmann equation (30) simplifies when rewritten in terms of H. Thus, we need to solve where f H := d f /dH, f HH := d 2 f /dH 2 , and ρ m and k has been neglected, as discussed before. The most general solution to the above differential equation is [88] being C 1 and C 2 integration constants and Ei the exponential integral function (see definition, e.g., in 5.1.2 of reference [92]). To obtain the final f (R) expression, the relation (41) must be inverted. This is At a first sight, it seems that this transformation may be non-univocally defined. However, since the phantom DE density (6) increases with the expansion of the universe towards the LR abrupt event, whereas the matter content is quickly diluted, then the corresponding Hubble rate (21) and scalar curvature (41) also grow. Hence, only the positive branch in the preceding expression applies. Therefore, represents the most general family of metric f (R) theories of gravity predicting the occurrence of that specific LR abrupt event [88] (see also [91]). We recall that β, which depends on the parameter A, is observationally constrained in Table 3.

The LSBR in f (R) Gravity
The selected DE model in this case corresponds to α = 0 in the EoS (5). That is where A takes the value showed in Table 3. In this scenario, see equation (24), the cosmic time derivative of the Hubble rate remains constant. Its value is given bẏ Hence, the scalar curvature and its cosmic time derivative read Consequently, the modified Friedmann equation (30) becomes whose most general solution reads [90] where c 1 and c 2 are integration constants and 1 F 1 is the confluent hypergeometric functions or Kummer's function, see definition in chapter 13 of reference [92]. See Table 3 for the possible values of A we shall consider here.

Viability and Local System Tests
It is a well-known fact that local tests pose rather tight constraints on the metric formulation of f (R) theories of gravity (see, for example, references [93][94][95][96][97][98][99]). Thus, any candidate for a reliable alternative to GR should pass or, somehow, evade these low-curvature-regime tests (see also [100][101][102][103][104][105] for an interesting discussion). However, the metric f (R) models presented here were expressly built to reproduce a high curvature regime very different from that of (an effective) ΛCDM. Please note that the contributions of the matter (dark and baryonic) and the spatial curvature k to equations (25) and (30) have been de facto ignored. Therefore, the resulting f (R) theories obtained here must be seen as useful asymptotic models for the theoretical evaluation of the quantum fate of classical singularities rather than complete proposals for viable alternatives to GR at all scales.

f (R) Quantum Geometrodynamics
Even though there is a lack of consensus on what is the correct quantum theory of gravity, the application of ordinary quantum mechanics to the universe as a whole leads to an interesting framework, known as quantum cosmology (for a review of the topic see, e.g., references [106,107]). Currently there are multiple proposals to quantize the cosmological background. A non-exhaustive listing of different approaches to quantum cosmology is string theory, loop quantum cosmology [108,109], causal dynamical triangulation [110][111][112] and quantum geometrodynamics [74], among other examples (see also reference [106]). Nonetheless, in this review we shall focus only on the latter approach, which corresponds, in fact, to one of the first attempts to quantize cosmological backgrounds [74]. In the DeWitt's pioneering work [74], a quantization procedure for a closed Friedmann universe was provided, leading to the first minisuperspace model in quantum cosmology. This quantum cosmology is based on a canonical quantization with the Wheeler-DeWitt equation for the wave function Ψ of the universe playing a central role [74,113,114]. The expression minisuperspace stands for a cosmological model truncated to a finite number of degrees of freedom. This nomenclature is derived from the usage of "superspace" to denote the full infinite-dimensional configuration space of GR and the prefix "mini" for the truncated versions. Moreover, DeWitt also proposed a criterion for the avoidance of classical singularities within the quantum regime, namely the DW criterion. This is, the classical singularity is potentially avoided if the wave function of the universe vanishes in the nearby configuration space. This criterion is based on a generalization of the interpretation of the wave function in ordinary quantum theory, where the wave function is the fundamental building block for any observable. Consequently, regions of the configuration space that lay outside of the support of Ψ are, therefore, irrelevant in practice 3 . It should be noted, however, that the non-vanishing of the wave function does not necessary entail a singularity. Therefore, the DW criterion can only be a sufficient but not necessary criterion for the avoidance of singularities. This criterion has been successfully applied in several cosmological scenarios, see, e.g., references [65,66,[71][72][73]89,90,[115][116][117] among others.

Modified Wheeler-DeWitt
Following the ideas presented in reference [75], the Wheeler-DeWitt equation must be adapted to the framework of f (R) theories of gravity. It is worthy to note that when investigating the canonical quantization of an f (R) theory we are interpreting that theory as a fundamental theory of gravity, rather than an effective framework coming from a quantum gravity proposal. As this classical theory of gravity implies the occurrence of singularities, one concludes that it must be quantized. In this section, we shall show how the formalism of quantum geometrodynamics can be adopted to perform such quantization. The resulting scheme is known as f (R) quantum geometrodynamics.
For the gravitational interaction being that provided by metric f (R) theories of gravity, cosmological models can be described by the action in the so-called Jordan frame. For a FLRW background metric, the above action reduces to where the point-like Lagrangian reads denoting by V (3) the spatial three-dimensional volume. Please note that metric f (R) theories of gravity have an additional degree of freedom when compared with GR (see Einstein's and Jordan's frame formulation in [118][119][120][121] and references therein). Consequently, a new variable can be introduced for the canonical quantization of these alternative theories of gravity. Furthermore, this new variable can be selected in such a way to remove the dependence of the action (54) on the second derivatives of the scale factor. Hence, following the line of reasoning presented in reference [75], we select the scalar curvature, R, to be the new variable. Then, the action (54) becomes However, since the scalar curvature and the scale factor are not independent (at the classical level), their relation needs to be properly introduced in the theory via a Lagrange multiplier, µ, for the classical constraint R = R(a,ȧ,ä) given in equation (25). Thence, The Lagrange multiplier can be determined by varying the action with respect to the scalar curvature. This leads to Accordingly, the point-like Lagrangian (57) can be reformulated as For the sake of the quantization procedure, the derivative part of the above point-like Lagrangian can be diagonalized by the introduction of a new set of variables [75]. These are being R a constant needed for the change of variables to be well-defined. There are different proposals among the existing literature for the value of this constant, see, for instance, references [75,[89][90][91]. In this work we adopt the convention discussed in reference [91], where R is defined as the value of the curvature scalar evaluated at some future scale factor a = a on which the description of the universe by means of DE only becomes appropriate. That is where R(a,ȧ,ä) obeys the classical relation (25). Furthermore, for the sake of concreteness, we can safely assume a = 100a 0 hereafter. Since at that moment in the expansion the matter content will be diluted by a factor of 10 −6 with respect to the present concentration and, therefore, Ω DE ≈ 1. Please note that from this definition it follows R > R , since the scalar curvature asymptotically approaches an increasing function with the expansion of the universe. For different definitions of R see references [75,89,90]. In these new variables, equation (59) transforms into 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 of the point-like Lagrangian has been diagonalized, the derivation of the corresponding Hamiltonian is straightforward. The conjugate momenta are Then, the corresponding Hamiltonian reads For the quantization of the theory, we assume the procedure Therefore, the Hamiltonian (65) of the classical system is promoted to a quantum operator. Consequently, the classical Hamiltonian constraint becomes the modified Wheeler-DeWitt (mWDW) equation for the wave function Ψ of the universe [74,75,106]. That iŝ After simple manipulations, the preceding expression can be cast in the form of the hyperbolic differential equation [75] where the factor-ordering parameter has been set equal to zero. The effective potential entering the preceding equation is given by being λ := R /(12V (3) f R ). Please note that for a given f (R) expression, the variables q and x are univocally fixed. Then, the relation in (60b) must be reversed to express f and R f R in terms of x. However, this may not always be possible analytically, limiting the non-numerical evaluation of the mWDW equation (69) to some suitable group of f (R) theories. In the next sections, we will review the results for the wave function Ψ presented in the literature when considering the f (R) expressions previously found by means of the reconstruction techniques (see Sections 3.1-3.3). We recall that those metric f (R) theories of gravity lead to the same asymptotic background expansion as their respective general relativistic phantom DE models and, therefore, they predict (at the classical level) singular cosmological behaviors. Accordingly, we will neglect the contribution of the spatial curvature k in the effective potential (70) when studying the asymptotic fate of those models, since it will be quickly redshifted with the superacceletared expansion. Before proceeding further, we want to address some comments on the structure of the mWDW equation (69). First, we want to emphasize the well-known ambiguity in the factor ordering in equation (69), i.e., we could have chosen a different factor ordering when applying the quantization procedure on the Hamiltonian of the classical theory, which could have led to a different wave function of the universe. According to reference [75], a variation of the factor ordering affects only the pre-exponential factor of the semiclassical wave function. Thus, the actual value of this parameter can be considered to be unimportant for the evaluation of the DW criterion, which is the main goal of this review. In fact, the DW criterion was found to be satisfied independently of the chosen factor ordering for some particular models related to the topic of this review, see references [72,73,122,123]. On the other hand, note that the mWDW equation (69) is a globally hyperbolic differential equation, i.e., the signature of the minisuperspace DeWitt metric, is (+, −). This is quite different to what happens in GR when the DE content is described by a minimally coupled phantom scalar field. In that case, the DeWitt metric has a positive signature and, therefore, the WDW equation is of an elliptic type. Examples of these models can be found, for instance, in references [66,[71][72][73]. Additionally, a change of signature in the WDW equation has also been noticed in the presence of non-minimally coupled scalar fields [124]. Moreover, the change in the signature has implications in the imposition of boundary conditions. Whereas the hyperbolic equation has a wave-like solution and, therefore, a wellposed initial value problem, a perturbation in the initial or boundary condition for an elliptic (or parabolic) equation will spread instantly over all the point in that domain.

The BR in f (R) Quantum Cosmology
In sight of the form of expression (37), only the negative branch gives the expected limit to a de Sitter universe when the parameter A shrinks to zero. (On the contrary, the exponent γ + diverge when A vanishes) Since small deviations from the EoS of a cosmological constant are the observationally preferred situation today [57] (see also references [37][38][39][40][41][42][43][44][45], among others), hereafter we consider only the negative branch of the solution presented in equation ( 36). This is, we assume c + = 0 in expression (36). Therefore, in this section, we quantize a subclass of the more general family of metric f (R) theories of gravity predicting the classical occurrence of a BR singularity. That is In the following of this section, we drop the subindex "-" for the sake of the notation but keeping in mind that the preceding expression corresponds only to one of the two independent solutions obtained from the reconstruction procedure. Additionally, note that cosmological constraints on the model [57] (see also references [37,38]) satisfy the condition A ≤ (10 − 4 √ 6)/3 ≈ 0.067. Therefore, those constraints favor a real valued exponent in equation (72) [see discussion below equation (37)]. Hence, hereafter we shall consider that the aforementioned exponent is real for values of A of physical interest. Additionally, note that the quantum fate of the BR singularity in f (R) gravity has already been studied in the literature for some particular values for A, finding that the DW criterion can be satisfied, see reference [89].
For the choice of the f (R) gravity in (72), the change of variables in equation (60) reads Moreover, from the evolution equations for the Hubble rate and its cosmic time derivative, and the definition adopted in expression (61), it follows that the constant R is Accordingly, the effective potential (70) entering the mWDW equation (69) becomes where, for the sake of the calculation, we have adopted the notation Thus, the mWDW equation (69) reads However, the dependence of the effective potential on the minisuperspace variables can be simplified when considering the following change of variables [89] (see also reference [73]) Please note that this implies [73] Then, for the choice of r(z) = e Cz/6 , the potential term will only depend on θ [89]. Accordingly, the mWDW equation (77) simplifies to where the factor 1 − C 2 /36 is different from zero at least for the range of values for A of our interest, i.e., A ≤ 0.067 [see discussion below equation (72)]. Therefore, to analyze whether the BR singularity is avoided in the quantum realm by means of the DW criterion, we focus on solving the preceding equation for the wave function Ψ in the configuration space near the singularity. As we do not expect the wave function to be peaked along the classical trajectory in this regime, R and a may take completely independent values. Therefore, to consider a region close to the BR singularity, we should assume either a → ∞ or R → ∞. Both choices imply θ → ∞, but in the former case z is arbitrary whereas in the latter one z → ∞. Please note that the divergence of the scalar curvature can be argued to be the dominant condition, from a geometric point of view, for the occurrence of the BR singularity. Thence, we shall consider both θ and z going to infinity as the main parametrization of the BR singularity in the configuration space. Nevertheless, the results and conclusions presented in this section are independent of this choice and still hold for θ → ∞ and z arbitrary. Additionally, further simplifications can be made when solving (81) close to the BR. By considering the third term containing the cross partial derivatives to be unimportant when θ → ∞, the above equation can be solved via the separation ansatz where b˜k gives the amplitude of each solution andk is related to the associated energy. Please do not confusek with the spatial curvature k, which has been set to zero. (The validity of this approximation is analyzed in Appendix A.) Under these approximations, equation (81) reduces to the following system of equations 4h The former equation can be straightforwardly worked out. The solutions are being d 1 and d 2 arbitrary constants. This corresponds to either exponential or trigonometric functions on z, depending whetherk 2 is positive or negative, respectively. On the other hand, equation (84) can be solved in an exact way by means of Bessel functions, cf. 9.1.53 of reference [92]. This solution can be expressed as where u 1 and u 2 are integration constants, J ν and Y ν are the Bessel functions of first and second order, respectively, and When near the BR singularity, i.e., for large θ, the χ˜k part of the wave function Ψ behaves asymptotically as whereũ 1 andũ 2 now depend onk, cf. 9.2.1-2 in reference [92]. Thence, the total wave function of the universe asymptotically reads As a result, for the condition of d 1 = 0 whenk 2 positive, the wave function vanishes at the BR singularity. Therefore, for this boundary condition, the DW criterion is satisfied. This result points towards the avoidance of this cosmological singularity in the quantum real of f (R) cosmology. Nevertheless, it should be noted that by setting d 1 = 0 we have dismissed a subgroup of solutions to the mWDW equation as unphysical. If future investigations show the importance of these ignored solutions, then it would be concluded that the DW criterion may not always be satisfied.

The LR in f (R) Quantum Cosmology
In this section, we address the quantum fate of the LR abrupt event in the framework of f (R) quantum geometrodynamics. Previously, in Section 3.2, it was shown that the alternative theory of gravity given by the function (46) is the most general expression for a metric f (R) theory of gravity which gives the same asymptotic expansion history as GR filled with phantom DE governed by equation (38). Moreover, since the presence of the Ei function in the term multiplying c 2 spoils the analytical inversion of x(R) in equation (60b), which is crucial for the exact derivation of the mWDW equation (69), hereafter we consider c 2 = 0 [91]. Thus, we focus on the simple, still general, group of alternative metric f (R) theories of gravity predicting the classical occurrence of a LR abrupt event given by [91] where we recall that β = √ 3A/2 and the parameter A is observationally constrained to approximately 2.75 × 10 −28 m −1 , see Table 3. Given the expression for f (R), the change of variables (60) reads being f R = 2c 1 75β 2 − 9β 9β 2 + 12R + 2R . The definition of the constant R needed for the above change of variables to be well-defined was previously discussed for the most general case, see equation (61). Consequently, from equations (41) and (21) it follows [91] On the other hand, for the purpose of computing the effective potential entering the mWDW equation, the inverse of the relation (92b) applies. This is However, only the positive branch is compatible with R being an increasing function of x and R > R . Thus, we choose the positive sign in the preceding expression. Thereafter, the effective potential in the mWDW equation becomes [91] V(q, x) = −U(x)q 6 , where beingŨ := f R /(48c 1 λ 2 R ) a constant. Contrary to the previous section, a change of variables such as (78) will no longer be able to make the potential one-variable-dependent and, therefore, a separation ansatz such as (82) is not appropriate in this case. Instead, note that the U(x) part of the effective potential converges very quickly to a constant value when the model is observationally constrained. This feature suggests an adiabatic semiseparability 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 [125] and first introduced in the framework of quantum cosmology in references [126][127][128]. Furthermore, in reference [91] we have argued that the scalar curvature can be considered more fundamental from a geometrical point of view than the scale factor, justifying that the following ansatz à la Born-Oppenheimer should apply where we recall that x depends only on R, whereas q contains both R and a, see the definitions in (60). Additionally, b˜k represents the amplitude of each solution andk is related to the associated energy. (Please do not confusek with the spatial curvature k which has been set to zero). As a result of this ansatz, the mWDW equation (69) becomes The contribution of the second and third terms in the above expression can be neglected by virtue of the adiabatic assumption 5 . Therefore, equation (98) implies The former equation can be directly solved. The solutions are exponential and trigonometric functions on x, depending on the sign ofk 2 . These are [91] ϕ˜k(x) = d 1 exp where d 1 and d 2 are integration constants. The solutions for χ˜k, on the other hand, are obtained assuming that the potential U(x) behaves as a quasiconstant. Please note that this approximation is based on the fact that U(x) converges very quickly to a constant value when observational constraints on the parameter A are taken into account. Hence, the solutions are [91] χ˜k(q, 5 The validity of this approximation is checked in Appendix B.
being u 1 and u 2 integration constants, cf. 9.1.53 of [92]. Thus, near the LR abrupt event, when both q and x diverge 6 , the resulting wave function of the universe becomes where the integration constantsũ 1 andũ 2 now depend onk 2 , see discussion in reference [91].
Since U(x) > 0 and tends to a constant value when x grows, the wave function cancels at the LR abrupt event when one of the integrations constants is set to zero. This is d 1 = 0 wheñ k 2 is positive. Hence, the DeWitt criterion can be, indeed, satisfied. This points towards the avoidance of the LR abrupt event in the quantum realm of metric f (R) theories of gravity. Nevertheless, as discussed for the case of the BR in the previous section, the fulfilment of the DW criterion is conditioned to the cancellation of one of the integration constants in Ψ. Accordingly, if future investigations claim for the physical importance of the dismissed solution, then it would have to be concluded that the DW criterion may not always be satisfied.

The LSBR in f (R) Quantum Cosmology
Finally, we present the quantum fate of the LSBR abrupt event predicted in the metric f (R) theories of gravity (52) obtained in Section 3.3. Moreover, following a line of reasoning similar to that presented in the previous sections, we consider c 2 = 0 in (52) as a necessary assumption to analytically obtain the corresponding mWDW equation. Therefore, we focus on the simple, still general, f (R) cosmological model exhibiting a future LSBR abrupt event given by This model has been already studied in the f (R) quantum geometrodynamics approach, see reference [90]. In the next part of this section, we not only review the conclusion there presented but also provide a different, less restrictive and more general, approximation when solving for the wave function. For the f (R) gravity model (104), the change of variables (60) reads being f R = 2c 1 (R − 9A). Following the spirit for a physically meaningful definition of the constant R given in equation (61), and in view of equations (23) and (49), we use where we recall that a = 100a 0 has been set, for the sake of concreteness, as the moment in the expansion history from which we can safely assume that DE is the only relevant component of the universe [see discussion below equation (61)].
Subsequently, the mWDW equation for the particular f (R) expression considered in equation (104) reads where W(x) is withW := f R /(24c 1 λ 2 R ) a constant. It should be noted that (107) resembles the form of the corresponding mWDW equation for the LR case presented in the previous section. In fact, the W(x) part of the effective potential also converges quickly to a constant value when A is observationally constrained. Therefore, we consider here the very same BO approximation discussed in the previous section 7 . That is ansatz (97), where b˜k represents the amplitude of each solution andk is related to the associated energy, do not confuse it with the spatial curvature k which has been neglected. Accordingly to the results presented in the previous section, mutatis mutandis, the total wave function of the universe near the LSBR abrupt event reads whereũ 1 andũ 2 now depend onk 2 . The validity of the BO approximation in this case is verified in Appendix B. Since W(x) > 0 and tends to a positive constant value when x grows, the wave function cancels at the LSBR abrupt event when one of the integrations constants is set to zero. This is d 1 = 0 whenk 2 is positive. Hence, the fulfilment of the DeWitt criterion points towards the avoidance of the LSBR abrupt event in the quantum realm of metric f (R) theories of gravity. Nevertheless, as stated before, if future investigations find out the importance of the dismissed solutions, then it would be concluded that the DW criterion may not always be satisfied.

Conclusions
The BR, LR and LSBR are cosmic curvature doomsdays predicted in some cosmological models where the superaccelerated expansion of the universe leads to the disintegration of all bounded structures and, ultimately, to the tear down of space-time itself. Within the context of GR, these cosmological catastrophes are intrinsic to phantom DE models, although a phantom fluid may also induce other cosmic singularities. For a FLRW background, these events can be characterized by the behavior of the scale factor a, the Hubble rate H and its cosmic time 7 Please note that in reference [90] we have used a different ansatz for solving the mWDW equation (107). Since the exponential terms appearing in (108) are strongly suppressed not only by the observational value of the parameter A, but also by the divergence of x, we have considered only the asymptotic value of the function W(x), i.e., we have approached the effective potential in (107) to asymptotically depend only on the variable q. However, as argued in [91], subdominant contributions to W(x) have imprints in the shape of Ψ near the abrupt event. Those subdominant contributions are, in fact, important when comparing different wave functions that share a common asymptotic regime.
derivativeḢ near the singular fate, see Table 1. More importantly, some of these models have been shown to be compatible with the current observational data, see, for instance, reference [57]. Therefore, our own universe may evolve towards some of these (classical) doomsdays. Nevertheless, quantum gravity effects can ultimately become significant and smooth out, or even avoid, the occurrence of these classically predicted singularities. In that sense, quantum cosmology is the natural framework for addressing the quantum fate of cosmological singularities. Among the different approaches to quantum cosmology, we have focused on the canonical quantization of cosmological background due to DeWitt's pioneering paper [74]. In this scheme, the DW criterion can be understood as a sufficient but not necessary condition for the avoidance of a classical singularity. Thus, the classical singularity is potentially avoided if the wave function of the universe vanishes in the nearby configuration space. This criterion has been successfully applied in GR for the phantom DE models considered in Section 2.1, which predict a classical fate à la BR, LR and LSBR, see references [66,72,73]. Consequently, this hints towards the avoidance of these cosmological doomsdays for those specific phantom DE models.
On the other hand, since the background late-time classical evolution can be equivalently described in the context of GR or within the framework of alternative theories of gravity, it is interesting to wonder whether the DW criterion is still fulfilled in the quantum realm of a different underlying theory of gravity. To find that out, reconstruction methods can be applied to obtain the general group of alternative theories of gravity able to reproduce the same expansion history as that of a given general relativistic model. Thus, in Section 3, we have collected the different f (R) theories of gravity found among the literature that produce the same asymptotic background expansion as the phantom DE models summarized in the preceding Section 2.1. Consequently, these f (R) models predict the classical occurrence of a BR, LR or a LSBR abrupt event. For the evaluation of the quantum fate of these classical models, we have followed the f (R) quantum geometrodynamics approach introduced by Vilenkin [75], where the Wheeler-DeWitt equation is adapted for the case of f (R) theories of gravity. Within this framework, the application of the DW criterion has been discussed, showing that it is possible to find solutions to the mWDW equation with a vanishing wave function for the aforementioned cosmic events. Therefore, as it happens when the gravitational interaction is that provided by GR, this result hints towards the avoidance of these cosmological doomsdays within the scheme of metric f (R) theories of gravity [89][90][91].
It should be noted, however, that the validity of the wave functions obtained in this review is subject to the fulfilment of the conditions discussed in the appendices. Thus, for the avoidance of the BR singularity, for instance, we found our approximations to the wave function of the universe to be legit only for very tiny values of the parameter A, indeed smaller even than those found in reference [57] where a strong ansatz was imposed on the possible values of w. Of course, by relaxing such an ansatz tinier values of A are compatible with cosmological observations; see, for example, references [37][38][39] among others. On the other hand, it should be also noted that we have applied the quantization procedure only to some particular solutions to the reconstruction method, i.e., we have quantized just some specific f (R) functions from the more general family of metric f (R) gravity found to classically predict a future fate à la BR, LR or LSBR. Furthermore, certain boundary conditions have also been imposed when solving for Ψ to obtain vanishing solutions at the singular points. These conditions typically involve the cancellation of one of the integration constants. Therefore, a whole subgroup of solutions to the mWDW equation has been disregarded 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.
The pioneering investigations presented in this paper have raised additional questions that should be addressed. From a classical point of view, the metric f (R) theories under consideration should be further analyzed, investigating their compatibility with different observational constraints, as those coming from Solar System tests. On the other hand, it would be also interesting to consider whether semiclassical effects may be important when approaching the quantum realm and, in that case, how they will affect the disintegration of bounded structures. Furthermore, we should highlight that strictly speaking, the DW criterion is not a sufficient condition to guarantee the avoidance of classical singularities. To have a complete analysis, one should obtain the expectation values of the relevant operators. Nonetheless, the calculation of expectation values and probability distributions is related to various open questions in quantum cosmology that still must be properly addressed in the framework of quantum geometrodynamics. Specifically, the correct boundary or initial conditions, the Hilbert space structure and the classical-quantum correspondence, see reference [106]. Finally, we hope that the topics summarized in this work and the presented open questions will motivated further investigations.
Author Contributions: All authors contributed equally to this review. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

ΛCDM
According to that approximation, the solutions for χ˜k and ϕ˜k are presented in equations (85) and (86). Therefore, the terms we kept in equation (81) However, the terms that we have neglected behave asymptotically as (A6) We recall that the integration constants d 1 must be set to zero to fulfil the DW criterion at the BR whenk 2 positive. Thus, to obtain the compliance region of the performed approximation we compare the largest of the neglected terms with the smallest of the saved ones. This is the ratio = B 3h 2 θ∂ θ χ˜k∂ z ϕ˜k This ratio keeps below one if γ − 2 is sufficiently small to compensate the increase of the variable θ towards the BR singularity. Therefore, for γ ≈ 2, this approximation is valid throughout the semiclassical regime towards the BR singularity, where θ increases but not sufficiently rapidly to compensate the small value of γ − 2. In the expression (36) for the metric f (R) theory of gravity predicting the BR singularity, this would correspond to having a small parameter A and c + = 0, since γ + diverge at A → 0. Consequently, this argument would favor small deviations from the ΛCDM model. It is worth noting that the values estimated for A in reference [57] are not small enough to ensure 1; however smaller values for A are compatible with the fits claimed in other references (see, for example, references [37][38][39]). Thus, we consider this approximation to Ψ valid for the appropriate γ value.
Additionally note that in this section, we have enhanced the discussion on the viability of the performed approximations for Ψ originally presented in reference [89], where the authors have focused only on the particular case ofk = 0 when addressing the validity of the wave function there found.
For the case of the LR abrupt event analyzed in Section 4.3 that ratio reads Consequently, the approximation is valid as long as ε 1. To evaluate this condition, note that when β is observationally constrained, see Table 3. Thus, in the configuration space near the LR cosmic event, we have where both q and x diverge. Finally, ε 1 near the LR if β is sufficiently small, i.e., for small value of A. Please note that this corresponds, in fact, to the observationally preferred situation [57]. (We recall that A is of order 10 −28 m −1 when observationally constrained, see Table 3). Therefore, when the parameters of the theory are observationally constrained, the approximation is valid throughout the semiclassical regime towards the LR abrupt event; where the variables q and x increase but not sufficiently rapidly to compensate the small value of β 2 . Hence, for the purpose of the present work, i.e., to analyze the fulfilment of the DW criterion in the configuration space close to the LR, this approximation is valid.
On the other hand, for the verification of the validity of the BO approximation performed in Section 4.4, U(x) must be exchanged for W(x) in the preceding formulas. Therefore, considering the expression for W(x) given in (108), the ratio ε becomes Following the same line of reasoning as that presented before, this ratio keeps below 1 since the parameter A for the LSBR is of order 10 −54 m −2 , see Table 3. Therefore, the BO approximation applied in Section 4.4 is valid throughout the semiclassical regime towards the LSBR abrupt event.