Extended matter bounce scenario in ghost free f ( R , G ) gravity compatible with GW 170817

E. Elizalde, S. D. Odintsov, V. K. Oikonomou, Tanmoy Paul 1) Institut de Ciències de l’Espai (ICE-CSIC/IEEC), Campus UAB, c. Can Magrans s/n, 08193, Barcelona, Spain. 2) ICREA, Passeig Luis Companys, 23, 08010 Barcelona, Spain 3) Institute of Space Sciences (IEEC-CSIC) C. Can Magrans s/n, 08193 Barcelona, Spain 4) Institute of Physics, Kazan Federal University, 420008 Kazan, Russia 5) Department of Physics, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece 6) International Laboratory for Theoretical Cosmology, Tomsk State University of Control Systems and Radioelectronics (TUSUR), 634050 Tomsk, Russia 7) Tomsk State Pedagogical University, 634061 Tomsk, Russia 8) Department of Physics, Chandernagore College, Hooghly 712 136. (9) Department of Theoretical Physics, Indian Association for the Cultivation of Science, 2A & 2B Raja S.C. Mullick Road, Kolkata 700 032, India

In the context of a ghost free f (R, G) model, an extended matter bounce scenario is studied where the form of the scale factor is given by a(t) = (a0t 2 + 1) n . The ghost free character of the model is ensured by the presence of a Lagrange multiplier, as developed in [1]. The conditions under which, in this model, the speed of gravitational waves becomes equal to the speed of light (equal to one, in natural units), thus becoming compatible with the striking event GW170817, is investigated. It is shown that this happens for a class of Gauss-Bonnet (GB) coupling functions (h(t)) which satisfies a constraint equation of the formḧ =ḣH, with H the Hubble parameter. This constraint is then imposed on the ghost free f (R, G) gravity theory to be consistent with the GW170817 event, subsequently, the corresponding non-singular bouncing cosmology with the aforementioned scale factor is extensively studied. The forms of the coupling function and Lagrange multiplier in the "low curvature limit" of the theory are reconstructed, yielding a viable approximation for n < 1/2. Correspondingly, by solving the cosmological perturbation equation, the main observable quantities, namely the spectral index, tensor to scalar ratio, and the running index are determined and confronted with the latest Planck 2018 data. Consistency with the data is proven for those parametric regimes that which correspond to n < 1/2. This makes the low curvature approximation a viable one for calculating the scalar and tensor power spectra.

