Cosmological perturbations and stability of nonsingular cosmologies with limiting curvature

TL;DR

This paper investigates the stability of nonsingular cosmologies with limiting curvature, identifying conditions to avoid instabilities, and introduces a new curvature invariant that enables stable early universe models.

Contribution

It proposes a novel curvature invariant scalar function that avoids Ostrogradski, ghost, and gradient instabilities, facilitating stable nonsingular cosmological scenarios.

Findings

Identifies the necessity of Weyl tensor squared in the action for stability.

Shows existing models suffer from ghost and gradient instabilities.

Introduces a new scalar invariant avoiding instabilities in most trajectories.

Abstract

We revisit nonsingular cosmologies in which the limiting curvature hypothesis is realized. We study the cosmological perturbations of the theory and determine the general criteria for stability. For the simplest model, we find generic Ostrogradski instabilities unless the action contains the Weyl tensor squared with the appropriate coefficient. When considering two specific nonsingular cosmological scenarios(one inflationary and one genesis model), we find ghost and gradient instabilities throughout most of the cosmic evolution. Furthermore, we show that the theory is equivalent to a theory of gravity where the action is a general function of the Ricci and Gauss-Bonnet scalars, and this type of theory is known to suffer from instabilities in anisotropic backgrounds. This leads us to construct a new type of curvature invariant scalar function. We show that it does not have Ostrogradski instabilities, and it avoids ghost and gradient instabilities for most of the interesting background inflationary and genesis trajectories. We further show that it does not possess additional new degrees of freedom in an anisotropic spacetime. This opens the door for studying stable alternative nonsingular very early universe cosmologies.

Equations

$$S=\frac{M_{\mathrm{Pl}}^2}{2}\int {\mathrm{d}}^{4}x\sqrt{-g}\left[R+\sum _i{\chi }_i{I}_i-V({\chi }_i)\right],$$

$$S=$$

$$\displaystyle+\chi_{2}I_{2}(\nabla,R^{\mu}{}_{\nu\rho\sigma})-V_{2}(\chi_{2})\Big{]}+S_{\mathrm{m}}~{},$$

$$I_1^{(0)}=12H^2,I_2^{(0)}=-6\dot{H}.$$

$$g_{\mu \nu }^{(0)}\mathrm{d}{x}^{\mu }\mathrm{d}{x}^{\nu }=-N(t{)}^2\mathrm{d}{t}^2+a(t{)}^2{\delta }_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j,$$

$$S^{(0)}=$$

$$+6(1-{\chi }_2)\frac{1}{N}\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\dot{a}}{Na}\right)-V_1-V_2]+S_{\mathrm{m}}^{(0)}.$$

$$12{\left(\frac{\dot{a}}{a}\right)}^2=\frac{\mathrm{d}{V}_1}{\mathrm{d}{\chi }_1}$$

$$-6\left[\frac{\ddot{a}}{a}-{\left(\frac{\dot{a}}{a}\right)}^2\right]=\frac{\mathrm{d}{V}_2}{\mathrm{d}{\chi }_2},$$

$$\frac{\epsilon }{3M_{\mathrm{Pl}}^2}=(1-2{\chi }_1-3{\chi }_2){\left(\frac{\dot{a}}{a}\right)}^2-{\dot{\chi }}_2\left(\frac{\dot{a}}{a}\right)-\frac{1}{6}(V_1+V_2),$$

$$-\frac{p}{M_{\mathrm{Pl}}^2}=$$

$$-2(2{\dot{\chi }}_1+3{\dot{\chi }}_2)\left(\frac{\dot{a}}{a}\right)-{\ddot{\chi }}_2-\frac{1}{2}(V_1+V_2).$$

$$12H^2=$$

$$-6\dot{H}=$$

$$\frac{\epsilon }{3M_{\mathrm{Pl}}^2}=$$

$$-\frac{1}{6}[V_1({\chi }_1)+V_2({\chi }_2)],$$

$$-\frac{p}{M_{\mathrm{Pl}}^2}=$$

$$-{\ddot{\chi }}_2-\frac{1}{2}[V_1({\chi }_1)+V_2({\chi }_2)].$$

$$\dot{a}$$

$${\dot{\chi }}_1$$

$${\dot{\chi }}_2$$

$$\frac{\mathrm{d}{\chi }_2}{\mathrm{d}{\chi }_1}$$

$$\frac{\mathrm{d}{\chi }_1}{\mathrm{d}a}$$

$$\frac{\mathrm{d}{\chi }_2}{\mathrm{d}a}$$

$$V_1({\chi }_1)=$$

$$V_2({\chi }_2)=$$

$$\frac{H^2}{H_{\max }^2}$$

$$\frac{\dot{H}}{H_{\max }^2}$$

$$V_1({\chi }_1)$$

$$V_2({\chi }_2)$$

$$\frac{H^2}{H_{\max }^2}$$

$$\frac{\dot{H}}{H_{\max }^2}$$

$$I_1$$

$$I_2$$

$$I_1^{(0)}=12H^2,I_2^{(0)}=-6{\dot{H}}^2,$$

$$g_{\mu \nu }^{(1)}\mathrm{d}{x}^{\mu }\mathrm{d}{x}^{\nu }=-\mathrm{d}{t}^2+a^2\left({\delta }_{ij}+h_{ij}+\frac{1}{2}{h}_{ik}{h}^k{}_j\right)\mathrm{d}{x}^i\mathrm{d}{x}^j,$$

$$h_{ij}=\int \frac{{\mathrm{d}}^3\vec{k}}{(2\pi {)}^3}\left[h_{\vec{k}}^{+}{e}_{ij}^{+}(\vec{k})+h_{\vec{k}}^{\times }{e}_{ij}^{\times }(\vec{k})\right]{\mathrm{e}}^{\mathrm{i}\vec{k}\cdot \vec{x}},$$

$$S_T^{(2)}=$$

$$\displaystyle-\mathcal{M}_{T}\frac{k^{2}}{a^{2}}h_{\vec{k}}^{I}h_{-\vec{k}}^{I}\Big{)}~{},$$

Full text

I Introduction

