Beyond the linear analysis of stability in higher derivative gravity with the Bianchi-I metric

TL;DR

This paper investigates the stability of higher derivative gravity with Bianchi-I metrics, confirming that linear stability extends to the nonperturbative level for small anisotropies, supporting the robustness of previous linear analyses.

Contribution

It provides a comparison between linear and nonperturbative stability solutions in higher derivative gravity with Bianchi-I metrics, validating the linear stability results.

Findings

Linear stability aligns with nonperturbative asymptotic behavior.

Small anisotropies do not lead to ghost instabilities.

Linear analysis is reliable for certain gravitational perturbations.

Abstract

The study of stability of gravitational perturbations in higher derivative gravity has shown that at the linear level the massive unphysical ghost is not generated from vacuum if the initial seed of metric perturbation has frequency essentially below the Planck threshold. The mathematical knowledge indicated that the linear stability is supposed to hold even at the nonperturbative level, but in such a complicated case it is important to perform a verification of this statement. We compare the asymptotic stability solutions at the linear and full nonperturbative levels for the Bianchi-I metric with small anisotropies, which can be regarded as an extreme, zero frequency limit of a gravitational wave. As one should expect from the combination of previous analysis and general mathematical theorems, there is a good correspondence between linear stability and the nonperturbative asymptotic behavior.

Equations

$$\displaystyle S=\int d^{4}x\,\Big{(}-\frac{M_{P}^{2}}{16\pi}\,R+a_{1}C^{2}+a_{2}R^{2}\,\Big{)}.$$

$$C^2=R_{\mu \nu \alpha \beta }^2-2R_{\alpha \beta }^2+\frac{1}{3}{R}^2.$$

$$d{s}^2=d{t}^2-a_1^2(t)d{x}^2-a_2^2(t)d{y}^2-a_3^2(t)d{z}^2.$$

$$a_1(t)$$

$$a_2(t)$$

$$a_3(t)$$

$$g_{\mu \nu }=e^{2\sigma (\eta )}{\bar{g}}_{\mu \nu },$$

$$\displaystyle R=e^{-2\sigma}\big{[}\bar{R}-6(\sigma^{\prime})^{2}-6\sigma^{\prime\prime}\big{]}.$$

$$\bar{R}$$

$${\bar{C}}^2$$

$$g_{ik}=-{\delta }_{ik}+h_{ik},$$

$$\displaystyle h_{ik}\,=\,-\,\,\mbox{diag}\,\big{(}\beta_{+}+\sqrt{3}\beta_{-},\,\,\beta_{+}-\sqrt{3}\beta_{-},\,\,-2\beta_{+}\big{)}.\mbox{\qquad}$$

$$\mathcal{L}$$

