Future dynamics in f(R) theories

The $f(R)$ gravity theories provide an alternative way to explain the current cosmic acceleration without invoking dark energy matter component. However, the freedom in the choice of the functional forms of $f(R)$ gives rise to the problem of how to constrain and break the degeneracy among these gravity theories on theoretical and/or observational grounds. In this paper to proceed further with the investigation on the potentialities, difficulties and limitations of $f(R)$ gravity, we examine the question as to whether the future dynamics can be used to break the degeneracy between $f(R)$ gravity theories by investigating the future dynamics of spatially homogeneous and isotropic dust flat models in two $f(R)$ gravity theories, namely the well known $f(R) = R + \alpha R^{n}$ gravity and another by A. Aviles et al., whose motivation comes from the cosmographic approach to $f(R)$ gravity. To this end we perform a detailed numerical study of the future dynamic of these flat model in these theories taking into account the recent constraints on the cosmological parameters made by the Planck team. We show that besides being powerful for discriminating between $f(R)$ gravity theories, the future dynamics technique can also be used to determine the fate of the Universe in the framework of these $f(R)$ gravity theories. Moreover, there emerges from our numerical analysis that if we do not invoke a dark energy component with equation-of-state parameter $\omega<-1$ one still has dust flat FLRW solution with a big rip, if gravity deviates from general relativity via $f(R) = R + \alpha R^n $. We also show that FLRW dust solutions with $f''<0$ do not necessarily lead to singularity.

The f (R) gravity theories provide an alternative way to explain the current cosmic acceleration without invoking dark energy matter component. However, the freedom in the choice of the functional forms of f (R) gives rise to the problem of how to constrain and break the degeneracy among these gravity theories on theoretical and/or observational grounds. In this paper to proceed further with the investigation on the potentialities, difficulties and limitations of f (R) gravity, we examine the question as to whether the future dynamics can be used to break the degeneracy between f (R) gravity theories by investigating the future dynamics of spatially homogeneous and isotropic dust flat models in two f (R) gravity theories, namely the well known f (R) = R + αR n gravity and another by A. Aviles et al., whose motivation comes from the cosmographic approach to f (R) gravity. To this end we perform a detailed numerical study of the future dynamic of these flat model in these theories taking into account the recent constraints on the cosmological parameters made by the Planck team. We show that besides being powerful for discriminating between f (R) gravity theories, the future dynamics technique can also be used to determine the fate of the Universe in the framework of these f (R) gravity theories. Moreover, there emerges from our numerical analysis that if we do not invoke a dark energy component with equation-of-state parameter ω < −1 one still has dust flat FLRW solution with a big rip, if gravity deviates from general relativity via f (R) = R + αR n . We also show that FLRW dust solutions with f < 0 do not necessarily lead to singularity.

