Bouncing cosmology in $f(R,\mathcal{G})$ gravity by order reduction
Bruno J. Barros Elsa M. Teixeira Daniele Vernieri
gr-qc astro-ph.CO hep-th








Instituto de Astrofísica e Ciências do Espaço, 
Faculdade de Ciências da Universidade de Lisboa, Campo Grande, PT1749-016 Lisboa, Portugal



A bouncing universe is a viable candidate to solve the initial singularity problem. Here we consider bouncing solutions in the context of f(R,𝒢) gravity by using an order reduction technique which allows one to find solutions that are perturbatively close to General Relativity. This procedure also acts as a model selection approach. Indeed, several covariant gravitational actions leading to a bounce are directly selected by demanding that the Friedmann equation derived within such gravity theories coincides with the one emerging from Loop Quantum Cosmology.


Bouncing cosmology in f(R,𝒢) gravity by order reduction
    Daniele Vernieri
    December 30, 2023
=======================================================





§ INTRODUCTION


An effective description of Loop Quantum Cosmology (LQC) <cit.> enables one to write cosmological dynamics in a classical fashion incorporating leading order quantum corrections. While at low energies the predictions of LQC are very close to General Relativity (GR), in high density regimes, where quantum effects become dominant, these might drastically depart from the classical description <cit.>. One such phenomenon is the Big Bang singularity, which in light of LQC is replaced by a quantum bounce <cit.>. This becomes evident when writing the modified Friedmann equation emerging from the effective description of LQC <cit.>:


    H^2 = κ/3ρ( 1-ρ/ρ_c),


where H := ȧ/a is the Hubble rate, a(t) is the scale factor of the universe, κ≡ 8π, and ρ_c = √(3)/(32π^2γ^3) gives the critical energy density for the bounce, with γ being the Barbero–Immirzi parameter. Notice that we are using physical units such that the Newton's gravitational constant G_N, the speed of light c and the reduced Planck constant ℏ have been set to one.

We will consider bouncing scenarios in the context of modified gravity theories <cit.> whose gravitational action contains higher-order derivatives and that generically propagate extra degrees of freedom. In the specific case of the generalized Gauss–Bonnet f(R,𝒢) gravity theories, that will be the main focus of this work, R being the Ricci scalar and 𝒢 the Gauss–Bonnet term, one can indeed write equivalently the corresponding higher-order gravitational action as a scalar–tensor theory with two extra dynamical scalar fields <cit.>. 

These generalized Gauss–Bonnet gravity models have been considered to describe cosmological phenomena, such as the inflationary epoch <cit.> and dark energy scenarios <cit.>. Future singularities within f(R,𝒢) gravity have been discussed in Ref. <cit.>, together with some proposals to avoid them, such as considering a Lagrangian with terms of the form 𝒢^α/R^β. The study of the stability conditions for scalar and tensor perturbations in f(R,𝒢) gravity was conducted in Refs. <cit.>. Primordial inflation was also studied in Ref. <cit.>, where it was shown that a double inflationary setting naturally arises in the context of such theories. Bouncing scenarios and cosmological applications of f(𝒢) models have been studied in Refs. <cit.>.

In order to avoid the presence of spurious degrees of freedom that are usually present in such theories, we will implement the so-called order reduction technique <cit.>. The latter consists in reducing the field equations by expressing the geometrical quantities in terms of the matter fields at the lowest order in the perturbative expansion. Hence, the gravity theory at hand is framed within an effective field theory approach and the reduced field equations contain derivatives up to second order. Then, only the solutions that are perturbatively close to GR are selected.

