Bounce and Stability in the Early Cosmology with Anomaly-Induced Corrections

Abstract

An extremely fast exponential expansion of the Universe is typical for the stable version of the inflationary model, based on the anomaly-induced action of gravity. The total amount of exponential e-folds could be very large, before the transition to the unstable version and the beginning of the Starobinsky inflation. Thus, the stable exponential expansion can be seen as a pre-inflationary semiclassical cosmological solution. We explore whether this stable phase could follow after the bounce, subsequent to the contraction of the Universe. Extending the previous consideration of the bounce, we explore both stable expansion and the bounce solutions in the models with non-zero cosmological constant and the presence of background radiation. The critical part of the analysis concerns stability for small perturbations of the Hubble parameter. It is shown that the stability is possible for the variations in the bounce region, but not in the sufficiently distant past in the contraction phase.

1. Introduction

The importance of singularities in general relativity is partially based on the fact that they are unavoidable [1,2,3,4] and, therefore, indicate to the limits of applicability of the theory. The common belief is that, by properly extending general relativity, one arrives at the singularity-free theory. In particular, the cosmological singularities (see, e.g., [5] and references therein) led to various theoretical developments, such as, e.g., exploration of $f\left(R\right)$-type models with simplified form of quantum corrections [6,7], string cosmology [8], different versions of cyclic universe (see, e.g., [9]), quantum cosmology (see, e.g., [10,11,12]) and, in general, the interest in the bouncing models (see [13,14] for reviews and many references). Typically, bouncing models require fundamentally new assumptions, such as the especially designed non-localities (see, e.g., [15,16]) or introducing a scalar field with specially adjusted features (or a modified $f\left(R\right)$-type gravity, see, e.g., [17]). One of the most economical ways to avoid singularity is the assume a relevant positive curvature of the space section of the spacetime manifold [18] in the very early Universe.

The natural question is whether one can achieve the cosmological solutions with bounce with the minimal set of, or even without any, additional assumptions. The proposal of this type was done in the paper [19]. In the present paper, we elaborate the same idea in more detail and consider it from the modern perspective of the effective field theory framework.

Even without quantum gravity, the consistent treatment near the singularity requires taking into account the back reaction of the quantum matter fields on the classical gravitational background. In general, the derivation of the semiclassical corrections to the gravitational action is an unsolved problem, except for an extreme UV, where at least the one-loop effects are controlled by the logarithmic divergences, admitting a simple and elegant description in terms of conformal anomaly [20,21,22] and anomaly-induced effective action [23,24]. It was suggested that these quantum corrections may be sufficient to provide a bouncing solution [19], without additional assumptions. In the present work, we elaborate this proposal in detail. First, we extend the information about the bounce solution by including the large cosmological constant and the background radiation. These aspects of the model are relevant because a huge cosmological constant may be generated in the UV by symmetry restoration, while an intensive creation of matter from the vacuum is expected in the region of bounce, when the Universe leaves the de Sitter-type contraction phase. The second important point is the generalization of the phase diagrams for unstable [25,26] and stable [27,28] phases, in the presence of the cosmological constant. Finally, we explore the stability of the approximately exponential expansion and contraction with respect to the small perturbations of the conformal factor of the metric. Taken together, these new results extend our general understanding of the possibility of the bounce due to the semiclassical contributions to the gravitational action.

The paper is organized as follows. In Section 2, we briefly review the anomaly and induced action, and describe the derivation of the 00 component of the gravitational equations in the presence of radiation from the conservation law. In Section 3, we derive the particular de Sitter-like solutions from the 00 equations with semiclassical correction and formulates the detailed conditions of their stability. These conditions can be applied to both expanding and contracting phases. In Section 4, we present the phase diagrams with the cosmological constant for stable and unstable versions of the theory, and also give new examples of the bounce, including with the background radiation or cosmological constant. Section 5 is devoted to the stability analysis of the contracting case, clarifying the status of the bounce solutions. Finally, in Section 6, we draw our conclusions and also discuss general perspectives of solving the problem of cosmological singularities using quantum corrections.

The definitions of geometric quantities, used below, included the signature $(+,-,-,-)$ of the metric, the curvature tensor $R_{\tau \alpha \beta }^{\lambda }={\partial }_{\alpha }{\mathsf{\Gamma }}_{\tau \beta }^{\lambda }+\dots$, and the scalar curvature $R=g^{\alpha \beta }{R}_{\alpha \lambda \beta }^{\lambda }$.

2. Anomaly-Induced Action and the Early Universe

In the very early epoch of the Universe, the typical energy scale is very high. Different models of particle physics give different estimates for the upper bounds of the masses of the particles, but usually these upper bounds are of the order of magnitude, or below, the GUT scale $M_X\propto {10}^{16} \mathrm{GeV}$, that is a few orders below the Planck scale $M_P$. The Friedmann equation, based on general relativity, predicts an unbounded growth of the Hubble parameter near the beginning of time $t=0$, and the singularity emerges with $H\sim t^{-1}$. Thus, in the region close to the spacetime singularity, the gravitational field is so intensive, that all real particles can be regarded massless. The same concerns the virtual particles, which are interacting with the strong gravitational background. Thus, the relevant quantum effects in the vicinity of cosmological singularity may be limited by the effects of massless fields, which is, in fact, the simplest case.

It is a natural hypothesis that the fundamental theory of particle physics possesses the asymptotic freedom in the UV, including in the vicinity of the singularity. Furthermore, for fermion and vector fields, the masslessness implies local conformal invariance. This is a natural feature from the perspective of quantum theory. The point is that, free massless particles have a vanishing trace of the energy-momentum tensor. Demanding the field-particles duality in the description of quantum matter, we arrive at the condition for the trace of the dynamical energy-momentum tensor,