$$\displaystyle 72a_{2}\,\Big{[}\sigma^{(4)}-2\,\sigma^{\prime\prime}\,\big{(}3\,\sigma^{\prime 2}+\,\beta_{-}^{\prime}{}^{2}+\beta_{+}^{\prime}{}^{2}\big{)}$$

$$\displaystyle 24a_{1}\big{(}8\beta_{\mp}^{\prime}{}^{2}\,\beta_{\pm}^{\prime\prime}+16\beta_{\pm}^{\prime}\,\beta_{\mp}^{\prime}\,\beta_{\mp}^{\prime\prime}+24\beta_{\mp}^{\prime}{}^{2}\,\beta_{\pm}^{\prime\prime}-\beta_{\pm}{}^{(4)}\big{)}$$

$$\displaystyle 72a_{2}\,\Big{[}\,\sigma^{(4)}+12\dot{\sigma}^{2}\,\ddot{\sigma}+4\,\ddot{\sigma}^{2}$$

$$\displaystyle 144a_{2}\Big{\{}\ddot{\beta}_{\pm}\,\big{(}\,2\,\dot{\sigma}^{2}+\ddot{\beta}_{\mp}^{2}+3\,\dot{\beta}_{\pm}^{2}+\ddot{\sigma}\,\big{)}$$

$$\dot{\sigma }={\sigma }^0+ϵ{\sigma }^1+\cdots$$

$${\dot{\beta }}_{\pm }=0+ϵ{\beta }_{\pm }^1+\cdots ,$$

$$d_t\mathbf{y}=\frac{d}{dt}\mathbf{y}=\mathbf{f}(\mathbf{y}),$$

$$\displaystyle d_{t}\big{[}\mathbf{y}^{0}+\epsilon\mathbf{y}^{1}+\cdots\big{]}\,=\,\mathbf{f}(\mathbf{y}^{0})+\epsilon\nabla\mathbf{f}(\mathbf{y}^{0})\mathbf{y}^{1}+\cdots,\mbox{\qquad}$$

$$d_t{\mathbf{y}}^1=\nabla \mathbf{f}({\mathbf{y}}^0){\mathbf{y}}^1.$$

$$\dot{\sigma }=H,\dot{H}=Q_1,\dot{Q_1}=Q_2,$$

$$\dot{{\beta }_{\pm }}=x_{\pm },\dot{x_{\pm }}=y_{\pm },\dot{y_{\pm }}=z_{\pm }$$

$$\displaystyle\dot{Q_{2}}\,=\,-\Big{[}\,12H^{2}\,Q_{1}+4\,Q_{1}^{2}$$

$$+H(6x_{+}{y}_{+}+6x_{-}{y}_{-}+7Q_2)$$

$$\displaystyle+2\big{(}\,y_{+}^{2}+y_{-}^{2}+x_{+}\,z_{+}+x_{-}\,z_{-}\big{)}\,\Big{]}$$

$$\displaystyle-\frac{M_{p}^{2}}{96\pi a_{2}H_{0}^{2}}\Big{[}2H^{2}+Q_{1}+x_{+}^{2}+x_{-}^{2}-2\Omega_{\Lambda}$$

$$\displaystyle-\frac{1}{2}\,\Omega_{0}\,e^{-3\sigma(1+\omega)}\,(1-3\omega)\Big{]},$$

$$\displaystyle\dot{z_{\pm}}\,=\,6\frac{a_{2}}{a_{1}}\,\Big{\{}x_{\pm}\,\Big{[}6\,H^{3}+3\,H\,(\,x_{+}^{2}+x_{-}^{2}\,)+7\,H\,Q_{1}$$

$$\displaystyle+2x_{\mp}\,y_{\mp}+Q_{2}\,\Big{]}+y_{\pm}\,\big{(}\,2\,H^{2}+x_{\mp}^{2}+3\,x_{\pm}^{2}+Q_{1}\,\big{)}\Big{\}}$$

$$\displaystyle-\Big{\{}x_{\pm}\big{[}6\,H^{3}-16x_{\mp}\,y_{\mp}+Q_{2}+7H\,Q_{1}-24H(x_{+}^{2}+x_{-}^{2})\big{]}$$

$$\displaystyle+6\,H\,z_{\pm}+y_{\pm}\,\big{(}11\,H^{2}-8\,x_{\mp}^{2}-24x_{\pm}^{2}+4\,Q_{1}\big{)}\Big{\}}$$

$$\displaystyle+\,\frac{3M_{p}^{2}}{4\pi a_{1}H_{0}^{2}}\big{(}y_{\pm}+3H\,x_{\pm}\big{)}.$$

$$\dot{{\beta }_{\pm }}=x_{\pm },\dot{x_{\pm }}=y_{\pm },\dot{y_{\pm }}=z_{\pm },$$

$$\displaystyle\dot{z_{\pm}}=\frac{6a_{2}}{a_{1}}\Big{\{}x_{\pm}\big{[}6\dot{\sigma_{0}}^{3}+7\dot{\sigma_{0}}\ddot{\sigma_{0}}+\sigma_{0}^{(3)}\big{]}+y_{\pm}\big{(}2\dot{\sigma_{0}}^{2}+\ddot{\sigma_{0}}\big{)}\Big{\}}$$

$$\displaystyle-\Big{\{}x_{\pm}\big{[}6\dot{\sigma_{0}}^{3}+\sigma_{0}^{(3)}+7\dot{\sigma_{0}}\ddot{\sigma_{0}}\big{]}+6\,\dot{\sigma_{0}}\,z_{\pm}$$

$$\displaystyle+y_{\pm}\,\big{(}11\,\dot{\sigma_{0}}^{2}+4\,\ddot{\sigma_{0}}\big{)}\Big{\}}\,+\,\frac{3}{4\pi a_{1}\,H_{0}^{2}}\,M_{p}^{2}\,\Big{(}\,y_{\pm}+3\,\dot{\sigma_{0}}\,x_{\pm}\Big{)}.$$

Full text

Bianchi-I solutions, higher derivative gravity, massive ghosts, stability, nonlinear analysis

I Introduction

There is well known controversy between renormalizability of quantum gravity and the problems which are caused by the introduction of higher derivatives, which are capable to provide this renormalizability Stelle . The theory with sufficiently general higher derivatives always has massive unphysical ghosts in the spectrum, making physical interpretation of such a theory problematic. In the presence of ghosts the vacuum state is not stable, and even Minkowski space may decay into Planck-mass ghost plus the gravitons with huge overall energy which is compensating the negative energy of the ghost.

Indeed, the presence of the ghost in the spectrum of the theory does not necessary mean that there should be such a particle “alive”. It might happen, e.g., that there is an unknown physical principle which forbids the concentration of gravitons with Planck energy density, resolving the mentioned puzzle with Minkowski space HD-Stab ; DG , and also providing the stability of a qualitatively similar, low curvature space-times. Certain arguments which support this expectation have been given in the recent papers GW-Stab ; HD-Stab ; GW-HD-MPLA . In a perfect agreement with the previous works on the evolution of gravitational waves on the deSitter background, star81 ; asta ; hhr , we have found that these waves do not have growing amplitudes, regardless of the presence of higher derivatives. The situation was analysed in the context of ghosts in HD-Stab , where it was shown that there are no growing modes also in other cosmological backgrounds, if the initial frequency of the gravitational wave is much smaller than the Planck scale. On the opposite, in case of Planck-scale frequencies there is an expected explosion of gravitational waves. Our interpretation of this situation in HD-Stab was that the presence of the ghost in the spectrum of the theory does not necessary means that there is a ghost as a real particle. For the low-energy frequencies of the gravitational waves the positive energy modes don’t form a Planck-density distribution and then the ghost can not be created from vacuum. This solution of the problem is certainly incomplete, because i) quantum gravity is supposed to work at all frequencies, even over-planckian ones;   ii) The linear stability guarantees non-linear perturbative stability from the mathematical point of view, but it does not look sufficient from the point of view of Physics, because the exponential instabilities are expected at the non-linear level Woodard .

The item i) has been addressed in PP , where we have shown that, at least for the cosmological background, if the cosmological solution corresponds to the rapidly expanding universe, the explosive behaviour of the gravitational waves does not last for long, and after that the metric perturbations get stabilized. The reason is that the wave equation includes the wave vector ${{\bf k}}$ only in the combination ${{\bf q}}=\textstyle{\frac{{\bf k}}{a(t)}}$, such that the physical frequency of the wave is decreasing as $1/a$. Of course this is not the complete solution of the problem, but just a useful hint on how the problem can be eventually solved. What is still needed is certainly the physical principle explaining why gravitons can not accumulate with the over-Planck energy density on a weak gravitational background, and how this principle may be violated by the fast expansion of the universe.