This paper addresses the following questions: Can the Friedmann equation emerging from LQC be derived from an effective modified gravity theory? And, if so, what is the most general form for the Lagrangian under these conditions? Indeed, this was already proved to be possible for f(R) <cit.> and f(𝒢) <cit.> gravity. In the present manuscript we approach this issue in the context of f(R,𝒢) gravity models. We find novel solutions for the form of the gravitational Lagrangian allowing for bouncing cosmological scenarios by means of the order reduction technique. Moreover this work extends the results reported in Refs. <cit.>. Notice that previous works have been conducted to reproduce the Friedmann equation from LQC but using different methods than the one employed in this work and within distinct gravity theories. In particular, in the framework of Palatini f(R) gravity <cit.>, by demanding that the Friedmann equation matches the one from LQC, the specific form of the f(R) Lagrangian was numerically reconstructed by means of interpolation techniques. Moreover, the covariant action for a class of higher-derivative scalar–tensor theories, referred to as mimetic gravity, was reconstructed in Ref. <cit.> such that the effective dynamics of LQC is reproduced within homogeneous and isotropic backgrounds, resorting to an Hamiltonian analysis. 

The structure of the paper is the following. In Sec. <ref> we write the full action for the gravitational sector and for the matter fields, the modified field equations, and then apply the order reduction technique. In Sec. <ref> we follow to specify the background cosmology and write the evolution equations accordingly. The results are presented and discussed in Sec. <ref>. Finally we conclude in Sec. <ref>.



§ F(R,𝒢) GRAVITY AND ORDER REDUCTION



We begin by writing down the action for f(R,𝒢) modified gravity theories in a 4-dimensional spacetime,


    𝒮= 1/2κ∫ d^4 x √(-g)[ R + f(R,𝒢)  ] + 𝒮_m (g_μν,ψ) ,


where g is the determinant of the metric tensor and the gravitational sector consists of the standard Einstein–Hilbert term plus an additional general function of the Ricci scalar curvature and the Gauss–Bonnet term, the latter defined as


    𝒢 = R^2 - 4R^μνR_μν + R^μναβR_μναβ ,


with R_μν and R_μναβ being the components of the Ricci and Riemann tensors, respectively. The deviations from GR are generically enclosed in the function f(R,𝒢). The object 𝒮_m in Eq. (<ref>) portrays the matter sector, with ψ standing for a collective representation of the matter fields. 

Following the standard procedure of varying the action in Eq. (<ref>) with respect to (wrt) the metric g^μν, we arrive at the following modified field equations:


    R_μν - 1/2g_μνR - 1/2g_μν f + f_,RR_μν+ ( g_μν □  - ∇_μ∇_ν)f_,R 
    + f_,𝒢( 2RR_μν - 4R_μ^    αR_να + 2R_μαβγR_ν^    αβγ  -4R_μανβR^αβ) 
    
    + (2 g_μνR □ -2R ∇_μ∇_ν +4R_μ^    α ∇_α∇_ν+4R_ν^    α ∇_α∇_μ - 4g_μνR^αβ ∇_α∇_β - 4R_μν □ + 4R_μ    ν^    α    β ∇_α∇_β)f_,𝒢
    =κ T_μν ,


where ∇_μ is the covariant derivative and □ = g^μν∇_μ∇_ν is the D'Alembertian operator. We have defined f_,R := ∂ f/ ∂ R, f_,𝒢 :=  ∂ f/ ∂𝒢, and


    T_μν= -2/√(-g)δ𝒮_m/δ g^μν


is the matter stress-energy tensor.

We rely on the fact that the Riemann tensor can be decomposed into a traceless, semi-traceless and a scalar part, respectively as:


    R_μναβ =         C_μναβ  + 1/2( g_μαR_βν+g_νβR_αμ- g_μβR_αν  -g_ναR_βμ) +1/6( g_μβg_αν - g_μαg_βν) R ,


with C_μναβ being the Weyl tensor components. As a first simplification, we require that the Weyl tensor vanishes: C_μναβ=0. This amounts to assuming that our spacetime is conformally flat, i.e., every point in our spacetime manifold contains a neighborhood which can be mapped to a flat space by means of a conformal transformation. This condition always holds for a Friedmann–Lemaître–Robertson–Walker (FLRW) background on which we will base our study throughout the next sections. Thus, the Riemann tensor components reduce to


    R_μναβ =        1/2( g_μαR_βν+g_νβR_αμ - g_μβR_αν -g_ναR_βμ) +1/6( g_μβg_αν - g_μαg_βν) R ,


and the Gauss–Bonnet term 𝒢 can be written as


    𝒢 = 2/3R^2 - 2R_μνR^μν .


