Unification of Constant-roll Inflation and Dark Energy with Logarithmic $R^2$-corrected and Exponential $F(R)$ Gravity

In this paper we investigate how to describe in a unified way a constant-roll inflationary era with a dark energy era, by using the theoretical framework of $F(R)$ gravity. To this end, we introduce some classes of appropriately chosen $F(R)$ gravity models, and we examine in detail how the unification of early and late-time acceleration eras can be achieved. We study in detail the inflationary era, and as we demonstrate it is possible to achieve a viable inflationary era, for which the spectral index of primordial curvature perturbations and the scalar-to-tensor ratio can be compatible with the latest observational data. In addition, the graceful exit issue is briefly discussed for a class of models. Finally, we discuss the dark energy oscillations issue, and we investigate which model from one of the classes we introduced, can produce oscillations with the smallest amplitude.


I. INTRODUCTION
Modified gravity in it's various forms, has a prominent role in describing the Universe's evolution [1][2][3][4][5]. Particularly, a vast number of phenomena related to various stages of the Universe's evolution can be explained by modified gravity theories, both at an astrophysical level and also at a galactic level, see for example the review [1]. Among the various theories of modified gravity, F (R) gravity is the most commonly used and, to our opinion, the most appealing modified gravity theory, due to the simplicity offered by F (R) gravity descriptions, but also due to the conceptual rigidity offered by this theory. Among other applications, the F (R) gravity theoretical framework, offers the possibility of a unified description of the early-time and late-time acceleration eras [6], see also Refs. [7] (for reviews on the unification of the early and late-time acceleration, see [1][2][3]). This unified description is one of the most demanding aims of a theoretical description, and various approaches towards this issue have been presented in the literature.
In the present work we aim to present some F (R) gravity models which can describe in a unified way a constant-roll inflationary era with a late-time acceleration era. The constant-roll inflation description, is an alternative approach to the standard slow-roll inflationary era, and it's implications have recently been studied in the context of scalar-tensor theories [8? -21], see also [22][23][24][25] and also in the context of F (R) gravity [26][27][28]. As it was shown in [26], it is possible that non-viable F (R) gravity models in the context of slow-roll inflation, may become viable in the context of constant-roll inflation. As we will demonstrate, in the context of the constant-roll evolution, it is possible to obtain observational indices that are compatible with the Planck constraints. We shall use two well-known F (R) gravity models, namely an R 2 -corrected logarithmic model and a curvature corrected exponential model, and we explicitly calculate the spectral index of primordial curvature perturbations and the scalar-to-tensor ratio, and as we show, these are compatible with the observational data. Also we shall briefly discuss some interesting features of the reheating era corresponding to the R 2 -corrected logarithmic model. In addition, we analyze the late-time evolution corresponding to both the aforementioned F (R) gravity models, and it is noteworthy that these models provide a unified description of the late and early-time acceleration era. Finally, we analyze the dark energy oscillations for the curvature-corrected exponential model, and we demonstrate that the most phenomenologically appealing case corresponds to nearly R 2 curvature corrections. This paper is organized as follows: In section II we introduce and study two models of F (R) gravity, and we investigate the constant-roll inflation implications on the inflationary era. For the R 2 -corrected logarithmic model of F (R) gravity we also briefly discuss the reheating era implications. In section III, we study the late-time evolution implications of the two models we introduced in the previous section, and for the curvature corrected exponential model, we investigate which curvature correction produces less dark energy oscillations during the last stages of the matter domination era and after the deceleration acceleration transition. Finally, the conclusions follow in the end of the paper.
Before we start, let us briefly present the conventions we shall assume. With regard to the geometric background, we assume that it is a flat Friedmann-Robertson-Walker (FRW), with the line element being, where a(t) is the scale factor. In addition, we assume that the metric connection is an affine, torsion-less, and metriccompatible connection, the Levi-Civita connection. Also in the following, the parameter κ stands for κ 2 = 16π/M 2 Pl , where M Pl is the Planck mass scale. The first model we shall consider has the following action, where M denotes the spacetime manifold, g is the determinant of the metric tensor g µν , L m is the Lagrangian of the matter and radiation perfect fluids and R is the Ricci scalar. The corrections to the Hilbert-Einstein Lagrangian of General Relativity (GR) assume the form of F (R)-gravity and are represented by the higher curvature R 2 -term for the early-time inflation, and a function of the Ricci scalar, f DE (R), for the dark energy sector. The running coefficient γ(R) in front of R 2 depends also on the Ricci scalar and has been introduced in order for the graceful exit from inflation to be possible. When γ(R) = γ is constant, then the Starobinsky inflationary scenario is obtained, where the early-time de Sitter expansion is governed by an R 2 gravity. The first Friedmann equation for the model at hand is the following, where the "dot" denotes the derivative with respect to the cosmic time and the "prime" denotes differentiation with respect to the Ricci scalar. The Ricci scalar reads, with H =ȧ(t)/a(t) being the Hubble parameter. In the equation above, ρ m denotes the energy density of matter which satisfies the following conservation law,ρ where p m is the matter pressure. In order to reproduce the early-time acceleration we introduce the following expression for the function γ(R), where R 0 is the curvature of the Universe at the end of inflation and γ 0 , γ 1 are positive dimensional constants. Note that logarithmic corrections of the form appearing in Eq. (II.5) are actually motivated by one-loop corrections to coupling constants in multiplicatively renormalizable higher-derivative quantum gravity, for a general review, see [29]. Also it has been demonstrated [30,31], that such an R 2 -corrected logarithmic gravity predicts viable inflation. Furthermore, these terms may also find a theoretical explanation in the framework of holographic renormalization group flow and the corresponding cosmological implications, in the context these issue were discussed in Refs. [32,33].
In fact, if one incorporates in a gravitational system an holographic surface term, located near to the Hubble horizon, one obtains higher derivative corrections to the FRW field equations. Thus, the solutions of these equations may be used to describe the early-time acceleration and, eventually, by making use of two holographic screens, the late-time accelerated expansion of our Universe today may be realized.
Returning to the model at hand, it turns out that γ(R 0 ) = γ 0 . Since we would like to avoid the effects of R 2 -gravity in the limit of small curvature, we require that the following general condition holds true, where R = 4Λ is the curvature of the Universe when the dark energy is dominant, and Λ is the Cosmological constant. In the following, we will assume that f DE (R) and L m in (II.1) are negligible in the limit of high curvatures. The de Sitter solution with constant curvature R dS = 12H dS follows from (II.2) and it reads, In the case of a constant value of γ(R) = γ 0 , namely for γ 1 = 0, the de Sitter solution is obtained as an asymptotic limit of the first Friedmann equation, when the R 2 -term is dominates the evolution. This is the so-called Starobinsky model, where the Hilbert-Einstein term guarantees a graceful exit from the accelerated phase. Since in the Starobinsky model 1/κ 2 ≪ R, we have a regime of super-Planckian curvature. Here, the R 2 term has the same order of magnitude as the Hilbert-Einstein term during the inflationary era. In this case, the running constant γ(R) in (II.5) with γ 1 = 0, determines the value of the de Sitter solution as in Eq. (II.7). If we perturb the de Sitter solution as follows, by keeping first order terms with respect to δH(t), we obtain from Eq. (II.2), In the limit R 0 ≪ R dS the solution of this equation reads, where h ± are constants depending on the sign of ∆ ± . When we choose the plus sign, the solution diverges and the de Sitter expansion is unstable. Thus, by Taylor expanding, the divergent solution is given by, 11) and in effect, we obtain, where t 0 is the time at the end of inflation when R ≃ R 0 and also h 0 , R 0 and N stand for, In order to study the behavior of the solution during the exit from inflation, we introduce the e-foldings number, (II.14) By using Eq. (II.12) we have, where we have assumed that N ≪ H dS (t−t 0 ), or equivalently N ≪ N . Thus, the Hubble parameter may be expressed as follows, (II. 16) At the beginning of inflation we have N ≪ N and H ≃ H dS , while at the end of the early-time acceleration, when N = 0, one recovers H = H 0 . During the quasi de Sitter expansion of inflation the Hubble parameter slowly decreases. The slow-roll parameters are defined as follows, where we assumed that the constant-roll condition holds true. At the beginning of the early-time acceleration the first slow-roll parameter ǫ is small, in which case the slow-roll approximation regime is realized. For the solution (II. 16) in the limit N ≪ N , we get, On the other hand, for the β parameter we obtain a constant value, namely, This means that the model at hand satisfies the condition for constant-roll inflation. This fact has an important consequence on the form of the spectral index of primordial curvature perturbations, which will be independent from the total number of the e-foldings during inflation. The inflationary paradigm offers two important predictions about the inhomogeneities of our Universe at a galactic scale. Particularly, the perturbations around the FRW metric lead to a non-flat spectral index n s and also lead to a non-zero scalar-to-tensor ratio r.
In the case of F (R)-gravity, the inflationary indices have the following form, (II.20) By calculating these, we obtain, We can see that in the computation of the spectral index n s we can omit the contribution of ǫ which tends to vanish for N ≪ N (like in the Einstein frame case of the Starobinsky inflation). Since the constant-roll inflationary condition is assumed, it turns out that this index is in fact independent on the total e-foldings number.
The latest Planck data [34] constrain the spectral index and the scalar-to-tensor ratio as follows, n s = 0.9644 ± 0.0049 , r < 0.10 . (II.22) As a consequence, we must require N ≃ 60 in order to obtain a viable inflationary scenario. This means that at the beginning of inflation we have 60 ≪ N , a condition which solves the problem of initial conditions of the Friedmann Universe model we study. We can compare this result with the one corresponding to the ordinary Starobinsky model with γ 1 = 0. In this case, the Hilbert-Einstein term of (II.2) evaluated on the background solution (specifically, 6H 2 dS /κ 2 ) can be inserted in the perturbed Friedmann equation and we obtain instead of Eq. (II.9), the following equation, where, as we have discussed above, we have to consider 1/κ 2 ≪ H 2 dS . This equation implies the following, where t i is appropriately chosen to be the time instance at the beginning of inflation. Note that now we can eliminate the non-linear term of the solution, such that δḦ(t) ≃ 0 and the slow-roll parameter β in (II.17) is approximately equal to zero. We also require that for t = t 0 , when inflation ends, the Hubble parameter vanishes, namely, (II. 26) For the e-foldings number (II.14) we get, As a consequence, the resulting Hubble rate solution has the following form, (II.28) When t = t i , the e-foldings number is given by N = 18γ 0 κ 2 H 2 dS , and it is easy to verify that H = H dS . The first slow-roll parameter ǫ for the Starobinsky inflation is equal to, and the spectral index and the scalar-to-tensor ratio are equal to, Therefore, we must require N ≃ 60 in order to have concordance with the Planck data and therefore we need to specify the boundary value of the Hubble parameter as follows H 2 dS ≃ 60/(18γ 0 κ 2 ). In the model at hand, the first slow-roll parameter ǫ appearing in Eq. (II.18) decreases in exponential way, a behavior which is different in comparison to Eq. (II.29). However, when N ≃ 60 we obtain the same spectral index with the Starobinsky inflation with N ≃ 60, but we need to point out that in our model, the total amount of inflation measured by the e-foldings at the beginning of inflation, can be larger than N ≃ 60. This fact contributes to easily suppress the value of the scalar-to-tensor ratio in Eq. (II.21), as it is strongly encouraged by the BICEP2/Keck-Array data [35].
By imposing N ≃ 60 in Eq. (II.13) we obtain, 31) and the expansion curvature rate during inflation is defined in this way. The characteristic curvature at the time of inflation is R dS ≃ 10 120 Λ, in which case one has R 0 ≃ 1.8 × 10 85 Λ and from Eq. (II.6) we must require γ 1 ≪ 0.005. Finally, the relation between γ 0 and γ 1 is fixed by Eq. (II.7) and we obtain, (II.32)

