Stable Singularity-free Cosmological Solutions in non-projectable Horava-Lifshitz Gravity
Mitsuhiro Fukushima Yosuke Misonoh Shoichiro Miyashita Seiga Sato


WU-AP/1809/18







m.fukushima"at"aoni.waseda.jp
Department of Physics, Waseda University, 
Okubo 3-4-1, Shinjuku, Tokyo 169-8555, Japan


y"underscore"misonou"at"moegi.waseda.jp
Department of Physics, Waseda University, 
Okubo 3-4-1, Shinjuku, Tokyo 169-8555, Japan


miyashita"at"gravity.phys.waseda.ac.jp
Department of Physics, Waseda University, 
Okubo 3-4-1, Shinjuku, Tokyo 169-8555, Japan


s.seiga"at"gravity.phys.waseda.ac.jp
Department of Physics, Waseda University, 
Okubo 3-4-1, Shinjuku, Tokyo 169-8555, Japan














We find stable singularity-free cosmological solutions 
in non-flat Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime in the context of Hořava-Lifshitz (HL) theory.
Although we encounter the negative squared effective masses of the scalar perturbations in the original HL theory,
the behaviors can be remedied by relaxing the projectability condition.
In our analysis, the effects from the background dynamics are taken into account 
as well as the sign of the coefficients in the quadratic action for perturbations. 
More specifically, we give further classification of the gradient stability/instability into five types.
These types are defined in terms of the effective squared masses of perturbations ℳ^2, 
the effective friction coefficients in perturbation equations ℋ and
these magnitude relations |ℳ^2|/ℋ^2.
Furthermore, 
we indicate that oscillating solutions possibly show 
a kind of resonance especially in open FLRW spacetime.
We find that the higher order spatial curvature terms with Lifshitz scaling z=3 are significant
to suppress the instabilities due to the background dynamics.








04.60.-m, 98.80.-k, 98.80.Cq 


Stable Singularity-free Cosmological Solutions in non-projectable Hořava-Lifshitz Gravity

    Seiga Sato
    December 30, 2023
===========================================================================================













§ INTRODUCTION


From the invention of Big Bang theory, 
resolving the cosmological initial singularity problem has been one of the most intriguing 
research frontier in theoretical physics.
According to the singularity theorem proved by Hawking and Penrose<cit.>, 
a spacetime singularity must exist at  a finite past of the Big Bang Universe 
based on the general relativity (GR) unless some unnatural conditions are imposed.
Even if the inflation which resolves the fine tuning problems in the very early stage 
of the Universe is introduced, the initial singularity cannot be remedied<cit.>.

Since the initial singularity spoils the predictability
at the beginning of the Universe, 
many researchers have been proposed cosmological scenarios to remove 
the singularity. 
For example, in the context of braneworld<cit.>,
string gas<cit.>, 
loop quantum gravity<cit.>,
Horndeski theory (or generalized Galileon)<cit.>,
and the non-local gravitational theory<cit.>.
In spite of these efforts, we have not achieved the complete singularity-free
cosmological scenario, yet.
Recently, striking studies have been conducted.
A quite wide class of singularity-free cosmological solutions 
in Friedmann-Lemaître-Robertson-Walker (FLRW)  spacetime
is proved to be unstable<cit.>.
Since the no-go theorem is established based on the Horndeski theory, 
in other words, the most generalized scalar-tensor theory 
whose equation of motion is up to second order<cit.>,
one may consider it is difficult to find stable cosmological solutions
without a singularity.
However, certain loopholes of the no-go theorem are known<cit.>.
One example of such loophole is to consider a Lorenz violating gravitational theory. 
The no-go theorem is possibly violated by introducing higher order spatial derivative
terms<cit.>.


Hořava-Lifshitz (HL) theory is known as a gravitational theory without Lorentz symmetry, 
which is a candidate for quantum gravity<cit.> (the recent progresses of the HL theory are reviewed in <cit.>).
As we have known, the spin-2 gravitational field cannot be quantized 
in perturbative approach, 
that is, there is infinite number of loop diagrams with ultraviolet divergence.
In contrast, the HL theory realizes renormalizability 
at least at power-counting level by introducing the Lifshitz scaling, 

    t→ b^-z t  and  x^i → b^-1x^i  ,

which is an anisotropic scaling between time t and space x^i.
A dynamical critical exponent z characterizes a degree of anisotropy.
If z is equal to or greater than the spatial dimension, 
ultraviolet divergence can be suppressed by finite number of counter terms<cit.>.


In the context of the HL theory, a number of attractive applications 
for cosmology have been conducted (see <cit.> for a review): 
for example, primordial perturbations<cit.>, 
gravitational waves<cit.> 
and other cosmological aspects<cit.>.
It is remarkable that the bouncing and oscillating solutions are discovered
as singularity-free cosmological models in HL theory<cit.>.
Lifshitz scaling terms up to z=3 which realizes
the power-counting renormalizability in four-dimensional spacetime 
derive squared and cubed three-dimensional Ricci curvatures in the action.
These terms behave as radiation-like and stiff matter-like components in non-flat 
FLRW spacetime.
As we will mention, these components can effectively violate the energy condition,  
which is one of the postulates in the singularity theorem.

In our previous paper<cit.>, we investigated the stability of the singularity-free
cosmological solutions based on the projectable HL theory.
Namely, we impose following the projectability condition:
 
    N(t, x^i) → N(t)  ,

 where, 
N is the lapse function which is one of the Arnowitt-Deser-Misner (ADM) variables.
It turns out that the HL theory with the projectability condition suffers from the instability 
when the effects of higher order curvature 
terms are irrelevant<cit.>.
More specifically, we encounter the negative squared effective masses 
of scalar perturbations in FLRW spacetime.
To suppress the instability, we must consider a strong Hubble friction 
by introducing a positive cosmological constant. 

In fact, the pathological behavior in infrared region is possibly remedied by considering a theory without the projectability condition,
that is the non-projectable HL theory<cit.>.

It is discovered that the gradient instability in flat FLRW spacetime can be avoided<cit.>.
Therefore it is expected that the infrared instabilities of the singularity-free solutions 
in non-flat FLRW spacetime are improved.


Although the non-projectable HL theory possesses 
the attractive feature<cit.>, 
the most general form of the gravitational action is extremely complicated.
Since what we would like to clarify is to examine whether the non-projectable HL theory can improve
the infrared behavior after singularity avoidance or not, 
we will focus only on an additional term which is dominant in infrared regime, 
which is the minimally extended HL action.

The rest of this paper is organized as follows.
In Sec. <ref>, 
we construct the non-projectable HL action with the minimal extension from the projectable one.
The background dynamics in non-flat FLRW spacetime and the classification of the singularity-free solutions are
reviewed in Sec. <ref>.
In Sec. <ref>, 
we perform the perturbative analysis by deriving the quadratic action.
The decomposition of the perturbation modes and the manner of the gauge fixing are also summarized.
Then, in Sec. <ref>, we discuss the stability of singularity-free solutions in non-flat FLRW spacetime.
We give further classification of the gradient stability/instability into five types
to consider the background effects. 
Sec. <ref> is devoted to summary and discussions.











§ MINIMALLY EXTENDED NON-PROJECTABLE HOŘAVA-LIFSHITZ THEORY

 
As we mentioned, Lifshitz scaling with z >1 induces the anisotropy 
between space and time, which means 
the general covariance in four-dimensional spacetime is no longer valid.
We instead find a fundamental symmetry called foliation-preserving diffeomorphism:

    t → t + f(t) , x^i → x^i + ζ^i (t,x^j)  ,

namely, a boost transformation is prohibited. 
Thus, it is clearly preferable to adopt the ADM quantities, the three-dimensional spatial metric g_ij,
the lapse function N and the shift vector N^i as the fundamental variables.
In order to preserve the invariance under (<ref>), the terms in action must be composed of the following quantities: 
The  extrinsic curvature which is defined by
 
    𝒦_ij := 1  2N( ġ_ij - ∇_i N_j - ∇_i N_j ) ,

 where, the dot symbol represents the derivative with respect to time coordinate and 
 ∇_i denotes the spatial covariant derivative. 
 The three dimensional Ricci curvature ℛ_ij associated with 
 the spatial covariant derivative. 
 Since we relax the projectability condition, the spatial dependence of the lapse function
 is restored. 
 Thus, the following vector quantity can be included in the action:
  
    Φ_i := ∇_i N  N ,
 
 which represents the acceleration of rest observer on three-dimensional hypersurface.
Furthermore, spatial covariant derivatives of these quantities also form 
the invariant action, i.e., ∇_i 𝒦_jk, ∇_i ℛ_jk, ∇_i Φ_j, ∇_i ∇_j 𝒦_kl
 and so on.
 
Although we can include every scalar quantities which are composed of these variables 
to construct the invariant action, our model is restricted as follows:

    S = ^2  2∫ dt  d^3x   ( ℒ_K
    +ℒ_P +ℒ_NP)  ,

 with
   
    ℒ_K   :=    N√(g)( 𝒦_ij𝒦^ij -λ𝒦^2  )  
      ,
       
    
    ℒ_P   :=    -N √(g)[ 
     𝒱_z=1 + ^-2𝒱_z=2 + ^-4𝒱_z=3]
      , 
     
    
    ℒ_NP   :=    N √(g) 𝒱 [Φ_i]  , 
    
    𝒱_z=1   :=    2Λ + g_1ℛ , 
    𝒱_z=2   :=    g_2 ℛ^2  + g_3 ℛ^i_ jℛ^j_ i , 
    𝒱_z=3   :=    g_4 ℛ^3 + g_5 ℛ ℛ^i_ jℛ^j_ i 
    +g_6 ℛ^i_ jℛ^j_ kℛ^k_ i
          +g_7 ℛ∇^2 ℛ +g_8 ∇_iℛ_jk∇^iℛ^jk ,

where, m_LV is a Lorentz violating scale which is expected to be around the Planck scale, 
λ and g_i (i=1-8) are dimension-less coupling constants 
and Λ is a cosmological constant.
ℒ_NP  is constructed by the terms including Φ_i field,
which is with no effect to FLRW background.
The reason for the restriction is to ensure the comparability of our previous result based on the SVW action<cit.>.
As we will see in the next section, 
the background dynamics are identical to those of the SVW action under specific condition.

We further impose a restriction on ℒ_NP.
Recall that the purpose of this paper is to remedy the infrared behavior of bouncing solutions
discovered in the projectable HL theory.
It is naturally expected that the terms with the lowest order of spatial derivative
is essential to stabilize the infrared region. 
In fact, such a term is uniquely determined, that is Φ^2 := Φ_i Φ^i.
This diagnosis seems not to be inconsequent.
According to the papers <cit.>, the stabilities of  Minkowski and flat FLRW spacetime 
are remedied by the non-projectable HL theory with Φ^2 term.
Therefore, we do minimally extend the theory by adding only Φ^2 term: 

    𝒱[Φ_i] := ς Φ^2  .

with a dimensionless coupling constant ς.


In the rest of this paper, the unit = 1 is adopted.
We additionally can set the value of the coupling constant by rescaling the
 time coordinate.
Thus, we set g_1 = -1, which is equivalent to take a coordinate in which
the propagating speed of the spin-2 gravitational wave in infrared limit is precisely unity<cit.>.








§ BACKGROUND SOLUTIONS IN FLRW SPACETIME


The one of the striking differences from the projectable HL theory is the structure of the basic equation system,
that is, the localness of the Hamiltonian constraint is recovered by relaxing the projectability condition.
In this section, we briefly summarize the structure of the basic equation system and
the classification of singularity-free solutions in FLRW spacetime.



 §.§ Basic equations

To unify the description, we consider the FLRW spacetime in a spatial coordinate system x^i = (χ, θ, ϕ):

    ds^2 
       =    g_ij dx^i dx^j 
       =     a^2 [ d χ ^2 + f(χ)^2 (dθ^2 + sin^2θ dϕ^2) ]  ,

with 

    f(χ) := 
    χ   for K=0 
    sinχ   for K=1 
    sinhχ   for K=-1 
     ,

where, a is the scale factor which depends only on the cosmic time, 
and K is related to the spatial Ricci curvature as ℛ = 6 K/ a^2.
The cases with K=0, 1 and -1 correspond to the flat, closed and open FLRW spacetime, respectively.
The domain of the coordinate χ is defined by 0 ≤χ < ∞ for K=0, -1 and
0 ≤χ≤π for K=1.
Furthermore, angular coordinates θ and ϕ take 0 ≤θ≤π 
and 0 ≤ϕ≤ 2π, respectively.