I. INTRODUCTION
A wide range of cosmological observations coming from different sources, including the supernovae type Ia (SNe Ia) [1], the cosmic microwave background radiation (CMBR) [2] and baryon acoustic oscillation (BAO) surveys [3], clearly indicate that the Universe is currently expanding with an accelerating rate. A fair number of models and frameworks have been proposed to account for this observed accelerated expansion. These approaches can be roughly grouped in two families. In the first, the so-called dark energy is invoked and the underlying framework of general relativity (GR) is kept unchanged. In this regard, the simplest way to account for the accelerating expansion of the Universe is through the introduction of a cosmological constant, Λ, into Einstein's field equations. This is entirely consistent with the available observational data, but it faces difficulties such as the order of magnitude of the cosmological constant and its microphysical origin. In the second family, modifications of Einstein's field equations are assumed as an alternative for describing the accelerated expansion. This latter group includes, for example, generalized theories of gravity based upon modifications of the Einstein-Hilbert action by taking nonlinear functions, f (R), of the Ricci scalar R or other curvature invariants (for reviews see Refs. [4]).
The fact that f (R) theories can potentially be used to explain the observed accelerating expansion have given birth to a number of articles on these gravity theories, in which several features of f (R) gravity have been discussed [5], including the stability conditions [6], compatibility with solar-system tests [7], energy conditions [8], nonlocal causal structure [9], and observational constraints from a diverse set of cosmological observations [10].
However, although the freedom in the choice of the functional forms of f (R) has motivated many different suggestions of f (R) gravity theories, which account for the accelerating expansion and are compatible with the solar-system tests, it also gives rise to the problem of how to constrain and break the degeneracy among these gravity theories on theoretical and/or observational grounds. In this regard, observational constraints on some f (R) gravity from a diverse set of observations have been placed [10], and tests of the cosmological viability of some specific forms of f (R) have been explored [11].
A pertinent question that arises here is whether the future dynamics can be used to break the degeneracy between f (R) gravity theories. In this article, to proceed further with the investigation on the potentialities, difficulties and limitations of f (R) gravity, we examine this question by investigating the future dynamics of Friedmann-Lemaître-Robertson-Walker (FLRW) dust flat model in two f (R) gravity theories, namely the well known f (R) = R + αR n gravity, for which many results are available in the literature [12], and another by A. Aviles et al. [13] (ABCL gravity for short), whose motivation comes from the cosmographic approach to f (R) gravity [14,15]. We show that besides being powerful for discriminating between f (R) gravity theories, the future dynamics technique can also be used to determine the fate of the Universe in f (R) gravity theory.
Until the discovery of the accelerating expansion virtually any textbook on cosmology describes the future dynamics of FLRW pressure-free dust models as follows. It expands forever if it has an Euclidean or hyperbolic spatial geometry, and expands and eventually recollapses if it has a spherical spatial geometry. However, the discovery of the accelerating expansion made apparent that these simple future forecasts had to be modified, since the negative-pressure dark energy component, invoked to account for the acceleration, plays a crucial role in the evolution of the Universe. Indeed, the dark energy (DE) is usually described by the equation-of-state parameter ω = p/ρ, which is the ratio of the DE pressure to its density. A value ω < −1/3 is required for cosmic acceleration. When −1 < ω < −1/3 the DE density decreases with the scale factor a(t). However, if ω < −1 the dark energy density becomes infinite in a finite-time, driving therefore the Universe to a future finite-time singularity, called big rip [16]. Afterwards it was determined that this is not the only possible doomsday of a dark energy dominated universe. It may, for example, come to an end in a sudden singularity [17], a big freeze doomsday [18], or a little rip [19]. In this paper we also determine that even if we do not invoke a dark energy component with ω < −1 one still has pressure-free dust solution with a big rip, if gravity deviates from general relativity via f (R) = R + αR n . Finally, by using the future dynamics scheme of this paper, we further present an example of FLRW dust solution in which the ghost-like regimes (f < 0) do not necessarily lead to singularity 1 .
Our paper is organized as follows. In Sec. II we give a brief review of f (R) gravity theories, derive the field equations for the flat FLRW metric with dust matter content, state the initial conditions for the dynamical evolutions, and present the future dynamics in the context of the general relativity theory. In Sec. III we introduce the ABCL [13] Lagrangian, develop the necessary technique for solving the dynamical equations, and derive our 1 We note that future dynamics in f (R) gravity was also considered in Ref. [20]. Their approaches are different from ours since they have considered perfect fluid matter source with ρ, p and ω = p/ρ, while we consider simply a pressure-free dust. Furthermore, our analysis is numerical while theirs is not.
numerical results regarding the ABCL f (R) gravity. In Sec. IV we use our method to study the polynomial Lagrangian f (R) gravity and make a comparative analysis of these theories. Final remarks and main conclusions are presented in Sec. V.

II. PREREQUISITES
In this section we briefly review the f (R) gravity, derive the field equations for the flat FLRW metric with dust matter content, state the initial conditions for all numerical analyses of this paper, and present the future dynamics for the particular Einstein's gravity theory for later comparison with the dynamics in other gravity theories.