One of the biggest problems with the classical theory of general relativity is the occurrence of singularities, which are inevitable under realistic assumptions [1, 2, 3] and which signify the breakdown and the incompleteness of the theory. The big bang singularity in cosmology and the singularity at the center of a black hole are two well-known instances of singularities in general relativity that one would like to resolve. These singularities often find themselves in regions of high density, high energy, and high curvature, where one may expect the breakdown of the classical theory and the emergence of quantum behavior. For this reason, there is hope that a quantum theory of gravity would provide the resolution to the otherwise pathological classical singularities.

Without a proper theory of quantum gravity, one may approach the problem with an effective theory that could mimic the low-energy behavior of the full quantum theory. The effective theory could be constructed with one or more new degrees of freedom, e.g. a new scalar field. This allows one to study nonsingular theories of the very early Universe within effective field theory (EFT) [4, 5, 6] as is done in, for instance, bouncing cosmologies [7, 8, 9, 10] or genesis scenarios [11, 12, 13, 14, 15, 16] with a generalized scalar field such as in Horndeski [17] and generalized Galileon [18] theories, whose equivalence was first proved in [19]. Alternatively, one may attempt to modify the Einstein-Hilbert gravity action to include higher-order curvature terms (e.g., [20, 21, 22]). Interestingly, this can serve as the basis for implementing the limiting curvature hypothesis, which seeks to incorporate the idea of a fundamental limiting length (that would be realized in the full theory of quantum gravity) into the effective theory for gravity.

In line with special relativity where the speed of light is bounded from above and quantum mechanics where the Heisenberg uncertainty relation holds, the idea of the limiting curvature hypothesis comes from the fact that one may expect quantum gravity to possess a fundamental length scale $\ell_{f}$ below which no measurement can be made and on which all physical observables must be smeared out. Presumably, this fundamental scale is at least of the order of the Planck length $\ell_{\mathrm{Pl}}=\sqrt{\hbar G/c^{3}}$, although it could be larger. Taking $\ell_{f}\sim\ell_{\mathrm{Pl}}$, it is straightforward to see that if all curvature invariants are bounded throughout the spacetime manifold ($|R|<\ell_{\mathrm{Pl}}^{-2}$, $|R_{\mu\nu}R^{\mu\nu}|<\ell_{\mathrm{Pl}}^{-4}$, $|\nabla_{\rho}R_{\mu\nu}\nabla^{\rho}R^{\mu\nu}|<\ell_{\mathrm{Pl}}^{-6}$, $|C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}|<\ell_{\mathrm{Pl}}^{-8}$, etc., where $R_{\mu\nu}$ is the Ricci tensor, $R$ is the Ricci scalar, $C_{\mu\nu\rho\sigma}$ is the Weyl tensor, and $\nabla$ is the covariant derivative), then the spacetime is nonsingular. Indeed, in the well-known cases of the big bang and black hole singularities, some of the physically measurable curvature invariants such as $R$, $R_{\mu\nu}R^{\mu\nu}$, and $C^{2}:=C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}$ blow up; hence finding theories in which all invariants are bounded is certainly a necessary condition for constructing nonsingular cosmologies.

Unfortunately, bounding an infinite number of curvature invariants is rather nontrivial. Indeed, there are well-known instances where low-order curvature invariants are bounded while higher-order invariants are still unbounded (e.g., $|R_{\mu\nu}R^{\mu\nu}|<\ell_{\mathrm{Pl}}^{-4}$ while $|\nabla_{\rho}R_{\mu\nu}\nabla^{\rho}R^{\mu\nu}|\rightarrow\infty$). This is where the limiting curvature hypothesis comes in handy. The hypothesis states that [23, 24, 25, 26, 27] if one finds a theory that allows a finite number of curvature invariants to be bounded by an explicit construction (e.g., $|R|\leq\ell_{\mathrm{Pl}}^{-2}$ and $|R_{\mu\nu}R^{\mu\nu}|\leq\ell_{\mathrm{Pl}}^{-4}$), and when these invariants take on their limiting values, then any solution of the field equations reduces to a definite nonsingular solution (e.g., de Sitter space), for which all curvature invariants are automatically bounded. We note, though, that the assumptions of the limiting curvature hypothesis generally do not ensure that solutions avoid singularities when curvature scalars are not on their limiting values.

The limiting curvature hypothesis has been used and tested in the context of black hole physics [27, 26, 28, 29, 30, 31, 32, 33, 34, 35]. In this context, the geometry outside the black hole horizon is described by the usual Schwarzschild metric, but inside the event horizon, the black hole singularity is replaced with a nonsingular de Sitter spacetime, which, in turn, could be the source of a new “baby” Friedmann universe. Similarly, in a cosmological context [36, 37, 38, 30, 31, 39], a nonsingular universe can be constructed, in vacuum, such that it is asymptotically de Sitter in the past and Minkowski in the future (or vice versa). This is in line with Penrose’s vanishing Weyl tensor conjecture [40] (see also the discussion in Ref. [41]), which states that the Weyl tensor should vanish at the beginning of the universe, since de Sitter space has $C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}=0$. With the addition of matter sources, one can obtain asymptotically de Sitter and Friedmann cosmologies, remaining nonsingular throughout cosmic time. Recent works also show that the ideas of limiting curvature could allow one to construct nonsingular bouncing cosmologies [42, 43, 44].

In this paper, we want to revisit nonsingular cosmological models that make use of an effective theory for gravity in which the limit curvature hypothesis is realized. It was shown in Refs. [36, 37] that interesting background cosmologies can be found within this framework by constructing a theory in which the curvature invariants $R$ and $4R_{\mu\nu}R^{\mu\nu}-R^{2}$ are bounded. However, these studies did not explore the cosmological perturbations [45] for the action containing the above curvature invariants. Recent developments in nonsingular cosmology within EFT [4, 5, 6] have shown that it is often rather difficult to avoid instabilities in the cosmological perturbations (e.g., see Refs. [46, 47, 48] for no-go theorems within Galileon and Horndeski theories; see, also, Refs. [49, 50, 51]). For this reason, one could tend to believe that nonsingular models constructed as in Refs. [36, 37] are going to be very unstable at the perturbation level, thus rendering the models unviable.