$$\begin{array}{c}\hfill T_{\mu }^{\mu }=-\frac{2}{\sqrt{-g}}{g}_{\mu \nu }\frac{\delta S}{\delta g_{\mu \nu }}=0.\end{array}$$

(1)

On the other hand, (1) is the Noether identity for the local conformal symmetry

$$\begin{array}{c}\hfill g_{\mu \nu }=e^{2\sigma }{\overline{g}}_{\mu \nu },{\mathsf{\Phi }}_i=e^{k_i\sigma }{\overline{\mathsf{\Phi }}}_i,\sigma =\sigma \left(x\right),\end{array}$$

(2)

where $k_i=(-1,-\frac{3}{2},0)$ is the conformal weight for the matter fields (scalars, fermions and vectors, ${\mathsf{\Phi }}_i=(\phi ,\psi ,A_{\mu })$, correspondingly. If the Noether identity (1) is satisfied, the action S does not depend on the conformal factor $\sigma \left(x\right)$. For the scalar field, the symmetry (2) imposes an additional condition on the value of the nonminimal parameter $\xi =\frac{1}{6}$ of interaction with the scalar curvature, such that the Lagrangian of this field (we write down the Lagrangian for a single real scalar, for simplicity), becomes

$$\begin{array}{c}\hfill {\mathcal{L}}_{scal}=\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi +\frac{1}{2}\xi R{\phi }^2.\end{array}$$

(3)

Furthermore, it is sufficient to choose the action of vacuum to satisfy the local conformal symmetry (2), without jeopardizing the renormalizability of the theory. Taking the power counting arguments into account, this means that the part of the vacuum action, which is the subject of an infinite renormalization in the UV, has the form

$$\begin{array}{c}\hfill S_{cv}=\int d^{4}x\sqrt{-g}\left\{a_1{C}^2+a_2{E}_4+a_3☐R\right\},\end{array}$$

(4)

where the integrand of the Gauss-Bonnet topological term and the square of the Weyl tensor are

$$\begin{array}{c}{E}_4=R_{\alpha \beta \mu \nu }{R}^{\alpha \beta \mu \nu }-4R_{\alpha \beta }{R}^{\alpha \beta }+R^2,\hfill \\ C^2=R_{\alpha \beta \mu \nu }{R}^{\alpha \beta \mu \nu }-2R_{\alpha \beta }{R}^{\alpha \beta }+\frac{1}{3}{R}^2,\hfill \end{array}$$

respectively. Action (4) satisfies the Noether identity (2) at the classical level. On the top of this, the classical action of vacuum may include the non-conformal part, e.g., the Einstein-Hilbert term with the cosmological constant $\mathsf{\Lambda }$, and the $R^2$ term,

$$\begin{array}{c}\hfill S_{ncv}=\int d^{4}x\sqrt{-g}\left\{-\frac{1}{16\pi G}(R+2\mathsf{\Lambda })+a_4{R}^2\right\}.\end{array}$$

(5)

Finally, $a_1,a_2,a_3,a_4$, $G$ and $\mathsf{\Lambda }$ are the parameters of the vacuum action. The renormalization of the parameters $a_1,a_2$ and $a_3$ generates the anomalous violation of the Noether identity (2). The main idea of the corresponding inflationary model is to use the anomaly, or the anomaly-induced action of gravity, to arrive at the cosmological solution which replaces the usual Friedmann-Lemaître solution of general relativity.

Consider the particle physics theory with $N_0$ real scalars (a complex scalar counts as two real scalars), $N_{1/2}$ Dirac spinors, and $N_1$ massless vector fields. The beta functions for the parameters $a_1,a_2,a_3$ are as follows:

$$\begin{array}{c}w={\beta }_1=\frac{1}{{\left(4\pi \right)}^2}\left(\frac{N_0}{120}+\frac{N_{1/2}}{20}+\frac{N_1}{10}\right),\hfill \end{array}$$

(6)

$$\begin{array}{c}b={\beta }_2=-\frac{1}{{\left(4\pi \right)}^2}\left(\frac{N_0}{360}+\frac{11N_{1/2}}{360}+\frac{31N_1}{180}\right),\hfill \end{array}$$

(7)

$$\begin{array}{c}c={\beta }_3=\frac{1}{{\left(4\pi \right)}^2}\left(\frac{N_0}{180}+\frac{N_{1/2}}{30}-\frac{N_1}{10}\right),\hfill \end{array}$$

(8)

where we introduced the useful notations for the coefficients $w$, $b$ and $c$. At the quantum level, for free matter fields, the Noether identity (2) is violated by the anomaly in the vacuum sector,

$$\begin{array}{c}\hfill T=\langle T_{\mu }^{\mu }\rangle =-\frac{2}{\sqrt{-g}}{g}_{\mu \nu }\frac{\delta {\overline{\mathsf{\Gamma }}}_{ind}}{\delta g_{\mu \nu }}=-\left(w{C}^2+b{E}_4+c☐R+\beta F^2\right).\end{array}$$

(9)

Here ${\overline{\mathsf{\Gamma }}}_{ind}$ is the renormalized quantum correction to the classical action (4) plus the contribution to the radiation term, $F^2=F_{\mu \nu }{F}^{\mu \nu }$. One has to include the background radiation term because in the early Universe there may be radiation. The radiation term is conformal invariant and can be kept together with $C^2$ term in the derivation of anomaly-effective action ${\overline{\mathsf{\Gamma }}}_{ind}$. The general form of the beta function is

$$\begin{array}{ccc}\hfill \beta & =& -{\displaystyle \frac{2e^2}{3{\left(4\pi \right)}^2}}{N}_f-{\displaystyle \frac{e^2}{6{\left(4\pi \right)}^2}}{N}_s,\hfill \end{array}$$

(10)

where $N_f$ and $N_s$ are multiplicities of charged fermions and scalars.

Solving Equation (9) with respect to ${\overline{\mathsf{\Gamma }}}_{ind}$ is a relatively simple exercise, described in many papers (e.g., in the first original works [23,24]), the review [29], book [30] and the forthcoming textbook [31]). Here we need just the simplest version of the solution, which can be obtained directly from the two relations