A. f (R) gravity and field equations
We begin by recalling that the action that defines an f (R) gravity theory can be written as where κ 2 ≡ 8πG/c 4 , g is the determinant of the metric g ab , f (R) is a function of the Ricci scalar R, and S m the standard action for the matter fields. Varying this action with respect to the metric we obtain the field equations where here and in what follows primes denotes differentiation with respect to R and ≡ g ab ∇ a ∇ b . Clearly, for f (R) = R+Λ these field equations reduce to the Einstein equations with the cosmological constant, Λ, term. Two important constraints, often used to simplify the calculations, come from the fact that the covariant divergence of both sides of Eq. (2), is null. This implies that the 0i and 00 components of the field equations give rise to the following constraints Clearly, Eq. (3) is identically satisfied, while the constraint given by Eq. (4) must be fulfilled throughout time evolution. We shall use this fact as a way of checking the accuracy of the numerical integration of the remaining equations of the dynamical system in the numerical analyses of the following sections.
In this work we focus on the flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric, which is supported by the recent observations [2], and is consistent with the standard inflationary models. Thus, the non vanishing component of the Ricci tensor and the Ricci scalar can be written in the form where H =ȧ/a and the over-dot denotes derivative with respect to the time t.
Since we are interested in dust of density ρ with zero pressure (p = 0) then we havė where here and in what follows the subscript zero denotes present-day values of the cosmological parameters. Taking into account equations (6)-(9) the field equations (2) reduce to One can easily show that Eq. (10) is nothing but the constraint Eq. (4), which for the dust flat FLRW models takes the form which is in a suitable form for checking the accuracy of the numerical integration of the dynamical Eq. (11) for the f (R) gravity theories we are concerned with in this paper.

B. Initial conditions
To study the future dynamics for the spatially flat FLRW dust models one needs to choose initial conditions for the numerical integration. In this work we use the numerical values of the cosmological parameters reported by the Planck Collaboration team [2] along with the values of the cosmographic parameters given in Ref. [22]. In Table I we collect together the values of the cosmological parameters we shall employ in our numerical analyses.   The values of the cosmological parameters, deceleration and jerk parameters, are taken from the Planck results [2] from Ref. [22]. Table I also contains details of the units and convention we have adopted in this paper.
To investigate the future dynamics in the following sections, we recall that for the flat FLRW models the dimensionless deceleration (q) and jerk (j) parameters [22] given in Table I are such that the relations hold.

C. Dynamics in General Relativity
For a later comparison with the dynamics in other gravity theories, we briefly present here the analysis for the spatially flat FLRW dust model in the Einstein theory with cosmological constant Λ, that is for f (R) = R + Λ. In this case the field equations (12) and Eq. (11) reduce, respectively, to The left panel in Figure 1 shows the evolution of H for spatially flat FLRW dust model in the Einstein theory, where initial conditions given in Table I were employed in the numerical calculations. This panel also shows a de Sitter asymptotic behavior for H for a non vanishing cosmological constant Λ. The right panel shows the constraint equation (12), providing an assess of the reliability of the numerical integration.  (12) making apparent the high level of accuracy in the numerical integration has been obtained.
The current values of the cosmological parameters as given by Planck team [2] were taken as initial conditions in the numerical integration.