In this work, we will show that the naive models of Refs. [36, 37] are indeed generically unstable. We will see that minimal extensions in which one also includes the Weyl tensor squared, $C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}$, in the curvature invariants are more robust, i.e. there are fairly large regions of parameter space that are stable. Yet, there do not seem to exist nonsingular cosmological solutions that remain stable throughout cosmic history, and moreover, the theory will be shown to be equivalent to a $f(R,\mathcal{G})$ theory of gravity (where $\mathcal{G}$ is the Gauss-Bonnet term), which has unavoidable ghosts [52]. We will then construct a completely new curvature-invariant function and show that it allows for stable nonsingular cosmological solutions throughout time. There will remain some difficulties though in constructing a physically relevant model in certain cases.

The outline of this paper is as follows. In Sec. II, we will review the construction of a nonsingular cosmology in which the limiting curvature hypothesis is realized, set up the action for the theory, and find the background equations of motion. In particular, we will discuss two specific scenarios: an inflationary scenario and a genesis scenario. We will then study the cosmological perturbations, determine the general stability conditions and check them for specific models. The equivalence with $f(R,\mathcal{G})$ gravity will also be demonstrated. In Sec. III, we will construct a new model with a new curvature scalar, derive the resulting cosmological perturbations and the stability conditions, and comment on the case of an anisotropic background. We will then discuss the recovery of vacuum Einstein gravity and Friedmann cosmology with the addition of matter sources. We will end with a summary of the results and a discussion in Sec. IV.

II Nonsingular cosmology with limiting curvature

II.1 Setup of the theory and background evolution

The approach taken in Refs. [36, 37] to implement the limiting curvature hypothesis consists of introducing a finite number of nondynamical scalar fields $\chi_{i}$, or Lagrange multipliers, such that the action takes the form

[TABLE]

where the $I_{i}$’s are functions of curvature invariants that can depend on $R$, $R_{\mu\nu}$, $R^{\mu}{}_{\nu\rho\sigma}$, $C_{\mu\nu\rho\sigma}$ and combinations and derivatives thereof. Accordingly, given an appropriately chosen potential $V(\chi_{i})$, one can rewrite this action into a general $F(R,R_{\mu\nu}R^{\mu\nu},...)$ effective theory of gravity. By virtue of the Lagrange multipliers, a given potential imposes constraints on the $I_{i}$’s, hence the idea that the right choice of $V(\chi_{i})$ can naturally bound the curvature invariants and satisfy the limiting curvature hypothesis asymptotically. We will give examples where this is realized below.

As is done in Refs. [36, 37], we are going to consider two nondynamical scalar fields and start with a general action of the form

[TABLE]

where $S_{\mathrm{m}}$ is the action for possible matter sources. At this point, we do not make any assumption on the functional form of $I_{1}$ and $I_{2}$, but we want them to scale as $R$, so let us require that we recover a certain limit at the background level:

[TABLE]

The superscript $(0)$ refers to the metric of a flat111It is straightforward to generalize this to include curvature (see Ref. [37]). Friedmann-Lemaître-Robertson-Walker (FLRW) universe,

[TABLE]

where $N$ is the lapse function and $a$ is the scale factor. Accordingly, $H:=\dot{a}/(Na)$ is the Hubble parameter, and a dot is a derivative with respect to physical time, $t$. At the background level, we can further ask that $\chi_{1}=\chi_{1}^{(0)}(t)$ and $\chi_{2}=\chi_{2}^{(0)}(t)$. The original action then becomes222We omit the superscript $(0)$ for $\chi_{1}$ and $\chi_{2}$ when it is clear that they represent background quantities.

[TABLE]

Varying $S^{(0)}$ with respect to $\chi_{1}$ and $\chi_{2}$ yields the equations of motion (EOMs)

[TABLE]

and

[TABLE]

respectively, where we set the lapse function to $N=1$. Letting $T^{\mu}{}_{\nu}=\mathrm{diag}(-\varepsilon(t),p(t)\delta^{i}_{\ j})$, where $T_{\mu\nu}$ is the stress-energy tensor associated with the matter action $S_{\mathrm{m}}$ and where $\varepsilon$ is the energy density and $p$ is the pressure, one can then vary the background action with respect to $N$ to find

[TABLE]

again setting $N=1$. Finally, varying with respect to $a$ gives (setting $N=1$ once more)

[TABLE]

Let us denote $V_{1}^{\prime}:=\mathrm{d}V_{1}/\mathrm{d}\chi_{1}$ and $V_{2}^{\prime}:=\mathrm{d}V_{2}/\mathrm{d}\chi_{2}$ for shorthand notation from here on. We can then summarize the set of EOMs as

[TABLE]

In the limit where $\chi_{1}=\chi_{2}=V_{1}=V_{2}=0$, we note that we recover the usual Friedmann equations, as expected. Also, in the limit where $\varepsilon=p=0$, one obtains the EOMs in vacuum. Demanding that $V_{1}^{\prime}(\chi_{1})>0$ for all $\chi_{1}$ values so that $H$ is real and looking at an expanding universe (so $H>0$; this could be generalized to a contracting universe with $H<0$, in which case a minus sign would appear in certain equations), we can write the EOMs as

[TABLE]

Furthermore, the equations can be written in the following form:

[TABLE]

II.1.1 Example of inflationary scenario

To solve the background EOMs, one needs to specify a form for the potentials $V_{1}(\chi_{1})$ and $V_{2}(\chi_{2})$. As a first example, one can consider

[TABLE]

which is inspired from Ref. [37], and as we will see, this gives rise to a nonsingular inflationary scenario. In vacuum ($\varepsilon=p=0$), the EOMs are given by

[TABLE]

We plot these functions in Fig. 1. As we can see in the left-hand plot, the Hubble parameter is finite everywhere: as $\chi_{1}\rightarrow 0$, the spacetime is asymptotically Minkowski ($H\rightarrow 0$; recall that we are in vacuum), whereas when $\chi_{1}\rightarrow\infty$, the spacetime is asymptotically de Sitter since $H\rightarrow H_{\mathrm{max}}$. Similarly, looking at the right-hand plot, $\dot{H}$ is finite everywhere, and it is asymptotically vanishing as $\chi_{2}\rightarrow\pm\infty$. Accordingly, this verifies the limiting curvature hypothesis as in Ref. [37].