In the present work we address the point ii) and check out whether the situation with stability changes when we go beyond the linear perturbations level. In fact, we are able to get even the non-perturbative results, but not for the usual gravitational waves. Instead, we shall consider the evolution of anisotropies in the framework of the Bianchi I cosmological metric. Since the pioneering work Kasner:1921zz , the Bianchi I metric have been extensively studied as a model of anisotropic homogeneous cosmology. For cosmologic solutions and stability in fourth derivative gravity, see recent works Barrow:2006xb ; Toporensky:2016kss ; Muller:2017nxg .

With respect to an arbitrary perturbations of the metric our approach means the following two restrictions: (a) small amplitude of the perturbations; (b) zero frequencies of the perturbations. In what follows we perform numerical analysis of the dynamics of anisotropies under these two assumptions.

The paper is organized as follows. In the next Sec. II the equations for the Bianchi-I metric in the fourth derivative gravity are derived in Misner parametrization Misner:1967uu ; Misner:1969hg . Before starting the numerical analysis of the full and linearized version of these equations, in Sec. III we present a brief survey of the mathematical knowledge on the subject of stability in the systems described by differential equations. Namely, we discuss to which extent the stability with respect to linear perturbations defines the behavior of the system at the nonperturbative. In Sec. IV we present the results of numerical analysis including comparison of linear and full versions of equations. Finally, in Sec. V we draw our conclusions and discuss possible extensions of the present work.

II Dynamical equations

The theory of our interest has the classical action

[TABLE]

Here $M_{P}$ is the Planck mass, while other parameters $a_{1}$ and $a_{2}$ are arbitrary dimensionless constants. $R$ and $C^{2}$ are, respectively, the Ricci scalar and the square of Weyl tensor,