In order to avoid the presence of spurious degrees of freedom that usually arise in these theories, we employ the order reduction technique <cit.> by introducing a parameter ϵ such that its vanishing promptly leads us back to GR. Thus, following the procedure presented in Refs. <cit.>, we parametrize the function f(R,𝒢) as


    f(R,𝒢) = 2Λ + ϵφ(R,𝒢) ,


where Λ is a cosmological constant term and the dimensionless parameter ϵ regulates the deviations from GR. This is consistent with the fact that we wish to frame f(R,𝒢) gravity within an effective field theory approach such that its solutions are perturbatively close to GR, namely that ϵφ≪ R. Accordingly, we will only consider the lowest order terms in ϵ. By doing so, we constrain the higher order field equations to contain up to second order derivatives of the metric, thus avoiding the propagation of additional degrees of freedom arising in such theories.

To apply the order reduction technique we simply replace all the 𝒪(ϵ) terms in the modified field equations by their lowest order (ϵ =0) counterpart, that ends up with expressing the geometrical quantities R, R_μν, R_μναβ, and 𝒢, in terms of the matter fields at the lowest order in the perturbative expansion. Indeed, evaluating Eq. (<ref>) and its trace, with f as defined in Eq. (<ref>), for ϵ=0, yields the following expressions for the order reduced curvature scalar and Ricci tensor components:


    R^T    =    -κ T -4Λ , 
     
    R^T_μν   =   κ T_μν-κ/2g_μνT - Λ g_μν ,


where the superscript T denotes the order reduced quantities. Based on the above relations, together with the assumption of vanishing Weyl tensor, the Riemann tensor components in Eq. (<ref>), and the Gauss–Bonnet invariant in Eq. (<ref>), can be expressed at the lowest order as follows:


    R^T_μναβ =         -κ/2( g_μβT_αν+g_ναT_βμ - g_μαT_βν-g_νβT_αμ) -1/3( g_μαg_νβ-g_μβg_αν) ( Λ+κ T ), 
    𝒢^T =         2/3κ^2 T^2 - 2 κ^2 T_μνT^μν + 8/3Λ^2 + 4/3Λκ T .


Finally, we can formally write the reduced field equations, with each order-ϵ term replaced by its corresponding order reduced analogue,


    R_μν - 1/2g_μνR- g_μνΛ +ϵ[ -1/2g_μνφ^T + φ^T_,RR^T_μν+ ( g_μν □  - ∇_μ∇_ν)φ^T_,R 
    + φ^T_,𝒢( 2R^T R^T_μν - 4R_μ^T   αR^T_να.  
    . + 2R^T_μαβγR_ν^T   αβγ -4R^T_μανβR^T  αβ)+ (2 g_μνR^T □ -2R^T ∇_μ∇_ν +4R_μ^T  α ∇_α∇_ν +4R_ν^T  α ∇_α∇_μ - 4g_μνR^T  αβ ∇_α∇_β. 
    . - 4R^T_μν □ + 4R_μ    ν^Tα    β ∇_α∇_β)φ_,𝒢] =κ T_μν ,


in which, accordingly, the derivatives are wrt the lowest order quantities, i.e., φ^T_,R := ∂φ^T / ∂ R^T, with φ^T := φ^T (R^T,𝒢^T). 

Equation (<ref>) is then simply written as GR plus order-ϵ terms which will dominate near the bounce. The assumption ϵφ≪ R is a necessary condition for the whole scheme to be valid <cit.>, guaranteeing that the corrective terms are treated within an effective field theory approach, where R≪ρ_c ∼ l_p^-2 (l_p being the Planck scale) as will be verified a posteriori through the bouncing solutions presented in Sec. <ref>.

We proceed with the implementation of this model to a particular geometrical setting, specifying a form for g_μν and writing the cosmological equations in this context.



§ ORDER REDUCTION WITHIN A FLRW GEOMETRY



We will now consider the Friedmann–Lemaître–Robertson–Walker (FLRW) cosmological background, with line element given by


    ds^2 = -dt^2 +a(t)^2 [ dr^2/1-kr^2 +r^2( dθ^2 + sin^2θ  dϕ^2 ) ].