We note at this point that since we regard our theory as a low-energy effective theory of a possible quantum theory of gravity, there should be a cutoff scale beyond which the EFT is no longer valid. The model here includes only two dimensionful parameters: $M_{\mathrm{Pl}}$ and $H_{\mathrm{max}}$. Therefore, the cutoff scale should naively be given by these parameters as $\Lambda_{\mathrm{cut}}=(M_{\mathrm{Pl}}H_{\mathrm{max}}^{n})^{1/(1+n)}$ for a given integer $n\neq-1$, and in particular, it should be at least of the order of $H_{\mathrm{max}}$. Determining the exact value for $\Lambda_{\mathrm{cut}}$ involves a rather nontrivial computation for the given theory, but since the energy scale of our cosmological solutions is always less than $H_{\text{max}}$ by construction, the validity of EFT is naturally ensured.

The phase-space diagram for the model is plotted in Fig. 2, where the arrows show the vectors with components $(\dot{\chi}_{1},\dot{\chi}_{2})$ computed from Eqs. (15) and (16) (in vacuum with $\varepsilon=0$) with the potentials (20) and (21). We highlight a specific trajectory in green for illustration. In this case, the universe starts asymptotically at $\chi_{1}\rightarrow\infty$ and $\chi_{2}\rightarrow 0$ and ends asymptotically at $\chi_{1}\rightarrow 0$ and $\chi_{2}\rightarrow 0$, so as we saw from Fig. 1, the universe starts in a de Sitter spacetime and ends in a Minkowski spacetime.333With the addition of matter sources, it would end in a FLRW spacetime as shown in Ref. [37].

At this point, one may wonder how the given scenario evades the singularity theorems of [53, 54, 55] regarding the past incompleteness of inflationary cosmology. First, it is important to recall that the singularity theorem for inflation [53, 54] is proved under the assumption that gravity is given by the Einstein-Hilbert action and that inflation is driven by matter obeying the null energy condition. Our higher derivative gravity terms, when taken to the matter side of the equations of motion, act as matter violating the null energy condition. Hence the theorem does not apply in our setup. We also note that some of the past directed geodesics would have finite affine length (in agreement with the situation in Ref. [55]), but this is simply due to the fact that the flat FLRW chart does not cover the entire de Sitter space. One can extend the spacetime so that all geodesics are complete as in the case of de Sitter space. Thus, our inflationary universe is free from initial singularities.

II.1.2 Example of genesis scenario

As another example, let us consider

[TABLE]

In vacuum, the EOMs become

[TABLE]

We plot these functions in Fig. 3. As we can see in the left- and right-hand plots, the Hubble parameter and its time derivative are again everywhere finite: as $\chi_{1}\rightarrow 0$ or $\chi_{1}\rightarrow\infty$, the spacetime is asymptotically Minkowski with $H\rightarrow 0$, and $\dot{H}\rightarrow 0$ as $\chi_{2}\rightarrow\pm\infty$. We note that $H_{\mathrm{max}}$ is now reached when $\chi_{1}\rightarrow 1$. Thus, this is another type of scenario that verifies the limiting curvature hypothesis, namely a genesis scenario, which starts in Minkowski space rather than de Sitter space.

The phase-space diagram for this model is plotted in Fig. 4, where again, the arrows show the vectors with components $(\dot{\chi}_{1},\dot{\chi}_{2})$ computed from Eqs. (15) and (16) (in vacuum with $\varepsilon=0$) with the genesis potentials (24) and (25). We highlight different trajectories in green, red, black, and purple for illustration. All of these curves either start or end at $\chi_{1}\rightarrow\infty$, which corresponds to Minkowski spacetime. However, the green and red curves are pathological trajectories since they either end or start at $\chi_{1}\rightarrow 1$, at which point it can be shown that $V_{1}^{\prime\prime}\rightarrow 0$. Accordingly, from Eq. (15), one finds $\dot{\chi}_{1}\rightarrow\pm\infty$ at that point. More interestingly, the black and purple curves start and end at $\chi_{1}\rightarrow\infty$, and they turn around at some minimal $\chi_{1}>1$ value, so they never reach the “singularity” at $\chi_{1}=1$. Also, these trajectories always have $\chi_{2}\ll 1$. When looking at the left-hand plot of Fig. 3, these trajectories suggest that the universe starts in the far right at $\chi_{1}\rightarrow\infty$ (Minkowski spacetime), rolls up to the left but does not reach $H_{\mathrm{max}}$, and rolls back down to Minkowski spacetime again. In light of a structure formation scenario for the very early universe, one would like to have some form of reheating444For instance, reheating could occur via gravitational particle production [56, 57] (see Refs. [16, 58] for examples of gravitational particle production in nonsingular cosmologies). near the maximal value that the Hubble parameter reaches. Thus, the universe would start as Minkowski spacetime, but it would end as a radiation- and then matter-dominated FLRW spacetime.

II.2 Cosmological perturbations and stability analysis

We now turn to the study of the cosmological perturbations for the action given by Eq. (II.1). At this point, one needs to specify the form of the curvature-invariant functions $I_{1}$ and $I_{2}$. Motivated by Refs. [36, 37], let us take

[TABLE]

where at this point $\kappa$ is just some real constant. In a flat FLRW background, these curvature invariants reduce to

[TABLE]

under the assumption555We note that the requirement $\dot{H}<0$ restricts our attention to the regions of phase space in which $\chi_{2}<0$ for the examples given in Secs. II.1.1 and II.1.2. that $\dot{H}<0$, which was the hypothesis of Eq. (3) that allowed us to find the general background EOMs in Sec. II.1. We note that the background expressions for the curvature invariants do not depend on the constant $\kappa$ since flat FLRW spacetime is conformally flat, so the term proportional to the Weyl tensor squared does not affect the dynamics of background spacetime.

II.2.1 Tensor modes

Let us begin by studying the tensor fluctuations. We start by perturbing the metric linearly as

[TABLE]

where the perturbation tensor $h_{ij}$ is transverse and traceless, i.e. $h^{i}{}_{i}=\partial_{i}h^{i}{}_{j}=0$ (adding the last term on the right-hand side does not change the linear equations but simplifies the derivation). We define the Fourier components of $h_{ij}$ by

[TABLE]

where $\{e_{ij}^{+},e_{ij}^{\times}\}$ represents the polarization basis. Given the curvature-invariant functions of Eq. (II.2), we can now perturb Eq. (II.1) to second order with the above metric to find

[TABLE]

where

[TABLE]