III. DYNAMICS IN ABCL GRAVITY
In this section we shall study the future dynamic of the spatially flat FLRW dust model in the f (R) gravity recently suggested by A. Aviles et al. [13], referred to in this paper as ABCL gravity theory. This gravity theory has been obtained through an optimal Monte-Carlo fitting of cosmographic results and is given by where a, b and c are free parameters, and R 0 is the present-day value of the Ricci scalar. Regardless of the values of these parameters for this theory one has where f 0 ≡ (∂f /∂R) R=R0 and similar notation is used for higher order derivatives. In addition to the constraints (18) and (19), which are required to ensure that both Einstein's theory and Newton's constant are recovered in the lowest order, we shall take into account that f ≥ 0, which is a condition to avoid the presence of ghosts [23]. We also note that the conditions (18)- (20) also insure that (10) and (11) reduce to the Friedmann equations in the lower order.
In order to study the dynamics we will need Before proceeding to the numerical analysis for this gravity theory, we note that due care ought to be taken in using the initial conditions since f = 0 at t = t 0 . Note that in equations (10) and (11) the higher derivatives of the scalar factor are multiplied by f . When f = 0, we have a set of differential equations describing the dynamics. When f = 0 locally, at t = t 0 we have an entirely different set of differential equations obtained by (10) and (11). In what follows, to deal with this difficulty we first assume that the solution possesses a Taylor expansion about t = t 0 up to second order. Then we substitute this expansion into the field equations (10) and (11) in an order by order manner, which results in a perturbative solution up to some order. Second, instead of assuming the initial condition exactly at t = t 0 , the initial condition is taken at t = t 0 + ( 2 ≪ 1) through this perturbative scheme, for which now f (t 0 + ) = 0.
To carry out the outlined perturbative procedure, it is required to distinguish two different regimes in establishing the initial condition (f = (t 0 + ) = 0), namely one when at t = t 0 , f = 0, and another when t = t 0 , f = 0. In the following we shall treat separately these two cases.
A. The case f 0 = 0 Since for the ABCL gravity theory f (t 0 ) = 0, in order to find a suitable form for the field equations within a perturbative scheme, we first assume that the solution can be expanded in a Taylor series about t = t 0 up to second order. Thus, to obtain the terms of (10) and (11) we have Then we substitute these terms into the field equations (10) and (11) in an order by order mode to have • Zero order in (t − t 0 ). Equations (10) and (11) give • First order in (t − t 0 ). Equations (10) and (11) yield From (22)- (25) one has where from equation (8) we havė As mentioned, if f = 0 the initial condition follows directly from eqs. (10) and (11). Otherwise, if f 0 = 0 at t = t 0 , then the initial condition is chosen for t near t 0 , so that we can be sure that f (t) = 0 and the dynamics is directly described by eqs. (10) and (11). In this case, instead of taking H 0 ,Ḣ 0 andḦ 0 the initial conditions are chosen at t − t 0 = whereḢ 0 ,Ḧ 0 , ... H0 satisfy the relations (26)-(28). We note that according to eqs. (27) the initial valueḢ 0 ought to obey the constraint The left panel of Figure 2 shows a representative numerical future dynamics of spatially flat FLRW dust model in the ABCL gravity theory with f 0 = 0. In this case, as indicated by H(t), the universe would present a noteworthy expanding phase after an initial future decelerating period. The right panel in this figure shows the constraint E 00 = 0, given in eq. (12), fluctuates randomly and increases but is always smaller than 2 × 10 −17 , which  Table I as given by Planck team [2] were taken as initial conditions in the numerical integration. Right panel: The behavior of the constraint E00 as given by (12), making apparent the degree of confidence of the numerical calculations. No singularities were found throughout the numerical evolution.
is a strong indication of the correctness of the numerical solution.
Regarding the choice of the a, b and c used to find the numerical future dynamics solution, it is important to point out that the independence of equations (18)- (20) allows some freedom in their choice. Although a more general phase space analysis, with its associated attractors, would be necessary to determine every possible future dynamics solution, here we have restricted our analysis to a set of values which are consistent with the present-day constraints on the cosmological parameters. This choice of values was also motivated by a similar procedure used in the study of the polynomial f (R) gravity, which we study in details in the Section IV.
B. The case f 0 = 0 The condition f 0 = 0 can be achieved in the ABCL f (R) gravity provided that the constraint holds. Now, similarly to the previous section we assume that the solution can be expanded in a Taylor series t = t 0 up to second order. Then we substitute Taylor series terms into the field equations (10) and (11) in an order by order manner. This gives the following: • Zero order in (t − t 0 ). Equations (10) and (11) give • First order in (t − t 0 ). Equations (10) and (11) yield • Second order in (t − t 0 ). Equations (10) and (11) furnish Now, sinceṘ is given by Eq. (29) it is clear from eq. (38) thatḦ 0 is given by a third order algebraic equation. Thus, from the above equations (35)-(40) one has  Table I as given by Planck team [2]. The dynamic evolution of the flat FLRW spacetime tends asymptotically to Sitter space. Right panel: The behavior of the constraint E00 as given by (12), which is smaller than 10 −13 for all time t, making apparent the precision of the numerical calculation throughout evolution of the model. the case f 0 = 0. The left panel of the figure shows the dynamic evolution of this dust flat FLRW model tends asymptotically to de Sitter space. Its asymptote approaches the de Sitter vacuum solution obtained by analytically solving Eq. (2) with no matter content, which returns a value H = 1.8476 × 10 −4 , i.e. the asymptote shown in Figure 3. By carrying out a stability analysis having this solution as background, it is possible to show this is an attractor solution with f > 0. The right panel of Figure 3 makes clear that the constraint on E 00 [Eq. (12)], is smaller than 2 × 10 −13 during the evolution of the FLRW flat model. This makes apparent the accuracy of the numerical calculations performed in the study of the future dynamics in this case.

