Stability of Singularity-free Cosmological Solutions in Ho\v{r}ava-Lifshitz Gravity

We study stability of singularity-free cosmological solutions with positive cosmological constant based on projectable Ho\v{r}ava-Lifshitz (HL) theory. In HL theory, the isotropic and homogeneous cosmological solutions with bounce can be realized if spacial curvature is non-zero. By performing perturbation analysis around non-flat Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetime, we derive a quadratic action and discuss the stability, i.e, ghost and tachyon-free conditions. Although the squared effective mass of scalar perturbation must be negative in infrared regime, we can avoid tachyon instability by considering strong Hubble friction. Additionally, we estimate the backreaction from the perturbations on background geometry, especially, against anisotropic perturbation in closed FLRW spacetime. It turns out that certain types of bouncing solution may be spoiled even if all perturbation modes are stable.


I. INTRODUCTION
Spacetime singularity at the beginning of Universe is a problem of great importance in standard cosmology. According to singularity theorem proved by Hawking and Penrose [1], a spacetime singularity must be appeared in finite past. Since the appearance of singularity means breakdown of classical gravitational theory, one may expect that fundamental theory beyond General Relativity (GR), i.e., quantum theory of gravity, resolves the problem of infiniteness. A lot of attempts to resolve the singularity at the beginning of Universe have been proposed based on extension of GR [2], e.g., superstring theory [3], loop quantum gravity [4], causal dynamical triangulation [5] and gravity with non-local operator [6]. Regardless those efforts, the dynamics of very early Universe is still unclear. Because we have not achieved complete theory of quantum gravity, yet.
Recently, a gravitational theory attracts attention as a candidate for quantum gravity, which is called Hořava-Lifshitz (HL) theory [7]. The theory is characterized by Lifshitz scaling [8] which is an anisotropic scaling of spacetime : t → b −1 t and x → b −z x with dynamical exponent z. If we set z = 3, all types of ultraviolet divergence via Feynman diagrams can be suppressed in fourdimensional spacetime [9]. It means that gravitational interaction can be renormalized by adding appropriate counterterms. Thus, the renormalizable gravitational action is composed of second order time derivatives and up to sixth order spacial derivatives.
Based on HL theory, ultraviolet spacetime structures have been discussed such as black hole solutions with universal horizon which is a causal boundary for superluminal propagating modes [10]. In particular, singularity avoidance is an intriguing subject for study. We expect that quantum gravitational theory resolves spacetime singularities in classical theory of gravity. The key point to avoid spacetime singularity is violating null energy condition. As indicated in [11], in Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime with non-zero spacial curvature, higher spacial curvature terms in action avoid evolving into singularity. Namely, z = 2 and z = 3 terms mimic "dark radiation" and "dark stiff matter", respectively. Since the energy densities of these effective matter components depend on coupling constants in the theory, null energy condition can be violated if the values of coupling constants are arbitrary. As a result, we can find bouncing universe and oscillating universe as singularity-free solutions [12][13][14].
Although the cosmological singularity avoidance can be realized via higher spacial curvatures as ultraviolet modification of gravity, one may consider such solutions show unstable behavior. More specifically, it may possible that the effective matter components derived from z > 1 Lifshitz scaling terms make spacetime unstable because of violation of the energy condition. It is reasonable that the spacetime around Planck scale is perturbed by quantum fluctuation of gravity. Thus, to examine the spacetime stabilities of these singularity-free solutions is indispensable in order to construct a cosmological scenario without initial singularity.
In the paper [14], dynamics of Bianchi type IX spacetime, i.e., spacially homogeneous closed cosmological model, is discussed. In other words, the effect of spacetime anisotropy to singularity-free solutions in closed FLRW spacetime is examined. From the result, it is found that the stability against small anisotropic perturbation depends on the coupling constants of the theory. Thus, we expect that the dynamics of the other types of perturbations are also affected by the coupling constants. Since HL theory is renormalizable, the values of coupling constants at Planck scale can be evaluated via beta functions from renormalization group, in principle. If we obtain the values of coupling constants in ultraviolet regime, we may predict the dynamics of our Universe at very early stage, which cannot directly observed. Thus, we set the goal of this paper to show the stability conditions for singularity-free cosmological solutions against linear order perturbation not only anisotropic modes but also inhomogeneity ones.
The perturbation analysis in FLRW spacetime in the context of HL theory had been discussed. Without spacial curvature, there is quite a lot of study, e.g., primordial perturbation [15], stability of scalar perturbation [16,17] and stability of de Sitter spacetime [18]. Turning our attention to the case including non-zero spacial curvature, the analyses of scalar perturbation have been performed. In [19], the authors show scalar perturbation in vacuum FLRW spacetime, and in [20], dynamics of scalar field in bouncing universe is discussed. Recalling our motivation to investigate spacetime stabilities of singularity-free solutions, it is necessary to see dynamics of tensor and vector degree of freedoms as well as scalar ones.
Thus, in this paper, we perform perturbation analysis regarding tensor, vector and scalar modes based on projectable HL theory. Although non-projectable HL theory which is infrared completion of projectable HL theory has been proposed [17], its renormalizablity is still unclear. Therefore we focus only on the projectable case because of the renormalizable characteristic.
The rest of this paper is organized as follows : In Section II, we briefly review projectable HL theory, especially in FLRW background with non-zero spacial curvature. The perturbation theory around non-flat FLRW spacetime in HL theory is discussed in Section III. In Section IV, we discuss stabilities of bouncing solution of HL theory in non-flat FLRW spacetime by investigating ghost and tachyon-free conditions. Additionally, we estimate the backreaction from the perturbation on background geometry, especially, against anisotropic perturbation in closed FLRW spacetime. Section V is devoted to conclusion of this paper.