The constant k∈{-1,0,1} determines the spatial curvature of the universe (hyperbolic, flat or spherical, respectively).

Furthermore, we assume a perfect fluid form for the matter stress-energy tensor, namely


    T_μν = ( ρ + p )u_μu_ν + p g_μν ,


where u_μ is the fluid's 4-velocity, defined in such a way that u_μ u^μ=-1, and ρ and p denote the energy density and pressure, respectively. We consider a barotropic equation of state (EoS) for the pressure, i.e.,


    p = wρ ,


with w∈[-1,1] being a constant.

The time component of the conservation equation, ∇_μT^μν=0, yields the continuity equation


    ρ̇= -3H ( 1+w )ρ,


where a dot denotes derivative wrt cosmic time.

Using Eqs. (<ref>)-(<ref>), the order reduced scalars R^T and 𝒢^T read:


    R^T    =   κρ( 1-3w ) -4Λ ,  
    𝒢^T    =   4/3( Λ - κρ) [ 2Λ + κρ( 1+3w ) ].


The modified Friedmann equation is computed from the 0–0 component of the reduced field equations in Eq. (<ref>):

    6H^2 + 6k/a^2 -2Λ+2κρ+ ϵ{φ^T +6 H φ̇^T_R +φ^T_R [2 Λ +κρ (1 +3  w) ]-8 H φ̇^T_𝒢 ( Λ -κρ) 
    
     -4/3φ^T_𝒢 (Λ -κρ) [2 Λ +κρ (1 +3  w)]} =0 .


It is convenient to use the chain rule, alongside the continuity equation, Eq. (<ref>), and Eqs. (<ref>)-(<ref>), to express the time derivative of φ^T in terms of derivatives wrt 𝒢^T and R^T: 


    φ̇^T = ∂φ^T/∂ R^T∂ R^T/∂ρρ̇ + ∂φ^T/∂𝒢^T∂𝒢^T/∂ρρ̇ = -3 φ^T_R H(1+w)(1-3w)κρ +4 φ^T_𝒢 H(1+w)[ Λ (1-3w)  + 2κρ (1+3w) ]κρ .


The latter can be substituted in Eq. (<ref>) in order to finally arrive at the expression of the reduced modified Friedmann equation, with a cosmological constant Λ, spatial curvature k and EoS parameter w:


    H^2=κ/3ρ - Λ/3 - k/a^2
    
    + ϵ/18{6k/a^2κρ  (1+w)  [16 φ^T_𝒢𝒢 (κρ -Λ) [Λ  (1-3w)+2 κρ  (1+3w)]+12 φ^T_R 𝒢[κρ  (1+9w)+2Λ (1-3 w)]-9 φ^T_R R
       (1-3w)]  
    
        -  3 φ^T +2 (Λ -κρ ) (2 φ^T_𝒢[2 Λ +κρ  (1+3 w)]+κρ  (1+w) [16 φ^T_𝒢𝒢 (κρ -Λ) [Λ  (1-3 w)+2
       κρ  (1+3 w)]  
    
       +12 φ^T_R 𝒢[ κρ  (1+9 w)+2 Λ  (1-3 w)]-9 φ^T_R R (1-3w)]) -3 φ^T_R[2 Λ +κρ  (1+3 w)] }.

In the equation above we have also replaced H^2 by its zero-th order counterpart, i.e., H^2=1/3κρ -k/a^2-Λ/3.



§ BOUNCING COSMOLOGY



Particularizing the cosmology to a flat spacetime (k=0) (in which case the second line of Eq. (<ref>) vanishes) without a cosmological constant (Λ=0), Eq. (<ref>) can be recast into:


    H^2=κρ/3 -ϵ/3[ 1/2φ^T +2/3(1+3w) κ ^2 ρ ^2  φ^T_𝒢 + 32/3(1+3w)(1+w) κ ^4 ρ ^4  φ^T_𝒢𝒢 +1/2(1+3w) κρφ^T_R
    
    -3(1-3w)(1+w)κ ^2 ρ ^2 φ^T_RR +4(1+9w)(1+w) κ ^3 ρ ^3 φ^T_R𝒢],