Then, the dynamical equation for the scale factor and the Friedmann equation are 
obtained by taking variation of the action with respect to a and N:  

    2 Ḣ +  3 H^2  = 
      2  3λ -1( Λ - K  a^2 - g_r K^2  3a^4 - g_s K^3  a^6)  ,  
    
    
             H^2  =  2   3 (3λ-1) [ Λ -3K  a^2 +g_r K^2  a^4 + g_s K^3  a^6]  ,

where, H:=ȧ/a is the Hubble parameter, 
g_r and g_s are the linear combinations 
of the coupling constants[
It should be noted that we adopt the different definitions of g_r and g_s from
our previous paper<cit.> to simplify the perturbed action.
The previous definitions are

    g_r := 6K^2 (3g_2 +g_3) ,  
    g_s := 12 K^3 (9g_4 +3g_5 +g_6)  ,

thus, the sign of g_s is flipped in open FLRW spacetime.
]:

    g_r    :=    6 (3g_2 +g_3) , 
     
    g_s    :=    12 (9g_4 +3g_5 +g_6)  .

We have already fixed the gauge as N=1 and N_i=0.
In what follows, we assume λ > 1/3.

As we know, the equation of the scale factor (<ref>) can be derived by
taking derivative of  (<ref>) with respect to the cosmic time.
It means the independent equation is only (<ref>). 
However, the situation with the projectability condition is quite different.
Due to the lack of the local lapse function, 
the Hamiltonian constraint is an integration over whole space. 
Therefore, we have to adopt the scale factor equation instead of 
the global Hamiltonian constraint in the projectable case. 
Friedmann-like equation is derived by integration with respect to the cosmic time:

    H^2  =  2   3 (3λ-1) [ Λ -3K  a^2 +g_r K^2  a^4 + g_s K^3  a^6 + 𝒞  a^3]  ,
 
with dust-like term with an integration constant 𝒞<cit.>.
Thus, if we consider 𝒞 = 0 case in the projectable HL theory, 
the same background FLRW solutions are realized<cit.>.



 §.§ Singularity-free background solutions


In order to investigate the background solutions, it is convenient to rewrite 
the Friedmann equation (<ref>) into

    1  2ȧ^2 + 𝒰(a) =0 ,
 
with

    𝒰 (a) = 1  3λ-1 [ K -Λ 3a^2  -g_r K^2  3a^2 -g_s K^3  3a^4]  .
 
Since the possible ranges for the scale factor are limited within the region in which 𝒰≤ 0, 
the background evolution is completely
determined by the coupling constants in ℒ_P.
It is remarkable that the g_r and the g_s terms simulate a radiation component
and a stiff matter component, respectively.
The important point is that these terms can effectively violate the energy condition
if either or both of g_r K^2 and g_sK^3 are negative, 
which may lead singularity-free solutions even if we do not introduce some exotic matters.

The classification of the possible solutions
is given in Ref.<cit.>.
In our analysis, we focus on the following two kinds of singularity-free solutions.
One is a bouncing solution 
denoted by ℬ^[ sgn(Λ); K],
where the function sgn(x) gives a sign of x. 
An initial contracting universe shifts to expanding phase at a=a_T
and keep expanding forever.
The other is oscillating solution 
denoted by 𝒪^[ sgn(Λ); K].
A universe shows periodic oscillating behavior, in other words, 
bounces at a=a_min and recollapses[
In this paper, a word of recollapse means that the expanding universe shifts to contracting phase.
Thus, a recollapcing point satisfies the conditions 𝒰=0 and 𝒰' >0.
] at a=a_max.
Therefore the scale factor is limited within 0<a_min≤ a ≤ a_max<∞.
As we will see, the typical size of the oscillating amplitude is expected to be 
the Lorentz violating scale, thus, it is difficult to represent 
the cyclic universe scenario whose maximum scale factor is macroscopic.



  §.§.§ Without cosmological constant

In our previous paper, we have found that the stable singularity-free solutions require 
a positive cosmological constant based on the projectable HL theory.
Once the projectability condition is relaxed, it is expected
to find stable bouncing solutions without a cosmological constant
as is the case in flat FLRW spacetime. 
Therefore, we summarize the singularity-free cosmological solutions 
without a cosmological constant in this part. 

Since the sign of 𝒰 determines the possible ranges for the scale factor, 
it is convenient to derive the roots of the equation 𝒰=0:

    a^[K]_± = √(K  6( g_r  ±√(g_r^2  + 12g_s ))  ) .
 
The points a=a^[K]_± correspond to the bouncing or recollapsing points
of the universe.
Of course, the corresponding a^[K]_± must be real and positive to find such points.
By examining the forms of the potential 𝒰, 
we find the following three types of  singularity-free solutions. 



  (a) 
𝒪^[0;1]— 
A universe shows oscillating behavior 
in closed FLRW spacetime (K=1), which we call 𝒪^[0;1].
This type of the solutions can be found if the following conditions are satisfied.
 
    g_s <0 , g_r >0   and   g_r^2 +12 g_s >0
                      with  a_max = a^[1]_+  and  a_min = a^[1]_- .

In open FLRW spacetime, this kind of solutions is never found. 


  (b) 
 ℬ^[0;-1]— 
An initially contracting universe shifts to expanding phase at bouncing point a=a_T, 
and keeps expanding forever in open FLRW spacetime (K=-1).
This type of the solutions is classified into the following two cases: 

    (i)      ∀ g_r  and  g_s >0   with  a_T = a^[-1]_- . 
    
    
      (ii)        g_r <0   and   g_s =0   with  a_T = a^[-1]_- .

The solutions satisfying the conditions (i) and (ii) are referred to as ℬ^[0;-1](i) and
ℬ^[0;-1](ii), respectively.
In closed FLRW spacetime, we cannot find this type of the solutions.
 
 
  (c) ℬ^[0;-1]_BC—
A universe shows bouncing behavior if the initial value of scale factor 
a_ini is larger than a_T,
whereas, falls down into the singularity if a_ini≤ a_BC.
Since this type of the solutions possibly induces a Big Crunch (BC),  
we add a subscript BC.
Note that the domain a_BC < a < a_T is prohibited 
because the Hamiltonian constraint is never satisfied.
If the following conditions are satisfied in open FLRW spacetime (K=-1), 
ℬ^[0;-1]_BC can be found: 
 
    g_s <0 , g_r <0  and  g_r^2 +12 g_s >0
                       with  a_T = a^[-1]_-  and  a_BC = a^[-1]_+ ,

 where, a_BC is a recollapsing point of the Big Crunch universe.  
 This type of bouncing solutions is found only in open FLRW spacetime.
 
We show the typical forms of potentials of 
𝒪^[0,1], ℬ^[0,-1] and ℬ^[0,-1]_BC
in FIG. <ref>
and the distribution of these types of solutions on (g_r, g_s) plane 
in FIG. <ref>.
 
 

As we will discuss later, 
we can construct the bouncing solutions 
ℬ^[0;-1]  and ℬ^[0;-1]_BC
whose squared effective masses of scalar perturbations are positive.
Therefore, there is a possibility to construct stable singularity-free solutions without
a cosmological constant. 
It is a striking difference from the projectable HL theory.



  §.§.§ With positive cosmological constant

Although the properties of singularity-free solutions for Λ >0 have already given 
in the papers <cit.>, 
we again summarize the solutions because the notation is slightly changed from previous ones.
When we consider the case with non-zero cosmological constant, 
it is convenient to introduce the following quantities normalized by a cosmological constant; 
ã := a/ℓ, g̃_r := g_r/ℓ^2 and g̃_s := g_s /ℓ^4 with ℓ := √( 3 / |Λ| ). 
Then, the potential 𝒰 is rewritten as
 
    𝒰̃ (ã) = 1  3λ -1[ K - εã^2 - g̃_r K^2 3 ã^2 - g̃_s K^3 3 ã^4]  ,

where ε:=±1 expresses the sign of the cosmological constant.
The three roots of the equation 𝒰̃ (ã) =0 are given by 
 
    ã^[ε; K]_I     :=   √(K  6  ε[ 2 + 4(1-g̃_r) ξ̃^[ε; K]_I  
    +  ξ̃^[ε; K]_I  ]  ) ,

with 

    ξ̃^[ε; K]_I    :=    2^2/3 (e^ 2 π i / 3 )^I  pv[ 2 -3εg̃_r  
    - 9g̃_s  
          +9 sgn(K)√((g̃_s -g̃^[ε; K]_s(+)) 
    (g̃_s -g̃^[ε; K]_s(-)) )]^1/3 , 
    
    g̃^[ε; K]_s(±)   :=   1  9[ 2 -3 εg̃_r  ± 2  sgn(K) ( 1 - εg̃_r )^3/2]   ,
 
where, pv means a principal value of cubic root
and I =1, 2 , 3. 
 If ã^[ε; K]_I takes real and positive value,  
𝒰̃ (ã) =0 possesses 
a corresponding root 
in the possible ranges for the scale factor. 
We further 
derive the roots of 𝒰̃ (ã) =0 with g̃_s =0: 
 
    ã_±^[ε; K] := √(K  2 ε[ 1 ±sgn(K) √( 1 -4 εg̃_r   3)] ) .
 

We classify the singularity-free solutions with a positive cosmological constant 
(ε=+1) into the following three types: 


  (a) ℬ^[1;K]—
A bouncing solution whose bouncing radius is given as ã_T.
We refer to this type of the solutions as ℬ^[1;K].
Unlike the case without a cosmological constant, 
the universe after the bounce approaches de Sitter spacetime.
 For closed FLRW spacetime (K=1), we can find ℬ^[1;1] as the following three cases: 

    (i)      g̃_s(+)^[1; 1]< g̃_s < 0  and  g̃_r<1   with  ã_T = ã^[1; 1]_1 . 
    
    
     (ii)      g̃_s <0   and  g̃_s< g̃_s(-)^[1; 1]   for  |2g̃_r-1|<1
    g̃_s <0    for  |2g̃_r-1| ≥ 1
    
                                                with ã_T = ã^[1; 1]_3 . 
    
    
    (iii)      g̃_s=0  and  g̃_r ≤ 0   with  ã_T = ã^[1; 1]_+ .

The solutions satisfying the conditions (i), (ii) and (iii) are referred to as ℬ^[1;1](i), ℬ^[1;1](ii) and
 ℬ^[1;1](iii), respectively.
 For open FLRW spacetime (K=-1), we find ℬ^[1; -1] as the following two cases: 

    (i)       g̃_s >0  with  ã_T = ã^[1; -1]_3  . 
    
    
     (ii)       g̃_s =0  and  g̃_r <0   with  ã_T = ã^[1; -1]_+  .

ℬ^[1;-1](i) and ℬ^[1;-1](ii) correspond to the conditions (i) and (ii), respectively.



  (b)  ℬ^[1; K]_BC—
 A universe with ã_ini≤ã_BC
 evolves into a Big Crunch, while it shows bouncing behavior 
 if ã_ini≥ã_T  > ã_BC.
 We classify this type of the solutions as ℬ^[1; K]_BC.
 Unlike ℬ^[0;K]_BC, the asymptotic behavior of spacetime after bounce is de Sitter spacetime.
 For closed FLRW spacetime (K=1), we find the following two cases: 

    (i)      0 < g̃_s < g̃_s(+)^[1; 1]
                    with  ã_BC = ã^[1; 1]_2  and  ã_T = ã^[1; 1]_3 . 
    
    
     (ii)      g̃_s=0  and  0<g̃_r< 3/4
                    with  ã_BC = ã^[1; 1]_-   and  ã_T = ã^[1; 1]_+ .

The solutions satisfying the conditions (i) and (ii) are referred to as ℬ^[1; 1]_BC(i) and
 ℬ^[1; 1]_BC(ii), respectively.
 For open FLRW spacetime (K=-1), the solutions ℬ^[1; -1]_BC can be seen 
when the following conditions are satisfied: 

    0 > g̃_s > g̃_s(+)^[1; -1]  and  g̃_r<0  
                     with  ã_BC= ã_2^[1; -1]  and  ã_T= ã_3^[1; -1] .

 

  (c) ℬ^[1;1]_O— 
