Bounce cosmology from $F(R)$ gravity and $F(R)$ bigravity

We reconstruct $F(R)$ gravity models with exponential and power-law forms of the scale factor in which bounce cosmology can be realized. We explore the stability of the reconstructed models with analyzing the perturbations from the background solutions. Furthermore, we study an $F(R)$ gravity model with a sum of exponentials form of the scale factor, where the bounce in the early universe as well as the late-time cosmic acceleration can be realized in a unified manner. As a result, we build a second order polynomial type model in terms of $R$ and show that it could be stable. Moreover, when the scale factor is expressed by an exponential form, we derive $F(R)$ gravity models of a polynomial type in case of the non-zero spatial curvature and that of a generic type in that of the zero spatial curvature. In addition, for an exponential form of the scale factor, an $F(R)$ bigravity model realizing the bouncing behavior is reconstructed. It is found that in both the physical and reference metrics the bouncing phenomenon can occur, although in general the contraction and expansion rates are different each other.

that big crunch singularities of negative-energy (Anti-de Sitter) bubbles in the multiverse can be removed. Also, the behavior of bounce in anisotropic cosmology in F (R) gravity [30] and a bounce in modified gravty theories [31] have recently been discussed. We further mention that the form of F (R) leading to non-singular bounce cosmologies has been derived in Ref. [32]. Furthermore, bounces in gravity inspired by string theories [33] and non-local gravities [34] have been studied.
In addition, it has recently been revealed that a massive graviton can lead to the current cosmic acceleration. At the early stage, the Fierz-Pauli (FP) action [35] was considered to describe a linearized or free theory of massive gravity (for reviews, see, for example, [36,37]). Recently, the de Rham, Gabadadze, Tolley (dRGT) theory [38,39] and the Hassan-Rosen (HR) theory [40], which are non-linear massive gravity theories, have been proposed. These theories have two desirable properties: One is there is not the Boulware-Deser (BD) ghost [41,42]. Another is in the massless limit of the mass of massive graviton the van Dam-Veltman-Zakharov (vDVZ) discontinuity [43] can be screened through the Vainstein mechanism [44]. The latter is a similar feature appearing in the Galileon models [45] due to the operation on the Dvali-Gabadadze-Porrati (DGP) brane world scenarios [46]. Currently, in various aspects, massive gravity and bi-metric gravity have extensively been studied in the literature .
However, thanks to recent works [70], it has been found that in such non-linear massive gravity theories, the flat homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) universe, which is supported by various cosmological observations, cannot be stable. Hence, massive gravity theories in the context of general relativity explained above, which is called "the massive general relativity (GR)" in the literature, have been extended, for instance, an extended version of the dRGT theory [71], a massive bi-metric F (R) gravity theory [72,73], a new massive F (R) gravity [74] proposed very recently, a scale invariant theory with a dilaton field, or which is called the "Quasi-Dilaton" massive gravity (QMG) [75,76] and its extended versions [77], and mass varying scenario in which a massive graviton mass depends on a dynamical scalar field [78].
In this paper, with the procedure proposed in Ref. [79], which corresponds to a kind of simpler and more useful reconstruction method made by developing the formulation in Ref. [80], we derive F (R) gravity models in which bounce cosmology can occur. In particular, in the flat FLRW universe we perform the analysis for two cases that the scale factor is described by exponential and power-law forms. We study the perturbations from the background solutions and explicitly explore the stability conditions for these models to be stable. In addition, we investigate an F (R) gravity model with the scale factor having a sum of exponentials form, where the unification of the bouncing behavior in the early universe and the late-time cosmic acceleration can be realized. Furthermore, in the FLRW universe with non-zero spatial curvature, for an exponential form of the scale factor, we reconstruct F (R) gravity models in which a function of F (R) is expressed by a polynomial in terms of R. Also, for the scale factor with an exponential form, in the flat FLRW universe, we build F (R) gravity models by using the reconstruction method [81] and explore the stability conditions. We also reconstruct an F (R) bigravity model realizing bounce cosmology. Incidentally, the bouncing behavior and cyclic cosmology in extended non-linear massive gravity [82] and bounce cosmology in bigravity [83] have been investigated.
Here, we clarify our purpose of this study. As the first step, in this work we reconstruct F (R) gravity and F (R) bigravity models with the bouncing behavior. In particular, for F (R) gravity, we build models in which not only the bounce in the early universe but also the late-time cosmic acceleration occurs and examine the stability of these models. At the current stage, these models are still toy models. However, we note that some of the considered models with an R 2 term are known to be viable models for the early-time inflation. Moreover, it should be emphasized that bounce cosmology may be a natural part of the complete and viable history of the universe. This is the reason to study better such cosmologies. Our final goal is to construct the so called viable F (R) gravity and F (R) gravity models, in which all the cosmological various processes of expansion history of the universe with a bounce can be realized. This developed subject should be executed as another separate work in the near future. We use units of k B = c l = = 1, where c l is the speed of light, and denote the gravitational constant 8πG N by GeV. The paper is organized as follows. In Sec. II, we explain a reconstruction method of F (R) gravity. With this procedure, we derive F (R) gravity models realizing bounce cosmology in Sec. III. Furthermore, in Sec. IV we examine the stability of the reconstructed F (R) gravity models. In Sec. V, we also build an F (R) model where both the bounce in the early universe and the late-time cosmic acceleration can occur in a unified manner. In Sec. VI, we investigate an exponential form of the scale factor for the non-zero spatial curvature, while in Sec. VII, we explore it for the zero spatial curvature. Moreover, in Sec. VIII we reconstruct F (R) bigravity models in which the bouncing phenomenon can happen. In Sec. IX, conclusions are presented.