and the scalar curvature invariants given by Eqs. (<ref>) and (<ref>) become, respectively:


    R^T    =   (1-3w ) κρ ,  
    𝒢^T    =    -4/3( 1+3w )κ^2ρ^2 .


As we wish to account for the existence of solutions within the framework of f(R,𝒢) gravity, it seems reasonable to look only for solutions where both the curvature scalar R and the invariant 𝒢 do not vanish, otherwise our treatment simply reproduces f(R) or f(𝒢) gravity models. Nonetheless, we will comment on how the general solutions at hand reduce to these particular scenarios. Thus, hereafter we will assume w ≠± 1/3. It is worth mentioning that radiation-like fluids are therefore automatically excluded.

The main purpose of this work is to find an analytical form for the function f(R,𝒢) such that we recover the effective Friedmann equation of a bouncing universe from LQC, which reads


    H^2 = κ/3ρ( 1-ρ/ρ_c).


Since the quantity ρ_c defines the characteristic scale of the problem, and also for dimensionality reasons, it is convenient to define the following dimensionless quantities <cit.>: 


    R̅^T=R^T/ρ_c ,   𝒢̅^T= 𝒢^T/ρ_c^2 ,   ρ̅= ρ/ρ_c  and  φ̅^T=φ^T/ρ_c .


Then, we can directly equate the right-hand side of the Friedmann Eq. (<ref>), to the one we wish to reproduce, that is Eq. (<ref>). This amounts to a differential equation for φ̅^T(R̅^T,𝒢̅^T) expressed in terms of the dimensionless variables:


    1/2φ̅^T +2/3(1+3w) κ ^2 ρ̅ ^2  φ̅^T_𝒢̅ + 32/3(1+3w)(1+w) κ ^4 ρ̅ ^4  φ̅^T_𝒢̅𝒢̅ +1/2(1+3w) κρ̅φ̅^T_R̅ -3(1-3w)(1+w)κ ^2 ρ̅ ^2 φ̅^T_R̅R̅
    
     +4(1+9w)(1+w) κ ^3 ρ̅ ^3 φ̅^T_R̅𝒢̅ = κρ̅^2/ϵ ,


which is a non-linear second-order partial differential equation for a two variable function. We proceed with an empirical approach, by proposing an ansatz for the f(R,𝒢) function (to be more rigorous, for φ̅^T) as a specific solution of Eq. (<ref>). Afterwards, we investigate the conditions under which the Lagrangian found is capable of reproducing a bouncing universe. We allow the EoS parameter w for the fluid to be arbitrary while in previous works regarding f(R) <cit.> and f(𝒢) <cit.> gravity only the case w=1 was considered.

In the following, we omit the order reduction superscript T for notational convenience. 



 §.§ Solution I



The simplest morphology for φ̅ as a candidate for the solution of the differential equation, Eq. (<ref>), is a power-law form:


    φ̅ (R̅,𝒢̅)= c_1 |𝒢̅|^α|R̅|^β ,


where α, β and c_1 are real constants. From Eqs. (<ref>) and (<ref>) we gather that, depending on the value of w, the scalars R and 𝒢 may become negative. Therefore, in order to avoid the emergence of complex solutions, and without loss of generality, we choose to write the ansatz in terms of the modulus of the scalars, with powers {α,β}∈ℝ (this however becomes superfluous for {α, β}∈ℤ). Particular cases of this form have been widely considered in the literature to study the stability of cosmological solutions <cit.>, and to describe both inflationary <cit.> and dark energy <cit.> scenarios. The particular case (α, β) = (0, 2) has already been shown <cit.> to successfully reproduce a cosmological bounce for an f(R) gravity scenario under the same order reduction technique employed in the present manuscript.

Inserting the expression for φ̅ given in Eq. (<ref>) directly into the left-hand side of Eq. (<ref>), and making use of the identities in Eqs. (<ref>)-(<ref>), it follows that:


    ρ̅^2α+β∝ρ̅^2 .