[TABLE]

According to the recent work Barrow:2006xb , every vacuum solution of Einstein field equations is also a solution of the theory (1). However, since there are higher derivatives, the theory (1) can develop strong instabilities which are not present in general relativity. These instabilities represent our main interest in what follows.

In a comoving and synchronous frame, the Bianchi-I anisotropic metric is

[TABLE]

One can switch to a more useful parametrization, introduced by Misner in Misner:1967uu ; Misner:1969hg , in which there is a separation between the functions of time responsible for expansion $\sigma(t)$ and shear of the universe $\beta_{\pm}(t)$ respectively,

[TABLE]

In what follows the term anisotropies will refer to the functions $\beta_{\pm}$. The trivial case $\beta_{\pm}=0$ corresponds to an isotropic metric. A usefulness of Misner parametrization resides in the possibility of perform a local conformal transformation

[TABLE]

where the conformal time $\eta$ is defined by the relation $\,dt=e^{\sigma(\eta)}d\eta$. The fiducial metric $\bar{g}_{\mu\nu}$ is given by (3) with $\,\sigma(t)\equiv 0$. Under a conformal transformation, the Weyl-squared part of the action (1) is expressed only in terms of the metric $\bar{g}_{\mu\nu}$, while Ricci scalar transforms as

[TABLE]

It is easy to check that $\sqrt{-\bar{g}}=1$ and the expressions for $\bar{R}$ and $\bar{C}^{2}$ are

[TABLE]

In these expressions the prime stands for the derivative with respect to conformal time.

Let us remember that we regard the anisotropy parameters as a truncated part of the gravitational wave, or the gravitational wave with zero frequency. The gravitational wave of our interest is supposed to be created by quantum fluctuations GW-Stab , and if it does not experience fast growth due to the presence of ghosts, its amplitude remains very small. This is our main assumption and we need to know whether it is violated by the dynamics of the gravitational wave or, in the truncated case, of the anisotropies. Thus, consider the physically most interesting case when the anisotropy parameters in Eq. (3) are small, $\,|\beta_{\pm}|\ll 1$. Then one can write the space components of the metric in the form

[TABLE]

It is easy to see that the trace of the last expression is zero, $\,\delta^{ik}h_{ik}=0$, exactly as in the case of the gravitational wave, also in both cases we have two degrees of freedom.

Another desired similarity would be a transverse nature of the wave. However, in the case of Bianchi-I metric this feature can not be verified, because the perturbation in (7) is dependent only on time, and there is no wave vector. Therefore there is no complete correspondence between (7) and the gravitational wave, and we can speak only about a qualitative similarity between the two types of the perturbations. At the same time, since the Ostrogradsky instabilities which are expected in the higher derivative theories Ostrog (see Woodard:2015zca for a recent review) appear due to the higher derivatives in time, we can expect that the data obtained by using Bianchi-I metric will provide a useful hint for the general situation with the stability of metric perturbations in the higher derivative theories. Since the wave vector is zero in the case of (7), we can expect that, according to the results of HD-Stab , the classical isotropic solutions will be stable in the linear approximation. The Bianchi-I metric offers a possibility to have an independent check of these results and, most relevant, to go beyond the linear approximation.

In terms of the new variables, discarding superficial terms and taking into account that in Bianchi-I case all metric components depend only on time and not on the spatial coordinates, the Lagrangian of the action (1) becomes

[TABLE]

It is worth noting that in the limit of general relativity $a_{1,2}\to 0$ and after a rescaling anisotropies, that doesn’t affect the dynamics of the conformal factor, we recover the conventional Lagrangian for the gravitational waves beyond the horizon Mukhanov:2005sc ; Gorbunov:2011zzc . This means that, at least in the linear order, the Bianchi-I model under consideration can be seen as a zero-frequency approximation of the equation for the gravitational waves. Thus we shall assume that this correspondence holds beyond the linear order and regard the Bianchi-I as a simplest version of the equation for the gravitational wave.

It is easy to see that that Lagrangian expression has terms which are second and fourth order in conformal time derivatives. It is useful to show explicitly the unit of time $\eta_{0}$. The dynamical equations can be obtained by taking the variational derivatives of the action with the Lagrangian (II). The presence of isotropically distributed matter, radiation or cosmological constant does not affect the equations for $\,\beta_{\pm}$ Reis:2017bjf ; hervik2007 , but only changes the equation for $\sigma$ through the trace of the energy-momentum tensor. We will only consider a perfect fluid with linear equation of state defined by the constant $\omega$ which is assuming the values $\frac{1}{3}$, [math] and $-1$ for radiation, dust and cosmological constant, respectively. Taking variational derivatives with respect to $\sigma(\eta)$ and $\beta_{\pm}(\eta)$ and adding the matter part, we arrive at the equations