IV. DYNAMICS IN f (R) = R + αR n
In this section we apply the same numerical scheme used in the previous section to study the well known f (R) gravity theory The above f (R) has the special feature of presenting equivalent results to the ΛCDM model, when applied to a FLRW setting within a range of choices for the constants α and n. In this sense no observationally relevant prediction would distinguish these cases as established in [24].
As referred to in the end of Sec. III A, a straightforward calculation to guide the choice of parameters α and n can be achieved by using the field equations (10) along with equations (43), (13) and (14). Indeed, these equations allow to write α as Where observed values of parameters H 0 , q 0 , j 0 and ρ 0 and their uncertainties are given in Table I. Now, the standard error deviation tainties of the parameters H 0 , q 0 , j 0 and ρ 0 (Table I) one can plot the curves in the panels of Figure 4 to illustrate these bounds 2 . They have been used to guide suitable choices of values of the parameters α and n below.
By applying the same analysis introduced in Section III and using again the main values at Table I as initial conditions for the FLRW dust flat models, we have that for n > 0, many parameter choices lead to asymptotic de Sitter solutions, which is consistent with the results in the literature [12]. As an example, Figure 5 depicts H for of n = 0.2 and α = −4.295 × 10 −6 , which are typical values between the above-mentioned bounds in the left panel of Figure 4. In Figure 5 we also show the behavior of the constraint E 00 , which exhibits fluctuations around zero with rather small amplitudes, making apparent the degree of confidence of the numerical calculations. For this particular value of α, n and initial conditions as given by Planck team, the de Sitter analytical vacuum solution obtained by solving Eq. (2) gives H = 1.8431 × 10 −4 , in rather good agreement with the asymptotic values shown in Figure 5. It is possible to show this solution is an attractor over this phase space region, with f > 0, as is also the case for the example depicted before in Figure 3.
As for the n < 0 case, many initial conditions lead to a big rip singularity. This is a curvature singularity in the sense it cannot be removed by a coordinate transformation. As an example, we show in Figure 6 the time evolution of the Hubble parameter of the dust flat FLRW model for n = −0.1 and α = −4.557 × 10 −8 , chosen by taking into account right panel of Figure 4. In this case, the analytical vacuum solution to the field equation (2) gives H = 1.8618 × 10 −4 . It is also possible to show this solution is a repellor in this phase space region, rendering the divergent evolution shown in Figure  6 an expected feature, complemented by our verification that it presents f < 0 asymptotically. As indicated by the right panel of Figure 6, the numerical constraint E 00 noticeably increases as the solution approaches the physical singularity, which is expected from the outset, and then the numerical solution must be truncated.
It is also interesting to compare the future dynamics shown in Figure 2, for ABCL gravity, with that of Figure 6, for f (R) = R + αR n with n < 0. The dynamic depicted by Figure 6 is a typical example of the run away evolution that suddenly ends in a big rip singularity. This evolution is usually understood in terms of the ghost-like behavior due to its transition to a regime where f < 0. On the other hand, Figure 2 presents another runaway solution, that also develops a ghost-driven regime, but with no sudden singularity within the accuracy of the numerical analysis. Thus, even when f < 0 the evolution in one theory (ABCL gravity) presents the remarkable property of smoothing out the expected divergence. These different behaviors make apparent the richness of possible evolutions in the context of f (R) gravity theories.
We can further examine the connection between big rip singularities and ghost-like regimes (f < 0) through another example. In fact, by taking the same values for the parameters n and α used for Figure 6 along with a slightly different value for the Hubble parameter (H 0 = 2.0 × 10 −4 , in units of Table I), we have calculated the future dynamics of the FLRW dust flat models shown in Figure 7. For this case, sinceḢ → 0 and H → 0, from Eq. (8) one has the Ricci scalar tends to zero, and thus f → −∞. This evolution illustrates a case where, in a limit strongly associated to ghost-dominated regimes (f 0), the solution actually evolves simply to Minkowski spacetime. To the best of our knowledge, this interesting dynamical behavior has not been highlighted so far in the literature. It also illustrates how a direct association of big rip singularities, or even run away solutions, with a ghost-like regime of the field equations can be misleading in the framework of f (R) theories.  Table I. Right panel: The behavior of the constraint E00 as given by (12), which is smaller than 10 −22 for all time t, making explicit the precision of the numerical calculation throughout evolution of the model.  Table I as given by Planck team [2]. Right panel: The behavior of the constraint E00 as given by (12). As expected, the constraint increases as the curvature singularity approaches, where the numerical solution must be halted.