and the coefficients $\mathcal{K}_{T}$ and $\mathcal{M}_{T}$ will be specified shortly. We pause here to note that, in general, since the second-order action in the tensor sector has nondegenerate higher-derivative terms ($\propto\ddot{h}^{2}$), there appear to be ghost degrees of freedom according to Ostrogradski’s theorem. The only way to avoid those Ostrogradski ghosts would be if $\mathcal{G}_{T}$ were to vanish identically ($\mathcal{G}_{T}\equiv 0$). We see that this is not possible for a generic real constant $\kappa$, but if one sets $\kappa=-2$, then the model is safe with regards to Ostrogradski instabilities. The original models of Refs. [36, 37] did not include $C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}$ in their curvature invariants at all, so they had $\kappa=0$. Accordingly, the above implies that these models are inherently unstable. Yet, the addition of the Weyl tensor squared in the invariants with a specific prefactor ($\kappa=-2$) avoids this conclusion while having no effect on the background evolution. Still, it does not mean that the theory is necessarily free of all types of instabilities as we will see shortly.

The other coefficients of Eq. (II.2.1) are given by

[TABLE]

The criteria to avoid ghost and gradient instabilities are $\mathcal{K}_{T}>0$ and $\mathcal{M}_{T}>0$, respectively. By using the background EOMs (see Eq. (14), which can be rewritten as $H=\sqrt{V_{1}^{\prime}/12}$, Eq. (15), and Eq. (16) in the case for vacuum with $\varepsilon=0$), the conditions can be written solely in terms of the fields $\chi_{1}$, $\chi_{2}$ and their potentials $V_{1}(\chi_{1})$, $V_{2}(\chi_{2})$ as

[TABLE]

These conditions will be tested for the different models of Secs. II.1.1 and II.1.2 shortly.

II.2.2 Vector modes

We shall consider vector fluctuations in the following gauge, where

[TABLE]

Here, the vector perturbations $\beta^{i}$ satisfy $\partial_{i}\beta^{i}=0$. The Fourier components of $\beta_{i}$ are then defined by

[TABLE]

where $\{e^{1}_{i},e^{2}_{i}\}$ are orthogonal spatial vectors perpendicular to $\vec{k}$. The second-order action for vector perturbations becomes

[TABLE]

where $\mathcal{G}_{V}=\mathcal{G}_{T}$ and $\mathcal{K}_{V}=\mathcal{K}_{T}$. Accordingly, when $\kappa=-2$, which sets $\mathcal{G}_{T}=0$ to avoid Ostrogradski ghosts, it turns out that $\mathcal{G}_{V}=0$ as well, and as a result, there are no dynamical vector modes.

II.2.3 Scalar modes

We shall then focus on the scalar fluctuations in the spatially flat gauge, where $\chi_{1}=\chi_{1}^{(0)}+\delta\chi_{1}$, $\chi_{2}=\chi_{2}^{(0)}+\delta\chi_{2}$, and

[TABLE]

The second-order action for scalar modes is then given by

[TABLE]

where $\mathcal{G}_{S}=\mathcal{G}_{T}$ and666We omit the subscript $\vec{k}$ from the perturbation variables $\Psi^{I}$ when it is clear that they represent the Fourier components. $\Psi^{I}:=(\Phi,B,\delta\chi_{1},\delta\chi_{2})$. The matrix $M_{IJ}$ is given by777As before, we omit the superscript $(0)$ for $\chi_{1}$ and $\chi_{2}$ when it is understood that they represent background quantities.

[TABLE]

where $*$ stands for symmetric components. Since no time derivatives of $\Phi$, $B$, and $\delta\chi_{1}$ appear in the second-order action with $\kappa=-2$, these modes are nondynamical. Then, these variables can be eliminated by their equations of motion. After removing the nondynamical modes, the resulting action can be written solely in term of $\delta\chi_{2}$ as follows:

[TABLE]

where

[TABLE]

with

[TABLE]

We do not write down the form of the other $C_{n}$ coefficients because they are not so relevant in the following stability analysis.

In the small scale limit ($k/a\rightarrow\infty$), one finds

[TABLE]

where $A_{S}$ and $B_{S}$ are simply two real positive constants. Thus, the ghost instability in the scalar sector is avoidable when it is absent in the tensor sector, i.e. when ${\cal K}_{T}>0$ is satisfied. In addition, the gradient instability is absent when ${\cal M}_{S}>0$, so when

[TABLE]

By using the vacuum background EOMs, this condition can be rewritten as

[TABLE]

Regarding gradient instabilities which can occur on sub-Hubble scales ($k/a\gg H$), the general procedure is to check for the stability of modes within the validity range of the EFT, i.e. for $H\ll k/a\leq\Lambda_{\mathrm{cut}}$ (see for instance Ref. [10]), where for our models $H\leq H_{\mathrm{max}}\leq\Lambda_{\mathrm{cut}}$. The condition given by Eq. (51) can be viewed as the one which ensures that the shortest wavelength modes ($k/a\sim\Lambda_{\mathrm{cut}}$) do not suffer from gradient instabilities. However, since the perturbed action exhibits a modified dispersion relation, i.e., since Eq. (46) has terms of order $k^{2}$, $k^{4}$, $k^{6}$, and $k^{8}$, longer wavelength modes (longer than $\Lambda_{\mathrm{cut}}^{-1}$, but still smaller than $H^{-1}$) could still suffer from gradient instabilities. As long as the duration for such gradient instabilities is not too long though, their amplification remains controllable in comparison to the smaller wavelength modes which easily blow up (within a time scale of the order of $\Lambda_{\text{cut}}^{-1}$). This is why we only focus on the stability of the shortest wavelength modes.

II.2.4 Stability analysis

In summary, with $\kappa=-2$, we saw that there is no Ostrogradski instability. Then, we derived three conditions given by Eq. (36), which determines when the model is free of ghost instabilities in the tensor and scalar sectors, and Eqs. (37) and (51), which determine when the model is free of gradient instabilities in the tensor and scalar sector, respectively. The conditions depend on the potentials $V_{1}(\chi_{1})$ and $V_{2}(\chi_{2})$ and on the field values $\chi_{1}$ and $\chi_{2}$, so we need to study specific models to comment on the stability of the given theory.