II. RECONSTRUCTION METHOD OF F (R) GRAVITY
In this section, we explain the reconstruction method of F (R) gravity [79]. The action of F (R) gravity with matter is expressed as with L M the matter Lagrangian and Ψ M matter fields.
In the flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe, the metric is given by with a the scale factor.
Here, we introduce the number of e-folds defined by N ≡ ln (a/a * ), where a is a scale factor and a * is a value of a at a time t * . When we take t * = t 0 with t 0 the present time and and a * = a 0 at t = t 0 , we can also define the redshift z as z ≡ a 0 /a − 1. Moreover, the Hubble parameter is given by H ≡ȧ/a, where the dot denotes the time derivative of ∂/∂t, and we describe it by using a function ofg(N ) as H =g(N = − ln (1 + z)). Furthermore, we write H 2 as H 2 ≡ G(N ) =g 2 (N ) with G(N ) a function of N . With the quantities defined above, in this background the Friedmann equation reads Here, ρ M is the sum of energy density of all matters assumed to be fluids with a constant equation of state w i defined as w i ≡ P M i /ρ M i , where the subscription "i" shows the label of the fluids and ρ M i and P M i are the energy density and pressure of the i-th fluid, respectively, ρ M i0 is a constant, and the prime denotes the derivative with respect to N as G ′ (N ) ≡ dG/dN and G ′′ (N ) ≡ d 2 G/dN 2 .

III. F (R) GRAVITY REALIZING BOUNCE COSMOLOGY
In this section, we study the cosmological background evolutions in the matter bounce cosmology and reconstruct F (R) gravity models realizing it.