$$\begin{array}{c}\frac{\delta A}{\delta \sigma \left(x\right)}=\int d^{4}y\frac{\delta g_{\mu \nu }\left(y\right)}{\delta \sigma \left(x\right)}\frac{\delta A}{\delta g_{\mu \nu }\left(y\right)}=2g_{\mu \nu }\frac{\delta A}{\delta g_{\mu \nu }}\hfill \end{array}$$

and

$$\begin{array}{c}\hfill \sqrt{-g}\left(E_4-\frac{2}{3}☐R\right)=\sqrt{-g}\left({\overline{E}}_4-\frac{2}{3}\overline{☐}\overline{R}+4{\overline{\Delta }}_4\sigma \right).\end{array}$$

(11)

In these relations, we use the conformal parametrization of the metric (2). ${\Delta }_4$ is the fourth-order, Hermitian, conformal invariant operator [32,33,34] acting on a conformal invariant, dimensionless scalar field,

$$\begin{array}{c}\hfill {\Delta }_4={☐}^2+2R^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }-\frac{2}{3}R☐+\frac{1}{3}\left({\nabla }^{\mu }R\right){\nabla }_{\mu }.\end{array}$$

(12)

The general solution for the anomaly-induced effective action is (In this paper, we deal only with the dynamics of the conformal factor. Thus, this simplest non-covariant solution is equivalent to the nonlocal covariant one [23,24] and to the local covariant with two auxiliary fields [35].) [23,24]

$$\begin{array}{c}{\mathsf{\Gamma }}_{ind}=S_c\left[{\overline{g}}_{\mu \nu }\right]+\int d^{4}x\sqrt{-\overline{g}}\left\{\sigma \left(w{\overline{C}}^2+\beta {\overline{F}}^2\right)+b\sigma \left(\overline{E}-\frac{2}{3}\overline{☐}\overline{R}\right)+2b\sigma {\overline{\Delta }}_4\sigma \right\}\hfill \\ -\frac{3c+2b}{36}\int d^{4}x\sqrt{-g}{R}^2,\hfill \end{array}$$

(13)

where $\sqrt{-g}=\sqrt{-\overline{g}}{e}^{4\sigma }$ and $R=e^{-2\sigma }\left[\overline{R}-6{\left(\overline{\nabla }\sigma \right)}^2-6\overline{☐}\sigma \right]$. The first term in the solution (13) is an arbitrary conformal invariant functional, $S_c\left[{\overline{g}}_{\mu \nu }\right]=S_c\left[g_{\mu \nu }\right]$, which plays the role of an integration constant for Equation (9). There is no regular method for deriving $S_c\left[{\overline{g}}_{\mu \nu }\right]$ in the explicit form. However, in the known physical applications it is usually considered insignificant. E.g., for the background cosmological solutions (not for perturbations), this functional is irrelevant because an isotropic and homogeneous metric has the form

$$\begin{array}{c}\hfill g_{\mu \nu }={\overline{g}}_{\mu \nu }·a^2\left(\eta \right).\end{array}$$

(14)

As usual, the conformal time $\eta$ is related to the cosmic time t by the formula $a\left(\eta \right)d\eta =dt$. Thus, $S_c\left[{\overline{g}}_{\mu \nu }\right]$ is obviously independent of the conformal factor $a\left(\eta \right)$.

The form of Equations (9) and (13) demonstrates that the $☐R$ term in the anomaly contributes to the covariant local $R^2$ term in the induced action. The coefficient of this term is regularization-dependent [22,36]. This dependence was explored in detail in [37], where it was shown, in particular, that it is equivalent to modifying the finite $a_4$ term in the classical non-conformal action (5).

In contrast to the $C^2$ term in the action (4), the $R^2$ term affects the equation for the conformal factor $a\left(\eta \right)$. Therefore, it is important to fix the value of $a_4$ by some physical conditions. In the literature, there are two main choices for such a condition. The most important one is $a_4\approx 5\times {10}^8$, providing the consistent inflationary model of Starobinsky [25,26,38]. The large magnitude of this coefficient, compared to the beta function (8), makes the quantum contributions to the $R^2$ term completely irrelevant. Another interesting choice is $a_4\approx 0$, such that the particle contents $N_0$, $N_{1/2}$, $N_1$ in (8) becomes an important issue. In what follows, we shall follow this option and discuss how the choice of the particle contents defines the stability of the cosmological solutions.

The basis of the cosmological model for the early Universe is the total action, including (5), (4), matter and the quantum corrections described by (13)

$$\begin{array}{c}\hfill S_t=S_{vacuum}+S_{matter}+{\mathsf{\Gamma }}_{ind}.\end{array}$$

(15)

