Cosmological perturbations in a class of fully covariant modified theories: Application to models with the same background as standard LQC

Bouncing cosmologies are obtained by adding to the Einstein-Hilbert action a term of the form $\sqrt{-g}f(\chi)$, with $\chi$ a scalar depending on the Hubble parameter only, not on its derivatives, and which is here shown to arise from the divergence of the unitary time-like eigenvector of the stress tensor. At background level, the dynamical equations for a given $f$-theory are calculated, showing that the simplest bouncing cosmology resulting leads to exactly the same equations as those for holonomy corrected Loop Quantum Cosmology (LQC). When dealing with perturbations, the equation for tensor ones is the same as in General Relativity (GR); for scalar perturbations, when one uses the $f$-theory which leads to the same background as the standard version of holonomy corrected LQC, one obtains similar equations (although a bit more elaborated) as those coming from LQC in the so-called deformed algebra approach.


INTRODUCTION
One of the most simple bouncing backgrounds (see [1][2][3] for recent reviews about bounces) is obtained from holonomy corrected Loop Quantum Cosmology (LQC), where the corresponding Friedmann equation depicts an ellipse on the plane (H, ρ) [4][5][6][7][8][9], being H the Hubble parameter and ρ the energy density. As shown in several papers, this simple background can be mimicked by modifying the Einstein-Hilbert action through the introduction of a term of the form √ −gf (χ) where g is the determinant of the metric, f is a well-known function [10][11][12][13] and χ is a scalar, only depending on H but not on its derivatives. The problem with this method is to actually find a scalar that for synchronous observers in the Friedmann-Lemaître-Robertsont-Walker (FLRW) spacetime will only depend on the Hubble parameter.
From our viewpoint the simplest scalar is the extrinsic curvature [14], which appears in a natural way when using the ADM formalism [15]. Disappointingly, this is not a covariant theory because the extrinsic curvature is not a true scalar in the sense that it depends on the slicing chosen and, thus, its use is only justified if there exists a preferred foliation of the spacetime. Following the spirit of Weyl's principle (see [16] for a historical review), one could choose a preferred slicing as follows. The time-like eigenvector of the stress tensor, which always exists for realistic matter (see pages 89-90 of [17]), generates a preferred non-crossing family of world-lines and one can construct, at any given time t, a family of hypersurfaces ortogonal to these world-lines, obtaining in this way the so-called co-moving slicing.
Another simple scalar could be obtained by working in the Weitzenböck spacetime (the usual Levi-Civita connection is replaced by the Weitzenböck one) [18], where torsion does not vanish. In this spacetime, the scalar torsion for synchronous observers in the flat FLRW geometry is equal to minus three times the Hubble parameter, thus satisfying the required property. Unfortunately, as has been shown in [19], this theory is not locally Lorentz invariant.
For this reason, people have kept looking for a really covariant invariant, and have dealt with the Carminati-McLenaghan invariants [20,21], with a general function of the Ricci and Gauss-Bonet scalars or with second derivatives of the Riemann tensor [22]. The problem with these purely gravitational invariants is that they are quite involved and lead to very complicated equations for cosmological perturbations. Moreover, in our approach -not using the principle of the limiting curvature hypothesis considered in [22]-they can easily lead to Ostrogradski or gradient instabilities, as well as to the appearance of ghost fields. All these problems led specialists to explore other ways, such as modified mimetic gravity [13,23,24], so as to find out such scalar.
Following these arguments, guided by Weyl's principle and taking into account that for co-moving observers in flat FLRW spacetime the divergence of a unitary time-like vector is equal to −3H, in the present paper we propose as our fully covariant scalar the divergence of the unitary time-like eigenvector of the stress tensor, which, in the case of a universe filled up with a scalar field φ minimally coupled to gravity, is equal to u µ = φ,µ √ φ,ν φ ,ν (throughout the work we will use the notation: φ ,µ ≡ ∂ µ φ = ∇ µ φ).
With this covariant scalar, we will show how to construct f -theories leading to bouncing backgrounds and calculate the perturbed equations for scalar and tensor perturbations for a given f -theory. When dealing with perturbations and working in the longitudinal gauge -where the Newtonian potential, namely Φ, and the variation of the scalar field, namely δφ, are the dynamical variables-the corresponding dynamical system turns out to be a coupled one. This is an essential difference with respect to theories such as General Relativity (GR) or LQC in the deformed algebra approach, where the dynamical equation for the potential Φ decouples (δφ does not appear, see for instance [30,31]).
Once these equations are obtained, we study some characteristics of the matter-ekpyrotic bounce scenario [34,35] as the calculation of some spectral quantities and the reheating temperature via the gravitational particle production of massless particles, in the contracting regime, during the phase transition from matter domination to the ekpyrotic regime.
The manuscript is organized as follows: In Section II we present our class of modified gravitational theories and obtain the corresponding dynamical equations. We study them at the background level and find which is the model that leads to the simple bounce predicted by holonomy corrected LQC. A Hamiltonian analysis of our theory is performed in Section III, which leads to the conclusion that the present theory, as in the case of mimetic gravity, has one more degree of freedom than GR. In section IV, we study scalar and tensor perturbations. For scalar perturbations, working in the longitudinal gauge, we obtain the equations for the Newtonian potential and for the perturbed part of the scalar field, showing that they are coupled. Moreover, we derive the corresponding Mukhanov-Sasaki equations for our theory. Then, dealing with tensor perturbations, we show that, since the modification of our theory is performed on the matter sector, the equations must be the same as for GR. Section V is devoted to the comparison, at the perturbative level, of our model which leads to the same background as holonomy corrected LQC, with other theories which also lead to the same background, as: LQC in the deformed algebra approach [31][32][33], teleparallel LQC [36,37], extrinsic curvature LQC [14] and mimetic LQC [12,13,24]. In Section VI, we study the matter-ekpyrotic scenario applied to the model that leads the same background as LQC. We calculate the spectral index and its running and show that they match the most recent observational data. Moreover, we study the reheating process via gravitational massless particle production. Finally, the last Section is devoted to conclusions.
The units used throughout the paper are = c = 1 = M pl = 1, where M pl is the reduced Planck mass with the convention that a temporal vector v µ satisfies v µ v µ < 0 and with the notation: 1. ϕ ,µ ≡ ∂ µ ϕ = ∇ µ ϕ for a given scalar ϕ.
2.ḡ is the unperturbed part of g.
3. g χ means derivative of g with respect to χ and g φ derivative of g with respect to φ, for a given function g.

