text:
# 1 Introduction.

## 1 Introduction.

The importance of inflation becomes more and more apparent in the last years (see, e.g.  for the review). One of the reasons is that the existing and future astrophysical data may give a chance to use the early universe as a laboratory for the high energy physics, providing the Theoretical Physics a new important information. However, in order to use this opportunity, one needs a model of inflation which should be free of any serious phenomenological input. Otherwise, the experimental data will mainly deal with this input and not with the underlying high energy physics.

In this paper, we are going to discuss the example of such a model - the new version  of the anomaly-induced inflation (Starobinsky model) . This model is based on the well known methods of quantum field theory in curved space-time (see, e.g. ) and on the effective quantum field theory approach. The advantage of the anomaly-induced inflation is that it does not contain usual elements of the conventional cosmological phenomenology such as inflaton, the precise form of its potential and/or the necessity of a fine-tuning for the initial data. The anomaly-induced inflation provides the explanation for both the start of inflation and the graceful exit to the post-inflationary FRW-like phase . Under a natural supposition about the supersymmetry breaking scale we arrive at the magnitude of the Hubble constant which does not generate unacceptable amplitudes of the gravitational waves .

The problems of stability play a great role in the anomaly-induced inflation. In particular, the transition between the stable and unstable regimes, which could take place due to the breaking of supersymmetry  and the quantum vacuum effects of massive fields , could be responsible for the graceful exit to the FRW evolution. In the next sections, we investigate the stability with respect to the small perturbations of the conformal factor for various regimes associated to the anomaly-induced inflation. One can distinguish three stages of the expansion:

$`i)`$ For the almost exponential inflation that takes place at the very beginning (when the effect of particle masses is negligible) we consider the spaces with zero, positive and negative space curvatures, and with zero, positive or negative cosmological constant. The main part of the results can be obtained analytically such that some very important general features are established.

$`ii)`$ At the post-inflationary stage, after the supersymmetry breaking, the universe is approaching the stable fixed point with a small Hubble constant. In this case the FRW solution is a very good approximation for the theory with vacuum quantum corrections, such that these corrections do not play an essential role at the later time. We shall verify that the FRW-like solution is an attractor for the theory with quantum corrections if the cosmological constant has positive sign.

$`iii)`$ The most complicated is the intermediate transitional epoch between the inflationary and the post-inflationary evolutions. This period is characterized by the vacuum quantum effects of both massless and massive fields. Indeed, the exact analytic methods for the stability analysis  do not apply in this case, and one has to use numerical and approximate analytical methods.

The paper is organized as follows. In the next section, we present basic concepts of the anomaly-induced inflation. In part, we go beyond the results of the previous publications . In particular, we take into account the effect of the renormalization of the cosmological constant (here, our approach is somehow similar to the one of ). In section 3 the stability under the small perturbations of the conformal factor is analyzed for the “massless” stage $`i)`$, mainly through the analytical methods. The main new aspect here is that we clarify the role of the cosmological constant. In section 4 we verify, using both approximate analytical and numerical methods, that the stability of the model with the supersymmetric particle content really holds until the scale of the supersymmetry breaking. In section 5 we discuss, using the approximate analytical methods tested in section 4 and also the numerical integrations, the stability under the metric perturbations. Finally, in section 6 we draw our conclusions and also present a general discussion about the possible theoretical significance of the model.

## 2 Brief survey of the anomaly-induced inflation

Consider the vacuum quantum effects in the Early Universe, below the Planck scale. Then, the appropriate framework is not a string theory but some effective quantum field theory on a classical curved background. Let us make some assumptions about a quantum field theory which is valid at the sub-Planck energies. It proves useful to concentrate on the asymptotically free theories, that helps to avoid the discussion of the higher-loop effects. Thus, in what follows we shall deal with the one-loop approximation and only mention possible effects of the higher loops.

Consider the theory which has $`N_0`$ real scalars, $`N_{1/2}`$ Dirac spinors, and $`N_1`$ vectors. One has to notice that the vacuum quantum effects originate from the virtual particles. Therefore, the numbers $`N_0,N_{1/2},N_1`$ correspond to the particle content of the quantum theory, but they do not describe the real matter which might fill the Universe. As we shall see later on, the real matter has little importance for the anomaly-induced inflation. For simplicity, one can assume that the initial universe is empty and that all matter content of the Universe is created during the reheating period .

The renormalizable quantum theory of matter fields on curved background includes the action of vacuum, and its minimal necessary form is

$`S_{vacuum}=S_{HD}+S_{EH}.`$ (1)
Here the first term contains higher derivatives of the metric

$`S_{HD}={\displaystyle d^4x\sqrt{g}\left\{a_1C^2+a_2E+a_3^2R+a_4R^2\right\}}`$ (2)
and the second is the Einstein-Hilbert term 444Our notations are $`\eta _{\mu \nu }=diag(+)`$ and $`R_{\mu \nu }=_\lambda \mathrm{\Gamma }_{\mu \nu }^\lambda \mathrm{}.`$ The sign of the cosmological term in (3) corresponds to the positive vacuum energy for positive $`\mathrm{\Lambda }`$.

$`S_{EH}={\displaystyle \frac{1}{16\pi G}}{\displaystyle d^4x\sqrt{g}(R+2\mathrm{\Lambda })},`$ (3)
where $`a_1,..,a_4`$, $`G`$ and $`\mathrm{\Lambda }`$ are the parameters of the vacuum action. The renormalization of these parameters and the quantum corrections to the action (1) are the basics of the inflationary model we are discussing here. The discussion in this section will perform in two stages: first we review the massless limit and then present a useful approximation  for the massive case.

### 2.1 Massless conformal case

At high energies the masses of the fields are negligible and massless conformal fields serve as a good approximation. The vacuum quantum effects correspond to the Feynman diagrams with a loop of the quantum matter field and external lines of the classical gravity. Then, in order to evaluate the relevance of the mass of the quantum field, one has to compare it to the energy corresponding to the external line. In this paper we shall follow Ref.  and associate the energy of a quanta of the gravitational field with the magnitude of the Hubble constant $`H`$. Therefore, the massless approximation is justified for the very large values of $`H`$. As we shall see in the next subsection, this corresponds to the initial stage of inflation.

The actions of the free real scalar, Dirac fermion and gauge vector fields are

$`S_s={\displaystyle \frac{1}{2}}{\displaystyle d^4x\sqrt{g}\left\{(\phi )^2+m_s^2\phi ^2+\xi R\phi ^2\right\}},`$ (4)
where $`(\phi )^2=g^{\mu \nu }_\mu \phi _\nu \phi `$;

$`S_{1/2}=i{\displaystyle d^4x\sqrt{g}\left\{\overline{\psi }[\gamma ^\mu _\mu +im_f]\psi \right\}};`$ (5)
and

$`S_1={\displaystyle d^4x\sqrt{g}\left\{\frac{1}{4}F_{\mu \nu }F^{\mu \nu }\right\}}.`$ (6)
The massless matter fields possess local conformal invariance in the case where scalars couple to gravity in a nonminimal conformal way, with $`\xi =1/6`$. In what follows, we consider only this value, because it provides the traceless stress tensor for the massless theory at the classical level. It is well known that at the one-loop level the condition $`\xi =1/6`$ holds under renormalization (see e.g. ), but at higher loop this is not the case . However, even at higher loops the deviation from $`1/6`$ can be very small , such that $`\xi =1/6`$ can serve as a good approximation. Furthermore, massless fermions are conformal invariant, also the gauge fields can be treated as massless and conformal at very high energies.

The renormalization of the Einstein-Hilbert, cosmological and $`\sqrt{g}R^2`$ terms in (1) is absent in the case of conformal massless fields, and from this point of view their presence is not necessary. Indeed, we have to introduce the Einstein-Hilbert term because it is proved to fit with all experimental and observational data at the cosmic scale. For the same reason we introduce the cosmological term which was detected in the recent observations . Regarding the $`\sqrt{g}R^2`$ term, we have to accept the presence of this term since it is important for the renormalizability at higher loops. But its coefficient $`a_4`$ may be safely taken very small  such that it is negligible compared to the quantum corrections of the same form (see Eq. (34)). For this reason, in what follows we shall simply take $`a_4=0`$. Under this choice the action $`S_{HD}`$ satisfies the Noether identity

$`{\displaystyle \frac{2}{\sqrt{g}}}g_{\mu \nu }{\displaystyle \frac{\delta S_{HD}}{\delta g_{\mu \nu }}}=\mathrm{\hspace{0.17em}0},`$ (7)
which is usually interpreted as zero trace for the stress tensor of vacuum $`T_\mu ^\mu =0`$. At quantum level, this identity is violated by the anomaly such that

$`T=<T_\mu ^\mu >={\displaystyle \frac{2}{\sqrt{g}}}g_{\mu \nu }{\displaystyle \frac{\delta \overline{\mathrm{\Gamma }}^{(1)}}{\delta g_{\mu \nu }}}=(wC^2+bE+c^2R).`$ (8)
Here $`\overline{\mathrm{\Gamma }}^{(1)}`$ is the quantum correction to the classical vacuum action (1), the coefficients $`w`$, $`b`$ and $`c`$ are the $`\beta `$-functions for the parameters $`a_1,a_2,a_3`$ correspondingly

$`w={\displaystyle \frac{1}{(4\pi )^2}}\left({\displaystyle \frac{N_0}{120}}+{\displaystyle \frac{N_{1/2}}{20}}+{\displaystyle \frac{N_1}{10}}\right),`$ (9)
$`b={\displaystyle \frac{1}{(4\pi )^2}}\left({\displaystyle \frac{N_0}{360}}+{\displaystyle \frac{11N_{1/2}}{360}}+{\displaystyle \frac{31N_1}{180}}\right),`$ (10)
$`c={\displaystyle \frac{1}{(4\pi )^2}}\left({\displaystyle \frac{N_0}{180}}+{\displaystyle \frac{N_{1/2}}{30}}{\displaystyle \frac{N_1}{10}}\right).`$ (11)
Some observation is in order. The direct application of the dimensional regularization  leads to the ambiguous result for the coefficient of the $`^2R`$-term in the anomaly. We accept the expression for $`c`$ which emerges from other regularizations such as point-splitting  and $`\zeta `$-regularization .

The general solution for the anomaly-induced effective action is