A. Exponential model
We examine the case that the scale factor is expressed by an exponential form. For instance, we consider a bouncing solution which behaves as (3.1) Here, α is a constant with the dimension of mass squared ([Mass] 2 ). In the following, we set N ≡ ln a(t)/a(t = 0), where a(t = 0) = 1 because we study the bouncing behavior around t = 0. We now use the reconstruction in Ref. [79]. From Eq. (2.3), we solve the following differential equation: Here, we have neglected a contribution from matters and G(N ) = H(N ) 2 and the scalar curvature R is given by For the model (3.1), we find 5) and therefore Then, Eq. (3.2) has the following form: A solution of (3.7) is given by In Fig. 1, we show the behavior of the Hubble parameter in the second relation in (3.4) for α = 1/2 around a bounce at t = 0. From this figure, we see that before the bounce (t < 0), H < 0, while after it (t > 0), H > 0. Thus, the bouncing behavior occurs.
It should clearly be mentioned that the metric with the scale factor (3.1) does not have a finite maximum in the Riemann curvature, whereas the Riemann curvature takes its minimum by modulus in the bounce epoch. Thus, this space-time is irrelevant to the thing necessary to remove a cosmological singularity. By the same reason as written above, in the model described by Eq. (3.8) the Starobinsky inflation [26] cannot be realized. When the Starobinsky inflation occurs, the scale factor at the slow-roll inflationary stage is given by a(t) ∝ exp H 1 t − M 2 t 2 /12 . Also, we explore the stability with respect to tensor perturbations, namely, the required condition F ′ (R) > 0. It follows from the second relation in (3.5) with the first one in (3.4) and Eq. (3.8) that we have F ′ (R) = (2/α) (R − 36α) = 48 2αt 2 − 1 . Hence, when a bounce occurs at t = 0, we find F ′ (R) < 0. We also see that F ′ (R) = 0 at R = 36α. As a result, the bounce of the scale factor in Eq. (3.1) occurs in the unphysical regime of a negative effective gravitational constant, so that graviton can become a ghost. In addition, for α > 0, F ′′ (R) = 2/α > 0, and therefore the stability condition for the cosmological perturbations [9,84] can be satisfied. Moreover, in Refs. [85] and [86] it has been found that even at the classical level, it is not able to pass the point in which F ′ (R) = 0 for a finite R because in a generic solution a strong anisotropic curvature singularity appears.

B. Power-law model
On the other hand, it is known that for the case in which the scale factor is expressed by a power-law model, given by whereā( = 0) a constant,t is a fiducial time, and q = 2n with n is an integer, a power-law model of F (R) gravity would be reconstructed. In this case, we acquire With Eqs. (3.3), (3.10), (3.11) and G = H 2 , we find Around the bounce behavior, we have N ≃ 0. Hence, by adopting an approximation e N ≃ 1 to Eq. (3.13), we obtain where we have also neglected the matter contributions. As a solution, we have withF ( = 0) a constant. It has first been shown in Ref. [87] that there exist power-law solutions for such a monomial form of F (R). In Fig. 2, forā = 1.0, q = 2 with n = 1, andt = 1, we depict the behavior of the Hubble parameter in Eq. (3.11) around a bounce at t = 0. It follows from this figure that before the bounce (t < 0), H < 0, whereas after it (t > 0), H > 0, similarly to that in Fig. 1. As a result, the bouncing behavior happens.
We also remark that in the matter bounce cosmology with two fields [19], for a ∝ (t −t) s , cosmological background evolutions consist of the following four phases: (i) matter contraction phase, (ii) the Ekpyrotic contraction phase, (iii) bounce phase, and (iv) fast-roll expansion phase. In the matter contraction, the Ekpyrotic contraction, and fast-roll expansion phases, the scale factor a behaves as power-law type in Eq. (3.9), while in the bounce phase, a evolves as exponential type in Eq. (3.1). For the matter contraction phase, we find s = 2/3, for the Ekpyrotic contraction phase, s would not be set to a specific value, whereas in the fast-roll expansion phase, we have s = 1/3. On the other hand, for the Ekpyrotic contraction phase, if the scale factor is described by Eq. (3.1), we see that α is determined by the detailed physics on micro scales of the bounce process. Finally, it should be emphasized that a specific case of the matter bounce scenario [11,19] investigated by Brandenberger et al. is able to be realized also in F (R) gravity.