I. INTRODUCTION
A major issue in theoretical cosmology is to ascertain whether the Universe was created from an initial singular or if, on the contrary, the Universe initiated its expansion from a non-singular, bounce-like stage. In other words, if the Universe evolved according to the standard cosmology evolution, for which the Universe emerges from an initial spacelike singularity, or if, in accordance with the bouncing theory, the Universe started its expansion with a non-zero value of the scalae factor. The inflationary picture [2][3][4][5][6][7][8][9][10][11][12][13] has earned a lot of support, since it is able to resolve the horizon and flatness problems, and also generate an almost scale invariant power spectrum which is perfectly consistent with the observational data. However even the most recent observations could not yet undoubtedly confirm that inflation indeed took place, a direct proof of this fact (as would be, e.g., the direct or indirect detection of associated primordial gravitational waves) is still missing.
One of the alternative theories different from inflation is the bouncing scenario . This theory is also able to generate a nearly scale invariant power spectrum, thus becoming compatible with the observational data available. Moreover, a definite advantage of the bouncing scenario is that it avoids any spacetime singularity. In fact, the Big-Bang singularity could just be a manifestation of the shortcomings of classical gravity, which is certainly unable to describe the physical evolution at such small scales. Quite possibly, a yet-to-be-built quantum theory of gravity might be able to resolve the Big Bang singularity, as turned out to be actually the case in classical electrodynamics, with the Coulomb potential singularities at the origin of the potential, which are resolved in the context of quantum electrodynamics. However, in the absence of a fully accepted quantum gravity theory, the bouncing scenario reveals as a most promising theory to deal with this issue.
Among the various bouncing models proposed so far, the matter bounce scenario (MBS) [19,28,29,[60][61][62][63][64][65][66][67][68][69][70][71][72][73] has attracted a lot of attention, since it produces a nearly scale invariant power spectrum. In the matter bounce theory, the Universe evolved from an epoch at large negative time in the contracting era where the primordial spacetime perturbations are generated deeply inside the Hubble radius, and is thus able to solve the horizon problem. After bouncing, the Universe starts to expand with an evolution that is symmetric to the contraction phase. Here, it may be mentioned that another popular bounce scenario is the Big-Bounce one, which is also able to give a nearly scale invariant power spectrum. However in this last scenario, the perturbations are produced near the bouncing point, unlike what happens in the matter bounce model, where the perturbations originate, as mentioned earlier, deeply in the contracting phase. This is the scenario we will consider in the present paper, in particular the extended matter bounce scenario about which we will discuss later.
Despite the mentioned successes, the MBS still faces some problems, firstly, in an exact MBS characterized by a single scalar field, the power spectrum turns out to be exactly scale invariant (i.e the spectral index of the curvature perturbation is exactly equal to one), what does not match the observational constraints. Such inconsistency was also confirmed in [70] from a slightly different viewpoint, namely from an F (R) gravity theory. It turns out that F (R) models can be equivalently mapped to scalar-tensor ones via conformal transformation of the metric [11,12,74] and, thus, the inconsistencies of the spectral index in the two different models are well justified. Secondly, according to the Planck 2018 data, the running of the spectral index is constrained to be −0.0085 ± 0.0073. However, for the MBS in the case of a single scalar field model, the running of the index becomes zero and hence it is not compatible with observations. A zero running index is a consequence of the first problem. If the spectral index becomes exactly scale invariant, then the variation of the spectral index with respect to the mode momentum, i.e the running index, will obviously vanish. Thirdly, in the simplest MBS model, the amplitude of scalar fluctuations is found to be comparable to that of tensor perturbations, which in turn makes the value of the tensor-to-scalar ratio to be of order one, again in conflict with the Planck constraints. Lastly, the null energy condition in most of the bouncing models is violated when the bouncing is realized.
In order to solve all these issues, theories of modified gravity [75,76] have been added to the picture, with success in many of the cases. For instance, in the so-called quasi-matter bounce scenario (which improves the exact matter bounce one), according to which the scale factor of the Universe evolves as t 3(1+w) (with w = 0), deeply into the contracting era, it is possible to recover the consistency of the spectral index and of the running index, even in a single scalar field model [19]. However, the tensor-to-scalar ratio still remains a problem. From a different perspective, it has recently been shown [20] that an F (R) gravity model with Lagrange multiplier is able to resolve most of the problems arising in the context of matter or quasi-matter bounce scenarios, albeit it fails to restore the energy conditions. In such Lagrange multiplier F (R) gravity model, the tensor power spectrum becomes smaller in comparison with that of scalar perturbation and thus the tensor to scalar ratio lies within Planck constraints, unlike in pure F (R) gravity, where the observable parameters are not simultaneously compatible with the observational constraints.
The recently observed striking event of neutron star merging GW170817 [77], reflects the fact that the propagation speed of the gravitational and electromagnetic waves are the same, i.e. equal to one in natural units. This observation has definitely narrowed down the viability of the modified gravitational theories, since every theory that predicts a gravitational wave speed different from one is no more a viable description (see Ref. [ [78]] for a list of theories that are ruled out by [77]).
Motivated by the above arguments, we will here study an extended matter bounce scenario in the context of a ghost free f (R, G) model, as developed in [1], and the conditions for such model to have gravitational waves speed (c 2 T ) equal to one and thus become consistent with GW170817. As it was demonstrated in [1], the constraint c 2 T = 1 restricts the form of the Gauss-Bonnet coupling function. After finding the suitable form which also realizes the considered bounce model, we address the bouncing phenomenology of the resulting theory.
The paper is organized as follows. In Sect.
[II], we briefly discuss the essential features of the ghost free f (R, G) model compatible with the event GW170817. Sections [III] and [IV] are respectively devoted to an extensive study of the primordial scalar and tensor perturbations and the observable parameters. Finally, the conclusions and remarks follow in the end of the paper.
In this section we shall recall the essential features of the ghost free f (R, G) gravity theory developed in Ref. [1]. We consider f (R, G) = R 2κ 2 + f (G) which, owing to the presence of f (G), contains ghosts with respect to perturbations of the spacetime metric. However the ghost modes may be eliminated by introducing a Lagrange multiplier λ in the standard f (G) gravity action [1], leading to a ghost-free action, as follows where µ is a constant with mass dimension [+1]. Upon variation with respect to the Lagrange multiplier λ, one obtains the following constraint equation Effectively, the kinetic term is a constant, so it can be absorbed in the scalar potential, as and the action of Eq. (1) can be rewritten as The scalar and gravitational equations of motion for the action (4) take the form where D τ η µν has the following form, with having in mind g µν D τ η µν = 4 − 1 2 g τ η R + R τ η . Upon multiplication of Eq. (6) with g µν , we get and solving Eq. (7) with respect to λ yields Eventually, we will deal with non-singular bounce cosmology in this ghost free f (R, G) model and, thus, the spatially flat Friedmann-Robertson-Walker (FRW) metric ansatz will fulfill our purpose. Let us now see how the equations of motion become if the metric background is a flat FRW one, with line element Assuming that the functions λ and χ are only cosmic time dependent, and also that no matter fluids are present, that is, that T matter µν = 0, then Eq. (2) admits the following simple solution Hence, the (t, t) and (i, j) components of Eq. (6) can be written as and, in addition, from Eq. (5) we get By solving Eq. (11) with respect to λ, we obtain It is easy to see that, by combining Eqs. (14) and (13), we obtain Eq. (12). Also by solving Eq. (12) with respect to the scalar potentialṼ µ 2 t , we get Hence, for an arbitrarily chosen function h(χ(t)), and with the potentialṼ (χ) being equal tõ then we can realize an arbitrary cosmology corresponding to a given Hubble rate H(t). Finally, the functional form of the Lagrange multiplier reads As mentioned in the introductory section, here we are interested on bouncing scenario followed from the theory with Lagrangian (4). The cosmological perturbation, in particular the tensor perturbation caused by the quantum vacuum fluctuation, indicates the gravitational wave. On other hand, recently the event GW170817 confirms the detection of the gravitational wave which may has a different source than the quantum vacuum fluctuation, and ensures that the speed of the gravitational wave is equal to the speed of light (in natural units, it is equal to unity). In the current paper, we consider that all the gravitational waves generated from different sources (i.e from black hole merging or from quantum vacuum fluctuations) have same propagation speed, and due to GW170817, the speed is equal to unity. With this consideration, we apply the result of GW170817 in the present context of bouncing scenario. Earlier, the implications of GW170817 on non-singular bounce has been investigated, but in a different context i.e from a degenerate higher order scalar-tensor (DHOST) theory [79], unlike to the present paper where we discuss the bouncing scenario from ghost free f (R, G) gravity model taking the constraint of GW170817 into account. The resulting theory with Lagrangian (4) is a form of the scalar Einstein-Gauss-Bonnet gravity along with a Lagrange multiplier, for which it is well known that the speed of gravitational waves is different from one, in particular, the gravitational wave speed in the present context can be expressed as [80][81][82], with H =ȧ a being the Hubble parameter. Eq.(18) apparently reflects the non-viability of the model with respect to GW170817 (which validates the fact that the gravitational and electromagnetic waves have the same propagation speed). However the gravitational wave speed in the ghost free f (R, G) model becomes one if the coupling function satisfies the following constraint equation [83],ḧ =ḣH.
Thereby, in order to make compatible our present model with GW170817, we consider such Gauss-Bonnet coupling functions which obey Eq. (19). Thus, a form of the Hubble parameter (later, we will consider the Hubble parameter in such a way that it gives an extended scenario of matter bounce cosmology) fixes in turn the coupling function via Eq. (19). At this stage it is worth mentioning that the new constraint on h(χ) also fits with the original equations of motion, this being clear from the fact that there are two independent equations, namely the (t, t) component of the gravitational equation and the equation for χ, and yet two unknown functions (λ(t), V (χ)) to determine. In the later section, we will consider a Hubble parameter which leads to an an extended matter bounce cosmology. In this context, it is worth mentioning that the presence of the potential termṼ (χ) in f (R, G) model to realize a non-singular bounce cosmology which is also consistent with c 2 T = 1. Let us demonstrate this in more detail by examining the equations of motion without the potential term. ForṼ = 0, eqn.(12) takes the following form- Recall from eqn. (19), the Gauss-Bonnet coupling satisfies a constraint equation in order to make the propagation speed of the gravitational wave unity. On applying this constraint in eqn. (20), we get, which immediately leads to the solution as -2Ḣ + 3H 2 = 0, OR ḧ + 1 8κ 2 = 0. The first possibility leads to the evolution of the scale factor as a(t) = t 2/3 i.e an eternal matter dominated universe, On other hand, the second possibility givesḧ = − 1 8κ 2 which along with the constraintḧ =ḣH makes the scale factor evolution as a(t) ∼ t i.e a Milne universe. Thus eqn. (21) is unable to realize a non-singular bounce cosmology. This indicates that in absence of potential, f (R, G) model can not realize a non-singular bounce cosmology which is also consistent with c 2 T = 1. However, in the present paper, as mentioned earlier, we will explore a non-singular bounce scenario in f (R, G) model compatible with GW170817 and thus the consideration of the potential term becomes obvious from its own right. In the next section, we shall extensively discuss the bouncing scenario of this f (R, G) model, which is both ghost free and compatible with the GW170817 event.

III. REALIZATION OF THE EXTENDED MATTER BOUNCE SCENARIO
In this section we will discover which functional forms for h(χ) and λ(t) can realize a bouncing Universe cosmological scenario with the following scale factor where a 0 and n are the free parameters of the model, a 0 having mass dimension [+2], while n is dimensionless. The above scale factor leads to a matter like bouncing universe. Here it may be mentioned that the scale factor in eqn. (22) is not the unique bouncing model that our universe may realize. There are other possibilities of bounce, as for example -(1) the exponential bouncing model where the scale factor evolves as a(t) = e αt 2 and the bounce occurs at t = 0, (2) the singular bounce scenario where the scale facor has the form a(t) = e 1 α+1 (t−ts) α+1 [36]. In the case of singular bounce, the perturbation mode(s) generate near the time of bouncing unlike to the usual bouncing model where the primordial perturbations generate in the past contracting phase deep inside the Hubble radius. However, as mentioned in the introductory section, the matter-like bounce scenario demonstrated by the scale factor in eqn. (22), faces some questions regarding the observational compatibility with Planck constraints, the instability of linear order perturbation theory in super horizon scale etc. Motivated by such questions, we investigate the matterlike bounce scenario in a modified gravity theory like f (R, G) model. Thus the calculations in the following sections of the paper are based on this choice of the extended matter bounce scale factor. In this sense, one may argue that the results we will obtain depend on the ansatz of the scale factor shown in eqn. (22). The Universe's evolution in a general bouncing cosmology, consists of two eras, an era of contraction and an era of expansion. The above scale factor describes a contracting era for the Universe, when t → −∞, then the Universe reaches a bouncing point, at t = 0, at which the Universe has a minimal size, and then the Universe starts to expand again, for cosmic times t > 0. Hence, the Universe in this scenario never develops a crushing type Big Bang singularity. It may be mentioned that for n = 1/3, the scale factor describes a matter bounce scenario. Eq. (22) leads to the following Hubble rate and its first derivative With the help of the above expressions, the Ricci scalar is found to be, Using Eq. (24), one can determine the cosmic time as a function of the Ricci scalar, that is the function t = t(R). As a result, the Hubble rate and its first derivative can be expressed in terms of R (this statement holds true for all analytic functions of t) and, also, the differential operator d dt can be written as d dt =Ṙ(R) d dR . However, for the purpose of determining t = t(R) as well as the observable quantities, we will consider the low curvature limit of the theory. Before proceeding further, let us comment on the viability of this approximation. We do it in the context of matter bounce cosmology, which is obtained by taking n = 1/3 in Eq. (22), the primordial perturbations of the comoving curvature, which originate from quantum vacuum fluctuations. At subhorizon scales during the contracting era in the low-curvature regime, their wavelength was much smaller than the comoving Hubble radius, which is defined by r h = 1 aH . In the matter bounce evolution, the Hubble horizon radius decreases in size, and this causes the perturbation modes to exit from the horizon, eventually, with this exit occurring when the contracting Hubble horizon becomes equal to the wavelength of these primordial modes. However, in the present context, we consider a larger class of bouncing models of the form a(t) = (a 0 t 2 + 1) n , always within Lagrange multiplier ghost free Gauss-Bonnet gravity. Thus, it will be important to check what are the possible values of n which make the low-curvature limit, that is, R/a 0 ≪ 1 a viable approximation in calculating the power spectrum for the bouncing model a(t) = (a 0 t 2 + 1) n . This expression of the scale factor immediately leads to the comoving Hubble radius Thereby r h diverges at t ≃ 0, as expected because the Hubble rate goes to zero at the bouncing point. Furthermore, the asymptotic behavior of r h is given by r h ∼ t 1−2n , thus r h (|t| → ∞) diverges for n < 1/2, otherwise r h goes to zero asymptotically. Hence, for n < 1/2, the comoving Hubble radius decreases initially in the contracting era and then diverges near the bouncing point; unlike in the case n > 1/2, where the Hubble radius increases from past infinity and gradually diverges at t = 0. As a result, the possible range of n which leads the perturbation modes to exit the horizon at large negative time and make the low-curvature limit a viable approximation in calculating the power spectrum, is given by 0 < n < 1/2, see Fig[1], where it is shown that for n = 0.30, the perturbation modes can generate inside the Hubble horizon at a large negative time. Moreover, we will show in the later sections that this range of n makes the observable quantities compatible with the Planck constraints and, thus, the "low curvature limit" comes as a viable approximation to calculate the power spectra of scalar and tensor perturbations.
In Fig[1], we show the spacetime sketch of our non-singular bouncing cosmology where the bounce time is taken at t = 0. We also plot the evolution of the physical length corresponding to a fixed co-moving scale. This scale is the wavelength of the fluctuation mode k, with k standing for the co-moving wavenumber. The wavelength begins at large negative time in the contracting phase on sub-Hubble scale, exits the Hubble radius during this phase at a time (say −t h (k)), and re-enters the Hubble radius during the low curvature regime in expanding phase at the time t h (k) (the exit and entry time are symmetric about the bouncing point as the scale factor is itself symmetric) and make the present observation useful. Thus the spacetime perturbation (or more explicitly the comoving curvature perturbation) has to be evolved from the contracting phase to the expanding one, followed by the bouncing phase, in order to get the power spectrum at late time. In the large scale limit (i.e in the super-Hubble scale k ≪ aH) of the contracting phase, the comoving curvature perturbation (ℜ(k, τ )) satisfies the cosmological perturbation equation where τ is the conformal time defined as dt = a(t)dτ and prime denotes the differentiation with respect to τ throughout the paper. The above equation is written in terms of the canonical variable : v(k, τ ) = zℜ(k, τ ), and the variable z(τ ) depends on the specific model, as for example -in scalar-tensor theory z(τ ) = aφ ′ H where φ is the scalar field and for our present considered model (i.e the Lagrange multiplier ghost free Gauss-Bonnet theory of gravity) the explicit form of z(τ ) is shown in the later section. However in terms of a general z(τ ), the solution of Eq.(26) is given by, where the suffix 'c' denotes the contracting phase and D c (k), S c (k) are independent of time and carry the information about the spectra of the two modes. The above solution of v(k, τ ) immediately leads to the curvature perturbation in the super-Hubble scale of the contracting phase as, As evident from the above expression that the D mode is a constant mode and generally the S mode appears as an increasing mode. Similarly in the large scale limit of the expanding phase, the curvature perturbation has the following solution, The D e mode of the curvature perturbation is constant in time, as is the D c mode in contracting phase. However, the role of the S mode is very different. In the expanding phase S e is the sub-dominant decreasing mode, whereas in the contracting phase, the behavior gets changed i.e S c mode is the decreasing mode. Therefore, the dominant mode of the curvature perturbation in the period of expansion is D e . This leads to the power spectrum of the curvature perturbation at late time (which is useful for the present observation) as At this stage it may be mentioned that a model will be a viable one if the power spectrum P e (k) becomes nearly scale invariant accordance to the observations of Planck 2018. Moreover the curvature perturbation should be continuous and thus ℜ c , ℜ e have to be matched through the bouncing point as explicitly performed in [62,63]. During this matching procedure, the D e mode may inherit the contribution from both D c and S c modes [62,63]. However as shown in [62], a nearly scale invariant power spectrum of the constant mode in the contracting phase (P c (k) ∼ k 3 |D c | 2 ) will eventually lead to a scale invariant power spectrum in the expanding phase. So in the present paper, we consider the bouncing scenario in Lagrange multiplier Gauss-Bonnet gravity model and concentrate in calculating the power spectrum for the constant mode in the contracting phase of the Universe, which will be discussed in detail in the next section.
Back to our bouncing model a(t) = (a 0 t 2 + 1) n , the Ricci scalar (R(t)) during the low-curvature regime (or at large negative time) can be written as R(t) ∼ 12n(4n−1) from Eq. (24). This helps to express the scale factor, the Hubble rate and its first derivative in terms of the Ricci scalar R, as follows the '+' and '-' signs in the expression of H(R) indicate the Hubble parameter in the expanding and contracting phases, respectively. Moreover, by using H(t) = 2n/t (which is valid in the low curvature regime), we immediately get the coupling function h(χ(t)) from Eq. (19), as where h 0 is an integration constant having mass dimension [2n+1]. Further, from the expression R = 1/t 2 , we obtain the coupling function in terms of the Ricci scalar, as Again the '+' and '-' signs indicate h(R) in the expanding and contracting phases, respectively. Plugging back these expressions of H(R) and h(R) into Eq. (17), we obtain the explicit form of the Lagrange multiplier λ(R) in the low curvature regime of the contracting phase of the Universe, as Eqs. (31), (33), and (34), which are all valid in the low curvature approximation, i.e. for R a0 ≪ 1, are the main ingredients to determine the observable quantities (recall that the observable parameters are eventually evaluated at the time of horizon exit, which in turn occurs at a large negative value of time in the present context) for our considered f (R, G) model with the scale factor depicted in Eq. (22).
This will be the subject of the next section. Later we will determine the exact expressions of the Hubble parameter, coupling function. and Lagrange multiplier, which are valid for all cosmic time beyond the low curvature approximation.

IV. COSMOLOGICAL PHENOMENOLOGY AND VIABILITY OF THE MODEL
In this section we study the first order metric perturbations of the theory, following Refs. [80][81][82], and where the scalar and tensor perturbations are calculated for various variants of the higher curvature models. Scalar, vector and tensor perturbations are decoupled, as in general relativity, so we can focus our attention to tensor and scalar perturbations separately. However before starting the cosmological perturbation calculations, we first determine the functions Q [80][81][82] in the context of the ghost free f (R, G) gravity model, which will be useful later. Such functions are defined as where we use the forms of h(t) and H(t) that we have derived previously. Moreover the function Q f becomes zero due to the fact that the coupling function h(χ) obeys the constraint Eq. (19).

A. Scalar perturbations
The scalar perturbation of FRW background metric is defined as follows, where Ψ(t, x) denotes the scalar perturbation. In principle, perturbations should always be expressed in terms of gauge invariant quantities, in our case the comoving curvature perturbation defined as ℜ = Ψ − aHv, where, v(t, x) is the velocity perturbation. However, we shall work in the comoving gauge, where the velocity perturbation is taken as zero, thus with such gauge fixing ℜ = Ψ. Thereby, we can work with the perturbed variable Ψ(t, x). The perturbed action up to Ψ 2 order is [80], where z(t) and c 2 s (the speed of the scalar perturbation wave) have the following expressions, and respectively, where F = 1 in our case. For more details on this we refer the reader to [80]. The definition of the wave speed is for the general Gauss-Bonnet corrected theory with F = ∂f ∂R , but in our case f = R and F = 1. Also the waves speed is affected from the Gauss-Bonnet coupling via the functions Q f and Q b which in our case have the form (42). With the above expressions, we concentrate on determining various observable quantities and specifically, the spectral index of the primordial curvature perturbations, the tensor-to-scalar ratio and the running of the spectral index, which are eventually determined at the time of horizon exit. For the scale factor we consider in the present paper, the horizon exit occurs during the low-curvature regime deeply in the contracting era. Thereby, for the purpose of finding the observable parameters, the condition R/a 0 ≪ 1 stands as a viable approximation.
In the low-curvature limit, we determine various terms present in the expressions of z(t) and c 2 s (see Eqs. (38) and (39)) as, . Consequently z(t) takes the following form, where P (R) and Q(R) are defined as follows, and Moreover c 2 s has the following form Therefore up to first order inR 1/2−n , the speed of the scalar perturbation wave is unity, however here we retain the terms up to the order of R 1−2n and thus c 2 s becomes different than unity in the low curvature regime of our Universe.
Moreover at the bouncing point, the functions Q a = Q b = 0 (as these functions are proportional to the Hubble parameter, see Eq. (35)) which lead to the speed of the scalar perturbation wave c 2 s = 1 at the bouncing phase. The positivity of c 2 s stops the exponential growth of the scalar perturbation near the bouncing regime, as the squeezing term becomes negligible with respect to c 2 s k 2 in the cosmological perturbation equation, and makes the first order perturbation calculation reliable. Eq. (37) clearly indicates that Ψ(t, x) is not canonically normalized and to this end we introduce the well-known variable as v = zℜ (= zΨ as we are working in the comoving gauge). The corresponding fourier mode of the variable z satisfies, where τ = dt/a(t) is the conformal time and v k (τ ) is the Fourier transformed variable of v(t, x) for the kth mode. Eq. (44) is quite hard to solve analytically in general, since the function z depends on the background dynamics and also c 2 s is not constant. However the equation can be solved analytically at super-Hubble scale (i.e when c 2 s k 2 can be neglected with respect to the squeezing term z ′′ /z) in the regime R/a 0 ≪ 1 as we now show. The conformal time (τ ) is related to the cosmic time (t) as τ = dt a(t) = 1 a n 0 (1−2n) t 1−2n for n = 1/2, however we will show that the observable quantities are compatible with Planck data [84] for n < 1/2 and thus we can safely work with the aforementioned expression of τ = τ (t). Using this, we can express the Ricci scalar as a function of the conformal time, Having this in mind, along with Eq. (40), we can express z in terms of τ as follows, The above expression of z = z(τ ) yields the expression of 1 z d 2 z dτ 2 , which is essential for the cosmological perturbation equation, with ξ = (2n) (1−2n) . In the above expression, we retain the terms up to (R/a 0 ) 1−2n as the cosmological parameters are eventually determined in the low curvature regime, as discussed earlier. Moreover without any loss of generality, the integration constant h 0 is replaced by h 0 =h 0a n 0 (2n+1) and κ 2h 0 is taken as κ 2h 0 = 1/ √ a 0 . Recall, the low curvature limit is valid for n < 1/2 which clearly indicate that 1 2 − n (or 1 − 2n) is a positive quantity. Thus the term within parenthesis in Eq. (47) can be safely considered to be small in the low-curvature regime R/a 0 ≪ 1. As a result, 1 z d 2 z dτ 2 becomes proportional to 1/τ 2 i.e., 1 z d 2 z which is approximately a constant in the era, when the primordial perturbation modes are generated deeply inside the Hubble radius. In effect, the cosmological perturbation equation in super-Hubble scale can be solved as follows, with ν = σ + 1 4 and the k-dependent multiplication factor is fixed by assuming the Bunch-Davies vacuum initially. Having the solution of v k (τ ) at hand, next we proceed to evaluate the power spectrum (defined for the Bunch-Davies vacuum state) corresponding to the k-th scalar perturbation mode, which is defined as follows, In the superhorizon limit, using the mode solution in Eq. (49), we have, At this stage it deserves mentioning that in calculating the power spectrum P Ψ (k, τ ), we use the linear order perturbation theory, which is reflected through the fact that the perturbed action is taken up to the quadratic order of the perturbed variable. However it is important to investigate the validity of the linear order perturbation in our considered f (R, G) model. For this purpose, we use the super-horizon solution of v(k, τ ) from eqn. (49) and determine the comoving curvature perturbation in super-Hubble scale as, where z(τ ) is given in eqn. (46). Plugging this expression of z(τ ) in eqn. (52), we get Using eqns. (41) and (42), we determine the following expressions of P (R) and Q(R) (in the limit κ 2 R < 1) as, and respectively. Plugging back the above expressions of P (R) and Q(R) into eqn.(53) yields the perturbed variable Ψ(k, τ ) as, where we use R(τ ) = 12n(4n−1) [a n 0 (1−2n)] 2/(1−2n) 1 τ 2/(1−2n) from eqn. (45) and Z(R) has the following form, The linear order perturbation theory is valid as long as k 3/2 Ψ(k, τ ) is less than unity i.e In order to find an explicit condition from the above inequality, we consider the exact matter bounce scenario (i.e n = 1/3, later we will show that the exact MBS is consistent with the Planck constraints in our considered model) and put ν = 3/2. With these values of n and ν, eqn.(58) becomes, where, as previous, h 0 is replaced by h 0 =h 0a n 0 (2n+1) and κ 2h 0 is taken as κ 2h 0 = 1/ √ a 0 . However the inequality in eqn. (59) is valid in the regime R > 0 i.e for entire cosmic time. This indicates that the linear order perturbation theory is valid in the entire range of cosmic time in the ghost free F (R, G) bouncing model. However in a scalar-tensor matter bounce model, P (R) and Q(R) take the form as P = 4/3 and Q = 2/3 respectively, actually the dependence of Ricci scalar in P (R) and Q(R) (see eqns. (54) and (55)) arise due to the presence of f (G) gravity in the present model. Therefore, in scalar-tensor matter bounce model, the inequality k 3/2 Ψ(k, τ ) < 1 leads to the following condition (in terms of conformal time) : Eqn. (60) clearly indicates that in a scalar-tensor matter bounce model, the perturbation is less than unity and makes the linear order perturbation valid in the regime |τ | > τ f . Thereby the perturbation theory in a scalar-tensor matter bounce scenario is ill-defined to address the evolution of the perturbed variable in the entire range of cosmic time. This is in contrary with our considered f (R, G) theory where the linear order perturbation theory is valid for the entire range of cosmic time, as explained earlier.
Actually the presence of f (G) term in the gravitational action modifies the solution of the comoving curvature perturbation, which in turn validates the linear order perturbation theory for entire cosmic time. By using Eq. (51), we can determine the observable quantities like spectral index of the primordial curvature perturbations and the running of spectral index. Before proceeding to calculate these observable quantities, we will consider first the tensor power spectrum, which is necessary for evaluating the tensor-to-scalar ratio.

B. Tensor Perturbations
Let us now focus on the tensor perturbations, which for the FRW metric background are defined as where h ij (t, x) is the perturbation. The tensor perturbation is itself a gauge invariant quantity, and the tensor perturbed action up to quadratic order is given by where z T (t) and c 2 T (the speed of the gravitational wave) are Recall that the coupling function in the present context satisfies the relationḧ =ḣH, which in turn makes the speed of the gravitational wave to be one, i.e. c 2 T = 1, and thus the model becomes compatible with the GW170817 event. Similar to the case of scalar perturbations, the variable z for tensor ones is defined as (v T ) ij = z T h ij which, upon performing the Fourier transformation, satisfies the equation By using Eq. (63), along with the condition R/a 0 ≪ 1, we evaluate z T (τ ) and 1 zT (τ ) respectively, where S(R(τ )) = and recall ξ = (2n) (1−2n) . The above expressions yield the tensor power spectrum, defined with the initial state of the Bunch-Davies vacuum, so we have, The factor 2 arises from the two polarization modes of the gravity wave, and ν T = σ T + 1 4 , where σ T is defined in Eq. (67). Now we can explicitly confront the model at hand with the latest Planck observational data [84], so we shall calculate the spectral index of the primordial curvature perturbations n s and the tensor-to-scalar ratio r, which are defined as Eqs. (51) and (68) immediately lead to the explicitly form of n s and r where σ, z(τ ) and z T (τ ) are given in Eqs. (48), (45), and (65), respectively. As it is evident from the above equations, n s and r are evaluated at the time of horizon exit, when k = aH, or equivalently at τ = τ h . It should be noticed that n s and r depend on the dimensionless parameters R h a0 and n, with R h = R(τ h ). We can now directly confront the spectral index and the tensor-to-scalar ratio with the Planck 2018 results [84], which constrain the observational indices to be n s = 0.9649 ± 0.0042 , r < 0.064 .
For the model at hand, n s and r are within the Planck constraints for the following ranges of parameter values: 0.01 ≤ R h a0 ≤ 0.05 and 0.30 n 0.40, and this behavior is depicted in Fig. 2. The viable range of R h /a 0 is in agreement with the low-curvature condition R/a 0 ≪ 1 that we have considered in our calculations. Moreover, the range of the parameter n clearly indicates that the matter bounce scenario, for which n = 1/3, is well described by the ghost free Einstein-Gauss-Bonnet gravity model with the Lagrange multiplier term. At this stage it is worth mentioning that Eq.(48) clearly reveals that in absence of the Gauss-Bonnet coupling term (i.e for h 0 = 0, for which the present model resembles a scalar-tensor model), the quantity σ becomes σ = 2 for the pure matter bounce scenario, which in turn makes the power spectrum completely scale invariant (i.e n s = 1). However, it is well known that in scalar-tensor theory the matter bounce scenario is not consistent with the Planck observations. Moreover the matter bounce scenario also does not fit well even in the standard F (R) gravity, as it was confirmed in a previous paper [20]. However, here we show that in the ghost free f (R, G) model which is also compatible with GW170817, the matter bounce can be considered as a good bouncing model, which allows the simultaneous compatibility of n s and r with the astronomical observations. Furthermore, the running of the spectral index is defined as follows, and this is constrained by Planck 2018 results as α = −0.0085 ± 0.0073. Thus, it is also important to calculate the running of spectral index before concluding the viability of a model. By using the expression of σ (see Eq. (48)) and R = R(τ ) (see Eq. (45)), we get To arrive to the above result, we use the horizon crossing relation of the k-th mode k = aH to determine d|τ | d ln k = −|τ |, i.e., the horizon exit time |τ | increases as the momentum of the perturbation mode decreases, as expected. Eq. (73) indicates that similarly to n s and r, the running index (α) also depends on the parameters R h /a 0 and n. Taking n = 1/3, in Fig. 3 we give a plot of α with respect to R h /a 0 . As it can be seen in Fig. 3, the parameter α lies within the Planck constraint for 0.01 R h /a 0 0.05, for n = 1/3. For the ghost free f (R, G) model, we have shown that both the pure matter bounce scenario as well as the quasi-matter bounce scenario are consistent with the Planck observations, what is not true in the case of scalar-tensor or pure F (R) models. Therefore, the Lagrange multiplier ghost free Gauss-Bonnet gravity model has a richer phenomenology as compared to the scalar-tensor or standard F (R) gravity models, which fail to describe in a viable way these two bounce cosmology scenarios. Before proceeding further, it is worth mentioning that in [85,86], the authors explored how to build a viable nonsingular cosmological models within the framework of effective field theory (EFT). They showed that a healthy nonsingular bounce model can be achieved by introducing an effective operator of R (3) δg 00 , in particular, the effective action looks like where Λ(t), c(t), M 2 (t), m 3 (t) andm 4 (t) are various time dependent coefficients. Having this effective action in hand, we like to investigate whether the f (R, G) model that we consider in the present paper can be embedded in or have some differences in respect to this EFT action (S ef f ). The demonstration goes as follows: The first line in the expression (74) describes the background model and the rest is for perturbations. Therefore the background Friedmann equations are given by, Moreover the effective action (74) leads to the quadratic action for scalar and tensor perturbations as [85], and respectively, where Ψ and h ij are scalar and tensor perturbations respectively. Moreover c i (s) have the following form, Comparing Eqns. (11), (12) with Eqn. (75), one can argue that the background equations of our considered f (R, G) model can be embedded within that of the effective field theory (EFT) action (74), if the coefficients Λ(t) and c(t) are related with various functions of the f (R, G) model as, and respectively, recall that χ(t) = µ 2 t. Moreover comparing Eqns. (37) and (76), we may conclude that the scalar perturbation of our considered f (R, G) model is also consistent with that of the effective action (74), provided that the coefficients c i (t) are given by, where the functions Q i are shown in Eqn. (35). However eqns. (62) and (77) clearly reveal that the tensor perturbation of the present ghost free f (R, G) model is not consistent with that of the effective action (74), because of the time dependence of the factor z T in eqn. (62), unlike in the expression of eqn. (77). At this stage it deserves mentioning that "z T " factor of the tensor perturbation in an effective field theory may be made time dependent if one adds the terms like −m 2 3) to the effective action i.e if S ef f takes the following form [85], The above effective action leads to the quadratic action of tensor perturbation as follows, It is evident from eqn.(83) that due to the time dependence of m 2 4 (t), the front factor of δS h (apart from a 3 (t)) becomes time dependent. Therefore the addition of the new term like −m 2 4 (t) δK 2 − δK µν δK µν in the effective action makes the "z T " factor of the tensor perturbation time dependent, which is consistent with our considered f (R, G) model. However eqn. (83) also indicates that due to the presence of the term −m 2 4 (t) δK 2 − δK µν δK µν , the gravitational wave speed in the EFT becomes different than unity, unlike to the current f (R, G) model where the gravitational wave speed is unity. This in turn clarifies that the present ghost free f (R, G) model is not consistent to the effective action (82). Thus the current Lagrange multiplier Gauss-Bonnet theory of gravity can neither be embedded within the EFT action (74) nor within the action (82). However with a non-trivial coefficient (other than unity) of Ricci scalar along with m 4 = 0 in the effective field theory action, in particular, if the EFT action has the following form [85], then the scalar perturbed action will remain the same as in Eqn. (76), however the tensor perturbed action comes as, Therefore the effective action in eqn. (84) makes the front factor of δS h (apart from a 3 (t)) time dependent and also allow an unit gravitational wave speed, which are consistent to the present Lagrange multiplier Gauss-Bonnet theory of gravity. Thus the model considered in the present paper is consistent with the effective action (84) if the parameters Λ(t), c(t) and f (t) satisfy the following expressions: and where Q b is shown in Eq. (35). The first two equations (i.e Eqs. (86) and (87)) are due to matching the background equations of motion where as Eq.(88) are due to matching the tensor perturbation equation of the EFT (described by action (84)) with that of the Lagrange multiplier Gauss-Bonnet theory of gravity.