Starting with the inflationary model of Sec. II.1.1, where the potentials are given by Eqs. (20) and (21), we plot the regions of phase space that satisfy the three conditions in Fig. 5. The individual conditions are shown in the left-hand plot, and we see that there are large regions of phase space that can avoid ghost or gradient instabilities in the tensor or scalar sectors. However, when we look at the right-hand plot, which shows the combined region where all stability conditions are met, we see that there is, actually, only a very small region of phase space that is not unstable. In particular, the trajectories that correspond to asymptotically de Sitter and Minkowski (e.g., the green curve in Fig. 2) are generally unstable throughout their evolution, except in a very small region of phase space.

We show the same types of plots in Fig. 6 for the genesis scenario of Sec. II.1.2. There as well, there are large regions of phase space that can avoid ghost and gradient instabilities in the tensor or scalar sectors, but it remains that only small portions of those can be stable with regards to all types of instabilities at the same time. In particular, the interesting trajectories (e.g., the black and purple curves of Fig. 4) are unstable throughout their evolution.

II.3 Equivalence with $f(R,\mathcal{G})$ gravity

We have shown that the Weyl square term with $\kappa=-2$ kills the dangerous degrees of freedom and can relax the instability around a flat FLRW background spacetime. This result can be understood because this theory is included by $f(R,\mathcal{G})$ gravity, whose cosmological solutions are perturbatively stable. However, it is shown by Ref. [52] that the perturbations of flat FLRW spacetime in $f(R,\mathcal{G})$ gravity lose a degree of freedom, which reappears as a ghost around anisotropic spacetimes. Therefore, our theory necessarily suffers from the same instability.

Let us demonstrate the equivalence between the above model and $f(R,{\cal G})$ gravity. The Gauss-Bonnet term ${\cal G}$ is defined by

[TABLE]

Plugging the relation

[TABLE]

into the definition of the Gauss-Bonnet term, we obtain

[TABLE]

Then, we note that the second curvature-invariant function of Eq. (II.2) can be written as (with $\kappa=-2$)

[TABLE]

Accordingly, we see that $I_{2}$ and $I_{1}$ (which are given by $I_{1}=I_{2}+R=\sqrt{R^{2}-6{\cal G}}+R$) are functions of $R$ and ${\cal G}$. Thus, the solutions of the EOMs for $\chi_{1}$ and $\chi_{2}$ can be written as

[TABLE]

and by plugging these solutions into Eq. (II.1), the original action becomes only a function of $R$ and the Gauss-Bonnet term ${\cal G}$, i.e.

[TABLE]

In conclusion, this theory is a specific model of $f(R,{\cal G})$ gravity.

III New nonsingular model with a new curvature scalar

III.1 Setup

Let us investigate other curvature scalars, which reduce to $H$ and $\dot{H}$ at the background level. In particular, we focus our attention on curvature scalars that are functions of the Ricci scalar and its derivatives.

Let us consider the following tensor constructed from the first derivative of $R$,

[TABLE]

For a flat FLRW background, this quantity reduces to

[TABLE]

The trace of the tensor $X^{\mu}{}_{\nu}$ is

[TABLE]

and at the background level, it reduces to

[TABLE]

Since the second time derivative of the Hubble parameter appears here, we need to consider another curvature scalar to remove $\ddot{H}$. Thus, let us consider the second derivative of $R$,

[TABLE]

whose expression in a FLRW spacetime reduces to

[TABLE]

Since $\dddot{H}$ appears only in the $(0,0)$ component, it can be removed by the tensor defined by

[TABLE]

whose FLRW limit now is $P^{\mu}{}_{\nu}^{(0)}=\mathrm{diag}(0,1,1,1)$. In fact, we can construct the scalar quantity $\mathrm{tr}[P\mathcal{R}P\mathcal{R}]$, which gives

[TABLE]

at the background level. Therefore, the following quantity is a good candidate as a limiting curvature-invariant function,

[TABLE]

which satisfies

[TABLE]

as required. In order to bound $\dot{H}$, we introduce as before

[TABLE]

hence

[TABLE]

again as required.

III.2 Cosmological perturbations and stability analysis

III.2.1 Tensor modes

By substituting the definition of tensor perturbations [Eq. (30)], the second-order action for tensor modes can be obtained as

[TABLE]

where

[TABLE]

Thus, the tensor sector does not include higher-derivative terms, so the number of degrees of freedom is 2, and no Ostrogradski instabilities appear. The stability conditions for ghost and gradient instabilities in the tensor sector are then given by

[TABLE]

III.2.2 Vector modes

Using Eq. (38), the second-order action for vector modes can be derived as

[TABLE]

with $\mathcal{K}_{V}=\mathcal{K}_{T}$. Therefore, no vector modes exist in this theory.

III.2.3 Scalar modes

Let us investigate the scalar perturbations defined in Eq. (41). Let us introduce a new perturbation variable $\varphi$ in terms of $\Phi$ and $B$ by the equation

[TABLE]

and we regard the components of $\Psi^{I}:=(B,\delta\chi_{1},\delta\chi_{2},\varphi)$ as independent variables.

The second-order action for scalar modes can then be calculated as

[TABLE]

where we omit writing down the specific form of the terms $K$, $L_{I}$, and $M_{IJ}$. Since the action does not include the time derivative of $B$, $\delta\chi_{1}$, and $\delta\chi_{2}$, we can remove these variables by the EOMs. The resulting action is given by

[TABLE]

where the coefficients, in the limit where $k/a\rightarrow\infty$, are given by

[TABLE]

with

[TABLE]

III.2.4 Stability analysis

In summary, we derived four conditions given by Eqs. (75)-(76) and Eqs. (III.2.3)-(84), which determine when the model is free of ghost and gradient instabilities in the tensor and scalar sectors, respectively. Once again, the conditions depend on the potentials $V_{1}(\chi_{1})$ and $V_{2}(\chi_{2})$ and on the field values $\chi_{1}$ and $\chi_{2}$, so we need to study specific models to comment on the stability of the given theory.

Starting with the inflationary model of Sec. II.1.1, we plot the regions of phase space that satisfy the four conditions in Fig. 7. Once again, there are large regions of phase space that can avoid ghost or gradient instabilities in the tensor or scalar sectors. The major difference with Fig. 5 is apparent in the right-hand plot of Fig. 7, which shows the combined region where all stability conditions are met. Indeed, whereas the previous theory had only a small strip of phase space that was stable, there is now a large portion of the phase space that is free of all types of instabilities. Furthermore, it is precisely in this region that we obtained the interesting background trajectories starting from de Sitter spacetime and ending in Minkowski spacetime (e.g., the green curve in Fig. 2). For this class of trajectories, there only seems to be a small region of phase space around $0\lesssim\chi_{1}\lesssim 2$ and $-0.5\lesssim\chi_{2}\lesssim-2.5$ where there still appears to be some instability. Yet, it seems clear that the curvature invariants given by Eqs. (66) and (70) lead to a much more stable theory compared to the curvature invariants of Eq. (II.2) that were analyzed in Sec. II.2.