Let us remember that the matter corresponds to the pure radiation, its action is conformal invariant and thus, it does not affect the dynamical equation We shall see how the radiation density appears in the first integral of this equation.) for $a\left(\eta \right)$.

Taking variational derivative with respect to $a\left(\eta \right)$, we arrive at the equation of motion in the form

$$\begin{array}{c}\hfill T_{EH}+T_{HD}=T_{radiation},\end{array}$$

(16)

where the traces of the equations of motion for the corresponding terms in the action, are

$$\begin{array}{c}{T}_{EH}=\frac{3}{4\pi G}\left(\frac{a^{″}}{a^3}+\frac{k}{a^2}-\frac{2\mathsf{\Lambda }}{3}\right),\hfill \end{array}$$

(17)

$$\begin{array}{ll}{\overline{T}}_{HD}& =6c\left[\frac{4a^{‴}{a}^{\prime }}{a^6}-\frac{a^{\left(4\right)}}{a^5}+\frac{3{a^{″}}^2}{a^6}-\frac{6a^{″}{a^{\prime }}^2}{a^7}+2k\left(\frac{a^{″}}{a^5}-\frac{{a^{\prime }}^2}{a^6}\right)\right]\\ & +24b\left[\frac{a^{″}{a^{\prime }}^2}{a^7}-\frac{{a^{\prime }}^4}{a^8}+k\left(\frac{a^{″}}{a^5}-\frac{{a^{\prime }}^2}{a^6}\right)\right],\end{array}$$

(18)

$$\begin{array}{ccc}& & T_{radiation}=\frac{\beta {\overline{F}}^2}{a^4},\hfill \end{array}$$

(19)

where $k=0$ or $k=\pm 1$ for different space geometries.

Equation (16) represents the trace of the equations of motion, and can be used for analytical and numerical study of the system. However, it is more useful to derive 00-component ${\rho }_i$ using the conservation law for each of the components $T_{EH}$, $T_{HD}$ and $T_{radiation}$,

$$\begin{array}{c}\hfill d\left({\rho }_i{a}^3\right)=-p_{i}d\left(a^3\right),\mathrm{where}{T}_i={\rho }_i-3p_i.\end{array}$$

(20)

In what follows, we shall use the notations and general approach of the previous works, e.g., [39,40]. Let us stress that this consideration implies that we omit the conformal invariant functional $S_c\left[g_{\mu \nu }\right]$ in the purely gravitational part of the expression (13). This approximation is necessary here, because all the information we have, comes from the anomaly and reflects only the trace of the equations of motion. From Equation (20), follows that

$$\begin{array}{c}\hfill \frac{d{\rho }_i}{d{a}^3}+\frac{4}{3}\frac{{\rho }_i}{a^3}=\frac{T_i}{3a^3},\end{array}$$

(21)

with the general solution defined by

$$\begin{array}{c}\hfill {\rho }_i\left(a\right)=C\left(a\right)a^{-4},\mathrm{with}\frac{dC}{d\eta }=T_i{a}^3{a}^{\prime }.\end{array}$$

(22)

Integrating (22) for each of the components, we arrive at the following results:

$$\begin{array}{ll}{\rho }_{EH}& =\frac{3}{8\pi G}\left(\frac{{a^{\prime }}^2}{a^4}+\frac{k}{a^2}-\frac{\mathsf{\Lambda }}{3}\right)=\frac{3M_P^2}{8\pi }\left(\frac{{\dot{a}}^2}{a^2}+\frac{k}{a^2}-\frac{\mathsf{\Lambda }}{3}\right)\hfill \\ & =\frac{3M_P^2}{8\pi }\left({\dot{\sigma }}^2+k{e}^{-2\sigma }-\frac{\mathsf{\Lambda }}{3}\right),\hfill \end{array}$$

(23)

$$\begin{array}{ll}{\overline{\rho }}_{HD}& =\frac{6b\left({a^{\prime }}^4+2k{a}^2{a^{\prime }}^2\right)}{a^8}+\frac{3c\left(a{a^{″}}^2-2a{a}^{‴}{a}^{\prime }+4a^{″}{a^{\prime }}^2+2ka{a^{\prime }}^2\right)}{a^7}\hfill \\ & =(6b+9c)\frac{{\dot{a}}^4}{a^4}+6k(b+c)\frac{{\dot{a}}^2}{a^4}+3c\left(\frac{{\ddot{a}}^2}{a^2}-\frac{2\ddot{a}{\stackrel{.}{a}}^2}{a^2}-\frac{2\stackrel{⃛}{a}\dot{a}}{a^2}\right)\hfill \\ & =6b{\dot{\sigma }}^4+6k(b+c){\dot{\sigma }}^2{e}^{-2\sigma }+3c\left({\ddot{\sigma }}^2-2\stackrel{⃛}{\sigma }\dot{\sigma }-6\ddot{\sigma }{\dot{\sigma }}^2\right),\hfill \end{array}$$

(24)

$$\begin{array}{c}{\rho }_{radiation}=\frac{{\rho }_{r0}+\beta {\overline{F}}^{2}lna}{a^4}=({\rho }_{r0}+\beta {\overline{F}}^2\sigma )e^{-4\sigma }.\hfill \end{array}$$

(25)

In these expressions we have changed from conformal time $\eta$ to the cosmological time t, from $a\left(t\right)$ to $\sigma \left(t\right)=lna\left(t\right)$ and used the definition of Planck mass $M_P^2=1/G$. One of the advantages of the 00-component of the equations of motion compared to the trace is that, in the expression (25), we can include the classical radiation term. Since the corresponding term is conformal, it was hidden in the trace of Equation (16).