$$\mathrm{\Gamma }_{ind}=S_c[\overline{g}_{\mu \nu }]+d^4x\sqrt{\overline{g}}\{w\sigma \overline{C}^2+b\sigma (\overline{E}\frac{2}{3}\overline{}^2\overline{R})+2b\sigma \overline{\mathrm{\Delta }}_4\sigma $$

$`{\displaystyle \frac{1}{12}}(c+{\displaystyle \frac{2}{3}}b)[\overline{R}6(\overline{}\sigma )^26(\overline{}^2\sigma )]^2)\},`$ (12)
where

$$g_{\mu \nu }=\overline{g}_{\mu \nu }e^{2\sigma },\sqrt{g}R^2=\sqrt{\overline{g}}[\overline{R}6(\overline{}\sigma )^26(\overline{}^2\sigma )]^2$$

and

$$\mathrm{\Delta }_4=^4+2R^{\mu \nu }_\mu _\nu \frac{2}{3}R^2+\frac{1}{3}(^\mu R)_\mu $$

is the fourth derivative, conformal operator acting on scalars. $`S_c[\overline{g}_{\mu \nu }]=S_c[g_{\mu \nu }]`$ is the “integration constant” for the equation (8). In general, there is no regular method for deriving $`S_c[\overline{g}_{\mu \nu }]`$. Fortunately, for some cosmological problems the form of this functional has no importance. Let us consider an isotropic and homogeneous metric

$`g_{\mu \nu }=\overline{g}_{\mu \nu }a^2(\eta ),`$ (13)
where $`\eta `$ is conformal time and

$`d\overline{s}^2=\overline{g}_{\mu \nu }dx^\mu dx^\nu =d\eta ^2{\displaystyle \frac{dr^2}{1kr^2}}r^2d\mathrm{\Omega }.`$ (14)
Then, the conformal functional $`S_c[\overline{g}_{\mu \nu }]`$ is constant, for it does not depend on $`a(\eta )`$ and does not contribute to the equations of motion. In this particular case, the solution (12) is an exact one-loop quantum correction. However, if one takes a more general background, for example

$`g_{\mu \nu }=\overline{g}_{\mu \nu }a^2(\eta )+h_{\mu \nu }(\eta ,𝐫),`$ (15)
where $`h_{\mu \nu }`$ is a metric perturbation. Then the functional $`S_c[\overline{g}_{\mu \nu }]`$ becomes relevant, because

$$S_c[g_{\mu \nu }]=S_c[\overline{g}_{\mu \nu }+a^2(\eta )h_{\mu \nu }(\eta ,𝐫)]$$

and the solution (12) is only an approximation to the unknown general expression. In this and the next two sections we shall consider the metrics of the form (14) and in the section 5 the metric (15).

The total action with the quantum corrections

$`S_t=S_{vacuum}+\mathrm{\Gamma }_{ind},`$ (16)
leads to the following equation for the metric (13):

$`{\displaystyle \frac{\stackrel{\mathrm{}.}{a}}{a}}+{\displaystyle \frac{3\stackrel{.}{a}\stackrel{\mathrm{}}{a}}{a^2}}+{\displaystyle \frac{\stackrel{..}{a}^2}{a^2}}\left(5+{\displaystyle \frac{4b}{c}}\right){\displaystyle \frac{\stackrel{..}{a}\stackrel{.}{a}^2}{a^3}}2k\left(1+{\displaystyle \frac{2b}{c}}\right){\displaystyle \frac{\stackrel{..}{a}}{a^3}}{\displaystyle \frac{M_P^2}{8\pi c}}\left({\displaystyle \frac{\stackrel{..}{a}}{a}}+{\displaystyle \frac{\stackrel{.}{a}^2}{a^2}}+{\displaystyle \frac{k}{a^2}}{\displaystyle \frac{2\mathrm{\Lambda }}{3}}\right)=\mathrm{\hspace{0.17em}0},`$ (17)
where we used the standard notation for the Planck mass $`M_P^2=1/G`$. The point over the variable designates the derivative with respect to the physical time $`t`$, the last is related to the conformal time $`\eta `$ by the relation $`dt=a(\eta )d\eta `$.

The equation (17) has remarkable particular solutions for all three cases $`k=0,\pm 1`$:

$`k=0a(t)=a_0\mathrm{exp}(Ht)`$ (18)
$`k=1a(t)={\displaystyle \frac{1}{H}}\mathrm{sinh}(Ht),`$ (19)
$`k=+1a(t)={\displaystyle \frac{1}{H}}\mathrm{cosh}(Ht),`$ (20)
where the Hubble constant $`H`$ has the form (we do not consider solutions with negative $`H`$)

$`H={\displaystyle \frac{M_P}{\sqrt{32\pi b}}}\left(\mathrm{\hspace{0.17em}1}\pm \sqrt{1+{\displaystyle \frac{64\pi b}{3}}{\displaystyle \frac{\mathrm{\Lambda }}{M_P^2}}}\right)^{1/2}.`$ (21)
For a negative $`\mathrm{\Lambda }`$, there is only a particular solution with positive sign in (21).

The solutions for the case of $`\mathrm{\Lambda }=0`$ have been found by Starobinsky  (see also ). According to (21) and (10), the two solutions with $`\pm `$ exist only for $`\mathrm{\Lambda }0`$. Let us remark that in the high energy region the sign of $`\mathrm{\Lambda }`$ is positive. Let us advocate this point. In the nowadays universe $`\mathrm{\Lambda }`$ is also positive but very small, due to the cancelation of the induced and vacuum energy densities (see, e.g.  and also  for the pedagogical introduction to the formal aspects of the cosmological constant problem).

The induced cosmological constant appears due to the spontaneous symmetry breaking in the Standard Model, as a vacuum energy corresponding to the minima of the Higgs potential

$`V_{cl}={\displaystyle \frac{1}{2}}m^2\varphi ^2+{\displaystyle \frac{f}{8}}\varphi ^4,`$ (22)
and has the classical value

$`\mathrm{\Lambda }_{ind}=<V_{cl}>={\displaystyle \frac{m^4}{2f}}<0`$ (23)
with the magnitude about 55 orders greater than the observable cosmological constant. Obviously, the vacuum counterpart must be positive in order the provide a cancelation at low energies. But, at high energies (temperatures) above $`300GeV`$ the symmetry in the potential (22) gets restored such that the induced part (23) disappears. After that the value of the cosmological constant will be given by the vacuum term with the corresponding quantum corrections . Similar increase of $`\mathrm{\Lambda }`$ may be expected in any other phase transition which could occur at high energies. In particular, the contribution to the cosmological constant may be expected from the supersymmetry breaking in case this breaking is spontaneous. The last observation is that after inflation, when the matter content of the Universe gets closer to the equilibrium state, the high-temperature radiation dominates over the cosmological term .

In what follows, except if this is indicated explicitly, we assume that the cosmological constant $`\mathrm{\Lambda }`$ in the high energy inflationary epoch satisfies the relation $`0<\mathrm{\Lambda }M_P^2`$. Then the expansion of (21) in the small parameter $`\mathrm{\Lambda }/M_P^2`$ gives the two values

$`H=\sqrt{{\displaystyle \frac{\mathrm{\Lambda }}{3}}}\mathrm{and}H=\sqrt{{\displaystyle \frac{M_P^2}{16\pi b}}{\displaystyle \frac{\mathrm{\Lambda }}{3}}},`$ (24)
where the first solution is exactly that we meet without the quantum corrections. The error in both cases has the order of magnitude $`\mathrm{\Lambda }/M_P<<H.`$ In the present-day universe this means that what we observe is really $`\mathrm{\Lambda }`$, without significant quantum corrections. In the inflationary epoch, if we compare the second solution (24) and the ones in the $`\mathrm{\Lambda }=0`$ case, the spectrum of $`H`$ is modified by the relatively small $`\mathrm{\Lambda }`$-term such that the two values of $`H`$ become closer to each other .

As we shall see in the next section, the inflationary solutions (18), (19) and (20) are stable for $`c>0`$ and unstable for $`c<0`$  (see also ). According to (11), $`c<0`$ is equivalent to the condition for a particle content of the theory

$`N_1>{\displaystyle \frac{1}{3}}N_{1/2}+{\displaystyle \frac{1}{18}}N_0.`$ (25)
The unique way to provide the non-stability of inflation in the case of $`c>0`$ is to introduce some additional term or terms into the classical action of vacuum. For instance, the $`\sqrt{g}R^2`$-term with the sufficiently large positive coefficient does this job. But, for the first sight, this term is not necessary. Let us notice, that the quantum correction (12) is valid not only in the Early Universe but also at the later stages of the evolution. However, the content $`N_0,N_{1/2},N_1`$ may be different, because for the small energies the loops of the fields with large masses decouple. Let us consider the present-day universe as an example. The value of the Hubble constant today is extremely small, about $`H10^{42}GeV`$, that is about 30 orders smaller than the mass of the lightest neutrino. Therefore, the unique particle which makes a nontrivial contribution to the anomaly is the photon 555If the neutrino would be massless, there could be an equality in (25).. Then, the relation (25) testifies a non-stability and the present-day universe is safe without the special $`\sqrt{g}R^2`$-term. We shall continue the discussion of this example in the section 3.

The condition (25) enables one to construct a very attractive inflationary scenario . The universe could start in the stable phase, such that inflation starts independent on the initial data. The simplest way to provide stability in (25) is to assume supersymmetry in the high energy region and also suppose that the supersymmetry is absent at a lower energy because the sparticles are heavy . If, in the course of inflation, the magnitude of the Hubble constant decreases, then at some point the loops of the sparticles decouple and the matter content $`N_{0,1/2,1}`$ gets modified. As a result, the sign of the inequality (25) changes to the opposite and the universe falls into the non-stable inflation with the eventual transition to the FRW evolution. The important problem is that why the Hubble parameter should decrease during the inflation? Indeed, the Starobinsky solutions (18), (19) and (20) are characterized by the constant $`H`$. The answer to this question is that the value of $`H`$ really decreases if we go beyond the massless approximation . We shall consider the effect of massive fields below.

Another interesting aspect is the possible role of matter. Here we have two different situations. The first is the possibility to have a set of heavy massive fields at the beginning of inflation. These fields can be remnants of the string phase transition. For our model it is important that this massive content does not destroy the stability which holds at the beginning of the stable inflation. Let us suppose that these fields create zero pressure. Then the equation (17) is modified by adding the term $`\rho _ma^3`$ to the r.h.s.. The investigation of the equation with this term and with $`\mathrm{\Lambda }=k=0`$

$`{\displaystyle \frac{\stackrel{\mathrm{}.}{a}}{a}}+{\displaystyle \frac{3\stackrel{.}{a}\stackrel{\mathrm{}}{a}}{a^2}}+{\displaystyle \frac{\stackrel{..}{a}^2}{a^2}}\left(5+{\displaystyle \frac{4b}{c}}\right){\displaystyle \frac{\stackrel{..}{a}\stackrel{.}{a}^2}{a^3}}2k\left(1+{\displaystyle \frac{2b}{c}}\right){\displaystyle \frac{\stackrel{..}{a}}{a^3}}{\displaystyle \frac{M_P^2}{8\pi c}}\left({\displaystyle \frac{\stackrel{..}{a}}{a}}+{\displaystyle \frac{\stackrel{.}{a}^2}{a^2}}\right)={\displaystyle \frac{\rho _m}{ca^3}},`$ (26)
has been performed numerically. The result is such that even for $`\rho _m=M_P^4`$ the system does not loose stability and the stabilization performs in a few Planck times. The plots are almost identical to the matterless case (see Fig. 2). The reason of this output is that the new term decreases very fast during the inflation and becomes irrelevant.

Second aspect is the transition to the FRW evolution after some of the massive fields decouple and the inflation becomes unstable. In order to understand the situation better, let us just substitute the corresponding FRW solution $`a(t)t^{2/3}`$ into the equation (26). It is easy to see that the “classical” terms - both Einstein and matter, both behave like $`1/t^2`$, while the quantum correction - with higher derivatives, behave like $`1/t^4`$. Indeed, $`H(t)0`$, as usual. Therefore, with the properly taken initial data for the unstable inflation, the FRW solution is a perfect approximation for the later epochs. In  we have checked that the “proper” initial data emerge naturally in the framework of the sharp cut-off approach to the decoupling.

### 2.2 Massive case

Now, let us assume that some of the scalar and fermion fields are massive. At high energies, the leading vacuum quantum effect of the masses is the renormalization of the Newton and cosmological terms. In  the renormalization of the Newton term has been taken into account using some special approach which gives direct link to the cosmological applications. Here we shall follow the same method, but perform the calculations in a more detailed and general form. In particular, we shall take into account the renormalization of the cosmological term.

The conformal invariance of the actions (4) and (5) is violated by the masses and we can not use the conformal anomaly to derive quantum corrections. But, this can be changed if we apply the conformal description of the massive theory in the framework of the cosmon model  (see also  for other applications of the cosmon model). Let us replace the massive parameters by the powers of a new auxiliary scalar field $`\chi `$ according to

$$m_s^2\frac{m_s^2}{M^2}\chi ^2,m_f\frac{m_f}{M}\chi ,$$

$`{\displaystyle \frac{1}{16\pi G}}R{\displaystyle \frac{M_P^2}{16\pi M^2}}\left[R\chi ^2+6(\chi )^2\right],\mathrm{\Lambda }{\displaystyle \frac{\mathrm{\Lambda }}{M^2}}\chi ^2,`$ (27)
where $`M`$ is some dimensional parameter. In the original version of the cosmon model , the introduction of the new scalar fields $`\chi `$ is called to provide the dilatation invariance with the consequent breaking of the dilatation symmetry. In our case we need a more general, local conformal invariance. Therefore, we postulate that the auxiliary field $`\chi `$ transforms as

$`\chi \chi e^\sigma ,\sigma =\sigma (x),`$ (28)
independent of the space-time dimension 666Other fields transform in a usual way: $`g_{\mu \nu }g_{\mu \nu }\mathrm{exp}(2\sigma ),\phi \phi \mathrm{exp}(\frac{n2}{2}\sigma ),\psi \psi \mathrm{exp}(\frac{n1}{2}\sigma ).`$ Here we use the dimensional regularization, because it helps to simplify the calculations.. After we apply the procedure (27), the massive terms in (4) and (5) are replaced by Yukawa and quartic scalar interactions between physical fermion and scalar fields and the new auxiliary scalar $`\chi `$. The advantage of this procedure is that the new actions in both matter and gravitational sectors become conformal invariant. This invariance is violated by a conformal anomaly and one can derive the effective action of the background fields $`g_{\mu \nu },\chi `$ using the anomaly-induced effective action scheme.

The divergences of the theory in the conformal representation have the form

$`\mathrm{\Gamma }_{div}^{(1)}={\displaystyle \frac{\mu ^{n4}}{(n4)}}{\displaystyle d^nx\sqrt{g}\left\{wC^2+bE+c^2R+\frac{f}{M^2}[R\chi ^2+6(\chi )^2]+\frac{g}{M^4}\chi ^4\right\}},`$ (29)
where the coefficients $`w,b`$ and $`c`$ are given by Eqs. (9), (10) and (11) correspondingly, and we introduced useful notations for the dimensional quantities

$`f={\displaystyle \frac{1}{3(4\pi )^2}}{\displaystyle \underset{f}{}}N_fm_f^2,`$ (30)
and

$`g={\displaystyle \frac{1}{2(4\pi )^2}}{\displaystyle \underset{s}{}}N_sm_s^4{\displaystyle \frac{2}{(4\pi )^2}}{\displaystyle \underset{f}{}}N_fm_f^4.`$ (31)
Here the sums are taken over all massive fermion $`f`$ and scalar $`s`$ fields with the masses $`m_f`$ and $`m_s`$ correspondingly. $`N_f`$ and $`N_s`$ are multiplicities of fermions and scalars, for example $`N_{lepton}=1`$ and $`N_{quark}=3`$. In order to derive the conformal anomaly, one needs to consider the classical Noether identity for the vacuum part of the effective action

$`𝒯={\displaystyle \frac{2}{\sqrt{g}}}g_{\mu \nu }{\displaystyle \frac{\delta S_{vac}}{\delta g_{\mu \nu }}}+{\displaystyle \frac{1}{\sqrt{g}}}\chi {\displaystyle \frac{\delta S_{vac}}{\delta \chi }}=\mathrm{\hspace{0.17em}0}.`$ (32)
It is easy to see that, in the theory under discussion, the conformal invariance does not mean that the stress tensor is traceless, but instead that the quantity $`𝒯`$ is zero. Correspondingly, the conformal anomaly means $`<𝒯>\mathrm{\hspace{0.17em}0}`$ instead of usual $`<T_\mu ^\mu >\mathrm{\hspace{0.17em}0}`$ . However, since both $`\chi `$-dependent terms in (29) have the same conformal properties as the square of the Weyl tensor, they can be absorbed into the Weyl term such that the problem of deriving anomaly reduces to the standard one. Finally, we arrive at the expression

$`<𝒯>=\left\{wC^2+bE+c^2R+{\displaystyle \frac{f}{M^2}}[R\chi ^2+6(\chi )^2]+{\displaystyle \frac{g}{M^4}}\chi ^4\right\}.`$ (33)
In the standard way, using the parametrization $`g_{\mu \nu }=e^{2\sigma }\overline{g}_{\mu \nu }`$ and $`\chi =e^\sigma \overline{\chi }`$, one can derive the anomaly-induced effective action of the background fields $`g_{\mu \nu }`$ and $`\chi `$

$$\mathrm{\Gamma }_{ind}=S_c[g_{\mu \nu },\chi ]+d^4x\sqrt{\overline{g}}\{w\overline{C}^2\sigma +b(\overline{E}\frac{2}{3}\overline{}^2\overline{R})\sigma +2b\sigma \overline{\mathrm{\Delta }}\sigma +$$

$`+{\displaystyle \frac{f}{M^2}}[\overline{R}\overline{\chi }^2+6(\overline{\chi })^2]\sigma +{\displaystyle \frac{g}{M^4}}\overline{\chi }^4\sigma \}{\displaystyle \frac{3c+2b}{36}}{\displaystyle }d^4x\sqrt{g}R^2.`$ (34)
The last step is to fix the conformal unitary gauge $`\chi =\overline{\chi }e^\sigma =M`$, such that the classical Einstein-Hilbert and cosmological terms acquire their standard form. The one-loop effective action becomes

$$\mathrm{\Gamma }=S_{vacuum}+\mathrm{\Gamma }_{ind}[\chi M]=S_{HD}+S_c[g_{\mu \nu },M]+$$

$$+d^4x\sqrt{\overline{g}}\{w\overline{C}^2\sigma +b(\overline{E}\frac{2}{3}\overline{}^2\overline{R})\sigma +\mathrm{\hspace{0.17em}2}b\sigma \overline{\mathrm{\Delta }}\sigma \}\frac{3c+2b}{36}d^4x\sqrt{g}R^2$$

$`{\displaystyle d^4x\sqrt{\overline{g}}e^{2\sigma }[\overline{R}+6(\overline{}\sigma )^2]\left[\frac{1}{16\pi G}f\sigma \right]}{\displaystyle d^4x\sqrt{\overline{g}}e^{4\sigma }\left[\frac{\mathrm{\Lambda }}{8\pi G}g\sigma \right]}.`$ (35)
It is worth noticing that the expression (35) does not depend on the magnitude of $`M`$, but only on the fields and parameters of the theory. This effective action differs from the massless counterpart (16) in two respects. First, there are new terms which emerge due to the renormalization of the Einstein-Hilbert and cosmological term. Second, the integration constant $`S_c[g_{\mu \nu },M]`$ is not irrelevant anymore, since we imposed the conformal unitary gauge. Therefore, the expression (35) can be regarded only as an approximation to the effective action of vacuum for the massive fields and we have to learn the limits of validity for this approximation.

The higher-derivative part of the Eq. (35) is identical to that for the massless fields, as it has to be (see, e.g. ). Furthermore, there is a strong link between the anomaly-induced effective action Eq.(35) and the quantum corrections coming from the renormalization group. The expansion of the homogeneous, isotropic universe means a conformal transformation of the metric $`g_{\mu \nu }(t)=\overline{g}_{\mu \nu }\mathrm{exp}[\sigma (\eta )]`$. On the other hand, the renormalization group in curved space-time corresponds to the scale transformation of the metric $`g_{\mu \nu }g_{\mu \nu }e^{2\tau }`$ simultaneously with the inverse transformation of all dimensional quantities . For any $`\mu `$ we have $`\mu \mu e^\tau `$. One can compare the $`\sigma `$-dependence of the anomaly-induced effective action (35) and the $`\tau `$-dependence of the renormalization-group improved classical action

$`S_{vac.imp}=S_{vacuum}[P(\tau )],`$ (36)
where $`P=\{a_{1,2,3,4},G,\mathrm{\Lambda }\}`$ denote vacuum parameters of the theory and $`P(\tau )=P_0+\beta _P\tau `$. The expression (36) is a leading-log approximation for the solution of the renormalization group equation for the effective action

$`\mathrm{\Gamma }[e^{2\tau }g_{\alpha \beta },\mathrm{\Phi }_i,P,\mu ]=\mathrm{\Gamma }[g_{\alpha \beta },\mathrm{\Phi }_i(\tau ),P(\tau ),\mu ],`$ (37)
where $`\mathrm{\Phi }_i`$ are matter fields. It is easy to see that (35) becomes completely equivalent to (36) if we set $`\sigma =const`$. The coefficient $`f`$ is a factor of the $`\beta `$-function for the inverse Newton constant $`\mathrm{\hspace{0.17em}1}/16\pi G`$ and the coefficient $`g`$ is a factor of the $`\beta `$-function for the cosmological term $`\mathrm{\Lambda }/8\pi G`$. Indeed, (37) can be considered as a generalization of the renormalization group improved classical action (36).

Therefore, we have strong reasons to regard Eq. (35) as a leading-log approximation for the effective action for the massive fields. This approximation picks up the logarithmic quantum corrections and is reliable in the high energy region where masses of the fields are much smaller than the Hubble constant $`H`$. On the other hand, at low energies the massive fields decouple and their quantum effects become negligible. In this paper we shall interpolate between these two regimes using the MS scheme and the “sharp cut off” approximation at some specific scale $`M_{}`$. This scale is defined such that at $`H=M_{}`$ a sufficient number of the sparticles decouple and the inequality (25) changes its sign. Within our approximation, the functional $`S_c`$ contains the sub-leading non-log corrections and can be disregarded. Let us remark that the $`M_{}`$ scale can be different from the supersymmetry breaking scale $`M_{SUSY}`$. In particular, we can suppose that for a GUT model $`M_{SUSY}`$ is about $`10^{16}GeV`$ while the upper bound for $`M_{}`$ is $`10^{14}GeV`$. This bound appears because the amplitude of the generated gravitational waves has an admissible range only for $`H/M_P10^5`$ in the last stage ($`65e`$-folds) of inflation.

The equation of motion for $`\sigma (t)`$ has the form

$$\stackrel{\mathrm{}.}{\sigma }+7\stackrel{\mathrm{}}{\sigma }\stackrel{.}{\sigma }+4\stackrel{..}{\sigma }^2+4(3\frac{b}{c})\stackrel{..}{\sigma }\stackrel{.}{\sigma }^24\frac{b}{c}\stackrel{.}{\sigma }^42k(1+\frac{2b}{c})(\stackrel{.}{\sigma }^2+\stackrel{..}{\sigma })e^{2\sigma }$$

$`{\displaystyle \frac{M_P^2}{8\pi c}}[(\stackrel{..}{\sigma }+2\stackrel{.}{\sigma }^2+ke^{2\sigma })(1\stackrel{~}{f}\sigma ){\displaystyle \frac{1}{2}}\stackrel{~}{f}\dot{\sigma }^2]+{\displaystyle \frac{M_P^2\mathrm{\Lambda }}{12\pi c}}(1\stackrel{~}{g}\sigma \stackrel{~}{g}/4)=\mathrm{\hspace{0.17em}0},`$ (38)
where we introduced useful notations

$`\stackrel{~}{f}={\displaystyle \frac{16\pi f}{M_P^2}}={\displaystyle \frac{1}{3\pi }}{\displaystyle \underset{f}{}}{\displaystyle \frac{N_fm_f^2}{M_P^2}};\stackrel{~}{g}={\displaystyle \frac{8\pi g}{M_P^2\mathrm{\Lambda }}}={\displaystyle \frac{1}{4\pi }}{\displaystyle \underset{s}{}}{\displaystyle \frac{N_sm_s^4}{M_P^2\mathrm{\Lambda }}}{\displaystyle \frac{1}{\pi }}{\displaystyle \underset{f}{}}{\displaystyle \frac{N_fm_f^4}{M_P^2\mathrm{\Lambda }}}.`$ (39)
It proves useful to rewrite this equation, also, in terms of $`a(t)`$777We correct a small mistake in the similar expression (without the cosmological constant term) in . This correction does not modify the behavior of $`a(t)`$.

$$\frac{\stackrel{\mathrm{}.}{a}}{a}+3\frac{\stackrel{.}{a}}{a}\frac{\stackrel{\mathrm{}}{a}}{a}+\frac{\stackrel{..}{a}^2}{a^2}\left(5+\frac{4b}{c}\right)\frac{\stackrel{..}{a}}{a}\frac{\stackrel{.}{a}^2}{a^2}2k\left(1+\frac{2b}{c}\right)\frac{\stackrel{..}{a}}{a^3}$$

$`{\displaystyle \frac{M_P^2}{8\pi c}}\left[\left({\displaystyle \frac{\stackrel{..}{a}}{a}}+{\displaystyle \frac{\dot{a}^2}{a^2}}+{\displaystyle \frac{k}{a^2}}\right)(1\stackrel{~}{f}\mathrm{ln}a){\displaystyle \frac{\stackrel{~}{f}}{2}}{\displaystyle \frac{\dot{a}^2}{a^2}}\right]+{\displaystyle \frac{M_P^2\mathrm{\Lambda }}{12\pi c}}(\mathrm{\hspace{0.17em}1}\stackrel{~}{g}\mathrm{ln}a\stackrel{~}{g}/4)=\mathrm{\hspace{0.17em}0}.`$ (40)
The solution of the equation (38) can not be performed analytically in the general form, but it is obvious that the approximate solution for small $`\mathrm{ln}a(t)`$ can be obtained by the replacement

$`M_P^2\stackrel{~}{M}_P^2=M_P^2\left[1\stackrel{~}{f}\mathrm{ln}a(t)\right];\mathrm{\Lambda }\stackrel{~}{\mathrm{\Lambda }}=\mathrm{\Lambda }\left\{1\stackrel{~}{g}\mathrm{ln}a(t)\right\}`$ (41)
in the expression for the Hubble parameter (21) corresponding to the solutions (18), (19) and (20). Another possibility is the numerical integration of the equation (38). We present the plots produced by the numerical integration at Figure 1 for the case $`\mathrm{\Lambda }=0`$. These plots were obtained using Mathematica program  and also FORTRAN and standard algorithms of . It is easy to see that the initial parts of these numerical solutions are identical to the ones which were obtained in . Let us make one more step in the understanding of these plots and the quality of the approximation (41). For the sake of simplicity, we shall consider only the value $`\mathrm{\Lambda }=0`$ 888The numerical analysis shows that this is a good approximation for $`\mathrm{\Lambda }\stackrel{~}{f}M_P^2`$.. In this case one can use (41) and easily integrate the equation

$`\dot{\sigma }=H=H_0\sqrt{1\stackrel{~}{f}\sigma },H_0={\displaystyle \frac{M_P}{\sqrt{16b}}},\sigma (0)=0.`$ (42)
and find its solution in the form

$`\sigma (t)=H_0t{\displaystyle \frac{H_0^2}{4}}\stackrel{~}{f}t^2.`$ (43)
One might expect that this simple formula should serve as a reasonable approximation only at the initial stage when $`\sigma (t)=\mathrm{ln}a(t)`$ is small. However, if we compare the plots (a) and (b) at the Figure 1 and the parabola Fig. 1c (43), the numerical difference is remarkably small - about $`0.1`$ at the maximal point with $`\sigma _{max}=10^4`$ (for the toy model with the MSSM field content and $`\stackrel{~}{f}=10^5`$). The two solutions are almost identical for all $`t`$. Remarkably, this is not a specific feature of this particular value of $`\stackrel{~}{f}`$. If we substitute (43) into the equation (38), the result is

$$1.210^8+(25.9\stackrel{~}{f}2.5910^8)t+(3.1410^{17}3.1410^8f)t^2+$$

$$+(1.2210^{26}+1.2210^{17}f)t^3+(1.4810^{36}1.4810^{27}f)t^4$$

that is indeed negligible for all situations of interest, when $`\stackrel{~}{f}`$ is very small. This unexpected feature of the solution is extremely useful, for we can sometimes use simple expression (43) instead of the complicated numerical solution of the equation (38).

In the last considerations we took into account that $`\stackrel{~}{f}`$ and $`\mathrm{\Lambda }`$ are very small quantities. This leads to the condition for the magnitude of the quantum corrections to the cosmological constant

$`\left|{\displaystyle \frac{1}{2}}{\displaystyle \underset{s}{}}m_s^42{\displaystyle \underset{f}{}}m_f^4\right|\mathrm{ln}a(t)8\pi \mathrm{\Lambda }M_P^2.`$ (44)
Furthermore, one can request that the cosmological constant plays smaller role than the higher derivative induced contributions:

$`\left|{\displaystyle \frac{1}{2}}{\displaystyle \underset{s}{}}m_s^42{\displaystyle \underset{f}{}}m_f^4\right|M_P^2\left[1\stackrel{~}{f}\mathrm{ln}a(t)\right]^2\mathrm{and}\mathrm{\Lambda }8\pi M_P^2\left[1\stackrel{~}{f}\mathrm{ln}a(t)\right]^2.`$ (45)
Indeed, these conditions are satisfied at the beginning of inflation when $`\mathrm{ln}a(t)`$ is small and the higher derivative terms dominate over the cosmological constant which must be, at most, of the order $`M_{SUSY}^4`$. It is convenient to assume the soft supersymmetry breaking scheme such that the magnitude of the cosmological constant is much smaller than $`M_{}^4/M_P^2`$. When $`\mathrm{ln}a(t)`$ becomes large and the universe is close to the transition from stable to unstable inflation, we can not rely much on Eq.(41) and have to use the numerical methods to verify both stability of the tempered inflationary solution and the graceful exit to the FRW phase after the supersymmetry breaks down.

## 3 The conditions of stability in the massless case

In this section we shall consider the stability of the relatively simple massless case (17). Instead of the Eq. (17) one can take the third order $`00`$-component of the equation

$`R_{\mu \nu }{\displaystyle \frac{1}{2}}(R2\mathrm{\Lambda })g_{\mu \nu }=<T_{\mu \nu }>.`$ (46)
Without the matter fields the $`00`$-equation is equivalent to the Eq. (17) which we used in the previous section . The $`00`$-equation can be solved analytically  and the phase diagram shows the existence of the single attractor corresponding to the inflationary solution in the stable case and various attractors, including $`H=0`$ in the non-stable case with $`k=\mathrm{\Lambda }=0`$. Unfortunately, our equations (38) or (40) are much more complicated than in the massless case and there it is not possible to solve them exactly. However, there is no real need to do that. At the beginning of inflation the stabilization of the exponential solution performs very fast . Therefore, at the latest stages it is sufficient to check the asymptotic stability of the inflationary solution with respect to the small perturbations in the sense of Lyapunov (see, e.g.  regarding the stability issues in higher derivative models).

In this section we establish some important general properties of the stability under small perturbations and will apply this knowledge for the approximate analytical and numerical analysis of the massive case in the next section.

With respect to the criterion (25) of stability , there are three relevant questions: i) which kind of perturbations one has to consider; ii) whether the choice of $`k=0,1,1`$ may affect (25); iii) whether the presence of the cosmological constant may affect this criterion.