A universe in closed FLRW spacetime shows oscillating behavior if the initial scale factor is 
in ã_min≤ã_ini≤ã_max, 
or bounces if ã_ini≥ã_T > ã_max.
We classify this type of solutions as ℬ^[1; 1]_O. 
The subscript represents the oscillating behavior.
This type of the solutions can be found in closed FLRW spacetime (K=1)
if the following conditions are satisfied.

    g̃_s <0   and  g̃^[1; 1]_s(-) < g̃_s < g̃^[1; 1]_s(+)
          with  ã_min = ã^[1; 1]_1 , ã_max = ã^[1; 1]_2  and  ã_T = ã^[1; 1]_3 .
 
For open FLRW spacetime, we never find this type of solutions.

The typical forms of potentials and the distribution of the singularity-free solutions on (g̃_r, g̃_s) plane 
are shown in FIG. <ref> and FIG. <ref>, respectively.
 
 



  §.§.§ With negative cosmological constant
 
In this paper, we discuss the stabilities of the oscillating solutions 
as well as the bouncing solutions.
Thus, we do not exclude the solutions with a negative cosmological constant (ε=-1),
which we cannot construct the bouncing solutions.
We find the following two kinds of singularity-free oscillating solutions:


 
  (a) 𝒪^[-1;K]— 
A universe which shows periodic oscillation whose oscillating amplitude is 
given by ã_min≤ã≤ã_max.
We classify this type of the solutions as 𝒪^[-1;K].
If the following conditions are satisfied in closed FLRW spacetime (K=1), 𝒪^[-1;1] is realized: 
 
    g̃_s(-)^[-1;1] <g̃_s <0   
                 with  ã_min=ã_2^[-1;1]  and  ã_max=ã_1^[-1;1] .

In open FLRW spacetime (K=-1), we find the two cases of 𝒪^[-1;-1] are obtained:
 
    (i)       0< g̃_s < g̃_s(-)^[-1; -1]  
           with  ã_min=ã_2^[-1;-1]  and  ã_max=ã_1^[-1;-1] . 
    
     (ii)      g̃_s =0  and  -3/4< g̃_r < 0 
           with  ã_min=ã_+^[-1;-1]  and  ã_max=ã_-^[-1;-1] .

 The solutions satisfying the conditions (i) and (ii) are referred to as 𝒪^[-1;-1](i) and 𝒪^[-1;-1](ii), respectively.

 
  (b) 𝒪_BC^[-1;-1]— 
A universe in open FLRW spacetime shows oscillating behavior if the initial radius of the universe is in
ã_min≤ã_ini≤ã_max, 
or falls into the singularity if ã_ini<ã_BC<ã_min.
We refer to this type of solutions as 𝒪_BC^[-1;-1].
For open FLRW spacetime (K=-1), we find the solutions under the following conditions:

    g̃_s<0   and  g̃_s(+)^[-1;-1] <g̃_s <g̃_s(-)^[-1;-1]
                 with  ã_BC=ã_3^[-1;-1] , ã_min=ã_2^[-1;-1]
                         and  ã_max=ã_1^[-1;-1] .

In closed FLRW spacetime, we cannot construct this type of solutions. 

The typical forms of potentials and the distribution of the singularity-free solutions on (g̃_r, g̃_s) plane 
are shown in FIG. <ref> and FIG. <ref>, respectively.
 
 












§ PERTURBATION ANALYSIS AROUND A NON-FLAT FLRW BACKGROUND


In this section, we derive the perturbed quadratic action of the minimally extended non-projectable HL theory.
The perturbed ADM variables are defined by 

    N    =   N̅ + δ N , 
    
    N_i    =   N̅_i + δ N_i , 
    
    g_ij   =   g̅_ij + δ g_ij ,

with

    δ N    =   α , 
    δ N_i    =   β_i  , 
    δ g_ij   =     h_ij +1  2g̅^ab h_a i h_b j ,

where, N̅, N̅_i and g̅_ij denote
the background lapse function, shift vector and three-dimensional induced metric, respectively.
Furthermore, α, β_i and h_ij mean first order 
perturbation of the above ADM variables.
Note that indices of these perturbed variables are raised by g̅^ij
as h^i_j:=g̅^iah_aj, h^ij:=g̅^iag̅^jbh_ab, 
h:=h^i_i=g̅^iah_ia and β^i:=g̅^iaβ_a.



 §.§ Harmonic expansion

To decompose the perturbations into scalar, vector and tensor modes, 
we perform the harmonic expansion by equipping the set of tensor harmonics: 

    Y   =   {Q , Q_i , Q_ij , P_ij , S_(o)i , S_(e)i , S_(o)ij , S_(e)ij , 
    
            G_(o)ij , G_(e)ij} ,

where, Q , Q_i , Q_ij , P_ij are scalar type harmonics. 
Q_ij and P_ij are trace and trace-less part, respectively. 
S_(o)i , S_(e)i , S_(o)ij , S_(e)ij are vector type harmonics. 
The symbols (o) and (e) correspond to the odd parity and the even parity, respectively.
G_(o)ij , G_(e)ij are tensor type harmonics with the odd and the even parity.
The explicit forms and these characteristics can be seen in Appendix A and B in Ref. <cit.>.

The scalar perturbations of the ADM variables are decomposed into 

    α^(scalar)   =   ∑_n,l,mα_(Q)^(n;lm)Q^(n;lm) , 
    β_i^(scalar)   =   ∑_n,l,m  a^2 [ β^(n;lm)_(Q) Q^(n;lm)_i]  , 
    
    h_ij^(scalar)   =   ∑_n,l,m a^2 [ h^(n;lm)_(Q) Q^(n;lm)_ij +h^(n;lm)_(P) P^(n;lm)_ij]  ,
 
the vector perturbations also can be expanded by 

    β_i^(vector)   =   ∑_n,l,m a^2 [ β^(n;lm)_(S;o) S^(n;lm)_(o) i+β^(n;lm)_(S;e) S^(n;lm)_(e) i]  , 
    
     
    h_ij^(vector)   =   ∑_n,l,m a^2 [ h^(n;lm)_(S;o) S^(n;lm)_(o)ij +h^(n;lm)_(S;e) S^(n;lm)_(e)ij]  ,
 
and the tensor part is given by

    h_ij^(tensor)   =   ∑_n,l,m a^2 [ h^(n;lm)_(G;o) G^(n;lm)_(o)ij +h^(n;lm)_(G;e) G^(n;lm)_(e)ij]  .
 
where, the degrees l, m ∈ℤ are constrained by 0 ≤ l ≤ n-1 and 0 ≤ |m| ≤ l.
n ≥ 1 is a continuous real number for K=0, -1, 
while a discrete natural number for K=1.

Since we relax the projectability condition, the localness of the lapse perturbation α
is recovered.
Thus, we expand α^(scalar) by the harmonic functions.



 §.§ Gauge fixing


Before calculating quadratic action, we have to remove the gauge degree of freedom 
from the perturbations.
Since the HL theory respects the foliation-preserving diffeomorphism (<ref>), 
the infinitesimal transformations of the perturbed ADM quantities are given by 

    α^(gauge)   =   ∂_t f , 
    β_i ^(gauge)   =   ∂_t ζ_i -2 H ζ_i  , 
    
    h_ij^(gauge)   =    2∇_(iζ_j) + 2 H f g̅_ij .

We stress that f does not depend on space, thus, 
only ζ^i can be expanded by the harmonics:

    ζ_i    =   ∑_n,l,m a^2 [ ζ^(n;lm)_(Q) Q_i^(n;lm)
              + ζ^(n;lm)_(S;o) S_(o)i^(n;lm)  + ζ^(n;lm)_(S;e) S_(e)i^(n;lm)] .

Then, the infinitesimal gauge transformations of the harmonic expansion
 (<ref>)-(<ref>) are given by

    ∑_n,l,mα^(n;lm)_(Q) Q^(n;lm)   →   ∑_n,l,mα^(n;lm)_(Q) Q^(n;lm) + ∂_t f  , 
    
    
    β^(n;lm)_(Q)   →   β^(n;lm)_(Q)+ ∂_t ζ_(Q)^(n;lm) ,  
    
    
    h^(n;lm)_(P)   →    h^(n;lm)_(P) +√(8 (ν^2 -3K)  3)ζ_(Q)^(n;lm) ,


    ∑_n,l,m a^2 h_(Q)^(n;lm) Q^(n;lm)_ij
          →∑_n,l,m a^2 [ h^(n;lm)_(Q)-2 ν√(3)ζ_(Q)^(n;lm)] Q^(n;lm)_ij 
    +2 H f g̅_ij ,


    β^(n;lm)_(S)   →   β^(n;lm)_(S) + ∂_t ζ_(S)^(n;lm) ,  
    
    
    h^(n;lm)_(S)   →    h^(n;lm)_(S) +√(2 (ν^2 -3K) )ζ_(S)^(n;lm) ,


    h^(n;lm)_(G)   →    h^(n;lm)_(G) .

where, ν^2 is a eigenvalue of the harmonics which is defined by the following regions (n≥1):

    ν^2 := 
    
    n^2 -1 ,  n ∈   for  K=0 
    
    n^2 -1 ,  n∈   for  K=1 
    
    n^2 +1 ,  n ∈   for  K=-1 
     .

Since both odd and  even parity modes obey the same transformation law, 
the parity subscripts are abbreviated in vector and tensor perturbations.

We shall eliminate h^(n;lm)_(P), h^(n;lm)_(S;o) and h^(n;lm)_(S;e) by choosing the following gauge 

    ζ^(n;lm)_(Q)   =    - [ 8 (ν^2 -3K)  3]^-1/2 h^(n;lm)_(P) ,
    ζ^(n;lm)_(S;o)   =    -[2(ν^2-3K) ]^-1/2 h_(S;o)^(n;lm) ,
    ζ^(n;lm)_(S;e)   =    -[2(ν^2-3K) ]^-1/2 h_(S;e)^(n;lm) .

Unlike the projectable case, we cannot eliminate α^(n;lm)_(Q).
In what follows, we have abbreviated the superscript (n;lm) 
because the perturbed quantities do not mix with different modes.



 §.§ Quadratic action

As is the case in the projectable HL theory, the quadratic action can be decomposed 
into the tensor part and the scalar one. 
The vector perturbations are not dynamical.



  §.§.§ Tensor perturbations

Since the additional perturbation terms coming from ℒ_NP are only scalar modes, 
the tensor part of the quadratic action is identical to that of the projectable case: 

    δ_(2)ℒ^(tensor)   =   a^3  2[ ℱ_(G)ḣ_(G)^2 - 𝒢_(G) h_(G)^2 ] ,
 
where, we introduced ℱ_(G) and 𝒢_(G) which can be regarded as the kinetic term and the mass term
of the tensor perturbations, respectively.
These variables, except for the total derivative terms, are given by
 
    ℱ_(G)   :=    1  , 
    
    𝒢_(G)   :=   ν^2  a^2  +ν^2  3a^4[ -2 g_r  K +3g_3 ν^2  ] 
    
          
    +ν^2  a^6[ 
    -g_s K^2 +6 g_56 K ν^2 
    +g_8 ν^2 (ν^2 -2K)
     ] ,

and we define g_56:=g_5+ g_6.
The tensor quadratic action is defined for the case with l ≥ 2, 
because the tensor harmonics G^(n;lm)_(o)ij and G^(n;lm)_(e)ij are vanished when l < 2.
 


  §.§.§ Scalar perturbations

The scalar perturbations are drastically changed since the lapse perturbation α_(Q) 
cannot be eliminated by the gauge condition.
Furthermore, ℒ_NP adds a new degree of freedom to scalar perturbations. 
The quadratic action of the scalar perturbations is given by


    δ_(2)ℒ^(scalar)   =   
    -a^3  2(3λ-1) ḣ_(Q)^2 
    +a/3(ν^2-3K) h_(Q)^2
    -1/27a(ν^2-3K) [ 2g_r(2ν^2-3K) + 3g_3ν^2 ]h_(Q)^2 
            
    -1/9a^3(ν^2-3K) [ g_s K (4ν^2-9K) +2(3g_56-4g_7)Kν^2 + (-8g_7+3g_8)ν^2(3ν^2-10K) ]h_(Q)^2 
    
            
    + 2√(3) a^3 (3λ-1)Hḣ_(Q)α_(Q) 
    +4a/√(3)(ν^2-3K)[ 1 - 2/3a^2g_r K - 1/a^4K(g_s K -6g_7 ν^2) ] h_(Q)α_(Q)
            
    +2a[ ςν^2 -3a^2(3λ-1)H^2  ] α_(Q)^2
    -(3λ-1)  2ν√(3)ḣ_(Q)β_(Q) 
    +4a^3(3λ-1)ν H α_(Q)β_(Q)
            
    -2a^3[ (λ-1)ν^2 +2K  ] β_(Q)^2  .

