Stable Singularity-free Cosmological Solutions in non-projectable Horava-Lifshitz Gravity

TL;DR

This paper discovers stable, singularity-free cosmological solutions within non-projectable Hořava-Lifshitz gravity, addressing scalar perturbation instabilities by relaxing projectability and analyzing background and perturbation dynamics.

Contribution

It introduces a classification of gradient stability types and highlights the role of higher order spatial curvature terms in stabilizing solutions in HL gravity.

Findings

Stable singularity-free solutions in non-projectable HL gravity.

Classification of gradient stability into five types based on perturbation parameters.

Higher order spatial curvature terms suppress instabilities.

Abstract

We find stable singularity-free cosmological solutions in non-flat Friedmann-Lema\^{i}tre-Robertson-Walker (FLRW) spacetime in the context of Ho\v{r}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 $\mathcal{M}^2$, the effective friction coefficients in perturbation equations $\mathcal{H}$ and these magnitude relations $|\mathcal{M}^2|/\mathcal{H}^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.

Tables

Table 1. Table 1: The classification of stabilities and instabilities which is determined by ℳ 2 superscript ℳ 2 \mathcal{M}^{2} and ℋ ℋ \mathcal{H} .

	Stability types	$\mathscr{H}$	${ℳ}^2$	$\left|{ℳ}^2\right|\stackrel{>}{\sim }{\mathscr{H}}^2$
	Negative-${ℳ}^2$ instability	$-$	$-$	$\times$
		$-$	$-$	$○$
		$+$	$-$	$○$
	$\mathscr{H}$-accelerated instability	$-$	$+$	$\times$
	${ℳ}^2$-dominated stability	$-$	$+$	$○$
	$\mathscr{H}$-suppressed stability	$+$	$-$	$\times$
	Complete stability	$+$	$+$	$○$
		$+$	$+$	$\times$

Table 2. Table 2: The examples for the stable singularity-free solutions, which mean that these solutions are under the ℳ 2 superscript ℳ 2 \mathcal{M}^{2} -dominated stability when ℋ < 0 ℋ 0 \mathcal{H}<0 and under the complete stability when ℋ > 0 ℋ 0 \mathcal{H}>0 . The types of the solutions are introduced in Sec.  III . N/A means that there is no corresponding variables.

		Type	$\Lambda$	$\lambda$	$g_2$	$g_3$	$g_4$	$g_5$	$g_6$	$g_7$	$g_8$	$\varsigma$	$g_r$	$g_s$	$a_{\mathrm{BC}}$	$a_{\min }$	$a_{\max }$	$a_T$
	(i)	${𝒪}^{[0;1]}$	$0$	$1$	$-\frac{1}{18}$	$\frac{1}{3}$	$-\frac{1}{108}$	$0$	$\frac{7}{90}$	$\frac{1}{2}$	$1$	$1$	$1$	$-\frac{1}{15}$	N/A	0.304	0.491	N/A
	(ii)	${ℬ}^{[0;-1]}$	$0$	$6$	$-\frac{1}{18}$	$\frac{1}{8}$	$0$	$\frac{1}{3}$	$-\frac{2}{3}$	$-1$	$3$	$\frac{3}{2}$	$-\frac{1}{4}$	$4$	N/A	N/A	N/A	$1.094$
	(iii)	${ℬ}_{\mathrm{BC}}^{[0;-1]}$	$0$	$\frac{5}{2}$	$-\frac{1}{18}$	$-\frac{1}{3}$	$\frac{1}{20}$	$-\frac{1}{4}$	$\frac{1}{4}$	$0$	$1$	$\frac{1}{2}$	$-3$	$-\frac{3}{5}$	$0.526$	N/A	N/A	$0.851$
	(iv)	${ℬ}^{[1;1]}$	$1$	$1$	$\frac{1}{6}$	$-\frac{1}{6}$	$-\frac{1}{12}$	$\frac{1}{3}$	$-\frac{1}{3}$	$-\frac{1}{2}$	$1$	$\frac{3}{2}$	$2$	$-1$	N/A	N/A	N/A	$1.525$
	(v)	${ℬ}_{\mathrm{BC}}^{[1;1]}$	$2$	$2$	$\frac{1}{9}$	$-\frac{1}{2}$	$-\frac{1}{12}$	$\frac{1}{3}$	$-\frac{1}{6}$	$-\frac{1}{2}$	$1$	$1$	$-1$	$1$	$0.707$	N/A	N/A	$1.272$
	(vi)	${ℬ}_O^{[1;1]}$	$2$	$1$	$\frac{2}{45}$	$0$	$-\frac{1}{108}$	$\frac{2}{75}$	$0$	$0$	$2$	$\frac{1}{2}$	$\frac{4}{5}$	$-\frac{1}{25}$	N/A	$0.256$	$0.511$	$1.083$
	(vii)	${ℬ}^{[1;-1]}$	$1$	$\frac{5}{2}$	$0$	$-\frac{1}{6}$	$0$	$\frac{1}{6}$	$-\frac{1}{3}$	$\frac{1}{2}$	$2$	$1$	$-1$	$2$	N/A	N/A	N/A	$0.928$
	(viii)	${ℬ}_{\mathrm{BC}}^{[1;-1]}$	$2$	$3$	$-\frac{1}{3}$	$\frac{1}{2}$	$-\frac{1}{30}$	$\frac{1}{5}$	$-\frac{1}{3}$	$-1$	$1$	$\frac{3}{2}$	$-3$	$-\frac{2}{5}$	$0.403$	N/A	N/A	$0.745$
	(ix)	${𝒪}^{[-1;1]}$	$-1$	$1$	$-\frac{1}{12}$	$\frac{1}{2}$	$-\frac{1}{108}$	$\frac{1}{40}$	$0$	$1$	$1$	$1$	$\frac{3}{2}$	$-\frac{1}{10}$	N/A	0.282	0.604	N/A
	(x)	${𝒪}^{[-1;-1]}$	$-1$	$3$	$\frac{1}{4}$	$-\frac{1}{2}$	$0$	$-\frac{1}{40}$	$\frac{1}{5}$	$0$	$3$	$1$	$\frac{3}{2}$	$\frac{3}{2}$	N/A	0.730	1.821	N/A
	(xi)	${𝒪}_{\mathrm{BC}}^{[-1;-1]}$	$-\frac{3}{2}$	$\frac{5}{2}$	$\frac{1}{8}$	$-\frac{1}{2}$	$\frac{1}{216}$	$-\frac{1}{60}$	$\frac{1}{150}$	$-1$	$2$	$\frac{1}{2}$	$-\frac{3}{4}$	$-\frac{1}{50}$	0.174	0.507	1.309	N/A