V. FINAL REMARKS AND CONCLUSIONS
There has been a great deal of recent papers on f (R) gravity motivated by the attempts to explain the current cosmic acceleration with no need of invoking a dark energy component. Despite the arbitrariness in the choice of different functional forms of f (R), which call for ways of constraining the possible f (R) gravity theories on physical grounds, several features of these gravity theories have been discussed in a number of recent articles. In this paper we have proceeded further with the investigation of potentialities and limitations of f (R) gravity theories by examining whether the future dynamics can be used to break the degeneracy between f (R) gravity theories. To this end, by taking the recent constraints on the cosmological parameters made by the Planck team, we have performed a detailed numerical study of the future dynamic of spatially homogeneous and isotropic dust flat models in the framework of two gravity theories. As a first result, we have shown that besides being powerful for discriminating between f (R) gravity theories, the future dynamics numerical technique introduced in this paper can also be used to determine the fate of the Universe in the framework of these f (R) gravity theories. Fig-ure 8 collects together the results of the future dynamics numerical analyses of the FLRW dust flat models in several different cases. Curve (a) shows the future evolution of this model in general relativity theory [f (R) = R)], while the other curves represent future dynamics of these FLRW models in several instances as follows. Curve (b) ABCL theory, Eq. (17)] with f = 0 ; Curve (c) ABCL theory with f = 0 ; Curve (d) f (R) = R + αR n with positive n ; Curve (e) f (R) = R + αR n with negative n ; Curve (f ) f (R) = R + αR n with negative n but now with a slightly different value of the Hubble parameter. Figure 8 shows that, although with differences for t 1.5M yr, the ultimate fate of the Universe is a de Sitter model (with slightly different Hubble constant) for the cases (a), (c) and (d). Clearly the cases e evolves to Minkowski space whereas the case f develops singularity. Thus, the future dynamics scheme developed in this paper is indeed a powerful tool to discriminate f (R) gravity theories.
The development of a big rip singularity as shown in curve (e) of Figure 8 is consistent with the results found in the literature for f (R) = R + αR n with negative n, and are generally associated with ghost-like regimes (f < 0). In this regard, an interesting outcome of our  Figure 6, the values of the parameters were taken to be n = −0.1 and α = −4.557 × 10 −8 , but now a slightly different value for H0. Other initial values are given by Table I. It can be seen that Minkowski space is obtained asymptotically when t → ∞, as H → 0 = const. Right panel: The behavior of the constraint E00 as given by (12), which clearly is smaller than 10 −20 for all time t.  Table I as given by Planck team [2]. (a) Standard FLRW dust flat solution in Einstein's equations, which tends asymptotically to a de Sitter universe. (b) ABCL's theory, eq. (17), when initially f = 0, as discussed in subsection III A. (c) A case when initial conditions satisfy f = 0 in ABCL theory, as discussed in subsection III B. (d) Future dynamics of FLRW dust flat model for f (R) = R + αR n with n > 0. (e) Same as (d) but now with n < 0. (f ) The same (α, n) as in (e), but now with a slightly different value of the Hubble parameter H0.
analyses is the evolution given by curve (b) of Figure 8, where another ghost-driven regime follows a smooth accelerated evolution with no associated singularity. Along the same lines, we also note the interesting case presented by curve (f) in the same figure. This is a solution evolving to Minkowski spacetime for f negative and unbounded. Thus, our numerical analyses suggest that ghost-like regimes (f < 0) do not necessarily lead to singularities.
The discovery of the cosmological expansion along with earlier theoretical investigations concerning spatially homogeneous and isotropic models by Friedmann and Lemaître sparked the first scientific studies on the future dynamics and the ultimate fate of the Universe in the framework of Einstein's theory of gravitation. It was shown that if the matter content of the universe is a pressure-free dust, then the future of the Universe would depend only on the sign of the spatial curvature. It would expand forever if it has an Euclidean or hyper-bolic spatial geometry, and would expand and eventually recollapse if it has a spherical spatial geometry. These predictions may be found in virtually any textbook on cosmology, but the discovery of the accelerating expansion, through type Ia supernovae observations, made apparent that this simple future forecast for the Universe does not work anymore, since the negative-pressure dark energy (DE) component, invoked to account for the acceleration, plays a crucial role in the evolution of the Universe. Indeed, the dark energy is usually described by the equation-of-state parameter ω which is the ratio of the DE pressure to its density (ω = p/ρ). A value ω < −1/3 is required for cosmic acceleration. When −1 < ω < −1/3 the DE density decreases with the scale factor a(t). However, it has been shown that if ω < −1 the dark energy density becomes infinite in a finite-time, t s (say), driving therefore the universe to a future finite-time singularity called Big Rip in which when t → t s , then ρ → ∞, a(t) → ∞ and |p| → ∞ [16].
Afterwards, it was understood that this is not the only possible doomsday of a dark energy dominated universe. It may, for example, come to an end in a sudden singularity [17], a big freeze doomsday [18], or a little rip [19]. In this paper we have examined numerically the future dynamics of the Universe but, instead of assuming a dark energy component, we have assumed that gravity is governed by f (R) gravity theories and have kept a pressurefree dust as matter content. As a result, we have also shown that even if we do not invoke a dark energy component with ω < −1 one still can have pressure-free dust FLRW flat solution with a big rip, if gravity deviates from general relativity via f (R) = R + αR n , as shown in Figure 6.
Finally, by using the future dynamics scheme of this paper, we also show an example in which the ghost-like regimes (f < 0) do not necessarily lead to singularity ( Figure 7). Thus, the future dynamics scheme we have developed in this paper is not only a powerful tool to discriminate between f (R) gravity theories, but it has also permitted to shed some light on the ghost-like regime and its connection with singularities in the context of f (R) gravity theories.