We eliminated Ḣ by applying the scale factor equation (<ref>).  
Since α_(Q) and β_(Q) are not dynamical, 
we can eliminate both of them with the following constraint equations: 

    β_(Q)   =    (3λ-1)ν/ (λ-1)ν^2 + 2K [ Hα_(Q) - 1/2√(3)ḣ_(Q)]  , 
    α_(Q)   =   (ν^2-3K)[ 
    ( (λ-1)ν^2 + 2K ) ( 3a^4 - 2 g_r K a^2 -3K(g_s K -6g_7 ν^2)  ) h_(Q)
    - 3(3λ-1)Ha^6 ḣ_(Q)] /3√(3)a^4 [ ςν^2( (λ-1)ν^2 +2K ) +2(3λ-1)(ν^2-3K)H^2 a^2 ] ,

Substituting, (<ref>) and (<ref>) into (<ref>), we achieve the quadratic action of  scalar perturbations:

    δ_(2)ℒ^(scalar)   =   a^3  2[ ℱ_(Q)ḣ_(Q)^2 - 𝒢_(Q) h_(Q)^2 ] ,
 
where,

    ℱ_(Q)   :=   2ς (3λ-1) ν^2 (ν^2 -3K)  3ςν^2 [ (λ-1)ν^2  +2K ] +6(3λ-1) (ν^2 -3K)H^2 a^2   , 
    
    
    𝒢_(Q)   :=    -2  3a^2(ν^2 -3K)  +2  27 a^4(ν^2 - 3K) [ 2 g_r   ( 2ν^2 - 3K )
    + 3g_3 ν^2  ] 
          +2  9a^6 (ν^2 - 3K) [ g_s K (4ν^2 -9K) +2 (3g_56 -4 g_7)  Kν^2 
    +(-8 g_7 + 3 g_8) ν^2 (3ν^2 - 10K) ] 
           +4 (ν^2 -3K)^2 [ (λ-1)ν^2  +2K ] (  -3g_s K^2 + 18g_7 K ν^2  -2 g_r K a^2 + 3a^4   )^2  27 a^10[ ςν^2 { (λ-1)ν^2  +2K }  +2(3λ-1) (ν^2 -3K)H^2 a^2  ] 
           +  8 K (ν^2 -3K)^2 (  3g_s K^2 +2 g_r K a^2 -3 a^4 )( 3g_s K^2 -18g_7 K ν^2  +2 g_r K a^2 -3a^4   )
      27 a^10[ ςν^2 { (λ-1)ν^2  +2K } +2(3λ-1) (ν^2 -3K)H^2 a^2  ] 
           +  4 (3λ-1) (ν^2 -3K)^2 H^2 
      ( 3g_s K^2 -18g_7Kν^2 -2g_r K a^2 +9 a^4 )
      9 a^4 [ ςν^2 { (λ-1)ν^2  +2K } +2(3λ-1) (ν^2 -3K)H^2 a^2  ] 
           +  32 (3λ-1) (ν^2 -3K)^3 H^2 (  -3g_s K^2 +18g_7 K ν^2  -2 g_r K a^2 +3a^4   ) (  2g_s K^3 +g_r K^2 a^2  -Λ a^6 ) 
      27 a^8 [ ςν^2 { (λ-1)ν^2  +2K } +2(3λ-1) (ν^2 -3K)H^2 a^2  ]^2  .
 

In closed FLRW spacetime, we focus only on the case with n ≥ 3.
Actually, the case with n=1 corresponds to a constant shift of the scale factor, which is less important, and n=2 is not dynamical mode.
Thus, we can regard (ν^2 - 3K) as a positive value.
It should be noted that, even if we take the limit ς→ 0,
we cannot replicate the quadratic scalar action in the projectable case
(see (3.26) and (3.27) in ref.<cit.>).
In fact, a difference is caused by the gauge structure.
As we have seen, α_(Q) is eliminated by applying constraint equation (<ref>).
Due to the recovering the local lapse function, it cannot be removed by 
global infinitesimal transformation f.








§ STABILITY ANALYSIS OF THE SINGULARITY-FREE SOLUTIONS


To examine whether a singularity-free cosmological solution is truly realized or not, 
it is essential to consider its stability.
If the background solution is unstable with respect to small perturbations, 
the possibility of singularity avoidance may be spoiled.
Thus, we examine the stabilities of singularity-free background solutions which are
shown in previous section. 

When we discuss the stability of a specific solution, 
the sign of ℱ and 𝒢 in the quadratic action provides guideposts
(Since the discussion holds for both scalar and tensor modes, the subscripts are abbreviated.).
The sign of ℱ is relevant to the ghost instability.
A perturbation mode with negative ℱ losses the lowest energy state
which leads the fatal collapse of the perturbative approach.
Thus, we consider the condition ℱ >0 as the absolute requirement. 

 
On the other hand, the case with 𝒢<0 is, in general, regarded as the gradient instability.
Such a situation corresponds to the negative squared mass, which seems to induce 
an exponential growth of perturbations.
However, that unstable behavior may be suppressed by the effect of 
background dynamics.
To clarify the effect from the background dynamics, 
we examine the perturbation equation of motions derived by
taking variation of the quadratic action: 
 

    ḧ + 3ℋḣ + ℳ^2 h = 0 ,

where, 

    ℋ := H  + ℱ̇ 3 ℱ , ℳ^2 := 𝒢ℱ .

From the equation, one can see that the perturbation dynamics is 
influenced  by  the effective friction coefficient ℋ as a background effect
as well as the effective squared mass ℳ^2.
The effective friction coefficients of the tensor perturbations 
ℋ_(G) coincide with the Hubble parameter H, 
however, those of the scalar perturbations are not: 

    ℋ_(Q) 
    = H (1  -4( Ḣ + H^2 ) a^2  3ςν^2 
    [ (λ-1)ν^2  +2K   (3λ-1) (ν^2 -3K) ] +6H^2 a^2  )   .

We would like to stress that the positive ℋ generates friction,
on the other hand, the perturbations feel acceleration if ℋ is negative. 
It is not completely determined by the Hubble parameter H unlike the projectable case.
In what follows, the forces caused by ℋ>0 and ℋ<0 are called as
a ℋ-friction and a ℋ-acceleration, respectively.


  Thus, we pay deep attention to the values of ℳ^2 and ℋ to 
discuss the stability of the perturbations.
In this paper, we classify the stabilities and instabilities into the following five types.


  (i)
Negative-ℳ^2 instability: 
As we mentioned, the negative ℳ^2 causes undesirable exponential 
growth unless it is not suppressed by  the ℋ-friction.
Therefore, the unstable behavior due to the negative ℳ^2
is observed in the case with 
(1) ℳ^2<0 and ℋ<0, 
(2) ℳ^2<0 and |ℳ^2| ℋ^2>0 
which means the effect of the third term in LHS of (<ref>) is relatively dominant than 
that of the second term. 
If any one of the perturbation modes satisfies this condition, 
we classify the solution is under a negative-ℳ^2 instability.

  (ii)
ℋ-accelerated instability: 
Even if the squared effective mass is positive, 
it is possible to enhance the oscillating amplitude of perturbations.
When  ℋ <0, 
the second term in LHS of (<ref>) 
reinforces the amplitude rather than friction,  thus  the perturbation experiences a ℋ-acceleration. 
Clearly, the accelerating effect is manifested 
if the positive squared effective mass is relatively  smaller than the effect of ℋ term.
Thus, we regard the solution is under a  ℋ-accelerated instability 
if at least one of the perturbation modes satisfies 0 < ℳ^2 ℋ^2 and ℋ<0.

  (iii)
   ℳ^2-dominated stability:
On the other hand, the negative  ℋ is not problematic 
if it is sufficiently suppressed by the heavy positive effective squared mass.
Thus, we call the solution is under a ℳ^2-dominated stability
if every perturbation mode satisfy ℳ^2 ℋ^2 and ℋ<0.

  (iv)
 ℋ-suppressed stability:
If ℋ>0, we experience a friction effect which suppresses the dynamics of the perturbations.
With this taken into consideration, 
 we can find a stable solution even if ℳ^2 <0.
That is, the effective squared mass is negative, however, the unstable behavior is suppressed by the positive ℋ.
We call the solution is under a 
 ℋ-suppressed stability 
 if every perturbation mode satisfy |ℳ^2| ℋ^2 and ℳ^2 < 0.

  (v)
Complete stability:
Clearly, there is no problematic perturbation dynamics 
if every perturbation mode satisfy ℳ^2 > 0 and ℋ>0.
We call this situation as a complete stability. 

In TABLE <ref>, we summarize the classification of the stabilities and instabilities.


In principle, the types of stabilities are entirely determined by the coupling constants in the action.
However, it is difficult to show the explicit conditions for these stabilities
due to the complicated forms of 𝒢 and ℱ for general cases.
Therefore, our main procedure is to numerically trace 𝒢 and ℱ
based on the target background solutions.
Before starting numerical analysis, we consider several simplified cases 
in which the stability conditions are given in explicit forms.





 §.§ Late-time universe after bounce

When the higher order curvature terms are neglected, the forms of 𝒢 and ℱ are quite simplified.
The tensor perturbations are approximated as

    ℱ_(G) =1 ,   𝒢_(G)≈ν^2  a^2  >0 ,

which means these are under the complete stability.
On the other hand, those of the scalar perturbations possibly take negative values.
Therefore, we investigate the asymptotic forms of 𝒢 and ℱ in the late-time universe after bounce. 



  §.§.§ Asymptotic flat spacetime

Although it is impossible to realize a bounce solution in flat spacetime, 
it is expected to be an approximate solution inside the Hubble radius
in the late-time of the universe.
Therefore, we examine the stability of the flat spacetime 
as an asymptotic spacetime after bounce.
We assume K=0, Λ =0 and the scale factor is large, 
then, ℱ_(Q) and 𝒢_(Q) are given by<cit.>

    ℱ_(Q)   =   2 ς (3λ -1)ν^2  3ς (λ -1)ν^2 + 6(3 λ -1) H^2 a^2  , 
    𝒢_(Q)   ≈   2(2-ς)ν^2  3ς a^2 .
 
Since 𝒢_(Q) cannot be positive if ς <0, 
we focus  only on the case with ς >0.
To preserve ℱ_(Q) >0 and 𝒢_(Q)>0,
we must impose 

    λ >1  and  0< ς < 2  .
 
Then, the solutions are under the complete stability.



  §.§.§ Asymptotic Milne spacetime

The Milne universe is also a solution of open FLRW spacetime without a
cosmological constant in deep infrared regime.
The scale factor evolves as
 
    a ∼√(2  3λ-1)  t  .

 Then, ℱ_(Q) and 𝒢_(Q) are approximated by

    ℱ_(Q)   ≈   2ς (3λ-1) ν^2 (ν^2 +3)/3[ ς(λ-1)ν^4 +2(2-ς)ν^2 +12 ] , 
    
    𝒢_(Q)   ≈    2ν^2 (ν^2 +3)  3 a^2 
           ×[ (2-ς)(λ-1)ν^2 +2ς +6(λ-1) ς(λ-1)ν^4 +2(2-ς)ν^2 +12 ]  .
 
Then, the conditions for the positive ℱ_(Q) and 𝒢_(Q)
are given by

    λ≥ 1 and  0 < ς < 2  .
 
The important point is we can stabilize asymptotic spacetime after bounce
without cosmological constant unlike the projectable case.

We also analyze the asymptotic dynamics of the perturbations.
The asymptotic behaviors of the effective friction coefficients are given by

    ℋ^2_(G)   ∼   1  t^2 , 
    ℋ^2_(Q)   ∼   1  t^2  ,

and those of the effective squared masses are reduced into 

    ℳ^2_(G) a^2    ∼   ( 3λ -1 ) ν^2  2  , 
    ℳ^2_(Q) a^2    ∼     (2 - ς)( λ - 1 ) ν^2 + 2ς + 6(λ -1)    2 ς .

