Bouncing cosmology in the limiting curvature theory of gravity

In this paper we discuss models satisfying the limiting curvature condition. For this purpose we modify the Einstein-Hilbert action by adding a term which restricts the growth of curvature. We analyze cosmological solutions in such models. Namely, we consider a closed contracting homogeneous isotropic universe filled with thermal radiation. We demonstrate that for properly chosen curvature constraints such a universe has a bounce. As a result its evolution is nonsingular and contains a"de Sitter-type"supercritical stage connecting contracting and expanding phases. Possible generalizations of these results are briefly discussed.


I. INTRODUCTION
The idea that the Universe can have a prehistory before the big bang is very old. Cyclic or oscillating cosmological models were considered almost 90 years ago. Such models were discussed in the famous book by Tolman published in 1934 [1]. He also demonstrated that the validity of the second law of thermodynamics applied to the Universe and increasing entropy make pure periodic models impossible: each of the successive cycles should be longer and larger than the previous one.
Even if one does not require the existence of an infinite number of cycles before the formation of the present Universe it is interesting to analyze an option that the Universe before the big bang had a phase of contraction, usually called a big crunch. In such models the Universe should experience a bounce, where its size takes some minimal value. We denote the scale factor that enters the Friedmann-Robertson-Walker metric for a homogeneous isotropic universe as a(t). Then, one of the Einstein equations implies thaẗ where ε and P are the matter energy density and pressure, respectively. Since at the bounce point whereȧ = 0 one hasä > 0 the equation of state should be such that ε + 3P < 0. The famous Penrose-Hawking singularity theorems imply that in a general case in order to escape a cosmological singularity some of the energy conditions for matter should be violated [2]. Singularities in standard cosmological models are connected with an infinite growth of the spacetime curvature. Markov [3,4] suggested that the existence of the limiting curvature should be considered as a new physical principle. He demonstrated that for a proper choice of the equation of state in the cosmology the limiting curvature condition is satisfied and solutions in such a model describe a bouncing universe. A bouncing universe was discussed by Gasperini and Veneziano in their * vfrolov@ualberta.ca † zelnikov@ualberta.ca pre-big-bang string cosmology [5,6]. Nonsingular cosmological models that are based on the use of nondynamical scalar fields to implement the limiting curvature hypothesis were studied some time ago by Brandenberger, et al. [7,8]. More recently, the interest in bouncing cosmological models has increased. This is mainly connected with the remarkable increase in the accuracy of cosmological observations. An interesting and intriguing question is: if there was of a big crunch phase, is it possible to find observational evidence of this? A variety of different proposed bouncing cosmological models have been discussed in several nice review articles, which also contain references to the original publications [9][10][11][12][13][14][15][16][17][18][19][20].
In this paper we discuss bouncing cosmological models in a new recently proposed limiting curvature gravity (LCG) theory [21]. The main idea of this approach is to modify the Einstein-Hilbert action by adding a constraint term which controls the curvature behavior and forbids its infinite growth. In fact, this is a realization of the Markov's old idea about the existence of a limiting curvature. A limiting curvature modification of a twodimensional dilaton gravity was considered in [21]. In this paper we discuss four-dimensional LCG models. In the Friedmann-Robertson-Walker metric for a homogeneous isotropic universe the Weyl tensor vanishes. Therefore, it is sufficient to restrict the growth of the Ricci tensor.
We shall discuss two types of models. We first introduce linear-in-curvature constraints. For this purpose we add to the action terms that are linear in the Ricci scalar and the eigenvalues of the Ricci tensor. After this, we discuss quadratic-in-curvature constraints. In both cases, we demonstrate that there exists a wide class of curvature constraints for which the curvature remains uniformly bounded during the evolution of the universe. A common property of such limiting curvature gravity models is that the cosmological solutions have a bounce. A contracting universe at some stage of its evolution, when its curvature reaches the critical value, enters a supercritical regime. If the initial size of the universe was large, then the corresponding supercritical solution is always close to the de Sitter solution. After passing the bounce point the universe expands. We demonstrate that at some moment of time it can leave its supercritical regime and one arXiv:2108.09927v2 [hep-th] 27 Nov 2021 gets an expanding universe filled with matter. After this it follows the standard Einstein equations.
The paper is organized as follows. In Sec. II we recall some well-known properties of the isotropic homogeneous cosmological models and introduce notations that are used later in the paper. LCG models and reduced actions for these models are discussed in Sec. III. Sections IV-VII discuss LCG models with linear-in-curvature constraints. Sections VIII-X are devoted to study LCG models with quadratic-in-curvature constraints. More general curvature constraints are discussed in Sec. XI. Finally, Sec XII contains a summary of the obtained results, a discussion of different aspects of LCG cosmological models and their possible generalizations. Some additional technical details and results used in the main text are collected in the Appendix.

II. ISOTROPIC HOMOGENEOUS COSMOLOGY
Let us consider the cosmological metric in the form (2.1) This metric is a direct sum of the one-dimensional metric b 2 (t)dt 2 and three-dimensional metric a 2 (t)dγ 2 , where dγ 2 = γ ij dx i dx j is a line element on a unit 3D sphere S 3 . The metric dγ 2 admits group O(4) of symmetries. It is well known that: • A scalar function on S 3 invariant under the action of this group is a constant.
• There does not exist a nonvanishing vector field invariant under the group of symmetries.
• A symmetric rank-two tensor field A ij on S 3 invariant under the group of symmetries is A ij = Aγ ij , where A is a constant.
Consider a symmetric tensor A in a spacetime with metric (2.1) which respects its symmetry. Then, it has the following form: It is easy to see that A(t) andÂ(t) are eigenvalues of the tensor A µ ν . We call them temporal and spatial eigenvalues, respectively.
In what follows we use similar notations for other symmetric rank-two tensors. For example, the Ricci tensor R µν has the form Then, the Ricci scalar is We keep the coefficient b 2 (t) of the metric (2.1) as an arbitrary function. This will allow us to obtain a complete set of the gravitational field equations from a reduced metric, but later, after the variations, we put b(t) = 1. This is nothing but a gauge-fixing condition corresponding to synchronous gauge. The eigenvalues R andR of the Ricci tensor can be expressed in terms of two structures p and q R = 6(q + p), R = 3q,R = q + 2p, (2.5) (2.6) The traceless part of the Ricci tensor has the following eigenvalues The Weyl tensor for the metric (2.1) vanishes, that is, all of the information about the spacetime curvature is encoded in the Ricci tensor. Our goal is to study cosmological models that obey the limiting curvature condition. A natural way to do this is to impose restrictions on the eigenvalues of the Ricci tensor. For example, one may try to restrict the value of the Ricci scalar (2.4) which is a linear combination of these eigenvalues. However, the form (2.5) of this invariant implies that this does not work. The reason is simple: the function p(t) is positive definite, while the function q(t) does not have a definite sign. Hence, the growth of p(t) for a contracting universe can be compensated by an increasing negative value of q, so that R remains bounded. The well-known example of a contracting universe filled with a thermal radiation clearly illustrates this. The Kretschmann invariant for this solution grows infinitely, so that the limiting curvature condition is violated. In what follows we discuss constraints that can be used to prevent infinite curvature growth.