II. HOŘAVA-LIFSHITZ THEORY IN FLRW SPACETIME
We briefly review the projectable HL theory, especially, in FLRW spacetime. The gravitational action is given by [21] 2) where, m LV is Lorentz violating mass scale which may be expected Planck mass. The extrinsic curvature K ij is defined in terms of the lapse function N , the shift vector N i and the three-dimensional induced metric g ij : where, ∇ i represents the three-dimensional covariant derivative. The potential terms are defined by where, Λ is the cosmological constant, λ and g n (n = 1-8) are dimensionless coupling constants. The potential terms include the higher order spacial curvatures R ij , R := R i i up to sixth order spacial derivatives. In what follows, we adopt the unit m LV = 1 unless otherwise noted.
We shall focus on the non-flat FLRW spacetime whose induced metric g ij is given by where, a is a scale factor which depends only on time.
The metrics whose spacial curvature K = 1 and −1 correspond to closed and open FLRW spacetime, respectively. For closed case, the domains of the variables are defined by 0 ≤ χ < π, 0 ≤ θ < π and 0 ≤ φ < 2π. For open case, 0 ≤ χ < ∞, the domains of θ and φ are as same as closed ones. Assuming above ansatz and taking variation with respect to a, N and N i , we obtain the dynamical equation of scale factor, the Hamiltonian constraint and the momentum constraint, respectively. Since the momentum constraint gives a trivial relation, we omit it. The dynamical equation of the scale factor is given by (2.9) g s := 12K 3 (9g 4 + 3g 5 + g 6 ) . (2.10) After taking variation, we have imposed the gauge condition so that N = 1 and N i = 0. It is notable that the terms derived from the higher spacial curvatures behave as virtual matter fields. More specifically, the g r and g s terms effectively work as "radiation" and "stiff matter", respectively. Note that the g 7 and g 8 terms do not affect to the background solution because of the spacetime symmetry. Turning our attention to the Friedmann equation which corresponds to the Hamiltonian constraint in FLRW spacetime. In our case, we cannot construct the Friedmann equation via taking variation of the action. Since we have imposed the projectability condition, i.e., N (t, x) → N (t), the Hamiltonian constraint takes the following form : namely, the Hamiltonian constraint turns to be a global condition instead of local one. It means that we have to know the information within the entire spacetime to construct the Hamiltonian constraint, and thus, the equation (2.11) is not viable without special assumption. Then, we derive the Friedmann-like equation by considering the structure of the basic equations in FLRW spacetime. Note that the Friedmann equation (and matter conservation law) basically generates the dynamical equation of the scale factor by differentiating with respect to time. Thus, we can obtain the following equation by performing the time integration of (2.8) : where, C is an integration constant. Note that the C term behaves as a dust whose energy density is proportional to spacial volume, i.e., a −3 [22]. For later convenience, we define the following quantities : (2.14) Then, the dynamical equation of the scale factor and Friedmann-like equation can be written by E = 0 and H = 6(3λ − 1)C/a 3 , respectively. Since C does not depend on time, the value of H is determined by the initial condition.
One may notice that the spacetime dynamics for λ > 1/3 and λ < 1/3 are completely different, i.e., the sign of time derivative terms are flipped. Since the limit to GR can be obtained by taking λ → 1, we exclude the λ ≤ 1/3 case, in what follows. Additionally, the value of g 1 must be negative to recover the result based on GR at least at a background level. Then, we set g 1 = −1 in the rest part of this paper by performing a suitable rescaling of time.

III. PERTURBATION ANALYSIS AROUND NON-FLAT FLRW BACKGROUND
Since the ADM formalism is employed in this paper, we define the perturbed quantities of ADM variables δN , δN i and δg ij as follows : where,N ,N i andḡ ij denote the background lapse function, shift vector and three-dimensional induced metric, respectively. Since we consider the quadratic gravitational action, we define δN , δN i and δg ij as follows : where, α, β i and h ij are perturbations of first order. Furthermore, we define the perturbed quantities with upper indices as h i j :=ḡ ia h aj , h ij :=ḡ iaḡjb h ab , h :=ḡ ab h ab and β i :=ḡ ia β a . Then, we shall decompose α, β i and h ij into the scalar, vector, tensor modes.

A. spherical and pseudo-spherical harmonics
In non-flat FLRW background, the scalar, vector and tensor perturbations can be decoupled by employing spherical or pseudo-spherical harmonics [23][24][25]. More specifically, (χ, θ, φ) dependences of perturbed ADM variables can be expanded by each modes of harmonics, which is similar to black hole perturbation theory. For example, a scalar function Θ can be expanded by scalar spherical harmonics Q (n;lm) in closed FLRW spacetime : where, Θ (n;lm) is a coefficient of each (n; l, m) modes which depend only on time. We summarize the definitions of the tensor spherical and pseudo-spherical harmonics in Appendix A.
Since we consider the four-dimensional spacetime, ten types of independent tensor harmonics must be equipped as a basis set. In this paper, we employ one of possible orthonormal basis set Y as follows : When the perturbation in closed (open) FLRW spacetime is considered, we refer the quantities with hatted (checked) superscript. In what follows, we abbreviate these superscripts to unify the discussions 1 . Turning our attention to the decomposition of ADM variables by harmonics. The scalar perturbation can be expanded by the scalar harmonics : Note that we do not have to expand the perturbed lapse function in terms of harmonic function. Since projectability condition is imposed, the lapse perturbation also depend only on time variable. The vector perturbation can be expanded by the vector harmonics : Note that the lapse and shift perturbation do not exist in tensor perturbation. In this paper, we eliminate the following perturbations by choosing gauge : In Appendix C, the gauge structure in FLRW spacetime is summarized.