Table 3. Table 3: The maximum ratios h ( G ) / h ( G ) ​ ini subscript ℎ 𝐺 subscript ℎ 𝐺 ini h_{(G)}/h_{(G)\mathrm{ini}} and h ( Q ) / h ( Q ) ​ ini subscript ℎ 𝑄 subscript ℎ 𝑄 ini h_{(Q)}/h_{(Q)\mathrm{ini}} throughout the evolutions for each perturbation modes based on TABLE  2 (ii).

	${\nu }^2$	$\max \left(h_{\left(G\right)}/h_{\left(G\right)\mathrm{ini}}\right)$	$\max \left(h_{\left(Q\right)}/h_{\left(Q\right)\mathrm{ini}}\right)$
	$2$	49.425	32.907
	$5$	33.761	20.835
	$10$	24.876	13.945
	$50$	11.612	5.933
	$100$	8.263	4.132
	$300$	4.795	2.308

Table 4. Table 4: An example of the coupling constants for the oscillating solution with resonance.

	type	$\Lambda$	$\lambda$	$g_2$	$g_3$	$g_4$	$g_5$	$g_6$	$g_7$	$g_8$	$\varsigma$	$g_r$	$g_s$	$a_{\min }$	$a_{\max }$
	${𝒪}^{[-1;-1]}$	$-1$	$3$	$\frac{1}{4}$	$-\frac{1}{2}$	$0$	$\frac{1}{20}$	$-\frac{1}{15}$	$-\frac{1}{5}$	$\frac{1}{5}$	$\frac{3}{2}$	$\frac{3}{2}$	$1$	0.629	1.832

Equations

$$t\to b^{-z}t\mathrm{and}{x}^i\to b^{-1}{x}^i,$$

$$N(t,x^i)\to N(t),$$

$$t\to t+f(t),x^i\to x^i+{\zeta }^i(t,x^j),$$

$${\mathcal{K}}_{ij}:=\frac{1}{2N}\left({\dot{g}}_{ij}-{\nabla }_i{N}_j-{\nabla }_i{N}_j\right),$$

$${\Phi }_i:=\frac{{\nabla }_{i}N}{N},$$

$$S=\frac{m_{\mathrm{LV}}^2}{2}\int dt{d}^{3}x\left({\mathcal{L}}_{\mathrm{K}}+{\mathcal{L}}_{\mathrm{P}}+{\mathcal{L}}_{\mathrm{NP}}\right),$$

$${\mathcal{L}}_{\mathrm{K}}$$

$${\mathcal{L}}_{\mathrm{P}}$$

$${\mathcal{L}}_{\mathrm{NP}}$$

$${\mathcal{V}}_{z=1}$$

$${\mathcal{V}}_{z=2}$$

$${\mathcal{V}}_{z=3}$$

$$\mathcal{V}[{\Phi }_i]:=\varsigma {\Phi }^2.$$

$$d{s}^2$$