[TABLE]

and

[TABLE]

Here the primes mean the derivative with respect to the conformal time measured in the units of $\,\eta_{0}$. Eq. (9) corresponds to the variation with respect to $\,\sigma$ with the perfect fluid contribution, where $\Omega_{0}$ is the relative energy density of matter or cosmological constant. The sum of $\Omega$ and the contribution of higher derivative terms is equal to one identically. The Eqs. (10) describe the nonlinear dynamics of anisotropies.

We can also express the dynamical equations in terms of physical time through the relation $\,dt=e^{\sigma(\eta)}d\eta$. The results are

[TABLE]

and

[TABLE]

Here the dots mean derivative with respect of dimensionless time $\tau=H_{0}\,t$, where $H_{0}$ is the Hubble - Lemaître parameter measured at some instant of time. The set of Eqs. (11) and (10) or (11) and (12) represent systems of three coupled ordinary differential equations of the fourth order.

Besides the Einstein space solutions in vacuum (with cosmological constant), there are no much chances to find an exact solution of this system, and this is not our purpose in the complicated case with higher derivative terms included. Instead, we shall explore the stability of the cosmological (homogeneous and isotropic) solutions, corresponding to $\,\beta_{\pm}=0\,$ and the $\,\sigma(t)=\sigma_{0}(t)\,$ given by some cosmological solutions.

An important point concerns the choice of the background solution $\,\sigma_{0}(t)$. Let us start from a few preliminary observations. The first one is that in the action (1) the term $a_{1}C^{2}$ (regardless being most relevant for the tensor perturbations and massive ghosts) does not affect the dynamics of the the conformal factor and therefore the one of $\,\sigma_{0}(t)$. Thus when we choose this functions, we do not need to take the Weyl-squared term into account. Second, our main target in the previous works on the cosmological stability in the presence of massive ghosts was the low-energy cosmological solutions. Thus, the canonical approach would be to ignore also the $a_{2}R^{2}$-term as being Planck-suppressed (the last means we consider such solutions for which $\,\left|a_{2}R^{2}\right|\ll\left|M_{P}^{2}R\right|$ in the action and the corresponding hierarchy in the equations of motion), and consider the classical radiation-dominated and dust-radiation solutions only. Let us stress that this hierarchy can be assumed only for the background $\,\sigma_{0}(t)$. For the perturbations such as gravitational waves, the run-away solutions are capable to easily destroy this hierarchy. The main subject of the present work is to explore the effect of non-linearities in this possible breaking in the framework of the simple Bianchi-I based model.

The third point is that we can easily extend the low-energy region for the background up to the inflation scale, just taking the $a_{2}R^{2}$-term into account. According to the available set of observational and experimental data, this term is the main ingredient of the Starobinsky model star , that is exactly the most successful phenomenologically model of inflation. In order to achieve this success star83 , the value of $a_{2}$ should be chosen at about $\,5\times 10^{8}$. Then the inflationary solution corresponds to the slowly decreasing Hubble parameter, with an approximately linear dependence $H(t)$. Then, since our ultimate interest is the dynamics of the gravitational waves with the initial frequencies much greater in magnitude than $H$ (and at the same time much below the Planck scale HD-Stab ), it is a very good approximation to regard $H$ as a constant. Thus, we can safely consider, instead of the linear $H(t)$, the constant $H$ and derive it from the classical cosmological constant. All in all, we arrive at the situation where the main features of our model can be explored taking the three simplest examples of $\,\sigma_{0}(t)$, namely cosmological constant-, radiation- and matter-dominated classical solutions.

Let us repeat that the main advantage of the Bianchi-I metrics is that the Eqs. (9) and (9) or (11) and (12) are relatively simple and can be explored numerically even at the non-perturbative level. Thus we get a chance to check by direct calculation whether the mathematical statements about the general relation between linear stability and the nonperturbative asymptotic behavior, which were used in HD-Stab and PP , are correct. However, before going to numerics we shall give a brief survey of the mentioned mathematical statements in the next section.