B. quadratic action
Turning our attention to the perturbed action at quadratic order by employing (pseudo-)spherical harmonics. Note that the formulae of harmonic functions we have applied in this part are summarized in Appendix B.
1 For simplicity, we denote to sum up n modes. To be precisely, it must be replaced by integration symbol when we consider pseudo-spherical ones because n turns to be continuous number.
Before performing perturbation analysis, we shall mention n = 1 case for closed FLRW. In this case, the perturbation is given by (3.14) One notice that this mode corresponds to just a shift of the scale factor, i.e., a(t) → a(t) + δa(t), and thus, we exclude this perturbation mode in what follows. We firstly consider the perturbation of kinetic terms to clarify the dynamical degree of freedom. The quadratic kinetic actions for scalar, vector and tensor perturbation are given by where, ν 2 is a eigenvalue of the harmonics which is defined in terms of the perturbation mode n ≥ 1 : Since the perturbations with different degrees do not mix, we have abbreviated the superscript (n; lm). Note that β (S) and h (G) are vanished when l = 0 and l ≤ 1, respectively.
Taking variation with respect to β (Q) and β (S) , we obtain the following constraint equations.
and plugging these relations into the actions, we obtain the simplified quadratic action as follows : 21) and the vector mode is vanished. Summing up the quadratic potential terms, we finally obtain where, F (G) and G (G) can be regarded as kinetic term and mass term of tensor perturbation, respectively, which are given by where, g 56 := g 5 + g 6 . On the other hand, F (Q) and G (Q) are regarded as kinetic term and mass term of scalar perturbation : The terms including Hubble parameter H and its time derivativeḢ have been eliminated by applying the background equation of motion E = 0. Note that there does not exist the integration constant C appeared in the Friedmann-like equation H = 6(3λ−1)C/a 3 . It is because the Hamiltonian constant basically arises as a coefficient of the lapse perturbation. In our case, we have fixed the gauge so that α = 0, then, there is no ambiguity in quadratic action. Additionally, one notice that there are no dynamical scalar modes in the closed FLRW spacetime with n = 2 because F (Q) is automatically vanished. Particular attention should be given to the fact the quadratic actions (3.22) and (3.23) cannot be reduced into the result based on GR even if we take the limit λ → 1 and the terms from the higher spacial curvature are neglected. Based on GR, the scalar degree of freedom is absent if we consider the vacuum FLRW spacetime. However, in our case, the scalar perturbations cannot be eliminated when we take such a limit. This fact is due to the gauge structure shown in Appendix C. More specifically, the time derivatives of scalar perturbatioṅ h (Q) must be appeared because of the Lorentz violation, i.e., the rotation of time direction is not allowed.
In order to clarify the stability conditions for the perturbations, we firstly focus on the coefficient of kinetic terms and mass terms.

ghost avoidance
We shall concentrate on the terms derived from L K . Since F (G) and F (Q) are the coefficients of the kinetic terms of each perturbation modes, the conditions F (G) ≥ 0 and F (Q) ≥ 0 are required to avoid ghost instabilities in tensor and scalar perturbation, respectively. Since F (G) = 1, the tensor perturbations do not show ghost instability for any choice of the coupling constants.
On the other hand, the condition for the scalar one is not trivial. In closed FLRW spacetime, namely K = 1, to satisfy F (Q) ≥ 0 for every n ≥ 2, we find λ ≥ 1 . (3.28) One may notice that (3.28) is almost same as the stability conditions in Minkowski spacetime [13]. In open FLRW spacetime, namely K = −1, the stability condition for every n ≥ 1 mode is given by which gives tighter condition than that of closed case.

tachyon avoidance
The terms derived from L P , namely, G (G) and G (Q) can be regarded as squared masses of tensor and scalar perturbations, respectively. Therefore, G (G) ≥ 0 and G (Q) ≥ 0 must be satisfied, otherwise tachyon instability appears. One can see that equalities G (G) = 0 and G (Q) = 0 give quadratic equations with respect to a 2 and ν 2 with the coefficients related with g i . It means that the ranges for the scale factor in which G (G) ≥ 0 and G (Q) ≥ 0 can be expressed in terms of the coupling constants, in principle.
In this part, we focus on the infrared stability, i.e., the case without effects of higher order spacial curvatures is considered. Then, the mass terms in the quadratic action are reduced into We see the tensor perturbation mode has positive squared mass, while, the scalar one shows opposite sign. Thus, in infrared regime, the negative squared mass of the scalar perturbations cannot be avoided. The similar situation has been found in perturbation around Minkowski spacetime [16,17]. Of course, the negative squared mass does not always lead the tachyon instability. If the growing time scale for the scalar perturbation is sufficiently small relative to cosmological time scale, the instability is suppressed. Another possibility to avoid infrared tachyon instabilities is to consider the extended version of HL theory, namely, relaxing projectability condition [17]. It is known that additional (∇ i ln N ) 2 term in extended action can stabilize the scalar perturbation at least in flat background.