The Reheating Era
At the end of the inflationary era, the model enters in a reheating phase. We will still assume that the contribution of the term f DE (R) in Eq. (II.2) is negligible, while the perfect fluid matter content of the Universe (namely, the energy density term ρ m ) has been shifted away, during the severe inflationary accelerated expansion that the Universe experienced, and has to be thermalized by the reheating process. During this phase, we will assume that the radiation/matter fields are bosonic scalar fields, which are described by the Klein-Gordon Lagrangian. By assuming that R ≃ R 0 and by taking the limit (II.6), Eq. (II.2) reads, The reheating solution takes place in an oscillatory phase, when the following holds true, Thus, if one completely neglects the right hand side term of Eq. (II.33), we have H(t) ∝ cos 2 [ωt] with ω = 1/ 24γ 0 κ 2 . In detail, one has to match the solution (II.12) with the new one in the form H(t) = f (t) cos 2 [ωt], where f (t) is a function of the time such thatḟ (t) 2 /f (t) ≃ 0 and f (t)ḟ (t) ≃ 0 (slowly-damping approximation). The form of the reheating (oscillatory) solution is [36], where t r is the time when the reheating starts. For t < t r , the solution is given by (II.12) and therefore, the transition to the reheating phase takes place when |Ḣ 2 /(2H)| ≃ |3HḢ|, namely, with N ≪ N being the total e-foldings of inflation. By equating (II.12) with (II.35) and by imposing h 0 ≃ 1, we can specify the frequency ω of the reheating solution as follows, where we have used Eq. (II.7) and we have introduced the constant-roll parameter β of Eq. (II. 19).
For large values of (t − t r ) the Hubble parameter behaves as follows, Since < H >≃ 2/(3(t − t r )), we get a matter dominated cosmological evolution. By taking into account that R ≃ 6Ḣ, we obtain, The Lagrangian of a scalar bosonic field χ with mass m χ , which is non-minimally coupled with gravity, is given by, where ξ is a coupling constant. Thus, the field equation is derived as follows, As a consequence, the Fourier modes of the field χ k ≡ χ k (t) with momentum k for a FRW spacetime, satisfy the following differential equation,χ By introducing the conformal time dη = dt/a(t) and also by redefining u k = a(t)χ k we have, As a result, the oscillations of the Ricci scalar change the effective mass term m 2 eff and the number of massive particles χ k increases with time. We note that, even in the case of minimal coupling with gravity ξ = 0, the effective mass m 2 eff still depends on the Ricci scalar R and the reheating mechanism takes place.
The energy density of radiation-ultra-relativistic matter is related to the temperature as follows, where γ is some constant. Since the energy density is connected to the square average of the Ricci scalar [36], from Eq. (II.40) we have ρ ∝ ω 4 and therefore we obtain, where we used (II.38).
In the Starobinsky inflationary scenario, the equation (II.33) with solution (II.35) is still valid, but during the reheating era we have, or equivalently, with N being the total number of e-foldings during the inflationary era. In this case, the frequency ω in (II.35) results to be, Thus, if we compare the reheating temperature of our model T with the reheating temperature corresponding to the Starobinsky inflation T St , we find, Since the β constant-roll parameter in (II. 19) for N ≃ 60 is small, we see that the temperature of the Universe after inflation predicted by our model is smaller respect to the one of the standard R 2 -scenario. This result shows that a constant-roll inflationary era has important effects on the reheating phase, as was also pointed out in Ref. [28].
B. Model II of Inflation: Curvature-corrected Exponential F (R) gravity with Constant-roll Evolution