$$f(\chi ):=\{\begin{array}{ll}\chi & \mathrm{for}K=0\\ \sin \chi & \mathrm{for}K=1\\ \sinh \chi & \mathrm{for}K=-1\end{array},$$

$$2\dot{H}+3H^2=\frac{2}{3\lambda -1}\left(\Lambda -\frac{K}{a^2}-\frac{g_r{K}^2}{3a^4}-\frac{g_s{K}^3}{a^6}\right),$$

$$H^2=\frac{2}{3(3\lambda -1)}\left[\Lambda -3\frac{K}{a^2}+\frac{g_r{K}^2}{a^4}+\frac{g_s{K}^3}{a^6}\right],$$

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

$$g_r$$

$$g_s$$

$$H^2=\frac{2}{3(3\lambda -1)}\left[\Lambda -3\frac{K}{a^2}+\frac{g_r{K}^2}{a^4}+\frac{g_s{K}^3}{a^6}+\frac{\mathcal{C}}{a^3}\right],$$

$$\frac{1}{2}{\dot{a}}^2+\mathcal{U}(a)=0,$$

$$\mathcal{U}(a)=\frac{1}{3\lambda -1}\left[K-\frac{\Lambda }{3}{a}^2-\frac{g_r{K}^2}{3a^2}-\frac{g_s{K}^3}{3a^4}\right].$$

$$a_{\pm }^{[K]}=\sqrt{\frac{K}{6}\left(g_r\pm \sqrt{g_r^2+12g_s}\right)}.$$

$$g_s<0,g_r>0\mathrm{and}{g}_r^2+12g_s>0$$

$$\mathrm{with}{a}_{\max }=a_{+}^{[1]}\mathrm{and}{a}_{\min }=a_{-}^{[1]}.$$

$$(\mathrm{i})$$

$$(\mathrm{ii})$$

$$g_s<0,g_r<0\mathrm{and}{g}_r^2+12g_s>0$$

$$\mathrm{with}{a}_T=a_{-}^{[-1]}\mathrm{and}{a}_{\mathrm{BC}}=a_{+}^{[-1]},$$

$$\tilde{\mathcal{U}}(\tilde{a})=\frac{1}{3\lambda -1}\left[K-\epsilon {\tilde{a}}^2-\frac{{\tilde{g}}_r{K}^2}{3{\tilde{a}}^2}-\frac{{\tilde{g}}_s{K}^3}{3{\tilde{a}}^4}\right],$$

$${\tilde{a}}_I^{[\epsilon ;K]}$$

$${\tilde{\xi }}_I^{[\epsilon ;K]}$$

$${\tilde{g}}_{s(\pm )}^{[\epsilon ;K]}$$

$${\tilde{a}}_{\pm }^{[\epsilon ;K]}:=\sqrt{\frac{K}{2\epsilon }\left[1\pm \mathrm{sgn}(K)\sqrt{1-\frac{4\epsilon {\tilde{g}}_r}{3}}\right]}.$$

$$(\mathrm{i})$$

Full text

††preprint: WU-AP/1809/18

I 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 Penrosesingularity_theorem , 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 remediedinflation_not_complete .

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 braneworldbraneworld_bounce , string gasstring_gas , loop quantum gravityLQG_bounce1 ; LQG_bounce2 ; LQG_bounce3 , Horndeski theory (or generalized Galileon)genesis1 ; genesis2 ; genesis3 ; genesis4 ; Galileon_bounce ; Galileon_bounce2 ; Galileon_bounce3 , and the non-local gravitational theorynon-local_bounce . 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 unstableno-go_Horndeski ; no-go_multi ; no-go_non-flat ; no-go_papers . 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 orderHorndeski_theory , one may consider it is difficult to find stable cosmological solutions without a singularity. However, certain loopholes of the no-go theorem are knownbounce_EFT ; no-go_loopholes1 ; no-go_loopholes2 ; no-go_loopholes3 . 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 termsbounce_EFT .

Hořava-Lifshitz (HL) theory is known as a gravitational theory without Lorentz symmetry, which is a candidate for quantum gravityHorava (the recent progresses of the HL theory are reviewed in HL-review2 ). 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,

[TABLE]

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 termsHorava ; ref_renormalize_HL ; power-counting_renormalizability .

In the context of the HL theory, a number of attractive applications for cosmology have been conducted (see HL-review1 for a review): for example, primordial perturbationsCP1 ; CP2 ; CP3 ; CP4 ; CP5 , gravitational wavesHLC1 ; GW-HL and other cosmological aspectsHLC2 ; DM_integrate ; HLC3 ; dS_projectable . It is remarkable that the bouncing and oscillating solutions are discovered as singularity-free cosmological models in HL theorybounce_HL1 ; bounce_HL2 ; bounce_HL3 ; bounce_HL4 ; previous_oscillating_universe ; Bianchi_IX . 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 paperprevious_stability_bounce , we investigated the stability of the singularity-free cosmological solutions based on the projectable HL theory. Namely, we impose following the projectability condition:

[TABLE]

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 irrelevantproblem_projectable1 ; problem_projectable2 ; problem_projectable3 ; problem_projectable4 . 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 theoryhealthy_extension . It is discovered that the gradient instability in flat FLRW spacetime can be avoidedNP_stability1 ; NP_stability2 . 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 featureHL-review2 ; EA ; HL_BHTS1 ; HL_BHTS2 ; HL_BHTS3 ; HL_BHTS4 ; HL_BHTS5 ; HL_BHTS6 ; inflation_Horava , 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. II, 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. III. In Sec. IV, 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. V, 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. VI is devoted to summary and discussions.

II 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:

[TABLE]

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 (3), the terms in action must be composed of the following quantities: The extrinsic curvature which is defined by

[TABLE]

where, the dot symbol represents the derivative with respect to time coordinate and $\nabla_{i}$ denotes the spatial covariant derivative. The three dimensional Ricci curvature $\mathcal{R}_{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:

[TABLE]

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., $\nabla_{i}\mathcal{K}_{jk},\nabla_{i}\mathcal{R}_{jk},\nabla_{i}\Phi_{j},\nabla_{i}\nabla_{j}\mathcal{K}_{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:

[TABLE]

with

[TABLE]

where, $m_{\mathrm{LV}}$ is a Lorentz violating scale which is expected to be around the Planck scale, $\lambda$ and $g_{i}$ ($i=1$-$8$) are dimension-less coupling constants and $\Lambda$ is a cosmological constant. $\mathcal{L}_{\mathrm{NP}}$ is constructed by the terms including $\Phi_{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 actionSVW_paper . 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 $\mathcal{L}_{\mathrm{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 $\Phi^{2}:=\Phi_{i}\Phi^{i}$. This diagnosis seems not to be inconsequent. According to the papers healthy_extension , the stabilities of Minkowski and flat FLRW spacetime are remedied by the non-projectable HL theory with $\Phi^{2}$ term. Therefore, we do minimally extend the theory by adding only $\Phi^{2}$ term:

[TABLE]

with a dimensionless coupling constant $\varsigma$.

In the rest of this paper, the unit $m_{\mathrm{LV}}=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 unityGW ; GW-HL .

III 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.

III.1 Basic equations

To unify the description, we consider the FLRW spacetime in a spatial coordinate system $x^{i}=(\chi,\theta,\phi)$:

[TABLE]

with

[TABLE]

where, $a$ is the scale factor which depends only on the cosmic time, and $K$ is related to the spatial Ricci curvature as $\mathcal{R}=6K/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 $\chi$ is defined by $0\leq\chi<\infty$ for $K=0,-1$ and $0\leq\chi\leq\pi$ for $K=1$. Furthermore, angular coordinates $\theta$ and $\phi$ take $0\leq\theta\leq\pi$ and $0\leq\phi\leq 2\pi$, 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$:

[TABLE]

where, $H:=\dot{a}/a$ is the Hubble parameter, $g_{r}$ and $g_{s}$ are the linear combinations of the coupling constants111 It should be noted that we adopt the different definitions of $g_{r}$ and $g_{s}$ from our previous paperprevious_stability_bounce to simplify the perturbed action. The previous definitions are

$\displaystyle g_{r}:=6K^{2}(3g_{2}+g_{3})\,,~{}g_{s}:=12K^{3}(9g_{4}+3g_{5}+g_{6})\,,$

thus, the sign of $g_{s}$ is flipped in open FLRW spacetime. :

[TABLE]

We have already fixed the gauge as $N=1$ and $N_{i}=0$. In what follows, we assume $\lambda>1/3$.

As we know, the equation of the scale factor (14) can be derived by taking derivative of (15) with respect to the cosmic time. It means the independent equation is only (15). 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:

[TABLE]

with dust-like term with an integration constant $\mathcal{C}$DM_integrate . Thus, if we consider $\mathcal{C}=0$ case in the projectable HL theory, the same background FLRW solutions are realizedprevious_oscillating_universe ; previous_stability_bounce .

III.2 Singularity-free background solutions

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

[TABLE]

with

[TABLE]

Since the possible ranges for the scale factor are limited within the region in which $\mathcal{U}\leq 0$, the background evolution is completely determined by the coupling constants in $\mathcal{L}_{\mathrm{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_{s}K^{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.previous_oscillating_universe . In our analysis, we focus on the following two kinds of singularity-free solutions. One is a bouncing solution denoted by $\mathcal{B}^{[\mathrm{sgn}(\Lambda);K]}$, where the function $\mathrm{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 $\mathcal{O}^{[\mathrm{sgn}(\Lambda);K]}$. A universe shows periodic oscillating behavior, in other words, bounces at $a=a_{\mathrm{min}}$ and recollapses222 In this paper, a word of recollapse means that the expanding universe shifts to contracting phase. Thus, a recollapcing point satisfies the conditions $\mathcal{U}=0$ and $\mathcal{U}^{\prime}>0$. at $a=a_{\mathrm{max}}$. Therefore the scale factor is limited within $0<a_{\mathrm{min}}\leq a\leq a_{\mathrm{max}}<\infty$. 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.

III.2.1 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 $\mathcal{U}$ determines the possible ranges for the scale factor, it is convenient to derive the roots of the equation $\mathcal{U}=0$:

[TABLE]

The points $a=a^{[K]}_{\pm}$ correspond to the bouncing or recollapsing points of the universe. Of course, the corresponding $a^{[K]}_{\pm}$ must be real and positive to find such points. By examining the forms of the potential $\mathcal{U}$, we find the following three types of singularity-free solutions.

1. (a)

$\mathcal{O}^{[0;1]}$— A universe shows oscillating behavior in closed FLRW spacetime ($K=1$), which we call $\mathcal{O}^{[0;1]}$. This type of the solutions can be found if the following conditions are satisfied.

[TABLE]

In open FLRW spacetime, this kind of solutions is never found. 2. (b)

$\mathcal{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:

[TABLE]

The solutions satisfying the conditions (i) and (ii) are referred to as $\mathcal{B}^{[0;-1]}$(i) and $\mathcal{B}^{[0;-1]}$(ii), respectively. In closed FLRW spacetime, we cannot find this type of the solutions. 3. (c)

$\mathcal{B}^{[0;-1]}_{\mathrm{BC}}$— A universe shows bouncing behavior if the initial value of scale factor $a_{\mathrm{ini}}$ is larger than $a_{T}$, whereas, falls down into the singularity if $a_{\mathrm{ini}}\leq a_{\mathrm{BC}}$. Since this type of the solutions possibly induces a Big Crunch (BC), we add a subscript BC. Note that the domain $a_{\mathrm{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$), $\mathcal{B}^{[0;-1]}_{\mathrm{BC}}$ can be found:

[TABLE]

where, $a_{\mathrm{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 $\mathcal{O}^{[0,1]}$, $\mathcal{B}^{[0,-1]}$ and $\mathcal{B}^{[0,-1]}_{\mathrm{BC}}$ in FIG. 1 and the distribution of these types of solutions on $(g_{r},g_{s})$ plane in FIG. 2.

As we will discuss later, we can construct the bouncing solutions $\mathcal{B}^{[0;-1]}$ and $\mathcal{B}^{[0;-1]}_{\mathrm{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.

III.2.2 With positive cosmological constant

Although the properties of singularity-free solutions for $\Lambda>0$ have already given in the papers previous_oscillating_universe ; previous_stability_bounce , 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; $\tilde{a}:=a/\ell$, $\tilde{g}_{r}:=g_{r}/\ell^{2}$ and $\tilde{g}_{s}:=g_{s}/\ell^{4}$ with $\ell:=\sqrt{3/|\Lambda|}$. Then, the potential $\mathcal{U}$ is rewritten as

[TABLE]

where $\varepsilon:=\pm 1$ expresses the sign of the cosmological constant. The three roots of the equation $\tilde{\mathcal{U}}(\tilde{a})=0$ are given by

[TABLE]

with

[TABLE]

where, $\mathrm{pv}$ means a principal value of cubic root and $I=1,2,3$. If $\tilde{a}^{[\varepsilon;K]}_{I}$ takes real and positive value, $\tilde{\mathcal{U}}(\tilde{a})=0$ possesses a corresponding root in the possible ranges for the scale factor. We further derive the roots of $\tilde{\mathcal{U}}(\tilde{a})=0$ with $\tilde{g}_{s}=0$:

[TABLE]

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

1. (a)

$\mathcal{B}^{[1;K]}$— A bouncing solution whose bouncing radius is given as $\tilde{a}_{T}$. We refer to this type of the solutions as $\mathcal{B}^{[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 $\mathcal{B}^{[1;1]}$ as the following three cases:

[TABLE]

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

[TABLE]

$\mathcal{B}^{[1;-1]}$(i) and $\mathcal{B}^{[1;-1]}$(ii) correspond to the conditions (i) and (ii), respectively. 2. (b)

$\mathcal{B}^{[1;K]}_{\mathrm{BC}}$— A universe with $\tilde{a}_{\mathrm{ini}}\leq\tilde{a}_{\mathrm{BC}}$ evolves into a Big Crunch, while it shows bouncing behavior if $\tilde{a}_{\mathrm{ini}}\geq\tilde{a}_{T}>\tilde{a}_{\mathrm{BC}}$. We classify this type of the solutions as $\mathcal{B}^{[1;K]}_{\mathrm{BC}}$. Unlike $\mathcal{B}^{[0;K]}_{\mathrm{BC}}$, the asymptotic behavior of spacetime after bounce is de Sitter spacetime. For closed FLRW spacetime ($K=1$), we find the following two cases:

[TABLE]

The solutions satisfying the conditions (i) and (ii) are referred to as $\mathcal{B}^{[1;1]}_{\mathrm{BC}}$(i) and $\mathcal{B}^{[1;1]}_{\mathrm{BC}}$(ii), respectively. For open FLRW spacetime ($K=-1$), the solutions $\mathcal{B}^{[1;-1]}_{\mathrm{BC}}$ can be seen when the following conditions are satisfied:

[TABLE] 3. (c)

$\mathcal{B}^{[1;1]}_{O}$— A universe in closed FLRW spacetime shows oscillating behavior if the initial scale factor is in $\tilde{a}_{\mathrm{min}}\leq\tilde{a}_{\mathrm{ini}}\leq\tilde{a}_{\mathrm{max}}$, or bounces if $\tilde{a}_{\mathrm{ini}}\geq\tilde{a}_{T}>\tilde{a}_{\mathrm{max}}$. We classify this type of solutions as $\mathcal{B}^{[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.

[TABLE]

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 $(\tilde{g}_{r},\tilde{g}_{s})$ plane are shown in FIG. 3 and FIG. 4, respectively.

III.2.3 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 $(\varepsilon=-1)$, which we cannot construct the bouncing solutions. We find the following two kinds of singularity-free oscillating solutions:

1. (a)

$\mathcal{O}^{[-1;K]}$— A universe which shows periodic oscillation whose oscillating amplitude is given by $\tilde{a}_{\mathrm{min}}\leq\tilde{a}\leq\tilde{a}_{\mathrm{max}}$. We classify this type of the solutions as $\mathcal{O}^{[-1;K]}$. If the following conditions are satisfied in closed FLRW spacetime ($K=1$), $\mathcal{O}^{[-1;1]}$ is realized:

[TABLE]

In open FLRW spacetime ($K=-1$), we find the two cases of $\mathcal{O}^{[-1;-1]}$ are obtained:

[TABLE]

The solutions satisfying the conditions (i) and (ii) are referred to as $\mathcal{O}^{[-1;-1]}$(i) and $\mathcal{O}^{[-1;-1]}$(ii), respectively. 2. (b)

$\mathcal{O}_{\mathrm{BC}}^{[-1;-1]}$— A universe in open FLRW spacetime shows oscillating behavior if the initial radius of the universe is in $\tilde{a}_{\mathrm{min}}\leq\tilde{a}_{\mathrm{ini}}\leq\tilde{a}_{\mathrm{max}}$, or falls into the singularity if $\tilde{a}_{\mathrm{ini}}<\tilde{a}_{\mathrm{BC}}<\tilde{a}_{\mathrm{min}}$. We refer to this type of solutions as $\mathcal{O}_{\mathrm{BC}}^{[-1;-1]}$. For open FLRW spacetime ($K=-1$), we find the solutions under the following conditions:

[TABLE]

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 $(\tilde{g}_{r},\tilde{g}_{s})$ plane are shown in FIG. 5 and FIG. 6, respectively.

IV 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

[TABLE]

with

[TABLE]

where, $\bar{N}$, $\bar{N}_{i}$ and $\bar{g}_{ij}$ denote the background lapse function, shift vector and three-dimensional induced metric, respectively. Furthermore, $\alpha$, $\beta_{i}$ and $h_{ij}$ mean first order perturbation of the above ADM variables. Note that indices of these perturbed variables are raised by $\bar{g}^{ij}$ as $h^{i}_{j}:=\bar{g}^{ia}h_{aj}$, $h^{ij}:=\bar{g}^{ia}\bar{g}^{jb}h_{ab}$, $h:=h^{i}_{i}=\bar{g}^{ia}h_{ia}$ and $\beta^{i}:=\bar{g}^{ia}\beta_{a}$.

IV.1 Harmonic expansion

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

[TABLE]

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. previous_stability_bounce .

The scalar perturbations of the ADM variables are decomposed into

[TABLE]

the vector perturbations also can be expanded by

[TABLE]

and the tensor part is given by

[TABLE]

where, the degrees $l,m\in\mathbb{Z}$ are constrained by $0\leq l\leq n-1$ and $0\leq|m|\leq l$. $n\geq 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 $\alpha$ is recovered. Thus, we expand $\alpha^{(\mathrm{scalar})}$ by the harmonic functions.

IV.2 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 (3), the infinitesimal transformations of the perturbed ADM quantities are given by

[TABLE]

We stress that $f$ does not depend on space, thus, only $\zeta^{i}$ can be expanded by the harmonics:

[TABLE]

Then, the infinitesimal gauge transformations of the harmonic expansion (32)-(IV.1) are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, $\nu^{2}$ is a eigenvalue of the harmonics which is defined by the following regions $(n\geq 1)$:

[TABLE]

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

[TABLE]

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

IV.3 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.

IV.3.1 Tensor perturbations

Since the additional perturbation terms coming from $\mathcal{L}_{\mathrm{NP}}$ are only scalar modes, the tensor part of the quadratic action is identical to that of the projectable case:

[TABLE]

where, we introduced $\mathcal{F}_{(G)}$ and $\mathcal{G}_{(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

[TABLE]

and we define $g_{56}:=g_{5}+g_{6}$. The tensor quadratic action is defined for the case with $l\geq 2$, because the tensor harmonics $G^{(n;lm)}_{(o)ij}$ and $G^{(n;lm)}_{(e)ij}$ are vanished when $l<2$.

IV.3.2 Scalar perturbations

The scalar perturbations are drastically changed since the lapse perturbation $\alpha_{(Q)}$ cannot be eliminated by the gauge condition. Furthermore, $\mathcal{L}_{\mathrm{NP}}$ adds a new degree of freedom to scalar perturbations. The quadratic action of the scalar perturbations is given by

[TABLE]

We eliminated $\dot{H}$ by applying the scale factor equation (14). Since $\alpha_{(Q)}$ and $\beta_{(Q)}$ are not dynamical, we can eliminate both of them with the following constraint equations:

[TABLE]

Substituting, (57) and (58) into (56), we achieve the quadratic action of scalar perturbations:

[TABLE]

where,

[TABLE]

In closed FLRW spacetime, we focus only on the case with $n\geq 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 ($\nu^{2}-3K$) as a positive value. It should be noted that, even if we take the limit $\varsigma\to 0$, we cannot replicate the quadratic scalar action in the projectable case (see (3.26) and (3.27) in ref.previous_stability_bounce ). In fact, a difference is caused by the gauge structure. As we have seen, $\alpha_{(Q)}$ is eliminated by applying constraint equation (57). Due to the recovering the local lapse function, it cannot be removed by global infinitesimal transformation $f$.

V 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 $\mathcal{F}$ and $\mathcal{G}$ in the quadratic action provides guideposts (Since the discussion holds for both scalar and tensor modes, the subscripts are abbreviated.). The sign of $\mathcal{F}$ is relevant to the ghost instability. A perturbation mode with negative $\mathcal{F}$ losses the lowest energy state which leads the fatal collapse of the perturbative approach. Thus, we consider the condition $\mathcal{F}>0$ as the absolute requirement.

On the other hand, the case with $\mathcal{G}<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:

[TABLE]

where,

[TABLE]

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

[TABLE]

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

Thus, we pay deep attention to the values of $\mathcal{M}^{2}$ and $\mathcal{H}$ to discuss the stability of the perturbations. In this paper, we classify the stabilities and instabilities into the following five types.

1. (i)

Negative-$\mathcal{M}^{2}$ instability: As we mentioned, the negative $\mathcal{M}^{2}$ causes undesirable exponential growth unless it is not suppressed by the $\mathcal{H}$-friction. Therefore, the unstable behavior due to the negative $\mathcal{M}^{2}$ is observed in the case with (1) $\mathcal{M}^{2}<0$ and $\mathcal{H}<0$, (2) $\mathcal{M}^{2}<0$ and |\mathcal{M}^{2}|\,\mbox{ \raisebox{-4.30554pt}{\stackrel{{\scriptstyle\textstyle>}}{{\textstyle\sim}}}}\,\mathcal{H}^{2}>0 which means the effect of the third term in LHS of (62) 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-$\mathcal{M}^{2}$ instability. 2. (ii)

$\mathcal{H}$-accelerated instability: Even if the squared effective mass is positive, it is possible to enhance the oscillating amplitude of perturbations. When $\mathcal{H}<0$, the second term in LHS of (62) reinforces the amplitude rather than friction, thus the perturbation experiences a $\mathcal{H}$-acceleration. Clearly, the accelerating effect is manifested if the positive squared effective mass is relatively smaller than the effect of $\mathcal{H}$ term. Thus, we regard the solution is under a $\mathcal{H}$-accelerated instability if at least one of the perturbation modes satisfies 0<\mathcal{M}^{2}\,\mbox{ \raisebox{-4.30554pt}{\stackrel{{\scriptstyle\textstyle<}}{{\textstyle\sim}}}}\,\mathcal{H}^{2} and $\mathcal{H}<0$. 3. (iii)

$\mathcal{M}^{2}$-dominated stability: On the other hand, the negative $\mathcal{H}$ is not problematic if it is sufficiently suppressed by the heavy positive effective squared mass. Thus, we call the solution is under a $\mathcal{M}^{2}$-dominated stability if every perturbation mode satisfy \mathcal{M}^{2}\,\mbox{ \raisebox{-4.30554pt}{\stackrel{{\scriptstyle\textstyle>}}{{\textstyle\sim}}}}\,\mathcal{H}^{2} and $\mathcal{H}<0$. 4. (iv)

$\mathcal{H}$-suppressed stability: If $\mathcal{H}>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 $\mathcal{M}^{2}<0$. That is, the effective squared mass is negative, however, the unstable behavior is suppressed by the positive $\mathcal{H}$. We call the solution is under a $\mathcal{H}$-suppressed stability if every perturbation mode satisfy |\mathcal{M}^{2}|\,\mbox{ \raisebox{-4.30554pt}{\stackrel{{\scriptstyle\textstyle<}}{{\textstyle\sim}}}}\,\mathcal{H}^{2} and $\mathcal{M}^{2}<0$. 5. (v)

Complete stability: Clearly, there is no problematic perturbation dynamics if every perturbation mode satisfy $\mathcal{M}^{2}>0$ and $\mathcal{H}>0$. We call this situation as a complete stability.

In TABLE 1, 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 $\mathcal{G}$ and $\mathcal{F}$ for general cases. Therefore, our main procedure is to numerically trace $\mathcal{G}$ and $\mathcal{F}$ 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.

V.1 Late-time universe after bounce

When the higher order curvature terms are neglected, the forms of $\mathcal{G}$ and $\mathcal{F}$ are quite simplified. The tensor perturbations are approximated as

[TABLE]

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 $\mathcal{G}$ and $\mathcal{F}$ in the late-time universe after bounce.

V.1.1 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$, $\Lambda=0$ and the scale factor is large, then, $\mathcal{F}_{(Q)}$ and $\mathcal{G}_{(Q)}$ are given byhealthy_extension

[TABLE]

Since $\mathcal{G}_{(Q)}$ cannot be positive if $\varsigma<0$, we focus only on the case with $\varsigma>0$. To preserve $\mathcal{F}_{(Q)}>0$ and $\mathcal{G}_{(Q)}>0$, we must impose

[TABLE]

Then, the solutions are under the complete stability.

V.1.2 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

[TABLE]

Then, $\mathcal{F}_{(Q)}$ and $\mathcal{G}_{(Q)}$ are approximated by

[TABLE]

Then, the conditions for the positive $\mathcal{F}_{(Q)}$ and $\mathcal{G}_{(Q)}$ are given by

[TABLE]

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

[TABLE]

and those of the effective squared masses are reduced into

[TABLE]

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

[TABLE]

with integration constants $C_{1}$ and $C_{2}$. One can see that the perturbation amplitude approaches to zero.

V.1.3 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

[TABLE]

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, $\mathcal{F}_{(Q)}$ and $\mathcal{G}_{(Q)}$ are approximated by333 It is natural to consider that there exists an ultraviolet cut-off momentum $p_{\mathrm{cut}}$ at which a non-perturbative quantum effect of gravity cannot be ignored. More precisely, let $\nu_{\mathrm{cut}}$ be an ultraviolet cut-off mode which is related with the cut-off momentum as $\nu^{2}_{\mathrm{cut}}/a_{T}^{2}\approx p^{2}_{\mathrm{cut}}$ at the early stage of the universe. Then, we find some finite time so that $\nu^{2}/a^{2}<H^{2}$ for any $\nu\leq\nu_{\mathrm{cut}}.$ Therefore the all perturbation modes are rapidly suppressed by the strong $\mathcal{H}$-friction after that time.

[TABLE]

The positivity of $\mathcal{F}_{(Q)}$ is satisfied if $\varsigma>0$ which is consistent with that of in the asymptotic flat spacetime. We show the explicit conditions for $\mathcal{G}_{(Q)}\geq 0$ in closed spacetime:

[TABLE]

and in open spacetime:

[TABLE]

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 $\mathcal{M}_{(Q)}^{2}$ approach asymptotically to zero because it is proportional to $a^{-2}$. Therefore the case with $\mathcal{G}_{(Q)}<0$ seems not to be quite problematic because the solution is under the $\mathcal{H}$-suppressed stability. To confirm it, we examine the asymptotic dynamics of the perturbations. The effective friction coefficients asymptotically behave as

[TABLE]

and those of the effective squared masses are

[TABLE]

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

[TABLE]

and those of the scalar perturbations are given by

[TABLE]

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)}\sim C_{4}$ and $h_{(Q)}\sim 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 $\mathcal{M}^{2}<0$, the dynamics of the tensor perturbations are approximated as

[TABLE]

and those of the scalar perturbations are given by

[TABLE]

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)}\sim C_{8}$ and $h_{(Q)}\sim C_{9}$, which means the $\mathcal{H}$-suppressed stability is truly stabilized state.

V.2 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.

V.2.1 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 $\mathcal{G}_{(G)}$. As we mentioned, the forms of $\mathcal{F}_{(G)}$ and $\mathcal{G}_{(G)}$ in tensor perturbation are identical to those of the projectable case. Thus, we summarize the positivity conditions of $\mathcal{G}_{(G)}$, which we have argued in previous paperprevious_stability_bounce .

To investigate the stability of the tensor perturbations, we focus on the case with large $\nu^{2}$. Then, $\mathcal{G}_{(G)}$ is approximated as

[TABLE]

One can see that the ultraviolet stability imposes $g_{8}\geq 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, $\mathcal{G}_{(G)}$ must be positive for any $a>0$ and $\nu^{2}\geq 8$:

[TABLE]

Similarly, in open spacetime ($K=-1$):

[TABLE]

which are equivalent to $\mathcal{G}_{(G)}\geq 0$ for any $a>0$ and $\nu^{2}\geq 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 previous_stability_bounce ).

V.2.2 Scalar perturbations

The ghost-free condition for scalar perturbations can be analytically discussed. Let the coupling constants be $\lambda>1/3$ and $\varsigma>0$, and adopting $\nu^{2}-3K>0$, the condition $\mathcal{F}_{(Q)}>0$ is equivalent to

[TABLE]

The most stringent conditions are imposed at the bouncing and the recollapse points, i.e., $H=0$:

[TABLE]

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 $\mathcal{G}_{(Q)}$ cannot be expressed without any approximation due to the extremely complicated form. Therefore, we restrict our analysis to the case with large $\nu^{2}$ modes at this stage. Then, $\mathcal{G}_{(Q)}$ is given by

[TABLE]

Since $\varsigma$ should be positive, we assume

[TABLE]

Then, $\mathcal{G}_{(Q)}\geq 0$ for large $\nu$ is guaranteed.

V.2.3 Numerical analysis

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

1. (1)

The cases with $K=1$ and $\Lambda\leq 0$: Since there exists a finite maximum value of scale factor $a_{\mathrm{max}}$, i.e., the solutions $\mathcal{O}^{[0;1]}$ and $\mathcal{O}^{[-1;1]}$, we assume the stability conditions around the bouncing and the recollapsing points and in the large $\nu^{2}$ region, thus,

[TABLE] 2. (2)

The cases with $K=1$ and $\Lambda>0$: In this cases, we observe $\mathcal{B}^{[1;1]}$, $\mathcal{B}^{[1;1]}_{\mathrm{BC}}$ and $\mathcal{B}^{[1;1]}_{\mathrm{O}}$ as the singularity-free solutions. To guarantee the stabilities around bouncing point (and the recollapsing point if $\mathcal{B}^{[1;1]}_{\mathrm{O}}$), in asymptotic de Sitter spacetime and in the large $\nu^{2}$ region, we assume

[TABLE]

If we focus only on the oscillating phase which is observed in the solution $\mathcal{B}^{[1;1]}_{\mathrm{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 $\Lambda\leq 0$, that is, (96). 3. (3)

The cases with $K=-1$ and $\Lambda=0$: We find $\mathcal{B}^{[0,-1]}$ and $\mathcal{B}^{[0,-1]}_{\mathrm{BC}}$ as the singularity-free solutions whose asymptotic behaviors are the Milne spacetime. Thus, we impose

[TABLE] 4. (4)

The cases with $K=-1$ and $\Lambda>0$: In this case, the solutions $\mathcal{B}^{[1,-1]}$ and $\mathcal{B}^{[1,-1]}_{\mathrm{BC}}$ are realized. Since both of the solutions approaches the de Sitter spacetime after the bounce, we impose the conditions as follows:

[TABLE] 5. (5)

The cases with $K=-1$ and $\Lambda<0$: The possible singularity-free solutions are $\mathcal{O}^{[-1;-1]}$ and $\mathcal{O}^{[-1;-1]}_{\mathrm{BC}}$. Therefore we impose the stability conditions around the bouncing and recollapsing points and in the large $\nu^{2}$ region:

[TABLE]

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 $\mathcal{M}^{2}$-dominated stability when $\mathcal{H}<0$ and under the complete stability when $\mathcal{H}>0$. The concrete examples are shown in TABLE 2.

Bouncing solutions with asymptotic Milne spacetime

The solutions $\mathcal{B}^{[0;-1]}$ and $\mathcal{B}_{\mathrm{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}\sim a^{-1}$, the effect of the squared effective masses do not suppressed by the Hubble friction, that is, the values of $\mathcal{M}^{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 $\mathcal{M}^{2}_{(Q)}$ to be positive during the whole evolution. In FIG. 7, a typical example of such a solution is shown.

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

One may wonder whether the squared effective masses truly suppress the effects of $\mathcal{H}$-acceleration under the $\mathcal{M}^{2}$-dominated stability, because we do not mention any explicit criterions of the ratio $\mathcal{M}^{2}/\mathcal{H}^{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 $\mathcal{H}$-acceleration.

To clarify the effects, we numerically solve the perturbation equations of motion (62). The initial conditions are given at $a_{\mathrm{ini}}=100a_{T}$ so that $\dot{h}_{(G)}=0=\dot{h}_{(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)\mathrm{ini}}$ and $h_{(Q)\mathrm{ini}}$. Therefore, we trace the ratios to each initial values, i.e., $h_{(G)}/h_{(G)\mathrm{ini}}$ and $h_{(Q)}/h_{(Q)\mathrm{ini}}$. In FIG 8, the time evolutions are shown.

This result is based on the solution $\mathcal{B}^{[0;-1]}$ whose coupling constants are given in TABLE 2 (ii), thus, also corresponds with FIG 7. From the figure, we find that the perturbation amplitudes are enhanced around the bouncing point even if the solution is under the $\mathcal{M}^{2}$-dominated stability. Thus, we should clarify the growth rates of the perturbations. In TABLE 3, we show the detail data of the maximum ratios $h_{(G)}/h_{(G)\mathrm{ini}}$ and $h_{(Q)}/h_{(Q)\mathrm{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 $\nu^{2}$ is lower. Since the large perturbation mode $\nu^{2}$ corresponds to the heavy effective squared mass $\mathcal{M}^{2}$, the result seems appropriate.

Bouncing solutions with asymptotic de Sitter spacetime

The solutions $\mathcal{B}^{[1;1]}$, $\mathcal{B}^{[1;1]}_{\mathrm{BC}}$, $\mathcal{B}^{[1;1]}_{\mathrm{O}}$ with $a\geq a_{T}$, $\mathcal{B}^{[1;-1]}$ and $\mathcal{B}^{[1;-1]}_{\mathrm{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 $\mathcal{H}$-suppressed stability due to the negative $\mathcal{M}_{(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 $\mathcal{M}^{2}/\mathcal{H}^{2}$ in FIG. 9.

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 $\nu^{2}$. Therefore we conclude this solution is under complete stability for any possible $\nu^{2}$.

The magnitude relationships between $\mathcal{M}^{2}$ and $\mathcal{H}^{2}$ are also analyzed. As we discussed, both $\mathcal{M}^{2}_{(G)}$ and $\mathcal{M}^{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 $\mathcal{M}^{2}$ is equal to $\mathcal{H}^{2}$. We define such a value of scale factor as a critical scale factor $a_{\mathrm{crit}}$. In our numerical example FIG. 9, the critical scale factor is $a_{\mathrm{crit}}=5.118$ at which the lowest $\nu^{2}$ mode of the tensor perturbation shows $\mathcal{M}_{(G)}^{2}=\mathcal{H}^{2}_{(G)}$ and $\mathcal{M}_{(G)}^{2}$ is weaker than $\mathcal{H}_{(G)}^{2}$ for $a>a_{\mathrm{crit}}$. Thus, in order to prevent $\mathcal{H}$-accelerated instability, the initial value of the scale factor $a_{\mathrm{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, $\mathcal{O}^{[0;1]}$, $\mathcal{B}^{[1;1]}_{\mathrm{O}}$ for $a_{\mathrm{min}}\leq a\leq a_{\mathrm{max}}$, $\mathcal{O}^{[-1;1]}$, $\mathcal{O}^{[-1;-1]}$ and $\mathcal{O}^{[-1;-1]}_{\mathrm{BC}}$ for $a_{\mathrm{min}}\leq a\leq a_{\mathrm{max}}$. In this case, we also construct solutions whose squared effective masses are positive and dominate the $\mathcal{H}$-term within the possible ranges of the scale factor. A typical example is shown in FIG. 10.

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

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

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.

VI 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 (6). 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 $\Phi^{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 $\mathcal{C}=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 $\mathcal{F}$ and $\mathcal{G}$ showed in (53) and (59). $\mathcal{F}<0$ corresponds to the ghost instability which is equivalent to a lack of the lowest energy state. The case with $\mathcal{G}<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 $\mathcal{G}<0$ is suppressed by the effect of background dynamics. As we showed in (62), the stability of the perturbations can be judged by the sign of the effective friction coefficients $\mathcal{H}$ and the effective squared masses $\mathcal{M}^{2}$, additionally the magnitude relationships between $\mathcal{H}^{2}$ and $|\mathcal{M}^{2}|$. Thus, we introduced five types of stabilities and instabilities (i) negative-$\mathcal{M}^{2}$ instability, (ii) $\mathcal{H}$-accelerated instability, (iii) $\mathcal{M}^{2}$-dominated stability, (iv) $\mathcal{H}$-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 $\mathcal{M}^{2}$-dominated stability condition when $\mathcal{H}<0$ and the complete stability condition when $\mathcal{H}>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. V.2.3 is the stable solution which satisfy the $\mathcal{M}^{2}$-dominated stability condition when $\mathcal{H}<0$ and the complete stability condition when $\mathcal{H}>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 $\mathcal{M}^{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 methodref_Hill_eq which is summarized in Appendix A.

We numerically show a example of the oscillating solution with resonance instability in FIG. 12 and the coupling constants of the solution is shown in TABLE 4. Note that this solution is satisfied both $\mathcal{M}^{2}$-stability condition and complete stability condition (see FIG. 13). 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 $\epsilon_{\pm}$ defined in (109). If $|\epsilon_{\pm}|$ exceeds 1, the amplitude of the corresponding perturbation mode exponentially grows. In fact, the solution includes certain tensor perturbation modes with $|\epsilon_{\pm}|>1$ as we show in FIG.14. 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 2 as far as our numerical analysis is concerned. One may notice that such a manipulation corresponds to consider the small effective squared mass. See (LABEL:Gg) and (61), it is found that ultraviolet dominant terms $g_{8}\nu^{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.14 that the degree of enhancement $|\epsilon_{\pm}|$ approaches unity as the perturbation mode $\nu^{2}$ increases. We would like to stress that the resonance instability is mainly the problem in open FLRW spacetime. Since the perturbation mode $\nu^{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 $\mathcal{M}^{2}$ is preferred to avoid the instabilities due to the background dynamics, i.e., $\mathcal{H}$-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 perturbationU1_PHL1 . 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 conditionsU1_PHL2 ; U1_NPHL1 ; U1_NPHL2 . 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. V.2.1 are satisfied. Then, we can construct the stable singularity-free solutions only by examining $\mathcal{M}^{2}_{\mathrm{(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 $|\mathcal{M}^{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 $\nu^{2}=8$ in closed FLRW spacetime include the homogeneous modes, which corresponds the Bianchi-type IX spacetime whose anisotropy is small (in detail, see previous_stability_bounce ). 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).

Appendix A 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

[TABLE]

Then, the perturbation equations are transformed as

[TABLE]

with

[TABLE]

Since the effective squared angular frequency $\omega^{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

[TABLE]

Here $\hat{h}(t)$ is an arbitrary perturbation function and $X(t)$ is a real solution vector for the Hill equation (102).

Let $T$ be a period which is characterized by the effective squared angular momentum $\omega^{2}$, then, the evolution of $X$ is expressed as

[TABLE]

where, $\mathcal{A}$ is a $2\times 2$ matrix which represents the time translation $t\to t+T$. Suppose {$\epsilon_{+},X_{\epsilon_{+}}$} and {$\epsilon_{-},X_{\epsilon_{-}}$} are two independent eigensystems of matrix $\mathcal{A}$, where, $\epsilon_{\pm}$ are eigenvalues and $X_{\epsilon_{\pm}}$ are corresponding eigenvectors. Then, (108) can be rewritten as some linear combination of the following relations:

[TABLE]

Therefore, we find that both $|\epsilon_{+}|\leq 1$ and $|\epsilon_{-}|\leq 1$ are required to suppress the resonance of the perturbations.

To specify the explicit conditions for $|\epsilon_{\pm}|\leq 1$, we introduce two independent complex vector $X_{1}$ and $X_{2}$ whose initial conditions are

[TABLE]

Any solution vector can be constructed by the linear combination of $X_{1}$ and $X_{2}$. Substituting (114) and (119) into (108), we find

[TABLE]

Solving the Hill equation (102) with these initial conditions, $\mathcal{A}$ can be estimated. The eigenvalues of $\mathcal{A}$ are obtained from the following equation.

[TABLE]

Since the Hill equation (102) and the initial conditions (114) and (119) give $\mathrm{det}\mathcal{A}=1$, the eigenvalues $\epsilon_{\pm}$ are given by

[TABLE]

One can see that $|\mathrm{tr}\mathcal{A}|\leq 2$ is equivalent to $|\epsilon_{\pm}|=1$. When $\epsilon_{\pm}$ is a complex number, the perturbations correspond to the real part of $X(t)$. For the case with $|\mathrm{tr}\mathcal{A}|>2$, any one of $|\epsilon_{\pm}|$ is greater than unity. Thus, we conclude that if $|\mathrm{tr}\mathcal{A}|\leq 2$, the perturbation keeps oscillating forever without growing or decaying. On the other hand, the perturbation shows an exponential instability if $|\mathrm{tr}\mathcal{A}|>2$.