A. Action and gravity equations
The limiting curvature gravity model is described by an action of the form where I g is the Einstein-Hilbert action,

2)
and I m is the matter action. The term I c is the constraint action depending on the metric and the Lagrange multipliers that generate constraints on the curvature. Later, we will specify the form of these constraints and the corresponding term of the action I c . Now we just mention that the imposed restriction on the curvature has the form of inequalities (see discussion in Ref. [21]). They have the following properties. Before the curvature reaches its critical value the corresponding subcritical metric coincides with a standard solution of the unmodified Einstein equations. After the curvature reaches the critical value the solution becomes supercritical and it follows the modified constraint equations which prevent further growth of the curvature. The variation of the total action (3.1) over the metric gives the following "gravity" equations: which have the form Here G µν is the Einstein tensor and are the stress-energy tensors of matter and constraints. Besides the gravity equations the action I also gives additional equations for matter and constraints which are obtained by its variation over the Lagrange multipliers, that are variables additional to the metric. If these equations are satisfied and the actions I c and I m are covariant, the following relations are valid: These conservation laws guarantee consistency of the gravitational equations.

B. Reduced action and reduced gravity equations
The tensor G for the metric (2.1) is Similarly, the tensors T and X, respecting the symmetry of the metric (2.1) have the form X µ ν = diag(X ,X ,X ,X ) . (3.10) Then the gravity equations (3.4) reduce to the following equations Equation (3.6) and the conservation property of the Einstein tensor, G µν ;ν = 0, givê 13) 14) In particular, these relations imply that if the temporal gravity equation (3.11) is valid, the spatial gravity equation (3.12) is also satisfied. It is convenient to use symmetries of the cosmological spacetimes and write down a reduced action for our gravitational system. Namely, it is easy to check that the 4D gravity equations (3.3) taken on the spacetime (2.1) can be equivalently derived from the dimensionally reduced action where L is the Lagrangian of the system evaluated on the metric (2.1) and 2π 2 is the volume of the unit sphere S 3 . For example the dimensionally reduced Einstein action is The variation of this action over the temporal b and spatial a components of the metric after imposing the gaugefixing condition b = 1 gives (3.18)

C. Thermal radiation
We choose the stress-energy tensor of the matter in the form where ε is the energy density and P is the pressure. In what follows we assume that the matter is hot thermal radiation with the equation of state P = 1 3 ε. The conservation law (3.15) is satisfied if where the factor C is defined by the temperature T of radiation and the number of massless degrees of freedom n. 1 For a closed homogeneous and isotopic universe filled with a thermal radiation the total energy E and entropy S of the universe are given by Here T is the radiation temperature and V = 2π 2 a 3 is the volume of the closed universe. In the case of only electromagnetic radiation n = 2. At high temperature many other fields become effectively massless. For example at the temperature corresponding to 300 GeV this number is about n ≈ 106. The relation (3.21) allows one to express the constant C in Eq. (3.19) in terms of the entropy S, which is a conserved quantity. One gets and, hence, For pure electromagnetic radiation n = 2 and we have ν ≈ 0.014686. Thus the value of the constant C is defined by the entropy of the thermal gas in the Universe and can be estimated from observations [22]. For example the contribution of photons to the entropy is S/k B ∼ 5.4 · 10 89 and, hence, C ∼ 6.5 · 10 117 c, which is a huge number. The other massless particles like neutrinos contribute similar amounts to the entropy and energy density. Let us note that at the stage of contraction the thermal radiation dominates. When the growing temperature becomes high enough particles with mass m T becomes ultrarelativistic and their contribution to the energy density is similar to the contribution of massless particles (photons), while the contribution of the particles with m > T is relatively small. For this reason, in what follows we assume that the contracting universe is radiation dominated.

A. General form of linear constraints
In order to control curvature growth one can impose a restriction on the eigenvalues of tensors constructed as a 1 Let us note that the stress-energy tensor for thermal radiation can be derived from the reduced action 20) linear combination of the Ricci tensor R µ ν and Rδ µ ν , where R is a Ricci scalar. We call such constraints linear in curvature. Let us discuss the case of the linear constraints first and return to the discussion of other constraints constructed from curvature invariants later.
As earlier, we denote by S a traceless part of the Ricci tensor, and Q = 1 6 RI. Here I is a unit tensor. We denote Z = c S S + c R Q .
We shall restrict the curvature by imposing the conditions Consider two-dimensional (p, q) plane (see Fig. 1). A domain in this plane where the inequalities (4.5) are satisfied is restricted by straight lines, We call of these lines as describing temporal and spatial constraints, respectively.