IV. STABILITY OF THE SOLUTIONS
In this section, with the procedure of the first reference in Ref. [10], we examine the stability of the solutions in F (R) gravity models obtained in Sec. III.
We suppose that a solution of Eq.
Here, we note that N (≥ 0) is equal to or larger than 0. (This is clearly seen from the first equation in (3.4) with α > 0 and Eq. (3.10) withā > 0.) By substituting this expression into Eq. (2.3), we find where the values of dF (R)/dR, d 2 F (R)/dR 2 and d 3 F (R)/dR 3 are the ones at R = 3G ′ b (N ) + 12G b (N ) following from Eq. (3.3). Thus, the stability conditions J 2 /J 1 > 0 and J 3 /J 1 > 0 can be written as A. Stability of the exponential model In the case that the scale factor is described by an exponential form in the exponential model, with Eqs. (3.3), (3.5) and (3.8) we see that G b = 4αN and therefore the first condition in (4.5) reads 6 − 1/N > 0. Moreover, regarding the second condition in (4.6), the quantity on the left-hand side is equal to zero. In other words, the quantity J 3 /J 1 is not negative. Consequently, if N < 0 or N > 1/6, the solution could be stable. The latter condition can be satisfied because N has to be much larger than unity. Thus, the exponential model of the scale factor could be stable.

B. Stability of the power-law model
When the scale factor has a power-law form, given by Eq. (3.9), the first stability condition (4.5) becomes whereas the second stability condition (4.5) reads Here, in deriving Eqs. (4.7) and (4.8), we have used e N ≃ 1 in those numerators. From Eq. (4.7), we see that if β (1 − 2/q) + 5/q − 4 > 0, the first stability condition can be satisfied. Furthermore, it follows from Eq. (4.8) that for q < 1/2 or q > 2, the second stability condition can be met.

V. UNIFIED F (R) MODEL OF BOUNCE AND THE LATE-TIME COSMIC ACCELERATED EXPANSION
In this section, we reconstruct an F (R) model where not only the bouncing behavior in the early universe but also the late-time accelerated expansion of the universe can be realized in a unified manner.
A. Sum of exponentials model As a concrete model, we investigate a sum of exponentials form for the scale factor Here, we again note thatt is a fiducial time. In this model, for the limit t/t → 0, i.e., in the early universe, we obtain a → e Y , which is equivalent to a = e αt 2 with α = 1/t 2 in Eq. (3.1), and hence the bouncing behavior can occur. While, in the limit t/t ≫ 1, we find a → e Y 2 and henceä = 4 (1/t) 2 Y 3 + 4Y 2 e Y 2 > 0. Consequently, the late-time accelerated expansion of the universe can be realized. In the following, we analyze cosmological quantities around t = 0 in order to examine bounce cosmology. With N = ln a and H =Ṅ , the form of a in Eq. (5.1) leads to Accordingly, by applying Eqs. (5.7) and (5.8) to Eq. (3.2) and providing that contributions from matter are negligible, we acquire We find a solution of this equation as Here, the reason why the solution in Eq. (5.10) includest is that the dimension of the F (R) form is adjusted to be mass squared ([Mass] 2 ). From Eq. (5.10) with Eqs. (5.5) and (5.8), we acquire F ′ (R) = 2t 2 R − 18/t 2 = 24 (t/t) 2 − 1 . Accordingly, when a bounce happens at t = 0, F ′ (R) < 0. We also see that F ′ (R) = 0 at R = 18/t 2 . Consequently, for the scale factor in Eq. (5.1), the bounce is realized in the regime when a effective gravitational constant is negative, namely, graviton is a ghost. On the other hand, since F ′′ (R) = 2t 2 > 0, the cosmological perturbations can be stable [9,84] . In Fig. 3, we display the behavior of the Hubble parameter in the first equality witht = 1 in Eq. (5.4) around a bounce at t = 0. In this figure, before the bounce (t < 0), we have H < 0, and after it (t > 0), we obtain H > 0. This is the same behavior as Figs. 1 and 2, and therefore the bouncing behavior emerges.