In order to address the issue i), let us make small perturbations in the equation (17) and set, for the sake of simplicity, $`\mathrm{\Lambda }=0`$ and $`k=0`$. First we try the perturbations of the form

$$a(t)=a_0(t)+x(t),\mathrm{where}a_0(t)=\mathrm{exp}(Ht)\mathrm{and}H=\frac{M_P}{\sqrt{16\pi b}}.$$

After simple calculations we arrive at the equation for the perturbations $`x(t)`$

$`\stackrel{\mathrm{}.}{x}+3H\stackrel{\mathrm{}}{x}\left(3+{\displaystyle \frac{2b}{c}}\right)H^2\stackrel{..}{x}\left(7+{\displaystyle \frac{4b}{c}}\right)H^3\stackrel{.}{x}+6\left(1+{\displaystyle \frac{b}{c}}\right)H^4x=0.`$ (47)
Indeed, the factors of $`H`$ can be easily removed by the simple rescaling of time, and then the standard Routh-Hurwitz conditions (see the Appendix) show that the stability can not be achieved for any values of $`b/c`$. Of course, this contradicts our previous statement (25), and we can conclude that the perturbations of $`a(t)`$ do not serve for our purposes.

Fortunately, from the theory of ordinary differential equations (see, e.g. ) we know that the stability may depend on the choice of the variables. Therefore, let us consider the same problem in terms of $`\sigma (t)`$. The corresponding equation can be obtained from (38) by setting $`\mathrm{\Lambda }=k=\stackrel{~}{f}=\stackrel{~}{g}=0`$ and making small perturbation of the form

