Viable Inflationary Models in a Ghost-free Gauss-Bonnet Theory of Gravity

In this work we investigate the inflationary phenomenological implications of a recently developed ghost-free Gauss-Bonnet theory of gravity. The resulting theory can be viewed as a scalar Einstein-Gauss-Bonnet theory of gravity, so by employing the formalism for cosmological perturbations for the latter theory, we calculate the slow-roll indices and the observational indices, and we compare these with the latest observational data. Due to the presence of a freely chosen function in the model, in principle any cosmological evolution can be realized, so we specify the Hubble rate and the freely chosen function and we examine the phenomenology of the model. Specifically we focus on de Sitter, quasi-de Sitter and a cosmological evolution in which the Hubble rate evolves exponentially, with the last two being more realistic choices for describing inflation. As we demonstrate, the ghost-free model can produce inflationary phenomenology compatible with the observational data. We also briefly address the stability of first order scalar and tensor cosmological perturbations, for the exponential Hubble rate, and as we demonstrate, stability is achieved for the same range of values of the free parameters that guarantee the phenomenological viability of the models.

In this work we investigate the inflationary phenomenological implications of a recently developed ghost-free Gauss-Bonnet theory of gravity. The resulting theory can be viewed as a scalar Einstein-Gauss-Bonnet theory of gravity, so by employing the formalism for cosmological perturbations for the latter theory, we calculate the slow-roll indices and the observational indices, and we compare these with the latest observational data. Due to the presence of a freely chosen function in the model, in principle any cosmological evolution can be realized, so we specify the Hubble rate and the freely chosen function and we examine the phenomenology of the model. Specifically we focus on de Sitter, quasi-de Sitter and a cosmological evolution in which the Hubble rate evolves exponentially, with the last two being more realistic choices for describing inflation. As we demonstrate, the ghostfree model can produce inflationary phenomenology compatible with the observational data. We also briefly address the stability of first order scalar and tensor cosmological perturbations, for the exponential Hubble rate, and as we demonstrate, stability is achieved for the same range of values of the free parameters that guarantee the phenomenological viability of the models. Nearly forty years ago, three of the major problems in contemporary cosmology, namely the Horizon Problem, the Flatness Problem and the Magnetic-Monopoles Problem, have been given a successful solution in the context of the inflationary scenario. This scenario was firstly proposed in Ref. [1] and was further developed in Ref. [2,3]. According to the inflationary scenario, merely fractions of seconds after the Big Bang, the spatial coordinates of the Universe expanded exponentially. An expansion of this sort is supposed to last from about 10 −36 sec to 10 −15 sec and the size of the Universe is increased by a factor of 10 26 . The nature of this scenario is rather bizarre for classical cosmology, since traditional Big Bang Friedmann-Robertson-Walker (FRW) models do not match the fast evolution of the Universe [4][5][6]. The first approximation is to consider the expansion as the de Sitter phase of the Universe. The standard approach to achieve the de Sitter inflationary phase in cosmology is to use scalar fields, and many of the initial models of inflation made use of the scalar field formalism. However, it is also possible to produce an inflationary phase of the Universe in the context of modified gravity, see Refs. [7][8][9][10][11][12][13] for reviews on this. In fact, the first model of f (R) which remains viable up to date is the Starobinsky model [14], and ever since many models have been developed in various forms of modified gravity [7][8][9][10][11][12][13]. In all the modified gravities the key element is that geometric terms are included in the gravitational Lagrangian, which are absent in the Einstein-Hilbert gravity. These terms may dominate the Universe's evolution at early times or even at late times. Such models may include additional curvature terms, namely the f (R) theories, torsional terms namely the teleparallel f (T ) theories, or the Gauss-Bonnet modified gravities f (G) theories, as well as the generalized f (R, G) theories (see [7][8][9][10][11][12][13]). Such theoretical formulations of gravity are able to model both the early-time expansion and In this section we shall recall the essential features of the ghost free f (G) gravity developed in Ref. [16]. The whole ghost-free construction scheme is based on introducing a Lagrange multiplier λ in the standard f (G) gravity action, so the ghost-free action is the following, where µ is a mass-dimension one constant. Upon variation with respect to the Lagrange multiplier λ, we obtain the following constraint equation, Effectively, the kinetic term is a constant, so it can be absorbed in the scalar potential in the following way, and in effect, the action of Eq. (1) is rewritten as, The equations of motion for the action (4), are (2) and the following, Upon multiplication of Eq. (6) with g µν , we get, By solving Eq. (7) with respect to λ, we get, Let us now see how the equations of motion become if the metric background is a flat Friedmann-Robertson-Walker (FRW), with line element, Assuming that the functions λ and χ are only cosmic time dependent, and also that no matter fluids are present, that is, T matter µν = 0, Eq. (2) has the following simple solution, Hence, the (t, t) and (i, j) components of Eq. (6) can be written, and in addition from Eq. (5) we get, 0 = µ 2λ + 3µ 2 Hλ + 24H 2 Ḣ + H 2 h ′ µ 2 t −Ṽ ′ µ 2 t .
By solving Eq. (11) with respect to λ we get, It is easy to see that by combining Eqs. (14) and (13), we easily 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 to, then we can realize an arbitrary cosmology corresponding to a given Hubble rate H(t). Finally, the functional form of the Lagrange multiplier is equal to, The resulting theory with Lagrangian (4) is a form of the scalar Einstein-Gauss-Bonnet gravity and in the next section we shall extensively discuss the inflationary dynamics of this model. The presence of the arbitrary function h(χ) provides us with the freedom of realizing several viable cosmologies.