B. Stability of the sum of exponentials model
For the double exponential model in Eq. (5.1), the stability condition (4.5) reads This is satisfied if N > ln (5/2). Since N has to be much larger than unity, this condition can be met. Moreover, for Eq. (5.1), the left-hand side of the inequality (4.6) becomes zero. Presumably, if we include higher order term in Y , the left-hand side of the inequality (4.6) might be non-zero, and therefore that we can have some conditions on N , although it might be quite difficult to execute the nvestigations analytically. Hence, it would be expected that such a condition could be satisfied because of the large value of N . It follows from the above considerations that the sum of exponentials model could be compatible with the stability conditions.

VI. EXPONENTIAL FORM OF THE SCALE FACTOR FOR THE NON-ZERO SPATIAL CURVATURE
In Secs. III A and IV A, we have seen that in the flat FLRW universe, an exponential form of the scale factor and the resultant second order polynomial model of F (R) gravity could be a stable theory realizing the bounce cosmology. In this section, we examine an exponential form of the scale factor for the non-zero spatial curvature, namely, in the non-flat FLRW universe.
From this set of parameters, we see that where w DE is the equation of state of the dark energy component defined by From this set of parameters, we find that where f is a function of the parameters k, σ, τ , λ, α 0 , α 1 , and α 2 .
• Type III From this set of parameters, we obtain It should be noted that for the form of the function F (R) in Eq. (6.9), there is no solution other than the de Sitter solution, if the cosmic curvature k is zero (and also f = 0). We also remark that as the scale factor a(t) satisfying the above solutions, a more general expression can be described by replacing t → t − t 1 with t 1 another fiducial time, i.e., Similarly, F (R) can be generalized as any function of the form where β j (j = 1, . . . , l) and α i (i = 0, . . . , m) are constants.

B. Model consisting of an inverse power-law term
Next, we investigate the function F (R) expressed as [8] F With the similar procedure developed in the preceding subsection, we obtain the following restrictions on the parameters In addition, the following condition has to be met It is easy to rewrite this equation in the following form From this equation, we obtain the restrictions on the parameters From this set of parameters, we see that It should be cautioned that in the model in Eq. (6.15), the late-time cosmic acceleration which is accepted from the quantum field theoretical point of view cannot be realized because its de Sitter solution exists in the unstable region where F ′′ (R) < 0.
We can consider a slightly different form of the function F (R) as With the similar procedure developed in the preceding subsection, we obtain the following restrictions on the parameters In addition, the following condition has to be met It is easy to rewrite this equation in the following form For this type of a function F (R), we have more complicated solutions, but we can impose additional restrictions and find a set of parameters for which F ′′ (R) > 0. For example, 20) or From this set of parameters, we see that In summary, in this section, for the FLRW universe with non-zero spatial curvature, when the scale factor is given by an exponential form in Eq. (6.6), we have reconstructed a second order polynomial F (R) model in terms of R and an F (R) model consisting of both a term proportional to R and an inverse power-law term. It has been found that the de Sitter solution can exist for the case with non-zero spatial curvature. Related to these consequences, as noted in Introduction, we again mention that if the spatial curvature is positive, i.e., k(> 0), and a massive scalar field exists, the scale factor as well as the Riemann curvature can perform the bouncing behaviors [25]. Moreover, when the spatial curvature has a non-zero value, namely, k( = 0), in the Starobinsky model [26] there is a solution where the scale factor can behave a bounce.

VII. EXPONENTIAL FORM OF THE SCALE FACTOR FOR THE ZERO SPATIAL CURVATURE
In the study of the bouncing behavior with an exponential form of the scale factor, it seems that another version of the reconstruction method (with an auxiliary scalar field) is more suitable. Hence, in this section we apply it to the derivation of F (R) gravity models realizing bounce cosmology.