$$\sigma (t)=\sigma _0(t)+y(t),\mathrm{where}\sigma _0(t)=Ht.$$

The equation for $`y(t)`$ has the form

$`\stackrel{\mathrm{}.}{y}+7H\stackrel{\mathrm{}}{y}+[4(3{\displaystyle \frac{b}{c}})H^2{\displaystyle \frac{M_P^2}{8\pi c}}]\stackrel{..}{y}[{\displaystyle \frac{16b}{c}}H^3+{\displaystyle \frac{M_P^2H}{2\pi c}}]\stackrel{.}{y}=\mathrm{\hspace{0.17em}0}.`$ (48)
or, taking the solution for $`H`$ into account,

$`\stackrel{\mathrm{}.}{y}+7H\stackrel{\mathrm{}}{y}+2(6{\displaystyle \frac{b}{c}})H^2\stackrel{..}{y}{\displaystyle \frac{8b}{c}}H^3\stackrel{.}{y}=\mathrm{\hspace{0.17em}0}.`$ (49)
It is easy to see that in this case the stability is achieved when $`c`$ is positive that immediately gives (25) 999Of course, the conditions of stability for $`H(t)`$ are exactly those for $`\sigma (t)`$.. The illusory contradiction between the two results can be clarified if we identify

$$y(t)=\mathrm{ln}\left[\mathrm{\hspace{0.17em}1}+\frac{x(t)}{a_0(t)}\right].$$