Constant-roll Evolution in F (R) Gravity
In most theories of inflation in the context of F (R) gravity, it is common to adopt the slow-roll approximation, in order to produce the right amount of inflation that may be compatible with the observational data. In this section we shall deviate from the standard slow-roll approach, and we shall assume that a constant-roll era occurs during the inflationary era. The constant-roll inflationary paradigm was firstly used in the context of scalar-tensor theories [8][9][10][11][12][13][14][15][16][17][18][19][20][21], see also [22][23][24][25] and was extended in the context of F (R) gravity in [26][27][28]. As it was shown in [26], the most natural generalization of the constant-roll condition in the Jordan frame is the following, where β is some real parameter, which can be either positive or negative. The condition (II.51) is the most natural generalization of the constant-roll condition used in scalar-tensor approaches, which is, since the condition (II.52) is nothing else but the second slow-roll index η, which in the most general case is equal to η ∼ −Ḧ 2HḢ . We shall assume that the theory is described by a vacuum F (R) gravity and also that the background is a flat FRW metric. Upon variation of the gravitational F (R) gravity action with respect to the metric, we get the following equations of motion, where F R stands for F R = ∂F ∂R and also the "dot" denotes differentiation with respect to t. The dynamics of inflation in the context of F (R) gravity are governed by four inflationary indices, ǫ i , i = 1, ...4, which are defined as follows [1,[37][38][39][40][41][42], with the function E being equal to, Also for the calculation of the scalar-to-tensor ratio r, the quantity Q s is needed, which is defined as follows, (II.57) The spectral index of primordial curvature perturbations n s , in the case thatǫ i ≃ 0, is equal to [37][38][39], with ν s being equal to, The above relation is quite general and holds true not only in the case that ǫ i ≪ 1, but also when ǫ i ∼ O (1). With regard to the scalar-to-tensor ratio, in the context of vacuum F (R) gravity theories, it is defined as follows, where the quantity Q s is given in Eq. (II.57) above, and for the specific case of a vacuum F (R) gravity, the scalar-totensor ratio is equal to, The constant-roll condition (II.51), affects the inflationary indices of inflation ǫ i , i = 1, ..., 4 appearing in Eq. (II.55), which can be written as follows [26], where F RR and F RRR stand for F RR = ∂ 2 F ∂R 2 and F RRR F RRR = ∂ 3 F ∂R 3 respectively. It is conceivable that the inflationary dynamics crucially depend on the functional form of the F (R) gravity. We shall study a curvature corrected exponential F (R) gravity, which as we show, perfectly describes the early-time acceleration era. Also as we show in a later section, the late-time acceleration era is also successfully described by this model.

Constant-roll Evolution with Curvature-corrected Exponential F (R) Gravity
The constant-roll condition utterly changes the viability of an F (R) gravity model, as it was shown in Ref. [26]. Particularly, it is possible that an F (R) gravity model is not compatible with the observational data, when this is studied in the context of the slow-roll condition, however if the constant-roll condition is assumed, the model may be compatible with the observations, for a large range of values of the free parameters of the model. In this section we shall present one model of this sort, which is a curvature-corrected exponential model, with the functional form of the F (R) gravity being [43,44], where Λ = 7.93m 2 ,γ = 1/1000, m = 1.57 × 10 −67 eV, b is an arbitrary parameter [43,44,56] and n is a positive real parameter. This model has quite appealing inflationary dynamics in the context of constant-roll inflation, as we now evince. Consider the first equation of Eq. (II. 53), which for the model (II.63) and by assuming that H 2 ≫Ḣ during the inflationary era, it can be approximated as follows, where we have set R 0 = bΛ and γ = Λ 1000 1 (3m) n for notational simplicity. It is conceivable that the dynamical evolution is not affected by the terms containing the exponentials, at least in the era where the approximation H 2 ≫Ḣ holds true. By using the constant-roll condition (II.51) and after some algebraic manipulations, Eq. (II.64) can be simplified as follows, (II.65) The differential equation (II.65) can be analytically solved, with the solution being, with C 1 being an arbitrary integration constant which plays no role in the dynamics of inflation, as we shall show. Also it is notable that the parameter γ appearing in Eq. (II.64), does not appear in the final differential equation that governs the dynamical evolution. Having the Hubble rate at hand, it is easy to calculate the slow-roll indices ǫ i , i = 1, .. and by using these, the parameter ν s defined in Eq. (II.59), is equal to, ν s = 1 2 ((β + 2) (6 n + 1260) n 2 − n (4β (6 n + 297) + 3 (6 n + 852)) + (3β + 1)6 n ) 2 (−432(β + 2)n 2 + n (432β − 6 n + 900) + 6 n ) 2 , (II.68) and therefore in this case, the spectral index of primordial curvature perturbations n s can be cast after some algebra in the following form, (II.69) Accordingly, by using Eq. (II.61), the scalar-to-tensor ratio is found to be, It is noteworthy that both the spectral index and the scalar-to-tensor ratio depend only on β or n. Having the final expressions for the spectral index and the scalar-to-tensor ratio at hand, we shall investigate the parameter space in order to see for which values of the free parameters, the compatibility with the observational data can be achieved.
The latest Planck data [34] constrain the spectral index and the scalar-to-tensor ratio as in Eq. (II.22), so now we investigate which values of the parameters (n, β) may render the curvature-corrected constant-roll model of Eq. (II.63) compatible with the Planck data. A detailed analysis reveals that there is a large range of parameter values that may render the model compatible with the observations. For example by choosing (n, β) = (2.1, −8.7), the spectral index becomes n s = 0.966239 and the corresponding scalar-to-tensor ratio becomes r = 0.0119893. Also for (n, β) = (0.9, −1.08), the spectral index becomes n s = 0.96742 and the corresponding scalar-to-tensor ratio becomes r = 0.0936944. Finally for (n, β) = (1.5, −0.4), the spectral index becomes n s = 0.960444 and the corresponding scalar-to-tensor ratio becomes r = 0.0669277. In Fig. 1 we plotted the functional dependence of the spectral index as a function of β, for n = 1.5 (left plot), and for n = 2.1 (right plot). In both cases, In the red and black straight lines correspond to the Planck data allowed values for the spectral index, namely, n s = 0.9693 and n s = 0.9595 respectively. The same applies for the scalar-to-tensor ratio, and in Fig. 2, we plotted the functional dependence of the scalar-to-tensor ratio as a function of n. In the left plot, the parameter n is assumed to be n = 1.5, while in the right plot n = 3. Concluding, we demonstrated that the curvature-corrected F (R) gravity model of Eq. (II.63) in the context of the constant-roll condition, provides an inflationary era which is compatible with the latest Planck data. An intriguing feature of the model (II.63), is that the parameter n plays a crucial role in the late-time era. As we show in a later section, certain values of n may produce an undesirable amount of dark energy oscillations, so the values of n we used into this section, will be further constrained by the dark energy oscillations study. Before we continue, it is worth discussing the resulting scalar-to-tensor ratio as it appears in Eq. (II.70). The resulting scalar-to-tensor ratio is quite an important result, that might be tested in the next generation of CMB experiments. So far, there have been several ground based CMB projects under design or already on going, which includes the CMB-S4 [46], BICEP3 [47], and AliCPT [48], and so on. It is interesting to note that, the operation of all these ground based experiments are able to cover the full CMB sky without the help of satellites such as WMAP or Planck, and correspondingly, the observational limit on the tensor-to-scalar ratio could be improved to about two orders of magnitude higher in comparison to the current one. This means that the parameter space of the models under consideration would be much constrained in the near future.
Before we close this section, it is worth discussing in brief the possibility of a graceful exit from the inflationary era, for the model (II.63). In general, the graceful exit issue can be a cumbersome task and also a conceptual challenge, especially in the context of F (R) gravity. In the literature, the graceful exit from inflation instance is identified with the time instance that the slow-roll approximation ceases to hold true, which occurs when the first slow-roll index becomes of the order O(1). However, in Ref. [49] an alternative viewpoint was provided for the graceful exit from inflation issue. Particularly, the graceful exit from inflation could be triggered by the occurrence of growing curvature perturbations at some point during the inflationary era, which may disturb the inflationary attractor of the theory. In effect, the dynamical attractor becomes unstable, and hence ceases to be the final attractor of the theory, and the cosmological dynamical system is no longer described by this attractor, and therefore the dynamical cosmological evolution of inflation is interrupted. We believe that this way of thinking is conceptually more accurate for the description of the graceful exit from inflation.
So let us now investigate if the inflationary solution (II.66) ceases to be the final attractor of the cosmological system, and in order to see this, we linearly perturb this solution, as follows, where we identify H 0 (t) with the solution (II.66). By inserting Eq. (II.71) in the first equation of Eq. (II.53), and after some tedious algebraic manipulations, we obtain the following differential equation at leading order, which can be solved and the solution is, A simple analysis can easily reveal the behavior of the perturbation ∆H(t), for the various values of the parameter β, and in fact it can be shown that for β < 0, the perturbations decay as t −1 , while for β > 0, the perturbations grow as the time evolves to larger values. Thus, when β > 0, the cosmological solution H 0 (t) of Eq. (II.66) is unstable towards linear perturbations, and therefore this can be viewed as a strong indication that the graceful exit from inflation actually happens in this case. However, the full study of this issue deserves an article focusing on this, since there are various questions that remain unanswered, for example, how many e-foldings occur before the end of inflation, what happens when non-linear terms are taken into account and so on. Some of these tasks may prove quite difficult to address, since the lack of analyticity makes the problem quite cumbersome. We aim though in a future work, to address the aforementioned tasks by using an autonomous dynamical system approach.

III. LATE-TIME ACCELERATION ERA
In this section we will investigate the behavior of the model I appearing in Eq. (II.1) during the late-time era. In order to do so, we need to introduce the dark energy part in the action, in terms of the function f DE (R) appearing in Eq. (II.1). We will use a modified version of exponential gravity, namely, where b is a positive parameter and Λ is the cosmological constant. The function of the Ricci scalar g(R) is necessary to stabilize the theory at large redshifts and we shall assume that it has the following form, where c is a real and positive parameter. We need to note that, the existence of a quintom scenario where the equation of state (EoS) parameter evolves crossing the phantom divide line, is strictly connected with the possibility of having a stable de Sitter epoch at late times. This scenario is supported by cosmological observations, which allow an oscillatory behavior of the dark energy EoS parameter around the line of the phantom divide, although the data are still far from being conclusive. In this respect, an exhaustive review about the quintom scenario can be found in Ref. [50], where several successful examples of quintom cosmological models and their corresponding observational consequences were analyzed. As a general feature of the model, we immediately see that, at R = 0, one has f DE (R) = 0 and we recover the Minkowski spacetime solution of Special Relativity. When 4Λ ≤ R, f DE (R) ≃ −2Λ/κ 2 we obtain the standard evolution of the ΛCDM model. Moreover, since |f DE (R)| ∼ 10 −120 M 4 P l , we have that the modification of gravity for the dark energy sector is completely negligible in the high curvature limit of the inflationary era, where R/κ 2 ∼ M 4 P l .
We will now verify that with the modification we introduced in Eq. (III.1), it is possible to pass the cosmological and local tests. In order to do so, we introduce the following form of the F (R) gravity, in which case the trace of the field equations results to the following equation, where T is the trace of the stress-energy tensor of the matter-radiation perfect fluids that quantify the matter content of the Universe. When g(R) ≃ 1, it is easy to see that the following conditions hold true, The first condition is necessary in order to obtain the correct value of the Newton constant and avoid anti-gravitational effects during the matter, radiation and dark energy eras, while the second condition guarantees the stability of the model with respect to the matter perturbations [6,51,52]. When T = 0, which corresponds to the vacuum F (R) gravity case, the non-trivial solution of Eq. (III.4) is given by the de Sitter spacetime with R = 4Λ. This is a consequence of the fact that g(R) = 1 and f DE (R) ≃ −2Λ when R = 4Λ. The perturbations around the de Sitter solution are governed by the equation, For our model, the effective mass m 2 eff turns out to be positive, rendering the solution stable. This means that the de Sitter spacetime is a final attractor of the cosmological system for large values of the cosmic time, when in the expanding Universe the contents of matter and radiation vanish.
During the matter and radiation domination eras, the model we used mimics an effective cosmological constant, if the function g(R) in Eq. (III.2) is close to unity, namely where recall that R 0 is the curvature of the Universe at the end of the inflationary era. For example, if c = 10 −5 , we obtain f DE ≃ 2Λ/κ 2 up to the value R ≃ 4Λ × 10 4 . For larger values of the curvature, matter and radiation dominate strongly the evolution. Moreover, the conditions (III.5) read, The second condition always holds true for positive values of c, but on the other hand, the first condition requires a careful investigation in order to avoid anti-gravitational effects in the large curvature limit. For example, assume that R 0 ≃ 10 86 Λ, then if c ≃ 10 −5 , we obtain (F R (R) − 1) ≃ 10 −3 and it is obvious that this condition is satisfied.

A. The Radiation, Matter Domination Eras and Transition to Late-time Acceleration Era
In order to investigate the behavior of our model during radiation and matter domination eras, but also during the transition to the late-time era, and until present-time, we need to introduce the following variable, which is known as the "scaled dark energy" [1,53]. This variable encompasses the ratio between the effective dark energy and the standard matter density, evaluated at the present time, with the matter density defined as follows, where m is the mass scale associated with the Planck mass. In the expression (III.9), the variable z = [1/a(t) − 1] denotes the redshift as usual, and also χ stands for χ ≡ ρ r(0) /ρ m(0) .
If one extends the expression appearing in Eq. (III.3) as follows, in order to include the entire form of the gravitational Lagrangian (II.1), it is possible to derive from the FRW field equations the following equation [1,53], where the functions J i , i = 1, 2, 3 stand for, . (III.13) The Ricci scalar is easily derived from Eq. (III.9) and it is equal to, (III.14) At the late time regime, where z ≪ 1, we can avoid the contribution of the matter and radiation fluids, in which case, the solution of Eq. (III.12) reads, with y 0 being an integration constant. Since for the exponential gravity ΛF RR (4Λ) ≪ 1, the argument of the square root is positive, in effect, dark energy oscillates around the phantom divide line w = −1. The frequency of the oscillation with respect to log[z + 1] is given by, Generally speaking, since ΛF RR (4Λ) ≃ 2b 2 exp[−4b], the oscillation frequency at late times does not diverge. Such a problem may emerge in the high curvature regime, that is for redshifts that satisfy 0 ≪ z, when the contribution of the dark energy in Eq. (III.14) is negligible. In this case, after the expanding Eqs. (III.12)-(III.13) with respect to y H /(z + 1) 3 ≪ 1, we find the following solution in the vicinity of a given redshift z [43,44], where y 0 is an integration constant. Here, the oscillation frequency of the dark energy is equal to, As a consequence, when RF RR (R) is very close to zero, like in the case of pure (vacuum) exponential gravity, the effective dark energy induced by modified gravity oscillates with high frequency and diverges, rendering the theory unstable. However in our model, due to the presence of the function g(R) chosen as in Eq. (III.2), one has, ν ≃ 2/c 2π(z + 1) .
(III. 19) This means that, back into the past, during the radiation and matter domination eras, the frequency of the effective dark energy oscillations, tend to decrease and the theory is protected against singularities.

B. Dark Energy Era for the Curvature-corrected Exponential Model
In this section we shall analyze certain features of the late-time evolution corresponding to the curvature-corrected exponential F (R) gravity model of Eq. (II.63), having to do with dark energy oscillations. It is known that certain viable models of dark energy produce dark energy density oscillations during the matter domination era, and the frequencies of these oscillations may diverge [43,44,54,55]. In the modified gravity theories and especially in the F (R) gravity theories, fourth order derivatives appear in the theory, and hence high-frequency oscillations occur and the background may oscillate in a rapid way. In effect, perturbation theory may break down, since non-linearities may occur. It is therefore vital for F (R) gravity theories to overcome this theoretical obstacle. In this section we shall discuss this problem and by using some numerical analysis for the curvature corrected exponential model, and we shall investigate in which cases the dark energy oscillations may occur.
A convenient quantity that can easily quantify the dark energy oscillations is the scaled dark energy y H (z) = ρDE ρ (0) m appearing in Eq. (III.12). Actually, by solving the differential equation (III.12) numerically we will be able to find how the dark energy oscillations behave at late-times, so for small redshifts. For our numerical analysis we shall use the following initial conditions [43,44,[54][55][56], with z f being z f = 10 and also Λ is equal to Λ ≃ 11.89eV 2 . In Fig. 3 we have plotted the scaled dark energy y H (z) = ρDE behaves in a very similar way. Having the scaled dark energy at hand, we can easily investigate how the dark energy oscillations behave by studying the dark energy equation of state parameter ω DE = P DE /ρ DE , which in terms of y H (z) is defined as follows, In Fig. 4 we present the behavior of the dark energy equation of state parameter for n = 1.5 (blue curve), n = 2.1 (red curve) and n = 0.9 (black curve). As it can be seen in Fig. 4, the oscillations have smaller amplitudes for n = 2.1, and these start during the matter domination era. In all cases, the dark energy oscillations stop to occur at small redshifts. It is noteworthy that the values of n around n ∼ 2, always provide the smaller amplitudes for the oscillations.
In order to have a clear picture for the physical behavior of the cosmological system, we shall investigate the behavior of the total equation of state parameter ω eff (z), as a function of the redshift z, which in terms of y H (z) is defined as follows, (III. 22) In Fig. (5) we plot the behavior of the total equation of state parameter as a function of the redshift z, for n = 1.5 (blue curve), n = 2.1 (red curve) and n = 0.9 (black curve). As it can be seen, the case n = 2.1 has the mildest oscillatory behavior. It is noteworthy that in all cases, for small redshifts the phantom divide line w = −1 is not crossed. Finally, in Fig. 6 we plot the Hubble rate as a function of the redshift for n = 1.5 (blue curve), n = 2.1 (red curve) and n = 0.9 (black curve). Note that in terms of the scaled dark energy parameter y H (z), the Hubble rate is defined as follows, H(z) = m 2 y H (z) + g(a(z)) + χ(z + 1) 4 , (III.23) The analysis of the dark energy equation of state parameter and also of the total equation of state parameter, indicates strongly that the exponential model with nearly R 2 curvature corrections, has quite appealing properties, since the dark energy oscillations are less pronounced. In this way, we obtain an optimal reproduction of the ΛCDM model, and the effects of dark energy remain negligible during the early and mid stages of the matter and radiation eras. Now we need to fix the boundary conditions of our cosmological dynamical system at large redshift z = z max . They can be inferred from the form of ρ DE in Eq. (III.9) for the case of F (R)-modified gravity, namely, where R ≡ R(z) and H ≡ H(z) are functions of the redshift. At large redshift, during the matter era, we have to take R = 3m 2 (z + 1) 3 and H = m(z + 1) 3/2 and the boundary conditions of the system are given by, where, For z max = 10, in which case χ(z max + 1) ≃ 0.00341 ≪ 1, and we effectively are in a matter dominated Universe, we obtain, These values can be compared with the corresponding ones for the ΛCDM model, where y H is a constant, namely y H = Λ/(3m 2 ) = 2.17857. We argue that our model is extremely close to the ΛCDM model at very high redshift.
Here we recall that the first observed galaxies correspond to a redshift z ≃ 6. Finally, the contributions of matter and radiation are determined by the values of m 2 and χ in (III.9). The cosmological data indicate that, m 2 ≃ 1.82 × 10 −67 eV 2 , χ ≃ 3.1 × 10 −4 . (III.32) By solving numerically the differential equation for y H (z) in the redshift range, −1 < z < z max we obtained the results which appear in Fig. (7) and Fig. (8), where we plotted the behavior of y H (z) and ω DE (z) as a function of the redshift, for −1 < z < 10. Despite of the fact that at high redshifts, the amplitude of the oscillations of the effective EoS parameter around the phantom divide line gradually grows, we see that their frequency decreases and thus, singularities are avoided. In order to measure the matter energy density ρ m (z) at a given redshift, we introduce the parameter y m (z) as In Fig. (9) the plot of y H (z) is compared with the graphic of y m (z) for −1 < z < 1. We see that y H (z) is nearly constant and it is dominant over y m (z), for z < 0.4, a feature that is in full agreement with the ΛCDM description.
The Ω DE (z) parameter, Ω DE (z) ≡ ρ DE ρ eff = y H (z) y H (z) + (z + 1) 3 + χ(z + 1) 4 , (III.34) is frequently used to express the ratio between the dark energy density ρ DE and the effective energy density ρ eff of The latest cosmological data [34] indicate that, Ω DE (z = 0) = 0.685 ± 0.013 and ω DE (z = 0) = −1.006 ± 0.045. Thus, our model fits the observational data at present time. Hence, we have shown that the model I can mimic the late-time features of the ΛCDM Model. The additional degree of freedom introduced by the F (R)-gravity, leads to an oscillatory behavior during the last stages of the matter domination era and the early stages of the dark energy era. However, the theory is stable and protected against singularities. We need to stress that in the simulation presented in this section, we have considered the whole form of the gravitational Lagrangian in (II.1). Thus, we confirm that the higher curvature corrections for inflation are negligible in the limit of small curvatures.

Conclusions
In this paper we demonstrated that it is possible to have a unified description of a constant-roll inflationary era with the dark energy era, by using two F (R) gravity models. Particularly, we used a R 2 -corrected logarithmic F (R) gravity model and a curvature-corrected exponential F (R) gravity model. With regard to the inflationary era, we demonstrated that it is possible to obtain observational indices compatible with the latest Planck data, in the context of constant-roll inflation. Particularly, the constant-roll condition broadens the parameter space, and this feature makes it more easy to have compatibility with the observational constraints.
In addition, we studied the late-time evolution of the F (R) gravity models. A noticeable feature related with the exponential model is that the dark energy oscillations have quite small amplitude when the curvature corrections are nearly of the R 2 form.
It would be interesting to analyze the possibility of having a unified description of constant-roll inflation with a dark energy era, in the context of other modified gravities, like for example Gauss-Bonnet or F (T ) gravities, where it is possible to unify inflation with dark energy, see Ref. [1] for a comprehensive review. This issue deserves a study which we defer in a future work.
With regard to the graceful exit from inflation issue, which we briefly addressed for one of the inflationary models we presented in this paper, it should be noted that it deserves a more focused analysis. With the approach we adopted in this work, we provided hints that the graceful exit may occur as a result of growing curvature perturbations, and with exit from inflation we mean that an inflationary attractor solution may become unstable as the time evolves. In a more focused work, one should in principle construct an autonomous system of differential equations, and prove numerically that the inflationary attractors are indeed unstable stationary points of the dynamical system. Work is in progress towards the aforementioned line of research.
Finally, when gravitational systems different from the Einstein gravity are considered, we can find new cosmological solutions with a large variety of new features. In this work, we have seen that a constant roll inflation can be obtained by introducing higher derivative corrections in the action, while a quintom scenario for the dark energy can be obtained at late times, in the context of F (R)-gravity. A modified theory of gravity may also allow for an alternative description with respect to the Big Bang theory as in the bounce scenario, where a cosmological contraction is followed by an expansion at a finite time. The idea that instead from an initial singularity the universe has emerged from a cosmological bounce has been largely analyzed in the literature (see for example Refs. [57][58][59][60][61][62] or Ref. [63] for bounce cosmology description in the context of teleparallel gravity).