Dynamics of Inflation and Dark Energy from $F(R, \mathcal{G})$ Gravity

In this work we study certain classes of $F(R,{\cal G})$ gravity which have appealing phenomenological features, with respect to the successful realization of the dark energy and of the inflationary era. Particularly, we discuss the general formalism and we demonstrate how several inflationary and dark energy evolutions can be described in the context of $F(R,{\cal G})$ gravity. Also we propose a unified model, in the context of which the early and late-time dynamics are controlled by the $F(R,{\cal G})$ gravity, thus producing inflation and the dark energy era, while the intermediate era is approximately identical with standard Einstein-Hilbert gravity. Also we calculate the power spectrum of the primordial curvature perturbations corresponding to the unified $F(R,{\cal G})$ gravity model we propose, which as we demonstrate is nearly scale invariant and compatible with the latest observational data constraints.


INTRODUCTION
The Universe at present time experiences an accelerating expansion, as standard observations based on Type IA supernovae indicate [1]. This accelerating expansion, dubbed dark energy era, is admittedly one of the most curious and surprising observations for our Universe, up to date. This is due to the fact that only a few people actually expected this late-time accelerating evolution, see for example [2]. A vast amount of works in theoretical cosmology is devoted towards describing the dark energy epoch, see for example [3][4][5][6][7][8][9][10][11][12], and various descriptions in the literature exist on this issue, varying from quintessential or alternative approaches [13][14][15][16][17], to modified gravity descriptions see the reviews [18][19][20][21][22][23]. Another quite intriguing epoch is the early-time era before the radiation dominated the expansion and evolution of the Universe. This is known as inflationary era, during which the Universe increased vastly its size in an exponential rate. The inflationary paradigm was introduced in the early 80's [25][26][27], as a theoretical proposal that solved the shortcomings of the standard Big Bang cosmological model. The latest observations coming from Planck, severely constrained the spectral index of primordial curvature perturbations and the scalar-to-tensor ratio corresponding to this early epoch, however no direct observational proof that inflation occurred is given up to date. It may be possible that the bouncing cosmology paradigm [28] may also be a viable description of the Universe at early times and even at late times.
In all the above cases, modified gravity [18][19][20][21][22][23] in its various aspects, may successfully describe the early and latetime acceleration eras in a unified way, see for example [24] for a pioneer work towards this aspect. Modified gravity is based on generalizing Einstein's gravity at large scales and for various curvature limits. Actually it seems that at strong curvature limits and at large scales, Einstein's gravity fails to comply with the observational data. With regard to the dark energy epoch, the only way to describe an accelerating expansion in standard Einstein-Hilbert gravity, is to use a phantom fluid or field, and intriguingly enough, this phantom accelerating evolution leads to a crushing type singularity known as Big Rip [29]. In standard modified gravity descriptions, higher order curvature corrections are included in the standard Einstein-Hilbert gravitational action, and many examples of modified gravity exist in the literature, see the reviews [18][19][20][21][22][23]. For the f (R) gravity case, a well-known and quite successful model, from an observational point of view, is the Starobinsky model first proposed in [30] which is consistent with both local gravity constraints as well as large scale constraints. Also a quite appealing model that unifies late and early-time acceleration in the context of f (R) gravity is the Nojiri-Odintsov model [24]. Alternative proposals to the f (R) gravity theory use Lovelock invariants, such as the Gauss-Bonnet scalar G, and several works exist that can successfully describe theoretically the dark energy epoch and also the inflationary epoch, see references [31][32][33][34][35][36][37][38][39][40][41][42][43] for an important stream of papers and recent works on the subject. Generalizations of the f (R) and f (G) gravity, are offered by higher order gravities [44][45][46][47], which use combinations of higher order curvature invariants constructed from the Ricci and Riemann tensors, that is R µν R µν and R µνσρ R µνσρ respectively. Also theories which combine the Ricci and Gauss-Bonnet scalars also exist in the literature, namely, F (R, G) gravity theories, and studies on this type of theories can be found in [35,[48][49][50][51][52]. In this paper, we aim to realize several inflationary evolutions of main interest in the context of some classes of F (R, G) gravity, and also we shall investigate how a late-time acceleration era can be realized from these theories. Our main interest will be on inflationary evolutions of physical interest, such as the quasi-de Sitter evolution, or variant forms of this. Finally, we discuss how the primordial curvature perturbations evolve in the context of F (R, G) gravity, and we propose a phenomenologically interesting model, for which we calculate the power spectrum of the primordial curvature perturbations, and the corresponding spectral index. As we demonstrate, the viability of the model is achieved if the free variables are suitably constrained. Also, with regard to the dark energy epoch, we propose several phenomenologically viable models that may describe successfully this epoch.
The paper is organized as follows: In section 2 we present the fundamental features of F (R, G) gravity in general, and we investigate how several inflationary scenarios can be generated by an appealing class of F (R, G) gravity models. Also we investigate how a quasi-de Sitter inflationary evolution can be realized by a more general F (R, G) gravity. In section 3, we present the realization of several dark energy scenarios in the context of a class of F (R, G) gravity models, and in section 4 we propose a unification model that can both describe the inflationary and the dark energy epochs, and also during the intermediate epochs it can reduce to an Einstein-Hilbert model. Finally, the conclusions follow at the end of the article.

F (R, G) GRAVITY AND COSMOLOGICAL SOLUTIONS
In this section we present the general formalism of F (R, G), and we will demonstrate how the gravitational equations of motion can be used as a general reconstruction technique in order to realize several inflationary evolutions. We shall be interested in a certain class of F (R, G) gravity, and we also present how a quasi-de Sitter evolution can be realized by a more general F (R, G) gravity. The gravitational action of vacuum F (R, G) gravity is equal to, where G stands for the Gauss-Bonnet invariant, defined as follows, By varying the gravitational action with respect to the metric tensor g µν , we obtain the gravitational equations of motion, which are, In the following we shall assume that the geometric background is a flat Friedmann-Robertson-Walker (FRW) metric, with line element, so for the FRW metric (5), the gravitational equations (3) become, where R stands for the Ricci scalar, and G is the Gauss-Bonnet invariant, which for the FRW metric (5) these become equal to, and also F R = ∂F ∂R and F G = ∂F ∂G . An appealing class of F (R, G) gravity, with phenomenological interest, as we demonstrate at a later section, is the following [32], where α is an even positive number and A is an arbitrary constant of the theory. We shall assume that the parameter A takes such values so that the term AG α is subdominant during the radiation domination epoch. Therefore we can assume that A = (1/G i ) α , where G i corresponds to the value of the Gauss-Bonnet curvature during the inflationary era. It is conceivable that the dominant term during inflation is also the Gauss-Bonnet curvature term AG α , and also if α > 1/2, we have approximately the relation G ≈ R 2 , so effectively a rational power of the Ricci scalar controls the evolution, which is valuable phenomenologically. At this point it is worth discussing the question how many degrees of freedom has the model (10). This question was addressed in Ref. [45], and we repeat the argument of that work, so for a general theory with action, where P = R µν R µν and Q = R µνσβ R µνσβ , the equation for the metric perturbations (gravitational wave modes) is, where, and also F 0 = ∂f ∂R R0 , and finally f Q0 = ∂f ∂Q R0 . For a general theory of the form (11), there are one massive scalar mode with mass m 2 s corresponding to equation (12), one massless spin-2 mode corresponding to the term k 2 in Eq. (13) and finally one massive spin-2 ghost mode corresponding to the k 4 term in Eq. (13). Hence the general theory (11) has two more degree of freedom compared to the Einstein-Hilbert case. However, the model (10) has only the massless spin-2 mode corresponding to the term k 2 in Eq. (13) and the massive scalar mode with mass m 2 s , due to the fact that in this case f Q0 = 1 and f P 0 = −4. Hence the model (10) has an additional degree of freedom in comparison to the Einstein-Hilbert gravity.
With regard to the inflationary era, we shall work assuming that the slow-roll conditions hold true, which in terms of the Hubble rate are expressed as follows,Ḣ Also, for later use, the first two slow-roll indices of the Hubble slow-roll expansion are defined as follows, and according to the slow-roll conditions, these two must be small during the inflationary era. Substituting Eq. (10) into the second Friedmann equation, and also by using Eq. (9), we get, The above form is highly complicated, hence in order to get an analytical solution, we apply the slow-roll conditioṅ H H 2 ≪ 1 during the inflationary era. As a result, we can expand the terms like (1 +Ḣ H 2 ) α up to first order. Simplifying the result, we finally get, The above equation can have three possible solutions, which are, (α − 1) = 0 (19) or/and, H 4α (1 + αḢ H 2 ) = 0 (20) or/and, and since, H = 0, we get or equivalently, where c is an arbitrary integration constant. The corresponding solution for the scale factor a is given by, where A 1 is an arbitrary integration constant. Finally the third solution can be found by solving, In order to obtain the evolution, we apply the slow-roll conditions, and the solution for H is, where b is again an integration constant. Accordingly, the scale factor is, where A ′ is an integration constant.

Reconstructing the Universe Evolution: The e-foldings Number Approach
As we have seen previously, the equations of motion are highly complicated, hence only approximate solutions can be obtained analytically. Therefore, in order to obtain exact solutions, without neglecting the Einstein term R, we shall employ the reconstruction technique firstly developed in [53,54] and then to F (R) gravity in [55,56]. Following the analysis of [55], we express all the variables as functions of the e-foldings number N = ln(a/a 0 ) (with a 0 being the scale factor at present time, which can be set equal to one for simplicity) instead of the cosmic time. Using the relation dN = Hdt the second Friedmann equation appearing in Eq. (6) can be rewritten as follows, where we have used Eq. (9), expressed as follows in terms of the e-foldings number, In the above equations, the prime indicates differentiation with respect to the e-foldings number. Eq. (30) can be rewritten as follows, Solving with respect to H 2 (N ) we get, Thus we see that H is now a function of R, G, F (R, G) and of higher derivatives. Using the expressions for slow-roll parameters defined previously, in terms of H as given in Eq. (16), one can now express all the slow-roll indices as functions of R, G, F (R, G) and of higher derivatives (all being functions of the e-foldings number N ).
Using the relation dN = Hdt, we can rewrite the slow roll parameters such that t is replaced by N in the following way, Also from Eq. (30), we have, Using this information and substituting the expression for H 2 (N ) from Eq. (32), we may express the slow-roll parameters as follows, Thus we can see that given the F (R, G) function and also the functional form of H 2 (N ), one can find the expressions for the slow-roll parameters. If assumed that the F (R, G) gravity acts as a perfect fluid, one may calculate the spectral index of the primordial curvature perturbations and of the scalar-to-tensor ratio by using the usual expressions corresponding to a canonical scalar field. However, we need to stress that a more formal approach requires the calculation of the spectral index directly by calculating the power spectrum of primordial curvature perturbations, as we formally do in a later section for a phenomenologically appealing F (R, G) gravity.
What we have at hand up to now in this section, is a reconstruction method that can enable us, given the Hubble rate as a function of the e-foldings number, to find the F (R, G) function that may realize such an evolution, as a function of the e-foldings of course. The reconstruction method consists of the following steps, firstly we assume a specifically chosen functional form of the Hubble rate, for example, Then by using Eq. (30), we obtain, The Friedmann equation can now be written as follows, Thus, by appropriately choosing the function P (N ) and then by using the functional relation between R, G and the number of e-foldings, the resulting equation becomes a second order differential equation for the function F (R, G), whose solution can be obtained in principle. On the other hand, we can also assume a particular form of F (R, G) and solve the above equation for P (N ), however in general this might be technically more difficult.
In order to demonstrate how the method works in practise, let us present some illustrative and relatively simple examples. We start off by assuming that the Hubble rate has the following form [55], which in terms of the cosmic time is an exact de Sitter evolution H(t) ∼ e P0t . We shall seek which F (R, G) gravity of the form, can realize the above evolution. Using Eqs. (41) and (39), we obtain, Substituting the above forms into Eqs. (40), we find that the resulting form of the differential equation is still too complicated. In order to simplify the resulting equations, we shall assume that the slow-roll condition holds true, thus by neglecting H ′ terms in comparison to H 2 terms, we obtain, Using the above form for the Gauss-Bonnet scalar, we get, Thus Eq. (40) is simplified as follows, Accordingly, by using Eq. (42), we get, The solution of the above equation is the following, Here F 0 and F 1 are integration constants chosen such that the Gauss-Bonnet term has minimum effect as the Universe enters the radiation domination era, and G p q m n (z) is the Meijer G function. It should be noted that the first term in Eq. (50) does not affect the dynamics of the cosmological system, since it contributes a total divergence. Actually, only the second and third terms affect the dynamics. A detailed analysis of the Meijer G function in the large G regime, indicates that the second term in the solution (50) is dominant, thus we have during the inflationary era, This practically is the result of the reconstruction method we presented, so the exact de Sitter evolution of Eq. (41) is realized during the slow-roll inflationary era by the F (R, Let us consider another illustrative and relatively simple example, in which case the Hubble rate as a function of the e-foldings number N is, with P 2 > 0. Such exponential behavior gives rise to a power law scale factor as a function of the cosmic time, which can give rise to an inflationary evolution. In this case we have, and in effect, the differential equation for the F (R, G) is, which can be analytically solved and the solution is, where, Here F 0 , F 1 are integration constants. As it can be seen from the functional form of the solution, the first term is ∼ G 1/2 , hence it behaves like general relativity. The dominant term in the large Gauss-Bonnet curvature regime is the second term, and for it to be the dominating one, we should have a < 1. For example, if a = 0.5, then we get the simplest modification with G 2 . In this case, β = −1.5 or 0.38. The reconstruction method we presented works if the Hubble rate H(N ) is given, but the method can work the other way round, that is to fix the F (R, G) gravity and from it to find the Hubble rate H(N ). Indeed, by choosing a particular form of F (G) as was done earlier, substituting it into Eq. (40) and using Eq. (39) we get a differential equation for P (N ), solving which yields the form of H 2 (N ). Let us demonstrate this by choosing the functional form of F (R, G) to be, Using the above form, we get the following second order differential equation for the function P (N ), where we used the approximation (1 + P ′ (N )/P (N )) α ≈ (1 + αP ′ (N )/P (N )), which are imposed by the slow-roll approximation. Solving the above differential equation with respect to P (N ), we get the form of H 2 (N ) which is, with β = 2α. Thus, constraining β, one can in turn constrain α.

Quasi-de Sitter Evolution in the Context of a General F (R, G) Gravity
In principle, in the context of F (R, G) gravity, any cosmological evolution can be realized, however the resulting differential equations to solve are not always easy to tackle. Nevertheless, simple cosmological solutions can easily be handled, and a particularly interesting and phenomenologically valuable evolution, is the quasi-de Sitter evolution, in which case the Hubble rate is, where H 0 and H i are dimensionful parameters of the theory. In general, the latter two parameters are controlled by the modified gravity governing the evolution, and are constrained by the observational data. For example, in the case of R 2 gravity, a concise treatment on the allowed values for H 0 and H i can be found in Ref. [57], but for the purposes of this paper we shall assume that H 0 ≫ H i , which is justified by the conditions of having a nearly de Sitter evolution. We shall be interested in finding which F (R, G) gravity may generate the quasi de Sitter evolution (61), so we shall use the reconstruction method of Ref. [48]. We rewrite the gravitational action as follows, and upon variation with respect to t, which is an auxiliary scalar field identified with the cosmic time, we get, where the "prime" denotes differentiation with respect to the auxiliary scalar "t". From Eq. (63) one can determine the relation t = t(R, G), then if someone finds the functions P (t), Q(t) and Z(t), then the F (R, G) function is found by substituting t = t(R, G) found from Eq. (63). By using the equation of motion (6), we find that the differential equations which determine the functions P (t), Q(t) and Z(t), are the following, Thus by determining the functions P (t), Q(t) and Z(t), which can be found by directly solving the differential equation (64), the explicit form of the F (R, G) gravity can be found. Let us apply this method for the quasi-de Sitter evolution (61), so we assume for simplicity that the function Z(t) has the following form, where b and β are arbitrary real parameters. Then, by substituting Z(t) in Eq. (64), we obtain the differential equation that will yield the function P (t). However, the resulting differential equation is very complicated to even quote it here, and also to solve it, so some leading order approximation is needed. Emphasizing at early times, the differential equation (64), becomes at leading order, where the function S(t) is equal to, The differential equation (67) is easy to solve, so the resulting solution is, where C 1 and C 2 are integration constants, and µ = H 2 0 + 8H i . Substituting Eqs. (66) and (69) into Eq. (65), we obtain the function Q(t), which at leading order is, (H0+µ) . Finally, by combining Eqs. (66), (69), (65) and (63), we can obtain the relation t = t(R, G). By substituting the resulting expression for t = t(R, G) in F (R, G) = P (t)R + Z(t)G + Q(t), we obtain the approximate resulting form of the F (R, G) gravity which realizes the quasi-de Sitter evolution (61), which is, where the parameters Γ i , i = 1, ..., 4 and γ can be found in the Appendix. We need to note that there is in principle much freedom in choosing the function Z(t) and determining the resulting F (R, G) gravity, therefore it is possible that functionally different F (R, G) gravities may realize the same cosmological evolution, at least in the vacuum case. However, the complicated equations require extensive simplifications, so from now on we shall consider F (R, G) gravities of the form F (R, G) = R + f (G).

DARK ENERGY FROM F (R, G) GRAVITY
Having discussed the realizations of inflationary evolutions in the context of F (R, G) gravity, in this section we shall present several viable models of F (R, G) which may successfully generate a dark energy era. Some simple but in general phenomenologically viable models are given below, where G 0 corresponds to the present value of the Gauss-Bonnet scalar and F 0 is a constant which will be chosen to be equal to the present value of the cosmological constant. In this work we will mainly be interested in reproducing the dark energy era, using exponential F (G) gravities, since this functional form of F (R, G) gravity is relatively simple and it can easily satisfy the solar system constraints.
Since G = 24H 2 (H 2 +Ḣ) ≡ 24H 2ä /a, the termä changes sign, that is, goes from negative to positive, as the Universe transits from the decelerating to the accelerating phase, therefore G also goes from negative to positive. Thus in the model I in Eq. (76), the change in sign of G affects the exponential part whereas in the model II, since we have |G|, the change in the sign has no direct effect.
Neglecting radiation and matter fluids, thus considering the models in vacuum, de-Sitter points correspond to a constant Hubble rate value, say H = H 0 (which will correspond to the present value of Hubble constant in case we are dealing with late-time acceleration), and the corresponding equation that determines the de Sitter points is given by, where G 0 = 24H 4 0 . Let us consider first the model I in Eq. (76), so by using Eq. (76), we get F 0 = 3H 2 0 (78) In order to study the stability of the de Sitter point, let us consider a linear perturbation around the de Sitter point H = H 0 + δH, so we get the perturbation equation which is [51], which can be solved to yield the solution, where c 1 and c 2 are integration constants [51]. Thus in order for the de Sitter points to be stable, the following conditions need to be satisfied, Hence the stability depends on the sign and the magnitude of F GG . For the model I F (G), one gets the form of λ = 3H 0 /2(−1 ± √ 2.4) and the quantity H 6 0 F GG (H 0 ) = 3e 1 /(24 2 ) > 1/384. Thus, this model has slightly unstable de Sitter points. This implies that the late-time de Sitter attractor is unstable, however from a phenomenological point of view, this is not a serious issue, since we do not know how the Universe will evolve in the future, that is, if the late-time de Sitter attractor is stable or not.
For the model II, we have, and in order to have stability, the following condition must hold true F GG > 0, so in this case we must have λ = 3H 0 /2(−1 ± √ 2) and also H 6 0 F GG (H 0 ) = 3e −1 /((1 − 2e −1 ) × 24 2 ) > 1/384. As a result, in this case too we have an unstable de Sitter point as a solution.

Study of the Equation of State
Let us investigate the behavior of the effective equation of state (EoS) for the above models, so we rewrite the equations of motion as follows, where we have included the presence of matter, and ρ G and p G are given by, Thus the EoS parameter corresponding to the geometric dark energy is w = p G /ρ G . At the de Sitter point, by substituting H = H 0 and G = G 0 as constants, we get ρ G = −p G , as was expected. Let us consider model I, so in order to get the evolution for the Hubble parameter, and in effect to obtain the EoS corresponding to late times (so for G < G 0 ) we need to solve the Friedmann equation for the model I of Eq. (76). Since the resulting equation is very complicated, we apply certain approximations by using the fact that we are interested for low curvature terms, so one can expand the exponential term and keep up to the first non-zero contributions of the Gauss-Bonnet term in the action we get F (R, . Then substituting the above form in the Friedmann equation and by neglecting the Einstein term one finds that the approximate solution for Hubble parameter is given by, which denotes an accelerating Universe with an approximate power law scale factor. The EoS for any modified gravity is equal to [18], which in the case at hand becomes equal to w = −1 + 2 3A . Therefore, depending on the sign of the parameter A, one may have quintessential acceleration (for A > 0 and for values which render w < −1/3) or a phantom evolution (for A < 0). Also, for large values of A, the evolution is approximately a de Sitter evolution. The same results are obtained for the model II in Eq. (76), so we omit the study of the model II for brevity.

Stability of Matter-Radiation Points
Let us now focus on the stability of matter and radiation points solutions. In order to study the stability, we consider the following equation [51], We shall study perturbations of linear order and of the following form, By keeping linear order perturbations, the equation that determines the stability of matter and radiation points is given below [51], where we have assumed a ∝ t p . Then the stability conditions are the following [51], In addition, the modified gravity term originating from the functional form of the F (R, G) gravity, should also vanish or become subdominant during the matter and radiation domination epochs, so F (R, G) −→ R/2 as |G| ≫ G 0 . Let us firstly consider the model I in Eq. (76), so we need to stress that G is negative during the matter and radiation dominated epochs, hence the signs of G and G 0 are opposite in the exponent of the model I in Eq. (76). In effect, an exponentially decreasing function is obtained in the regime |G| >> G 0 . As can be seen from Eq. (76), and for the model I, for G >> G 0 , F GG → 0+. Also in the large curvature regime, F (R, G) = R/2 − F 0 ≈ R/2, thus we obtain a behavior similar to that of Einstein-Hilbert gravity. Let us now consider model II in Eq. (76), and since the absolute value of G appears in the exponent, we have a decreasing function of G, irrespective of the sign of G, thus satisfying all the above conditions. Therefore, as a result, both the models appearing in Eq. (76) can describe an accelerating late time evolution, but also can ensure an Einstein-Hilbert-like evolution during the radiation and matter domination epochs.

Reconstruction of Several Cosmic Scenarios
Before closing this section, let us demonstrate how to realize several cosmological scenarios by using the models (76). Firstly, for a de Sitter Universe, the Gauss-Bonnet and the Ricci scalars are constant (and can be set equal to their present value), so we have, Following the analysis of Ref. [34], the resulting gravitational equations of motion take the following form (after taking the trace and writing F (R, G) = R + F (G)), In the ΛCDM model, we have R 0 = 4Λ, thus the effective cosmological constant in this theory is equal to, Hence, for the models I we have, while for the model II, we have, In order to obtain several cosmic evolutions realizations, we shall apply the reconstruction techniques we presented in earlier sections. In this way, we shall be able to realize several cosmic scenarios so by following the analysis of [58], the Friedmann equation is, 6P (N (G)) +24(48(P (N (G))) 3 P ′ (N (G)) + 12(P (N (G))) 2 (P ′ (N (G))) 2 + 12(P (N (G))) 3 P ′′ (N (G)))F GG − (24(P (N (G))) 2 + 12P (N (G))P ′ (N (G) Thus considering different choices of the Hubble rate H, we can get the corresponding solution F (G).

ΛCDM Evolution
Let us firstly realize the ΛCDM evolution, so consider, where H 0 and ρ 0 are constants and these correspond to the value of the Hubble rate at present time, and to the cold dark matter respectively. By using Eq. (30), and upon keeping the leading order terms, we may obtain the expression of the e-foldings number N in terms of G as follows, Using the above expression, like in the cases for inflation, one may convert the above equation into a differential equation for the F (G) gravity. The resulting form of the equation is very complicated, however, since we are interested in the low curvature regime, we can keep low order terms in G. By doing this, we obtain F (G) = −c(1 − e aG ), which is identical to the model I of Eq. (76) at leading order. In effect, the model I can describe the ΛCDM successfully, without the need for a cosmological constant term.

Phantom Evolution
Let us now demonstrate how to obtain a phantom behavior, without the actual presence of any phantom fluid, so we assume that the function P (N ) is [58], Following the above analysis, the e-foldings number N as a function of G is equal to, so by solving the above equation, we get, The differential equation above can be solved to yield, where α, β are functions of H 0 and a, b are integration constants. Thus we see that a phantom dark energy era can also be obtained by some vacuum F (R, G) gravity without the need for a phantom fluid.

UNIFYING DARK ENERGY WITH INFLATION IN THE CONTEXT OF F (R, G) GRAVITY
Based on the previous sections, we can propose a phenomenological viable model that can describe in a unified way the dark energy and the inflationary epochs, but also that reduces to the standard Einstein-Hilbert description during the matter and radiation domination eras. The functional form of the proposed F (R, G) gravity is the following, and based on the previous sections analysis, the last term in Eq. (104) is dominant only in the low curvature regime, if α is sufficiently larger than unity, for example α = O(10), so at late times, the F (R, G) gravity mimics the ΛCDM model, since we have approximately F (R, G) ∼ R + F 0 . Also in the high curvature regime, the F (R, G) gravity is Gi α , hence the inflationary era is dominated by a power law Gauss-Bonnet term. The model of F (R, G) gravity appearing in Eq. (104) has great similarities with the inflationary exponential f (R) gravity models used in Refs. [59,60], to successfully describe the unification of inflation with the dark energy era. For example, such an exponential f (R) gravity model has the form [59,60], Now we shall investigate whether the predicted spectral index of the primordial curvature perturbations for this theory can be compatible with the latest Planck observational data. In order to calculate the spectral index, we shall firstly calculate the evolution of the power spectrum curvature perturbations. Also we shall investigate the behavior of the scalar primordial curvature perturbations, with respect to their evolution after the horizon crossing during the inflationary era. For the study of the scalar perturbations during the inflationary era, we assumed that the perturbed metric is the following [18], where ψ, β, γ and δ represent the smooth perturbation functions. We shall quantify our study of perturbation evolution in F (R, G) gravity, by using the comoving scalar curvature perturbation, which is defined as follows, which is a gauge invariant quantity. It can be shown (see for example [18]), that the differential equation which governs the evolution of scalar perturbations is the following, where Q, B 1 , and B 2 are functions of the F (R, G) function, of the Hubble rate H and their higher derivatives. We need to note that Eq. (108) refers to a general F (R, G) gravity and it has a complicated form and also the last term proportional to k 4 is only present for a general F (R, G). As we will show, in the case of an F (R, G) gravity of the form F (R, G) = R + f (G), the above equation is significantly simplified. For the model (104), the master differential equation (108) simplifies since B 2 is equal to zero, and also Q(t) and B 1 (t) are equal to, In effect, the master equation (108) can be rewritten as follows, As we already mentioned, the F (R, G) gravity during the inflationary era is approximately equal to ∼ AG α , which as we showed in previous sections, it generates an approximate power law evolution with scale factor a(t) ∼ t β . Hence, the functions Q(t) and B 1 (t) are equal to, Solving the differential equation (111) analytically is a formidable task, without using an approximation, so we can finding the leading order behavior for small cosmic times, which is, where C 0 , C 1 , C 2 are equal to, Upon solving the above differential equation, we obtain the following solution, where the parameters µ, F 1 and F 2 are given below, and P 1 and P 2 are (scale) k-dependent constants to be determined later on. By taking the small argument limit of the Bessel function J µ (z) in the solution (115), we have approximately, The power spectrum of the primordial scalar curvature perturbations at horizon crossing is given by, and we shall reveal the k-dependence of each term in the expression of the power spectrum. From Eqs. (114), (118) and (119), the parameters C 2 , F 2 and m have a scale dependence, which is, therefore we have, Also the cosmic time can be expressed in terms of the wavenumber k, by using the horizon crossing condition k = aH, so at leading order we have, It is important to define the value of the scalar field to have some initial vacuum state value. To this end, we introduce the variable u = z s Φ, where, z s = Q(t)a(t). Then the scalar perturbations action is written in terms of u in the following way, where the prime indicates differentiation with respect to the variable τ = a −1 (t)t. So we assume a Bunch-Davies initial state for the vacuum of the scalar field prior to the inflationary era so u ∼ e −ikτ √ k . The imaginary phase which depends on the conformal time will be eliminated since eventually we will be interested in the expression |Φ(t = t s )| 2 . By performing the calculation we obtain at the initial time t i prior to inflation that, so by using an approximate functional form for the function Q(t), we have, and eventually, P 2 reads, By combining the above resulting behaviors for the parameters, the power spectrum as a function of the wavenumber k is equal to, and we can easily extract the spectral index for the primordial scalar curvature perturbations, which is, Thus, by choosing for example β = 177.41 we obtain approximately n s ≃ 0.96, which is compatible with the 2015 Planck data [61]. Before closing this section, it is worth investigating how the cosmological perturbations evolve after the horizon crossing during inflation, and the differential equation that governs this evolution is [18], which can be solved to yield, and the function Q(t) is given in Eq. (126) and a(t) ∼ t β . Thus, by evaluating the second term in the above equation we get approximately, Hence, as the cosmic time increases, the perturbations after the horizon crossing remain approximately constant, which is a good feature, since this will render them relevant to present day observations, after they reenter the horizon during the radiation domination era.
In conclusion, the proposed model in Eq. (104), apart from providing a theoretical unified framework which can qualitatively describe inflation and the dark energy era, it also provides a quantitatively viable inflationary era, compatible with the observations.

CONCLUSIONS
In this paper we studied certain classes of phenomenologically appealing F (R, G) gravities. We focused on the successful realization of the inflationary and dark energy eras, but also we investigated how the intermediate matter and radiation domination eras can be realized. As we discussed, the most appealing class of F (G) is types of theories which have the form F (R, G) = R + f (G), since for these no superluminal modes appear in the evolution of scalar perturbations, which would be proportional to k 4 , only k 2 modes appear. In our study we presented how various inflationary and dark energy scenarios can be realized with this type of F (R, G) gravity, and also we investigated how a quasi-de Sitter evolution can be realized with a more general form of F (R, G). Finally, we proposed a unification model, which at early times produces a nearly scale invariant power spectrum of primordial curvature perturbations, at intermediate times, and particularly, during the matter and radiation domination eras, it behaves like ordinary Einstein-Hilbert gravity, and at late-time it produces a dark energy era. This model should be accompanied by two ordinary perfect fluids, describing the matter and radiation domination eras, and these produce insignificant effects at early and late times. The study of more general F (R, G) gravities in detail is a demanding task, and to some extend problematic, due to the fact that the calculation of the scalar primordial curvature perturbations is highly non-trivial and too complicated, plus superluminal modes proportional to k 4 powers of the wavenumber appear. However simpler models of F (R, G) gravity like the one appearing in Eq. (104) can provide a fertile ground for cosmological phenomenology, which is also free from k 4 instabilities in the perturbations.