B. Constraint action
In order to provide the inequality constraints (4.5) we add the following expressions to the reduced action (3.16): (4.7) As we shall see later, the constraint Z = −Λ does not give any restrictions on the physically interesting solutions. That is why we did not include the corresponding term in the action.
The variation of this action over the Lagrange multipliers gives the following equations These equations imply that the system has two different regimes. In the subcritical regime where χ =χ ± = 0, the nonvanishing parameters ζ andζ ± are defined in terms of Z andẐ ± , respectively. In this regime the action I c does not contribute to the gravity equations, so the evolution of the universe follows its standard solutions of the unmodified Einstein equations.
In the supercritical regime, when one of the constraint equations is saturated and the corresponding Lagrange multiplier ζ orζ ± becomes zero. This means that one of the constraint equations Z − Λ = 0 , (4.10) is valid. The corresponding control function χ orχ ± becomes "dynamical" and its evolution in the supercritical regime is defined by the gravity equations.
The contributions X andX of the constraint action I c to the gravity equations can be obtained as follows. Since the constraint functions Z andẐ ± are linear combinations of the functions p and q given by Eq.(2.6), it is sufficient to find the variations over the metric of the following reduced actions: I q = 2π 2 dt a 3 b uq , (4.13) (4.14) Here u = u(t) stands for one of the control functions. The variations of these reduced actions over the temporal b and spatial a coefficients of the metric after imposing the gauge-fixing condition b = 1 give (4.15)

V. SUBCRITICAL SOLUTIONS
At this stage the control functions χ orχ ± vanish and the functions ζ andζ ± drop out of the equations. As a result the standard Einstein equations govern the dynamics of the radiation-dominated universe. The temporal Einstein equation is In an explicit form this equation readṡ Here we fixed the gauge and put b = 1. In this gauge we have The spatial Einstein equation is the consequence of Eq.(5.2) and it reduces tö Note that the scalar curvature during this stage of evolution vanishes. This means that a point representing the state of the universe in the (p, q) plane moves along the line Γ − where q = −p (see Fig. 1). A solution of the equation Eq.(5.2) is well known (see, e.g., Ref. [23]). It has the form where the integration constant is fixed by the condition a(t = 0) = 0. Here a max is the maximal value of the scale factor of the Friedmann universe. Note that t < 0 during the contraction stage. During the collapse of the universe its scale parameter a(t) decreases and the energy density of matter grows. This subcritical regime continues until the curvature reaches the critical value at which the supercritical regime starts. It happens when one of the constraint equations (4.10)-(4.11) is satisfied. We denote the corresponding critical value of p by λ − . Thus, the supercritical solution starts at the point (p = λ − , q = −λ − ). Let us find the parameters a 0 andȧ 0 of the contracting universe at this point. For p = λ − Eqs.(5.2) and (5.3) give Using the definition of p one getṡ In the latter relation we choose a minus sign since we assume that the universe is initially contracting. The moment of transition to the supercritical stage is We are looking for constraints that restrict the curvature, so that during all of the subsequent evolution of the universe after it enters a supercritical regime the curvature remains finite and restricted by a chosen universal value. To characterize the value of the curvature one can use, for example, the Kretschmann invariant K Eq.(2.10).
Let a(t) be a solution for a scale function which determines the size of the universe. For a general constraint after the solution enters the supercritical regime it may terminate at some finite time t s . This may happen if the differential equation for a(t) determined by the constraint has a singular point which prevents an extension of the solution beyond time t s . We call such a constraint, which does not allow a complete description of the evolution of the universe, a singular one. In what follows we shall not consider such constraints. Namely, we assume the following.
• The supercritical solution is not terminated at finite time t.
• The constraint guarantees that during the supercritical regime the Kretschmann invariant is uniformly restricted by some universal value which does not depend on the parameters of the solution.
We call such a constraint a regular one. For this type of constraint, the corresponding supercritical solution can either be continued to t → ∞ or slip back to its subcritical phase. In principal, if there exist several constraints, the supercritical solution can also slip between them. We assume that the constraint line C intersects Γ − at p = λ − where the solution enters the supercritical regime. For the contracting universeȧ < 0 at this point. Using the definition (5.3) for p and q, one can obtain the following equation: While a point on C representing a supercritical solution is located in the domain where q < 0 the negative value ofȧ can only increase. Thus in this domainȧ < 0. One also has q − p < 0 in the domain of the (p, q) plane below Γ − . Equation (6.1) shows that under these conditions the point representing the supercritical contracting universe in the (p, q) plane can move only with an increase of the parameter p.

B. Regular linear constraints
To describe the evolution of the universe we use the two-dimensional (p, q) plane. Let Γ ± be two lines on this plane defined by the equations q = ±p, respectively. The subcritical evolution of a contracting radiationdominated universe is represented by the interval on Γ − (see Fig. 2).
Let C be a straight line representing the linear-incurvature constraint. We assume that this line intersects Γ − at a point 0 and write its equation in the form The parameter Λ has dimensions of [length] −2 and characterizes the value of the limiting curvature. Since the left-hand side of Eq.(6.2) is positive at the point 0, Λ is chosen to be positive as well. In the presence of the constraint (6.2) a point representing the evolution of the universe after it reaches the point 0 starts its motion along a constraint C. Let us discuss the corresponding supercritical solution.

Negative-µ case
Let us first assume that the parameter µ in Eq.(6.2) is negative. Then, dp/dq = µ < 0. If C does not meet another constraint, both p and |q| along C grow infinitely and as a result the Kretschmann invariant grows as well. This means that such a constraint is not regular. For this reason we assume that µ > 0.

