Unification of a Bounce with a Viable Dark Energy Era in Gauss-Bonnet Gravity

In this work we shall demonstrate that it is possible to describe in a unified way a primordial bounce with the dark energy era, in the context of Gauss-Bonnet modified gravity. Particularly, the early time bounce has a nearly scale invariant power spectrum of primordial scalar curvature perturbations, while the dark energy era is a viable one, meaning that it mimics the $\Lambda$-Cold-Dark-Matter model and also is compatible with the Planck 2018 data on cosmological parameters. In addition, our analysis indicates that the dark energy era is free from dark energy oscillations, which occur in the context of $f(R)$ gravity. We further addressed the later issue by examining $f(R)$ extensions of Gauss-Bonnet models, and we showed that the $f(R)$ gravity part of the action actually produces the dark energy oscillations at redshifts $z\sim 4$.


I. INTRODUCTION
The dark sector of the Universe constitutes the most mysterious problems in theoretical physics and cosmology, since these two sectors control the evolution of the Universe to an 96% extent. The dark sector consists of two parts, the dark matter and dark energy part, and both still urge for a consistent explanation. With regard to dark matter, it is still a question whether it is controlled by a weakly interactive massive particle [1][2][3][4][5][6], or it is simply some manifestation of modified version of general relativity [7]. On the other hand, dark energy is the name with which the late-time acceleration of the Universe, firstly observed in the late 90's [8], is now known, and this mysterious dark energy era has attracted a lot of attention in the literature [9][10][11][12][13][14][15][16][17][18][19][20][21][22]. In all the theoretical approaches towards consistently describing the dark energy era, modified gravity is to date the most promising description, see for example the reviews [23][24][25][26][27][28].
Apart from the mysterious dark sector of the Universe, another major issue which hopefully in the next two decades will be explained, is the primordial post-quantum gravity era of our Universe. To date there are two candidate descriptions for this primordial era, the inflationary scenario [29][30][31] and the bouncing cosmology scenario [32][33][34][35][36][37][38][39][40]. Both descriptions produce a nearly scale-invariant power spectrum for the primordial scalar curvature perturbations, with the cosmological bounces having the attribute of also producing a cosmological evolution free from the initial singularity.
In the context of modified gravity it is often possible to describe in an unified way the early and late-time eras of our Universe, see for example [41,42]. In fact, modified gravity might serve as the only consistent description of dark energy beyond general relativity. The reason for this is simple, since general relativity can describe late-time acceleration in a restricted way, by using a scalar field which produces either quintessential or phantom evolution, and also a simple cosmological constant may describe a de Sitter evolution at late-times. But phantom scalar fields are not necessarily the best description for the late-time era, since these inevitably drive the Universe towards a Big Rip singularity [43], and also phantom fields can be sources of instabilities. Modified gravity can successfully provide a consistent late-time era, that can mimic a quintessential or de Sitter or even a phantom dark energy era, see for example the reviews [23][24][25][26][27][28] for more details on these issues.
One promising sector of modified gravity theories, is the Gauss-Bonnet gravity [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59], in the context of which, the Gauss-Bonnet invariant appears in the Lagrangian in a non-linear way. Also, extensions of general relativity which contain higher orders of the Riemann and Ricci tensors can be found in Refs. [60][61][62]. The focus in this paper is in general on f (R, G) theories of gravity [48,59,63,64,69], and specifically we mainly focus on theories of the form R + f (G), in order to avoid having primordial superluminal perturbation modes, but for reasons that will be explained shortly, we also study the late-time behavior of an f (R) + g(G) model. Our aim is two fold: firstly to find appropriate model of Gauss-Bonnet gravity that may describe in a consistent way the dark energy era, and secondly to demonstrate that in some Gauss-Bonnet models it is possible to provide a unified description of the primordial post-quantum era and the dark energy era with the same model. One of the models which we shall present in this paper, is capable of describing primordially a Type IV singular bounce, and at late-times a dark energy era, which is mimics the Λ-Cold-Dark-Matter (ΛCDM) model and produces values for the cosmological quantities of interest that are compatible with the Planck 2018 data on cosmological parameters [66]. With regard to the primordial Type IV singular bounce, the singularity is Type IV type, so it is a smooth singularity which does not affect the evolution of the Universe in an extreme way, such as the crushing types of singularities. In addition, this particular singular bounce was shown in an earlier work that it generates a nearly scale invariant power spectrum of the primordial scalar curvature perturbations, compatible with the latest Planck 2018 constraints on inflation. Apart from the fact that our Gauss-Bonnet model of the form R + f (G), can both produce a singular bounce primordially and a dark energy era compatible with the ΛCDM model and the latest Planck observations [66], one major outcome of our work is that this specific type of models produce a dark energy era free from dark energy oscillations, known to be present in f (R) gravity models at large redshifts z ∼ 4. In fact, in order to verify this issue, we also examined the late-time phenomenology of a model of the form f (R)+g(G), and as we demonstrate, this type of models can also be compatible with both the Planck 2018 observations and the ΛCDM model, but it is not free from dark energy oscillations. Thus our work indicates the fact that the dark energy oscillations are possibly due to the f (R) gravity sector.