Then, the asymptotic dynamics of the perturbations are approximated by
the following form: 

    h    ∼   C_1  t cos[ √(ℳ^2   a^2  -1) ln t   ] 
           
    +  C_2  t sin[ √(ℳ^2   a^2  -1) ln t   ]  ,

with integration constants C_1 and C_2.
One can see that the perturbation amplitude approaches to zero.




  §.§.§ Asymptotic de Sitter spacetime

If  a positive cosmological constant is present, 
it may possible that the universe experiences an accelerating expansion 
whose expanding law is given by
 
    a ∼exp( √(2Λ 3( 3 λ-1 ))  t )  .


Once the universe enters the accelerating expansion phase, 
the spatial curvature turns to be irrelevant within a few Hubble time, 
then, the spacetime asymptotically approaches to de Sitter spacetime.
In that case, ℱ_(Q) and 𝒢_(Q) are approximated by[
It is natural to consider that
there exists an ultraviolet cut-off momentum p_cut 
at which a non-perturbative quantum effect of gravity cannot be ignored.
More precisely, let ν_cut be an ultraviolet cut-off mode 
which is related with the cut-off momentum as 
ν^2_cut/a_T^2 ≈  p^2_cut
at the early stage of the universe.
Then, we find some finite time so that ν^2/a^2 < H^2 for any ν≤ν_cut.
Therefore the all perturbation modes are rapidly suppressed 
by the strong ℋ-friction after that time.
]

    ℱ_(Q)   ≈   ςν^2  3H^2 a^2 , 
    𝒢_(Q)   ≈   2ν^2/27a^4[ 27(2+ς)/4Λ( (λ-1)ν^2+2K ) 
    
               
    + (ν^2-3K)(4g_r+3g_3) ]   .


The positivity of ℱ_(Q) is satisfied if ς >0 which is 
consistent with that of in the asymptotic flat spacetime.
We show the explicit conditions for 𝒢_(Q)≥ 0 in closed 
spacetime: 

    3g_3 + 4 g_r + 27 (λ -1)(2+ς)  4 Λ≥ 0 .

and in open spacetime: 

    3g_3 + 4 g_r + 27 (λ -2)(2+ς)  10 Λ≥ 0 .

If the above condition is satisfied, the solutions are under the complete stability 
in the asymptotic de Sitter era.

It is clear that the squared effective masses ℳ_(Q)^2
approach asymptotically to zero because it is proportional to a^-2.
Therefore the case with 𝒢_(Q)<0 seems not to be quite problematic because 
the solution is under the  ℋ-suppressed stability.
To confirm it, we examine the asymptotic dynamics of the perturbations.
 The effective friction coefficients asymptotically  behave as 

    ℋ_(G)^2    ∼    2Λ 3(3λ -1) , 
    ℋ_(Q)^2    ∼    2Λ 27(3λ -1) ,

 and those of the effective squared masses are

    ℳ^2_(G) a^2     ∼   ν^2  , 
    ℳ^2_(Q) a^2    ∼   2 H^2 /9 ς[ 27(2+ς)/4Λ( (λ-1)ν^2+2K ) 
    
               
    + (ν^2-3K)(4g_r+3g_3) ]   .



Substituting them into the perturbation equations of motion, we obtain the asymptotic behaviors.
For the case with ℳ^2 >0, which means the perturbations are under the complete stability,  
the asymptotic dynamics of the tensor perturbations are

    h_(G)   ∼   [ C_3 √(ℳ_(G)^2   a^2  H^2 ) e^- H t  + C_4  ] 
                   ×cos( √(ℳ_(G)^2   a^2  H^2 ) e^- H t) 
           - [ C_3 - C_4 √(ℳ_(G)^2   a^2  H^2 ) e^- H t] 
                   ×sin(  √(ℳ_(G)^2   a^2  H^2 ) e^- H t)  ,

and those of the scalar perturbations are given by

    h_(Q)   ∼     C_5 cos( √(ℳ_(Q)^2   a^2  H^2 ) e^- H t) 
                - C_6 sin(  √(ℳ_(Q)^2   a^2  H^2 ) e^- H t)  ,
 
 with integration constants C_3, C_4, C_5 and C_6.
 We find that the perturbations are rapidly suppressed and settled into some constants 
 as h_(G)∼ C_4 and h_(Q)∼ C_5.
 Since the constants C_4 and C_5 are expected to be a typical amplitude 
 of the perturbations when the universe enters the de Sitter era, 
 we generally observe non zero value of perturbation amplitude after entering de Sitter phase.
 
 For the case with ℳ^2 <0,  the dynamics of the tensor perturbations are approximated as

    h_(G)   ∼   [ C_7 √( |ℳ_(G)^2   a^2 |  H^2 ) e^- H t  + C_8  ]  
                ×cosh( √( |ℳ_(G)^2   a^2 |  H^2 ) e^- H t) 
           - [ C_7 - C_8 √( | ℳ_(G)^2   a^2 |  H^2 ) e^- H t] 
                ×sinh(  √(| ℳ_(G)^2   a^2 | H^2 ) e^- H t)  ,

and those of the scalar perturbations are given by 
 
    h_(Q)   ∼     C_9 cosh( √( |ℳ_(Q)^2   a^2 |  H^2 ) e^- H t) 
             -C_10sinh(  √(| ℳ_(Q)^2   a^2 | H^2 ) e^- H t)  ,

 with integration constants C_7, C_8, C_9 and C_10.
 Then, we also find that the perturbations are rapidly settled into some constant as h_(G)∼ C_8 and h_(Q)∼ C_9, 
 which means the ℋ-suppressed stability is truly stabilized state.



 §.§ Whole history of the universe

Let us turn our attention to the whole history of the universe including bouncing phase.
Although our main analysis depends on the numerical approach, 
we, in advance, consider some simple specific cases again in which explicit conditions can be derived.



  §.§.§ Tensor perturbations


Since the coefficients of kinetic terms in tensor perturbations are unity, 
these are free from the ghost instability.  
Therefore the stabilities are determined only by 𝒢_(G).
As we mentioned, the forms of ℱ_(G) and 𝒢_(G) in tensor perturbation
are identical to those of the projectable case.
Thus, we summarize the positivity conditions of 𝒢_(G), 
which we have argued in previous paper<cit.>.

To investigate the stability of the tensor perturbations, 
we focus on the case with large ν^2.
Then, 𝒢_(G) is approximated as 

    𝒢_(G)∼ν^2  a^2  + g_3 ν^4  a^4 + g_8 ν^6  a^6 .

One can see that the ultraviolet stability imposes g_8 ≥ 0. 

Although the explicit condition for general case is quite complicated, 
we show the conditions with g_3 = 0.
In closed spacetime (K=1), if one of the following conditions is satisfied, 
𝒢_(G) must be positive for any a>0 and ν^2 ≥ 8: 

    (i)        0<g_8 ≤ -3g_56/7 , g_r < 0 , g_s ≤  -(g_8-3g_56)^2  g_8 ; 
    
    (ii)        0<g_8 ≤ -3g_56/7 , g_r ≥ 0 , 
               g_s ≤  -(g_8-3g_56)^2  g_8 - g_r^2  9 ; 
    
    (iii)        0≤ g_8 ,  -3g_56/7 < g_8 , g_r < 0 , g_s ≤ 48(g_56 + g_8)  ; 
    
    (iv)         0≤ g_8 ,  -3g_56/7 < g_8 , g_r ≥ 0 , 
               g_s ≤ 48(g_56 + g_8) - g_r^2  9 ;
 
Similarly, in open spacetime (K=-1): 

    (i)        0 < g_8 ≤ g_56 , g_r >0   , g_s ≤ -(g_8-3g_56)^2  g_8 ; 
    
    (ii)        0 < g_8 ≤ g_56 , g_r ≤ 0  , g_s ≤ -(g_8-3g_56)^2  g_8 -g_r^2  9 ; 
    
    (iii)        0≤ g_8 , g_56 < g_8 , g_r >0  , g_s ≤ 4(2g_8 -3g_56)  ; 
    
    (iv)        0≤ g_8 , g_56 < g_8 , g_r ≤ 0  , g_s ≤ 4(2g_8 -3g_56) - g_r^2  9 ,

which are equivalent to 𝒢_(G)≥ 0 for any a>0 and ν^2 ≥ 2. 
It is worth mentioning that these stability conditions in open FLRW spacetime are not completely 
conflict, however slightly difficult to be compatible with the bouncing conditions
(see <cit.>).



  §.§.§ Scalar perturbations

The ghost-free condition for scalar perturbations can be analytically discussed.
 Let the coupling constants be λ >1/3 and ς >0, and  adopting ν^2 -3K >0,
the condition ℱ_(Q) >0 is equivalent to 

    ςν^2 [ (λ-1)ν^2 +2K ]  2(3λ-1)(ν^2 -3K)a^2  > -H^2  .
 
The most stringent conditions are imposed at the bouncing and the recollapse points, i.e., H=0:

    ς >0 and λ > 1  for K=0 
    λ≥ 1  for K=1 
    λ > 2  for K=-1
     .

If the above conditions are satisfied, the ghost instability can be eliminated 
throughout the evolution of the universe.

On the other hand, the positivity conditions of 𝒢_(Q) 
cannot be expressed without any approximation due to the extremely complicated form.
Therefore, we restrict our analysis to the case with large ν^2 modes at this stage. 
Then, 𝒢_(Q) is given by

    𝒢_(Q)≈[ (-8g_7+3g_8) + 72g_7^2K^2/ς a^4] 2ν^6/3a^6 ,

Since ς should be positive,  we assume 

    3g_8≥8g_7  .

Then, 𝒢_(Q)≥ 0 for large ν is guaranteed.



  §.§.§ Numerical analysis





 
Before performing the numerical analysis, we summarize the stability conditions 
for each cases.


  (1) The cases with K=1 and Λ≤ 0: 
Since there exists a finite maximum value of scale factor a_max,
i.e., the solutions 𝒪^[0; 1] and 𝒪^[-1; 1], 
we assume the stability conditions around the bouncing and the recollapsing points  
and in the large ν^2 region, thus, 

    ς >0 , λ≥ 1 ,  g_8≥ 0 and 3g_8 ≥ 8g_7  .
  

  (2) The cases with K=1 and Λ >0:
In this cases, we observe ℬ^[1;1], ℬ^[1;1]_BC
 and ℬ^[1;1]_O as the singularity-free solutions.  
To guarantee the stabilities around bouncing point (and the recollapsing point if ℬ^[1;1]_O), 
in asymptotic de Sitter spacetime 
and in the large ν^2 region, we assume

    ς >0 , λ≥ 1 , 
    3g_3 +4g_r + 27 (λ -1) (2+ς) 4 Λ≥ 0 , 
          g_8≥0 and 3g_8 ≥ 8g_7  .

If we focus only on the oscillating phase which is observed in the solution ℬ^[1;1]_O, 
the stability conditions in the asymptotic de Sitter spacetime  
can be relaxed. Then, we impose the same conditions as
the cases with K=1 and Λ≤ 0, that is, (<ref>).

  (3) The cases with K=-1 and Λ = 0: 
We find ℬ^[0,-1] and ℬ^[0,-1]_BC
as the singularity-free solutions whose asymptotic behaviors are the Milne spacetime.
Thus, we impose 

    0 < ς <2 , λ > 2 , g_8≥0 and 3g_8 ≥ 7 g_7 .


  (4) The cases with K=-1 and Λ > 0: 
In this case, the solutions ℬ^[1,-1] and ℬ^[1,-1]_BC
are realized.
Since both of the solutions approaches the de Sitter spacetime after the bounce, 
we impose the conditions as follows:  
 
    ς >0 , λ > 2 , 
    3g_3 +4g_r + 27 (λ -2) (2+ς) 10 Λ≥ 0 , 
          g_8≥0 and 3g_8 ≥ 8g_7  .


  (5) The cases with K=-1 and Λ < 0: 
The possible singularity-free solutions are 𝒪^[-1; -1] and 𝒪^[-1; -1]_BC.
Therefore we impose the stability conditions around the bouncing and recollapsing points and 
in the large ν^2 region:
 
    ς >0 , λ > 2 , g_8≥0 and 3g_8 ≥ 8g_7  .