µ > 1 case
Let us consider the case where µ > 1. At the point 0 where the supercritical regime starts p − µq > 0 and one After this, a representative point which is moving with the increase of p enters the domain above Γ − and remains there since the corresponding line C cannot intersect Γ + . To find how the scale factor behaves in this case we rewrite Eq.(6.1) in the form dp d ln(a 2 /a 2 0 ) Here a 0 is the size of the universe at the beginning of the supercritical regime when p = λ − . Integrating this equation with the imposed initial conditions we get .
Here F is the expansion factor and the function q(p) is defined by the Eq.(6.2). The integral can be easily calculated and one has The integration constant is chosen so that F | p=λ− = 0. Thus the relation between a and p takes the form . (6.7) Since µ > 1 and p grows monotonically the scale function a monotonically decreases. The Kretschmann invariant grows infinitely along the constraint while the size of the universe shrinks. Thus such a constraint is not regular.
Let us consider the last case where 0 < µ < 1. In this case the constraint line C crosses Γ + . At the point of the One also has Let us introduce the dimensionless quantities Then one has Equation (6.5) can be used to find a relation betweenp and α. It is sufficient to substitute into the expression for F the function q(p) defined by Eq.(6.11). The integral can be easily calculated and one has Thus the relation between α andp takes the form . (6.13) By inverting this relation we findp as a function of α, and then by using Eq.(6.11) we also computeq as a function of α,p Let us demonstrate that for a chosen linear constraint a contracting universe always has a bounce. Let us assume that such a bounce exists. At this point, whereȧ = 0, the universe has a minimal size, which we denote by α b . Letp b be the corresponding value of the parameterp at this point. Using Eq.(5.3), one has After using Eq.(6.13), this condition takes the form For every 0 < µ < 1 the function on the left-hand side of this relation grows infinitely whenp b → 1. This means that for an arbitrarily large α 0 Eq.(6.17) has a solution. In other words the universe has a bounce. For large α 0 , this happens whenp b is close to 1. In this case one can omit the termp in Eq.(6.16). Equation (6.17) allows one to expressp as a function of α. After substituting this expression into the relation one obtains the equation that determines the evolution of the universe in the supercritical regime. Here a prime is a derivative with respect to τ . After the size of the universe reaches the minimal value a b it expands again. The point representing it in the (p, q) plane moves again along the line C but now in the opposite direction with the decreasing value ofp. At the point where the solution intersects Γ − it can leave the supercritical phase. Such a solution describes an expanding universe filled with thermal radiation. Let us emphasize that during its evolution in the supercritical regime the value of the Kretschmann invariant remains uniformly restricted. Thus the linear constraint (6.2) with 0 < µ < 1 is regular.