A CLASS OF MODIFIED GRAVITATIONAL THEORIES
All models to be considered here come from a simple action, which consists in adding a term of the form f (χ) to the Einstein-Hilbert action. Here f is a given function which, in order to recover GR, vanishes at low energy densities, while χ is a fully covariant scalar built from the scalar field that fills the whole universe and whose value in any co-moving slicing of the flat FLRW spacetime is proportional to the Hubble parameter. More precisely, we consider R being the scalar curvature. We have assumed that the matter sector of the universe is described by a scalar field φ with a potential, V (φ) which is minimally coupled to gravity and whose Lagrangian is given by and χ ≡ −∇ µ u µ , being the vector u µ the time-like unitary eigenvector of the stress tensor that is, The dynamical equations are obtained by performing the variation of the action with respect to g µν , leading to where G µν = R µν − 1 2 g µν R is the well-known Einstein tensor, while the tensorT µ,ν , coming from the term f (χ), is given byT On the other hand, the variation of the action with respect to the scalar field φ leads to the following conservation equation, which differs from the usual one − φ+V φ = 0. However, for the FLRW geometry it leads to the standard conservation equation,φ

The background
Considering backgrounds (i.e., solutions of (8)) satisfyingφ(t) > 0 all the time, the simplest way to obtain the dynamical equations is as follows: We consider the flat metric ds 2 = −N 2 (t)dt 2 + a 2 (t)δ ij dx i dx j , which leads to χ = 3H N , where N (t) is the lapse function. Hence, after integration by parts, the action becomes where V is the volume of the fixed elementary spatial cell, where all spatial integrations are performed, and where the matter Lagrangian simplifies to Performing the variation with respect to N and taking at the end N = 1, one obtains the modified Friedmann equation for synchronous observers which depicts a curve on the plane (H, ρ). Finally, taking its temporal derivative and using the conservation equation (8) in the formρ = −3H(P + ρ), one finds the Raychaudhuri equation Remark 2.1 A different way to obtain such equations is to directly consider the metric for synchronous observers ds 2 = −dt 2 + a 2 δ ij dx i dx j and use the equations 0 − 0 and i − i, which also leads to the modified Friedmann and Raychaudhuri equations.
Therefore, given a curve ρ =ḡ(3H) =ḡ(χ) on the plane (H, ρ), so as to obtain the correspondingf (χ) theory one has to solve the first-order differential equation whose solution isf Note that, in order to obtain a bouncing background, a necessary condition is to choose a curve on the plane (H, ρ) cutting two or more times the axis H = 0. One of those points has to be (0, 0) because at low energy densities a viable bouncing background has to approach to GR, i.e., the chosen curve has to approach to the parabola H 2 = ρ 3 and the other points cutting the axis H = 0 have the form (0, ρ i ) with ρ i > 0.
The simplest example is the ellipse coming from the holonomy corrected Friedmann equation in standard LQC [38][39][40] (Note that here we are not dealing with the recent model of LQC proposed by Dapor and Liegener in [41]) where the so-called critical density has the value ρ c ∼ 0.4ρ pl = 0.4 × 64π 2 ∼ = 252 (see for instance [7]). For this background, using equation (14), the corresponding f -theory is given bȳ where s ≡ 2 √ 3ρcχ and the functions √ 1 − s 2 and arcsin s are bi-valued [21]. For the theory to be well defined, since the ellipse has two branches -the upper part corresponding to ρ = ρc Note also that, a simple calculation shows that, for this particular f -theory, the modified Friedmann and Raychaudhuri equations read and coincide with the ones obtained in holonomy corrected LQC [42,43].
A final remark is in order. Since in our case the matter sector is depicted by a non-phantom scalar fieldφ (see equation (10)), we will have infinitely many backgrounds because the conservation equationφ + 3Hφ + V φ (φ) = 0 is a second order differential equation. Effectively, for any initial conditionφ(0) = α 0 andφ(0) = α 1 , one has a different background. Therefore, if we choose a curve ρ = g(χ) in the plane (H, ρ) containing a point of the form (0,ρ) with ρ > 0 and initial conditions satisfying its corresponding solution will lead to a bouncing background provided that V φ (α 0 ) = 0.
As an example one can consider the matter bounce scenario in LQC, which is given by the potential [50]. In this case, the conservation equation has the following analytic solution leading to the following background which is the same nonsigular bouncing background obtained solving the equation (17) when one considers a pressureless fluid and whose dynamics is depicted in Figure 1. All the other solutions, which lead to different backgrounds, are obtained choosing different initial conditions (see Figure 3 of [14], where it is showed that nearly all solutions lead to nonsingular bouncing backgrounds).
Moreover, since we are dealing with theories beyond GR, in order to have a bouncing background it is not needed to violate the null energy condition ρ + P ≥ 0 near the bounce using quintom or Lee-Wick matter (see [44] and references therein) because the bounce occurs when the value of the energy density is strictly positive. For example, looking at (19) one can see that at the bouncing time t = 0 the energy density is given by ρ c and the pressure, which could be obtained from the Raychaudhuri equation (17), is zero. So, the null energy condition is fulfilled at the bounce.