Then, we shall investigate the spacetime stabilities in each cases: 
bouncing solutions with asymptotic de Sitter spacetime, 
bouncing solutions with asymptotic Milne spacetime and
oscillating solutions.
In our numerical analysis, we have found stable  singularity-free solutions 
throughout the whole evolutions, which means 
these solutions are under the  ℳ^2-dominated stability  when ℋ<0
and under the complete stability  when ℋ>0.
The concrete examples are shown in TABLE <ref>.





  §.§.§ Bouncing solutions with asymptotic Milne spacetime
 
The solutions ℬ^[0;-1] and ℬ_BC^[0;-1] correspond to this case. 
Recall that these types of singularity-free solutions cannot be stabilized in the projectable HL theory.
Since the Hubble parameter is dropped as t^-1∼ a^-1, 
the effect of the squared effective masses do not suppressed by the Hubble friction, 
that is, the values of ℳ^2/H^2 are negative constants and do not 
approach to zero at asymptotic Milne regime.

On the other hand, we can discover the stable solutions in the non-projectable HL theory 
because it is possible to keep the values of ℳ^2_(Q) to be positive during the whole evolution.
In FIG. <ref>, a typical example of such a solution is shown.

The figure shows that ℳ^2/ℋ^2 of both tensor and scalar perturbations keep the positive values
and always greater than unity.
Thus, the  ℋ-accelerated instability does not occur for any initial value of the scale factor a_ini.
As far as our numerical analysis is concerned, such behavior can also be seen for any possible perturbation mode ν.
Therefore we regard our example as a solution under complete stability  when ℋ>0 
and under  ℳ^2-dominated stability  when ℋ<0 for any a_ini.

One may wonder whether the squared effective masses truly suppress 
the effects of  ℋ-acceleration
under the  ℳ^2-dominated stability, because we do not mention any explicit criterions
of the ratio ℳ^2/ℋ^2 so far.
Thus, we additionally discuss the dynamical evolutions of each perturbation modes
to investigate the growth of the perturbation amplitudes caused by the  ℋ-acceleration. 

To clarify the effects, we numerically solve the perturbation equations of motion
(<ref>). 
The initial conditions are given at a_ini=100 a_T so that
ḣ_(G) =0 = ḣ_(Q). 
Since the equations of motion are linear with respect to h_(G) and 
h_(Q), 
the whole evolutions are proportional to the values of the initial conditions  h_(G) ini and h_(Q) ini.
Therefore, we trace the ratios to each initial values, i.e., h_(G) / h_(G) ini and h_(Q) / h_(Q)ini. 
In FIG <ref>, the time evolutions are shown.

 
This result is based on the solution ℬ^[0; -1] whose coupling constants
are given in TABLE <ref> (ii), thus,
also corresponds with FIG <ref>.
From the figure, we find that the perturbation amplitudes are enhanced around the bouncing point
even if the solution is under the  ℳ^2-dominated stability.
Thus, we should clarify the growth rates of the perturbations.
In TABLE <ref>, we show the detail data of the maximum 
ratios h_(G) / h_(G) ini and h_(Q) / h_(Q)ini 
throughout the evolutions  for each perturbation modes.
From the data, the amplitude possibly grows up to about 50 times and 
tends to be higher as the perturbation mode ν^2 is lower.
Since the large perturbation mode ν^2 corresponds to
the heavy effective squared mass ℳ^2, 
the result seems appropriate. 

 


  §.§.§ Bouncing solutions with asymptotic de Sitter spacetime
 
The solutions ℬ^[1;1], ℬ^[1;1]_BC, ℬ^[1;1]_O
with a ≥ a_T, ℬ^[1;-1] and ℬ^[1;-1]_BC correspond to this case. 
Actually, the previous work based on the projectable HL theory has shown that we cannot construct a solution 
under the complete stability. 
Instead of this, it is possible to find solutions under the ℋ-suppressed stability
due to the negative ℳ_(Q)^2 in infrared regime.

On the other hand, it turns out that solutions under the complete stability can be realized in the non-projectable case.
As an example, we show the evolutions of the ℳ^2/ℋ^2 in FIG. <ref>.

As can be seen from the figure, both tensor and scalar perturbations keep the positive values during whole evolution
of the solution. 
As far as our analysis is concerned, the squared effective masses of tensor and scalar perturbations are monotonically 
increasing functions with respect to the perturbation mode ν^2.
Therefore we conclude this solution is under complete stability for any possible ν^2.

The magnitude relationships between ℳ^2 and ℋ^2 are also analyzed.
As we discussed, both ℳ^2_(G) and ℳ^2_(Q) drop as a^-2, 
in contrast, the  effective friction coefficients are settled to a constant.
Thus, there must exist a certain value of the scale factor such that any one of ℳ^2 
is equal to ℋ^2.
We define such a value of scale factor as a critical scale factor a_crit.
In our numerical example FIG. <ref>, the critical scale factor is a_crit=5.118
at which the lowest ν^2 mode of the tensor perturbation shows
 ℳ_(G)^2=ℋ^2_(G)
and ℳ_(G)^2 is weaker than ℋ_(G)^2 for a>a_crit.
Thus, in order to prevent ℋ-accelerated instability, the initial value of the scale factor a_ini needs to be 
smaller than that critical value.


We also mention the growth of the perturbations around the bouncing point.
Actually, it is found that the behavior is not quite different from the previous case with asymptotic Milne spacetime.
It is reasonable because the effect of a cosmological constant is expected to be relatively weaker than those of
higher order curvature terms.



  §.§.§ Oscillating solutions
 
We examine the stabilities of the oscillating solutions, specifically, 𝒪^[0;1], ℬ^[1;1]_O for 
a_min≤ a ≤ a_max, 
𝒪^[-1;1], 𝒪^[-1;-1] and 𝒪^[-1;-1]_BC for 
a_min≤ a ≤ a_max. 
In this case, we also construct solutions whose squared effective masses are positive and 
dominate the ℋ-term  within the possible ranges of the scale factor. 
A typical example is shown in FIG. <ref>.

In our analysis, ℳ^2 of all perturbation modes are positive and greater than ℋ^2.
It means that the solution is under the ℳ^2-dominated stability  when ℋ<0
and under the complete stability  when ℋ>0.
It is remarkable that, in that case, the scalar perturbations with low ν^2 show Hℋ_(Q) <0 
around the bouncing point, which means that the scalar perturbations feel ℋ-friction slightly before the bounce and
receive ℋ-acceleration slightly after the bounce.

Additionally, we show the dynamics of the tensor and scalar perturbations 
in the oscillating solution in FIG. <ref>.

From these figures, the typical amplitudes of the perturbations are almost constant 
even if the universe experiences a number of oscillations.
Thus, we conclude the solution is stable with respect to the linear perturbations.








§ SUMMARY AND DISCUSSIONS


We have investigated the stability of the singularity-free cosmological solutions based on the non-projectable HL theory
whose action is given by (<ref>).
Since our aim is to remedy the infrared behaviors of the singularity-free solutions in the projectable HL theory,
we introduced only the single additional term Φ^2 which is expected to be dominant in infrared limit.
It is remarkable that our gravitational action realizes the identical background solutions in FLRW spacetime
based on the projectable HL theory with 𝒞=0.
Therefore, the bouncing solutions and the oscillating solutions which are induced by the higher order 
spatial derivative terms in the action
are also found as in the projectable case.

By considering the quadratic action, we discuss the stability of the singularity-free
solutions with respect to tensor and scalar perturbations.
The stabilities can be estimated by the sign of the coefficients ℱ and 𝒢 showed in (<ref>) and (<ref>). 
ℱ<0 corresponds to the ghost instability which is equivalent to a lack of the lowest energy state.
The case with 𝒢<0 is known as the gradient instability which is interpreted as
a unstable behavior due to the negative squared effective mass.
However, it is possible to consider the case that the unstable behavior caused by 𝒢<0 is suppressed by the effect of background dynamics.
As we showed in (<ref>), the stability of the perturbations can be
judged by the sign of the effective friction coefficients ℋ and 
the effective squared masses ℳ^2, additionally 
the magnitude relationships between ℋ^2 and |ℳ^2|.
Thus, we introduced five types of stabilities and instabilities 
(i) negative-ℳ^2 instability, (ii) ℋ-accelerated instability, (iii) ℳ^2-dominated stability, 
(iv) ℋ-suppressed stability, and (v) complete stability.

The novel feature of the singularity-free cosmological solutions in the non-projectable HL theory
is that we can find the bouncing solutions which satisfy the 
ℳ^2-dominated stability condition when ℋ<0
and the complete stability condition when ℋ>0.
Such solutions cannot be constructed in the projectable HL theory,
that is, the squared effective masses must be negative in infrared region. 

We additionally investigate the stability of the oscillating solutions.
The solution (x) we have shown in Sec. <ref> is 
the stable solution which satisfy the ℳ^2-dominated stability condition when ℋ<0
and the complete stability condition when ℋ>0.
In fact, the typical amplitudes of the perturbations stay almost constant. 


However, we would like to indicate that it is not impossible to construct 
an oscillating solution whose typical perturbation amplitude is exponentially enhanced 
even if either ℳ^2-dominated stability condition or 
the complete stability condition is satisfied.
That instability is caused by a resonance.
Whether the resonance is induced or not can be investigated 
by the Hill's method<cit.> which is summarized in Appendix <ref>.



We numerically show a example of the oscillating solution with resonance 
instability in FIG. <ref>
and the coupling constants of the solution is shown in TABLE <ref>.
Note that this solution is satisfied both ℳ^2-stability condition 
and complete stability condition (see FIG. <ref>).
In that case, certain modes of the tensor perturbation show unstable behavior.
The Hill's method indicates that the degree of enhancement is characterized by ϵ_± defined in (<ref>).
If |ϵ_±| exceeds 1, the amplitude of the corresponding perturbation mode exponentially grows.
In fact, the solution includes certain tensor perturbation modes with |ϵ_±|>1 
as we show in FIG.<ref>.
It shows that even if all stability conditions we have discussed are satisfied, 
the oscillating solution is possibly unstable due to the resonance.

 

The oscillating solutions with resonance instability can be discovered 
if we set the coupling constant g_8 to be positive, however, relatively smaller than
that of the stable solutions listed in TABLE <ref>
as far as our numerical analysis is concerned.
One may notice that such a manipulation corresponds to consider the small effective squared mass.
See (<ref>) and (<ref>), 
it is found that ultraviolet dominant terms g_8 ν^6/a^6 decrease.
Thus, it is natural to consider that the heavy effective squared masses prevent 
the oscillating solution from the resonance instability. 
Actually, we can see from FIG.<ref> that the degree of enhancement 
|ϵ_±| approaches unity as the perturbation mode ν^2 increases. 
We would like to stress that the resonance instability is mainly the problem in open FLRW spacetime.
Since the perturbation mode ν^2 takes discrete number greater than or equal to eight in closed FLRW spacetime, 
the resonance instability is not quite problematic.
In general terms, the positive heavy squared effective masses 
ℳ^2 is preferred to avoid the instabilities due to the background dynamics, 
i.e., ℋ-accelerated instability and the resonance instability.
Thus, the higher order curvature terms with z=3 are significant to stabilize the spacetime 
as well as inducing the bounce.


 

 

We also mention the stability of singularity-free solutions 
based on the general covariant HL theory
which is an alternative modification  to remedy the behavior of scalar perturbation<cit.>.
In this theory, the scalar propagating degree of freedom is eliminated by 
the additional fields, 
whereas the tensor perturbations are not modified.
It is notable that the identical background solutions in FLRW spacetime to 
those of projectable and non-projectable HL theory can be reproduced 
under certain conditions<cit.>.
Therefore the stability analyses we have performed are also valid 
in the general covariant HL theory.
For example, the background solutions with g_3 =0 hold positive squared effective
masses if the conditions we have explicitly shown in Sec. <ref> are satisfied.
Then, we can construct the stable singularity-free solutions 
only by examining ℳ^2_(G)/H^2.

The perturbative approach we have performed is based on a postulate.
The effect of the backreaction can be ignored.
In other words,  the background dynamics are never affected by the perturbations.
Of course, that assumption is not trivial.
If the perturbation amplitudes are much enhanced, it is possible that the perturbation fields 
impinge on the background dynamics.
Our numerical analysis shows that the typical amplitudes of the perturbations are enhanced around bouncing point.
The degree of enhancement tends to be larger if | ℳ^2 | is small.
Therefore one may wonder how much the perturbation amplitude is allowed not to affect to the background dynamics.
For certain perturbation mode, we can give an estimation. 
In fact, the tensor perturbations with ν^2 =8 in closed FLRW spacetime include the homogeneous modes, 
which corresponds the Bianchi-type IX spacetime whose anisotropy is small (in detail,  see <cit.>).
Thus it might be possible to estimate the effect of backreaction by considering more general background 
including inhomogeneity and anisotropy.