IV. STABILITY ANALYSIS OF SINGULARITY-FREE SOLUTIONS
In this section, we analyze the stabilities of singularityfree solutions in non-flat FLRW spacetime. Since the typical scale of the singularity avoidance is expected to be Planck scale, we focus only on the solutions which potentially connect to macroscopic universe. More specifically, we consider bouncing cosmological solutions with positive cosmological constant, which show bouncing behavior at a = a T > 0 and turn to accelerating expanding phase.
One may consider the bouncing solution without cosmological constant in open FLRW spacetime can also evolve to macroscopic universe with asymptotic Milne expansion. However, such solutions seem to be unstable because of weak Hubble friction during expanding phase. As we mentioned in previous, scaler perturbation possesses negative G G in infrared regime, and thus, we exclude the case with non-positive cosmological constant.

A. background solutions
The classification of the solutions in vacuum FLRW background is performed in the paper [13]. In that paper, it is found that there are two types of singularity-free solutions. One is bouncing universe, that is, the initial contracting universe turns into the expanding phase at a = a T , and the universe keeps expansion without finite upper bound of the scale factor. The other is oscillating universe whose scale factor is bounded in the range of 0 < a min ≤ a ≤ a max < ∞. Then, the universe shows periodic oscillatory behavior without singularity.
The dynamics of the background spacetime can be examined via rewritten Friedmann-like equation : For simplicity, we have taken the integration constant C = 0. Since the first term of the left hand side of (4.1) is not negative, the possible ranges of the scale factor are where U ≤ 0. Note that g r and g s are related with the coupling constants in L P , and then, these values can take both plus and minus sign if the coupling constants are arbitrary. Therefore, we can consider the situation in which the scale factor is bounded below by some non-zero minimum value. It means the universe is forbade to fall down into the singularity. We would like to stress that the singularity avoidance is induced because the energy condition is effectively violated.
To analyze the background dynamics with Λ > 0, it is convenient to rewrite the potential U by rescaled variables with respect to ℓ := 3/Λ: where,ã := a/ℓ,g r := g r /ℓ 2 andg s := g s /ℓ 4 . Sincẽ U = 0 is essentially a cubic equation ofã 2 , we obtain three analytic solutions as follows : where, I = 1, 2, 3. Note that the points at whichŨ(ã) = 0 can be found when correspondingã [K] I takes real and positive values.
Additionally, the above roots cannot be applied to the special case withg s = 0. In this instance, we obtain two analytic solutions as follows : 1. closed FLRW (K = 1) In closed FLRW universe, namely, K = 1 case, three types of singularity-free solutions are found. We show the typical potentialsŨ for these solutions in FIG. 1. (a) B [1] BC : A universe which shows bouncing behavior for initial scale factorã ini ≥ã T , however evolves into big crunch forã ini ≤ã BC . We classify this type of the solutions as B [1] BC . The typical potential  BC (red dashed curve), B [1] (blue solid curve) and B [1] O (green dotted curve). The coupling constants are chosen asgr = 3/10 andgs = 2/10 for B [1] BC ,gr = 17/20 andgs = −1/8 for B [1] ,gr = 17/20 andgs = −7/100 for B [1] O . The bouncing radii and the maximum radii of big crunch solution are denoted byãT andãBC , respectively. Additionally, the maximum and minimum radii of oscillation are given byãmax andãmin.
is given by the dashed red curve in FIG. 1. Note that the domainã BC <ã <ã T is forbidden. The solutions B [1] BC can be found in the following two cases : (b) B [1] : A universe which bounce atã =ã T without big-bang singularity for any possible initial scale factorã ini ≥ã T . We classify this type of the solutions as B [1] . The typical potential is given by the solid blue curve in FIG. 1. The solutions B [1] can be found in the following three cases : (c) B [1] O : A universe which shows bouncing behavior for initial scale factorã ini ≥ã T , on the other hand, oscillates ifã min ≤ã ini ≤ã max . The other choice of the initial scale factor is forbidden. We classify this type of the solutions as B [1] O . The typical potential is given by the dotted green curve in FIG. 1. The solutions B [1] O can be found if the following condition is satisfied : In open FLRW, namely, K = −1 case, two types of singularity-free solutions are found. We show the typical potentialsŨ for these solutions in BC , gr = 4/5 andgs = −3/20 for B [1] . The bouncing radii and the maximum radii of big crunch solution are denoted byãT andãBC , respectively. BC in closed case. Namely, this type of the solutions shows bouncing behavior for initial scale factorã ini ≥ã T , on the other hand, evolves into big-bang singularity forã ini ≤ã BC . We classify this type of the solutions as B BC can be found if the following condition is satisfied : : As is the case with B [1] in closed FLRW, this type of the solutions also shows bouncing behavior for any possible initial scale factorã ini ≥ã T . We classify this type of the solutions as B [−1] . The typical potential is given by the solid blue curve in 1 ≤ã ≤ã [1] 2 ,ã ≥ã  BC can be found in the following two cases : (ii)g s = 0 ,g r < 0 withã T =ã We show the properties of singularity-free solutions in TABLE I and the distribution of the singularity-free solutions in (g r ,g s ) plane in FIG. 3.

B. perturbation analysis
Turning our attention to the ultraviolet stability. Since the conditions for avoiding ghost instabilities have been discussed in previous section, we concentrate on the positivities of G (Q) and G (G) including higher order spacial curvatures. Our aim is to clarify the stability conditions for all perturbation modes with finite values of the coupling constants λ and g i (i = 1-8). Otherwise the asymptotic safety is violated, which may causes the divergence of gravitational force in ultraviolet regime.

