Viable non-singular cosmic bounce in holonomy improved F(R) gravity endowed with a Lagrange multiplier

Matter and quasi-matter bounce scenarios are studied for an F(R) gravity model with holonomy corrections and a Lagrange multiplier, with a scale factor $a(t) = \left(a_0t^2 + 1 \right)^n$, where the Hubble parameter squared has a linear and a quadratic dependence on the effective energy density. Provided $n<1/2$, it is shown that the primordial curvature perturbations are generated deeply into the contracting era, at large negative time, which makes the low-curvature limit a good approximation for calculating the perturbation power spectrum. Moreover, it is shown that, for $n$ within this range, the obtained cosmological quantities are fully compatible with the Planck constraints, and that the"low curvature limit"comes as a viable approximation to calculate the power spectra of both scalar and tensor perturbations. Using reconstruction techniques for F(R) gravity with the Lagrange multiplier, the precise form of the effective F(R) gravity is found, from which one determines the power spectra of scalar and tensor perturbations in such bouncing scenario. Correspondingly, the spectral index for curvature perturbations and the tensor to scalar ratio are obtained, and these values are successfully confronted with the latest Planck observations. Further, it is shown that both the weak and the null energy conditions are satisfied, thanks to the holonomy corrections performed in the theory--which are then proven to be necessary for achieving this goal. In fact, when approaching the bouncing era, the holonomy corrections become significant and play a crucial role in order to restore the energy conditions. Summing up, a cosmological bouncing scenario with the scale factor above and fulfilling the energy conditions can be adequately described by the F(R) model with a Lagrange multiplier and holonomy corrections, which prove to be very important.


