Is exponential gravity a viable description for the whole cosmological history?

Here we analyse a particular type of F(R) gravity, the so-called exponential gravity which includes an exponential function of the Ricci scalar in the action. Such a term represents a correction to the usual Hilbert–Einstein action. By using Supernovae Ia, Barionic Acoustic Oscillations, Cosmic Microwave Background and H(z) data, the free parameters of the model are well constrained. The results show that such corrections to General Relativity become important at cosmological scales and at late times, providing an alternative to the dark energy problem. In addition, the fits do not determine any significant difference statistically with respect to the Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}CDM model. Finally, such model is extended to include the inflationary epoch in the same gravitational Lagrangian. As shown in the paper, the additional terms can reproduce the inflationary epoch and satisfy the constraints from Planck data.


Introduction
Over the last decade, the study of some modifications of General Relativity have drawn a lot of attention, particularly in the framework of cosmology as an attempt to provide a more natural explanation to the accelerating expansion at early times (inflation) and at late times (dark energy epoch). In this sense, the most simple and natural extension of GR arises as the generalisation of the Hilbert-Einstein action by assuming a non-linear function of the Ricci scalar, which is commonly called f (R) gravity (for a review see [1][2][3][4][5][6][7][8]). Other extensions include curvature invariants such as the Gauss-Bonnet a e-mail: odintsov@ice.csic.es b e-mail: saez@ice.csic.es c e-mail: sharov.gs@tversu.ru gravity [9][10][11][12][13] or generalisations of the so-called Teleparallel gravity, an equivalent theory to GR constructed as a gauge theory of the translation group leading to a null-curvature theory with non-null torsion (see Ref. [14][15][16][17]). Nevertheless, f (R) gravities have been by far the most analysed extension of GR over the last years, also due to the motivation as regards more fundamental theories as string theory [18]. This extensive study has provided a very deep knowledge and comprehension of this type of theories, whose field equations turn out fourth order differential equations instead of second order as in GR. Nevertheless, f (R) gravities can easily be reduced to a type of scalar-tensor theory, i.e. f (R) gravity basically implies the appearance of an extra scalar mode [19][20][21]. As every theory with extra propagating modes, this may imply the existence of ghosts. Fortunately, this is not the case in f (R) gravities. However, the extra scalar mode may imply violations and deformations of well-known and tested predictions of GR. In order to avoid large corrections at scales where GR is very well tested, f (R) gravities can hide such extra mode through a mechanism known as chameleon mechanism, proposed initially in the framework of scalartensor theories [22,23], but rapidly extended to f (R) gravities [24,25].
In addition, the versatility of f (R) gravities allows one to reconstruct any cosmological solution with the suitable evolution [26][27][28][29][30][31][32][33][34][35]. Then late-time acceleration may arise in a natural way as a consequence of the gravitational theory instead of being the aftermath of any extra unknown field. Moreover, simultaneously f (R) gravities may contribute to the compensation of the large value predicted by quantum field theories for the vacuum energy density, and particularly may play an essential role in the framework of the so-called unimodular gravity theories. In this regard, f (R) gravity scenarios as an alternative to the CDM-cosmology are inter-esting and attractive, since they are able to describe simultaneously the early-time inflation and the late-time acceleration in the expansion of our universe [19][20][21]25,[36][37][38][39]. Particularly, some of the most promising inflationary models are constructed within the f (R) gravity scenario, since some of these models can easily reproduce slow-roll inflation by mimicking a cosmological constant at early times and then decaying, leading to a power spectrum for scalar perturbations nearly invariant and a negligible scalar-to-tensor ratio, coinciding with the last data released by the Planck collaboration [40,41]. This is the case for instance of Starobinsky inflation [42], a quadratic Lagrangian on the Ricci scalar that predicts the correct values for the spectral index and the scalarto-tensor ratio. Actually, some analysis suggest that any deviation from Starobinsky inflation should be small enough to avoid deviations from its well-established predictions [43]. Keeping this in mind, over the last years some efforts have been focussed on the attempt to unify inflation and dark energy epoch in the framework of f (R) gravities, and particularly within the so-called viable f (R) gravity models [25]. As mentioned above, these viable models pass the wellknown local tests, where the scalar mode acquires a large mass through the chameleon mechanism avoiding large corrections with respect to GR. Hence, the local tests or the Solar System tests for viable f (R) theories include correct Newtonian and post-Newtonian limit [24,25]. In addition, this type of models are capable of reproducing the correct late-time acceleration, in general by simulating an effective cosmological constant that becomes important at late times, while CDM behaviour is recovered at high redshift. Moreover, these models provide good fits when compared with observational data, being almost indistinguishable from CDM [44]. However, viable f (R) gravities contain a type of future cosmological singularity, the so-called sudden singularity, a consequence directly related to the mass of the scalar field that avoids corrections at local scales [45], although such singularity occurs in the future when the right parameters are set and can be avoided by adding some extra terms. Moreover, some extensions of such models are also capable of reproducing inflation at early times, when tends asymptotically to a power Lagrangian, leading to a Starobinsky-like inflation keeping the right predictions [25].
In this paper, we focus on the analysis of a type of viable f (R) models that reproduces late-time acceleration by mimicking a cosmological constant but where corrections may have some distinguishable effects. This class of f (R) models are given by a negative exponential of the Ricci scalar in the action, which turns out negligible at large redshifts but becomes important at late times, an effect easily controlled with a free parameter related to current Hubble parameter. Exponential gravity has been previously analysed in Refs. [39,[46][47][48] as a reliable alternative to other viable f (R) gravities, since GR results are recovered at local scales but reproduce dark energy behaviour at cosmological ones. In addition, the previous analysis has shown the existence of an asymptotically stable de Sitter solution in such exponential Lagrangians, leading to an approximated CDM behaviour at the present time [47,48]. Moreover, some recent analysis of such type of exponential gravities suggest that observational constraints can be well satisfied from the cosmological point of view, in such a way that f (R) gravity and CDM model turn out nearly indistinguishable, as suggested by previous analysis [49][50][51]. In addition, exponential gravity can be extended to cover the inflationary stage as well. To do so, an additional exponential is considered in the gravitational action becoming important at large curvature when the inflationary period occurs, and turning out to be negligible as curvature decreases [39,47]. Hence, in this paper we analyse such a type of f (R) gravities, firstly by fitting the free parameters of the model by using data from Type Ia supernovae, baryon acoustic oscillations (BAO), estimations of the Hubble parameter H (z) and parameters of the cosmic microwave background radiation (CMB) [52][53][54][55], and also considering different approaches. Then we analyse how the full gravitational Lagrangian can cover also the inflationary epoch, obtaining the spectral index for scalar perturbations and the tensor-to-scalar ratio.
The paper is organised as follows: Sect. 2 reviews the basics of f (R) gravities, while Sect. 3 is devoted to the introduction of the exponential f (R) gravity model and its dynamical equations. In Sect. 4 the observational data considered in the paper is shown, this includes Union 2.1 observations of Type Ia supernovae, BAO effects, the latest measurements of the Hubble parameter H (z) and CMB parameters. In Sect. 5 we estimate the constraints on the exponential F(R) model from the aforementioned data. In Sect. 6 we investigate the variant of the exponential model with inflation terms in the Lagrangian. Finally Sect. 7 is devoted to the conclusions of the paper.

F(R) gravity
Modified F(R) gravities are described by the following generalisation of the Einstein-Hilbert action [26][27][28][29][30][31][32][33][34][35]: where κ 2 = 8π G and S m is the matter action. Einstein General Relativity is very well understood and tested at many scales, so that one should assume the action (2.1) to contain slightly deviations from GR, such that we can rewrite the action in the following way: Here, the function f (R) accounts for the gravitational modifications and should become negligible at scales where GR is well tested. By varying the action (2.1) with respect to the metric tensor g μν , the field equations are obtained: where R and R μν are the Ricci scalar and Ricci tensor, respectively, whereas F R ≡ F (R) and T μν is the energymomentum tensor of matter. Note that the F(R) field equations are fourth order to be compared to the second order of General Relativity. However, the action (2.1) hides an additional scalar mode, such that it can be expressed by the Lagrangian of a type of scalar-tensor theory as follows: where the following relations are found: Hence, in order to avoid large deviations from GR, this additional degree of freedom should be hidden at the appropriate scale, a mechanism commonly known as the chameleon mechanism [22,23]. In this sense, some F(R) actions which accomplish this requirement have been proposed in the literature [24,25], particularly some of them with the form of a negative exponential, the type of Lagrangians we are exploring in this manuscript. Nevertheless, let us first analyse the general properties of F(R) gravities, and in particular in the cosmology framework. By assuming a spatially flat Friedman-Lemaître-Robertson-Walker (FLRW) spacetime with the metric where a(t) is the scale factor of the universe, c = 1, and the Ricci scalar is expressed as Here, the Hubble parameter is defined as usual by H =ȧ/a, where the dot denotes derivatives with respect to the cosmic time. By assuming an energy-momentum tensor T μ ν = diag (−ρ, p, p, p) as a perfect fluid, where ρ and p are the matter energy density and pressure, the field equations (2.3) turn out to be [26][27][28][29][30][31][32][33][34][35] While the divergence of the field equations lead to the energy conservation equation ∇ μ T μν = 0, which in a FLRW metric becomeṡ The FLRW equations (2.7) can be expressed in terms of other independent variables instead of the cosmic time for convenience. Here we use the number of e-folds, given by x = log a = − log(z + 1) with a(t 0 ) = 1 at the present time t 0 , to express the above equations in the form of a dynamical system as follows: Here F R R ≡ F (R) are derivatives with respect to x and we have used the Ricci scalar definition Eq. (2.6) and the continuity equation (2.8). Hence, the analysis of the above system can provide all the information as regards the dynamics produced by a particular action F(R).

Exponential gravity
Let us now introduce the type of exponential F(R) gravity, we are considering in this manuscript [39,[46][47][48][49][50][51] The model contains just two free parameters and R 0 , which may be expressed in a more convenient way as R 0 = 2 /β, where β is dimensionless [49][50][51], Note that in principle the model (3.1) can well describe the universe evolution for z < 10 4 , including the recombination epoch, the matter-dominated era and the late-time acceleration. This is true as far as β ≥ 0, since the exponential becomes negligible and the action (3.1) recovers the usual CDM model at large redshifts, where the curvature becomes much larger than . In Sect. 6 we will also consider corrections to this model such that early-time inflation is also described.
Here we focus on the epoch for 0 ≤ z ≤ 10 4 , when the content of the universe includes pressureless (nonrelativistic) matter and radiation (relativistic particles): ρ = ρ m +ρ r , such that the continuity equation (2.8) can be solved and yields where ρ 0 m and ρ 0 r are the present time values of these components, which can be normalised over the critical density as follows: Here H 0 is the Hubble parameter today. Let us first explore the behaviour of model (3.1) during the early universe, when the curvature becomes large R → ∞ as z → ∞ (or for z ≥ 10 4 in practical applications). Then the model (3.1) transforms into the CDM model with F(R) = R − 2 , so the solutions of the system (2.9) tend asymptotically to CDM at large redshifts, leading to Here the index " CDM" refers to quantities calculated within the CDM model, where CDM = 3(H CDM 0 ) 2 and H CDM 0 the Hubble parameter today as predicted by the CDM model, while X = r / m . However, despite the model (3.1) recovering CDM at large redshifts, late-time evolution deviates from CDM, such that the above quantities as measured today, t = t 0 , would differ from CDM unless initial conditions are fixed at z = 0, which is not the case of our paper. Note that other viable f (R) models shows a similar behaviour when they are analysed asymptotically [24]. Hence, we have in general where we have denoted by 0 those magnitudes measured today as predicted by our model (3.1). Nevertheless, we can connect the two models through the relation of the physical matter density [24] 0 As will be shown below, this remark is important when performing the fitting analysis for the observable parameters in Sect. 4. Moreover, note that the first FLRW Eq.(2.7) for the CDM model is a constraint equation which evaluated at t = t 0 can be expressed as follows: This expression is very well known in standard cosmology when GR is assumed but breaks down when other gravitational actions beyond GR are considered, like F(R) gravity.
In such a case, the first FLRW equation becomes a dynamical equation, since it contains second derivatives of the Hubble parameter. By evaluating the FLRW equation in (2.7) at z = 0, the above equation can be expressed as Note that here we have defined 0 = 3H 2 0 , which refers to the cosmological constant term in the action (3.1), while 0 f (R) includes the exponential function in (3.1). The smaller 0 f (R) is, the closer our model is to CDM at the present time, where Eq. (3.7) is evaluated. Nevertheless, note that our model recovers CDM asymptotically at high redshifts (z > 10) such that the differences among the relative densities i (z) become negligible in both models at high redshifts. Let us now, for convenience in the calculations, introduce the following dimensionless variables: Hence, the gravitational action (3.1) becomes while the system of equations (2.9) takes the form Recall that the variable x = log a = − log(z + 1) refers to the number of e-folds. This system can be solved numerically by setting the appropriate initial conditions. As is natural for the model (3.1), we assume initial conditions that match the CDM model at a particular high redshift: which corresponds to the CDM asymptotic solution (3.4) at an initial redshift z i , or alternatively at x i . The value of x i is determined by assuming the following condition: where ε is a small number in the range 10 −10 < ε < 10 −7 , such that our model mimics the CDM solution (3.4) at x < x i , and the corresponding solutions practically do not depend on ε or x i for all x.
Alternatively, we can also use the following variable [24,47,[49][50][51]: Then the equation for H in (2.9) can be rewritten as follows: (3.14) where the second expression corresponds to the initial condition (3.11). The advantage of the function y H and its derivative (3.14) lies on their finiteness at z → ∞ (x → −∞). However, the same does not apply for Eq. (3.10), since is not finite for every redshift. Hence, numerical integration of the system Eq. (3.14) or the corresponding second order differential equation for y H shows similar difficulties to the system (3.9) and (3.10).

Observational data
Let us now present the data we are using here to fit the free parameters of our model. Besides the late-time evolution data from SNe Ia, BAO and H (z), we are also considering the CMB parameters. Then, to do so, we have to include radiation in our equations, or in other words, assuming Eq. (3.11), or alternatively Eq. (3.14). As our model mimics CDM at high redshifts, we can reduce the number of free parameters by fixing the radiation-matter ratio as provided by Planck [55]: can be considered as a nuisance parameter, so that can be marginalized for all fits.

Supernovae Ia data
The Union 2.1 compilation provides [52] N SN = 580 Sne Ia with their observed (estimated) distance moduli μ i = μ obs i for redshifts z i in the interval 0 ≤ z i ≤ 1.41. In order to fit the free parameters of our model, we compare μ obs i with the theoretical value μ th (z i ), where the distance moduli is given by Here D L (z) is the luminosity distance. The corresponding χ 2 function is calculated by computing the differences between the SNe Ia observational data and the predictions of a particular model with parameters p 1 , p 2 , . . . , [52]. The marginalisation over the nuisance parameter H 0 is widely described in the literature (see Refs. [81][82][83]).

BAO data
Baryon acoustic oscillations (BAO) are obtained from galaxy clustering analysis and include measurements of two cosmological parameters [53] where r s (z d ) is the sound horizon at the decoupling epoch and D V (z) is given by The values (4.5) were estimated for redshifts z = z i (and redshift ranges) of galaxies from a peak in the correlation function of the galaxy distribution at the comoving sound horizon scale r s (z d ), which corresponds to the decoupling of the photons z d . In this paper we use the BAO data from Refs. [56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74] for the parameters (4.5), which provides N BAO = 17 data points for d z (z) and 7 data points for A(z), both shown in Table 1. We use the covariance matrices C d and C A for correlated data from Refs. [66,69] described in detail in Ref. [81]. So the χ 2 function for the values (4.5) yields where d and A are vector columns with As pointed out above, the Hubble parameter today H 0 as predicted by our model All these considerations have to be carefully studied in order to choose the appropriate approach to calculating the sound horizon r s (z d ) from different fitting formulae [50,74,84,85]. Here we are considering the following simple fitting formula [81]:  (2) making H (z) estimations from differential ages t of galaxies [75][76][77][78][79][80] via the following relation: To avoid additional correlation with the BAO data from Table 1, we use in this paper only N H = 30 values H (z) estimated from differential ages of galaxies, shown in Table 2.

CMB data
Unlike the SNe Ia, BAO and H (z) data described above, corresponding to the late-time era 0 < z ≤ 2.36, cosmological observations associated with CMB radiation [74,84,85] include parameters at the photon-decoupling epoch z * 1090 (z * = 1089.90 ± 0.30 [55]), particularly the comoving sound horizon r s (z * ) and the transverse comoving distance, In the present manuscript, we use the CMB parameters in the following form [84,85]: Here 0 b is the present time baryon fraction. The distance priors (4.11) with their errors σ i and the covariance matrix were derived in Ref. [85] from the Planck collaboration data [55] with free amplitude of the lensing power spectrum. For the value z * we use the fitting formula from Refs. [84,85,87]; the sound horizon r s (z * ) is estimated from Eq. (4.9) as the correction r s = dr s dz z. Hence, the χ 2 function corresponding to the data (4.10)-(4.11) is obtained as follows: which is minimised by marginalizing over the additional parameter ω b = 0 b h 2 , which should be considered as a nuisance parameter, as well as over H 0 or H CDM 0 . However, for the joint analysis of H (z) and CMB data, the marginalisation over H 0 is calculated simultaneously: Let us now present the results for the f (R) model considered here.

Testing exponential F(R) gravity
By considering the SNe Ia, H (z), BAO and CMB data illustrated in the previous section, the above exponential model is well constrained. Here, we calculate these limitations and the best-fitted values of the parameters for the exponential F(R) model (3.1). After marginalizing over H 0 (and over ω b for the CMB data in χ 2 CMB ), the F(R) model (3.1) has three free parameters: β, This means that the model is assumed to be close to CDM. This assumption relaxes the difficulties to fit the free parameters, since the free parameters of the model can be automatically reduced, leading to two free parameters: β and  Fig. 1   . So in order to calculate the corresponding χ 2 , the value of the three parameters has to be given before solving numerically the system (3.9), (3.10) with the initial conditions (3.11), as described above. In this case every χ 2 function (after marginalisation over H CDM 0 and ω b for χ 2 H and χ 2 CMB ) will depend on β,   The same panels and notations of Fig. 1 are used in Fig. 2, but the blue contours corresponds to χ 2 3 ( CDM m , β) in Fig. 2 while the blue lines refer to the one-dimensional distributions χ 2 3 min ( CDM m ) and χ 2 3 min (β) in the bottom panels of Fig. 2. In order to compare these results with those obtained under the approximation (5.1), the curves of Fig. 1 are depicted as well, denoted by magenta lines for χ 2 3 min and by thin black dash-dotted lines for χ 2 tot min (in the bottom panels).
The black stars in Fig. 2 denote the minimum points of the two-dimensional distributions χ 2 tot ( CDM m , β) and χ 2 tot ( 0 m , β). Their coordinates (the optimal values of parameters) are tabulated in Table 3. In the same way, the minimum points for χ 2 3 are shown as the blue circles.  Fig. 2. These enlarged domains of suitable model parameters include the above-mentioned "islands" or local minima of χ 2 3 and χ 2 tot functions existed under the restriction (5.1). This effect is hidden in the top panels of Fig. 2 (it is beyond the 3σ confidence level), but it is shown in the bottom-right panel, where the local minima at β 0.4 from Fig. 1 (for the magenta and black dash-dotted lines) are naturally included in the general behaviour of χ 2 3 (the blue line) and χ 2 tot (the solid black line). These one-dimensional distributions determine the optimal values and 1σ errors of the parameter β; this information is included in Table 3, where the absolute minima of χ 2 and the mean of the model parameters are provided.
The optimal values and 1σ errors for CDM m in Table 3 are deduced from the one-dimensional distributions χ 2 min ( CDM m ). They are shown in the bottom-left panel of Fig. 2 as the solid black and blue curves in comparison with the corresponding plots for the case CDM m + CDM = 1 (the dash-dotted and magenta lines; they are taken from Fig. 1). The latter distributions coincide with the predictions of the CDM model. One can conclude from Table 3 (and Fig. 2) that the absolute minima for the F(R) model (3.1) χ 2 3 572.07 and χ 2 tot 575.51 are smaller than the ones for CDM model (572.93 and 583.24 respectively). Such a difference lies in the existence of degrees of freedom for the model (3.1).
In Fig. 3 we demonstrate how the numerical solution of the system (3.9), (3.10) behaves (the black solid lines in the phase space diagram) in comparison with the CDM model (the red dashed lines). We show the plots for the Hubble parameter E = H/H CDM 0 and the Ricci scalar R = R/(2 ) (in the right panels) depending on a in the top panels with the logarithmic scale and the same plots E(z), R(z) (with z instead of a) in the usual scale in the bottom panels. The model parameters for both models are taken from Table 3; they are optimal for χ 2 tot . One can see that for the optimal parameters the F(R) and CDM models demonstrate rather close dynamics of For the parameter CDM , the optimal values and 1σ errors in Table 3 were found after preliminary calculation of two-dimensional distributions χ 2 ( CDM m , CDM ) for χ 2 3 and χ 2 tot . The contour plots for these two-dimensional distributions and the corresponding one-dimensional plot for  Fig. 5 shows the evolution of the luminos-ity distance and the Hubble rate for the exponential gravity model and for CDM for their best fits. As shown, both curves fit the data similarly, such that the two models become indistinguishable.

Exponential model and inflation
The exponential F(R) model (3.1) considered in the previous sections, describes all observational manifestations of the late-time acceleration. However, such a model can also explain the early-time inflation when introducing some suitable modifications in the form of F(R) as follows [47]: These additional terms can generate the expected cosmological constant i during the inflationary era, when R is near or larger than R i . The natural number n > 1 helps to avoid the effects of inflation during the matter era when R R i and the last term γ R α in Eq. (6.1) is necessary for a successful exit from inflation. As pointed out in Ref. [47], the model (6.1) has the following properties: a de Sitter phase naturally arises during inflation in the high-curvature regime, the inflationary terms become negligible in the small curvature era R R i , inflation ends successfully and avoids anti-gravity effects and instabilities during the matter era. To satisfy all these properties, the following requirements are obtained over the free parameters: Here the constant γ in Eq. (6.1) is expressed as γ = ( i ) 1−α . In Ref. [47] the de Sitter solution with parameters n = 4, γ = (4 i ) 1−α , α = 5/2, R d S = 4 i was considered. The corresponding values here are = 4, R d S = i with α = 5/2 satisfy the condition (6.5). However, below we analyse a wider set of inflationary solutions which obey the observational limitations for the slow-roll parameters.
Let us now focus on the realisation of slow-roll inflation within the model (6.1) and its predictions. In order to do so, the scalar-tensor counterpart of F(R) gravities is more convenient than its original action, such that the F(R) gravities can be expressed in terms of a scalar field as shown in (2.4) and (2.5). By varying the action (2.4) with respect to the scalar field φ, it yields Here recall that the relations of Eq. (2.5) for the scalar field and its potential hold: In order to analyse slow-roll inflation for the model (6.1), the action (2.4) can be transformed into the Einstein frame by the following conformal transformation: The action (2.4) turns out to bẽ Here, we have redefined the scalar field to keep the kinetic term in a canonical form: The FLRW equations for the action (6.8) become simpler than working directly on the F(R) action: whereas the scalar field equation is given bÿ During slow-roll inflation the scalar field mimics an effective cosmological constant, what basically means that Hφ φ andṼ φ2 . Both conditions can also be expressed through the so-called slow-roll parameters . (6.12) These quantities remain very small when inflation occurs such that 1 and η < 1, while at the end of inflation 1. In addition, the slow-roll parameters (6.12) are related to the amplitude and scale dependence of the perturbations originated during inflation, such that the spectral index of the perturbations and the tensor-to-scalar ratio, are given by Since both values are very well constrained by the last data from Planck and Bicep2 collaborations [40,41], which give n s = 0.968 ± 0.006, r < 0.07. (6.14) Then we can test whether the model (6.1) is capable to satisfy such constraints. Firstly, let us analyse the action (6.1), as during inflation, R i , the action (6.1) can be approximated as follows: Then, by Eq. (6.9), the relation among the scalar field and the curvature is obtained: (6.16) Here recall that γ = ( i ) 1−α . During inflation R i , and consequently κφ 1, so Eq. (6.16) can be approximated by (6.17) The scalar potential yields Applying the slow-roll approximation to the above relation (6.19), the number of e-foldings yields As shown above in (6.14), these values lie within the allowed ranges provided by Planck, such that the model (6.1) can reproduce well inflation and then recover late-time acceleration, leading to a unified description of the universe evolution.
Let us now analyse the system of equations (2.9) during the inflationary epoch, which are reduced to Here λ i = i /(2 ), r i = R/R i , e i = e −r 4 i ,γ = ( λ i ) 1−α . In the top-left panel of Fig. 6, the F(R) function (6.1) is depicted for the values n = 4, α = 5/2, ψ = 0.883, = 3.871, and compared to the exponential gravity (3.1) and the CDM model. As shown, the two F(R) models are not distinguishable at low curvature regimes but they are when the curvature becomes very large, as during inflation. Moreover, the model (6.1) practically coincides with the CDM model in the range < R < R i (or 1 < R < i / ), while differs for R > R i because of the γ R α term. Finally for R ∼ H 2 0 , both F(R) models (3.1) and (6.1) behave similarly and deviate from CDM.
One can conclude that the difference between the models (3.1) and (6.1) is negligible during the radiation/matter era, ΛCDM: The inflationary solution is unstable [47]: after N 55 e-foldings the de Sitter solution decays and the evolution transforms into the CDM behaviour.

Conclusions
Exponential gravity may be considered as an alternative to the so-called concordance model. The main gravitational action studied along this manuscript and described in (3.1) represents an slightly correction to the usual Hilbert-Einstein action with a cosmological constant. Such correction is modelled by a free parameter, which has been called β, such that it controls the scale at which corrections to GR become important. As shown in some previous works [39,47], such a type of F(R) models can reproduce a late-time acceleration epoch. However, the aim of this work was to show in an accurate way whether such type of gravitational actions fulfil the necessary cosmological constraints leading to a reasonable bound on the crucial parameter β. Note that an equivalent analysis was performed in [50] and more recently in [51]. In our paper we have updated the constraints obtained in previous works by assuming new released data from the last years. In addition, we have also analysed the exponential model at late times by following two approaches: the first one by assuming the condition (5.1) and the second one by assuming a more general approach. The former provides a more restricted case as we are forcing the model to mimic CDM at the present time, while the latter keeps the model free, except for the initial conditions which are the same for the two cases. As expected, the more restrictive on the theoretical conditions are, the better the constraints on the free parameter turn out, as shown in the bottom panels of Fig. 2. On the other hand, whereas the β parameter is not upper bounded (recall that GR is recovered for β → ∞), the fits realised in the paper, where SNe Ia, BAO, CMB and H (z) data were used, provides a sufficiently small lower bound on β, which may have consequences at the perturbation level, an aspect that should be studied in the future. In addition, the values for the matter density m do not differ too much among the one given by the CDM model and the one from the exponential gravity model, either when some approximations are assumed or when the general case is considered. Moreover, the χ 2 min is a bit smaller for the exponential gravity case than in CDM, such that the best fit does not correspond to CDM, although the difference is not statistically significant. From a qualitative point of view, the tiny differences among exponential gravity and CDM can be shown by looking at the form of the action in (3.1). By a sufficient small exponent, i.e. a β parameter large enough, the Lagrangian mimics quite well the Hilbert-Einstein action with a cosmological constant, such that there is not significant difference on the cosmological evolution among both models, as shown in Fig. 4.
In addition, the exponential gravity action (3.1) can be extended in such a way that the new action (6.1) can reproduce inflation as well. Note that inflationary models within the framework of F(R) theories have been widely studied in the literature, as shown by one of the most successful inflationary models, the so-called Starobinsky inflation [42]. In this sense, we have proposed here a model where the exponential term responsible for the late-time acceleration is suppressed at the large curvature regime and consequently its induced cosmological constant while two additional terms may become important: a different effective cosmological constant i (much larger than ) and a power term R α . Recalling that Starobinsky inflation is described by a R 2 term, the inclusion of R α just generalised the latter and ensures a successful exit from inflation [47]. However, as shown in some previous papers, such exponent has to be 2 < α < 3 in order to avoid instabilities [47]. Here, we have extended the previous analysis by using the usual scalar-tensor counterpart and obtaining the explicit form of the potential for the scalaron. Then we computed explicitly the spectral index of curvature fluctuations during inflation. An example fulfilling all the requirements provides an spectral index that leads to a nearly invariant power spectrum and an almost null amplitude for the tensor modes, predictions in agreement with the last data released by the Planck collaboration. Note that the constant i establishes the energy scale at which the last terms in the action (6.1) become important, such that then the action also provides a quasi-de Sitter inflationary expansion, similar to Starobinsky model, where the R α term guarantees a successful exit from inflation. As shown, the values for the free parameters which provide the correct values for the spectral index and the scalar-to-tensor ratio, also avoid further corrections when inflation ends. Such model is then able to reproduce inflation and successfully exit.
Hence, we can conclude that the full gravitational Lagrangian (6.1) is capable of reproducing inflation and late-time acceleration in such a way that no other fields are required. Recently, one more extension of this type of exponential gravity with log-corrected R 2 term was proposed to explain the unified universe history (see Ref. [89]). As shown here, the Lagrangian satisfies the observational constraints, with a no statistical significant difference with respect to the CDM model, what means that they cannot be distinguished. Simultaneously, the F(R) model (6.1) provides the right predictions for the inflationary epoch. The next steps should be focussed on the analysis of cosmological perturbations and the possible effects of such Lagrangian at the astrophysical level.