III Asymptotic series expansion for singular perturbation

Since our intention is to compare the linear approximation for the anisotropies with the nonperturbative numerical solution, it makes sense to briefly review the general mathematical theorems which cover the relation between first order stability and nonperturbative behavior in the systems described by differential equations.

In the zero-order case functions $\sigma$ and $\beta_{\pm}$ are approximated by $\sigma_{0}(t)$ and zero, because in the background solutions there are no anisotropies, by assumption. This fact motivates to explore the general solution of the system of equations Eqs. (11) and (12) in the form of asymptotic series expansion

[TABLE]

where $\epsilon$ is a small parameter, which one can easily implemented into the perturbations (7).

Eqs. (11) and (12) can be rewritten in the mathematically standard form as a system of twelve autonomous ordinary differential equations

[TABLE]

where the vector $\mathbf{y}$ includes $\sigma$, $\beta_{\pm}$ and also first, second and third derivatives of these functions. Substituting into this system the expansion (13), we arrive at the equations for the power series

[TABLE]

where $\nabla\mathbf{f}(\mathbf{y}^{0})$ is a Jacobian of the function $\mathbf{f}$ calculated on the background (unperturbed) solution $\mathbf{y}^{0}$. In order to solve this system we equate terms with the same order in $\epsilon$. This procedure is well known in Singular Perturbation Theory wasow87 .

Let us note that the order zero in $\epsilon$ corresponds to the equation $d_{t}\mathbf{y}^{0}=\mathbf{f}(\mathbf{y}^{0})$, that is satisfied for the background under consideration. Then the first order approximation corresponds to the linear differential equation

[TABLE]

Our main purpose is to compare the solution of this equation with the one for the complete version (15). For instance, let us assume that for the certain choice of initial conditions (small deviations from the background, as we explained above), linear system (16) does not show growing modes, but only those which asymptotically vanish or oscillate without growing amplitude in the limit $t\to\infty$. Then, under smoothness hypotheses on the dependence on the small parameter $\epsilon$, the first order approximation $\mathbf{y}^{0}+\epsilon\mathbf{y}^{1}$ is of the order $\epsilon$ close to the solution of the complete system $d_{t}\mathbf{y}=\mathbf{f}(\mathbf{y})$ wasow87 .

Finally, we can quote the following two theorems concerning sink equilibrium points, which can be found in the well-known book on differential equations hirsch1974differential :

Theorem 1.   Assume that the system $d_{t}\mathbf{y}=\mathbf{f}(\mathbf{y})$ possesses a sink in the point $\tilde{\mathbf{y}}$, i.e., there exists a constant $c>0$, such that all eigenvalues $\lambda_{i}$ of the Jacobian $\mathbf{f}(\tilde{\mathbf{y}})$ satisfy $Re(\lambda_{i})<-c$. Then all the solutions starting in some neighborhood of the point $\tilde{\mathbf{y}}$ converge to $\tilde{\mathbf{y}}$ exponentially.

Theorem 2.   If the system $d_{t}\mathbf{y}=\mathbf{f}(\mathbf{y})$ possesses a stable equilibrium in $\tilde{\mathbf{y}}$, then all eigenvalues $\lambda_{i}$ of the Jacobean $\mathbf{f}(\tilde{\mathbf{y}})$ have non positive real part of the eigenvalues $Re(\lambda_{i})\leq 0$.

Coming back to our problem of exploring Eqs. (11) and (12), we know that in the linear approximation there are no growing modes for the frequencies below the Planck-order threshold HD-Stab ; GW-HD-MPLA . This is certainly true for the zero frequency modes, which correspond to the Bianchi-I model. Thus we can claim that the condition of the Theorem 2 are satisfied and, therefore, the conditions of the Theorem 1 are also satisfied. Hence we can expect a good qualitative correspondence between the dynamics of anisotropies in the linear approximation and within the full nonperturbative consideration. In the next section we check this conclusion by using numerical methods and now let us discuss how these well-known theorems can be applied to evaluate the regions where one can expect the validity of the linear approximation.

First of all, let us construct the presentation (14) for the non-linear system formed by (11) and (12). For this end we introduce the new variables

[TABLE]

Then the first order equations equations include (17),

[TABLE]

and

[TABLE]

The first order version of linearized system consists from

[TABLE]

and

[TABLE]