HAMILTONIAN ANALYSIS
In this section we perform a Hamiltonian analysis in order to find the degrees of fredoom of our model. To this end, we will use the ADM formalism [15], where the line element acquires the form the matter Lagrangian becomes and the action is given by where R is the intrinsic curvature, i.e. the scalar curvature of Σ t , and the extrinsic curvature tensor, with D the induced Levi-Civita connection in the slicing Σ t .
Introducing a Lagrangian multiplier field, namely β, the action becomes with the canonical momenta being and the constraints Note that, from the equation , we obtainφ as a function of P β . In fact, taking the square of this expression one gets −φ ,ν φ ,ν = γγ ij φ,iφ,j P 2 β −γ and, thus, And, from , we readily obtainβ as a function of P β and P φ . What is important is that, after the Legendre transformation, one can get the Lagrangian as follows, where H matt is the matter part of the Hamiltonian; α N , α N i and α χ are Lagrange multipliers, and the functions H and H i lead to the hamiltonian and diffeomorphism constraints, which come from imposing stability under time evolution of the constraints P N ≈ 0 and Now we examine the stability of the constraint P χ . Looking for the only term in the Hamiltonian where χ appears and for the term where it appears in the action (24) is ensured by fixing the Lagrange multiplier α χ as α χ = 0.
Summing up, we have obtained the constraints H ∼ = 0, H i ∼ = 0, P χ ∼ = 0 and C χ ∼ = 0, and the canonical pairs (q ij , P ij ), (φ, P φ ), (χ, P χ ) and (β, P β ). Then, from the constraints P χ ∼ = 0 and C χ ∼ = 0 one may remove two variables, for example P χ ∼ = 0 and β ≈ −f (χ), thus obtaining, as in mimetic gravity [12], one more degree of freedom than in the case of GR. However, as we will show in next Section, when dealing with perturbations in longitudinal gauge, the degrees of freedom are the Newtonian potential Φ and the perturbation of the scalar field δφ for scalar perturbations and two degrees for the tensor ones, exactly the same as in GR. This is the same as what happens in mimetic gravity [45], where the degrees of freedom are the perturbed part of the mimetic field and δφ.