I. INTRODUCTION
There seems to be no doubt that, at present, our universe expansion is accelerating. Its expansion rate is quantified by the Hubble parameter, defined as H =ȧ/a, where a(t) is a scale factor of the universe at cosmic time t. When we go back in time, there are namely two possibilities: (i) the scale factor started from a value zero, what leads to the divergence of the Kretschmann scalar, which in turn ensures the singularity in the spacetime curvature, known as the Big Bang singularity. It may be possible that this singularity is just a manifestation of the shortcomings of classical gravity, unable to describe such small scales (in other words, extremely high energy ones). Quite possibly, a yet-to-be-built quantum theory of gravity will be able to resolve the Big Bang singularity, as turns out to be the case in classical electrodynamics with the Coulomb potential singularities at the origin of the potential, which are resolved in the context of quantum electrodynamics. (ii) In the absence of a fully accepted quantum gravity theory, there is at present another possibility to deal with this issue within the domain of classical gravity. This is known as the bouncing scenario , in which the spacetime curvature singularity is absent. In bouncing cosmology, the universe starts from a contracting era until it reaches a minimal size, then it bounces off at some specific cosmic time and starts to expand. Thereby, the scale factor of the bounce universe does never hit the value zero, which makes the spacetime curvature free from any singularity. Moreover, bounce cosmology is also appealing since it can be obtained as a cosmological solution of the theory of Loop Quantum Cosmology [46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61].
Among the various bouncing models proposed so far, the matter bounce scenario (MBS) [6,14,15,51,[61][62][63][64][65][66][67][68][69][70][71][72][73] has earned special attention because it generates a nearly scale invariant power spectrum. The MBS is characterized by a universe depicted by a matter dominated epoch at very large negative time in the contracting phase, where the primordial curvature perturbations are generated deeply inside the Hubble radius, and is thus able to solve the horizon problem. After it bounces off, the universe enters a regular expanding phase (symmetric to the contracting phase), arXiv:1912.05138v1 [gr-qc] 11 Dec 2019 in which it matches the behavior of the standard Big Bang cosmology. However, in order to obtain a viable matter bounce scenario, the observable parameters of the underlying model have to fulfill a number of stringent constraints, coming from the latest Planck 2018 and other astronomical observations. Along this direction, there are still questions to be answered, within the framework of the matter bounce scenario. Firstly, in an exact MBS, characterized by a single scalar field, the power spectrum turns out to be exactly scale invariant (i.e the spectral index of the curvature perturbation is exactly equal to one), what is in tension with the observational constraints. Such inconsistency of the spectral index in the context of a matter bounce scenario was also confirmed in [70] from a slightly different point of view, namely from an F(R) theory of gravity. F(R) models can be equivalently mapped to scalar-tensor ones via conformal transformation of the metric [75][76][77] and, thus, the inconsistencies of the spectral index from two different models are well justified. Secondly, according to the Planck 2018 data, the running of the spectral index is constrained to be −0.0085 ± 0.0073. However, for the MBS in the case of a single scalar field model, the running of the index becomes zero and hence it is not compatible with observations. Thirdly, in the simplest MBS model the amplitude of scalar fluctuation is found to be comparable to that of tensor perturbations, which in turn makes the value of the tensor-to-scalar ratio to be of order one, again in conflict with the Planck constraints.
However, in a so-called quasi-matter bounce scenario (improving the exact matter bounce one), according to which the scale factor of the Universe evolves as t 3(1+w) (with w = 0), deeply in the contracting era it is possible to recover the consistency of the spectral index and of the running index, even in a single scalar field model. However, the tensor-to-scalar ratio still remains problematic. Fourthly, a crucial drawback of the matter bounce scenario (just as it happens in most of the bouncing models) is the violation of the null energy condition by which the bouncing can be realized.
From a different perspective, it has been shown [74] that an F(R) gravity model with Lagrange multiplier is able to resolve most of the problems arising in the context of matter or quasi-matter bounce scenarios, albeit it fails to restore the energy conditions. Motivated by all these arguments, we will here study matter and quasi-matter bounce scenarios in an F(R) gravity model with a Lagrange multiplier and with the addition of holonomy corrections, which will be proven to be crucial to solve the remaining problems in the description above. The Hubble parameter squared (H 2 ) will be proportional to a linear as well as to a quadratic power of the effective energy density (ρ eff ), unlike in the usual Friedmann cosmology, where H 2 is proportional to the linear power of ρ eff , only. We will explore in this framework the viability of the matter and quasi-matter bounce scenarios, including, in special, the investigation of the energy conditions. The paper is organized as follows. In Sects.
[II], [III], and [IV] we discuss how the holonomy corrections modify the Hamiltonian expression for Einstein's gravity, F(R) gravity, and the Lagrange multiplier F(R) gravity model, respectively. Sects.
[V], [VI], and [VII] are devoted to the explicit calculation of the power spectrum, the observational indices, and the investigation of the energy conditions, in the holonomy corrected Lagrange multiplier F (R) gravity model. Conclusions follow at the end of the paper.
A technical observation is here in order. Before considering the holonomy modifications of the F(R) gravity model with a Lagrange multiplier (an essential ingredient of the present work), it is convenient to discuss the issue of the holonomy improvement in the more standard cases of the Einstein and pure F(R) gravity models. It should be mentioned that holonomy corrections can be introduced in a variety of ways, which are connected to one another by a canonical transformation (see [78], for the details). We will here restrict our analysis to a particular procedure for introducing the holonomy improvement (noted by unbarred quantities in [78]).

A. Ordinary case without holonomy corrections
In the flat Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetime, the Lagrangian of Einstein's gravity along with a scalar field can be written as, where V = a 3 is the volume, R = 2V V − 2V 2 3V 2 is the scalar curvature and Φ the scalar field endowed with the potential U (Φ). It should be noticed that the above Lagrangian contains higher derivatives (a second derivative) term of V and thus, to get the Hamiltonian from this Lagrangian, it is useful to introduce a Lagrange multiplier, namely µ, as follows The equation of motion for R is given by µ = 1 2 V . In order to remove the second derivativeV , we subtract a total derivative term dV dt from the above Lagrangian as (note that the subtraction of a total derivative from a Lagrangian does not change its dynamics)L where we have used µ = V /2. We may note thatL depends on (V, R) and its first derivatives. The corresponding conjugate momenta are With these expressions for the momenta along with Eq. (3), we get the Hamiltonian for Einstein's gravity with a scalar field, asH The above Hamiltonian immediately leads to the Hamiltonian equations, namely (i) the Hamiltonian constraintH = 0 gives Note that the Hamiltonians for Einstein's gravity with and without the bholonomy improvement match with each other in the limit kγ → 0. Using Eq. (5), we get the corresponding Hamiltonian equations as followṡ where we have use Eq.
Eqs. (8) and (9) are the Friedmann equations for holonomy improved Einstein's gravity, while Eq. (10) corresponds to the dynamics of the scalar field with potential U (Φ). It should be noticed that the above Friedmann equations are indeed modified due to the holonomy corrections, as compared to the case without holonomy corrections. Moreover, the above equations of motion have a term containing the generalized momenta p V . However, this term, sin − kγ 2 p V , can be replaced with the matter field energy density, by using Eqs. (6) and (8), which immediately lead to sin 2 − kγ represents the energy density of the scalar field, and ρ c = 3 k 2 γ 2 is known as the critical energy density. With this expression of p V , the gravitational equations can be written as with p = 1 2Φ 2 − U (Φ) is the pressure of the scalar field. It is evident that for kγ → 0 (or equivalently ρ c → ∞), Eqs. (11) and (12) converge to the usual Friedmann equations for usual Einstein's gravity without the holonomy corrections. However it is expected, as for kγ → 0, that the Hamiltonians with and without holonomy corrections match with each other (as mentioned earlier). It is clear, thereby, that the introduction of holonomy corrections modify the expressions of H 2 andḢ by the terms 1 − ρ ρc and 1 − 2ρ ρc , respectively, in comparison with the usual Friedmann equations. Thus, in the holonomy corrected scenario, the squared Hubble parameter is proportional to ρ as well as ρ 2 (apart from some coefficients), which is not the case in the usual FLRW cosmology where H 2 is proportional to the linear power of ρ, only.

III. F(R) GRAVITY WITH A SCALAR FIELD
This section is devoted to the calculation of the Hamiltonian and its corresponding equations for the F(R) model (along with a scalar field) without and with the holonomy improvement (for a general review of F (R) gravity, see e.g., [79][80][81]).

A. Previous to holonomy corrections
The Lagrangian of the F(R) model with a scalar field in the FLRW spacetime geometry is given by where V = a 3 is the volume and R is the scalar curvature. It may be mentioned that a F(R) model can be equivalently mapped to a scalar-tensor theory by using a conformal transformation of the metric, where the scalar field potential of the ST theory depends on the form of F(R). However, here we are interested in obtaining the Hamiltonian in the Jordan frame F(R) model and thus we stick to the Lagrangian shown in Eq. (13) (see [78], for the Hamiltonian formalism of the F(R) model in the Einstein frame). Similar to Einstein's gravity, here we introduce the Lagrange multiplier µ in the above Lagrangian L, to get Variation of R gives µ = 1 2 V F (R), where the prime denotes differentiation with respect to R. Further, to remove the second derivative of V we subtract a total derivative term d dt F (R)V from L 1 , as follows Note that the final LagrangianL depends on the variables (V, R, Φ) and their first derivatives. The corresponding conjugate momenta can be expressed as respectively. These expressions of canonical momenta along with the LagrangianL yield the Hamiltonian for ordinary F(R) gravity (without holonomy corrections) as The Hamiltonian constraintH = 0 leads to the well known Friedmann equation in the F(R) model, namelỹ The other Friedmann equation and the scalar field equation are obtained from the other Hamiltonian equations, as followsṗ respectively. In the following subsection, we will explore the Hamiltonian formalism of F(R) gravity in the presence of holonomy corrections, and examine how such corrections modify the Friedmann equations, in comparison with the usual F(R) case (without holonomy corrections).

B. Holonomy improvement
Performing the holonomy corrections p V → − 2 kγ sin − kγ 2 p V in Eq. (16), the Hamiltonian of F(R) gravity takes the following form, The above expression forH immediately leads to the corresponding Hamiltonian equations, aṡ where b = − kγ 2 p V and, moreover, we use Eqn. (21) The first three equations yield the generalized momenta, while Eqns. (23), (24) are the Friedmann equations of F(R) gravity in the holonomy corrected scenario, and Eqn. (25) is the conservation equation of the scalar field. Note that in the limit kγ → 0, all the above expressions match the equations of the F(R) model without holonomy corrections. Similarly to Einstein's gravity, here we also represent the generalized momenta p V in terms of the scalar field energy density (ρ = 1 2Φ 2 + U ) and the form of F(R). Following [78], we easily get Using this expression, Eqns. (23) and (24) take the form and respectively. It is clear that the holonomy corrections vanish in the limit kγ → 0, as expected.
In a matter or qausi-matter bouncing universe, the holonomy corrections may have significant imprints near the bouncing point, however in the deep contracting era, i.e., at large negative time when the spacetime perturbations are generated, the scalar curvature is large as compared to that of the bouncing era, and thus the holonomy corrections may be safely disregarded in the matter dominated epoch. In most of the previous models, the energy conditions have to be violated, in order to get a non-singular bounce. However in the present paper, we are mainly interested in whether the holonomy corrections may restore the energy conditions in the background of matter or quasi-matter bounce scenario. To this purpose, we consider the Lagrange multiplier F(R) gravity model as it yields a viable phenomenology in a matter (or quasi-matter) bouncing universe, unlike the case of pure F(R) gravity, as explored in our earlier paper [74]. Keeping this in mind, in the next section we shall determine the Hamiltonian and the corresponding equations of the F(R) model with Lagrange multiplier, without and with holonomy corrections, in order to explicitly characterize the modifications generated by the holonomy improvement.

A. No holonomy correction
Let us briefly recall the formalism of F(R) gravity with Lagrange multiplier, developed in Ref. [82]. The action of the model is, where κ 2 = 1 M 2 , with M the four dimensional Planck mass ∼ 10 19 GeV. Here, F (R) and G(R) are two differentiable functions of the Ricci scalar R, Φ is a scalar field with a self coupling kinetic term and the coupling is determined by the function λ, known as the Lagrange multiplier, in the action (28). It was shown in [82], that this variant of F (R) gravity with a Lagrange multiplier term is free of ghosts. However, the Lagrangian contains a higher derivative of V = a 3 and thus the suitable form of the Lagrangian for determining the Hamiltonian is given bỹ with f (R, λ) = F (R) + λG(R), and where we denote f R = df /dR (the same notation for G(R)). It should be observed that the above Lagrangian depends on the variables (V, R, λ, Φ), thus the canonical momenta turns out to be Consequently the Hamiltonian is given by, The Hamiltonian constraintH = 0 and the other Hamiltonian equationṗ V = − ∂H ∂V lead to the well-known Friedmann equations for Lagrange multiplier F(R) gravity without the holonomy corrections [82], as and respectively. The dynamics corresponding to the fields Φ and λ are obtained fromṗ Φ = − ∂H ∂Φ andṗ λ = − ∂H ∂λ , as followsṗ respectively, with E being an integration constant. Taking eqn. (34) into account, the Lagrange multiplier can be written as λ = E a 3 √ G . As a result, the gravitational equations take the form Eqn. (35) denotes the final Friedmann equations for Lagrange multiplier F(R) gravity in absence of holonomy corrections, which also match with [82].

B. Holonomy improvement
Similarly as in previous cases, the holonomy corrections can be introduced by the replacement p V → − 2 kγ sin − kγ 2 p V in eqn. (30) and thus the modified Hamiltonian takes the form Eqn. (36) clearly evidences thatH depends on the variables (V, R, λ, Φ) and their conjugate momenta. The dynamics corresponding to (V, R, λ, Φ) can be obtained aṡ Moreover, the Hamiltonian constraintH = 0 and the other Hamiltonian equationṗ V = − ∂H ∂V yield the Friedmann equations for the Lagrange multiplier F(R) gravity model in presence of holonomy corrections, as respectively, where we have used Eq. (37). The other two equations of motion, To complete the equations of motion, we need to represent sin b in terms of the scalar curvature and the scalar field energy density. Using the above equations, we get (see Appendix-I) This expression of sin b leads to the final form of the gravitational equations, as folllows and respectively. Comparing Eqs. (41) and (42) with Eq. (35), we immediately identify the modification to the Friedmann equations generated by the holonomy improvement. As mentioned earlier, such modifications become significant near the bouncing era in a bouncing universe. We investigate whether the holonomy modifications can restore the energy conditions in the backdrop of a matter or quasi-matter bouncing scenario. However, before the investigation of the energy condition, it is crucial to determine the observable parameters (like the spectral index and tensor to scalar ratio) and directly confront them with the values coming from the most recent Planck observations, which is the subject of next section.

V. REALIZATION OF THE BOUNCING COSMOLOGY
In the present section we will determine the observable quantities for the bouncing universe described by the scale factor where a 0 and n are the free parameters of the model, with a 0 having mass dimension [+2], while n is dimensionless. It must be mentioned that, for n = 1/3, the scale factor describes a matter bounce scenario. The Universe's evolution in a general bouncing cosmology consists of two characteristic eras: an era of contraction and one of expansion. It is obvious that the above scale factor describes a contracting era for the Universe, when t → −∞; then the Universe reaches a bouncing point, at t = 0, at which the Universe has a minimal size, and afterwards the Universe starts to expand again, for cosmic times t > 0. Hence, the Universe in this scenario never develops a crushing type Big Bang singularity.
For the purpose of determining the observable quantities, we will consider the low curvature limit of the theory. So, before proceeding further, let us comment on the viability of this approximation. We do it in the context of matter bounce cosmology, which is obtained by taking n = 1/3 in Eq. (43), the primordial perturbations of the comoving curvature, which originate from quantum vacuum fluctuations, were at subhorizon scales during the contracting era in the low-curvature regime, that is, their wavelength was much smaller than the comoving Hubble radius, which is defined by r h = 1 aH . In the matter bounce evolution, the Hubble horizon radius decreases in size, and this causes the perturbation modes to exit from the horizon eventually, with this exit occurring when the contracting Hubble horizon becomes equal to the wavelength of these primordial modes. However, in the present context, we consider a larger class of bouncing models of the form a(t) = (a 0 t 2 + 1) n , in the presence of a generalized Lagrange multiplier F (R) gravity. Thus, it will be important to check what are the possible values of n which make the low-curvature limit, that is, R/a 0 1 a viable approximation in calculating the power spectrum for the bouncing model a(t) = (a 0 t 2 + 1) n . This expression of the scale factor immediately leads to the comoving Hubble radius Thereby r h diverges at t 0, as expected because the Hubble rate goes to zero at the bouncing point. Furthermore, the asymptotic behavior of r h is given by r h ∼ t 1−2n , thus r h (|t| → ∞) diverges for n < 1/2, otherwise r h goes to zero asymptotically. Hence, for n < 1/2, the comoving Hubble radius decreases initially in the contracting era and then diverges near the bouncing point, unlike in the case n > 1/2, where the Hubble radius increases from the past infinite and gradually diverges at t = 0. As a result, the possible range of n which leads the perturbation modes to exit the horizon at large negative time and make the low-curvature limit a viable approximation in calculating the power spectrum, is given by 0 < n < 1/2. Moreover we will show in the later sections that this range of n makes the observable quantities compatible with the Planck constraints and, thus, the "low curvature limit" comes as a viable approximation to calculate the power spectra of scalar and tensor perturbations. Further, in such low curvature regime (i.e in the deep contracting era), the holonomy corrections may be safely disregarded, which in turn makes the Friedmann equations (see Eqs. (41) and (42), as follows and respectively. Thereby, the Lagrange multiplier F(R) gravity model with/without the holonomy corrections behave similarly in respect to the observable quantities, in the backdrop of a matter or quasi-matter bounce scenario. With these equations, we shall investigate which functional forms of F (R) and G(R) can realize a bouncing Universe cosmological scenario, with the scale factor shown in Eq. (43), which leads to the following Hubble rate and its first derivative With the help of the above expressions, the Ricci scalar is found to be Using Eq. (48), one can determine the cosmic time as a function of the Ricci scalar, that is, the function t = t(R). As a result, the Hubble rate and its first derivative can be expressed in terms of R (however, this statement holds for all analytic functions of t), and also the differential operator d dt can be written as d dt =Ṙ(R) d dR . By plugging the resulting expressions into Eqs. (45) and (46), we obtain differential equations which determine the functional forms of F (R), G(R) fully in terms of R, and by solving those differential equations, the forms of F (R) and G(R) can be found.
During the low-curvature regime ( R a0 1, or for a large negative time), R(t) can be written as R(t) ∼ 12n(4n−1) from Eq. (48). This helps to express the scale factor, the Hubble rate, its first derivative, and the differential operators d/dt, d 2 /dt 2 , in terms of the Ricci scalar R, as follows and respectively. By plugging back these expressions into Eqs. (45) and (46), and by introducing J(R) = F (R)+ 2E √ G(R) a 3 (R) , we get the following differential equations, and Eq. (51) has the following solution where ρ = 1 4 3 − 2n − 1 + 4n(5 + n) and δ = 1 4 3 − 2n + 1 + 4n(5 + n) , and also A and B are integration constants having mass dimensions [2 − 2ρ] and [2 − 2δ], respectively. This solution of J(R) along with Eq. (52), lead to the following functional form of F (R) where C and D are the corresponding coefficients of R ρ and R δ , respectively. With these solutions, the effective f (R) can be written as, where we have used the solution of λ(t) = E a 3 √ G . Eqs. (53), (54), and (55) constitute the main results of the present section. It is interesting to note that due to well-known equivalence of F(R) gravity with generalized fluid models (see [83]), the same scenario may be induced by a specific generalized fluid. In the next section we address concretely the cosmological perturbations issue and we shall confront the theory with the observational data.

VI. COSMOLOGICAL PERTURBATION: OBSERVABLE QUANTITIES
In this section we shall study the first order metric perturbations of the theory at hand, following Refs. [84][85][86], where the scalar and tensor perturbations are calculated for various variants of higher curvature gravity models. Scalar, vector and tensor perturbations are decoupled, as in general relativity, so that we can focus our attention to tensor and scalar perturbations separately.

A. Scalar perturbations
A scalar perturbation of the FRW background can be written as where Ψ(t, x) denotes the perturbation. In principle, perturbations should always be expressed in terms of gauge invariant quantities, namely in our case the comoving curvature perturbation, defined as = Ψ − aHq, where v(t, x) is the velocity perturbation. However, we shall work in the comoving gauge condition, where the velocity perturbation is taken to be zero; with such gauge fixing = Ψ. Thereby, we can work with the perturbed variable Ψ(t, x). The perturbed action up to Ψ 2 order is [84], where z(t) has the following expression It is plain, from Eq. (57), that c 2 s = 1, which guarantees the absence of ghost modes quantified in terms of superluminal propagating modes of the model. Recall that the low curvature approximation stands as a viable check in calculating the observable quantities. Thus, in the low-curvature limit, we determine various terms present in the expression of z(t) (see Eq. (58)), as a(t) Consequently, z(t) takes the following form where P (R) and Q(R) are defined as follows and Before moving further, we check at this stage whether Q(R) goes to zero or, equivalently, z(t) → ∞ at some time value. It is important to examine this, because as we will see, the Mukhanov-Sasaki equation (which is essential to determine the observable quantities) has a term containing 1/z(t) and, moreover, the Mukhanov variable (v = zΨ) diverges at the point where z(t) goes to infinity. As mentioned earlier, perturbations are generated in the low curvature regime deeply in the contracting era and thus the above expression of Q can be simplified, as follows where ρ = 1 4 3 − 2n − 1 + 4n(5 + n) , δ = 1 4 3 − 2n + 1 + 4n(5 + n) and A, B are two integration constants. Further, recall, the explicit expressions of C and D (see Eq. (54)) are Putting these expressions of C and D into Eq. (62) and using the forms of ρ and δ (in terms of n), it can be checked that the above expression for Q is a strictly positive definite quantity for n > 1 4 . Moreover, we will show in the later sections that the observable quantities are compatible with Planck observations for the parametric regime 0.27 n 0.40 (i.e for n > 1/4). Therefore, Q(t) does not hit zero , or equivalently z(t) does not tend to infinity, for the parametric values which are consistent with the Planck observations.
Eq. (57) clearly indicates that Ψ(t, x) is not canonically normalized and, to this end, we introduce the well-known Mukhanov-Sasaki variable, as v = z (= zΨ, since we are working in the comoving gauge). The corresponding Fourier mode of the Mukhanov-Sasaki variable satisfies where τ = dt/a(t) is the conformal time and v k (τ ) the Fourier transform of the variable of v(t, x) for the kth mode. Eq. (64) is quite hard to solve analytically, in general, since the function z depends on the background dynamics. However, the equation can be solved analytically in the regime R/a 0 1, as we now show. The conformal time (τ ) is related to the cosmic time (t) as τ = dt a(t) = 1 a n 0 (1−2n) t 1−2n , for n = 1/2, however we will show that the observable quantities are compatible with Planck data [87] for n < 1/2 and thus we can safely work with the aforementioned expression of τ = τ (t). Using this, we can express the Ricci scalar as a function of the conformal time Having this in mind, along with Eq. (59), we can express z in terms of τ , as follows The above form of z = z(τ ) yields the expression of 1 z d 2 z dτ 2 that is essential for the Mukhanov equation with ξ = (2n+1−ρ) (1−2n) . Recall ρ = 1 4 3 − 2n − 1 + 4n(5 + n) and δ = 1 4 3 − 2n + 1 + 4n(5 + n) , which clearly indicate that δ − ρ is a positive quantity. Thus, the term within parenthesis in Eq. (67) can be safely considered to be small in the low-curvature regime R/a 0 1. As a result, 1 z d 2 z dτ 2 becomes proportional to 1/τ 2 i.e., 1 z d 2 z which is approximately a constant in this era, when the primordial perturbation modes are generated deeply inside the Hubble radius. In fact, and in conjunction with the fact that c 2 s = 1, the Mukhanov equation can be solved, as follows with ν = σ + 1 4 , and c 1 and c 2 are integration constants. Assuming the Bunch-Davies vacuum initially, these integration constants become c 1 = 0 and c 2 = 1, respectively. Having the solution of v k (τ ) at hand, next we proceed to evaluate the power spectrum (defined for the Bunch-Davies vacuum state) corresponding to the k-th scalar perturbation mode, which is defined as In the superhorizon limit, using the mode solution in Eq. (69), we get By using Eq. (71), we can determine the observable quantities, as the spectral index of the primordial curvature perturbations and the running of the spectral index. Before proceeding to calculate these observable quantities, we will consider first the tensor power spectrum, which is necessary for evaluating the tensor-to-scalar ratio.

B. Tensor perturbations
Let us now focus on the tensor perturbations, which for the FRW metric background are noted by h ij (t, x) and defined as A tensor perturbation is itself a gauge invariant quantity, and the tensor perturbed action up to quadratic order is given by where z T (t) is Therefore, the speed of the tensor perturbation propagation is c 2 T = 1. Similar to scalar perturbations, the Mukhanov-Sasaki variable for tensor perturbation is defined as (v T ) ij = z T h ij , which, upon performing the Fourier transformation, satisfies the equation Using Eq. (74), along with the condition R/a 0 1, we evaluate z T (τ ) and 1 z T (τ ) d 2 z T dτ 2 , and these read and respectively, where S(R(τ )) = ρ(A+C) and also we used R = R(τ ) from Eq. (65). Due to the fact that δ − ρ is positive, the variation of the term in parenthesis in Eq. (77) can be regarded to be small in the low-curvature regime, and thus 1 and recall that ξ = (2n+1−ρ) (1−2n) . The above expressions yield the tensor power spectrum, defined with the initial state being the Bunch-Davies vacuum, so we have The factor 2 arises due to the two polarization modes of the gravity wave, and ν T = σ T + 1 4 , where σ T is defined in Eq. (78). Now, we can explicitly confront the model at hand with the latest Planck observational data [87], so we shall calculate the spectral index of the primordial curvature perturbations n s and the tensor-to-scalar ratio r, which are defined as follows Eqs. (71) and (79) immediately lead to the explicitly form of n s and r, as where σ, z(τ ) and z T (τ ) are given in Eqs. (68), (65), and (76), respectively. As it is evident from the above equations, n s and r are evaluated at the time of exit from the horizon, when k = aH or, equivalently, at τ = τ h . It may be noticed that n s and r depend on the dimensionless parameters R h a0 and n with R h = R(τ h ). We can now directly confront the spectral index and the tensor-to-scalar ratio with the Planck 2018 constraints [87], which constrain the observational indices as follows n s = 0.9649 ± 0.0042 , r < 0.064 .
For the model at hand, n s and r are within the Planck constraints for the following ranges of parameter values: 0.01 ≤ R h a0 ≤ 0.07 and 0.27 n 0.40. This behavior is depicted in Fig. 1. The viable range of R h /a 0 is in agreement with the low-curvature condition R/a 0 1 that we have considered in our calculations. Moreover, the range of the parameter n clearly indicates that the matter bounce scenario, for which n = 1/3, is well described by the Lagrange multiplier F (R) gravity model with the holonomy corrections (though the holonomy corrections may be disregarded in calculating the observable quantities of the low curvature regime, as mentioned earlier). At this stage it is worth mentioning that in scalar-tensor theory (with an exponential scalar potential), the matter bounce scenario is not consistent with the Planck observations. Moreover, the matter bounce scenario also does not fit well even in the standard F (R) gravity, as we confirmed in an earlier paper [74]. However, here we show that for the Lagrange multiplier generalized F (R) gravity, the matter bounce may indeed be considered as a good bouncing model, which allows the simultaneous compatibility of n s and r with observations. Here it may be mentioned that an analogue study of matter bounce cosmology can be found in [62] within another type of modified gravity, in particular F (T ) gravity, where the power spectrum becomes nearly scale invariant but suffers from an over large tensor-to-scalar ratio. However, in the present paper, we showed that within holonomy corrected Lagrange multiplier F (R) gravity, the matter bounce scenario gives rise to a nearly scale invariant power spectrum and also the tensor-to-scalar ratio lies within observational bound. The running of the spectral index is defined as and this is constrained by the Planck 2018 results, as α = −0.0085 ± 0.0073. Thus, it is also important to calculate the running of the spectral index before concluding the viability of a model. By using the expression of σ (see Eq. (68)) and R = R(τ ) (see Eq. (65)), we get To arrive at the above result, we use the horizon crossing relation of the k-th mode k = aH to determine d|τ | d ln k = −|τ | i.e., the horizon exit time |τ | increases as the momentum of the perturbation mode decreases, as expected. Eq. (84) indicates that, similarly to n s and r, the running index (α) also depends on the parameters R h /a 0 and n. Taking R h /a 0 = 0.05, we have produced a plot of α with respect to n, in Fig. 2. As can be seen in Fig. 2 takes negative values, crossing zero about n 0.30. Thus α lies within the Planck constraint for 0.30 n 0.40, which includes the matter bounce scenario. For the Lagrange multiplier generalized F (R) gravity model, we showed that the pure matter bounce scenario, as well as the quasi matter bounce scenario, are both consistent with the Planck observations. Therefore, generalized F (R) gravity in terms of the Lagrange multiplier has a richer phenomenology, in comparison with scalar-tensor or standard F (R) gravity models, which fail to describe in a viable way these two bouncing cosmology scenarios.

VII. ENERGY CONDITIONS
A crucial drawback in most bouncing models is the violation of the null energy condition near the bounce at which the bouncing phenomena is realized. However in the present context, we deal with the matter or quasi-matter bounce in the backdrop of the Lagrange multiplier F (R) gravity model in presence of holonomy corrections. And it does happen that such holonomy corrections may become significant about the bouncing point (where the spacetime curvature is large compared to the one at present) and thus may play an important role in restoring the energy conditions. Keeping this in mind, here we check the energy conditions in Lagrange multiplier F(R) gravity with holonomy corrections, where the effective energy density ρ eff and pressure p eff can be determined from anḋ respectively, with f R = F R + EG R a 3 √ G . Eq. (85) can be simplified to ρ 2 eff − ρ eff ρ c + 3H 2 ρ c = 0, which can be solved . At the bouncing point, the Hubble parameter becomes zero, H(t = 0) = 0, which immediately leads to ρ eff (t = 0) = 0 (we call + and − the branch solutions of ρ eff ). Recall that, without the holonomy improvement, the effective energy density goes to zero at the bouncing point, whereby the + branch solution comes just due to the presence of holonomy corrections. Considering that such holonomy corrections have an effect on the evolution of ρ eff near the bounce, we take the + branch solution; otherwise the holonomy corrections would have no effect on the effective energy density even at the bouncing point. Thus the evolution of ρ eff is given by Putting this solution into Eq. (86), we get and using the expressions of the Hubble parameter in the present context (H = dota a with a(t) = a 0 t 2 + 1 n ), we give the plots of ρ eff /ρ c and ρ eff + p eff (with respect to the cosmic time) as the left and right plots of Fig. 3, respectively, for n = 1/3, kγ = 1, a 0 = 1 (in reduced Planck units). As it can be observed from Fig. 3 that ρ eff + p eff asymptotically goes to zero, as expected sinceḢ vanishes for t → ±∞. Moreover Fig. 3 clearly demonstrates that both the weak and the null energy conditions are satisfied in the holonomy corrected version of Lagrange multiplier F (R) gravity, unlike in the usual Friedmann case where the energy conditions are generally violated near the bouncing point. It is worth mentioning that in the usual Friedmann cosmology, the gravitational equations come as H 2 = ρ eff 3 andḢ = − 1 2 ρ eff + p eff respectively and thus at the bouncing point (where H = 0 andḢ > 0), ρ eff + p eff becomes less than zero, which implies the violation of the energy conditions. However in presence of the holonomy improvements, there are extra terms in the equation of motion, as for exampleḢ = − 1 2 ρ eff + p eff 1 − 2ρ eff /ρ c , and thus, at the bouncing point (where ρ eff = ρ c andḢ > 0), ρ eff + p eff becomes positive, so that the energy conditions are restored. Therefore, it is clear that the holonomy improvement plays a significant role to rescue the energy conditions for the matter (or quasi-matter) bounce scenario.
Before concluding, we determine the form of the effective f(R) in the whole curvature regime. This has to be done numerically, owing to the complicated equations of motion. In a previous Section [], we noticed that in the lowcurvature regime, f (R) goes as f (R) ∝ R ρ . Recall, ρ = 1 4 3 − 2n − 1 + 4n(5 + n) which is negative for n > 0.25, and as shown earlier, the present model is consistent with the Planck results for 0.27 n 0.4, hence ρ is negative, in order to ensure the viability of the model. Therefore, it is clear that in the low-curvature regime f (R) is proportional to an inverse power of Ricci scalar ∝ R −|ρ| , using such form of f (R) as boundary condition along with the expression R(t) = 12n (4n−1)t 2 +1/a0 (t 2 +1/a0) 2 , we solve Eqs. (41) and (42) numerically, with the cosmic time t being the independent variable. Moreover, a 0 , n and kγ are taken as a 0 = 1 (in reduced Planck units), n = 1/3 and kγ = 1, respectively; so in effect we consider the matter bounce scenario. However, it may be mentioned that the n = 1/3 case makes the model consistent with the Planck 2018 constraints, as confirmed in the previous section. The numerical solution of f (R) in terms of R is obtained by using the expression R(t) = 12n (4n−1)t 2 +1/a0 (t 2 +1/a0) 2 and is presented in Fig. 4. We also give a plot of the Ricci scalar in the same diagram, in order to make a comparison of the effective f(R) in the present context with the one corresponding to Einstein's gravity. Having such form of f(R), it is important to explore whether this form of f (R) passes the astrophysical tests in the low curvature regime and, on the other hand, gives rise to an accelerating universe at late time. An example of the tests of infrared instability is matter instabilities, which are related to the fact that the spherical body solution in general relativity may not be a solution in the modified gravity theory considered. Matter instabilities may appear when the energy density or the curvature is large compared with the average density or curvature in the universe, as is the case inside of a planet. Following [80], we immediately write the potential (U (R b ), with R b the perturbed Ricci scalar) for the perturbed Ricci curvature over Einstein gravity, as where we consider f (R) ∼ R ρ in low curvature regime and we denote d k f (R)/dR k = f (k) (R). If U (R b ) > 0, the perturbation grows with time and the system becomes unstable. Recall ρ < 0 and δ > 0 and thus the term R ρ dominates over R δ in the low curvature regime. Thus we can approximate F (R) ∼ R ρ in the low curvature regime which immediately leads to the potential, as In the low curvature regime, the first term dominates in the above expression of U (R b ), and thus U (R b ) becomes negative. This indicates that the model considered here passes the matter instability test. Moreover, here it may be mentioned that for |ρ| = 1, the above potential U (R b ) becomes positive, which indicates that f (R) ∼ 1/R (in low curvature regime) leads to an infrared instability, as also confirmed in [80]. However for the purpose of investigating the late time acceleration in our present model, one may check that f (R) ∼ R ρ leads to a power solution of the scale factor as a(t) = a 0 t − (ρ−1)(2ρ−1) (see [80]), recall ρ = 1 4 3 − 4n − 1 + 4n(5 + n) . We showed in the earlier section that the viable range of n for which the observable parameters (n s , r and α) are simultaneously compatible with Planck observations is given by 0.30 n 0.40. However this range of n leads to the exponent in the solution of a(t) as 0 < − (ρ−1)(2ρ−1) (ρ−2) < 1, which in turn indicates a decelerating expanding universe at late time. The deceleration of the late time scale factor is also confirmed from the evolution of Hubble radius (the Hubble radius is increasing at the late time of the expanding phase), as explained in Section [V]. Thereby as per the evidences from supernova data [88,89], which points out that the present universe undergoes through an accelerating phase, our present model is unable to describe the behaviour of our universe at late time. However our results indicate that the present model must be combined with another cosmological scenario in the low curvature regime of the expanding phase, or should be modified appropriately, in order that it leads to a viable cosmology at late time. Such unification of bouncing with late time acceleration is an interesting issue and is expected to study in a near future work.

VIII. CONCLUSIONS
We have discussed a bouncing scenario which incorporates holonomy corrections in an F(R) gravity model with a Lagrange multiplier and where the Hubble parameter squared is proportional (aside from some coefficients) to a linear as well as to a quadratic power of the effective energy density; this differs from the usual Friedmann case, where H 2 is proportional to the linear power of energy density, only. We have specified our study to matter or quasi-matter bounce scenarios, with the scale factor being a(t) = a 0 t 2 + 1 n , where n is a parameter of the model and t the cosmic time. For such bouncing scale factor, it was shown that, for n < 1/2, the primordial curvature perturbation modes exit the Hubble horizon (defined by r h = 1/aH) at a large negative timevalue deeply in the contracting era (where the spacetime curvature is low as compared to the bouncing point), which in turn makes the "low curvature approximation" a viable one to calculate the power spectra of the scalar and tensor perturbations. We have determined the form of the effective f(R) theory in the low curvature regime, by using a reconstruction technique for the present model, which realizes the bouncing with the aforementioned scale factor. The form of f(R) leads to the explicit expression of the well known Mukhanov-Sasaki equation, by solving which we have determined the power spectrums of primordial perturbations and, correspondingly, have calculated the fundamental cosmological parameters, such as the spectral index of scalar perturbation, the tensor to scalar ratio, and the running spectral index. Such observational indexes are found to depend on the dimensionless parameters of the model, as R h /a 0 (with R h the Ricci curvature at the time of horizon exit) and n. It turned out that the observable quantities are simultaneously compatible with the Planck2018 constraints for the parameteric range 0.01 R h a0 0.07 and 0.27 n 0.40. It should be noticed that this range of n is also supported by the range 0 < n < 1 2 , which makes the low-curvature approximation a perfectly reliable one for calculating the power spectra.
Further, we have determined the explicit expressions for the effective energy density (ρ eff ) and effective pressure (p eff ) in the present holonomy improved scenario, from which we have investigated the problem of the energy conditions. Concerning this issue, we have obtained a quite remarkable result, namely that both the weak and the null energy conditions are fulfilled, owing precisely to the presence of the holonomy modifications in our Lagrange multiplier F(R) gravity model. This is a significant difference with respect to the usual Friedmann case, where the energy conditions are generically violated near the bouncing point. Actually in usual Friedmann cosmology, the gravitational equations have the form H 2 = ρ eff 3 andḢ = − 1 2 ρ eff + p eff , respectively, and thus at the bouncing point (where H = 0 anḋ H > 0), ρ eff + p eff becomes negative, what leads to the violation of the energy conditions. However, because of the holonomy corrections, extra terms appear in the equation of motion, as for instanceḢ = − 1 2 ρ eff + p eff 1 − 2ρ eff /ρ c , and thus at the bouncing point (where ρ eff = ρ c andḢ > 0), ρ eff + p eff becomes positive, what shows that the energy conditions are restored. Thus, the holonomy corrected Lagrange multiplier F(R) gravity model discussed here is fully vindicated as a viable description for the bouncing scenario with the aforementioned scale factor.

Appendix-1
Here we give the details of the calculation of Eq. (40). The demonstration goes as follows.

Introducep as
With this new defined quantity, the first and second expressions of eqn.(37) takes the following form