In order to estimate the radius of the region where the linearization procedure is valid for the ordinary differential equations written in the form (14), one needs to go into details of the proofs of the theorems mentioned above. In both cases the proofs are based on the Taylor expansions around the equilibrium point $\mathbf{y_{0}}$ in the form

[TABLE]

where $\mathbf{J}$ and $\mathbf{H}$ are the Jacobian and Hessian operators of the function $\mathbf{f}$ evaluates on the background solution $\mathbf{y_{0}}$. Remember that at the equilibrium point $\mathbf{f}(\mathbf{y_{0}})=\mathbf{0}$ by definition. Thus, the equation above can be rewritten for the perturbations as

[TABLE]

The theorems cited above are valid in the region where the terms of the higher order are negligible (or possibly vanish under certain change of variables) in a small neighborhood of $\mathbf{y_{0}}$. This means that the linear approximation ceases validity when linear and quadratic terms are of the same order of magnitude. A rough estimate for the region where the linear approximation is valid is

[TABLE]

where the Euclidean norm $|\cdot|$ is used for vectors and the operator norm $\|\cdot\|$ follows the standard definition and can be calculated using Riesz representation (see e.g. kreyszig78 ) as

[TABLE]

One can note that it is not possible to apply this formula to the linearized model as the Hessian tensor will be singular. This is a natural situation, because the criteria (III) and (25) are intended to compare linear and non-linear cases.

In order to calculate the radius given in 24 for the system (19) we need just to evaluate the norms of Jacobian and Hessian operators at the equilibrium point. The numerical simulations based on Eqs. (24) with (III) and (25) has been performed in the radiation and dust models, using the dimensionless units with $M_{P}=1$. The results were equal for the tested versions with $a_{1}=\pm 1$ and $a_{2}=5\times 10^{8}$. In both models we met the radius $R=(1/3)\cdot 10^{-9}$. It is interesting that the sign of $a_{1}$ did not make any difference for the radius of validity of the linear approximation $R$, regardless of the critical importance of the same sign for the asymptotic stability, as we will discuss in the next section.

An interesting observation is in order. After we submitted the first version of this work to arXiv we learned about a similar investigation Salvio-19 . The results of numerical analysis in this work concern the non-linear case and are qualitatively the same as ours, that are also close to those of the earlier paper Toporensky:2016kss , which did not link the study of the dynamics of anisotropies with the problem of massive ghosts in higher derivative gravity. The correspondence between the three independent investigations are certainly adding an extra safety to our conclusions. At the same time the results of Salvio-19 include the growing solution for the the initial conditions with relatively large first derivatives of the anisotropies. This output may look as a contradiction with our interpretation of Bianchi-I perturbations as a zero frequency gravitational wave. The analysis presented in this section shows that this is not a correct interpretation. The frequency is still zero, but in this case we have the situation when the initial conditions correspond to the point which is out of the region satisfying the condition (24). As we have discussed, out of this region we can not expect correspondence between linear and non-linear approximations.

IV Linear and non-linear numerical solutions

In this section we present the numerical solutions of differential equations (11) and (12) in both linear and full version. The first part requires the linearization. Let us note that in this section we exclusively work with set of Eqs. (11) and (12) in terms of dimensionless physical time.

As we have explained above, the linearization is performed around isotropic cosmological solutions, which means null values for anisotropies and the well-known cosmological solutions of general relativity $\,\sigma_{0}(\tau)\,$. It is easy to check that at the linear level the perturbations for $\,\sigma(\tau)\,$ and anisotropies completely decouple. Thus in the linearized case one can restrict consideration by the equations for anisotropies, which have the form

[TABLE]

The free parameters of the systems are Hubble - Lemaître parameter at the reference time instant $H_{0}$ and the coefficients $\,a_{1}\,$ and $\,a_{2}$. The theory with $\,a_{1}\,>\,0\,$ manifest instabilities for anisotropies as we know from the more general gravitational wave solutions HD-Stab (see also more detailed discussion in GWprT ). For the same of completeness we present the corresponding plots in Figs. 1 for the cosmological constant - dominated and matter - dominated backgrounds.

The radiation case is very similar to these two and hence will not be included here. In what follows, we consider only negative values of $\,a_{1}$.