§ ACKNOWLEDGMENTS


The authors would like to thank K. Maeda and K. Aoki for 
their insightful comments and discussions.
MF and SS are grateful to the Early Bird Program from Waseda Research Institute for Science and Engineering, 
Grant-in-Aid for Young Scientists.
The work of S. M. was supported in a part by Grants-in-Aid from the Scientific Research Fund of the Japan Society 
for the Promotion of Science (No.18J11983).














§ RESONANCE IN THE OSCILLATING SOLUTION



In the oscillating solutions, the perturbations feel periodical external forces, 
which may causes a resonance and the rapid growth of the perturbations.
To clarify whether the resonance is induced or not, we examine the perturbation equations by 
adopting Hill's method.
Since the procedures are completely the same between tensor and scalar perturbations,
the subscripts which represent tensor and scalar modes are abbreviated in what follows.

To rewrite the equations into the Hill's form, we consider a field redefinition as

    ĥ := √( a^3 ℱ)  h.

Then,  the perturbation equations are transformed as 

    d^2ĥ/dt^2   +ω^2(t)ĥ=0  ,

with

    ω^2(t):=1/4[ℱ̇^2/ℱ^2-2(3Hℱ̇+ℱ̈)/ℱ
    -9H^2 -6Ḣ + 4ℳ^2  ]  ,

Since the effective squared angular frequency ω^2(t) is essentially a function which depend only on the scale factor, 
the oscillating period is expected to be the same as that of the scale factor.
We further introduce the following equation in matrix expression as

    d/dtX(t)
        =   (
        [       0       1; -ω^2(t)       0 ])
      X(t) , 
    
    
    X(t)    :=   [  ĥ(t); ĥ̇(t) ] ,

Here  ĥ(t) is an arbitrary perturbation function and X(t) is a real solution vector for the Hill equation (<ref>).

Let T be a period which is characterized by the effective squared angular momentum ω^2, 
then, the evolution of X is expressed as 

    X(t+T)=𝒜 X(t)  ,

where, 𝒜 is a 2 × 2 matrix which represents the time translation t → t + T.
Suppose {ϵ_+, X_ϵ_+} and {ϵ_-, X_ϵ_-}  are two independent 
eigensystems of matrix 𝒜, 
where, ϵ_± are eigenvalues  and X_ϵ_± are corresponding eigenvectors.
Then, (<ref>) can be rewritten as some linear combination of the following relations: 

    X_ϵ_±  (t+T) = ϵ_± X_ϵ_± (t)  .

 Therefore, we find that both |ϵ_+| ≤ 1 and |ϵ_-| ≤ 1 are required 
 to suppress the resonance of the perturbations.

 To specify the explicit conditions for |ϵ_±| ≤ 1, 
 we introduce two independent complex vector X_1 and X_2
 whose initial conditions are

    X_1(0)   = 
      (
        [  ĥ_1(0); ĥ̇_1(0) ])
      = (
        [ 1; 0 ]) , 
    
    X_2(0)   =
       (
        [  ĥ_2(0); ĥ̇_2(0) ])
    =  (
        [ 0; 1 ]).

Any solution vector can be constructed by the linear combination of X_1 and X_2.
Substituting (<ref>) and (<ref>) into (<ref>), we find

    𝒜=(
        [  ĥ_1(T)  ĥ_2(T); ĥ̇_1(T) ĥ̇_2(T) ]).

 Solving the Hill equation (<ref>) with these initial conditions, 𝒜 can be estimated.
The eigenvalues of 𝒜 are obtained from the following equation.

    ϵ^2 -(tr𝒜) ϵ +det𝒜=0.

Since the Hill equation (<ref>) and the initial conditions (<ref>) and (<ref>) give det𝒜=1,
the eigenvalues ϵ_± are given by

    ϵ_±=tr𝒜±√(( tr𝒜)^2-4)/2.

One can see that | tr𝒜 | ≤ 2 is equivalent to |ϵ_±| = 1.
When ϵ_± is a complex number, the perturbations correspond to the real part of X(t).
For the case with | tr𝒜 | > 2, any one of |ϵ_±| is greater than unity.
Thus, we conclude that if | tr𝒜 | ≤ 2, the perturbation keeps oscillating forever without growing or decaying.
On the other hand,  the perturbation shows an exponential instability if | tr𝒜 | > 2.







99

 singularity_theorem
 R. Penrose, Phys. Rev. Lett. 14, 57 (1965); 
 S. W. Hawking, Proc. R. Soc., A 300, 187 (1967); 
 S. W. Hawking and R. Penrose, Proc. R. Soc., A 314, 529 (1970); 
 S. W. Hawking and G. F. R. Ellis, The Large Scale Structure of Space-Time (Cambridge Univ., Cambridge, 1973).

inflation_not_complete
A. Borde and A. Vilenkin,
 “Eternal inflation and the initial singularity",
Phys. Rev. Lett. 72, 3305 (1994);
A. Borde, A. H. Guth and A. Vilenkin,
 “Inflationary Spacetimes Are Incomplete in Past Directions",
Phys. Rev. Lett. 90, 151301 (2003) [arXiv:0110012[gr-qc]].

braneworld_bounce
J. Khoury, B. A. Ovrut, P. J. Steinhardt and N. Turok,
“Ekpyrotic universe: Colliding branes and the origin of the hot big bang", 
Phys. Rev. D 64, 123522 (2001) [arXiv:0103239[hep-th]];
J. Khoury, B. A. Ovrut, N. Seiberg, P. J. Steinhardt and N. Turok,
“From big crunch to big bang", 
Phys. Rev. D 65, 086007 (2002) [arXiv:0108187[hep-th]].

string_gas
R. Brandenberger and C. Vafa,
“Superstrings in the early universe", 
Nucl. Phys. B 316, 391 (1989);
A. Nayeri, R. Brandenberger and C. Vafa,
“Producing a Scale-Invariant Spectrum of Perturbations in a Hagedorn Phase of String Cosmology", 
Phys. Rev. Lett. 97, 021302 (2006) [arXiv:0511140[hep-th]];
R. Brandenberger, R. Costa, G. Franzmann and A. Weltman,
“Point particle motion in double field theory and a singularity-free cosmological solution", 
Phys. Rev. D 97, 063530 (2018) [arXiv:1710.02412[hep-th]].

LQG_bounce1
D. J. Mulryne, R. Tavakol, J. E. Lidsey and G. F. R. Ellis,
“An emergent universe from a loop", 
Phys. Rev. D 71, 123512 (2005) [arXiv:0502589[astro-ph]].
LQG_bounce2
I. Agullo, A. Ashtekar and W. Nelson,
“The pre-inflationary dynamics of loop quantum cosmology: confronting quantum gravity with observations", 
Class. Quant. Grav. 30, 085014 (2013) [arXiv:1302.0254[gr-qc]].
LQG_bounce3
E. W. Ewing,
“Ekpyrotic loop quantum cosmology", 
JCAP 1308, 015 (2013) [arXiv:1306.6582[gr-qc]].


genesis1
P. Creminelli, A. Nicolis and E. Trincherini,
“Galilean genesis: an alternative to inflation", 
JCAP 1011, 021 (2010) [arXiv:1007.0027[hep-th]].
genesis2
D. A. Easson, I. Sawicki, and A. Vikman,
“G-Bounce", 
JCAP 1111, 021 (2011) [arXiv:1109.1047[hep-th]].
genesis3
S. Nishi and T. Kobayashi, 
“Generalized galilean genesis", 
JCAP 1503, 057 (2015) [arXiv:1501.02553[hep-th]].
genesis4
T. Kobayashi, M. Yamaguchi and J. Yokoyama,
“Galilean creation of the inflationary universe", 
JCAP 1507, 017 (2015) [arXiv:1504.05710[hep-th]].

Galileon_bounce
A. Ijjas and P. J. Steinhardt,
“Classically Stable Nonsingular Cosmological Bounces", 
Phys. Rev. Lett. 117, 121304 (2016) [arXiv:1606.08880[gr-qc]];
A. Ijjas and P. J. Steinhardt,
“Fully stable cosmological solutions with a non-singular classical bounce", 
Phys. Lett. B 764, 289 (2017) [arXiv:1609.01253[gr-qc]].

Galileon_bounce2
A. Ijjas,
“Space-time slicing in Horndeski theories and its implications for non-singular bouncing solutions",
JCAP 1802, 007 (2018) [arXiv:1710.05990[gr-qc]].

Galileon_bounce3
D. A. Dobre, A. V. Frolov, J. T. G. Ghersi, S. Ramazanov, and A. Vikman,
“Unbraiding the Bounce: Superluminality around the Corner",
JCAP 1803, 020 (2018) [arXiv:1712.10272[gr-qc]].

non-local_bounce
 T. Biswas, E. Gerwick, T. Koivisto and A. Mazumdar, 
 “Towards Singularity and Ghost-free Theories of Gravity"
 Phys. Rev. Lett. 108, 031101 (2012) [arXiv:1110.5249[gr-qc]];
 T. Biswas, A. S. Koshelev, A. Mazumdar and S. Y. Vernov,
“Stable bounce and inflation in non-local higher derivative cosmology"
JCAP 08, 024 (2012) [arXiv:1206.6374[astro-ph.CO]];
Á. de la Cruz-Dombriz, F. J. M. Torralba, and A. Mazumdar,
“Non-singular and ghost-free infinite derivative gravity with torsion"
arXiv:1812.04037[gr-qc].


no-go_Horndeski
T. Kobayashi,
“Generic instabilities of non-singular cosmologies in Horndeski theory: a no-go theorem",
Phys. Rev. D 94, 043511 (2016) [arXiv:1606.05831[hep-th]].

no-go_multi
S. Akama and T. Kobayashi,
“Generalized multi-Galileons, covariantized new terms, and the no-go theorem for non-singular cosmologies"
Phys. Rev. D 95, 064011 (2017) [arXiv:1701.02926[hep-th]]

no-go_non-flat
S. Akama and T. Kobayashi,
“General theory of cosmological perturbations in open and closed universes from the Horndeski action",
arXiv:1810.01863[gr-qc].

no-go_papers
R. Kolevatov, and S. Mironov,
“Cosmological bounces and Lorentzian wormholes in Galileon theories with an extra scalar field",
Phys. Rev. D 94, 123516 (2016) [arXiv:1607.04099[hep-th]].

Horndeski_theory
G. W. Horndeski, 
“Second-order scalar-tensor field equations in a four-dimensional space", 
Int. J. Theor. Phys. 10, 363 (1974); 
T. Kobayashi, M. Yamaguchi and J. Yokoyama, 
“Generalized G-inflation: Inflation with the most general second-order field equations", 
Prog. Theor. Phys. 126, 511 (2011) [arXiv:1105.5723 [hep-th]].

bounce_EFT
Y. Cai, Y. Wan, H-G. Li, T. Qiu and Y-S. Piao, 
“The Effective Field Theory of nonsingular cosmology", 
JHEP 1701, 090 (2017) [arXiv:1610.03400[gr-qc]];
Y. Cai, Y. Wan, H-G. Li, T. Qiu and Y-S. Piao, 
“The Effective Field Theory of nonsingular cosmology:II", 
Eur. Phys. J. C 77, 369 (2017) [arXiv:1701.04330[gr-qc]].


no-go_loopholes1
P. Creminelli, D. Pirtskhalava, L. Santoni, and E. Trincherini,
“Stability of Geodesically Complete Cosmologies",
JCAP 1611, 047 (2016) [arXiv:1610.04207[hep-th]].

no-go_loopholes2
S. Banerjee, and E. N. Saridakis,
“Bounce and cyclic cosmology in weakly broken galileon theories",
Phys. Rev. D 95, 063523 (2017) [arXiv:1604.06932[gr-qc]].

no-go_loopholes3
R. Kolevatov, S. Mironov, N. Sukhov, and V. Volkova, 
“Cosmological bounce and Genesis beyond Horndeski",
JCAP 1708, 038 (2017) [arXiv:1705.06626[hep-th]].

Horava
P. Hořava, 
“Quantum gravity at a Lifshitz point", 
Phys. Rev. D 79, 084008 (2009) [arXiv:0901.3775[hep-th]].