III. INFLATIONARY DYNAMICS OF THE GHOST-FREE f (G) MODEL
As we already mentioned, the ghost-free f (G) model of Eq. (4) is a sort of scalar Einstein-Gauss-Bonnet model [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36], the cosmological perturbations of which were studied in Ref. [37]. In this section we shall use the formalism, notation and results of Ref. [37], and we shall calculate the spectral index of primordial curvature perturbations and the tensor-to-scalar ratio for the model (4), by specifying the functional form of h(χ) and the Hubble rate. Then, by replacing the cosmic time with the e-foldings number, we shall express all the observational and slow-roll indices as functions of the e-foldings number, and we shall put the phenomenology of the model into test by confronting the resulting theory with the latest observational data.
We begin by defining the functions Q i (χ) (see [37] for more details), as follows, where H is the Hubble rate, H ≡ȧ/a. In addition, the wave speeds c A and c T become, , ∂ 2 f ∂X 2 = 0 and F = 1 in our case. Note that c A is the wave speed of the perturbed field in the context of the perturbed FRW metric, and c T is the sound speed. For more details on this we refer the reader to [37]. The definition of the wave speeds is for the general Gauss-Bonnet corrected f (R, χ) theory with F = ∂f ∂R , but in our case f (R, χ) = R and F = 1. Also the waves speeds are affected from the Gauss-Bonnet coupling via the functions Q f and Q b which in our case have the form (18). As a result, the two wave speeds are further simplified with the wave speed of the perturbed field c A being trivial as in the classical case, while the wave speed of the gravitational waves in non-trivial, In order to calculate the slow-roll parameters, we first need to determine the function E(R, χ, X) which is defined as follows [37], The slow-roll parameters are defined as follows [37] (R, χ, X) HE(R, χ, X) , The two spectral indices, for scalar and for tensor perturbations in the inflationary era respectively, are defined using the slow-roll parameters [37], Finally, the tensor-to-scalar ratio is equal to [37], The above expressions of the parameters for the slow-roll inflationary dynamics, are in fact functions of the cosmic time, t. However, such a description is not sufficient for our study, since the preferable variable to perfectly quantify the evolution during the inflationary era is the e-foldings number, N . So we need to transform the above relations with respect to the e-foldings numbers. At first, we consider a given Hubble expansion rate for the inflationary era, as a function of time, H = H(t). The e-foldings number is defined as where t i is the initial and t f the final moments of inflation. Considering a given initial moment for inflation, t i ∈ [0, 10 −36 ], and an unspecified final moment, t, the e-foldings number is obtained via Eq. 25 as a function of time, N = N (t). Supposing this function is reversible, time is also given as a function of the e-foldings number, t = t(N ). Consequently, the first-and the second-order derivatives with respect to time, are transformed into first-and secondorder derivatives with respect to the e-foldings number, as follows, Since the scalar field, χ = χ(t) is a function of time, its potential,Ṽ (χ) =Ṽ (χ(t)), and the Lagrange multiplier, λ = λ(t), the coupling function, h(χ) = h(χ(t)), as well as the Ricci scalar, the Gauss-Bonnet invariant and the function E(R, χ) = E (R(t), χ(t)) are also functions of the cosmic time. As a result, they can all be rewritten with respect to the e-foldings number. Furthermore, the functions Q i (χ) are also transformed, taking the following forms, where the prime denotes differentiation with respect to the e-foldings number. In the same manner, we may redefine the wave speed for the gravitational waves, The next step is to express the slow-roll parameters, ǫ i , with respect to the e-foldings number, and the resulting expressions are, Through these, the spectral indices and the tensor-to-scalar ratio are directly calculated with respect to the e-foldings number, using Eqs. (23) and (24). What remains is to define a specific coupling function, h(χ), as well as the Hubble rate for the cosmological FRW background, and also to calculate the spectral indices and the tensor-to-scalar ratio and compare our results with that of the latest Planck [17] and BICEP2/Keck-Array [18] observations. With regard to the coupling function, we shall assume that it has either exponential or power-law forms, while with regard to the Hubble rate, we shall firstly assume the de Sitter evolution for a warm up study, and finally we shall assume the quasi-de Sitter evolution.