3. De Sitter-Like Solutions and Their Stability

Now we are in the position to write the dynamical equation for $\sigma \left(t\right)$ in the form

$$\begin{array}{c}\hfill {\rho }_{EH}+{\rho }_{HD}={\rho }_{radiation},\end{array}$$

(26)

where the ${\rho }_i$ are defined in Equations (23)–(25). Most of the further discussion will be based on this equation. For the sake of brevity, we restrict our attention to the $k=0$ case.

It is most economical to write the purely gravitational part of the equation in terms of the Hubble parameter $H=\dot{\sigma }$. In this way we get

$$\begin{array}{c}\hfill \frac{3M_P^2}{8\pi }{H}^2-\frac{\mathsf{\Lambda }{M}_P^2}{8\pi }+6b{H}^4+3c\left({\dot{H}}^2-2H\ddot{H}-6H^2\dot{H}\right)=({\rho }_{r0}+\beta {\overline{F}}^2\sigma )e^{-4\sigma }.\end{array}$$

(27)

As the first step, consider the equation without the radiation term. Looking for the solution with a constant H, we meet a biquadratic equation

$$\begin{array}{c}\hfill 6b{H}^4+\frac{3M_P^2}{8\pi }{H}^2-\frac{\mathsf{\Lambda }{M}_P^2}{8\pi }=0,\end{array}$$

(28)

with the following four solutions [41]:

$$\begin{array}{c}\hfill H=\pm \frac{M_P}{\sqrt{-32\pi b}}{\left(1\pm \sqrt{1+\frac{64\pi b}{3}\frac{\mathsf{\Lambda }}{M_P^2}}\right)}^{1/2}.\end{array}$$

(29)

It is clear that the solutions with the positive H correspond to expanding, and the ones with the negative H, to contracting Universes. From Equation (7) follows that $b<0$ for any particle contents of the underlying theory, such that the $\sqrt{-32\pi b}$ is always real. Consider the sign in the parenthesis, taking $H>0$ and $\mathsf{\Lambda }>0$ for the definiteness. Taking into account that $\mathsf{\Lambda }\ll M_P^2$, we get the two approximate solutions:

$$\begin{array}{c}\hfill H_1=\sqrt{\frac{\mathsf{\Lambda }}{3}}\mathrm{and}{H}_2=\frac{M_P}{\sqrt{-16\pi b}}.\end{array}$$

(30)

The first one, $H_1$, is the classical de Sitter solution where the anomaly-induced terms play no role. The second solution, $H_2$, is the solution corresponding to the equilibrium between quantum, anomaly-induced part and the classical part, of Equation (27). This solution cannot be obtained in the approach when the quantum part is treated a small perturbation. On the other hand, it can be derived from different approaches to quantum corrections in curved spacetime [42,43].

Let us consider the stability of these solutions. Our purpose is to explore the asymptotic stability, and, owing to the exponential time dependence of $a\left(t\right)$, it is important to perform variation of $\sigma \left(t\right)$ or $H\left(t\right)$, instead of $a\left(t\right)$ [41]. For the variation of $H\left(t\right)\to H\left(t\right)+X\left(t\right)$, without the radiation term in (27) and in the linear approximation, we obtain

$$\begin{array}{c}\hfill 24b{H}^{3}X+6c\dot{H}\dot{X}-6c\ddot{H}X-6c\ddot{X}H-18c\dot{X}{H}^2-36c\dot{H}HX+\frac{3M_P^2}{4\pi }HX=0.\end{array}$$

(31)

For a constant $H=H_0$, such as the ones in (30), this equation boils down to

$$\begin{array}{c}\hfill \ddot{X}+3H_0\dot{X}-\left(\frac{4b{H}_0^2}{c}+\frac{M_P^2}{8\pi c}\right)X=0.\end{array}$$

(32)

In the case of extreme inflationary solution with $H_0=H_2$ in (30), we have $-\frac{4b{H}_0^2}{c}=\frac{M_P^2}{4\pi c}$ and then the solution of (32) has the form

$$\begin{array}{c}\hfill X=C_1{e}^{{\lambda }_{1}t}+C_2{e}^{{\lambda }_{2}t},\end{array}$$

(33)

where

$$\begin{array}{c}\hfill {\lambda }_{1/2}=-\frac{3H_0}{2}\pm \frac{3H_0}{2}\sqrt{1+\frac{8b}{9c}}.\end{array}$$

(34)

Since we always have $b<0$, for $c>0$ there are no eigenvalues with the positive real part hence there are no growing modes and the solution with $H_0=H_2$ is stable [19,25,26]. Using the equation for the trace, one can show that in the model with $c<0$ the same solution is unstable.

It is instructive to consider the stability of the classical solution with $H_0=H_1$. In this case the equation for linear perturbation has the form

$$\begin{array}{c}\hfill \ddot{X}+3H_0\dot{X}-\frac{M_P^2}{8\pi c}X=0\end{array}$$

(35)

and the eigenvalues in (33) are

$$\begin{array}{c}\hfill {\lambda }_{1/2}=-\frac{\sqrt{3\mathsf{\Lambda }}}{2}\pm \frac{M_P}{\sqrt{8\pi c}}{\left(1+\frac{6\pi c\mathsf{\Lambda }}{M_P^2}\right)}^{1/2}\approx -\frac{\sqrt{3\mathsf{\Lambda }}}{2}\pm \frac{M_P}{\sqrt{8\pi c}}.\end{array}$$

(36)