C. Temporal and spatial constraints
Both temporal and spatial curvature constraints can be written in a form similar to Eq.(6.2). Let us first apply the results of the previous section to the temporal constraint. In order to present it in the form (6.2) it is sufficient to choose the coefficients c R and c S in Eq. (4.3) in the form Then one has We choose 0 < µ < 1. Then the temporal constraint . The evolution of the universe is represented in the (p, q) plane by two intervals: one is the interval along Γ − until the point 0 wherep = λ − , and the other is the interval on the constraint line C from 0 until the turning point (λ For a large initial size of the universe α 0 1 the positive quantity ∆p is small. After the turning point the universe moves back along C up to the point 0, where it can slip to the subcritical solution describing an expanding universe. Let us show that for this motion the spatial constraints (6.21) are always satisfied. The spatial constraints define a domain in the (p-q) plane, where the corresponding functions of curvatures are restricted. This domain is a strip located between the straight linesẐ = Λ and Z = −Λ. We call these lines the upper and lower bounds, respectively. We denote byp ± the coordinates p of the points where the spatial constraintẐ intersects Γ ± lines. At these points one haŝ One can check that In these relations the upper signs stand for the upper bound constraint and the lower signs stand for the lower bound constraint. It is easy to check that the curve representing the evolution of the universe obeying the temporal constraint always lies inside the domain restricted by upper and lower bound lines. In other words, the spatial constraints do not impose any restriction on the evolution of the universe and hence can be ignored.

D. Evolution of the control function χ
Let us now discuss the gravity equations (3.11)-(3.12). As we already mentioned, as a result of the conservation law the second of these equations (the spatial equation) is satisfied if the first (temporal) equation is valid. We rewrite the latter in the form Using expressions for G and T , one gets The control function χ vanishes in the subcritical regime where X is also zero. Equation (6.25) determines the evolution of the control function in the supercritical regime. In such a case one can put ζ = 0 and use the reduced action Taking the variation of the reduced action I c over b and putting b = 1, one gets Using Eq.(4.15), one obtains Combining Eqs.(6.25), (6.26), and (6.29), one can write the equation for the control function χ in the following dimensionless form: (6.30) Hereχ = χ/λ + . This equation determines the time dependence of the control functionχ in the supercritical regime. For a given α(τ ) this is a first-order linear inhomogeneous ordinary differential equation (ODE). This equation can be written in such a form that the control functionχ explicitly depends only on α,

E. Phase diagram
In the previous discussion we focused on the description of the evolution of the universe by using (p, q) planes. Let us now describe this evolution by using the phasespace diagrams. Let us consider a two-dimensional space with coordinates (α, α ). Equation (6.11) can be written in the form This second-order ODE is equivalent to the following system of two first-order equations: Phase diagrams for the system (6.33)-(6.34) are shown in Fig. 3. A dashed line represents the de Sitter solution which approximates a general solution near the turning points.
The dynamics of the universe is described by the system (6.33)-(6.34) with the initial condition  The parameter which has dimensions of length is the critical length corresponding to the limiting curvature Λ.
Then, by using Eq.(3.23) one can write α 0 in the form Here l Pl is the Planck length. An effective curvature radius during inflation that is consistent with observations is usually considered to be in the range 10 5 -10 9 l Pl . If one chooses the critical length to be of the same order of magnitude, then (l Pl / ) 1/2 ≈ 3 · 10 −3 -3 · 10 −5 . Since the entropy of our Universe is large (S/k ≈ 10 90 ), the parameter α 0 is also very large α 0 ≈ 10 25 -10 27 .
The minimal value of dimensionless radius is achieved at the bounce point α b . For every choice of parameters of the system µ and α 0 the bounce point α b can be found from the equation Because α 0 is assumed to be very large, the condition is satisfied for all 0 < µ < µ max , where µ max is close to 1.
In this case the bounce happens very close to unity For example for µ = 0.97 one has α b ≈ 1.01. For smaller values of µ the bounce radius becomes exponentially close to 1. In the range of α 0 ≈ 10 25 -10 27 the corresponding number of e-folds N = ln(α 0 /α b ) is about 57 < N < 62.
Recall that during the supercritical stage the universe first contracts from a 0 to a b . Then, the inflationary stage begins and it expands back to a 0 with the e-fold number N . After that, the Friedmann big bang expansion governed by the standard Einstein equations continues, as depicted in Fig. 4.
In the vicinity of the bounce point the trajectory of the supercritical evolution is very close to the de Sitter spacetime (see Fig. 3). For very large values of α 0 , the supercritical trajectory spends most of its time close to the de Sitter solution. Qualitatively the de Sitter-like behavior happens when the accelerationä changes sign from negative to positive. This is because the effective positive cosmological constant corresponds to repulsive gravity effects. Thus, the criterion of closeness of a supercritical solution to the de Sitter metric is that q > 0. The scale factor a dS when q = 0 can be estimated as follows. At this pointp = 1 − µ, and using Eq.(6.13) one gets For all values of 0 ≤ µ ≤ 1 we have 0.778 < a dS /a 0 ≤ 1, i.e., the de Sitter-like stage always happens very soon after the beginning of the supercritical regime.

F. Effective Lagrangian
Let us note that the Eq.(6.32) coincides with the Euler-Lagrange equation for the following Lagrangian (6.45) Figure 5 shows the potential V as a function of its argument α. Since the Lagrangian (6.43) does not contain an explicit dependence on time τ , the "energy" is conserved. Using the initial conditions one can find The motion with negative energy E in the potential V (α) is bound. In particular, α always has a "left" turning point where it takes the minimal value α b . This conclusion is in agreement with the above general analysis of the evolution of the scale factor in the theory of limiting curvature with linear-in-curvature constraints. Let us notice that the solution also has a "right" turning point where the scale factor reaches its maximal value, . (6.50) If the coefficient µ is not very close to 1, then α * is of order of α 0 and larger than it. It should be noted that before the scale factor reaches α * the solution crosses the line Γ − . If at this point the control function χ vanishes, the solution can slip to the subcritical regime. In the Appendix it is shown that such a solution for χ exists.
In such a case the solution for α leaves its supercritical phase and one gets an expanding Friedmann-Robertson-Walker universe filled with thermal radiation.

VII. A SPECIAL CASE: EINSTEIN CONSTRAINT
In the previous discussion we assumed that the parameter µ was positive. Let us discuss the supercritical solutions in the limiting case where this parameter tends to zero. Using Eq.(6.11), we rewrite Eq.(6.13) in the form The supercritical evolution starts at the point whereq = λ − = −(1 − µ)/(1 + µ) and continues its motion along the constraint line C until it reaches a bounce point in a close vicinity of a pointq = 1. During practically the entirety of this evolution the ratio a/a 0 is of the order of 1. Essential change of the scale factor α occurs only whenq becomes close to 1, so that 1 −q ∼ exp(−2/µ) . This means that such a limiting constraint is equivalent to putting restrictions on the eigenvalues of the Einstein tensor G. We call these restrictions the Einstein constraint. The temporal Einstein constraint is p = Λ/3 =const. The conservation law (3.13) implies that the spatial constraintĜ = Λ is satisfied. The constraint line C is vertical so that λ + = λ − . In the limit µ → 0 the parameterq "jumps" along this line from −1 to 1. The solution of the constraint equation After a bounce at the moment τ = 0 the universe begins to expand.

A. General remarks
In our discussion of the linear-in-curvature constraints we imposed restrictions on the eigenvalues of the linear combinations of the Ricci tensor and the diagonal tensor proportional to the Ricci scalar. Let us now discuss a more general approach where the constraints are composed of functions of scalar invariants constructed from the Ricci tensor 3 . The corresponding constraint can be written in the form This equation establishes a relation between the quantities p and q, defined by Eq. (5.3), and determines a corresponding constraint line C in the (p, q) plane. Let us discuss some general properties of such constraints. Let us assume that ∂f ∂q = 0, so that the equation for the curve C (at least over some its interval) can be written in the form q = q(p). This is nothing but a second-order (nonlinear) equation which is resolved with respect to the second derivative,ä = A(a,ȧ) , (8.2) where the function A(a,ȧ) is determined by the constraint equation (8.1). It may happen that this nonlinear equation has a singular point at which the solution terminates. In such a case the corresponding constraint is singular.
To illustrate this let us use the relatioṅ which directly follows from Eq.(5.3) and which is equivalent to Eq.(8.2). As earlier, we consider the evolution of a radiationdominated universe at the state of contraction which is represented (see Fig. 2) by the interval of line Γ − where q = −p. It starts at some point 1 and continues until it meats the constraint line C at point 0. After this, the solution becomes supercritical and moves along the constraint line C. Since initially the universe contracts, a < 0 at point 0. Let us assume that in its further motion along the constraint a point representing the universe enters the P − domain (see Fig. 6). Since q is negative thereȧ can only decrease and hence remains negative. A turning point of a(t), if it exists, can only be located in the domain P + where q > 0.
Let us assume that the constraint equation does not allow the parameter p to be bigger than p and a solution of the equation f (p, q) = 0 near the point O with Consider a constraint that has a point O with the maximal value of p located in the P + domain (see Fig. 7). Iḟ a is negative at O, then using the above given arguments one can conclude that such a constraint is singular. Let us assume now thatȧ at O is positive. Then, a point representing a solution moves away from O while p is decreasing. This means that if the motion along the constraint starts at point 0 on Γ − , then the solution cannot reach the point O. This happens because before the solution reaches O whereȧ > 0 it first reaches a point on the constraint line C whereȧ = 0. This is a turning point of the solution. At this point a(t) reaches its minimal value. After this the scale factor a(t) increases and a point representing the solution moves back along the constraint curve with decreasing parameter p. In other words a point of the constraint O whereȧ > 0 is not dangerous and the supercritical solution never reaches it.
In what follows we shall not consider singular constraints that do not allow a complete description of the evolution of the universe. Let us note that a "natural" quadratic-in-curvature constraint in which one restricts the Kretschmann invariant belongs to a class of singular constraints. This can be easily seen since the corresponding constraint function is f = p 2 + q 2 − Λ, and p reaches its maximum when q = 0. We shall focus on nonsingular constraints. We shall demonstrate that for a wide class of quadratic-in-curvature constraints there exists a turning point of a(t) located in P + domain, which for a "large" initial size of the scale factor is always very close to the line Γ + where q = p.

B. Quadratic-in-curvature constraints
Let us now discuss a limiting curvature gravity model with quadratic-in-curvature constraints. We denote The most general square-in-curvature expression can be written in the form As it will be explained later it is sufficient to use this constraint in the domain below Γ + where it can be written in the form The equation determines a second-order curve C in the (p, q) plane. We assume that this curve is an ellipse. The general ellipse can be parametrized by its two semiaxes A and B ≤ A, and the angle θ between the large semimajor axis and coordinate axis p. In this parametrization its equation is The coefficients A and B and the angle θ can be expressed in terms of the coefficients c SS , c SR , and c RR and Λ.
In these variables the restriction on the curvature (8.7) implies a restriction on the size of the ellipse and, in particular, on the "length" of its major semiaxis A. A relation between the limiting curvature Λ and A can be easily found. Instead of this it is more convenient to choose from the very beginning the scale defined by A as a limiting curvature parameter and to use A in order to introduce dimensionless quantities that describe our system. Namely, we set We also denote At the moment when the radiation-dominated Friedmann stage matches the evolution along the constraint, Using this initial condition and the constraint (8.11) we get the relation between λ − and the parameters γ and θ The point where dp/dq = 0 (see Fig. 8) has coordinates (p ,q )p = cos 2 θ + γ sin 2 θ, We impose a conditionq >p , that is the point is located above Γ + . This is possible if For γ = 0 the angle θ at which the point (p ,q ) lies on When 0 < γ ≤ 3 − 2 √ 2, the range of θ such that q > p is θ min < θ < θ max , (8.18) The ellipse intersects Γ + at a pointq =p = λ + where . (8.20) Note that for a fixed γ, λ 2 + as a function of θ gets its maximum and minimum values at θ max and θ min , respectively, For small γ,

IX. EVOLUTION ALONG THE CONSTRAINT AND BIG BOUNCE
Let us now discuss the evolution of the universe in the supercritical regime for the quadratic constraint described in the previous section. A point representing the unverse in the (p,q) plane starts its motion atp =q = λ − and moves with increasingp along the ellipse wherē q =q(p). We now demonstrate that this monotonic motion continues until the point reaches the vicinity of p =q = λ + where the scale factor α(τ ) has a turning point. For this purpose, we again use the following the relation, which follows from the definition of the quantitiesp andq It gives .
Here α 0 is the dimensional value of the scale function α at the beginning of the supercritical evolution, that is at p = λ − . According to our assumption α 0 1. In the next section we show that for the quadratic-incurvature constraint the integral in Eq.(9.2) can be calculated explicitly. Now we demonstrate that the general form of the relation Eq.(9.2) allows one to prove that the supercritical solution always has a bounce. Let us note that the integrand in the expression for F is positive in the domain below Γ + wherep < λ + and F is a monotonically increasing function ofp which is logarithmically divergent atp = λ + . One can use Eq.(9.2) to findp as a function of α. Using the definition ofp one has (α ) 2 =pα 2 − 1 .
(9.3) Substitutingp(α) into this relation one obtains an equation that determines the evolution of the scale factor α as a function of time τ . If this function has a minimum α b , then the following condition should be satisfied: For λ − λ + the parameterp is of order of one. Since F grows infinitely nearp = λ + , Eq.(9.4) for large α 0 always has a solutionp b which is located near λ + . This solution determines the size of the universe α b at the turning point.
If we denote by ∆p b < 0 a position of the turning point, then .
The main contribution to this integral comes from the vicinity of its upper limit. This gives In the turning point Then, using Eqs.(9.5) and (9.10), one finds Hence at the turning point α b is always larger than λ −1/2 + and for large α 0 its deflection from this value is small. The existence of the turning point means that the universe has a bounce where its contraction is changed to the expansion.

X. EXACT SOLUTION
Let us demonstrate that for the quadratic-in-curvature constraint one can obtain an explicit expression relatinḡ p and α. For this purpose we again use the equation Let us use Eq.(8.11) to expressq in terms ofp where u, v, w are the following constants (10.4) The integration constant was fixed by the matching condition at the moment when this supercritical solutions starts and wherep = λ − , where . Let us summarize: the evolution of the universe for models with quadratic-in-curvature constraints is quite similar to the case of linear constraints. Namely, there exists a range of free parameters that specify a model for which these constraints are not singular. For such constraints there always exists a bounce and the supercritical solution is symmetric with respect to this moment of time. The size of the universe at the bounce is close to λ −1/2 + and always larger than this parameter. During this supercritical phase the contraction of the universe is replaced by its expansion. The results presented in the Appendix imply that there always exists a solution of the gravity equations for the control function χ that is time symmetric with respect to the turning point. For this solution at the moment when the point representing the expanding supercritical universe crosses Γ − it can slip to the subcritical regime. This subcritical solution describes an expanding Friedmann universe filled with thermal radiation which contains the same entropy as the original collapsing world. The parameter N b = ln α0 α b for the expansion from the bounce point to the beginning of the Friedmann phase is nothing but the e-fold number (for a review of the restrictions on the e-fold number in inflationary models from observations and their cosmological implications see, e.g., Ref. [24] and references therein).

XI. GENERAL CASE
We have discussed cases of linear and quadratic-incurvature constraints which admit rather complete analysis. In this section we demonstrate that under quite general assumptions many of the features of these models are also valid for a wider class of curvature constraints. Namely, we consider constraints constructed from the Ricci tensor invariants. As earlier, we do not include invariants containing derivatives of this tensor. Such a general constraint can be written in the form f (p, q, Λ) = 0. (11.1) Here p and q are defined by Eq.(5.3) and Λ is a parameter defining the limiting curvature. This constraint defines a line C in the (p, q) plane (see Fig. 9). We make the following assumptions: 1. The constraint curve C intersects both lines Γ − and Γ + . We denote the coordinate p at the intersection points by λ − and λ + , respectively.
2. We assume that ∂f ∂q = 0 on a segment of C between Γ − and Γ + , so that on the interval p ∈ [λ − , λ + ] one can express q as a function of p and Λ, q = q(p, Λ).
The last condition implies that λ + > λ − . We denote by ξ the minimal value of dq dp on the interval p ∈ [λ − , λ + ]. Then one has We assume that ξ is not very close to 1, so that the parameters λ + and λ − are of the same order. One can use Eq.(6.5) for the scale factor evolution for a supercritical solution with the constraint (11.1) . (11.3) Here q(p) is defined by the constraint equation. We denote Then, near p = λ + one has Thus, the integrant in F has a pole at p = λ + , so that this integral is logarithmically divergent at this point and Using the same arguments as in the earlier discussion of linear and quadratic-in-curvature constraints, one can conclude that there exists a turning point where a(t) has the minimal values. The scale factor at this point can be found by using Eq. (11.3). After the bounce a point representing the supercritical solution moves back along the line C. To summarize, the supercritical evolution of the universe for the general constraint (11.1) satisfying the conditions 1-3 is qualitatively the same as for linear and quadratic constraints. The control function χ(t) for such a solution is discussed in the Appendix.

XII. DISCUSSIONS
In this paper we studied the evolution of an initially contracting isotropic homogeneous closed universe in the limiting curvature gravity theory. For this purpose, we modified a standard Einstein-Hilbert action by adding terms that restrict the curvature invariants. This was done in such a way that when the curvature is less than the critical one the evolution of the universe follows the standard (unconstrained) cosmological equations. We called this regime a subcritical one. For such a solution, the control function χ in the action vanishes.
After the spacetime curvature reaches its critical value a solution follows along the constraint and the control function χ becomes a nonvanishing function of time. The solution can leave its supercritical regime if the control function becomes zero. To make discussion more concrete we assumed that the contacting universe is initially filled with a thermal gas of radiation and a transition from sub-to supercritical regime occurs when the size of the universe is large in the following sense. If the critical value of the curvature is ∼ 1/ 2 , then we required that the size a 0 of the contracting universe at the moment when its curvature reaches this critical value obeys the condition a 0 / 1. There is a freedom in the choice of the term in the action which controls and restricts the growth of the curvature. We studied two types of constraints. The first class are constraints that are linear in eigenvalues of the Ricci tensor. Such constraints are represented by a straight line in (p, q) plane. If this line crosses Γ ± where q = ±p and 0 < dp/dq < 1 for it, then we demonstrated that the evolution of the universe with such an inequality constraint is the following. After the contracting universe reaches the critical curvature and the evolution becomes supercritical, its acceleration parameter q =ä/a quite soon becomes positive. If the scale factor at the transition point α 0 is large, then its further motion is very close to the motion of a contracting de Sitter universe. It has a bounce where the scale factor a(t) has the minimal value a b close to and begins expanding. Both the scale function a(t) and the control function χ(t) are symmetric under reflection with respect to the turning point. The function χ(t) becomes zero again when the size of the expanding universe becomes equal to a 0 . After this point the solution leaves the constraint and it describes an expanding universe filled with thermal radiation. The entropy of this radiation is the same as that during the contraction phase.
The second class of constraints that we discussed in this paper are quadratic in curvature. We demonstrated that there exists a wide variety of such (nonsingular) constraints that guarantee that solutions are complete, that is, they do not break at a finite time. In the (p, q) plane the constraint curves are ellipses with two parameters the angle θ characterizing the orientation of the ellipse, and the ratio of its semiminor and semimajor axes γ. The size of the major semiaxis characterizes the limiting curvature value. We showed that if θ and γ obey some inequalities the corresponding constraint is nonsingular. The evolution of the universe for such nonsingular quadratic constraints is similar to the case of linear-in-curvature) constraints. After the universe reaches the point where its curvature becomes critical, the solution evolves along the constraint. During this supercritical phase it reaches a point of bounce after which the scale function grows. The control function χ can become zero again at this phase and the universe can leave its supercritical regime. After this one has an expanding universe filled with thermal radiation which follows the corresponding solution of the Einstein equations.
Let us emphasize that these results allow a natural generalization. In Sec. XI we demonstrated that they can be easily extended to the case when a constraint is not linear or quadratic in the curvature but is described by a quite general function of it. 4 The qualitative behavior of the supercritical solutions in such models remains qualitatively the same. The solutions predict a bouncing point, when the contracting universe transitions to expansion.
In our discussion we assumed that a contracting universe is filled with thermal radiation. This simplified our analysis at one point, where we calculated the value of the scale parameter a 0 at the moment of transition of the solution from the sub-to supercritical regime. The value of a 0 can be easily found for any other choice of the equations of state. This changes nothing in the further supercritical evolution of the universe, which was the main point of our discussion. Another assumption was that our contracting universe is closed. The other two cases, k = 0 and k = −1 where the universe is open , can be analyzed similarly. The main difference is that a supercritical solution for a(t) does not have a turning point but it can reach zero value. However, this does not mean that one has a physical singularity at this point. The curvature invariants remain finite and bounded and the singularity of the solution is a reflection of a "bad" choice of the coordinates. The situation here is similar to the case of de Sitter model when coordinates with open space slices are chosen. One can expect that by using proper coordinates one can further trace the evolution of the supercritical solution. It would be interesting to study these cases in detail and to confirm (or disprove) that there is also a bounce for these universes.
Many gravity models have been proposed in the literature that describe an inflationary stage of the Universe. Some of these models involve either higher-order-incurvature terms or higher derivatives, or both. As a consequence, these models are typically prone to instabilities [19]. Some of these instabilities are related to the presence of ghosts (see however the discussion in Ref. [25]). Complications with ghosts can be avoided in some versions of nonlocal higher-derivative theories of gravity, and cosmologically viable models admitting nonsingular bouncing solutions can be constructed [16,26,27]. The analysis of the stability of cosmological solutions is a nontrivial problem in both ghost-free higher-derivative theories of gravity and systems with constraints [19].
In the models of limiting curvature gravity discussed in the present paper a set of pairs of Lagrange multipliers χ i and ζ i entering in a specific combination was introduced. As a result, as soon as some function of curvature invariants does not accede its limiting value, the gravity equations are exactly those of the pure Einstein theory. This means that during a subcritical stage all of the degrees of freedom and physical effects are exactly the same as in general relativity. No extra instabilities and ghosts appear. During the supercritical stage the metric evolution 4 Here we do not consider more general invariants constructed from the curvature and its derivatives.
is governed by the constraints. Matching conditions provide us with the initial data for the evolution of the metric and the Lagrange multipliers. Further evolution is unambiguous and respects the property of limiting curvature. If one includes a constraint involving the Weyl tensor, then the growth of all relevant curvature invariants can be bounded even if instability modes appear. Of course, in more realistic models one has to constrain all kinds of curvatures and take into account anisotropic and other deviations from the background geometry [19,28]. The control fields χ are very special. They identically vanish in the subcritical regime and this property allows one to obtain the uniquely specified initial conditions for them at the beginning of the supercritical regime. In the latter case the control field obeys the linear inhomogeneous equation. The initial conditions and the inhomogeneous term completely fix the solution for χ. The control fields do not bring extra degrees of freedom to the system. In this sense they are not dynamical. There is a well-known generic problem of bouncing cosmological models: the growth of the anisotropy at the stage of contraction. Even if the anisotropy is initially small and it can be described as a perturbation of isotropic homogeneous space, its amplitude for a physically reasonable equation of state grows fast so that during the contraction at some its stage the anisotropy would become to affect the dynamics of the contracting universe. One can expect that in the models with limiting curvature this anisotropy growth could be suppressed by a proper choice of the constraints that contain not only Ricci-tensor but also Weyl-tensor invariants. When properly included, the corresponding constraint would not allow infinite anisotropy growth. It is interesting and important to check whether this is really so.
One might interpret the obtained results as follows. After the curvature of the contracting universe reaches its maximal value the matter does not contribute to the growth of the curvature. Instead, the further growth of its stress-energy tensor is compensated by the generation of the control field χ. After passing the bounce and reaching the point of slipping back to the subcritical phase, the hidden thermal radiation (with its entropy) simply reappears. In this sense, the thermal state of the inflating universe arises without an additional reheating. Certainly this and other interesting features of the bouncing cosmologies in the limiting curvature gravity models require further detailed analysis.
Appendix: Evolution of the control function χ Let us consider the general curvature constraint that was discussed in Sec. XI. We now discuss the gravity equations that are obtained by variation of the action including this constraint over the metric. During the supercritical stage the metric evolution is governed by the constraint. The gravity equations, in fact, describe evolution of the control function χ. As we discussed earlier there is only one independent equation which can be obtained by variation of the dimensionally reduced action over the lapse function b. The constraint (11.1) can be obtained from the action (A.1) by its variation over the control function χ. Let us emphasize that in order to derive a complete set of gravity equations one should substitute general expressions (2.6) for p and q which contain the gauge function b.
In what follows we assume that conditions 1-3 formulated in Sec. XI are satisfied and one can use Eq.(11.3) to find p on the constraint line C as a function of the scale parameter a.
Let us consider the gravity equation we get the first-order linear inhomogeneous differential equation where y = ln a and y 0 = ln a 0 . At the matching point the control function vanishes. This condition fixes the integration constant Ω 0 = 0. The solution (A.12)-(A.13) is finite for all a between a 0 and the bounce radius a b . This fact is not evident and needs special analysis because both functions U and W have a pole near the bounce point. This happens because the function p − 1 a 2 =ȧ 2 a 2 in their denominators vanishes at the bounce, whereȧ = 0. In order to prove that these poles do not lead to singularities for ω at a b , we analyze Eq.(A.8) in its vicinity. Let us analyze the asymptotic of this equation when a → a b . Let x = ln(a/a b ) be a small parameter. Then, using the expansion near the turning point, one gets (A.14) Here q b = q(a b , a 0 , Λ) and p b = p(a b , a 0 , Λ). As soon as we know the evolution of the scale factor a(t) along the constraint, we know p b and q b and therefore can determine the evolution for the control functions ω and χ.
For small x one has where the constant C 1 is defined by the asymptotic behavior of W at the bounce point Because at the bounce point p b ∼ q b ∼ λ + > 0 and a 2 0 λ + is a very big number, the constant C 1 is positive and big too. Equation (A.17) shows that ω = C 1 + C 2 |x| + O(x). (A.18) The integration constant C 2 can be determined using the matching condition at a 0 . If t b is the moment of bounce, then near this point one has a ≈ a b [1 + q b 2 (t − t b ) 2 ]. So that The second derivative of Y is d 2 Y dp 2 = 2 (Y + p)(1 + dq(p) dp ) (p − q) 2 .
(A. 25) Since dq/dp is positive on the interval p ∈ [λ − , λ + ] the second derivative of Y is positive at the matching point. This means that in the vicinity of λ − for p > λ − both Y and dY /dp are positive. To prove that Y (p) is positive on the whole interval (λ − , λ + ) we assume the opposite. Namely, we assume that Y becomes negative at some point on this interval. This means that there exists point a p 1 ∈ (λ − , λ + ), where Y vanishes again. This may happen only if Y reaches its maximum at p 2 ∈ (λ − , p 1 ). In this case dY /dp| p2 = 0 and, hence, there exists p 3 ∈ (λ − , p 2 ) where dY /dp has a maximum. At this point d 2 Y /dp 2 | p3 = 0, dY /dp| p3 > 0, and Y (p 3 ) > 0. This conclusion is in contradiction with Eq.(A.25) evaluated at the point p 3 . Thus, our assumption that Y may vanish inside the interval (λ − , λ + ) leads to a contradiction and we must conclude that both dY /dp and Y are non-negative during the whole supercritical stage. Hence, W is negative for a < a 0 and W = 0 at a = a 0 . Using this property, one can show that ω is negative during the whole evolution along the constraint. The function χ differs from ω only by a factor ∂f /∂q and, hence, it does not vanish as well. This means that the control function χ does not vanish during a supercritical stage except for the initial and final matching points corresponding to the scale factor a 0 .