IV. THE CASE OF DE SITTER BACKGROUND EVOLUTION
In the de Sitter case, the Hubble rate is constant as a function of the cosmic time, therefore the e-foldings number and the cosmic time are related as follows, As a result, the Ricci scalar and the Gauss-Bonnet invariant are both constant, Finally, the scalar field given by Eq. (10), takes the following form, Using, Eqs. (30), (31) and (33) and in addition a specific form for the function h(χ), we can calculate the slow-roll indices and the observational indices for the de Sitter evolution cosmology.
A. A power-law coupling function, h(χ) = γχ b Let us assume that the coupling function is a simple power law, where γ and b are real constants, to be used as free parameters later. Using Eqs. (33) and (31), we can write the coupling function first as function of time, and then as a function of the e-foldings number, Using Eq. (16), we may derive the potential as a function of the e-foldings number, as well as the Lagrange multiplier, From the equations in (27), we can write the Q i functions with respect to the e-foldings number, as follows, while the wave-speeds appearing in Eqs. (19) and (28) are, The function E(R, χ) is written with respect to the e-foldings number in the following way, Using Eqs. (29), (30), (36) and (41), we obtain the slow-roll parameters of the de Sitter evolution case, which are, Using the above results, we can proceed in calculating the spectral indices, from Eqs. (23), and the tensor-to-scalar ratio, from Eq. (24), r = 16 Having these at hand, we can compare them directly to the Planck [17] and the BICEP2/Keck-Array data [18], which indicate that n S = 0.9649 ± 0.0042 and r < 0.064. It can be shown that the viability of the theory is achieved for a restricted range of values of the free parameters. Actually, if we set N = 50 (or N = 60) to indicate the end of the inflationary era, it is easy to see that the values of H 0 , γ and µ do not affect the resulting values. In effect, we choose γ = 1 and µ = 1 sec −1 for simplicity and H 0 = 10 26 sec −1 (or H 0 = 10 27 sec −1 ). The tensor-to-scalar ratio is constantly close to zero, while the spectral index coincides with the Planck data only for µ ∼ 4 sec −1 . Namely, n S = 0.9644 only for b = 3.78 when N = 50, or b = 4.136 for N = 60 for the same values, r ∈ [10 −50 , 10 −20 ]. In Fig. 1 we present the plots of the spectral index and of the tensor-to-scalar ratio as a function of b. As a result, a power-law coupling function for the de Sitter background evolution, may generate a viable inflationary model, only under the strict assumption of h(χ) ∼ χ 4 .
B. An Exponential Coupling Function, h(χ) = γe bχ In this case, we assume that the coupling function h(χ) has the following exponential form, where γ and b are real constants, to be used as free parameters later. Using Eqs. (33) and (31), we can write the coupling function first as function of time, and then as a function of the e-foldings number, At this point, by using Eq. (16), we may derive the potential as a function of the e-foldings number, as well as the Lagrange multiplier, Accordingly from Eqs. (27), we derive the Q i functions with respect to the e-foldings number, as follows, while the wave-speeds appearing in Eqs. (19) and (28), take the following form, The function E(R, χ) is written with respect to the e-foldings number as follows, Using Eqs. (29), (30), (47) and (52), we obtain the slow-roll parameters for the de Sitter evolution case with an exponential coupling function, which is, By using the above results, we can proceed in calculating the spectral indices, from Eqs. (23), and the tensor-to-scalar ratio, from Eq. (24), In order to examine the viability of the model, we need to calculate the numerical values for the spectral index n S and the tensor-to-scalar ratio r, for various values of the parameters H 0 , γ, b and µ at the end of inflation (for N ∈ [50, 60]) and compare these values to the observational results of the Planck collaboration [17] and the BICEP2/Keck-Array [18]. However in this case, no simultaneous compatibility with the observations can be obtained, and more specifically, the values of n S and r do not depend on the choice of γ, so we set it equal to one for simplicity. They also do not depend on the number of e-foldings, so N = 50 and N = 60 are used in the same manner. They depend on H 0 , b and µ, though, thus assuming that H 0 ∼ 10 27 sec −1 and setting µ = 10 12 sec −1 , we get b = 35.6 so that n S = 0.9644 (Planck's previous result) however, the resulting value of the tensor-to-scalar ratio is excluded. This can also be seen in Fig. 2.

V. A FLAT QUASI-DE SITTER VACUUM AS BACKGROUND
Now we assume that the Universe's evolution is described by the quasi-de Sitter Hubble rate, Integrating Eq. (56) with respect to the cosmic time, we obtain and solving with respect to time, we may write the latter with respect to the e-foldings number as follows, As a result, the Hubble rate with respect to the e-foldings number becomes, while the Ricci scalar and the Gauss-Bonnet scalar are equal to, Finally, we also express the scalar field of Eq. (10) with respect to the e-foldings number as follows, As in section IV, the Eqs. (30), (31) and (33) and a coupling function allow us to reveal the phenomenological implications of the model by calculating the observational indices of inflation.
A. An exponential coupling function, h(χ) = γe bχ At first, we shall assume that the function h(χ) has the functional form given in Eq. (45), which in the case at hand is written in terms of the e-foldings number as follows, By using Eq. (16), we may derive the potential as a function of the e-foldings number, as well as the Lagrange multiplier, The functions Q i with respect to the e-foldings number are derived from the Eqs. (27), while the wave-speeds are Interestingly, both the Q i functions and the wave-speeds have a trivial form in the case of the quasi-de Sitter expansion. This triviality is independent of the coupling function, as we see later on, and should be attributed to this specific FRW background.
The function E(R, χ) with respect to the e-foldings number takes the form, Using Eqs. (29), (30), (47) and (52), we obtain the slow-roll parameters of the flat quasi-de Sitter case with an exponential coupling function. Interestingly, the five of them take the following trivial form, that seems independent of the coupling function, while the fourth has a long and complex form depending on the coupling function, where ε exp is some notation for the complicated functional form of the slow-roll index ǫ 4 . Similarly, the spectral indices and the tensor-to-scalar ratio are also long and complex functions of the e-foldings number, the mass µ and the model parameters, H 0 and H 1 due to the expansion rate and γ and b due to the coupling function, thus we do not present them in close form. What is interesting to note is that the spectral indices and the tensor-to-scalar ratio yield the same values independently of which case of Eq. (57) we will use. Again, we perform comparisons using the observable values for n S and r obtained by the Planck with their latest data [17], along with [18]. As we stated before, the spectral index of the scalar modes must be within the interval [0.9607, 0.9691] and mean n S = 0.9649; the tensor-to-scalar mode, on the other hand, is restricted below 0.1 by [18], while [17] restricts further as r < 0.064. In our case, the parameters γ and b, as well as the mass µ of the scalar field seem not to affect the numerical values of the spectral index or the tensor-to-scalar ratio. As a result, we consider them equal to unity (γ = b = 1 and µ = 1 sec −1 ), so that the analysis is simplified and focused on the rest of the parameters. The e-foldings number is chosen N = 50 and N = 60, so as to indicate the end of inflation, but this also does not alter the results. As for the expansion rate, given that H 0 ≥ 10 14 sec −1 for H 1 ≈ 10 26 sec −2 (or that H 0 ≥ 5 × 10 14 sec −1 for H 1 ≈ 10 27 sec −2 ), the spectral index approaches unity, restricting our choices. We consider H 0 to be in the interval [10 12 , 10 15 ] sec −1 and H 1 in the respective interval [10 26 , 10 29 ] sec −2 , where the spectral index of our model equals to the observable value, as we can see in Fig. 3. For the majority of these cases, the tensor-to-scalar ratio is close to zero, as we can see in Fig. 4. As an example, choosing N = 50 (or N = 60) and H 1 = 10 27 sec −2 , then for H 0 = 4.91375 × 10 14 sec −1 , we have n S = 0.9644 and r = 0.0282787, which comply with the latest data of the Planck collaboration. What we need to notice is that these two parameters (H 0 and H 1 ) need careful fine-tuning and cannot differ significantly for the set of values we gave, otherwise the model collapses before the data. Now let us assume that the function h(χ) takes the form given in Eq. (34), which in terms of the e-foldings number is expressed as follows, From here, using Eq. (16), we may derive the potential as a function of the e-foldings number, as well as the Lagrange multiplier, The Q i functions with respect to the e-foldings number have the same trivial form given in Eqs. (64) and (65). The function E(R, χ) with respect to the e-foldings number takes the form, Using Eqs. (29), (30), (68) and (71), we obtain the slow-roll parameters of the flat quasi-de Sitter case with an exponential coupling function. Except from the fourth one, which has a long and complex expression, (N, H 0 , H 1 , a, b, µ) , the rest are given in Eqs. (67). The spectral indices and the tensor-to-scalar ratio have the same form as in the case of the exponential coupling function, presented above. Again, we perform comparisons using the observable values for n S and r obtained by the Planck with their latest data [17], along with [18]. We assume that the parameters γ and b, as well as the mass µ are equal to unity (γ = b = 1 and µ = 1 sec −1 ), so that the analysis is simplified and focused on the rest of the parameters. The e-foldings number is chosen N = 50 and N = 60, and as for the expansion rate, given that H 0 ≥ 10 14 sec −1 for H 1 ≈ 10 26 sec −2 (or that H 0 ≥ 5 × 10 14 sec −1 for H 1 ≈ 10 27 sec −2 ), the spectral index approaches unity, restricting our choices. We consider H 0 to be in the interval [10 12 , 10 15 ]sec −1 and H 1 in the respective interval [10 26 , 10 29 ]sec −2 , where the spectral index value of our model becomes equal to the observable value, as we can see in Fig. 5. For the majority of these cases, the tensor-to-scalar ratio is close to zero, as we can see we have n S = 0.9644 and r = 0.0400002 (or r = 0.0333335), that match the latest data of the Planck collaboration.
Again, these two parameters (H 0 and H 1 ) need careful fine-tuning and cannot differ significantly from the above values.

VI. THE CASE OF AN EXPONENTIAL HUBBLE EVOLUTION
Finally, let us assume that the evolution of our Universe is described by the following Hubble rate, where H 0 and Ω are model parameters with both having mass dimension [+1]. The Hubble rate of Eq. (73) becomes approximately a quasi de-Sitter like evolution at early times, when t → 0, that is, and also the exit from the inflationary epoch occurs at a finite time t f , which is, Moreover such exponential type Hubble parameter has been used in previous works, in the context of f (R) gravity [38,39] as well as in different theoretical frameworks [40,41]. Motivated by such properties of H = H 0 e −Ωt , here we use it in the context of ghost free f (G) gravity to describe the inflationary phase of our Universe we will test the viability of the model by confronting it with the Planck 2018 constraints. Also, we can express the cosmic time as a function of the e-foldings number N , by using the definition of the latter, where t h is the horizon crossing time instance. Inverting Eq. (76), we get t h in terms of N as follows, This expression of t h is important, since the inflationary parameters will be calculated at the horizon crossing time instance. In the following we will calculate the slow-roll indices and the observational indices of inflation by specifying the function h(χ).

A. Exponential coupling : h(χ) = e −αχ
Let us assume that h(χ) = e −αχ where α is a model parameter having mass dimension [-1]. For this exponential function h(χ), and also for the Hubble rate chosen as in Eq. (73), the scalar potential is equal to, while the Lagrange multiplier is equal to, Accordingly, the function E defined in Eq. (21) evaluated at the horizon crossing time instance, so by expressing it in terms of the e-foldings number, this reads, We also need to evaluate the expression ofĖ (= dE/dt) as it will be needed for the calculation of the slow-roll indices. In terms of the e-foldings number, this reads, with T = Ω(1+N )

H0
. Furthermore, by using Eq. (18), we explicitly determine the functions Q i in terms of e-foldings number, Having the above expressions at hand, we can easily calculate the spectral index n s and tensor-to-scalar ratio, which are, and where A 1 and B 1 are defined as follows, It may be noticed that n s and r depend on the parameters Ω/H 0 , αµ 2 /Ω, κH 0 and N . We can now directly confront the spectral index and the tensor-to-scalar ratio with the Planck 2018 constraints and the BICEP-2 Keck-Array data, which recall that constraint the observational indices as: n s = 0.9649 ± 0.0042 and r < 0.064, as shown earlier.
For the model at hand, n s and r lie within the Planck constraints for the following ranges of parameter values: 0 Ω/H 0 ≤ 0.035 , 0 αµ 2 /Ω ≤ 1.5 with κH 0 ∼ 0.01 and N = 60 and this behavior is depicted in Fig. 7. At this stage it deserves mentioning that an exponential coupling function in a scalar GB theory (without scalar field potential) admits, at early times, slowly expanding solutions of the form a(t) = (At + B) 1/5 (see [42]) and thus exhibits an epoch of deceleration. However here, we show that in the presence of ghost free f (G) gravity, the exponential coupling function may be considered as a "good inflationary" model, which allows an early acceleration and also it is compatible with observations. Before closing, we can also notice that if some sort of slow-roll conditions are employed in the model, viability with the observational data can also be achieved. The slow-roll conditions in the ghost free Gauss-Bonnet scenario are the following, The first condition carries the information about the slow-evolution of the Hubble rate, while the last two demand a slowly evolving of the function h(χ). These conditions, and especially the last two can significantly constrain the and the Lagrange multiplier is, Furthermore, the function E(R) and consequently its derivative, evaluated initially at the horizon crossing time instance, and expressed eventually in terms of the e-foldings number, are equal to, with S = ln H0 Ω (1+N ) . Accordingly, the functions Q i , in terms of e-folding number, are equal to, Hence, the spectral index becomes in this case, and the tensor-to-scalar ratio is equal to, respectively, with A 2 and B 2 being defined as follows, and B 2 (Ω/H 0 ,µ 2 /(ΩM ), n, κH 0 , N ) Now we shall confront the resulting theory with the observational constraints, by assuming two different values for the parameter n, namely, n = 2 and n = 3. For n = 2, the tensor-to-scalar ratio acquires a minimum value r min = − 4 (1+N ) which is equal to r min = 0.065 for N = 60 ( and to r min = 0.078 for N = 50 ). The behavior of the tensor-to-scalar ratio as a function of the free parameters, is given in Fig. 8. As it can be seen in Fig. 8 the  Fig. 9 we can see that the spectral index and the tensor-to-scalar ratio can be simultaneously compatible with the observational data, for a wide range of values of the free parameters. Before closing this subsection, we need to comment that it was shown in [42] that a quadratic coupling function in a scalar GB theory (without scalar field potential) gives either a pure de-Sitter evolution of our Universe or a de-Sitter solution at early times connected by a Milne phase at late times, while the cubic and higher order coupling functions describe contracting cosmological solutions with a final singularity at asymptotically infinite time. Thus none of the power law coupling function corresponding to n ∈ [2, 3] admits a successful inflationary model in scalar GB theory in the absence of scalar potential. However in the context of ghost free f (G) gravity, we demonstrated that h(χ) ∼ χ n with n ∈ [2,3] can realize an accelerating Universe at early times, although only the cubic coupling function h(χ) ∼ χ 3 produces a viable inflationary phenomenology, in contrast to the models studied in [42].

C. A Different Reconstruction Approach
In this subsection we shall consider an alternative approach in comparison to the previous cases, by providing the scalar potential and the Hubble rate, and we seek for the function h(χ) that may realize the cosmology with Hubble rate (73). We shall consider two types of potentials, namely exponential and power law potentials and we shall confront the resulting theories with the observational data.
and the Lagrange multiplier is, With the above expressions of h(χ) and λ, we get the function E(R) as well as its derivative, which are, respectively, with T = Ω(N +1)

H0
, defined earlier. In addition, the functions Q i as functions of the e-foldings number become in this case, Having the above expressions in hand, we determine the explicit expressions of the spectral index and of the tensorto-scalar ratio, which are, respectively, where we took V 0 = Ω 4 . Moreover C 1 , D 1 appearing in Eq. (97)) are defined as follows, From Eqs. (97) and (98), it easily follows that the spectral index of scalar perturbation and the tensor-to-scalar ratio depend on the dimensionless parameters : Ω/H 0 , βµ 2 /Ω, κH 0 and N . These theoretical expressions of n s and r should be confronted with the latest Planck constraints in order to check the viability of the model. As a consequence, it is found that the compatibility with the observational data occurs for a narrow range of values of the free parameters, and particularly for 0.001 ≤ Ω/H 0 ≤ 0.02, 82 ≤ βµ 2 /Ω ≤ 83, κH 0 ∼ 0.01 and N = 60. This can also be seen in Fig. 10 where we present the parametric plot of n s and r. With regard to the exponential potential, the classical single scalar theory has no inherent mechanism to trigger the graceful exit from inflation, since the slow-roll indices are constant and field-independent. However the ghost free f (G) theory has the slow-roll index ǫ 4 which is field dependent, and thus the slow-roll phase ends when this index becomes of order O(1). Moreover we have already shown that the model with V = V 0 e −βχ in f (G) gravity, is also in agreement with Planck observational constraints. Hence the ghost free f (G) gravity can make the exponential scalar potential a phenomenologically appealing model for inflation, in contrast to the single scalar canonical exponential theory.