II. MODIFYING THE GAUSS-BONNET GRAVITY THEORETICAL FRAMEWORK FOR THE DARK ENERGY ERA STUDY
The focus in this work is, as we already mentioned, the unification of a singular bounce with the dark energy era, and to our knowledge this is the first time that this proposal is quantitatively materialized. In this section we shall present the theoretical framework of a general f (R, G) gravity and we shall appropriately modify the Friedmann equation by using appropriate statefinder quantities, in order we study in an optimal way the late-time era. The starting point of our work is obviously the gravitational action and we shall assume an f (R, G) model accompanied by the presence of perfect matter fluids, with the following gravitational action, with R being the Ricci scalar, κ = 1 MP is the gravitational constant, where M P denotes the reduced Planck mass, and G signifies the Gauss-Bonnet invariant defined as G = R 2 − 4R αβ R αβ + R αβγδ R αβγδ with R αβ and R αβγδ being the Ricci and Riemann tensor respectively. Lastly, L (m) is the Lagrangian density of the perfect matter fluids, which contains all the information for non-relativistic matter, that is Cold Dark-Matter (CDM) and relativistic matter, so radiation. Furthermore, we shall assume that the cosmological background corresponds to that of a flat Friedman-Robertson-Walker (FRW) metric, with the line element being, where a(t) denotes the scale factor. As a result, the Ricci and Gauss-Bonnet scalar are reduced to simpler forms, which read, where H signifies Hubble's parameter defined as H =ȧ a and as usual, and the "dot" implies differentiation with respect to cosmic time t. Thus, by varying the gravitational action (1) with respect to the metric tensor g µν , the gravitational field equation is derived. Here, we shall separate our equations in space and time components, hence the equations of motion are, where for simplicity, we denote differentiation with respect to a scalar function with a subscript. Furthermore, as stated before, the matter density is comprised of both relativistic and non-relativistic particles and consequently is written as, where ρ d0 signifies the current value of the non-relativistic density and χ = ρr0 ρ d0 is the ratio of the current values of relativistic and non-relativistic matter. In addition, P (m) denotes the corresponding pressure which is connected to the matter density as, with ω i being the equation of state parameter for a specific kind of matter and i = d, r, either non-relativistic or relativistic matter perfect fluids. Both kinds are treated as a barotropic perfect fluid with continuity equations, The aim of our study is to derive a functional form for Hubble's parameter, hence only a single equation of motion is necessary. In the following, we shall utilize Eq. (5) which we aim to solve numerically, and to extract the form of Hubble rate during the dark energy era. Before we continue however, we shall perform certain changes which will facilitate our study. Specifically, as a dynamical variable we shall use the redshift, and we shall also introduce a statefinder variable y H (z) which we define shortly, in order to make the late-time study more concrete and easy to tackle numerically. So the cosmic time will be replaced by a more convenient variable which is the redshift z. From the definition of the redshift, where we assumed that at present time the scale factor is equal to unity, a new differential operator can be constructed by simply performing a differentiation on this particular relation, which in turn reads, where now, Hubble's parameter depends solely on the redshift, that is H = H(z). This operator is of paramount importance as each object in the equations of motion which is differentiated with respect to the cosmic time, can be transformed to a redshift-dependent quantity by using the above transformation. Below we quote some important quantities that will be used frequently in this paper, and these are transformed as, where the "prime" implies differentiation with respect to the redshift. Also we have, where X, Y take the values R, G. This is because, since R = R(t), G = G(t) and f = f (R, G). Both approaches are valid so the choice is up to the reader, however, even for the first case, a similar relation for the differential operator with respect to redshift applies. The second change which shall be made is a function replacement, and specifically, instead of using the Hubble rate, we shall use an appropriate statefinder function related to it. But before we continue, it is worth making certain changes in the equations of motion. Recalling Eq. (5), we shall treat each geometric term derived from the expression f (G) in the gravitational action (1) as a fluid, corresponding to dark energy, which will turn out to be a perfect fluid as well. Assuming that, then the continuity equation reads,ρ This continuity equation, as mentioned before, implies that the dark energy fluid is perfect, as is the case with the rest of the perfect fluids present. Consequently, the equations of motion obtain the familiar form of Friedman's and Raychaudhuri's equations, With these equations at hand, we define the new statefinder function y H (z) as, This is the new dimensionless function which will participate in the equations of motion instead of Hubble's parameter.
Since the dark energy density was defined in Eq. (18), this particular function is related to the Hubble rate, as follows, where m s is a mass scale defined as m 2 with H 0 being the current value of Hubble's parameter and Ω 0 (m) the current value of the matter density parameter. Their values will be assumed to be equal to H 0 = 67.4 ± 0.5 km sec×Mpc and Ω 0 (m) = 0.3153 which are both based on the latest Planck 2018 data [66]. This statefinder function will be used as a replacement for Hubble and its derivatives, as shown below, An observant reader might notice that Hubble's derivatives participate in the previous equations with these exact forms. Indeed the Ricci scalar and Gauss-Bonnet invariant time derivative contain the above expressions, so with this designation, all the previous equations can be rewritten easily. Furthermore, recalling equations (7) and (11), the density of matter is rewritten with respect to the redshift variable as follows, where χ ≃ 3.1 · 10 −4 , and with χ being the fraction of the present day energy density of the radiation over the cold dark matter fluids. Finally, we define the following parameters which we shall evaluate during the late-time era.
Concerning the dark energy Equation of State (EoS) parameter, this is equal to, while the dark energy density parameter is, Thus the aim in this paper, is to study some appropriate f (R, G) models and compare their behavior with the ΛCDM model. In order to extract the late-time behavior of each f (R, G) model chosen, we shall solve numerically the differential equation (5) with respect to the statefinder y H (z), and for appropriately chosen physically motivated initial conditions.