tensor perturbations
The stability condition for the tensor perturbation is given by G (G) ≥ 0. The positivity of F (G) is automatically satisfied, namely, there is no ghost tensor mode for any choice of coupling constants. Since G (G) = 0 gives a quadratic equation with respect to a 2 and ν 2 , we must require g 8 ≥ 0. Otherwise, tensor perturbations with large ν 2 show unstable behavior.
To illustrate the stability of each perturbation mode, we firstly analyze the lowest order of the tensor perturbation in closed FLRW background, namely, K = 1 and n = 3 (ν 2 = 8) case. Then, we find In order to stabilize this perturbation mode for any a > 0, the following condition must be satisfied : for g r < 12g 3 . (4.9) One may notice that the stability condition (4.9) reproduces the result shown in [14]. In fact, the tensor perturbations with n = 3 in closed FLRW background include homogeneous and anisotropic perturbation, namely, Bianchi type IX spacetime with small anisotropy (see Section IV D).
It is obvious that we should impose the positivity of G (G) for any perturbation mode ν 2 to ensure the stability of spacetime against tensor perturbation. Of course, we can express the stability conditions for tensor perturbation for any viable a > 0 and perturbation mode ν 2 in terms of the coupling constant g i . Since G (G) = 0 gives a quadratic equation in terms of a 2 and ν 2 , we can solve it, in principle. However, it is found that the explicit form is quite complicated. Especially, in closed background, ν 2 takes discrete value and this fact complicates the analysis. Thus, instead of considering general case, we show a special case, namely, g 3 = 0 in open FLRW spacetime.
Then, we find that every mode of tensor perturbation can be stabilized for any a ≥ 0, if one of the following four condition is satisfied : Note that there is some difficulty in reconciling above conditions with the conditions for bouncing solutions in open FLRW spacetime. Referring TABLE I  (ii) Since g s is constrained to be positive number under the condition (ii), there is no stable bouncing solution B [−1] . Additionally, one may notice that the lower bound of g s is given in terms of g 2 r . Since |g r | ≫ |g s | should be satisfied to stabilize the solution B

[−1]
BC , the condition (ii) is hard to compatible unless small Λ > 0 is set. As we shown later, small Λ > 0 is not preferable to stabilize scalar perturbation at infrared regime. . Conversely, in closed FLRW spacetime, the tensor perturbations can be stabilized without special tuning. Naively, g r < 0 and g s < 0 are preferred to satisfy the positivity of G (G) for any a > 0 and ν 2 . Referring TABLE I and FIG. 3, we find B [1] (ii) solution is under such a condition.

scalar perturbation
Although the negative sign of G (Q) in infrared regime cannot be avoided, G (Q) may take positive value in the deep ultraviolet region in which the effects of higher spacial curvatures are predominant. To realize stable bouncing phase, we require the positivity of G (Q) , at least in the range [0, a ini ] with a ini > a T . Since G (Q) = 0 gives a quadratic equation in terms of a 2 with upward convex, the positivity of G (Q) at least in [0, a ini ] is guaranteed if both of the following conditions are satisfied : Note that the case with g 8 < 8g 7 /3 must be excluded, otherwise the scalar perturbations with large ν 2 shows unstable behavior.
In our analysis, we can find bouncing solutions B  One may wonder about the appropriate value of a ini . The initial scale factor a ini seems to relate with a quantum creation of the universe. Therefore, it is natural to consider the typical energy scale is estimated at Planck scale. On the other hand, the typical scale of the bouncing radius is also expected to be around at Planck scale. Because, the bouncing behavior is induced by higher spacial curvature terms, namely, quantum gravitational corrections. More specifically, we assume the three-Ricci curvature represents the energy scale, namely, m ∼ √ R ∝ a −1 . Then, the ratio of quantum creation scale m ini to bouncing scale m T is given by It is natural to consider the ratio is order one. Additionally, we consider the upper limit of a ini to satisfy the positivity of G (Q) during bouncing phase. As we noted, G (Q) must be negative for large a. Thus, there exists a critical value of the scale factor a crit . Namely, any scalar perturbation mode possess positive squared mass for a > a crit , however, any one of scalar perturbation mode turns to be zero at a = a crit . Obviously, a ini must be in (a T , a crit ). Further constraint for a ini can be imposed by considering the dynamics of the perturbations.