For the genesis model of Sec. II.1.2, the stable phase-space regions are shown in Fig. 8. There as well, there now are larger regions of phase space that can avoid all types of instabilities. In particular, the interesting trajectories (e.g., the black and purple curves of Fig. 4) are stable throughout their evolution. This reinforces our conclusion that the limiting curvature theory with curvature invariants (66) and (70) is more stable and leads to interesting nonsingular cosmological scenarios that can remain stable throughout their evolution.

III.3 Stability around an anisotropic background

At this point, the new curvature-invariant function leads to generally stable nonsingular solutions around isotropic backgrounds. Still, one may worry that this new theory might still possess undesired ghosts. Indeed, it appears that the new curvature scalar of this section is not included in the ghost free theories found in Ref. [59], which can be further mapped into multifield extensions of Horndeski theories [60, 61, 62].

One situation where new ghosts may appear is around anisotropic backgrounds. In the model of Sec. II, which we showed to be equivalent to a $f(R,{\cal G})$ theory of gravity, a degree of freedom is lost around a FLRW background and such a mode reappears as a ghost once the background is deformed to an anisotropic Bianchi type I universe. The purpose of this subsection is to investigate the stability of our new theory against such a nonperturbative deformation of the background spacetime. We shall investigate the perturbations around a Bianchi type I universe with the rotational symmetry in the $y-z$ plane, so the background metric is given by

[TABLE]

At the background level, the curvature scalar invariants now include shear terms,

[TABLE]

where the shear $\sigma$ is defined by

[TABLE]

and so, the background dynamics is different from that analyzed in Sec. II.1.

We focus on whether or not additional degrees of freedom of perturbations appear independently of the background dynamics. One would usually start with the action given by Eq. (II.1) with curvature invariants given by Eqs. (66) and (70) and perturb the action according to the scalar-vector-tensor decomposition of the perturbed metric given by Eqs. (30), (38), and (41). Following the methodology of [63] for perturbations around an anisotropic background, thanks to the axisymmetry of the background spacetime given by Eq. (85), the linear perturbations can actually be decomposed into scalar and vector modes only, with respect to rotations of the $y-z$ plane,888Since there is no spherical symmetry in an anisotropic background, the usual (three-dimensional) scalar-vector-tensor modes cannot develop independently. However, thanks to one rotational symmetry, the one in the $y-z$ plane, a (two-dimensional) scalar-vector decomposition works well. For example, the $x$ component of a three-dimensional vector is just a scalar, and the $y$ and $z$ components can be decomposed into one scalar (gradient) mode and one vector (transverse) mode. The reason for the absence of (two-dimensional) tensor modes is just that $2\times 2$ symmetric tensors can be written as two scalar modes and one vector mode only. In this sense, the ten components of the metric tensor can be decomposed into seven scalar modes and three vector modes generally. Since three scalar modes and one vector mode can be killed by gauge symmetry, we have four scalar modes and two vector modes in the metric tensor.

[TABLE]

The Fourier components of $\delta g_{\mu\nu}^{\mathrm{scalar}}$ and $\delta g_{\mu\nu}^{\mathrm{vector}}$ are given by

[TABLE]

[TABLE]

where $k$ is defined by

[TABLE]

Here we have fixed the arbitrariness of rotations in the $y-z$ plane so that $\vec{k}=(k_{x},k_{y},0)$. We have also fixed the gauge degrees of freedom of perturbations. Our gauge choice corresponds to the spatially flat gauge [defined by Eqs. (30), (38), and (41)] in the limit where $\delta\rightarrow 0$ because the perturbed metric reduces to

[TABLE]

with

[TABLE]

The rotation matrix $U(\vec{k})$ is defined by

[TABLE]

and it transforms the vector $(0,k_{x},k_{y},0)$ into $(0,k,0,0)$.

Since the time derivative of $B_{\vec{k}}$ only appears through the following combination in the second-order action,

[TABLE]

it is useful to regard $\varphi$ as a dynamical variable instead of $\Phi$, analogous to Eq. (78). Thus, we now have six scalar mode perturbations ($\varphi$, $B$, $\beta_{1}$, $h_{+}$, $\delta\chi_{1}$, and $\delta\chi_{2}$), and two vector mode perturbations ($\beta_{2}$ and $h_{\times}$). By a straightforward calculation, one can show that the second-order action does not include any time derivatives of $B$, $\beta_{1}$, $\delta\chi_{1}$, $\delta\chi_{2}$, and $\beta_{2}$ after integration by parts. Thus, these variables are nondynamical and the remaining dynamical degrees of freedom are $\varphi$ and $h_{+}$, which are scalar modes, and $h_{\times}$, which is a vector mode. This result shows that no additional degrees of freedom appear at least in this anisotropic background. Therefore, our new model is not disturbed by the deformation of the background spacetime given by Eq. (85). This is a crucial difference compared to the model of Sec. II or $f(R,{\cal G})$ gravity.

III.4 Recovering Einstein gravity and the addition of matter sources

At this point, it looks like the action (II.1) with curvature invariants (66) and (70) can lead to interesting nonsingular cosmological background models such as an inflationary model and a genesis model that would remain stable against all types of instabilities throughout most of their evolution. These models start in de Sitter or Minkowski spacetime, respectively, and both end up in Minkowski spacetime. To be viable structure formation scenarios for the very early universe, as we already pointed out, one needs a reheating mechanism that would produce radiation and matter in a sufficient amount after the universe has acquired its adiabatic, scale-invariant curvature perturbation power spectrum. At this point, let us suppose that such a reheating mechanism exists and that it produces matter and radiation in large amounts. Then, at the background level, one is left with the nonsingular theory described in Sec. II.1 with nonzero energy density, and possibly, nonzero pressure as well. For the theory to successfully describe our universe, it is then necessary to recover the usual Einstein equations, i.e. the Friedmann equations in our context.