Since $`a_0(t)`$ grows up very fast, the growing modes in $`x(t)`$ may be negligible because they are suppressed by the factor of $`1/a_0(t)`$. Of course, this is the case when $`c>0`$, while for $`c<0`$ there is a mode in $`x(t)`$ which is growing up faster than $`a_0(t)`$. Finally, we learned a very important lesson for the consequent study of the complicated massive case: one has to perturb $`\sigma (t)`$, but not $`a(t)`$. The calculations for the $`k0`$ and $`\mathrm{\Lambda }0`$ cases are a bit more complicated, and we will not bother the reader with the details. The condition of stability is always the same as in the simplest situation described above. But, it is important that this condition must be the same for $`k=0,1,1`$ and must be independent on $`\mathrm{\Lambda }`$. Let us formulate two simple theorems:

Theorem 1. The value of $`\mathrm{\Lambda }`$ does not influence the asymptotic stability.

Proof. The equation for the perturbation $`x(t)`$ or $`y(t)`$ does not change if the original equation for $`a(t)`$ or $`\sigma (t)`$ is multiplied by the factor of $`a^n`$ or $`\mathrm{exp}(n\sigma )`$. Therefore, we can always isolate the cosmological term such that its contribution to the perturbation equation is zero. The dependence on $`\mathrm{\Lambda }`$ can be also seen in $`H`$, but the magnitude of $`H`$ is irrelevant for stability, because we can always renormalize time such that in the new equation $`H1`$.

Observation. This theorem is valid only for $`\mathrm{\Lambda }M_P^2`$.

Theorem 2. The choice of $`k=0,1,1`$ does not influence the asymptotic stability and the condition (25) is valid for all three cases.

Proof. It is sufficient to apply the corresponding theorem about the asymptotic stability in the limit $`t\mathrm{}`$ (see Appendix). In the equation for the perturbations one can always replace the coefficient functions by their asymptotic values. Since both $`sinh(Ht)`$ and $`cosh(Ht)`$ behave like an exponential at the late times, the conditions of stability for all three cases are the same.

Illustration. The equation for the perturbation in the general massless case is

$`\stackrel{\mathrm{}.}{y}[t]+7\underset{0}{\overset{.}{\sigma }}\stackrel{\mathrm{}}{y}[t]+[(1+{\displaystyle \frac{2b}{c}})(\underset{0}{\overset{2}{\stackrel{.}{\sigma }}}+\underset{0}{\overset{..}{\sigma }})+{\displaystyle \frac{M_P^2}{16\pi c}}]4ke^{2\sigma _0}h[t]`$
$`+\left[8\underset{0}{\overset{..}{\sigma }}+4\left(3{\displaystyle \frac{b}{c}}\right)\underset{0}{\overset{2}{\stackrel{.}{\sigma }}}2k\left(1+{\displaystyle \frac{2b}{c}}\right)e^{2\sigma _0}{\displaystyle \frac{M_P^2}{8\pi c}}\right]\stackrel{..}{y}[t]`$
$`+[7\underset{0}{\overset{\mathrm{}}{\sigma }}+8(3{\displaystyle \frac{b}{c}})\underset{0}{\overset{..}{\sigma }}\underset{0}{\overset{.}{\sigma }}{\displaystyle \frac{16b}{c}}\underset{0}{\overset{3}{\stackrel{.}{\sigma }}}4k(1+{\displaystyle \frac{2b}{c}})e^{2\sigma _0}\underset{0}{\overset{.}{\sigma }}{\displaystyle \frac{M_P^2}{2\pi c}}\underset{0}{\overset{.}{\sigma }}]\stackrel{.}{y}[t]=0.`$
If we neglect those terms in the coefficients of the perturbations which tend to zero at $`t\mathrm{}`$, the $`k`$ \- dependence disappears.

Thus, we have learned two lessons from the analysis of the massless case. It is sufficient to consider the case $`\mathrm{\Lambda }=0`$ and $`k=0`$, for this does not modify the stability of the solution. These statements can be used in the more realistic theory of inflation based on the effective action for the massive fields, but with a proper caution. In fact, the negligible role of $`\mathrm{\Lambda }`$ depends on that it is a constant. As we shall see in the next section, this is not true for the effective action of massive fields.

Until now, we have considered only the stable inflationary solution with $`c>0`$. In other words, we were interested in the high energy region where the supersymmetry is unbroken. Obviously, another end of the energy scale also represents a great interest, for it enables one to perform a simple and efficient test of the model. Let us consider, again, a nowaday universe with $`H=H_010^{42}GeV`$. According to the recent data , this magnitude of the Hubble constant is mainly due to the contribution of the cosmological constant. As we have already mentioned in the last section, the value of $`c`$ in the present-day universe is negative due to the contribution of the single massless particle - photon. Therefore, we have to check whether the higher derivative terms do not destroy the stability of the first solution in Eq. (24). Using (48), we arrive at the following characteristic equation for the perturbation of $`HH+const\mathrm{exp}(\lambda t)`$.

$`\lambda ^3+7H\lambda ^2+\left[4\left(3{\displaystyle \frac{b}{c}}\right)H^2{\displaystyle \frac{M_P^2}{8\pi c}}\right]\lambda \left[{\displaystyle \frac{16b}{c}}H^3+{\displaystyle \frac{M_P^2}{2\pi c}}H\right]=\mathrm{\hspace{0.17em}0}.`$ (50)
Since the explicit solution of this equation is tedious, let us start from the $`\mathrm{\Lambda }=0`$ case. Then $`H=0`$ and the roots of the equation (50) are

$`\lambda _1^{(0)}=0,\lambda _{2/3}^{(0)}=\pm {\displaystyle \frac{M_P}{\sqrt{8\pi |c|}}}i.`$ (51)
At this level, we do not have definite answer to our question about stability. Now, let us look for the solution of Eq. (50) making perturbations in (51). Due to the huge difference between $`H`$ and $`M_P`$, this approximation is perfect. The result is

$`\lambda _1=4H,\lambda _{2/3}={\displaystyle \frac{3}{2}}H\pm {\displaystyle \frac{M_P}{\sqrt{8\pi |c|}}}i,\mathrm{where}H=\sqrt{{\displaystyle \frac{\mathrm{\Lambda }}{3}}}>0.`$ (52)
It is very nice to see that the cosmological constant $`\mathrm{\Lambda }>0`$ really stabilizes the solution in the low-energy region, exactly as we should optimistically expect. The stability can be verified also for the earlier post-inflationary epoch, when the energy density of vacuum played smaller role than the density of radiation . For the later epoch, as we have already noticed, when we substitute the FRW solution into Eq. (17), the “quantum” terms decrease as $`1/a^2`$ compared to the Einstein-Hilbert and matter terms. Therefore, for the later epoch of the expanding post-inflationary universe the FRW is a very good approximation.

## 4 The stability in the interpolation regime

At the very beginning of the inflationary period, the evolution goes exactly as in the massless case and the condition of stability must be the same (25). Later on, when $`\sigma (t)`$ becomes large enough, the deviation from the exponential inflation becomes greater. Finally, the inflation slows down such that the Hubble parameter achieves the scale $`M_{}`$, and the massive sparticles decouple. Then the condition (25) breaks down. Our purpose is to see whether the Eq. (25) remains the condition of stability for the modified equations (38) until the point $`HM_{}`$.

Before starting the numerical analysis, let us perform some analytic investigation of the stability in a useful approximation framework. It is sufficient to investigate the behavior of the perturbations in the physically important region corresponding to the last 65 $`e`$-folds before $`H`$ achieves the value $`M_{}`$. Looking at the plot at Fig. 1, we see that the dependence $`\sigma (t)`$ is rather smooth. Hence, we can divide this region into small intervals $`(t_i,t_i+\mathrm{\Delta }t_i)`$ and consider $`H(t)=H_i=const`$ inside each of the intervals. Indeed, $`H`$ will change from one interval to another, and this may be the source of new perturbations. Of course, all this consideration is done for the stable case $`c>0`$, while for the non-stable case this approach has no sense. For the purpose of analytic study, we accept the approximation (41). Then, the direct calculations give the following equation for the perturbations $`\sigma \sigma +y(t)`$:

$`b_0\stackrel{\mathrm{}.}{y}+b_1\stackrel{\mathrm{}}{y}+b_2\stackrel{..}{y}+b_3\stackrel{.}{y}+b_4y=\mathrm{\hspace{0.17em}0},`$ (53)
where

$$b_0=1,b_1=7\stackrel{~}{H},b_2=4\left(3\frac{b}{c}\right)\stackrel{~}{H}^2\frac{\stackrel{~}{M}_P^2}{8\pi c},$$

$`b_3={\displaystyle \frac{16b}{c}}\stackrel{~}{H}^3{\displaystyle \frac{\stackrel{~}{M}_P^2\stackrel{~}{H}}{2\pi c}},b_4={\displaystyle \frac{\stackrel{~}{M}_P^2}{4\pi c}}\left(\stackrel{~}{H}^2\stackrel{~}{f}{\displaystyle \frac{1}{3}}\stackrel{~}{\mathrm{\Lambda }}\stackrel{~}{g}\right)`$ (54)
and