In order for the relation above to be consistent we find that 2α+β = 2 must hold, or equivalently, β=2-2α. Moreover, the Friedmann equation of a bouncing universe emerging from LQC, that is Eq. (<ref>), is reproduced from the requirement that the constant c_1 reads:


    c_1 = 2/93^α/4^α1/(1-α)ϵκ|1-3w|^2α|1+3w|^-α/(w+1)(3w-1) .


Hence, we find a family of solutions of the form,


    ℒ_g = R + C_1/ρ_c |𝒢|^α|R|^2-2α ,


with C_1=ϵ  c_1. The particular cases of α=1 or w=-1 are excluded as they do not yield a cosmological bounce. When α=0, the result found for f(R) gravity (with w=1) in Ref. <cit.> is extended to arbitrary values of w:


    ℒ_g = R+2/9 (w+1)(3w-1) κρ_c  R^2 .

Recall that the values w=±1/3 were ruled out a priori. The particular case w=1/3 yields R^T = 0, from which we require that α=1 so that φ≠ 0. However, this does not feature a cosmological bounce since it leads to φ∝ G, which is just a topological term. On the other hand, w=-1/3 (𝒢^T = 0) is indeed valid provided that α=0 and, according to Eq. (<ref>), the correction to the Friedmann equation stems from an R^2 term.
 
It is possible to verify that by taking the ansatz φ̅(R̅,𝒢̅)=R̅Ψ (𝒢̅), and solving Eq. (<ref>) for Ψ, we obtain a solution where the bounce is induced by Eq. (<ref>) with α=1/2.



 §.§ Solution II



Having attained a solution which encapsulates the model presented in Ref. <cit.>, for an arbitrary w, at this point we aim at constructing an ansatz which enclosures the solution found in the framework of f(𝒢) gravity <cit.>. This suggests looking for φ̅ functions of the following form:


    φ̅ (R̅,𝒢̅)=c_2 𝒢̅ln(|𝒢̅|^β|R̅|^γ) ,


where β, γ and c_2 are real-valued constants. 

Following the same reasoning as in the previous case, we have opted for defining φ̅ in terms of the modulus of the scalars in order to ensure that the solutions are real-valued and well-defined. This scenario is successful in reproducing the effective Friedmann equation from LQC. We have verified that the term multiplying the logarithm in Eq. (<ref>) must be exactly 𝒢 (for w≠ -1) and not some combination of different powers of R and 𝒢. This is intricately related to the form of the differential Eq. (<ref>) and the factors needed in order to cancel out the dependence on the logarithm in the Friedmann equation. The powers inside the logarithm have no influence whatsoever on the exponent of ρ in the modified Friedmann equation since Eq. (<ref>) can be manipulated such that they appear merely as multiplicative factors.

Directly inserting the expression for φ̅, Eq. (<ref>), in the differential equation (<ref>) yields:


    c_2 = 3/21/ϵκ1-3w/β(-11+24w+27w^2)+γ(-1+30w+27w^2) .


The parameters β, γ and w must be such that the denominator of the last expression is different from zero. The gravitational Lagrangian found is:


    ℒ_g = R+C_2/ρ_c 𝒢ln|𝒢|^β|R|^γ/ρ_c^2β+γ  ,


with C_2 = ϵ  c_2. Solving Eq. (<ref>) for the case φ̅ (R̅,𝒢̅)=𝒢̅Ψ(R̅) yields a solution where the bounce is driven by the corrective term in Eq. (<ref>) with β=0. It is also worth mentioning that the choice γ=0 corresponds to a generalization of the solution found for f(𝒢) gravity in the case of w=1 <cit.> to arbitrary values of w:


    ℒ_g = R-3/2 (11+9w) κρ_c𝒢ln|𝒢|/ρ_c^2 .