PERTURBATIONS
In this section we will calculate, for a given f -theory, the scalar and tensor perturbations using for scalar perturbations the longitudinal gauge (see for a review of the used setup [28,29]).

Scalar perturbations
In Newtonian gauge the the line element is given by [30] where the potentials Φ and Ψ coincide with the gauge invariant ones.
A simple calculation leads to δu 0 = Φ and δu k = ∂ k δφ φ , where δφ, in this gauge, coincides with the δφ gi (the gauge invariant perturbation of the scalar field). Then, at linear order, we have To obtain the dynamical equations for perturbations, note first of all that the perturbed i − j equation of (5) for i = j leads to the identity Ψ = Φ. Thus, perturbing the i − 0, i − i and 0 − 0 equations of (5) one gets, respectively, Adding equations 0 − 0 and i − i and using i − 0, one gets where we have introduced the notation Ω ≡ We may write equation i − 0 as:φ which leads to the following equation for the potential Φ, which, for the case of the choice of f given in (16), differs from the corresponding equation of LQC in the deformed algebra approach [31,32] only in the right hand side term, which in the last approach vanishes.
On the other hand, the equation for δφ is obtained from the linearization of the conservation equation (7) δφ + 3Hδφ − 1 From these equations, we will calculate the Mukhanov-Sasaki (M-S) equation for scalar perturbations in our approach. First of all, note that equation i − 0 can be written as On the other hand, equation 0 − 0 takes the form and, after a somewhat cumbersome calculation, one can see that it is equivalent to the following one, Introducing gauge-invariant variables (recall we are working in the Newtonian gauge) and using the conformal time v = a δφ +φ where is the curvature fluctuation in co-moving coordinates, one obtains the following M-S equations which, after inserting the second into the first one, lead to the equation for the potential (36) in the simple form while the equation for the variable v, after taking the Laplacian from the second equation and using the first one, becomes Note that in our approach, the right hand side of Eq. (44) does not vanish, meaning that the variable v, which encodes the scalar perturbations (it depends on Φ and δφ), is not independent and one needs another equation in order to calculate the evolution of the scalar perturbations. Fortunately in the matter (or matter-ekpyrotic) scenario, in the contracting phase, the pivot scale leaves the Hubble radius at rather low energy densities, as compared to the Planck one, so the corrections due to f can be safely disregarded and, thus, v satisfies approximately the usual equation v − ∆v − z z v = 0. This finally means that the calculation of the spectral quantities, such as the spectral index, its running, and the ratio of tensor to scalar perturbations, can safely be done using GR in the contracting phase, as we will show in next section.
On the other hand, to calculate the evolution of the scalar perturbations through time, we need two equations, since the equations for Φ and δφ are coupled. In this case, solving the system of Eqs. (36) and (37), we will obtain such evolution. The best suited variables for that are u and δσ ≡ aδφ, with the corresponding dynamical equations being We now have to look for the initial conditions for a matter or matter-ekpyrotic bouncing scenario [34,35,49]. Using that at very early times we are in the framework of GR, in Fourier space, we will have v k → e −ikτ one obtains the well-known result δσ k → e −ikτ √ 2k when τ → −∞. Thus, the asymptotic conditions at very early times are u k = i e −ikτ k √ 2k and δσ k = e −ikτ √ 2k .
A final remark is in order. In a bouncing scenario, as for instance the matter-ekpyrotic one, at early times GR holds, meaning thatf χχ could be disregarded and thus obtaining for the variable v, in Fourier space, the classical equation v k + (k 2 − a a )v k = 0. Then, when the pivot scale leaves the Hubble radius, i.e. in the long wavelength approximation (k 2 a 2 H 2 ), one can safely disregard the Laplacian terms which appear on the right hand side of Eq. (44) (see the end of section 4 in [50]). And maybe the same happens with the right hand side of our Eq. (44) because it contains a Laplacian. If so, Eq. (44), in the long wavelength approximation, will become, as usual, v k − z z v k = 0 and, for the f -theory given by (16), we will recover the results obtained in LQC using the deformed algebra approach [34,35,49]. Anyway, this has to be properly checked by solving numerically the system of Eqs. (45).
However, we want to stress that the system (45) can in fact be solved iteratively, taking the right hand side term as a perturbation (when the pivot scale is well inside the Hubble radius, which happens at very early times, the right hand side term can be dismissed, since GR holds, and after the pivot scale leaves the Hubble radius, i.e. in the long wavelength approximation, the Laplacian terms can be dismissed too). In fact, at very early times, disregarding the right hand term, since the universe is matter-dominated, the equation for u is, in Fourier space, Its solution satisfying the asymptotic condition where H Note that, in the case of the matter bounce scenario, given by the potential V = 2ρ c e − √ 3φ (1+e − √ 3φ ) 2 , and choosing the solution (18) which leads to the background (19) mimicking the same background as a pressureless fluid, one can see that the function θ θ , which is symmetric with respect to the contracting and expanding phase, satisfies that in the contracting regime is increasing from ρ = 0, where it vanishes, to ρ = ρc 2 , and decreasing from ρ = ρc 2 to the bounce, which occurs at ρ = ρ c , and where one has θ θ = 5 4 ρ c ∼ = 315, which, as we will see in Section VI, is bigger than the pivot scale k * ∼ 5 × 10 −33 a E ρ c . Then, the pivot scale, which occurs at a very early time, namely t = −t H , reenters to the Hubble radius at time t = t H , so in this case the equation (49) holds from −t H to t H and the equation (47) holds at very early and late times.
The solution of (49) is, then, where the sub-index I refers to an early time where GR holds, and the coefficients C 1 (k) and C 2 (k) are obtained by matching both expressions of u k in (48) and (50) when k 2 τ 2 1.
Once this solution has been obtained, we consider the i − 0 equation, which can be written as follows, Inserting Φ (0) k on it, one gets δσ (0) k . Then, to find the new iteration of u (1) k , one has to solve the equation using the well-known method of variation of constants for second-order differential equations (see for example Chapts. 13 − 17 of [51]). When one has u (1) k one calculates Φ (1) k and, inserting the result in (51), one gets δσ (1) k . And, thus, the new iterations are iterativelly obtained and a fully fledged method is constructed.
Alternatively, one could also deal directly with Eq. (44), considering the right hand side as a perturbation. In Fourier space, at very early times this equation becomes whose solution satisfying the corresponding asymptotic is When the pivot scale leaves the Hubble radius one can use the longwavelenght approximation (for the matter bounce scenario using the background (19) the function z z is symmetric with respect to the expanding and contracting phase, being an increasing function in the whole contracting period. Thus, if the pivot scale leaves the Hubble radius at whose solution is [37] v (0) At early times, if one deals with the matter bounce scenario with the background (19), this expression is equal to v τ , which one has to match with (54) so as to obtain When one has the expression of v k . Then, inserting this expression into Eq. i − 0, one obtains a solvable first-order differential equation in δσ (0) k . To obtain the next iteration one has to solve, using the method of variation of constants for second order differential equations, the equation Once v (1) k has been obtained, one has to use the definition of v and Eq. i − 0 to find δσ (1) k , and so continue successively to obtain the next iteration.
For example, we take the background (19). Since from the equation (51) one can see that the function δσ to regularize the integral one has to add an imaginary part to the poles obtaining t + + i and t − − i , then after an integration by parts one will have v (1) due, as we will show at the end of Section VI, to the small value of the pivot scale. One can also see it taking into account that Then, since the pivot scale k * leaves the Hubble radius at time −t * , when holonomy corrections could be disregarded, that is, when −t * 1 √ ρc , which is equivalent to k * = a(t * )H(t * ) √ ρ c , one can see that the dominant term in (59) is . Therefore, we have obtained approximately the same value as in the first iteration, which coincides with the result of LQC [37,49], meaning that in this iteration one will have

Tensor perturbations
The metric for tensor perturbations is [30] where h ij is a symmetric, traceless and transverse tensor ( Since the field φ does not affect the tensor perturbations, it is clear that in this theory this equation will coincide with the corresponding one in GR, that is, [30] Denoting h j i by h and introducing the variable v T = ah, in Fourier space this equation becomes

COMPARISON WITH OTHER MODELS
We will now compare our approach with a number of other different models, all of them sharing as a background the same as for holonomy corrected LQC. These models are: LQC in the deformed algebra approach, the teleparallel LQC approach, the intrinsic curvature LQC approach and the mimetic LQC approach. First, we shall review the dynamical equations for each of these theories: 1. LQC in the deformed algebra approach [31,32].
In this case the equation for the potential Φ decouples. It is given bÿ The variable M-S variable v adquires the usual form v = a(δφ +φ H Φ) and the M-S equations decouple in the simple form In particular the equation v − Ω∆v − z z v = 0 can be solved when the pivot scale is well inside and outside of the Hubble radius, obtaining the whole evolution of the variable v, which encodes all the information about scalar perturbations and, thus, the knowledge of the power spectrum for the curvature fluctuation in co-moving coordinates.
On the other hand, for tensor perturbation the corresponding M-S equation is [33] h where z T ≡ a √ Ω .
In teleparallel LQC the equation for the Newtonian potential also decouples where the square of the velocity of sound is c 2 , which, contrary to what occurs in LQC in the deformed algebra approach, is always positive, meaning that in this approach there are no gradient inestabilities. In the same way, the M-S also decouple and the only difference with the ones of LQC in the deformed algebra approach is that the velocity of sound is the same that appears in (68). As in LQC in the deformed algebra approach, the equation for the variable v decouples, which allows us to know the power spectrum of the curvature fluctuations in co-moving coordinates.

For tensor perturbations the velocity of sound is equal to 1, and the M-S equation is given by
where z T ≡ a ρc 12H 2 arcsin  3. Intrinsic curvature LQC [14].
In this approach the equations for scalar perturbations are the same as in LQC in the deformed algebra approach. However, for tensor perturbations, the corresponding M-S equation is where z T is the same as in teleparallel LQC, and in this case the square of the velocity of sound is given by

Mimetic LQC [45].
This is also a fully covariant approach where, at the perturbative level, the dynamical variables are the perturbed mimetic field, namely δϕ, and the perturbed scalar field δφ. In the Newtonian gauge, the potential is related with the mimetic field, as follows Φ = δφ, and the dynamical equation, as in our approach, becomes coupled [45] δφ + Hδφ − where the square of the velocity of sound is given, as in [46], by which exhibits the well-known gradient instability of the mimetic gravity case [46,47] (see also [48] for the study of perturbations in specific mimetic matter models).
Dealing with tensor perturbations, since the mimetic field does not alter the gravitational sector, as in our approach, the equations for tensor perturbations are the same as in GR.
Consequently, we can see that for non-fully covariant theories, such as LQC in the deformed algebra approach, teleparallel LQC, or intrinsic curvature LQC, the equations for scalar perturbations decouple, which actually simplifies the theory a lot, allowing us to calculate the corresponding power spectrum. On the contrary, for the fully covariant theories, i.e. our approach and mimetic LQC, the equations of scalar perturbations do not decople, which makes their analytic study more difficult and only a numerical analysis seems to be viable in order to understand their evolution. However, the clear advantage of the covariant theories is that the equation for tensor perturbations is the simplest one because it coincides with the one for tensor perturbations in GR.

REHEATING AND THE CALCULATION OF THE SPECTRAL PARAMETERS IN THE BOUNCING MATTER-EKPYROTIC SCENARIO
In this section we consider the background given by holonomy corrected LQC. In other words, we consider the function f given by (16), which means that the universe bounces when its energy density is ρ c . We will show how to calculate the reheating temperature of the universe via gravitational particle production due to a phase transition from the matter domination to an ekpyrotic era, and how the theoretical values of the spectral index and its running match well with the corresponding observational data.
As we will immediately show, in a viable bouncing scenario the pivot scale leaves the Hubble radius in the contracting phase, when GR does hold. Correspondingly, the spectral index and its running are given, respectively, by [52,53] where w is the effective Equation of State (EoS) parameter and the "star" means that the quantities are evaluated when the pivot scale leaves the Hubble radius.
First of all, note that in order to match the theoretical value of the spectral index with the observational data, w * has to be negative and close to zero. Now, assuming that, at very early times, we have a quasi-matter domination epoch (φ 2 ∼ = 2V →φ ∼ = V φ ), the background equations become In this case and, thus, Note that for a potential corresponding exactly to matter domination, i.e. for V = λe √ 3φ , one obtains an exactly flat spectrum (n s = 1 and α s = 0).
We consider here a phase transition, in the contracting phase, from matter domination to an ekpyrotic regime. To this end, we choose a model with potential given by where k is a positive parameter and λ andλ satisfy since, at the end of the phase transition, we assume |λ| ρ c ∼ = 252. This models depicts for φ → −∞ a matter dominated universe, and for φ ≥ 0 an ekpyrotic universe with EoS parameter w = 2 [54].
Then, for φ negative satisfying e φ/k 1, one has Using now the BICEP2/Keck Array and Planck observational data at 1σ C.L., one has n s = 0.968 ± 0.006 and α s = −0.003 ± 0.007 [55]. Then, choosing for instance k = 10, one can see that for −42.38 ≤ φ * ≤ −24.49 the theoretical values of the spectral index and its running belong to the 1-dimensional marginalized 2σ C.L..

Reheating
Here we will consider a reheating process due to the gravitational particle production of massless particles minimally coupled with gravity during a phase transition from a matter dominated regime to an ekpyrotic one in the contracting phase. Recall that in our model (77) we have an ekpyrotic phase with EoS parameter w = 2. Let H E be the value of the Hubble parameter at the beginning of this phase. Then, from the relations V = 3H 2 +Ḣ andḢ = − 9 2 H 2 , we obtain λ = 3 2 H 2 E . On the other hand, in holonomy corrected LQC a viable value of H E is approximately −10 −3 [34,35], which justifies our choice |λ| ρ c .
The energy density of the produced particles during this phase transition is given by [56] where a E is the value of the scale factor at the beginning of the ekpyrotic phase and R ∼ 10 −2 N s , being N s the number of different scalar fields.
Remark 6.1 The equation (80) was obtained considering a phase transition, in the expanding phase, from the de Sitter phase to another one with constant EoS parameter w > 1/3 [57,58]. In the case we consider a phase transition in the contracting phase from the matter domination phase to another one with constant EoS parameter w > 1, this formula is also valid due to the duality, pointed out in [59], between the de Sitter regime in the expanding phase and the matter domination in the contracting one.
On the other hand, as has been showed numerically in [34,35] (see the figures 11, 13 and 15 of [35]), in the matterekpyrotic bounce scenario, after the bounce the universe enters in a kination regime, i.e., the effective EoS parameter is equal to 1. Then, to clarify ideas we consider the background equations corresponding to LQC (17) and we consider a fluid with the following EoS, in the contracting phase, and P (ρ) = ρ in the expanding one.
Since for the EoS P = wρ in LQC the Hubble parameter evolves as , for our EoS (81) we will have where t E is the phase transition time and, thus, it has to satisfy 3 2 ρct E 27 4 t 2 E +1 = H E , andt has to be chosen imposing continuity at t = t E , that is, it has to satisfy Remark 6.2 Since in LQC when considering a fluid with EoS P = wρ the reconstruction method (see Section 3 of [35] for a detailed discussion) leads to the potential then the conservation equation with the following potential where φ E and φ B satisfy respectively 2V (φ E ) = ρ E and −2V (φ B ) = ρ c , has a solution which leads to the background (82). a E a 6 . Then, the universe will become reheated when both energy densities are of the same order, i.e. ρ r (a reh ) ∼ ρ b (a reh ). Since in the contracting phase the energy density of the background increases faster than the one of the produced particles, the reheating will occur in the expanding , and it will be

Now, since for a fluid with the linear Equation of State
Then, for N s ∼ 1 (GUT theories), ρ c ∼ = 252 and H E ∼ −10 −3 one obtains T reh ∼ 8 × 10 −10 . Finally, since M pl ∼ = 2.4 × 10 18 GeV, in natural units one has T reh ∼ 2 × 10 9 GeV. Remark 6.3 It is possible to obtain a lower reheating temperature by increasing the EoS parameter in the ekpyrotic phase. For example, if in the ekpyrotic phase one takes w = 5, then the reheating temperature is reduced by one order.
Once we have calculated the reheating temperature we can show that the pivot scale, in the contracting phase, leaves the Hubble scale when GR holds. The pivot scale is related with its physical value by k * = a 0 k phys (t 0 ), where the sub-index 0 means present time, and we choose, as usual, k phys (t 0 ) ∼ 10 2 H 0 ∼ 10 −59 .
On the other hand, as we have showed at the reheating time, i.e., when both energy densities are of the same order, one has a reh ∼ ρc 3H 2 E ρc 2 3 RH 4 E a E ∼ 5 × 10 4 a E . Now, from the conservation of the entropy we have the adiabatic relation a 0 ∼ T reh T0 a reh [60] and using that the current and reheating temperature are respectively T 0 ∼ 8 × 10 −32 and T reh ∼ 8 × 10 −10 , one gets a 0 ∼ 5 × 10 26 a E . As a consequence, k * ∼ 5 × 10 −33 a E , which means, since |H E | ∼ 10 −3 , that k * |H E |a E . In other words, the pivot scale leaves in the contracting phase the Hubble radius well after the phase transition, more precisely when H * ∼ 5 × 10 −33 a E a * ≤ 5 × 10 −33 (in the contracting phase a E < a * ) and, thus, since ρ c ∼ = 252, one can safely disregard the effects of the f theory when the pivot scale leaves the Hubble radius.

CONCLUSIONS
We have constructed a class of modified gravitational theories based on the addition to the Einstein-Hilbert action of a function f , which depends on the divergence of the unitary time-like eigenvector of the stress tensor. We have obtained in this way a fully covariant theory, which, as in the case of mimetic gravity, has one more degree of freedom than GR. The main advantage, at the background level, of our class of models is that, for the FLRW geometry, this divergence is minus three times the Hubble parameter, which allows, by choosing the function f appropriately, to obtain very simple bouncing backgrounds, as the one obtained in holonomy corrected LQC.
At the level of cosmological perturbations, working in the Newtonian gauge, the equations for scalar perturbations exhibit some of the same interesting features as those appearing in LQC in the so-called deformed algebra approach. However, they are not exactly the same, owing to the fact that our theory is fully covariant, in contrast with LQC in the deformed algebra approach [61][62][63]. In fact, contrary to what happens with LQC and other non-covariant approaches such as teleparallelism or modified theories using the extrinsic curvature, in our fully covariant approach the equations do not decouple, which is an added difficulty and ammounts to the fact that, in practice, the equations cannot be solved analytically but only by standard numerical methods.
A very positive feature is, however, that for tensor perturbations our model leads to the same equations as GR because the modification of the action does not affect the gravity sector.
Finally, we have studied the matter-ekpyrotic bouncing scenario for the LQC background, when reheating is a consequence of the production of massless particles minimally coupled with gravity, during the phase transition from matter domination to the ekpyrotic regime in the contracting phase. We have obtained a viable reheating temperature of the order of 10 9 GeV and have shown that the observable modes leave the Hubble radius in the contracting phase, when the holonomy correction can be disregarded. This permits a very simple calculation of the theoretical values of the spectral index and of its running, both of which, as it turns out, perfectly match the current observational data at the 2σ C.L..