C. The Energy Conditions Violation Issue for the Bounce
At this stage it deserves mentioning that apart from the computation of spacetime perturbation, the investigation of energy conditions is also important in a bouncing model from its own right as most of the bouncing models fail to rescue the null energy condition. As for example, the energy condition is violated in the bouncing scenario of Lagrange multiplier F (R) gravity model [20], however the holonomy generalization of such F(R) model is able to rescue the null energy condition [87]. Here, we want to investigate the energy conditions in the present context i.e for the bouncing scenario in the Lagrange multiplier Gauss-Bonnet theory of gravity. As we already mentioned, a crucial drawback in most of the bouncing models is the violation of the null energy condition. Here we will check the fulfillment of the energy conditions in the context of the ghost free f (R, G) gravity model. For this purpose, we first determine the effective energy density ρ eff + p eff from Eqs. (11) and (12), However, as we mentioned earlier, the energy condition has to be checked for all cosmic times, including the bouncing point which occurs at t = 0, where the low-curvature approximation no longer holds true. Thus, it will not be justified if we use the forms for the coupling function and Lagrange multiplier obtained in Eqs. (32) and (34), to check the energy condition near the bouncing point. Thereby, the best way to investigate the energy condition is to determine the forms of h(χ(t)) and λ(t) for the whole range of time −∞ < t < ∞, and then use such forms of h(t) and λ(t) in the expression of ρ eff + p eff . For this purpose, we reconstruct the GB coupling function and the Lagrange multiplier term from Eqs. (19) and (17), respectively, by using the exact form of H(t) = 2nt t 2 +1/a0 (which is valid for the entire range of the cosmic time)ḣ (t) =h 0 a 0 t 2 + 1 n µ 4 λ(t) = 4a 0 n κ 2 (1 − a 0 t 2 ) 1 + a 0 t t + 16nh 0 κ 2 (t 2 + 1/a 0 ) n 1 + a 0 t 2 3 , whereh 0 is an integration constant with mass dimension [+1], and we should recall thath 0 is related to h 0 by h 0 =h 0 a n 0 (2n+1) . It should be observed from the above equation that we stop at the level of the first derivative of the coupling function and we do not determine the explicit form of h(t). However, this approximation should be sufficient, because the right hand side of Eq.(89) does not contain any term that is without a derivative of h(t), which is also expected from the fact that the Gauss-Bonnet term without the coupling function becomes a surface integral and is eliminated from the action. Plugging back the expressions of h(t) and λ(t) into Eq.(89), in Fig. 4 we give the plot of ρ eff + p eff (with respect to the cosmic time), for n = 1/3,h 0 = 1, a 0 = 1 κ 2 = 1 (in reduced Planck units). As can be seen in Fig. 4, ρ eff + p eff becomes negative near the bouncing point. Furthermore, Eq.(89) clearly indicates that both the scalar field (χ) and the Gauss-Bonnet term contribute to ρ eff + p eff , their individual contributions being and respectively. Using these expressions we give the plot of ρ + p χ , ρ + p GB and ρ eff + p eff (with respect to the cosmic time, for n = 1/3,h 0 = 1, a 0 = 1 κ 2 = 1 i.e in reduced Planck units), see Fig. [5], in order to better understand the energy flow of the universe. Figs. [4] and [5] show that the null energy condition for the whole universe, as well as for the individual contributors, are violated, which further implies that the weak energy condition is necessarily violated. At this stage, we want to mention that the holonomy improvement is able to rescue the energy condition in the context of Lagrange multiplier F(R) gravity, as proven in [87]. Hopefully, therefore, the presence of holonomy modifications even in the f (R, G) gravity model, or the presence of extra spatial dimensions, where H 2 is proportional to linear powers as well as quadratic powers of the energy density, may play a significant role to rescue the null energy condition for a non-singular bounce. This investigation is expected to be carried out soon in a future work.