C. dynamics of perturbations
To construct a scenario for the non-singular cosmological evolution, we have to pay attention to the dynamics of the perturbations. Taking variation of the quadratic actions with respect to h (G) and h (Q) , we obtain the equations of motion for tensor and scalar perturbations, respectively :ḧ where, we define effective squared masses of the tensor and scalar perturbation as We firstly consider the contracting phase before bounce. In this era, the perturbations feel a Hubble acceleration which is derived from the second terms in (4.13) and (4.14) because of negative Hubble parameter H < 0. Since the Hubble acceleration enhances the both perturbation modes, this effect should be suppressed by effective mass terms.
Intuitively, the magnitudes of the Hubble accelerations for tensor and scalar perturbation are given by H 2 . Thus, to suppress the unstable behavior, we require (4.16) throughout contracting phase. This condition gives a further constraint on possible value of the initial scale factor. Namely, a ini should be in the range of (a T , a H ), where, a H is a value of scale factor in which any one of effective mass of perturbation turns to be M 2 = H 2 . Namely, for a > a H , every perturbation modes shows positive effective squared mass which is larger than H 2 . Then, the condition (4.16) is ensured in contracting era for a T < a ini < a H . Note that the possible range for the initial scale factor can be broadened by tuning the value of coupling constant λ. Namely, the dependence of λ in effective mass-Hubble parameter ratios can be evaluated as follows : then, one can see that large value of λ weakens the effect of Hubble acceleration in both cases.
After the bounce at a T , the universe turns to expand. As we mentioned, the effective squared mass of the scalar perturbation must be negative in infrared regime. Thus, to stabilize the perturbation, Hubble friction with H > 0 must overcome the effects of the negative squared mass of scalar modes. Namely, we require  One may wonder about the growth of scalar perturbation during tachyon instability. Seeing equation of motion for scalar perturbation, it is natural to speculate that the growth rate is related with the minimum value of M 2 (Q) /H 2 . Thus, we firstly clarify the minimum value of the squared effective mass M 2 (Q) . For simplicity, we consider the asymptotic region, i.e., for large perturbation mode ν 2 . In this limit, the kinetic and mass terms of scalar perturbations given by (3.26) and (3.27) are reduced into the following forms : . Then, the minimum value of M 2 (Q) is given by with η := 4g r + 3g 3 K 2 After taking above value, M 2 (Q) monotonically increases with time and approaches to zero. Then, the scalar perturbation is stabilized by Hubble friction. The important point is that the minimum value does not depend on the perturbation mode ν 2 in this limit. Thus, we can conclude that the minimum value of M 2 (Q) can be bounded in finite value. Namely, if the ratio to Hubble parameter for large ν 2 , is sufficiently suppressed, it is expected that there is no serious instability at least at the classical level. Note that the large value of the positive cosmological constant Λ and/or the small value of λ decreases above value.
In above discussion, we limited our analysis to the case with large ν 2 . However, we can investigate the case with intermediate value of ν 2 in the same manner, and we also find the minimum value of M 2 (Q) /H 2 is also affected by the values of Λ and λ. Namely, large Λ > 0 and small λ are preferred.
Then, we demonstrate the growth of the scalar perturbation without condition  Furthermore, we mention the relation between the scalar growth rate and the minimum value (4.22). To evaluate the relation, we perform numerical calculation by setting various λ with fixed perturbation mode. Then, the evolutions of the scalar perturbation is shown in FIG.  7. Additionally, the detailed data of the asymptotic values of r (Q) in terms of λ are shown in TABLE III. From this result, we can conclude that the large value of M 2 (Q) /H 2 with M 2 (Q) < 0 enhances the growth of scalar perturbation. In other words, the scalar growth rate is amplified by choosing large value of λ and small value of Λ > 0.  When the accelerating expansion caused by a cosmological constant persists, the effect of the spacial curvature K turns to be irrelevant. Then, the analysis can be simplified into the case with K = 0 and Λ > 0. In other words, the spacetime can be approximated as the de Sitter solution at the late time of the evolution after bounce (c.f. cosmic no-hair theorem [26]). It should be noted that the detailed analysis of de Sitter spacetime stability has been already performed in the papers [18]. In those paper, the authors indicated that the instability of scalar perturbation may be cured if the background spacetime is the de Sitter solution. It is worth mention-ing that the quadratic actions in flat FLRW spacetime can be reproduced by simply taking a limit K → 0 in (3.22) and (3.23) 2 .