III. R + f (G) GRAVITY: A SINGULAR BOUNCE AT EARLY-TIMES AND A DARK ENERGY ERA AT LATE-TIMES
For an arbitrary f (R, G), the possibility of ghosts being present is nonzero. For the first model, we shall assume a ghost free case where function f (R, G) is replaced by R + f (G). As a result, equation Eq. (5) is rewritten as This is the general differential equation that must be solved in the interval [-0.9,10] for redshift with respect to y H (z) defined in Eq. (24). Let us assume now that the Gauss-Bonnet function is given by the following expression where c 1 and c 2 are auxiliary constants with mass dimensions [m] 6 and [m] 2− 4α 3α−1 respectively for consistency while α is an additional parameter which we shall assume it satisfies the condition α > 1. This is an interesting model due to the fact that for small values of G, the inverse term c1 G becomes dominant during the late time whereas c 2 G α 3α−1 is dominant during the early time where the exponent α 3α−1 is lesser than unity given that the condition α > 1 holds true. Therefore, this model is capable of describing both early and late time era and thus unifying them smoothly.
The singular bounce cosmology is basically realized by an R + f (G) gravity of the form f (G) ∼ G 1−2α 3α−1 as it was shown in Ref. [52]. This particular f (G) gravity was able to realize a Type IV singular bounce with Hubble rate H(t) ∼ (t − t s ) α , with α strictly greater than unity, that is α > 1, and t s is the cosmic time instance that the singular bounce occurs. From Eq. (32) it is apparent that for α > 1 the term ∼ G α 3α−1 is dominant at early times, and the term ∼ G −1 is subdominant during the early-time era. We shall quantify this shortly, but let us discuss in short the singular bounce generation by the term ∼ G α 3α−1 . Following Ref. [52], the primordial curvature perturbations are generated near the bouncing point t = t s , and exit the Hubble horizon after the singular bouncing point. The Hubble rate near the bouncing point could be of the order H I ∼ 10 13 GeV (borrowing the value of the Hubble rate from low-scale inflation studies), thus, the term ∼ G α 3α−1 for α = 3.30579 and c 2 = 1eV 2− 4α 3α−1 , which are the values of the parameters α and c 2 we shall use in the following, is of the order c 2 G α 3α−1 ∼ 1.11887 × 10 30 eV 2 , while the term c 1 G −1 is of the order c 1 G −1 ∼ 10 −88 eV 2 , for c 1 = 1eV 6 . Thus indeed, the term ∼ G −1 is significantly subdominant during the early-time era, near the singular bouncing point. In Ref. [52], we calculated in detail the power spectrum of the primordial scalar curvature perturbations, and it was found that it is equal to, From the above expression we easily derived the spectral index of the primordial curvature perturbations, which is equal to, It is easy to see that for α taken in the range α = [3.2999, 3.31], the spectral index becomes compatible with the latest Planck 2018 data [67], which constrain the spectral index to be n s = 0.9649 ± 0.0042. Also for α < −1 the spectral index can also be compatible with the Planck data, however the α < −1 case corresponds to a Big Rip singularity, thus it is not physically acceptable. The singular bounce which has a Type IV singularity, is more physically appealing, since the Type IV singularity is quite smooth and does not affect any physical quantities that can be defined on the three dimensional spacelike hypersurface which is defined at the time instance that the singularity occurs.
Having settled that the term practically generates the singular bounce at early times, and is dominant during this primordial era, let us see how things are modified at late times. Apparently, the late-time era is controlled by the term ∼ G −1 , but let us see this explicitly in a quantitative way for the moment. Let us use the current value of the Hubble rate which is H 0 ∼ 10 −33 eV, so the term G α 3α−1 is approximately of the order c 2 G α 3α−1 ∼ 8.44947 × 10 −46 eV 2 , while the term c 1 G −1 ∼ 10 132 eV 2 . Thus it is quantitatively apparent that the term ∼ G −1 is quite dominant at late times and controls the evolution.
In this section we shall demonstrate numerically that the term G −1 indeed dominates the late time era, generating a viable dark energy era. Our results are robust towards changes of the parameter α for all values larger than unity, but the analysis will be focused on those values of the parameter α that yield a spectral index of the primordial curvature perturbations compatible with the Planck 2018 data, so α = 3.30579. So let us proceed with the analysis of the model and express all the differential equations and the physical quantities in terms of the statefinder quantity y H (z). Before going to the details of our analysis, let us quote some useful expressions, so for the f (G) gravity chosen as in Eq. (32), we have, and thus Eq. (5) reads, 1000 for z f = 10 and there is a strong physical motivation for using these initial conditions, see for example Ref. [42]. In addition, the choice for Λ is Λ = 1.1895 · 10 −66 eV, while the mass scale m s is m s = 4.32552 · 10 −34 eV and in addition in the following we shall take c 1 = 1eV 6 , c 2 = 1eV 2− 4α 3α−1 , α = 3.30579. Then by solving numerically the differential equation (37), by using the aforementioned initial conditions and values for the free parameters, we shall analyze several statefinder quantities of cosmological interest for the late-time era, and we shall compare the results with the ΛCDM model and an f (R) gravity model which is known to produce viability for the late-time era. The comparison of the results of the Gauss-Bonnet model with the ΛCDM is obvious, since the latter is the cornerstone model of late-time phenomenology, since it is highly compatible with the CMB. However we need to discuss the comparison of the Gauss-Bonnet model with the f (R) gravity theory, since these two are apparently two distinct and phenomenologically competing theories. Our motivation is simply to investigate whether the dark energy oscillations at large redshifts (z ≥ 4) persist in the Gauss-Bonnet theory case. Our results are quite interesting, since in the Gauss-Bonnet case the oscillations do not occur. Also we shall calculate the predicted values for some quantities of cosmological interest and compare these values with the latest constraints of the Planck collaboration on these cosmological parameters [66]. For our analysis, we shall use the CMB based value for the Hubble rate, which is, [66], so H 0 = 67.4km/sec/M pc which is H 0 = 1.37187 × 10 −33 eV. Also let us discuss in brief the cosmological quantities and statefinders which we shall analyze and compare their behavior in this paper, for the ΛCDM, the Gauss-Bonnet Basically, the above quantity is itself a statefinder quantity, since it depends on the geometry through its explicit dependence on the Hubble rate derivatives. Another important cosmological quantity is the dark energy density parameter Ω DE (z), defined as Ω DE = ρDE ρtot , which in terms of the statefinder quantity y H (z) is written as, Now, for the comparisons with the ΛCDM model, the Hubble rate for the ΛCDM model is equal to, where again H 0 is the value of the Hubble rate at present time, namely, H 0 ≃ 1.37187 × 10 −33 eV [66], while Ω Λ ≃ 0.681369 and finally Ω M ∼ 0.3153 [66]. Also Ω r /Ω M ≃ χ, where we defined the parameter χ below Eq. (28). Finally, regarding the statefinder quantities we shall examine in this late-time study, our focus will be on four statefinder quantities, the deceleration parameter q, the jerk j, the snap parameter s, and finally the parameter Om(z), which in terms of the Hubble rate are defined as follows, Of course, the first three must be expressed in terms of the redshift variable, so for simplicity we quote the simplest expressions of the three, and these are, For the ΛCDM model, the statefinders j, s and Om(z) have the following simple values, s = 0, j = 1 and Om(z) = Ω M ≃ 0.3153. Let us now proceed to the results of our numerical analysis. Firstly we shall compare the Gauss-Bonnet model with the ΛCDM model, and in Fig. 1 we present the comparisons of the deceleration parameter (left upper plot), the jerk (right upper plot), the snap (bottom left plot) and the parameter Om(z) (bottom right plot), for the Gauss-Bonnet model (blue curves) and the ΛCDM model (red curves). As it is obvious, in the case of the deceleration parameter, the two curves are indistinguishable, while differences can be found for the rest three statefinders. Also it is mentionable that for the snap parameter, for redshifts z ∼ 4 and smaller, the ΛCDM model and the Gauss-Bonnet model are indistinguishable. Now in order to investigate also the dark energy oscillations issue, known to affect the f (R) gravity theories, we shall present the results of our numerical analysis for the Gauss-Bonnet theory, focusing on the dark energy EoS ω DE and the dark energy density parameter Ω DE (z), the behavior of which is plotted in Fig. 2 for the pure Gauss-Bonnet theory. As it can be seen in Fig. 2, no dark energy oscillations occur, but in order to make this result more clear, we shall compare the Gauss-Bonnet theory with a viable f (R) gravity theory, the functional form of which is [42], where, M is an auxiliary parameter with mass dimensions [m] and is given by the expression M = 1.5 · 10 −5 50 N M P , with N being the e-foldings number referring to the inflationary era (N = 60) and M P is the reduced Planck mass.
Essentially, the R 2 term contributes to the inflationary era and early time whereas R 3m 2 s δ for δ < 1 becomes dominant in the late-time era. In Fig. 3 we plot the behavior of the dark energy EoS ω DE (left) and the dark energy parameter Ω DE (right) for the f (R) (red curves) and the Gauss-Bonnet gravity (blue curves). As it is apparent, the behavior of the dark energy density parameter is indistinguishable between the models, however, the dark energy EoS for the f (R) gravity model has strong oscillations for z ≥ 4, which are absolutely absent from the Gauss-Bonnet model. Thus the dark energy oscillations plague that haunted the f (R) gravity models, is absent from the Gauss-Bonnet models. This claim is further supported by the plots appearing in Figs. 4 and 5. In Fig. 4 we compare the deceleration parameter q and the statefinder y H (z) for the f (R) gravity model (red curves) and the Gauss-Bonnet model (blue). The absence of oscillations in both cases are obvious, for the Gauss-Bonnet case, and the same conclusions can be derived if we look at Fig. 5 where we plot the statefinder parameters j(z) (upper left), s(z) (upper right) and Om(z) (bottom) for the f (R) gravity (red curves) and the Gauss-Bonnet gravity (blue curves). In the upper diagrams, it becomes apparent that the main difference lies in the dark energy oscillations which as expected is enhanced in higher order derivatives of y H , for the f (R) gravity case. Finally, in order to have a concrete idea of how well the Gauss-Bonnet model behaves, we shall compare the values of several statefinders, for the ΛCDM model and the Gauss-Bonnet gravity, and also we directly compare the dark energy EoS and the dark energy density parameters at present time for the Gauss-Bonnet model with the latest Planck data. Our results are summarized in Table III. As it can be seen in Table III, the statefinders values for the Gauss-Bonnet model are quite close to the corresponding ΛCDM value, and in addition, the Gauss-Bonnet model value for the dark energy density parameter is Ω DE (0) = 0.679553, which is compatible with the latest Planck constraints Ω DE = 0.6847 ± 0.0073 [66]. In addition, the dark energy EoS parameter value for the Gauss-Bonnet model at present time is ω DE (0) = −0.999667 which is in good agreement with the corresponding Planck constraint ω DE = −1.018±0.031 [66]. In conclusion, the Gauss-Bonnet model of Eq. (32) is able to produce a phenomenologically  viable early-time bounce, while at late-times it produces a viable dark energy era, which is in addition free from dark energy oscillations. In the next section, we shall also discuss a more complicated f (R, G) model, however less appealing in comparison to the model we presented in this section, due to the primordial ghost modes in those f (R, G) models. Let us now investigate the stability of the f (G) gravity solutions and discuss the stability of the cosmological perturbations as the general relativistic limit is approached. The FRW equations for the f (G) gravity, constitute a dynamical system, and the stability of the solutions can be examined if we perturb this dynamical system. We shall consider linear perturbations of the cosmological solutions, as functions of the Hubble rate, and the presence of an instability would indicate that the solution is not the final attractor for the theory at hand. For the era near the bouncing point, this is somewhat expected because the evolution continues after the bouncing point, but for the latetime era, the expected behavior is rather vague. This is due to the fact that the late-time era seems to be a de Sitter like solution by looking at the dark energy EoS parameter, but eventually it is not an exact de Sitter solution. So let us check explicitly the stability of the dynamical system towards linear perturbations of the cosmological solutions, in order to shed some light on this issue. The Hubble rate at the vicinity of the bouncing point is, and the f (G) gravity which approximately realizes the Hubble rate (49) is, Following the strategy of Ref. [68], we linearly perturb the solution of the Friedmann equation g(N ) = H 2 , in the following way, where g(N ) is a solution of the Friedmann equation, and N is the e-foldings number. By using the function g(N ), the Friedmann equation can be cast in the following form, where J 1 stands for, while J 2 is, and finally J 3 is, Now let us calculate explicitly the parameters J 1 , J 2 and J 3 for the era corresponding to the bouncing point, so the f (G) gravity is given in Eq. (50). By expressing the Hubble rate (49) as a function of the e-foldings number N , with the latter being equal N = ln a, the function g(N ) becomes, with ζ = 2α 1+α . Now let us proceed with the stability conditions, and by calculating the parameters J 1 , J 2 and J 3 , we get, and J 3 /J 1 , Focusing at the vicinity of the bouncing point t → t s , which correspond to N → 0, we get, where the parameter A is, and A is obviously positive. Hence J 2 /J 1 > 0 and J 3 /J 1 < 0, therefore the dynamical system of the cosmological equations is unstable near the bouncing point, as we anticipated. Now let us consider the stability of the cosmological solution for the late-time era, focusing at redshifts z ∼ 0. Since the late-time evolution results to an EoS for the dark energy approximately equal to −1, the late-time evolution for redshifts z ∼ 0, can be approximated by a de Sitter evolution, and it is realized by the approximate f (G) gravity of the form, Hence assuming that the late-time Hubble rate has the de Sitter form, H(t) ∼ H 0 , in effect, the function g(N ) in this case has the form g(N ) = H 2 0 . Let us repeat the procedure we followed for the case of the bounce near the bouncing point, so in this case, by perturbing the de Sitter vacuum solution g(N ) = H 2 0 using a linear perturbation of the form (51), the variable J 1 becomes, while the parameter J 2 reads, and finally the parameter J 3 reads, Now the fraction J 2 /J 1 is equal to J 2 /J 1 = 3, while the fraction J 3 /J 1 , is equal to, and since for the late-time analysis we took H 0 ≃ 1.37187×10 −33 eV and c 1 = 1eV 6 , this means that J 3 /J 1 is obviously negative. This means that the dynamical system corresponding to the linear cosmological perturbations is unstable for the late-time solutions. This means that the late-time dynamical system has an unstable de Sitter attractor (fixed point), thus in the era that lies beyond z = 0, the system will not possibly remain to the de Sitter vacuum, and physically this could mean that the late-time dark energy EoS will not remain to the fixed −1 value, but might evolve to the phantom or quintessence regime. In any case, this result is somewhat interesting, since it deviates from the ΛCDM description, in which case the de Sitter state is stable and unchanged, the cosmological constant is constant. Before closing, let us briefly discuss the stability of the cosmological perturbations for the R + f (G) gravity, as the general relativistic limit is approached. This issue was covered in detail in Ref. [69], so we shall report their result, which is that the condition f ′′ (G) > 0 suffices to ensure the stability of any solution in the General Relativistic limit. In our case, at the late-time era, where the f (G) is approximately f (G) ∼ c 1 /G, the quantity f ′′ (G) is equal to f ′′ ((G)) = 2c1 G 3 , hence the stability is ensured. However, for the early-time era, in which case the f (G) gravity has the form f (G) ∼ c 2 G α 3α−1 , the quantity f ′′ (G) reads, which is obviously negative for the values of α we chose, hence this solution is unstable towards to the general relativistic limit. However, this is not a problem in this case, because near the bounce era, no consistent general relativistic limit exists, since general relativity cannot describe the early-time era, while at the late-time, there can be some overlap between the general relativistic description and the f (G) description, because the late-time era is approximately described by a nearly de Sitter evolution.
IV. GENERALIZED f (R, G) GRAVITY LATE-TIME PHENOMENOLOGY: THE CASE OF f (R) + g(G) GRAVITY In the previous section we compared the results of the Gauss-Bonnet late-time phenomenology with the f (R) model of Eq. (48), and we demonstrated that even though both f (R) and R + f (G) models are capable of uniting early and late-time eras, there exist many differences, in particular the dark energy oscillations spotted on redshifts z ≥ 4, which are present in the f (R) gravity model while these are absent in the Gauss-Bonnet model. It is therefore sensible to try and find out whether combining the aforementioned f (R) model with an appropriately chosen g(G) function makes the dark energy oscillations disappear. Therefore, let us assume that the f (R, G) model is written as f (R, G) = f (R) + g(G) with, and where M , γ and δ are the same parameters as in the previous section, whereas λ and β are extra auxiliary parameters, the first being dimensionless while the latter has mass dimensions [m] 4 for consistency and moreover β is assumed to be positive. The model g(G) seems quite convenient for our approach since it has a linear part in terms of G and moreover the derivative of tan −1 (G) produces a term 1 1+G 2 . Since we assume an f (R) + g(G) gravity, then Eq. (5) is written as follows, This is the Friedmann equation, and this particular differential equation shall be solved numerically in the same redshift interval z = [−0.9, 10], by expressing the above differential equation in terms of the statefinder quantity y H . In the case at hand, we have, Suppose now that the values for M , m s , γ, δ, Λ remain the same as in the previous section, hence the reason they were not relabelled, and in addition λ = 1, β = 1. Then by using exactly the same initial conditions, and by solving numerically the differential equation (70), we present the results of our analysis in Figs. 6, 7 and 8. Particularly, in Fig. 6 we plot the deceleration parameter and and statefinder y H as functions of the redshift, and in Fig. 7 we plot the dark energy EoS parameter and the dark energy density parameter as functions of the redshift. Finally, in Fig. 8 we plot several statefinders for the combined f (R) + g(G) model (blue curves) and we compare these to the ΛCDM model results (red curves). Also in Table IV, we compare the values of several quantities of cosmological interest at present time for the combined f (R) + g(G) model, and compare these to the pure f (R) gravity model. An overall apparent result derived from the analysis, is that the presence of the f (R) gravity utterly affects the late-time phenomenology, since it brings along the dark energy oscillations in several statefinder quantities and cosmological quantities. This result seems to be model-independent, thus in conclusion, the f (R) gravity dark energy oscillations cannot be remedied by adding a Gauss-Bonnet term in the Lagrangian.  In this paper we quantitatively addressed the late-time phenomenology of Gauss-Bonnet theories, and we investigated how a singular bouncing cosmology occurring primordially and a viable dark energy era can be realized by Gauss-Bonnet gravity. This unification scheme was possible to be realized by using a Gauss-Bonnet gravity of the form R + f (G), which is free from primordial superluminal modes. The Type IV singular bounce at early times, realized by an appropriate Gauss-Bonnet gravity, generates a nearly scale invariant power spectrum of the primordial scalar curvature perturbations. In addition, the late-time driving part of the R + f (G) gravity generates a viable dark energy era, which mimics the ΛCDM model for some statefinder quantities, and is compatible at present time with the latest Planck data on cosmological parameters. This specific model, has another appealing attribute, the fact that the dark energy era is free from large redshift (z ∼ 4) dark energy oscillations, known to occur in f (R) theories. For demonstrative reasons, we compared the R + f (G) with the ΛCDM model and an appropriate f (R) gravity model which is known to generate a viable dark energy, and we provided qualitative evidence for the presence of dark energy oscillations only in the f (R) gravity case. In order to further analyze the impact of a non-trivial f (R) gravity term in the dark energy oscillations, we also studied the late-time phenomenology of an f (R) + g(G), and as we showed, the dark energy oscillations are also present in this case. Thus the f (R) gravity part in the gravitational action seems to always lead to dark energy oscillations at large redshifts.