V. CONCLUSIONS
In the present paper we have considered an extended matter bounce scenario in a ghost free f (R, G) gravity, in particular, for the f (R, G) = R 2κ 2 + f (G) model. The idea for making the model ghost free consists in introducing a Lagrange multiplier, as discussed in [1]. In such gravity theory, we used the results of [83] which indicated that the model at hand is compatible with GW170817 that this happens for a class of Gauss-Bonnet (GB) coupling functions (h(t)) which satisfy the constraint equationḧ =ḣH, H(t) being the Hubble parameter. Thus, in order to make our model compatible with the event GW170817, we considered only such GB coupling functions which obey this constraint. At this stage, it is worth mentioning that this new constraint on the coupling function also fits with the equations of motion due to the presence of the scalar field potential.
In such ghost free f (R, G) model compatible with GW170817, we considered a non-singular bounce scenario with the scale factor being expressed as a(t) = (a 0 t 2 + 1) n , where n is a dimensionless parameter of the model and t the cosmic time. For this scale factor, it was shown that, for n < 1/2, the spacetime perturbation modes are generated deeply in the contracting era, at large negative values of time, where the Ricci curvature is low as compared to the one in the near-bouncing era. This, in turn, makes the "low curvature limit" a viable approximation in calculating the observable quantities, which are eventually determined at the time of horizon exit. We have determined the forms of the coupling function (h(t)) and Lagrange multiplier (λ(t)) in the low curvature regime, by using a reconstruction technique for the present model, which realizes the bouncing with the aforementioned scale factor. Such forms of h(t) and λ(t) led to an explicit expression of the cosmological perturbation equation, by solving which we have determined the power spectra of the primordial perturbations and, correspondingly, have calculated the fundamental cosmological parameters, namely the spectral index of scalar perturbations, the tensor to scalar ratio, and the running spectral index, respectively.
Such observational indexes are found to depend on the dimensionless parameters of the model, as R h /a 0 (with R h the Ricci curvature at the time of horizon exit) and n. It turned out that the observable quantities are simultaneously compatible with the Planck 2018 constraints for the parametric range 0.01 R h a0 0.05 and 0.30 n 0.40. It should be noticed that this range of n is also supported by the range 0 < n < 1/2, which ensures the low-curvature approximation to be a perfectly reliable one for calculating the power spectra. The expression of the spectral index we obtained clearly demonstrates that, in the absence of the GB term (the model resembling a scalar tensor one), the power spectrum becomes completely scale invariant for n = 1/3. This was expected, as it is well known that for a scalar tensor theory, the spectral index becomes one in a pure matter bounce scenario, which in turn makes the spectrum completely scale invariant. We recall that the pure matter bounce scenario is not consistent with observations, even in pure vacuum F (R) gravity. Here, we showed that for a ghost free f (R, G) model compatible with the GW170817 event, the matter as well as quasi-matter bounces may be considered as good bouncing theories, which confirms the richer bouncing phenomenology of our general model here. We further determined the effective energy density and pressure in the present context, in order to investigate the energy conditions. The energy condition has to be checked for all cosmic times, including at the bouncing point, where the low-curvature approximation no longer holds true. Thus, it will not be justified if we use the forms of coupling function and the Lagrange multiplier obtained from the low curvature approximation to check the energy condition near the bouncing point. Keeping this in mind, we have reconstructed h(t) and λ(t) for the entire range of the cosmic time (by relaxing the low curvature approximation) and then used such forms of h(t), λ(t) to determine the effective energy density and pressure. As a consequence, we found that the null energy condition is violated near the bouncing point, which further implies that the weak energy condition is necessarily violated. At this stage, we want to mention that the holonomy improvement has been proven to rescue the energy conditions in the context of Lagrange multiplier F (R) gravity [87]. Thus, hopefully, the presence of holonomy modifications, even in the f (R, G) gravity model, may similarly play a significant role to rescue the null energy conditions for a non-singular bounce. This investigation should be interesting and insightful, and we expect to carry it out in future work, soon.