$`\stackrel{~}{H}^2={\displaystyle \frac{\stackrel{~}{M}_P}{32\pi b}}\left[\mathrm{\hspace{0.17em}1}+\left(1+{\displaystyle \frac{64\pi b}{3}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Lambda }}}{\stackrel{~}{M}_P^2}}\right)^{1/2}\right].`$ (55)
If we treat the slowly-varying $`\sigma (t)`$ as a constant, then the stability depends on the signs of the Routh-Hurwitz determinants (see Appendix). Obviously, $`D_1=b_1`$ is positive. After some algebra we arrive at the expression

$`D_2={\displaystyle \frac{21\stackrel{~}{M}_P^2\stackrel{~}{H}}{8\pi b}}\left[\mathrm{\hspace{0.17em}1}+\left(1{\displaystyle \frac{b}{7c}}\right)\left(1+{\displaystyle \frac{64\pi b}{3}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Lambda }}}{\stackrel{~}{M}_P^2}}\right)^{1/2}\right]>0.`$ (56)
The general expression for $`D_3`$ is

$$D_3=\frac{21\stackrel{~}{M}_P^4\stackrel{~}{H}^2}{16\pi ^2bc}\left[\mathrm{\hspace{0.17em}1}+\frac{64\pi b}{3}\frac{\stackrel{~}{\mathrm{\Lambda }}}{\stackrel{~}{M}_P^2}+\left(1\frac{b}{7c}\right)\left(1+\frac{64\pi b}{3}\frac{\stackrel{~}{\mathrm{\Lambda }}}{\stackrel{~}{M}_P^2}\right)^{1/2}\right]$$

$`{\displaystyle \frac{49\stackrel{~}{M}_P^2\stackrel{~}{H}^2}{4\pi c}}\left[\stackrel{~}{H}^2\stackrel{~}{f}{\displaystyle \frac{1}{3}}\stackrel{~}{\mathrm{\Lambda }}\stackrel{~}{g}\right].`$ (57)
Here only the last term may be negative. Let us notice that the values of $`\stackrel{~}{f}`$ and $`\stackrel{~}{g}`$ are small, because they are proportional to the ratios of the particle masses and the Planck mass (39). Hence, we can safely conclude that $`D_3`$ is positive if the cosmological constant is not too large. Furthermore, $`D_4=b_4D_3`$ and the stability of the tempered inflation, within our framework, depends on the sign of the same expression which we met in (57)

$`\stackrel{~}{H}^2\stackrel{~}{f}{\displaystyle \frac{1}{3}}\stackrel{~}{\mathrm{\Lambda }}\stackrel{~}{g}.`$ (58)
The positivity of the last expression is a criterion of stability in the last stage of the interpolation regime \- just before the decoupling of sparticles and consequent transition to the unstable regime. The sign of (58) depends on the relation between the masses of the fermion and boson components and on the way of the supersymmetry breaking. The last defines the magnitude of the cosmological constant. If the supersymmetry breaking is spontaneous, the vacuum energy density $`\stackrel{~}{\mathrm{\Lambda }}\stackrel{~}{M}_P^2`$ may be as large as $`M_{SUSY}^4`$, that could be, in principle, much greater than $`M_{}^4`$. Then, in the case of a positive value of $`\stackrel{~}{g}`$ in (39), the inflation would be destabilized by the cosmological constant earlier than at the $`M_{}`$ scale. For the GUT-scale supersymmetry breaking this might be bad, because then the magnitude of the Hubble parameter in the last $`65`$ $`e`$-folds is too big 101010For the lower scale supersymmetry breaking this does not pose a problem.. Thus, for the sake of simplicity we shall assume that $`\stackrel{~}{g}`$ is negative or that the supersymmetry breaking is not spontaneous such that the vacuum energy density does not exceed $`M_{}^4`$. In each of these two cases the transition from the stable to unstable inflation performs, qualitatively, in the same way as for $`\stackrel{~}{\mathrm{\Lambda }}=0`$ and the stability of inflation holds until the supersymmetry breaking at the $`M_{}`$ scale.

In the same framework we can investigate the oscillations of the Hubble parameter during inflation. This is potentially important issue, since it is related to the possibility of a reheating. Let us look again at the equation for the perturbations (53). Typically, the characteristic equation (A4) has one or two couples of complex conjugate solutions

$`\lambda =\alpha \pm i\beta ,`$ (59)
where $`\alpha `$ are positive in the stable case. Also, there may be real solutions. For example, in the case of the MSSM, $`\mathrm{\Lambda }=0`$ and $`\stackrel{~}{f}=10^5`$, we meet one couple of complex solutions $`\lambda _{1/2}1.46\pm 3.47i=\alpha \pm i\beta `$ and two real solutions $`\lambda _33.88`$ and $`\lambda _44.8510^6`$. The solution for the perturbations is

$`y(t)=C_0e^{\alpha t}\mathrm{cos}(\beta t+\phi )+C_1e^{\lambda _3t}+C_2e^{\lambda _4t}.`$ (60)
We will be interested in the effect of $`\beta `$. Of course, for the constant coefficient $`b_k`$ in (54) the oscillations in the solution are invisible because they are suppressed by the exponentials. However, let us remember that the coefficients $`b_k`$ are not constants, for they depend on $`\stackrel{~}{H}=\stackrel{~}{H}(t)`$. We can consider $`\stackrel{~}{H}`$ as a constant only within a small interval $`(t_i,t_i+\mathrm{\Delta }t_i)`$. When we pass from one such interval to another, the value of $`\stackrel{~}{H}`$ changes. The evolution of $`\stackrel{~}{H}`$ will result in the oscillations (60) and they may become significant when $`\stackrel{~}{H}`$ changes faster. Let us evaluate this effect. If we define $`\overline{b}_k=\stackrel{~}{H}^ib_k`$ and $`\overline{\lambda }=\overline{\alpha }+i\overline{\beta },`$ $`\overline{\alpha }_{1/2}`$ to be the solutions of the equation

$`\overline{b}_0\overline{\lambda }^4+\overline{b}_1\overline{\lambda }^3+\overline{b}_2\overline{\lambda }^2+\mathrm{}+\overline{b}_3\overline{\lambda }+\overline{b}_4=\mathrm{\hspace{0.17em}0}.`$ (61)
Then $`\overline{\alpha }_{1/2}`$, $`\overline{\alpha }`$ and $`\overline{\beta }`$ are constants and moreover

$`\alpha _{1/2}=\stackrel{~}{H}\overline{\alpha }_{1/2}H(1\stackrel{~}{f}\sigma )\overline{\alpha }_{1/2},\beta =\stackrel{~}{H}\overline{\beta }H(1\stackrel{~}{f}\sigma )\overline{\beta },\alpha =\stackrel{~}{H}\overline{\beta }H(1\stackrel{~}{f}\sigma )\overline{\alpha }.`$ (62)
The last formula shows that the frequency of the oscillations of $`\sigma (t)`$ and $`H(t)`$ is greatly decreasing during inflation. But, the same concerns also the damping parameter $`\alpha `$. Now, let us look at the plot Fig. 1a. We can divide it into three regions with distinct properties: 1) The initial stage when masses of the particles is irrelevant. Then, the $`H(t)`$ is almost constant and the perturbations are not visible. 2) The intermediate stage of the ”tempered” inflation with a significant effect of the particle masses. In this region $`H(t)`$ changes faster and one can expect oscillations with a relevant amplitude. It is easy to see that the amplitude of the oscillations of $`H(t)`$ has the order of magnitude $`\sqrt{\alpha ^2+\beta ^2}C_0`$. 3) The final region when $`\sigma (t)`$ is approaching a plateau and $`H(t)0`$, where our approximation related to the “sharp cut-off” decoupling obviously fails. In this region the sparticles start to decouple and the universe performs the transition to the non-stable inflation and eventually into the FRW-like evolution.

The numerical analysis of stability has been performed using the numerical simulations method, and we found that all curves rapidly converge to the non-perturbed one in both the first and the last 65 $`e`$-folds of inflation. We have chosen the initial data $`\stackrel{\mathrm{}}{\sigma }(t_i)`$, $`\stackrel{..}{\sigma }(t_i)`$, $`\stackrel{.}{\sigma }(t_i)`$, $`\sigma (t_i)`$ for the equation (38) using the generator of random numbers independent for each derivative and then integrated this equation numerically. The choice of $`t_i`$ has been done in such a way that the corresponding $`\sigma _i=\sigma (t_i)=\sigma _f65`$, where $`\sigma _f`$ corresponds to $`H=M_{}`$. In all cases we took the same initial value $`\sigma (t_i)=1`$. The measure of the deviation from the non-perturbed solution $`H_{np}(t)`$ can be taken, for example, in the form

$$D=\left(\frac{\stackrel{.}{H}\stackrel{.}{H}_{np}}{H_{np}^2}\right)^2+\left(\frac{\stackrel{..}{H}\stackrel{..}{H}_{np}}{H_{np}^3}\right)^2+\left(\frac{\stackrel{\mathrm{}}{H}\stackrel{\mathrm{}}{H}_{np}}{H_{np}^4}\right)^2.$$

The result of the time-consuming calculations is that, independent on the choice of the measure for the deviation, the solution $`\sigma _0`$ presented at Fig. 1 is stable. Some plots, illustrating the behavior of the perturbed solutions are presented at the Fig. 2. The plot of $`H(t)`$ for the last $`65`$ $`e`$-folds at Fig. 3 enables one to see the oscillations 111111Similar plot with oscillations has been obtained using Mathematica program in the course of preparation of the previous article . This plot has not been included into  due to the size limits. which are not visible in the scale of Fig. 1

Finally, the numerical solutions and approximate analytic investigation of the stability issues in this model completely agree between each other. We can observe, at Figure 2, an illustration of the stability of the inflationary solution, and at Figure 3 the oscillations of $`H`$. We can conclude, that our approximate analytic approach is sufficiently robust and is applicable to other problems, like the analysis of the gravitational waves and (hopefully) of the density perturbations.

## 5 On the amplitude of the gravitational waves

The most difficult test for a high derivative cosmological model is the behavior of the cosmological perturbations. In this section we shall present the preliminary results concerning metric perturbations (gravitational waves) during the last 65 $`e`$-folds of the anomaly-induced inflation. Similar problem was previously studied in . In both papers, the analysis was restricted to the anomaly produced by the massless quantum fields, and the small difference in the equations for the perturbations is (in our terms) due to the different choice of the conformal invariant functional $`S_c[g_{\mu \nu }]`$ in Eq. (12). In the Ref. , the derivation of equations for the perturbations has been based on the anomaly-induced effective action and, in principle, this derivation can be generalized for the effective actions Eq. (34) and (35). However, this requires a formulation of the local covariant version of the effective actions (34) and (35), considerable calculations of perturbations for a variable $`H(t)`$ and consequent numerical analysis. At the same time, there is a possibility to evaluate the stability under the metric perturbations in a more economic way, and that is what we are going to discuss here. The detailed study of the gravitational waves will be performed separately.