Power law scalar potential
As a final consideration, we shall assume that the scalar field potential has the form, where n is a positive integer. For such power law potential, the function h(χ) and Lagrange multiplier are equal to, and respectively. Accordingly the function E(R) expressed in terms of the e-foldings number is equal to, and also its derivative is, For the Hubble rate given in Eq. (73) and with the expression of h(χ) we found above, we can easily find the Q i functions expressed in terms of the e-foldings number, Let us use the above results in order to investigate the viability of a power-law class of potentials. According to the latest Planck data, the cubic and quartic potentials are not compatible with the Planck data, so let us investigate whether compatibility with the observations is obtained if the ghost free f (G) theory is used. Let us first assume that n = 3 so we consider the cubic potential first. Using V (χ) = V 0 χ 3 along with the explicit expressions of Q i functions (see the equations in 103, we determine the spectral index and tensor to scalar ratio in terms of the model parameters as follows, and r = 3N (N + 1) where we assumed that V 0 = H 0 (for the cubic potential, V 0 has mass dimension [+1]). Moreover C 2 and D 2 have the following form, parameters, and in particular for 0.001 ≤ Ω H0 ≤ 0.003, 50 µ 2 Ω 2 (κH 0 ) 2/3 52 and N = 60, as shown in Fig. 11. We should note that the single canonical scalar field model with cubic potential without the Gauss-Bonnet coupling yields n s ≃ 0.9089 and r ≃ 0.01, so the spectral index is not compatible with the Planck data. Hence, the presence of the ghost free f (G) gravity can make the cubic potential scalar field class of models to be compatible with the observations. This kind of result is also shown in a different context [36]. Let us now consider the n = 4 case, in which case the potential is V = V 0 χ 4 . In this case, the spectral index of the primordial scalar curvature perturbations and the tensor-to-scalar ratio are equal to, inflationary observational indices lie within 0.960 ≤ n s ≤ 0.970 and 0.049 ≤ r ≤ 0.065 respectively. Thus the model becomes viable (with respect to the Planck 2018 constraints) for such narrow parameter space. However as may be noticed that µ 2 Ω 2 (κH 0 ) 1/2 must be fine-tuned within the values 109 and 110.9 to keep the model compatible with Planck constraints. The simultaneous compatibility of n s and r is illustrated in Fig. 12. However the single canonical scalar field theory with a quartic potential yields n s ≃ 0.8677 and r ≃ 0.066 for 60 e-foldings, so the spectral index of the corresponding canonical scalar field theory is excluded by the latest observational data. Hence, the presence of ghost free f (G) theory modifies the quartic scalar field theory and enhances the phenomenological viability of the model.

D. Stability of First Order Perturbations for the Exponential Cosmological Evolution
In this subsection we shall study stability of the first order perturbation cosmological perturbations, following the work of [43][44][45][46][47][48] where scalar, vector and tensor perturbations are calculated in the context of Gauss-Bonnet theory. Scalar, vector and tensor perturbations are decoupled, as in general relativity, so that we can focus our attention to tensor and scalar perturbations separately, as discussed in what follows. Let us consider first tensor perturbations, which the flat FRW perturbed line element has the form, where f ij (t, x) is the tensorial perturbation satisfying f i i = f ij ,j = 0. Plugging back the above metric into the original action and expanding, keeping terms up to O(f 2 ) (to obtain the first order equations), we get the following perturbed action [43,47,48], where we use the background equations of motion. With the Fourier decomposition as f ij (t, x) = dkf ij (t)e i k. x , the above perturbed action takes the following form, Thereby, the tensor perturbation is ghost free and stable if the following two conditions hold true, 1 + 8κ 2ḣ H >0 , and are satisfied simultaneously. If we assume that the slow-roll conditions of Eq. (85) hold true, the coupling function h(χ) rolls slowly if it obeys ḣ H ≪ 1/κ 2 and ḧ ≪ 1/κ 2 . We have shown in previous sections that a phenomenologically viable cosmological evolution also satisfies these constraints if the free parameters are chosen appropriately, so in view of Eq. (111) we may conclude that the tensor perturbations are ghost free and stable, at least at first order. So the theory is compatible with the observational data and stable up to first order cosmological tensor perturbations. Now let us turn our focus to scalar perturbations on the FRW background spacetime, in which case the line element is, with Ψ(t, x) being the scalar perturbation. Following [43], the perturbed action up to order O(Ψ 2 ) is equal to, where Z 1 and Z 2 are defined as follows, Clearly the scalar perturbation is ghost free and stable if Z 1 and Z 2 are both positive. With the slow-roll criteria taken into account, and the corresponding field equations, the positivity of Z 1 , Z 2 is guaranteed if the following two conditions hold true,Ḣ Now let us proceed to explore whether, for our considered choice of coupling or potential function, the above two conditions are in agreement with the Planck 2018 constraints. The first condition is satisfied for Ω > 0 ( asḢ = −ΩH 0 e −Ωt ) and simply gives the information that the Hubble parameter must decrease with cosmic time at the early universe, which is also expected in an inflationary scenario. On the other hand, it is shown that all the previous four cases (see Sections from [VI A] to [VI C 2]) need Ω > 0 in order to be compatible with Planck constraints and thus one of the stability condition of scalar perturbation is ensured. Now let us investigate the second condition case by case: for the exponential coupling i.e. h(χ) = e −αχ ,ḧ −ḣH becomes positive for α > 0 which is also needed to make the model observationally viable (as shown in Section [VI A]). To investigate what happens in the power-law case of h(χ), we provide the plot ofḣ Ḧ h as a function of the e-foldings number in Fig. 13. As it can be seen in Fig. 13, the ratioḣ Ḧ h remains less than unity for all the parameter values that render the theory compatible with the latest Planck data. Thus this ensures numerically the stability of the scalar perturbations. E. Reheating mechanism for the exponential cosmological evolution Before moving to the conclusion section, here we discuss the phenomenological implications of the theory we studied in the reheating era, and the possible effects of Gauss-Bonnet coupling on it. Needless to say that reheating describes the production of Standard Model matter after the period of accelerated expansion. For this purpose, we assume that the inflaton field (i.e the field χ) is coupled to another scalar field ζ, given by the interaction Lagrangian, where g is a dimensionless coupling constant and λ is a mass scale. The scalar field ζ quantifies Standard Model particles in our case study. With this interaction Lagrangian, the decay rate of the inflaton into ζ particles becomes, where m denotes the mass of the inflaton field and can be obtained from the effective potential V ef f (χ) = V (χ) − 24H 2 (H 2 +Ḣ)h(χ) through which the Gauss-Bonnet coupling function (h(χ)) affects indirectly the reheating mechanism. Moreover, the presence of Gauss-Bonnet term also affects the self-potential functionṼ (χ) as may be noticed in Eq. (16) ( see the terms dependent on h(χ) in the right hand side of Eq. (16)). Generally during the reheating epoch, the inflaton losses energy due to the expansion of the Universe, and due to transfer of energy to the ζ particles, controlled by the Hubble parameter and the decay rate respectively. As a result the production of ζ particles becomes effective when the Hubble parameter becomes less or comparable to Γ, otherwise the energy loss into particles is negligible compared to the energy loss due to the expansion of space as occurred during the early phases of the inflation. Therefore, the time scale t h (let us call it the reheating time) after when the production of ζ becomes effective is given by, where we used the form of the functionṼ (χ) as obtained in Eq. (78). Consequently the stable point (< χ > (ec) , where the notation "ec" stands for "exponential coupling") of V ef f can be determined by the following algebraic equation, The presence of h 0 in the above expression entails that the Gauss-Bonnet coupling indeed affects the stability of the inflaton. In order to understand the effect of the GB coupling more clearly, we write < χ > (ec) =< χ > 0 + < δχ > (ec) , where < χ > 0 is the stable point of χ in absence of GB term (h 0 = 0) i.e., Thus < δχ > (ec) is the deviation of stable point from χ > 0 solely due to the presence of the Gauss-Bonnet term. Expanding Eq. (120) in terms of < χ > (ec) =< χ > 0 + < δχ > (ec) , we get the following expression for < δχ > (ec) , where we kept terms up to first order in < δχ > (ec) and we also assumed Ω αµ 2 > 1, which is also consistent with the Planck observations, as mentioned earlier in Section[VI A]. Clearly < δχ > (ec) becomes zero as h 0 → 0, as expected. Eqns. (121) and (122) immediately lead to the stable point of V ef f in presence of Gauss-Bonnet coupling, which is, Using the above expression for < χ > (ec) , we determine the mass squared of the inflaton (m 2 (ec) ) for the case of exponential coupling function, which is, Thus in the absence of the Gauss-Bonnet term (i.e for h 0 = 0), m 2 (ec) becomes m 2 (ec) = 2Ω 4 µ 4 κ 2 which is also consistent with Eq. (121). However the presence of exponential coupling affects the inflaton mass by the factor proportional to h 0 , as is evident from Eq. (124), in particular the mass increases due to the presence of the Gauss-Bonnet term, compared to the case where h 0 = 0. Having the explicit expression of m 2 (ec) (see Eq. (124)) at hand, now we can determine the reheating time by using Eq. (118), which is, Thus we can argue that the presence of the exponential GB coupling function, enhances the mass of the inflaton which in turn makes the reheating time larger compared to the situation where the Gauss-Bonnet term is absent. For the cubic coupling (h(χ) = h 0 χ/M 3 ), the effective potential of the inflaton is equal to, Following the same procedure as above, we determine the stable point of the effective potential and the mass of the inflaton field, in the case of cubic coupling, which are, < χ > (cc) = µ 2 Ω ln 3H 0 /Ω − respectively, where x 0 = ln 3H 0 /Ω and the notation 'cc' stands for "cubic coupling". Thereby, the presence of cubic GB coupling function makes the inflaton mass larger relative to the situation where the GB term is absent.
As a consequence, the reheating time t (cc) h = 1 Ω ln 8πm (cc) H0 g 2 λ 2 also increases due to the effect of the cubic coupling function, similar to the case of the exponential coupling we discussed earlier.
Before closing, let us comment on an interesting issue, related to previous works in the field. In Ref. [49], the authors calculated the observational indices of inflation for a generalized Galileon theories, however these theories are quantitatively different to a great extent from the theory we developed in this paper. Particularly, the theory at hand with action (4) can be treated at a quantitative level as an generalized Einstein-Gauss-Bonnet theory of gravity, which is entirely different from the Galileon models studied in Ref. [49]. At a quantitative level, the theories developed in Ref. [49], allow the derivation of general forms of the observational indices, however in our case, and in Einstein Gauss-Bonnet models, it is hard to derive general relations for the observational indices. This is because the latter depend strongly on the choice of the Gauss-Bonnet scalar coupling function h(φ). Thus the results are strongly model dependent, as we evinced in the previous sections, for example in Section VI A with h(χ) = e −αχ or in Section VI B with h(χ) = χ/M n . As we have shown, the quadratic coupling is not viable although the exponential one in section VI A is viable. We have further extended our discussion to investigate the possible effects of GB coupling function on reheating mechanism, unlike to [49].

VII. CONCLUSIONS
In this work we studied the inflationary phenomenology of a recently developed ghost-free f (G) model of gravity. Particularly, the form of the model mimics the scalar Einstein-Gauss-Bonnet theory, so we employed the formalism of cosmological perturbations for the latter theory, in order to calculate the slow-roll indices and the corresponding observational indices for the theory at hand. The model has rich phenomenology due to the presence of a freely chosen function h(χ), in which case by choosing this function and the Hubble rate, the observational indices can be calculated easily. We examined three types of inflationary cosmic evolution and functional forms of the function h(χ), and as we demonstrated it is possible to have a viable inflationary era, compatible with the latest observational data. Particularly we used de Sitter, quasi-de Sitter and exponential cosmological evolutions, and also exponential and power-law functional forms for the function h(χ). The simple de Sitter evolution leads in some cases to problematic phenomenology, however no realistic cosmology gives the exact de Sitter evolution, so we investigated the quasi-de Sitter case, in which case the viability of the theory with the observational data comes more easily. The same applies for the exponential cosmological evolution. For the exponential Hubble rate case, we also tested the stability of the first order scalar and tensor cosmological perturbations, and as we demonstrated these are stable for the same range of values of the free parameters, for which the phenomenological viability of the model is ensured. Finally we explore the reheating mechanism and the possible effects of Gauss-Bonnet term on it for the case of exponential Hubble rate. As a result we found that the presence of GB coupling, in particular the exponential and cubic coupling function, enhance the mass of the inflaton which in turn makes the reheating time larger compared to the situation where the Gauss-Bonnet term is absent. In this work we mainly focused on realizing inflationary evolutions, however it is also possible to realize non-singular cosmological evolutions, such as cosmological bounces, however we defer this task to future work.