Obviously, for $c>0$ there is a Planck order positive eigenvalue and hence a very fast growing mode. For $c<0$ there are oscillations with the frequency of the Planck order of magnitude, suppressed by a relatively slow damping. Qualitatively, the results for the stability are the same as derived earlier from the trace Equation (16) in [41].

The anomaly-induced action (13) includes the semiclassical corrections to the gravitational action which are certain, in the sense they do not depend on the way we intend to quantize gravity, use string theory, etc. Thus, we have to correctly interpret the stability results, formulated above. The Starobinsky model of inflation [25,36] is based on the unstable version and requires especially chosen initial conditions. As we have mentioned above, the phenomenological consistency requires adding a classical $R^2$ term with $a_4\approx 5\times {10}^8$ in the action (5). On the other hand, the stable version is free from the restrictions on the initial conditions.

As we just discussed, without the classical $a_4$ term, the stability is defined by the sign of c and this sign depends on the particle contents, according to Equation (8). The stability of $H_1$ solution is achieved for $N_1>18N_0+3N_{1/2}$. This inequality holds for the minimal standard model of particle physics (MSM) with $N_1=12$, $N_0=4$ and $N_{1/2}=24$, but does not hold for supersymmetric extensions, such as MSSM, or for the “doubled” standard model with $N_1=12$, $N_0=8$ and $N_{1/2}=48$. In what follows, we shall use the last version for the numerical analysis, but one has to remember that there is no real difference between the behavior of different models with $c>0$. Let us note that the $c=0$ case, which occurs for the $N=4$ supersymmetric Yang-Mills theory [44], gives an unstable solution with $H_2$.

It is clear that, in the late Universe, there should be either $c<0$, or a sufficiently large classical $a_4$, as otherwise flat space would be unstable. (As we already mentioned above, these two possibilities are equivalent [22,36,37].) Thus, using the advantage of the stable version at the beginning of inflation, requires the scheme of transition between stable and unstable particle contents. One of the possibilities for such a transition is related to the low-energy decoupling of the heavy degrees of freedom beyond the MSM, e.g., the s-particles in the case of MSSM [45]. A possible explanation of this decoupling is based on the slowing down exponential expansion with $a\left(t\right)\sim exp\left\{H_{2}t\right\}$ owing to the masses of the quantum s-particles [41,46]. It turns out that the process of slowing down inflation takes a long while and, as a result, the observable effects of inflation come from the very last 60–70 e-folds out of much greater rate of inflation.

The next question regarding this transition is how to provide the sufficiently large coefficient $a_4$ at the physically relevant, last 60–70, e-folds. This problem has been addressed in [28] and will not be discussed here. Instead of this, in the rest of the paper, we shall discuss the beginning of the exponential phase and, in particular, whether it could follow after the bounce from a contracting phase.

4. Phase Diagrams

Before starting the analysis of stability of the special solutions with respect to the small perturbations, let us complete the previous section by considering the phase diagrams for the Equation (27). Without the radiation term the equation is

$$\begin{array}{c}\hfill 2H\ddot{H}+6H^2\dot{H}-{\dot{H}}^2-\frac{2b}{c}{H}^4-\frac{M_P^2}{8\pi c}\left(H^2-\frac{\mathsf{\Lambda }}{3}\right)=0.\end{array}$$

(37)

We shall follow the analysis of the original paper by Starobinsky [25,26], and also the ones in [27,28] for the stable case. The main difference with these publications is that we include the cosmological constant term. As before, we restrict the main attention to the $k=0$ flat-space case. Indeed, there are no strong arguments of why the space curvature could be important in the very early Universe.

On the other hand, the cosmological constant may be relevant in the vicinity of the cosmological singularity because its value may be dramatically different from the one in the late Universe. The reason is that, in the last case, the observed value of the vacuum energy density ${\rho }_{\mathsf{\Lambda }}^{\left(obs\right)}=\frac{{\mathsf{\Lambda }}^{\left(obs\right)}}{16\pi G^{\left(obs\right)}}$ is a sum of the vacuum term ${\rho }_{\mathsf{\Lambda }}^{\left(vac\right)}$ and the contribution ${\rho }_{\mathsf{\Lambda }}^{\left(ind\right)}$, induced from the symmetry breaking in the matter fields sector, e.g., from the spontaneous breaking in the Higgs potential [47] (see also [48]). In the early epoch, this symmetry can be (and most likely is) restored owing to the high temperature or/and high curvature effects. As a result, the ${\rho }_{\mathsf{\Lambda }}^{\left(obs\right)}$ acquires the value which is typical for the corresponding phase transition. E.g., in case of the electroweak phase transition, this value is about 56 orders of magnitude greater than the one in the present Universe, but still negligibly small compared to the Planck scale in the vicinity of the singularity. However, in the case of the possible GUT-scale phase transition the difference between ${\rho }_{\mathsf{\Lambda }}^{\left(obs\right)}\propto M_X^4$ and $M_P^4$, may be just a few orders, and in the case of the Planck-scale phase transition the two quantities may be of the same order of magnitude. In what follows we consider the last two possibilities.

The order of Equations (37) can be reduced using the following change of variables [25,26]:

$$\begin{array}{c}\hfill x={\left(\frac{H}{H_0}\right)}^{\frac{3}{2}}\mathrm{and}y=\frac{\dot{H}}{2H_0^{\frac{3}{2}}}{H}^{-\frac{1}{2}}.\end{array}$$

(38)

Indeed, the reduction can be achieved by a simpler transformation $\dot{H}=x\left(H\right)$, but the resulting equation is more complicated. Equation (37) can be elaborated using (38) as