HL-review2
A. Wang, 
“Hořava Gravity at a Lifshitz Point: A Progress Report", 
Int. J. Mod. Phys. D 26, 1730014 (2017) [arXiv:1701.06087 [gr-qc]].

 
 power-counting_renormalizability
T. Fujimori, T. Inami, K. Izumi and T. Kitamura,
“Power-counting and renormalizability in Lifshitz scalar theory", 
Phys. Rev. D 91, 125007 (2015) [arXiv:1502.01820[hep-th]];
T. Fujimori, T. Inami, K. Izumi and T. Kitamura,
“Tree-Level Unitarity and Renormalizability in Lifshitz Scalar Theory", 
Prog. Theor. Exp. Phys. 2016, 013B08 (2016) [arXiv:1510.07237[hep-th]].



ref_renormalize_HL
 A. O. Barvinsky, D. Blas, M. Herrero-Valea, S. M. Sibiryakov and  C. F. Steinwachs,
 “Renormalization of Horava gravity", 
 Phys. Rev. D 93, 064022 (2016) [arXiv:1512.02250[hep-th]].


HL-review1
S. Mukhoyama,
“Hořava-Lifshitz cosmology: a review", 
Class. Quant. Grav. 27, 223101 (2010) [arXiv:1007.5199[hep-th]].



CP1
S. Mukhoyama,
 “Scale-invariant cosmological perturbations from Hořava-Lifshitz gravity without inflation", 
JCAP 090, 001 (2009) [arXiv:0904.2190[hep-th]].

CP2
 R. Cai, B. Hu and H. Zhang, 
 “Dynamical Scalar Degree of Freedom in Horava-Lifshitz Gravity", 
Phys. Rev. D 80, 041501 (2009) [arXiv:0905.0255[hep-th]].

CP3
T. Kobayashi, Y. Urakawa and M. Yamaguchi, 
“Large scale evolution of the curvature perturbation in Horava-Lifshitz cosmology", 
JCAP 0911, 015 (2009) [arXiv:0908.1005[astro-ph.CO]].

CP4
A. Wang, D. Wands and R. Maartens,
 “Scalar field perturbations in Hořava-Lifshitz cosmology", 
JCAP 1003 013 (2010) [arXiv:0909.5167[hep-th]];
A. Wang,
 “Vector and tensor perturbations in Horava-Lifshitz cosmology", 
 Phys. Rev. D 82, 124063 (2010) [arXiv:1008.3637[hep-th]].

CP5
K. Izumi, T. Kobayashi and S. Mukohyama, 
“Non-Gaussianity from Lifshitz Scalar", 
JCAP 1010, 031 (2010) [arXiv:1008.1406 [hep-th]].


HLC1
T. Takahashi, and J. Soda,
 “Chiral Primordial Gravitational Waves from a Lifshitz Point", 
Phys. Rev. Lett. 102, 231301 (2009) [arXiv:0904.0554[hep-th]].


GW-HL
A. E. Gümrükçüoğlu, M. Saravani and T. P. Sotiriou,
“Hořava Gravity after GW170817", 
Phys. Rev. D 97, 024032 (2018) [arXiv:1711.08845[gr-qc]].



HLC2
S. Mukohyama, K. Nakayama, F. Takahashi, S. Yokoyama,
 “Phenomenological aspects of Hořava-Lifshitz cosmology", 
Phys .Lett. B 679, 6-9 (2009) [arXiv:0905.0055[hep-th]].


 DM_integrate
 S. Mukohyama, 
 “Dark matter as integration constant in Horava-Lifshitz gravity", 
 Phys. Rev. D 80, 064005 (2009) [arXiv:0905.3563[hep-th]].

HLC3
A. Wang and Y. Wu, 
“Thermodynamics and classification of cosmological models in the Horava-Lifshitz theory of gravity"
JCAP 0907 012 (2009) [arXiv:0905.4117[hep-th]].

dS_projectable
Y. Q. Huang, A. Wang and Q. Wu, 
“Stability of the de Sitter spacetime in Horava-Lifshitz theory",
Mod. Phys. Lett. A25 2267 (2010) [arXiv:1003.2003[hep-th]];
A. Wang and Q. Wu,
“Stability of spin-0 graviton and strong coupling in Horava-Lifshitz theory of gravity",
Phys. Rev. D 83 044025 (2011) [arXiv:1009.0268[hep-th]].

bounce_HL1
G. Calcagni,
 “Cosmology of the Lifshitz universe", 
JHEP 0909, 112 (2009) [arXiv:0904.0829[hep-th]].
bounce_HL2 
E. Kiritsis and G. Kofinas, 
“Horava-Lifshitz Cosmology", 
Nucl. Phys. B821: 467-480 (2009) [arXiv:0904.1334[hep-th]].
bounce_HL3
R. H. Brandenberger, 
“Matter Bounce in Horava-Lifshitz Cosmology", 
Phys. Rev. D80, 043516 (2009) [arXiv:0904.2835[hep-th]].
bounce_HL4
Y. F. Cai, and E. N. Saridakis,
“Non-singular cosmology in a model of non-relativistic gravity ", 
JCAP 0910 020 (2009) [arXiv:0906.1789[hep-th]].

previous_oscillating_universe
K. Maeda, Y. Misonoh and T. Kobayashi, 
“Oscillating Universe in Horava-Lifshitz Gravity"
Phys. Rev. D 82, 064024 (2010) [arXiv:1006.2739[hep-th]].

Bianchi_IX
Y. Misonoh, K. Maeda and T. Kobayashi, 
“Oscillating Bianchi IX Universe in Horava-Lifshitz Gravity"
Phys. Rev. D 84, 064030 (2011) [arXiv:1104.3978[hep-th]].

previous_stability_bounce
Y. Misonoh, M. Fukushima, and S. Miyashita, 
“Stability of singularity-free cosmological solutions in Hořava-Lifshitz gravity", 
Phys. Rev. D 95, 044044 (2017) [arXiv:1612.09077[gr-qc]].

problem_projectable1
C. Charmousis, G. Niz, A. Padilla and P. M. Saffin,
 “Strong coupling in Hořava gravity", 
JHEP 0908, 070 (2009) [arXiv:0905.2579[hep-th]].

problem_projectable2
M. Li and Y. Pang,
“A trouble with Hořava-Lifshitz gravity", 
JHEP 0908, 015 (2009) [arXiv:0905.2751[hep-th]].

problem_projectable3
D. Blas, O. Pujolas and S. Sibiryakov,
“On the extra mode and inconsistency of Hořava gravity", 
JHEP 0910, 029 (2009) [arXiv:0906.3046[hep-th]].

problem_projectable4
K. Koyama and F. Arroja,
 “Pathological behaviour of the scalar graviton in Hořava-Lifshitz gravity", 
JHEP 1003, 061 (2010) [arXiv:0910.1998[hep-th]].

healthy_extension
D. Blas, O. Pujolas, S. Sibiryakov, 
“Consistent Extension of Hořava Gravity", 
Phys. Rev. Lett. 104, 181302 (2010) [arXiv:0909.3525[hep-th]];
D. Blas, O. Pujolas, S. Sibiryakov, 
“Comment on “Strong coupling in extended Hořava-Lifshitz gravity” [Phys. Lett. B 685 (2010) 197]", 
Phys. Lett. B 688, 350 (2010) [arXiv:0912.0550[hep-th]];
D. Blas, O. Pujolas, S. Sibiryakov, 
“Models of non-relativistic quantum gravity: the good, the bad and the healthy", 
JHEP 1104, 018 (2011) [arXiv:1007.3503[hep-th]].

NP_stability1
T. Kobayashi, Y. Urakawa, and M. Yamaguchi,
“Cosmological perturbations in a healthy extension of Hořava gravity", 
JCAP 1004, 025 (2010) [arXiv:1002.3101[hep-th]].

NP_stability2
R. G. Cai, B. Hu, and H. B. Zhang,
“Scalar graviton in the healthy extension of Hořava-Lifshitz theory", 
Phys. Rev. D 83, 084009 (2011) [arXiv:1008.5048[hep-th]].

EA
T. Jacobson,
 “Extended Hořava gravity and einstein-aether theory",
Phys. Rev. D 81, 101502 (2010) [Phys. Rev. D 82, 129901 (2010)] [arXiv:1001.4823[hep-th]].

 HL_BHTS1
 D. Blas and S. Sibiryakov, 
 “Hořava gravity versus thermodynamics: The black hole case”, 
 Phys. Rev. D 84, 124043 (2011) [arXiv:1110.2195[hep-th]].

 HL_BHTS2
 E. Barausse, T. Jacobson and T. P. Sotiriou, 
 “Black holes in Einstein-aether and Hořava-Lifshitz gravity”, 
 Phys. Rev. D 83, 124043 (2011) [arXiv:1104.2889[gr-qc]].

 HL_BHTS3
 P. Berglund, J. Bhattacharyya and D. Mattingly, 
 “Mechanics of universal horizons", 
 Phys. Rev. D 85, 124019 (2012) [arXiv:1202.4497[hep-th]].

HL_BHTS4
K. Lin, E. Abdalla, R. G. Cai, A. Wang,
 “Universal horizons and black holes in gravitational theories with broken Lorentz symmetry", 
Int. J. Mod. Phys. D 23, 1443004 (2014) [arXiv:1408.5976[gr-qc]].

HL_BHTS5
K. Lin, V. H. Satheeshkumar, A. Wang,
 “Static and rotating universal horizons and black holes in gravitational theories with broken Lorentz invariance", 
Phys. Rev. D 93, 124025 (2016) [arXiv:1603.05708[gr-qc]].

 HL_BHTS6
 Y. Misonoh and K. Maeda, 
 “Black Holes and Thunderbolt Singularities with Lifshitz Scaling Terms", 
 Phys. Rev. D 92, 084049 (2015) [arXiv:1509.01378[gr-qc]].

inflation_Horava
E. G. M. Ferreira and R. Brandenberger, 
“Trans-Planckian problem in the healthy extension of Hořava-Lifshitz gravity", 
Phys. Rev. D 86, 043514 (2012) [arXiv:1204.5239[hep-th]].

SVW_paper
T. P. Sotiriou, M. Visser, S. Weinfurtner, 
“Phenomenologically viable Lorentz-violating quantum gravity", 
Phys.Rev. Lett. 102, 251601 (2009) [arXiv:0904.4464[hep-th]]; 
T. P. Sotiriou, M. Visser, S. Weinfurtner, 
“Quantum gravity without Lorentz invariance", 
JHEP 0910, 033 (2009) [arXiv:0905.2798[hep-th]].

GW
B. P. Abbott et al. [LIGO Scientific Collaboration, Virgo Collaboration], 
“GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral", 
Phys. Rev. Lett. 119, 161101 (2017) [arXiv:1710.05832[gr-qc]];
B. P. Abbott et al. [LIGO Scientific Collaboration, Virgo Collaboration, Fermi Gamma-Ray Burst Monitor, INTEGRAL], 
“Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A", 
Astrophys. J. 848 no.2, L13 (2017) [arXiv:1710.05834 [astro-ph.HE]].

ref_Hill_eq
G. W. Hill, 
“On the part of the motion of the lunar perigee which is a function of the mean motions of the sun and moon", 
Acta Math. 8, 1 (1886).

U1_PHL1
P. Horava, and C. M. Melby-Thompson, 
“General Covariance in Quantum Gravity at a Lifshitz Point", 
Phys.Rev. D 82, 064027 (2010) [arXiv:1007.2410[hep-th]].

U1_PHL2
A. Wang, and Y. Wu, 
“Cosmology in nonrelativistic general covariant theory of gravity", 
Phys.Rev. D 83, 044031 (2011) [arXiv:1009.2089[hep-th]].

U1_NPHL1
T. Zhu, Q. Wu, A. Wang, and F. Shu,
“U(1) symmetry and elimination of spin-0 gravitons in Horava-Lifshitz gravity without the projectability condition", 
Phys. Rev. D 84, 101502 (2011) [arXiv:1108.1237[hep-th]].

U1_NPHL2
T. Zhu, F.-W. Shu, Q. Wu, and A. Wang,
“General covariant Horava-Lifshitz gravity without projectability condition and its applications to cosmology", 
Phys. Rev. D 85, 044053 (2012) [arXiv:1110.5106[hep-th]].


 