Note that the case w=1/3 (R^T = 0) is valid for γ = 0 which leads to f=f(𝒢) governed by Eq. (<ref>). On the other hand, the choice w=-1/3 (𝒢^T = 0) is obviously not allowed by the ansatz in Eq. (<ref>). As a final remark we briefly mention that, for the particular choice of w=-1, it is possible to have a general combination of powers of R and 𝒢 multiplying the logarithm and still obtain solutions for the differential equation (<ref>). Thus, exclusively for the case w=-1, the gravitational Lagrangian can be extended to:


    ℒ_g = R-6^δ/41/κρ_c(2β+γ) |𝒢|^δ|R|^2-2δln|𝒢|^β|R|^γ/ρ_c^2β+γ ,


where {β, γ,  δ}∈ℝ such that 2β+γ≠0 .



§ CONCLUSIONS



In this work we have reconstructed gravitational Lagrangians for f(R, 𝒢) gravity that reproduce the Friedmann equation of a bouncing universe, emerging from an effective description of LQC. 

We have formally derived covariant actions for f(R,𝒢) gravity on the basis of an effective field theory approach, which has been implemented by virtue of the order reduction technique. As an artifact of this procedure we were able to find solutions which are perturbatively close to GR while avoiding the, possibly ghost-like, spurious degrees of freedom, generically arising in these theories. Such an approach can lead to much different and richer phenomenological implications than GR such as, in this particular case, to bouncing solutions found in LQC formulations, at early times and near to the Planck scale. 

In this framework, the most general gravitational Lagrangian found can be written as:


    ℒ_g    =    R+𝒜/ρ_c |𝒢|^α|R|^2-2α+ℬ/ρ_c 𝒢ln|𝒢|^β|R|^γ/ρ_c^2β + γ ,


where {𝒜,ℬ,α,β,γ}∈ℝ are free parameters. The Friedmann equation emerging from the Lagrangian in Eq. (<ref>) via the order reduction technique is:


    H^2 = κ/3ρ[ 1-( 𝒜/C_1+ℬ/C_2)ρ/ρ_c].


In order for the Friedmann equation (<ref>) to match the one emerging from LQC, the following constraint has to be satisfied: 𝒜/C_1+ℬ/C_2=1 (hence, only four parameters are indeed free).

Notice that, in order to reproduce the additional term ∝ρ^2 in the Friedmann equation, all our solutions satisfy ϵφ∼ R^2/ρ_c (since 𝒢∼ R^2). Thus, as previously mentioned in Sec. <ref>, it is now straightforward to verify that the necessary condition for the validity of the order reduction approach, i.e. ϵφ≪ R, leads to R ≪ρ_c, since ϵφ∼ R^2/ρ_c ≪ R  ⇒  R ≪ρ_c.

The Lagrangian in Eq. (<ref>) extends previous works concerning f(R) <cit.> and f(𝒢) <cit.> gravity, for a fluid with an arbitrary EoS parameter w. Let us stress that we have not assumed any specific form for the f(R,𝒢) gravitational action a priori, but we have found it just as a result of the order reduction technique employed. As a further extension of this work, it would be interesting to investigate if the generalized Friedmann equations obtained in LQC (see Ref. <cit.> and references therein) leading to a bounce, can be obtained in the context of f(R,𝒢) gravity by means of the same order reduction approach. 


Finally, it should be noted that since the dynamics of this model deviates from GR at high energies, the primordial power spectrum of both scalar and tensor perturbations is expected to be altered, since it directly depends on the Hubble expansion rate. Such phenomena have also been studied in the context of other theories which present deviations from the standard Friedmann expansion at early times, e.g. in Randall Sundrum II braneworld model <cit.>. Hence, in order to test the robustness of the theory it would be interesting to explore the effects stemming from Eq. (<ref>) on the power spectrum of primordial fluctuations and compare them with current observational data.

This research was supported by Fundação para a Ciência e a Tecnologia (FCT) through the research grant UID/FIS/04434/2019, by the projects PTDC/FIS-OUT/29048/2017, COMPETE2020: POCI-01-0145-FEDER-028987 & FCT: PTDC/FIS-AST/28987/2017. BJB was supported by the grant PD/BD/128018/2016 from FCT. EMT was supported by the grant BI-MESTRE-IF/00852/2015 from FCT. DV was supported by the grant BI-DOUTOR-IF/00852/2015 from FCT.