A. Reconstruction method of F (R) gravity
When the scale factor is given by an exponential form in Eq. (6.6), with the reconstruction method [81], we find F (R) gravity models with realizing the bounce cosmology. By using proper functions P (t) and Q(t) of a scalar field t which we identify with the cosmic time, the action in Eq. (2.1) can be represented as The variation with respect to t yields from which it is possible to solve t in terms of R as t = t(R). By substituting t = t(R) into Eq. (7.1), F (R) can be written as F (R) = P (t(R))R + Q(t(R)) .
With Eq. (6.2), Q(t) is given by Taking into account Eq. (7.4), from Eq. (6.3) we have the differential equation where we have used the expression of the Hubble parameter H =ȧ/a of the first equality in (6.7). There are two different cases.
The general solution of Eq. (7.5) is given by where c 1 and c 2 are constants. From Eq. (7.4), we have It follows from Eq. (7.2) that By solving Eq. (7.3), we find the most general form of F (R)

Case 2:
The general solution of Eq. (7.5) is given by From Eq. (7.4), we obtain From Eq. (7.2), we obtain By solving Eq. (7.3), we acquire the most general form of F (R) ± (R) + B (2) We caution that F Hence, we have executed a reconstruction of F (R) gravity for the scale factor in Eq. (6.6), so that we have been able to build several types of F (R) gravity theories realizing bounce cosmology. In Fig. 4, we plot F as a function of R with the parameters c 1 = 1, c 2 = 0; 1; 2; 3, λ = 1, σ = 1, and τ = −1.

B. Stability of the solutions
Next, we explore the stability of the obtained models. However, there are several problems associated with the large arbitrariness in the choice of the coefficients and the unwieldy of emerging relations. As an example, we study the stability of solutions F For the scale factor in Eq. (7.7), the stability conditions (4.5) and (4.6) can be written as follows.
• Case I (c 1 = 0, c 2 = 0) As a result, for case I, if N > 0.251224, whereas for case II, when N > 0.0701889, both stability conditions can be met. Since the value of N is much larger than unity, these stability conditions can be satisfied. Thus, we find that for the scale factor in Eq. (6.6), the model F (1) − (R) is stable. In Fig. 5, we illustrate the behavior of the Hubble parameter in the second relation with λ = 1 in (7.8) around a bounce at t = 0. From this figure, it is observed that before the bounce (t < 0), H < 0, and after it (t > 0), H > 0. This behavior is the same as Figs. 1-3. Accordingly, the bouncing behavior is realized.