For the inflationary scenario, this is not a problem. Once matter is included, and as $\chi_{1}$ and $\chi_{2}$ (and their time derivatives) go to 0, one notes from Eqs. (20) and (21) that $V_{1}\rightarrow 0$ and $V_{2}\rightarrow 0$. Thus, we see that Eqs. (II.1) and (II.1) reduce to the Friedmann equations and can lead to the expected radiation- and matter-dominated era of our universe.

In the context of the genesis scenario, one runs into difficulty though. Once reheating has occurred and matter is included, one remains in the regime $\chi_{2}\ll 1$, but $\chi_{1}\rightarrow\infty$. From Eqs. (24) and (25), this still implies $V_{1}\rightarrow 0$ and $V_{2}\rightarrow 0$. Then, simply at the level of the action (II.1), one notices that recovering the Hilbert-Einstein term alone is only possible if $\chi_{1}I_{1}\rightarrow 0$, which is to say that $I_{1}$ vanishes faster than $\chi_{1}\rightarrow\infty$. However, since $I_{1}\sim\mathcal{O}(R)$, it is not possible to have a nonzero Hilbert-Einstein term while $I_{1}\rightarrow 0$. One can also see that, for Eq. (II.1) to be valid as $\chi_{1}\rightarrow\infty$, the only possibility is that $H$ and $\varepsilon$ vanish faster than $\chi_{1}\rightarrow\infty$. Once again, this implies an empty Minkowski spacetime rather than a FLRW spacetime. Accordingly, it seems impossible, in the context of this genesis scenario, to have the higher-derivative terms from the curvature invariants vanish, i.e. to recover the Einstein equations, and be left with a nonempty Friedmann universe. Therefore, the genesis scenario within this theory remains at the level of a toy model.

It should be noted that nontrivial couplings between matter fields and $\chi_{1}$ or $\chi_{2}$ may relax this problem [64]. In this case, the bounds on the curvature are, however, weakened because $I_{1}$ and $I_{2}$ include matter fields. Such matter couplings must be subdominant at high energy scales in order to ensure the avoidance of curvature singularity, but they must be dominant at low energy scales in order to recover Einstein gravity.

IV Conclusions and discussion

In this paper, we revisited the nonsingular cosmologies of Refs. [36, 37], which implement the limiting curvature hypothesis. We extended the analysis beyond the background cosmology to include the linear cosmological perturbations and determined the criteria for stability. This showed that the original models of Refs. [36, 37] appear to have, generically, undesired additional degrees of freedom leading to Ostrogradski instabilities. These instabilities could be killed with the addition of the Weyl tensor squared in the curvature-invariant functions given the appropriate coefficient. Still, by exploring two nonsingular cosmological scenarios in which the limiting curvature hypothesis is realized (one inflationary and one genesis scenario), it appeared that the cosmologies inevitably possess either ghost or gradient instabilities through large portions of their evolution. Furthermore, we showed that the theory could be rewritten as a $f(R,\mathcal{G})$ theory of gravity, which is known to suffer from instabilities in anisotropic backgrounds.

We then constructed a new curvature-invariant function by taking a specific combination of covariant derivatives of the Ricci scalar. Given the same inflationary and genesis scenarios at the background level as before, we showed that the new curvature scalar could lead to stable cosmologies with respect to Ostrogradski ghosts, as well as ghost and gradient instabilities throughout most of their evolution. Furthermore, the theory does not possess additional new degrees of freedom around anisotropic backgrounds, contrary to $f(R,\mathcal{G})$ gravity.

In light of constructing a nonsingular theory for the very early universe, there remain some challenges though. If one starts in a vacuum universe (either de Sitter or Minkowski), one would need to provide some form of reheating mechanism, possibly via gravitational particle production, so that the universe can contain matter and radiation after the early epoch. Furthermore, one would need to assure that the theory reduces to the Einstein limit for gravity fast enough. This appears to be satisfied in our inflationary scenario, but it remains an issue in the genesis scenario. Finally, given a successful scenario at the background level and which is stable perturbatively, it would be straightforward to solve the perturbation equations to find the power spectra of these perturbations and compare with observations to validate the theory.

The analysis performed in this paper opens the window to construct and study other nonsingular cosmologies, e.g., bouncing cosmologies. As mentioned before, this has been explored in Ref. [42] (see also Refs. [43, 44]). However, other nonsingular models with limiting curvature [35, 42] could also suffer from instabilities as suggested by Ref. [65]. This might be due to the fact that Refs. [35, 42] implement the limiting curvature hypothesis within a mimetic theory [66, 67, 68, 69], whose stability (or instability) does not appear to have been settled yet (see, e.g., Refs. [70, 71, 72, 73, 69, 74]).

Finally, it would be interesting to study how the approach to implement the limiting curvature hypothesis used in this paper fits in the grand picture of general scalar-tensor theories of gravity (e.g., [73, 75]). In particular, it would be interesting to find general classes of curvature-invariant functions in which the limiting curvature hypothesis can be realized and where solutions are stable.

Acknowledgements.

We are grateful to Lavinia Heisenberg and David Langlois for valuable discussions. D. Y. is supported by the Japan Society for the Promotion of Science (JSPS) Postdoctoral Fellowships for Research Abroad. J. Q. wishes to thank the Tokyo Institute of Technology and RESCEU at the University of Tokyo for kind hospitality while this work was prepared and acknowledges financial support from the Walter C. Sumner Memorial Fellowship and from the Vanier Canada Graduate Scholarship administered by the Natural Sciences and Engineering Research Council of Canada (NSERC). M. Y. wishes to thank the Institute for Theoretical Studies (ITS) of the ETH Zürich and McGill University for kind hospitality while this work was initiated and prepared. M. Y. is supported in part by the JSPS Grant-in-Aid for Scientific Research Nos. 25287054 and 26610062 and by the MEXT Grant-in-Aid for Scientific Research on Innovative Areas “Cosmic Acceleration” No. 15H05888. R. B. wishes to thank the ITS of the ETH Zürich for a Senior Fellowship during the 2015/2016 academic year when this project was initiated and acknowledges financial support from an NSERC Discovery Grant, from the Canada Research Chair program, from a Simons Foundation fellowship, and (at the ITS) from Dr. Max Rössler, the Walter Haefner Foundation and the ETH Zürich Foundation. This work was also partially supported by the Open Partnership Joint Projects Grant of JSPS.