The examples of the results of numerical analysis can be seen in the figures presented below. The qualitative behavior is pretty much the same for any choice of initial data which we tried. The values for the plots which we selected are specified at the Captions of the figures. In all cases, the initial conditions for $\,\beta_{+}(\tau)\,$ for both linear and non-linear equations which we show in the plots are $\,\beta_{\pm}(0)\,=\,0$,  $\dot{\beta_{\pm}}(0)\,=\,0.01$,  $\ddot{\beta_{\pm}}(0)\,=\,-0.001$,  $\beta_{\pm}^{(3)}(0)\,=\,0.0001$. Furthermore in order to shorten the numerical procedure, the value of Hubble - Lemaître parameter has been taken as $\,H_{0}\,=\,10^{-2}\,M_{p}$. In the figures we present the plots of numerical solutions for $\,\sigma(\tau)\,$ and anisotropies. In the last case we show only $\,\beta_{+}(\tau)\,$ solutions, because it turns out that both anisotropies $\,\beta_{\pm}(\tau)\,$ have similar behaviour, which may differ only due to the choice of initial conditions and do not define the asymptotic behaviour. The time $\tau$ is measured in units of $1/H_{0}$, where we choose $H_{0}=0.01M_{P}$ for the sake of convenience of numerical analysis and plotting the figures.

In the first set, illustrated in Figs. 3, 4 and 5 the system of nonlinear equations have initial conditions for $\,\sigma(\tau)\,$ which are the same as for isotropic radiation - dominated universe in general relativity. Linearization is done around $\,\sigma_{0}(\tau)\,$ of isotropic radiation dominated universe.

The second set of Figs. 6,  7 and 8 illustrates the solutions for the background of $\,\sigma_{0}(\tau)\,$ corresponding to the matter - dominated universe.

The last cases are shown in Figs. 9, 10 and 11, they correspond to equations for the variation of conformal factor and anisotropies on the background of isotropic solution in the universe dominated by cosmological constant.

Let us conclude this section by repeating that we have also checked other choices of initial data and the results are always qualitatively the same as in the plots shown above. In general there is a very good correspondence between linearized Bianchi-I system and the dynamics of gravitational waves with low frequencies from one side, and the linearized and non-perturbative treatments from another side.

V Conclusions

We have explored the time dependence of anisotropies in the Bianchi-I model with fourth derivatives, which can be seen as a zero-frequency approximation for the gravitational waves in the model (1). Qualitatively we observe from the plots presented in the Figures that in all cases there is no qualitative difference between the behaviour of linearized and non-perturbative systems, exactly as it should be in accordance with the standard mathematical results cited in Sec. III.

In all cases which we were analysed, the dynamics of both linearized and general systems does not show instabilities related to the presence of higher derivatives, exactly as one should expect from the previous considerations of the gravitational waves from one side HD-Stab and the mentioned mathematical theorems from another side. Since Bianchi-I can be regarded as a zero-frequency approximation to the gravitational waves dynamics, we gain a strong reasons to expect the absence of explosive exponential type instabilities for the gravitational waves, even in the nonperturbative regime.

For the cases of radiation-dominated and dust-dominated background solutions the numerical results confirm show that for the values $\,a_{1}\,=\,-1\,$ and $\,a_{2}\,=\,1\,$ the numerical solutions of $\,\sigma(\tau)\,$ asymptotically tend to the isotropic ones with the same matter contents. At the same time, for larger value $\,a_{2}\,=\,100\,$ we can note stronger deviation between linear and nonperturbative regimes. This effect should be expected much stronger for the phenomenologically optimized value $\,a_{2}\,\approx\,5\times 10^{8}$, required for the successful Starobinsky inflation star ; star83 .

In general, we confirmed the expectations of HD-Stab and PP concerning the correspondence between linear and general nonlinear results. It would be certainly interesting to extend the analysis in several directions. For instance, to include the cases of the background cosmological metrics with strong curvature, such that the effect of higher derivatives on the background should be taken into account. Regardless of that this case is not expected to give great surprises (the reason is that the large $a_{2}$ is known to increase the value of $H_{0}$, in the first approximation), this check has to be done. In fact, the solutions for more complicated cosmological backgrounds would be an interesting issue to explore. A much more challenging problem is to consider more complicated anisotropic solutions, with a non-zero frequencies. Such an investigation would require more serious calculation, but in some cases it does not look impossible. Anyway, the results of the present work show that we have strong reasons to believe to the validity of the first-order perturbations if they show the strong signs of asymptotic stability.

Acknowledgements

S.C.R. is grateful to CAPES for supporting his Ph.D. project. I.Sh. was partially supported by CNPq (grant 303893/2014-1) and FAPEMIG (project APQ-01205-16).