D. backreaction of the perturbation on background geometry
We further discuss the stability of bouncing solution by considering a backreaction of the perturbation. Especially, against anisotropic and homogeneous perturbation in closed FLRW spacetime. Such perturbation modes can be derived by considering Bianchi type IX spacetime whose three-dimensional space is homogeneous, however, isotropy is not always hold. The metric is given by where, a represents the scale factor, ω i (i = 1, 2, 3) is an invariant basis which is given by The traceless symmetric tensor β ij represents anisotropy. Since Bianchi type IX space belongs to Bianchi class A, β ij can be diagonalized without loss of generality as follows [31] : When β ± = 0, the spacial isotopy is restored, namely, closed FLRW spacetime. The basic equations are given byḢ where, V (a, β ± ) := −a 3 L P /128 is a potential which is given by the spacial Ricci curvature terms (in [14], the explicit form is shown). The equation (4.26) and (4.27) correspond to the dynamical equations of the scale factor a and the anisotropy β ± , respectively. As is the case with the case of FLRW spacetime, we obtain a Friedmannlike equation (4.28). Due to integration with respect to time variable, a integration constant C IX appears. For simplicity, we consider the case with C IX = 0. Assuming |β ± | ≪ 1, the potential V is reduced into where, U 0 (a) and U 2 (a) are defined by It is worth mentioning that U 0 and U 2 are related to the potential in FLRW spacetime (4.2) and the squared effective mass M 2 := G/F of n = 3 tensor mode as follows :  One may notice that the equation of motion for the tensor perturbation with n = 3 is reproduced if β ± is replaced into h (G) in (4.27) (see (4.8) and (4.13)). In other words, this perturbation modes include the homogeneous and anisotropic perturbation in closed FLRW spacetime. Then, the equation (4.28) can be rewritten as follows : where, Since we impose the positivity of the tensor squared mass M 2 (G) to stabilize the perturbation, E β± always takes positive value. From (4.34), one can see that the scale factor is regarded as a particle with energy 2a 2 (E β+ + E β− )/(3λ − 1) in potential U(a). Namely, the possible range for the scale factor is broaden due to an anisotropic energy E β± .
Then, we examine the backreaction on the singularityfree solutions, especially B [1] BC . This type of bouncing solutions realizes singularity avoidance due to the potential barrier U ≥ 0 between a BC and a T (see FIG.1). However, the cosmological bounce at a T may be spoiled if the backreaction from anisotropic perturbation is considered. Namely, the energy for the scale factor is lifted up to 2a 2 (E β+ + E β− )/(3λ − 1) due to anisotropic effect, and then, the potential barrier can be overleaped if the anisotropic energy exceeds the local maximum value of the potential U. In FIG. 8, we show the typical example based on the solution (i) listed in TABLE II. In this analysis, we set λ = 1 and the initial conditions are given by a ini = 1.5a T ≈ 2.784, β ± = β ini andβ ± = 0. In the top figure, the evolutions of the anisotropic energies are shown. The anisotropic energy 2a 2 (E β+ + E β− )/(3λ − 1) takes maximum value at the bouncing time t = t T ≈ 1.261. If the initial anisotropy exceeds a critical value β crit ≈ 0.0569, the universe results in big crunch by overleaping the potential barrier (the purple curves in the top figure of FIG. 8).
We further mention the dynamics of the anisotropy β ± andβ ± (the bottom figure of FIG. 8). One may notice that the oscillating amplitudes of β ± are almost invariant throughout bounce, while those ofβ ± are enhanced whose maximum amplitudes reach up to 10 −1 order. It is not quite unnatural because the dynamics of anisotropic perturbation is approximately governed by the following equation :β Around the bouncing point, the term including Hubble parameter can be ignored. Then, the oscillating frequency is naively given by ω ≈ M (G) = G (G) /F (G) . As a result, the amplitudes ofβ ± is approximately estimated by |β ± | ≈ M (G) |β ± |. Since M 2 (G) ∝ a −6 for small scale factor, the large anisotropic energy at bouncing point is induced if we consider small bouncing radii. Thus, we conclude that the backreaction on the bouncing universe with large bouncing radii tend to be small.
The backreactions from the other perturbation mode are unclear. However, it is natural to consider that the most symmetric spacetime corresponds to the lowest energy state. Then, one may speculate that the other perturbation modes also lift up the energy for the scale factor like the case of the anisotropic perturbation as (4.34) and (4.35).
We also would like to point out that the oscillating universe obtained by B [1] O with a min ≤ a ini ≤ a min can possibly evolve into macroscopic universe. In this type of potential U(a), the oscillating and bouncing solutions are separated by a potential barrier between a max and a T (see FIG.1). However, if the potential barrier is sufficiently small, it is possible that the initial oscillating era shifts into accelerating expanding phase via energy induced by perturbations. In fact, the similar evolution of the universe is found when the spacetime anisotropy is large [14].