$$\begin{array}{c}\hfill \dot{H}=\frac{2}{3}{H}_0{x}^{-\frac{1}{3}}\dot{x},\mathrm{and}\ddot{H}=-\frac{2}{9}{H}_0{x}^{-\frac{4}{3}}{\dot{x}}^2+\frac{2}{3}{H}_0{x}^{-\frac{1}{3}}\ddot{x}.\end{array}$$

(39)

In this way, we arrive at

$$\begin{array}{c}\hfill \ddot{x}-\frac{2}{3}{x}^{-1}{\dot{x}}^2+3H_0{x}^{\frac{2}{3}}\dot{x}-\frac{3b}{2c}{H}_0^2{x}^{\frac{7}{3}}-\frac{3M_P^2}{32\pi c}\left(x-\frac{\mathsf{\Lambda }}{3}{H}_0^{-2}{x}^{-\frac{1}{3}}\right)=0.\end{array}$$

(40)

The second change of variables in (38) gives

$$\begin{array}{c}\hfill \dot{x}=3H_0{x}^{\frac{2}{3}}y\mathrm{and}\ddot{x}=2H_0{x}^{-\frac{1}{3}}\dot{x}y+3H_0{x}^{\frac{2}{3}}\dot{y},\end{array}$$

(41)

such that the last relation becomes

$$\begin{array}{c}\hfill \ddot{x}=6H_0^2{x}^{\frac{1}{3}}{y}^2+3H_0{x}^{\frac{2}{3}}\frac{dx}{dt}\frac{dy}{dx}=6H_0^2{x}^{\frac{1}{3}}{y}^2+9H_0^2{x}^{\frac{4}{3}}\frac{ydy}{dx}.\end{array}$$

(42)

Inserting these relations in (40), we obtain

$$\begin{array}{c}\hfill \frac{dy}{dx}=\frac{b}{6yc}\left[x+\frac{M_P^2{H}_0^{-2}}{16\pi b}\left(x^{-\frac{1}{3}}-\frac{\mathsf{\Lambda }}{3}{H}_0^{-2}{x}^{-\frac{5}{3}}\right)\right]-1.\end{array}$$

(43)

Taking $H_0=H_{+}$, defined in (29), we can write

$$\begin{array}{c}\hfill \frac{dy}{dx}=\frac{b}{6yc}\left[x-\frac{2}{\alpha }{x}^{-\frac{1}{3}}-\frac{64\pi b}{3{\alpha }^2{M}_P^2}\mathsf{\Lambda }{x}^{-\frac{5}{3}}\right]-1,\end{array}$$

(44)

where we used an abbreviation $\alpha \equiv \left(1+\sqrt{1+\frac{64b}{3}\frac{\mathsf{\Lambda }}{M_P^2}}\right)$.

The analogous formula from [28] can be reproduced in the limit $\mathsf{\Lambda }=0$. Using (44), we arrive to the phase diagrams for the two different values of $\mathsf{\Lambda }$, as shown in Figure 1 for the stable and in Figure 2 for the unstable versions of the theory (All numerical calculations and plots were done using Wolfram’s Mathematica [49].).

One can easily note that small cosmological constant does not produce essential modifications compared to the previously known $\mathsf{\Lambda }=0$ diagrams [25,26,27]. At the same time, the cosmological constant of the Planck order of magnitude, greatly transforms the phase diagram of the stable case making it qualitatively similar to the one of the unstable case with a small cosmological constant. This feature suggests a new mechanism of transition between the stable and unstable phases (different from the one described in [28,45,46]) and perhaps deserves further investigation in the future.

Imposing the special initial conditions at $t=0$, one can arrive at the bounce solutions, with the exponential [for a(t)] expansion after the exponential contraction, and the smooth transition between these two phases. The basic example of this sort is shown in Figure 3 for the model with zero cosmological constant. It is important that we cannot set $\dot{\sigma }\left(0\right)=0$, as otherwise the solution is $\sigma \left(0\right)\equiv 0$.

Let us present more details concerning the bounce solutions, based on the numerical analysis. Figure 4 demonstrates the result of the small variations of initial conditions. Here we still keep zero cosmological constant. These plots clearly show that the bounce survives such a small changes. However, we have found that larger modifications may exclude the bounce solutions.

Another two relevant aspects concern the choice of initial conditions. As we already know, changing these conditions, one can produce the solutions without bounce. It is especially important, that these conditions are imposed in the vicinity of the minimum $t=0$ of the curve. This point is in the past for the expansion region $t>0$, but, at the same time, it is in the future for the region $t<0$ where we meet a contraction. In the next section, we shall elaborate on the qualitative sense of this restriction.

The next step is to switch on the cosmological constant. As we have mentioned above, in the region close to singularity the magnitude of the corresponding density ${\rho }_{\mathsf{\Lambda }}^{\left(obs\right)}$ may be either a few (even many, actually) orders of magnitude smaller than $M_P^4$, or even of the same order of magnitude as $M_P^4$, if the Universe is approaching singularity after the corresponding phase transition. The example of the corresponding plots is shown in Figure 5. We can observe that the bounce may hold even for a relatively large cosmological constant. In fact, the cosmological constant of the given range does not change, qualitatively, the behavior near the singularity, even though the phase diagram gets modified. Let us note, without going into detail, that qualitatively the same situation with bounce holds for the $k=+1$ space geometry.