A. F (R) bigravity
We start with reviewing F (R) bigravity proposed in Ref. [72]. The consistent model of bimetric gravity, which includes two metric tensors g µν and f µν , was proposed in Ref. [89]. It contains the massless spin-two field, corresponding to graviton, and massive spin-two field. It has been shown that the Boulware-Deser ghost [42] does not appear in such a theory.
We consider the following action: Here, R (g) is the scalar curvature for g µν , R (f ) is the scalar curvature for f µν , m is constant mass of a massive graviton, M eff is defined by 1 with M g and M f constants, andβ j (j = 0, . . . , 4) are constants. Moreover, ϕ and ξ are scalar fields, and V (ϕ) and U (ξ) are the potential of ϕ and ξ, respectively. Furthermore, a tensor g −1 f is defined by the square root of g µρ f ρν , that is, For a general tensor X µ ν , e n (X)'s are defined by where [X] expresses the trace of arbitrary tensor X µ ν : [X] = X µ µ . By the conformal transformations g µν → e −ϕ g J µν and f µν → e −ξ f J µν , the action (8.1) is transformed as Note that the kinetic terms for ϕ and ξ vanish. By the variations with respect to ϕ and ξ as in the case of convenient F (R) gravity [9], we obtain We should mention, however, that it is difficult to solve Eqs. (8.4) and (8.5) with respect to ϕ and ξ explicitly. Therefore, it might be easier to define the model in terms of the auxiliary scalars ϕ and ξ as in (8.3). We now explore the cosmological reconstruction program following Ref. [72] but in a slightly extended form. For simplicity, we start from the minimal case:β 0 = 3,β 1 = −1,β 2 =β 3 = 0, andβ 4 = 24. In order to evaluate δ g −1 f , we examine two matrices M and N , which satisfy the relation M 2 = N . Since δM M + M δM = δN , we find tr δM = 1 2 tr M −1 δN . For a while, we investigate the Einstein frame action (8.1) in the minimal case and we neglect the contributions from matters. By the variation with respect to g µν , we have On the other hand, by the variation with respect to f µν , we find We should note that det √ g det g −1 f = √ f in general. The variations of the scalar fields ϕ and ξ are given by Here, g ( f ) is the d'Alembertian with respect to the metric g (f ), and the prime means the derivative of the potential in terms of the argument as V ′ (ϕ) ≡ ∂V (ϕ)/∂ϕ and U ′ (ξ) ≡ ∂U (ξ)/∂ξ. By multiplying the covariant derivative ∇ µ g with respect to the metric g with Eq. (8.6) and using the Bianchi identity 0 = ∇ µ g 1 and Eq. (8.8), we obtain Similarly, by using the covariant derivative ∇ µ f with respect to the metric f , from (8.7) we find In case of the Einstein gravity, the conservation law of the energy-momentum tensor corresponds to the Bianchi identity. In case of bigravity, however, the conservation laws of the energy-momentum tensor of the scalar fields are independent of the Einstein equation. The Bianchi identities give Eqs. (8.9) and (8.10) independent of the Einstein equation.

B. Cosmological bouncing models
Next, we construct cosmological bouncing models. The physical metric, where the scalar does not directly couple with matter, is given by multiplying the scalar field to the metric in the Einstein frame in (8.1): g J µν = e ϕ g µν . In the bigravity model, there appears another (reference) metric tensor f µν besides g µν . In our model, since the matter only couples with g µν , the physical metric could be given by g J µν . In our formulation, it is convenient to use the conformal time description. The conformally flat FLRW universe metric is given by This equation (8.22) with g J µν = e ϕ g µν shows e ϕ(t) a(t) 2 =ã(t) 2 , that is, ϕ = −2 ln a(t) + lnã(t). By using (8.21), we find (8.23) Here,H ≡ 1 a dã dt . In the following, by making the choice a(t) = b(t) = 1, we explicitly construct the model generating the bouncing behavior. We should remark that the choice a(t) = b(t) = 1 satisfies the constraint (8.16).
When a(t) = b(t) = 1, the Einstein frame metric g µν expresses the flat Minkowski space, although the metric we observe is given by g J µν . Equations (8.17), (8.18), (8.19), and (8.20) with (8.23) are simplified as follows We should note that both ω(t) and σ(t) are positive and hence there does not appear any ghost in the theory. We now study the bouncing solutionã withᾱ a positive constant. SinceH we find and ω(η) = 12ᾱM 2 g η 2 ,Ṽ (η) = −12ᾱ 2 η 2 , σ(ζ) = Consequently, for Eq. (8.29), the solutions in (8.32) can be obtained. Moreover, the exponential form of the scale factor in Eq. (8.29) is equivalent to that in Eq. (3.1), which can lead to the bouncing behavior. This means that in the flat FLRW universe, for an exponential form of the scale factor in F (R) bigravity, bounce cosmology can be realized, similarly to that in F (R) gravity, as demonstrated in Sec. III A. For the case that the scale factor has an exponential form in Eq. (8.29), in terms of the physical metric, the bouncing behavior in F (R) bigravity is the same as that in F (R) gravity. On the other hand, for this case, in the reference metric, i.e., the fiducial metric existing only in F (R) bigravity, it is clearly seen from Eqs. (8.29) and (8.31) that also in this reference metric, the bouncing behavior can occur, but the contraction and expansion rates are different each other. In the physical metric,H =ȧ/ã ∼ 2ᾱt as given by Eq. (8.30), while in the reference metric,ċ/c = 2It/ 1 + It 2 with I ≡ 12ᾱ 2 M 2 g / m 2 M 2 eff . The ratio ofH toċ/c reads R ≡H/ (ċ/c) ≃ᾱI −1 1 + It 2 . Thus, when It 2 ≫ 1, the contraction and expansion rates in the physical metric are much larger than those in the reference metric, while for It 2 = O(1), namely, around the bouncing epoch, the ratio defined above becomes R ∼ m 2 M 2 eff / ᾱ 2 M 2 g . This implies that whether the contraction and expansion rates in the physical metric is larger or smaller than those in the reference metric depends on the model parameters.
Furthermore, since the form of the scale factorã(t) in Eq. (8.29) in the physical metric is equivalent to that of a(t) in Eq. (3.1), it is considered that the same consequences as in Sec. III A in terms of the cosmological evolution and values of F ′ (R) [85,86] and F ′′ (R) would be obtained.