Let us make some helpful observations about what we really need to know about the metric perturbations. It is well known, that the only phenomenologically important region is the last 65 $`e`$-folds of inflation. All the gravitational waves which were generated before that are not visible in the present-day universe and hence they have no great importance. We shall follow the same procedure which proved reliable in the last section when we investigated the perturbations of $`\sigma (t)`$. Thus, we divide the time interval corresponding to the last $`65`$ $`e`$-folds into the small pieces $`(t_i,t_i+\mathrm{\Delta }t_i)`$. Within each interval $`(t_i,t_i+\mathrm{\Delta }t_i)`$, the time derivatives of the Hubble parameter are negligible compared to its own magnitude, and we can safely take constant $`H`$. In fact, this means that the last 65 $`e`$-folds can be approximated by the massless model, but with a greatly reduced magnitude of the Hubble constant and with much greater magnitude of $`\sigma =\mathrm{ln}a`$ 121212At this instant, we have to admit that the choice of the vacuum for the perturbations can be done in the same way for the effective actions generated by both massless and massive fields. It would be quite interesting to check its validity using some other approach to the derivation of the effective action for massive fields.. As we shall see below, the last point leads to the dramatic difference with the original Starobinsky model .

For the sake of simplicity, let us present the calculations for the $`\mathrm{\Lambda }=0`$ case only. Using the results of , we arrive at the following equations for the perturbations in the massive case

$$b_0\stackrel{\mathrm{}.}{h}+b_1\stackrel{\mathrm{}}{h}+b_2\stackrel{..}{h}+b_3\stackrel{.}{h}+b_4h+$$

$`+n_1e^{2\sigma }^2\stackrel{.}{h}+n_2e^{2\sigma }^2\stackrel{..}{h}+n_3e^{4\sigma }^4h=\mathrm{\hspace{0.17em}0},`$ (63)
where $`hh(t,𝐱)`$ and

$$b_0=a_1+w\sigma (t),b_1=6H[a_1+w\sigma (t)]+2wH,$$

$$b_2=11H^2[a_1+w\sigma (t)]+H^2(cb/2+7w),$$

$`b_3=6H^3[a_1+w\sigma (t)]+H^3(3c3b/2+5w),b_4=12H^4b,`$ (64)
$`n_1=2H[a_1+w+w\sigma (t)],n_2=2[a_1w\sigma (t)],n_3=a_1+w\sigma (t).`$ (65)
We are mainly interested in the dynamics of the perturbations in the last 65 $`e`$-folds of inflation when the magnitude of $`\sigma `$ is at least of the order $`10^4`$. Then, the first observation regarding Eq. (63) is that the terms with the space derivatives are suppressed by the factor of $`e^{2\sigma }`$ and can be safely neglected such that the equation does not depend on the magnitude $`n`$ of the wavenumber. This does not guarantee that the spectrum is flat, because the dependence on $`n`$ may appear from the initial data. For example, in the case of the perturbations of quantum origin the initial data are given by

$$h_0\frac{1}{\sqrt{2n}},\dot{h}_0\sqrt{\frac{n}{2}},\ddot{h}_0\frac{n^{3/2}}{\sqrt{2}},\stackrel{\mathrm{}}{h}_0\frac{n^{5/2}}{\sqrt{2}}.$$
(66)
Even without the $`n`$-dependent terms, the equation (63) is still complicated, but its analysis can be essentially simplified if we remember that $`\sigma 1`$. Taking into account that the $`\beta `$-functions (9), (10), (11) all have the same order of magnitude, one can disregard the $`\sigma `$-independent terms in (64). If we omit also the small term $`b_4`$, the coefficients of the Eq. (63) become constants and we can easily check whether the growing modes are present using the standard methods. Let us make a weaker approximation and regard $`\sigma `$ as a big constant. Then, the Routh-Hurwitz determinants are

$$D_0=b_0,D_1=b_1,D_2=60H^3[a_1+w\sigma (t)]^2,$$

which are positive for any parameters of the theory;

$$D_3=360H^6[a_1+w\sigma (t)]^3+432H^6wb[a_1+w\sigma (t)]^2,$$

which is indeed positive for $`\sigma (t)+a_1/w>6b/5`$. Since $`b`$ is a negative number of the order (at most) 10, this condition is satisfied if we do not include too large negative coefficient $`a_1`$ into the classical action of vacuum (2). Finally,

$$D_4=12bH^4D_3$$

is positive if $`D_3`$ is positive. As we shall see below, for the values $`n<n_c`$, where $`n_c`$ is a cut-off depending on $`M_{}`$, this means that the anomaly-induced inflation is stable under metric perturbations. This fact is very important for the whole model. If we admit that the transition from stable to unstable inflation occurs at $`H=M_{}10^5M_P`$, then this scale will define the maximal amplitude of the gravitational waves and our model fits with the standard restrictions coming from the anisotropy of the CMBR. This conclusion is achieved without fine-tuning of the parameters of the model and looks quite robust. In the consideration above we have neglected the derivative $`\dot{H}`$. This can be justified, for example, using Eq. (41). Remember that, by dimensional reason, $`\dot{H}`$ must be compared to $`H^2`$. A simple calculation gives

$$\frac{\dot{H}}{H^2}\frac{\stackrel{~}{f}}{2(1\stackrel{~}{f}\sigma )}.$$

For the SUSY GUT the magnitude of $`\stackrel{~}{f}`$ is at most $`10^5`$, while the expression $`(1\stackrel{~}{f}\sigma )`$ becomes comparable to $`10^5`$ only at the very end of the stable inflation when the expression (41) and the whole approach based on the “sharp-cut-off” needs to be replaced by a detailed description of the decoupling mechanism for massive sparticles in curved space-time. Therefore, our approximation $`\dot{H}/H^20`$ is reliable in the region where the whole approach is reliable, and can be considered as a qualitative hint for the decoupling region.

The numerical integration of the equation (63) have been performed for various magnitudes of the wavenumber vector $`n`$ and the initial conditions (66). The plots illustrating the stability of $`h(t)`$ are presented at Fig. 4a, where the plots correspond to $`nn_c=0.1`$. As we see, there is a very small difference between these plots such that the spectrum is flat. For the larger values of $`n`$ the plots have greater growth at the initial period (see Fig. 4b) After this initial growth the perturbations decrease exponentially according to our analytical consideration. But, this decrease is too slow and for the large $`n`$ the waves leave the horizon with the amplitude much greater than the initial one. The effect is essentially stronger for smaller values of the parameter $`\stackrel{~}{f}`$. The important observation is that, while the exponential decrease of the perturbations (after the initial growth) does not depend on the choice of the vacuum and the approximation $`\stackrel{~}{H}const`$, the dangerous initial growth may change a lot if we take these aspects into account. Therefore, the final conclusion concerning the amplitude of the perturbations with large $`n`$ can be achieved only after a more complete analysis. The positive result for $`n<n_c`$ give a hope to meet an acceptable behavior of the metric perturbations.

The main difference with the previous analysis of the perturbations  is that here we have to disregard all the $`\sigma `$-independent terms in (63). This changes the structure of the equation and the growing modes disappear. Let us remark that the numerical analysis of the metric perturbations in the first $`65e`$-folds of inflation shows that, say, the perturbations of quantum origin grow up very fast, such that the amplitude increases about 15 orders of magnitude. Indeed, those huge waves can not be observed and their existence does not contradict the CMBR data. However, it is very important that the magnitude of the perturbations is indeed decreasing with respect to the inflationary background metric. As a result, the universe enters into the last $`65e`$-folds of inflation in a very isotropic initial state. All in all, within our framework the anomaly-induced inflation does not meet unsolvable problems with the gravitational waves.

## 6 Conclusions and discussions

We have investigated the stability of the anomaly-induced inflation  with respect to the small perturbations. The stability holds from the initial stage when the quantum fields may be approximately considered massless, until the scale $`M_{}`$, when the most of the sparticles decouple and the inflation becomes unstable. Thus, this model of inflation does not need a fine tuning neither for the initial data nor for a graceful exit to the FRW stage. The equations for the gravitational waves are essentially different for the first and the last $`65e`$-folds of the inflation. In particular, during the last $`65e`$-folds the amplitude of the perturbations with $`n<n_c`$ does not grow up and the spectrum of these perturbations is flat. Thus, for the choice of the scale $`M_{}10^5M_P`$, the anomaly-induced inflation may be consistent with the CMBR constraints.

In the period of inflation previous to the transition to the unstable stage, we met a rapid oscillations of the curvature, very similar to those which were discovered earlier by Suen and Anderson for a similar (technically) higher-derivative model  (see also ). Therefore, the anomaly-induced inflation provides an excellent possibility to describe the reheating during the last stage of the stable inflation and in the transitional period. Finally, let us remark that the next important steps in the development of this model should be the quantitative description of the decoupling process and the corresponding cosmological epoch, and also derivation and analysis of the density perturbations. After these two principal problems will be solved, one can start to investigate different supersymmetric models and their cosmological consequences.

At the present level of understanding the model of anomaly-induced inflation looks vert attractive from both phenomenological (no fine tuning) and theoretical point of view. Let us discuss the last aspect. In the modern high energy theoretical physics there are several drawbacks. One example is the supersymmetry and its low-energy breaking. Even if this breaking is a high energy effect such that $`M_{SUSY}M_F`$, it looks unnatural that there is not a single couple of the observed particles which could be superpartners of each other. This problem becomes even more explicit for the low-energy supersymmetry model like MSSM. But, in the anomaly-induced inflation we see that exactly this spectrum of supersymmetry is favored by the scheme suggested in . Moreover, even the “less natural” soft supersymmetry breaking is favorable, because for the spontaneous supersymmetry breaking there could be problems with the decoupling of heavy sparticles. Another example is the sign universality of the vacuum $`\beta `$-functions (9) and (10) and the variability of sign for (11). As we saw in sections 4 and 5, the conditions $`b<0,w>0`$ are absolutely necessary for the nice stability properties of the inflationary model, while the variability of $`c`$ is very helpful for the successful realization of the transition between stable and unstable regimes . It might happen that all this indicates to a much more general theory which can give us a correct and unified understanding of both vacuum quantum effects and the supersymmetry breaking. From the effective quantum field theory point of view, the anomaly-induced inflation shows that there are chances to link these different issues in a new way and finally approach this general theory.

Acknowledgments. One of the authors (I.Sh.) is very grateful to J. Fabris, I.I. Kogan, V.N. Lukash, J. Sola and A.A. Starobinsky for useful discussions, and especially to J. Sola for the collaboration on the early stage of the work. He is also thankful to his colleagues and friends from the Department of Theoretical Physics at the University of Zaragoza (Spain) for warm hospitality and stimulating interest to his studies. The work of the authors has been supported by the scholarships from CNPq, by the research grant from FAPEMIG and (I.Sh.) by the grant from DGU-MEC (Spain).

Appendix.

This appendix contains a brief summary of the stability issues and especially of the Routh-Hurwitz conditions, which were extensively used in the text of the paper. We refer to the book  for the details. Suppose we know a particular solution $`x_0(t)`$ of the equation