The last example of the numerical analysis is included to illustrate the bounce in the presence of radiation. This is an important aspect to explore, because in the region close to the bounce, the solution is rapidly deviating from de Sitter and, as a consequence, there is an intensive creation of particles and radiation from the vacuum (see, e.g., [50,51,52,53] and references therein). Thus, there may be a large amount of radiation in the vicinity of the bounce, and this should be taken into account. However, since we use the equation for the trace (16), the classical radiation is not “visible” for the equation. On the other hand, we can still model the presence of radiation with the anomalous contribution, and the result of this is shown in Figure 6.

5. Bounce and the Stability Analysis

It is interesting to explore the stability of the bounce solutions with respect to small perturbations of initial conditions. The analysis can be easily done analytically using the previous results (33) and (34) for the exponential behavior.

According to Equation (29), there are two pairs of solutions, corresponding to each of the values (30) with positive and negative signs. The “pre-inflationary” exponential expansion corresponds to the value of $H_2>0$. Then, from the diagrams in Figure 4, Figure 5 and Figure 6, follows that the contraction phase corresponds to the change of sign $H_2\to -H_2$.

Starting from this point, we can address the issue of stability in the linear contracting regime $\sigma \left(t\right)=-H_{2}t$ using the eigenvalues (34) of the solution (33) for small perturbations. Since we assume that $a_4$ in Equation (5) is zero, the value of c should be defined from quantum contribution (8). It is easy to see that $1+\frac{8b}{c}<0$ for the MSSM particle contents, when $c>0$. Then the real parts of both roots ${\lambda }_{1/2}$ in Equation (34) are defined by the sign of $H_0$. In the expansion phase, $H_0=H_2>0$, and both roots have negative real parts. On the contrary, in case of contraction, we have $H_0=-H_2<0$ and the solution $\sigma \left(t\right)\approx -H_{2}t$ is unstable under small perturbations of $\sigma \left(t\right)$. Thus, the stable bounce due to the anomaly-induced semiclassical corrections is impossible.

A natural and important question is how it happened that we meet the bounce solutions as shown in the plots above, if the bounce can not be stable? As we discussed in the previous section, there is no strong dependence on the initial conditions in these solutions. The numerical analysis shows that the bounce survives small variations of these conditions, which is typical for the stable case. We know that typically, for the linearly unstable solutions, one has to expect that the solution of the desired type can be achieved only for an absolutely precise fine-tuning of the initial conditions.

The answer to this question is as follows. Indeed, the stability or instability depend on the way we interpret the dynamics of the system or, in other words, depends on our viewpoint. Defining the initial conditions at the bounce point $t=0$ and choosing $|\dot{\sigma }\left(0\right)|$ sufficiently small, we are exploring the contraction phase by looking backward in time t. In this case, for the contraction there is a negative t and therefore a positive product $H_{0}t$. Then Equation (33) tells us that the contraction phase is stable. However, if we look forward in time, that is choose the initial point at some point $t_0<0$ near the same decreasing line, the solution is unstable.

It is worth noting that the situation in many other bounce models (see, e.g., the review [13]) is different. In these models there is a typical relation $H\sim \frac{1}{t}$, such that the sign of the product $Ht$ does not change in the transition between the expansion and contraction phases, when we consider the dynamics forward in time. (We are grateful to Nelson P. Neto for this explanation.)

6. Conclusions and Discussions

It is quite natural to expect that the anomaly-induced semiclassical corrections to the gravitational action cure the cosmological singularity. Qualitatively, the reason for this expectation is that the theory with fourth and higher derivatives has a natural cut-off. In the fourth-derivative theory [54] this cut-off is the Planck mass divided by the dimensionless parameter of the relevant coefficient of the fourth-derivative term. From the dimensional arguments follows that, in the the cosmological case, the higher derivative terms should modify the dynamics of the conformal factor of the metric and remove the singularity.

For example, in the Newtonian singularity case, i.e., for a point-like static mass source, many types of higher derivative terms remove singularity in the modified Newtonian potential. This occurs starting from the simplest fourth-derivative theory with the action composed by (5) and (4) terms [55]. The same is true for the local (see, e.g., [56,57]) and nonlocal [58] higher derivative models. It is worthwhile to note that the detailed analysis of the Newtonian and especially black hole cases is more complicated (see, e.g., [59,60,61,62] and further references therein).

In cosmology, the basic four-derivative action (5) and (4) does not remove singularity because the $C^2$ term, contributing with the ghost degree of freedom to the spin-2 sector of the theory, equals zero on the cosmological background. Then, the $R^2$ term alone cannot help removing singularity because it does not produce a massive unphysical ghost.

At this point, it becomes clear that the cosmological case is more complicated, especially because the situation with singularity can not be reduced to the analysis of the flat space propagator, as it is in the Newtonian singularity case. Taking these considerations into account, one can expect that the anomaly-induced term (13) may remove singularity and promote bounce since it introduces a more complete form of dynamics to the conformal factor $\sigma \left(t\right)$. This sort of intuitive arguments should be carefully verified by direct calculations. In the previous sections, we saw that the bounce solutions in this model really take place. On the other hand, these solutions are unstable if we define the initial conditions sufficiently far backward in the contraction regime and look forward in time. In this respect, the anomaly-induced model is different from the models with nonlocal form factors, where the bounce solutions can be found [15,16].

It would be certainly interesting to obtain a stable bounce solutions from the semiclassical or quantum gravity corrections, derived in a consistent way. One of the options is to explore higher-loop semiclassical corrections in the initially conformal model, which are expected to give higher than linear powers of logarithms in the UV form factors and, therefore, stronger non-linear $\sigma$-dependence in the effective action (see, e.g., [63]). This possibility represents an interesting proposal for a possible future work.