IX. CONCLUSIONS
In the present paper, we have reconstructed F (R) gravity models where bounce cosmology can occur. As concrete models, we have demonstrated the cases that in the flat FLRW universe, the scale factor has exponential and powerlaw forms in Eqs. (3.1) and (3.9), respectively. For an exponential form of the scale factor in Eq. (3.1), an F (R) gravity model with the second order polynomial in terms of R is reconstructed, whereas for the power-law form, the resultant F (R) function is proportional to R, equivalent to that in general relativity. In addition, we have investigated the perturbations from the background solutions and examined the explicit stability conditions for these reconstructed models. As a result, it has been found that these models could be stable because the stability conditions can be satisfied. It has to be stressed that the matter bounce scenario [11,19] (for a specific case) proposed by Brandenberger et al. is able to be reproduced also in F (R) gravity.
Also, we have explored a sum of exponentials form of the scale factor in Eq. (5.1) in order to derive an F (R) gravity model in which the bounce in the early universe and the late-time accelerated expansion of the universe can be realized in a unified manner. In this case, a second order polynomial F (R) gravity model is derived as in a model where the scale factor consists of a single exponential term. For this model, we have analyzed the stability condition and confirmed that it can be met. Accordingly, it is considered that the model with the sum of exponentials form of the scale factor could be stable. It is remarkable that the R 2 -gravity theory of the same type as the one realizing inflation occurs as the theory which gives rise to bounce cosmology does.
Furthermore, in the FLRW universe with non-zero spatial curvature, for the scale factor with an exponential form in Eq. (6.6), we have reconstructed a second order polynomial F (R) gravity model and an F (R) gravity model with a term proportional to R and that proportional to 1/R [8]. As a consequence, it has been seen that only in the non-flat FLRW universe with non-zero spatial curvature, a solution can exist, and that if the cosmic curvature vanishes, we can obtain only the de Sitter solution and hence bounce cosmology cannot be realized.
Therefore, when the scale factor is given by an exponential form in Eq. (6.6), by using the reconstruction method, we have derived F (R) gravity models realizing bounce cosmology. Regarding one model leading to bounce cosmology, we have also analyzed the stability conditions and confirmed that these conditions can be satisfied and thus this model can be stable.
Moreover, we have reconstructed an F (R) bigravity model in which bounce cosmology can be realized. It has been verified that in F (R) bigravity, for an exponential form of the scale factor in Eq. (8.29), in the flat FLRW universe bounce cosmology can be realized. It is interesting to emphasize that not only in the physical metric but also in the reference metric the bouncing behavior can happen. Also, if the cosmic time is very far past or future from the bouncing epoch, the contraction and expansion rates in the physical metric are much larger than those in the reference metric. On the other hand, around the bouncing epoch, if the values of the model parameters are determined, we can see which contraction and expansion rates in the physical or reference metric are larger or smaller.