$$x^{(n)}=f(x^{(n1)},x^{(n2)},\mathrm{},x^{},y)$$
$`(A1)`$
where $`x=x(t)`$ and $`x^{(n)}=\frac{d^nx}{dt^t}`$. In order to check the asymptotic stability of this solution under small perturbations we consider $`x(t)=x_0(t)+y(t)`$ and expand (A1) in the first order in $`y(t)`$.

$$b_0y^{(n)}+b_1y^{(n1)}+b_2y^{(n2)}+\mathrm{}+b_{n1}y^{}+b_ny=0,$$
$`(A2)`$
where $`b_0>0`$ for definiteness. Let us consider two possibilities:

i) All $`b_i`$ are constants. Then the problem is reduced to the algebraic characteristic equation

$$b_0\lambda ^n+b_1\lambda ^{n1}+b_2\lambda ^{n2}+\mathrm{}+b_{n1}\lambda +b_n=0,$$
$`(A3)`$
because a general solution of (A2) is a linear combination of the terms $`\mathrm{exp}(\lambda _it)`$ where $`\lambda _i`$ are the roots of Eq. (A3) and the coefficients are polynomials in $`t`$. Hence, if all these roots have negative real part

$$\text{Re}(\lambda _i)<0,$$
$`(A4)`$
the solution $`y_0(t)`$ is asymptotically stable. A useful criterion of (A4) is called the Routh-Hurwitz conditions. The conditions are that the determinants $`D_i`$ must be positive, where $`D_1=b_1`$ and

$$D_k=det\left(\begin{array}{ccccc}b_1& b_3& b_5& \mathrm{}& b_{2k1}\\ b_0& b_2& b_4& \mathrm{}& b_{2k2}\\ 0& b_1& b_3& \mathrm{}& b_{2k3}\\ 0& b_0& b_2& \mathrm{}& b_{2k4}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right),k>1,$$
$`(A5)`$
where $`k=2,3,\mathrm{},n`$ and $`b_j=0,j>n`$. For $`n=4`$ we meet four determinants $`D_{1,2,3,4}`$.

ii) In the cases $`k=1,1`$ of the subsection 3.1 we met the case of the non-constant coefficients $`b_i=b_i(t)`$. But, there is a general theorem (see e.g. ) that the problem of asymptotic stability for such case can be reduced to the one with the constant coefficients $`\overline{b}_i`$, provided that there are finite limits $`b_i(t)\overline{b}_i`$ when $`t\mathrm{}`$. As we saw, this is exactly what happens for $`k=1,1`$, where $`\overline{b}_i`$ correspond to the $`k=0`$ case.

Figure captions

> Figure 1. The plots representing numerical solutions $`\sigma (t^{})`$ (where $`t^{}=t/16\pi `$) of the equation (38) for the toy model with $`\stackrel{~}{f}=10^6`$ and MSSM particle content. The plot (a) is the same as the one obtained in . In the region $`\dot{\sigma }<M_{}`$ the two plots (a) and (b) are equivalent. The plot (b) is almost identical to the parabola (43). At the top point of the plot (b) the value of $`\sigma `$ equals 1000000, while in case of parabola $`\sigma `$ equals 1000001.

> Figure 2. The graphics illustrating the behavior of the perturbation measure $`D`$ for the first (a) and last (c) 65 $`e`$-folds of inflation. The plot (b) slightly differs from (a), because of the presence of the matter (dust) with the density $`M_P^4`$ at the initial point of inflation. One can see that the stabilization takes a very short time in all cases.

> Figure 3. The plot of numerical solution for $`H(t)`$, in the region where the oscillations are explicitly visible.

> Figure 4. The samples of the plots of the numerical solutions $`h(t)`$ in the last $`65`$ $`e`$-folds of inflation and initial spectrum (66) with the normalization $`h(0)=1`$ and $`\stackrel{~}{f}=10^5`$.


format: text

text:
# A Test of The Source Galerkin Method

## 1 Introduction.

The Source Galerkin method is being developed as a flexible alternative to Monte Carlo approaches to solving Quantum Field Theories.

### 1.1 Overview of the Method

To illustrate the application of the Source Galerkin approach to solving Quantum

Field Theories we consider a theory defined by a Lagrangian $``$. All information of the theory is contained in the generating functional

$$Z(j)=𝒟\varphi \mathrm{exp}\left[(+j\varphi )d^dx\right].$$
(1)
This satisfies a Schwinger-Dyson set of differential equations. A procedure to determine an approximate solution of these equations is specified by our SG technique. To start, a trial solution or an ansatz must be introduced. In general, it is constructed from a predefined set of trial functions. The proper choice of trial functions allows us to take advantage of symmetries of the theory as well as of other known analytical properties of the generating functional. A simple example of an approximate solution is:

$$Z_a=\mathrm{exp}\left[j_xG_{xy}j_y+j_wj_xH_{wxyz}j_yj_z+\mathrm{}\right].$$
(2)
A set of residues is obtained by substituting this ansatz into the Schwinger-Dyson equations. The undefined coefficients of the trial solution are determined by solving a system of nonlinear equations obtained by projecting the residues on the trial functions and requiring these projections to be equal to zero. The procedure outlined above guarantees that the error associated with the approximate solution converges to zero in the mean as the number of members in the set of trial functions goes to infinity.

### 1.2 Motivation

In order to demonstrate the validity of the method, it was applied to the $`O(3)`$ nonlinear $`\sigma `$ model . This model is asymptotically free and is a useful toy model for approaching Non-Abelian gauge theories. First order solutions were obtained, but during our attempts to extend this work to higher orders it became evident that a better understanding of the general properties of the SG method is necessary. Theorems guarantee that Galerkin methods produce an approximate solution, which converges to the exact solution as the number of terms goes to infinity. In practice, the approximate solution can only have a fairly small number of parameters. Therefore, we must ensure that high accuracy and rapid convergence can reasonably be achieved with only a few terms in the ansatz. The fact that the error approaches zero in the mean implies that the accuracy of the result is heavily affected by the choice of inner product. We attempt to investigate the effectiveness of the Source Galerkin approach by using it to solve $`\varphi ^4`$ theory in zero dimensions. We study its performance for several choices of trial functions and scalar products.

## 2 Ultralocal Model

After introducing an external source $`j`$, we define the ultralocal model by the following lagrangian

$$=\frac{g}{4}\varphi ^4+\frac{\mu }{2}\varphi ^2j\varphi .$$
(3)
Theoretical solutions for this theory can be obtained by expressing generating functional as a power series

$$Z_a(j)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{1}{k!}a_kj^k.$$
(4)
The coefficients of this series are given by

$`a_{2n}`$ $`=`$ $`\left(2n1\right)!!{\displaystyle \frac{U(n,t)+\left(1\right)^n\rho U(n,t)}{U(0,t)+\rho U(0,t)}},`$
$`a_{2n+1}`$ $`=`$ $`{\displaystyle \frac{2n!!}{n!}}\left(t\right)^n{\displaystyle \frac{V(n+\frac{1}{2},t)}{V(\frac{1}{2},t)}}\alpha {\displaystyle \frac{t^{\frac{1}{2}}e^{\frac{t^2}{4}}}{U(0,t)+\rho U(0,t)}}.`$
Here $`U`$ and $`V`$ are parabolic cylinder functions, coefficients $`\rho `$ and $`\alpha `$ fix boundary conditions.

### 2.1 Numerical Solution

In analogy with the theoretical solution, a truncated power series of order $`N`$ can be chosen as a trial solution. After obtaining the residues by substituting this expression into the Schwinger-Dyson equations, we need to eliminate source dependence by projecting them on the trial functions. Two possible ways to define scalar product are presented below.

$$\underset{c}{\overset{c}{}}j^nR\left(j\right)dj,n=0..N3$$
(5)
and

$$\underset{\mathrm{}}{\overset{\mathrm{}}{}}j^nR\left(j\right)\mathrm{exp}[\frac{j^2}{ϵ^2}]dj,n=0..N3.$$
(6)
Both definitions are localized in the region close to $`j=0`$ since eventually we want to set the source to zero in order to compute physical values. Using these definitions produces similar results as shown in figure 1. However, the second expression is more general and it is easier to extend to calculations with space-time dimensions.

Inspection of these results shows that the error in determination of coefficients $`a_i`$ increases rapidly at high orders. In spite of this, the generating functional can be determined with accuracy as high as several parts in $`10^7`$ for a significant range of j. However, it must be noted that correct determination of the proper range of integration or a value of the parameter $`ϵ`$ is crucial for achieving high accuracy. The SG approach by itself does not provide an algorithm for setting these parameters. In a real problem when the exact result is not available it is necessary to have some external source of information which allows us to determine the correct values of the parameters.

### 2.2 Solution in Terms of Hermite Polynomials

A set of Hermite polynomials can be defined by

$$H_n\left(\xi \right)=\left(1\right)^n\mathrm{exp}\left\{\frac{\xi ^2}{ϵ^2}\right\}\frac{d^n\mathrm{exp}\left\{\frac{\xi ^2}{ϵ^2}\right\}}{d\xi ^n}.$$
(7)
These polynomials are orthonormal under the following inner product

$$\underset{\mathrm{}}{\overset{\mathrm{}}{}}H_k\left(x\right)H_l\left(x\right)\mathrm{exp}\left(\frac{x^2}{ϵ^2}\right)𝑑x=\delta ^{kl}.$$
(8)
This suggests a way to improve our implementation of the Source Galerkin procedure for the ultralocal theory. Trial functions can be constructed as a linear combination of Hermite polynomials of different orders, while the expression shown above is used to define the scalar product. This choice leads to significant simplification of the system of equations which give the coefficients in the ansatz. Figure 2 shows the dependence of relative error in determination of values $`a_i`$ on the parameter $`ϵ`$. From this graph we observe that there is a significant range of values of the parameters in which numerical results remain stable. We can use this fact to resolve the problem outlined in the end of previous section. In a real problem, it would be possible to set all the parameters without any knowledge of theoretical solution just by finding solutions for all possible values of parameters and then choosing a solution which is stable within the required accuracy for some range of the parameters. The accuracy of determination of the generating functional using this method is consistent with the power series ansatz approach. This implies that it is limited only by error introduced by the method used to solve the system of equations produced by the Source Galerkin procedure and by the numerical precision of the hardware used for computation.

## 3 Conclusion

This investigation demonstrates that very high precision computations can be performed using the Source Galerkin procedure. On the other hand, the method is very sensitive to the choice of parameters. Application of this method to a more general set of problems will require development of some consistent method of determination of proper values of parameters. Current results suggest that use of orthogonal trial functions may simplify this task. Source Galerkin technique is still under development. Additional work is needed to fully investigate it’s properties and to apply it to various field theories.

## 4 Acknowledgements:

We wish to thank R. Easther, D. Ferrante, Z. Guralnik and S. Hahn for helpful discussions. This research was supported in part by DOE grant DOE grant DE-FG02/19ER40688-(Task D).