V. CONCLUSION
To avoid big-bang singularity at the beginning of universe, it is essential to consider the situation which null energy condition is violated, i.e., p + ρ = (1 + w)ρ < 0. Based on HL theory, the higher spacial curvatures in action possibly behave as such exotic matters. Thus, one find singularity-free cosmological solutions, such as bouncing universe in non-flat FLRW spacetime. However, it is natural to consider the effective exotic matters which violate null energy condition destabilize spacetime.
In this paper, we investigate the stabilities of bouncing solutions via perturbation analysis around non-flat FLRW spacetime. Employing (pseudo-)spherical harmonic functions, both tensor and scalar perturbations can be decomposed into each (n; l, m) modes. Then, perturbed actions for tensor and scalar modes at quadratic order are reduced into (3.22) and (3.23), respectively. Note that the integration constant C induced by the lack of local Hamiltonian constraint does not affect to the quadratic action, however, background dynamics is influenced, i.e., dust-like additional term is joined in the Friedmann-like equations. In our analysis, the integration constant is set to zero, for simplicity. Thus, the result may be slightly changed if we consider non-zero C, i.e., the spacetime stability around bouncing point.
In order to avoid ghost instabilities, we must require the coefficients of kinetic terms in quadratic action to be positive, namely, F (G) ≥ 0 and F (Q) ≥ 0 for any a > 0 and perturbation mode ν 2 . Since F (G) = 1, tensor perturbations do not show ghost instability. On the other hand, the condition for ghost avoidance in scalar perturbation is expressed in terms of λ : Note that in flat FLRW spacetime or Minkowski spacetime, the stability condition for scalar perturbation is given by λ > 1 and scalar degree of freedom vanishes when λ = 1. It is known that there is no smooth connection between λ > 1 and λ = 1 because of strong coupling problem. However, in closed FLRW case, we can take smooth limit λ → +1 without any singular behavior, and then, the scalar perturbation can propagate even if the case with λ = 1 is considered. Thus, the limit λ → 1 does not mean GR is restored. The dissimilarity is due to the gauge structure mentioned in Appendix C.
We further consider the positivities of G (G) and G (Q) for any perturbation mode (n; l, m). In order to stabilize the perturbation modes with large n, the following conditions must be satisfied : Although it is possible that G (G) ≥ 0 for any a > 0 and viable ν 2 , the negativity of G (Q) cannot be avoided in infrared regime, i.e., there must be exist a crit at which any one of scalar perturbation mode turns from G (Q) > 0 to G (Q) = 0. This result is consistent with infrared instability of scalar graviton in flat background. Note that the negative value of G (Q) does not always mean the instability of the scalar perturbation. Then we have investigated the dynamics of perturbations in bouncing universe via equations of motion. In contracting phase, perturbations are possibly amplified by the negative sign of Hubble term. To suppress the instabilities, the following condition must be satisfied : In infrared regime, the squared effective masses of scalar perturbations must be negative for any choice of the coupling constants g i . Thus, we impose the following condition in order to overcome the effect of M 2 (Q) < 0 : Thus, the stable bouncing solutions are limited to the case with Λ > 0. If the condition (5.4) is violated, scalar perturbation is amplified, which means the tachyon instability is occurred. Since M 2 (Q) decreases as a −2 at infrared regime, the period for the scalar instability is temporal. The growth rates of the scalar perturbations are related with the minimum value of M 2 (Q) /H 2 . Then, small λ − 1 > 0 and/or large Λ > 0 are preferred to suppress the growth of scalar perturbations. It may be interesting to estimate permissible growth rate of scalar perturbation by referring observational cosmological data. Then, one can derive further constraints to the values of coupling constants.
Additionally, we have investigated backreaction from perturbation on background geometry, especially, against anisotropic perturbation in B [1] BC type solution. Considering Bianchi type IX spacetime with small anisotropy, the modified Friedmann-like equation including backreaction from the anisotropic perturbation can be derived as (4.34). It is found that the energy for the scale factor is lifted up by anisotropic perturbation. Thus, the universe can evolve into the singularity if the potential barrier between a BC and a T is sufficiently small. We also pointed out that the anisotropic energy tend to be large if the bouncing radius is small, because the oscillating amplitudes ofβ ± are enhanced.
We would like to stress that the bouncing solutions in open FLRW spacetime tend to be unstable for the following reasons : (i) the stability conditions for tensor perturbation basically contradicts to bouncing conditions. Intuitively, g r > 0 and g s > 0 are preferred to satisfy G (G) ≥ 0 for any a > 0 and viable ν 2 . However, we cannot find any types of bouncing solution in open FLRW spacetime under such a condition. Thus, a certain level of tuning of the other coupling constants is required. (ii) in infrared regime, the Hubble friction which is significant to suppress tachyon instabilities in scalar perturbation tend to be weak. In closed FLRW spacetime, the minimum value of M 2 (Q) /H 2 can be converged to zero if λ approaches to unity via renormalization group flow.
However, in open FLRW spacetime, we have to constrain λ > 2 to avoid ghost instabilities.
Our conclusion is that we have shown that non-singular cosmological solutions in non-flat FLRW spacetime can be stable against tensor and scalar perturbation, at least at the linear level. Since projectable HL theory is proved to be truly renormalizable in perturbation approach, we can calculate the values of the coupling constants via beta functions from renormalization group, in principle. Thus, it may be possible that the beginning of our Universe can be predicted based on well-known perturbative quantization approach.
In our case, i.e., HL theory under projectability condition, it is indispensable to consider accelerating expanding phase after bounce in order to suppress the effect of negative squared mass of scalar perturbation in infrared regime. We speculate that the scalar instabilities in infrared regime is the nature of projectable HL gravity theory. This infrared pathological behavior is conceiv-ably resolved by considering extended theory, i.e., nonprojectable HL gravity whose scalar graviton can be stable at least in Minkowski spacetime. We additionally mention that infrared limit of nonprojectable HL theory (i.e., without higher spacial curvatures) is included within the framework of Horndeski theory which is the general theory of ghost-free scalar-tensor gravity [27]. Based on Horndeski theory, it turns out that any non-singular cosmological solution is unstable [28]. On the other hand, it is worth mentioning that the nogo theorem for the stable non-singular solution can be violated if the extended theory including higher spacial curvatures is considered [29]. Thus, in view of the situation, it should be interesting to investigate the stability of bouncing solutions with higher spacial curvatures based on non-projectable HL theory as a special case of extended Horndeski scalar-tensor gravity [30].
The vector spherical harmonic functions are classified into two classes. Since the vector quantity can be decomposed into gradient part and rotational part, we define the gradient of the scalar harmonics ψ (lm) A and the dual of the gradient φ where, D A denotes a covariant derivative on two sphere, ǫ AB is Levi-Civita tensor on two-sphere. Note that ψ A and φ A posses even and odd parity, respectively. The tensoral ones are classified into three-types : η AB is a trace part, which is proportional to two-metric. ψ AB and φ AB are traceless with even and odd parity, respectively. The explicit forms are given by (A8)
In what follows, we abbreviate the superscriptsˆandˇif not otherwise specified.

a. scalar type
We introduce the scalar type harmonics which contribute to the scalar perturbation. Since the scalar quantities have already been introduced in previous part, we focus only on the vector and tensor quantities.
The vector quantities C i are defined by C (n;lm) i = D i Y (n;lm) .
Namely, C ij and D ij assume the traceless and trace parts, respectively. Considering the internal products, the normalized scalar type harmonics are defined by When we consider spherical (pseudo-spherical) case, the spacial curvature takes K = 1 (K = −1).

b. vector type
The vector type harmonics contribute to the transverse modes of the metric perturbation. Namely, the divergences of these harmonics are vanished.
The vector quantities includes two types of harmonics. One is odd parity mode A i whose explicit form is given by A (n;lm) i = 0, f (χ) X (n;l) φ (lm) A . (A29) where, the function f (χ) is defined in (2.7). Then, the eigenvalues are given by The other is even parity mode B i whose explicit form is given by where, ǫ ijk is Levi-Civita symbol associated with γ ij . The eigenvalues are as same as odd ones. The tensor quantities can be constructed by taking symmetrized gradient of each vector quantities :  To recap, we can eliminate the following perturbation modes by choosing gauge : (C17)