Anisotropic cosmological solutions in $R + R^2$ gravity

In this paper we investigate the past evolution of an anisotropic Bianchi I universe in $R+R^2$ gravity. Using the dynamical system approach we show that there exists a new two-parameters set of solutions that includes both an isotropic"false radiation"solution and an anisotropic generalized Kasner solution, which is stable. We derive the analytic behaviour of the shear from a specific property of $f(R)$ gravity and the analytic asymptotic form of the Ricci scalar when approaching the initial singularity. Finally we numerically checked our results.


Introduction
Inflation is a generic intermediate attractor in the direction of expansion of the universe, and in the case of pure R 2 gravity it is an exact attractor [1]. However, it is not an attractor in the opposite direction in time. Thus, if we are interested in the most generic behavior before inflation, more general anisotropic and inhomogeneous solutions should be considered. We know from General Relativity (GR) that already anisotropic homogeneous solutions help us much in understanding the structure of a generic space-like curvature singularity. Thus, a natural question is to investigate anisotropic solutions in the R + R 2 gravity, too. In light of the latest Cosmic Microwave Background constraints by PLANCK [2], the pioneer inflationary model based on the modified R + R 2 gravity (with small one-loop corrections from quantum fields) [3] represents one of the most favourable models. It lies among the simplest ones from all viable inflationary models since it contains only one free ajustable parameter taken from observations. Also it provides a graceful exit from inflation and a natural mechanism for creation of known matter after its end, which is actually the same as that used to generate scalar and tensor perturbations during inflation. This theory can be read as a particular form of f (R)-theories of gravity which, in turn, is a limiting case of scalar-tensor gravity when the Brans-Dicke parameter ω BH → 0, and it contains an additional scalar degree of freedom (scalar particles, or quasi-particles, in quantum language) compared to GR which is purely geometrical. The existence of a scalar degree of freedom (an effective scalar field) is needed if we want to generate scalar (matter) inhomogeneities in the universe from "vacuum" fluctuations of some quantum field [4,5]. Such generalizations of the familiar Einstein-Hilbert action have been also studied as an explanation for dark energy and late-time acceleration of the universe's expansion [6,7,8,9] and to include quantum behaviour in the gravitational theory [10]. As is already very well known, through the Gauss-Bonnet term which in four dimensions is a surface term, the most general theory up to quadratic in curvature terms is of the type R + R 2 + C abcd C abcd , where C abcd is the Weyl tensor. The investigations of this type of models began with [11,12,13,14]. After them, many authors have been analyzed the cosmological evolutions of such a model [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]. A particular attention is given to the asymptotic behavior in [30,31,32,33]. The addition of a Ricci square term creates a richer set of solutions, and in particular, in [1,34,35], has been shown that the cosmic no-hair theorems no longer hold. Quadratic theory like R + R 2 gravity, is a particular case of the more general quadratic type, and has higher order time derivatives in the equations of motion and this leads to appearance of solutions which have no analogs in GR. One of such solutions corresponds to the scale factor a ∼ √ t behavior, that coincides with a radiation dominated solution in GR. In quadratic gravity, this solution represents a vacuum solution which is also a stable past attractor for Bianchi I model and probably for all Bianchi models [35]. The other solution, being an analog of GR solution with matter with equation of state p = (γ − 1)ρ, is a ∼ t 4/3γ (instead of usual a ∼ t 2/3γ in GR), and it can be naively considered as a solution which would describe the last stages of a collapsing universe when quadratic terms dominate. However, this solution appears to be a saddle, so a collapsing universe for a general initial condition ends up with a vacuum "false radiation" regime, in principle not possible in GR. The a(t) ∝ √ t behavior near singularity in the R + R 2 model does not mean that the R 2 term behaves as radiation generically. This behavior is specific only for the purely isotropic case, as it will be shown in the present paper (and even in the isotropic case the late-time behavior is different: a(t) ∝ t 2/3 modulated by small high-frequency oscillations). Neither does it behave as an ideal fluid in the anisotropic case.
When shear is taken into account the situation becomes more complicated. Vacuum solutions exist in GR also, and in the simplest case of a flat anisotropic metrics (that is the case analyzed in the present paper) this is the famous Kasner solution [36]. On the other hand, studies of cosmological evolution in a general quadratic gravity (which includes apart from R 2 the Weyl tensor square term in the quadratic part of the action) indicate that the isotropic vacuum "false radiation" still exists and, moreover, it is an attractor [35]. Kasner solution is also a solution in quadratic gravity. Due to complicated nature of dynamics near the Kasner solution a generic trajectory could end up in Kasner or isotropic solution depending on initial conditions [37].
So, for these reasons, the stability and the full behavior of quadratic theories of gravity is still subject of investigation and one powerful tool to address these problems is the dynamical system approach [38] which allows to find exact solutions of the theory through the determination of fixed points and gives a description of the evolution of the system, at least at qualitative level.
Despite the obvious fact that the general quadratic gravity at the level of the action includes R + R 2 gravity (which can be obtained setting the coefficient before Ricci square term to zero) it is not so at the level of corresponding equations of motion for the universe model in question. The reason is that the number of degrees of freedom in a general quadratic gravity is bigger than in R + R 2 theory (which, on its turn is bigger than in GR). That is why we can not simply put corresponding constant to zero in cosmological equations of motion. This means that cosmological evolution of a flat anisotropic universe in R + R 2 gravity needs a special investigation which is the matter of the present paper.
The paper is organized as follows: In section 2, we present the basic equations of the model. In section 3 we describe schematically the strategy adopted to obtain the correct degrees of freedom and then we analyze the dynamics of R 2 -gravity both in the vacuum case and in the case with matter; we find exact solutions and determine their stability. In section 4, we derive the analytic behavior of the shear using a general line element. Finally, section 5 contains a summary of the results and conclusions. 2

System under consideration
The gravitational action considered in our analysis is the following where g is the determinant of the metric, G the Newton constant and β a parameter. This theory can be interpreted as a particular form of f (R) gravity. Observations tell us that the dimensionless coefficient β 16πG is very large, ≈ 5 × 10 8 . This follows from the fact that its expression in terms of observable quantities, in the leading order of the slow-roll approximation, is β 16πG = where P ζ (k) is the power spectrum of primordial scalar (adiabatic) perturbations, N is both ln k f /k and the number of e-folds from the end of inflation, k f being the wave vector corresponding to the comoving scale equal to the Hubble radius at the end of inflation (k f /a(t 0 ) is slightly less than the CMB temperature now): see e.g. [39]. For the R + R 2 inflationary model, P ζ ∝ N 2 , so β is a constant indeed. Note also that β = 1 6M 2 where M is the scalaron mass after the end of inflaton (and in flat space-time, too). On the other hand, the coefficient of the Weyl square term (that is present in general quadratic model of gravity) in the Lagrangian density generated by one-loop quantum gravitational effects is not expected to be so large. Typically it is of the order of unity (or even significantly less due to small numerical factors) multiplied by the number of elementary quantum fields. Thus, there exists a large range of the Riemann and Ricci curvature where the R 2 term dominates while the contribution from the Weyl square term is still small. For this reason, anisotropic solutions preceding the inflationary stage may be studied using the same R + R 2 model up to curvatures of the order of the Planck one.
Metric variation of the theory in (1) gives the following fields equations where H (1) ab is the contribution coming from the variation of R 2 term. Let us emphasize that every Einstein metric satisfying R ab = g ab Λ is an exact solution of (2). This implies that all vacuum solutions of GR are also exact solutions of the quadratic theory (1). Any source that satisfies ∇ c T ca = 0, can be consistently added to the right hand side of (2). As anticipated, a powerful tool to provide exact solutions of quadratic theory of gravity, is the dynamical system approach which allows for the determination of fixed point and for a qualitative description of the global dynamics of the system. It is particularly suited for the study of the dynamics of anisotropic spacetimes [38], like spatially homogeneous Bianchi metrics. In this case we can write the line-element as where the i, j indices refer to the spatial part and ω j is a triad of one-forms where C a bc are the spatial structure constants of the Bianchi group, and depend only on time. They are usually defined as where the values of the symmetric matrix n ab and the vector a b define the various Bianchi models. In our case, where we focus on Bianchi I metric, we have n ab = 0 and a b = 0. Defining the time-like vector u a = (H, 0, 0, 0), which satisfies the normalization condition u a u a = −1, and the projection tensor h ab = g ab + u a u b , we can define the relevant kinematical quantities where σ ab is the symmetric shear tensor (σ ab = σ (ab) , σ ab u b = 0, σ a a = 0), ω ab is the vorticity tensor (ω ab = ω (ab) , ω ab u b = 0) andu a is the acceleration vector. θ is the volume expansion, and it is related to the Hubble parameter by In our analysis we consider a cosmological model where the shear is diagonal and is defined as and since we will consider spatially homogeneous spacetimes, the time-like vector is geodesicu a = 0 with zero vorticity ω ab = 0, being normal to the time slices.
We consider a perfect fluid source with no anisotropic pressures so the energy-momentum tensor is and it can be decomposed schematically as where w is the equation of state (EoS) parameter. In order to have a system of autonomous first order differential equations we divide the shear σ ± , given in (9), and density parameters by appropriate powers of H, defining in this way the new dimensionless expansion-normalized variables (ENV) where ρ/H 2 = T 00 is the energy density. The rest of the ENV are zero since we are restricting to the Bianchi I case. The time evolution of the sources follow directly from the conservation of the energy momentum tensor (∇ b T ab = 0) and from the definition itself of (12),Ω 4 The fact that we have higher order theory of gravity, requires the introduction of additional ENVs, which reflect the higher order time derivatives in the equations of motion, as firstly done in [35] Q 1 =Ḣ H , According to their own definitions, these ENV must satisfy the following differential equationṡ So now we have all the ingredients the compute the evolution of our theory. The complete dynamical system is given by the equations (13), (15) and by the differential equations shown in the Appendix A.

Generalized anisotropic solutions
Now we start from this particular line element for the spacetimes For vanishing cosmological constant, and near the singularity when τ → 0, the Einstein tensor for the above line element goes like G ab ∼ 1/τ 2 , so it becomes negligible in comparison to the H ab ∼ 1/τ 4 given in (3). By direct substituting the line element (16) into the field equations (2) for vacuum source, a purely algebraic equation is obtained when B = 1 3βH 2 → 0 The set of solutions (17) can be parametrized using two angles ψ and φ where ψ = [0, 2π] and 0 < φ < π, and lies in the surface of an ellipsoid shown in figure 1. In this surface are contained both generalized Kasner solution and the isotropic solution.
The expansion normalized variables for the line element (16) read .
5 Figure 1: Ellipsoid in the parameter space p 1 , p 2 and p 3 given by equation (18).
The solution space given by eq. (17) can be written in a more compact form using the variables u = p 2 1 + p 2 2 + p 2 3 and s = p 1 + p 2 + p 3 , with τ → 0 or with respect to ENV with B → 0, as The solution of (20) is easily obtained as This compact way of writing the equation, is particular suitable to check the solutions: in fact it can be easily seen that Kasner (s = 1 and u = 1) and the isotropic vacuum (s = 3/2 and u = 3/4) are both particular solutions of this equation 1 . We also remember that, in terms of the ENV, generalized Kasner's solution is given by Q 1 = −3, Σ 2 + + Σ 2 − = 1, and the isotropic vacuum solution by Q 1 = −2 and Σ + = Σ − = 0. And that both of these solutions belong to the solution set given by (20).
The Ricci scalar can be written in terms of the ENV, and in terms of variables s and u, like such that the solution given in (22) , (21) as long as B = 0, results in By looking at (3), it is not difficult to convince ourselves that zero Ricci scalar (R = 0) is in fact the asymptotic solution of eq. (2). If near the singularity when τ → 0, the most important contribution to (2) comes from H ab given in (3), the field equation can be approximated by On the other hand for vanishing cosmological constant and absence of classical sources, in a non perturbative picture in the sense that the Einstein tensor G ab is not disregarded, it behaves as an effective source for field equations. And even though it diverges at the singularity (20), the following ratios of the effective pressures to energy densities obtained directly from (16) do have a well defined limit. The trace indicates that at the singularity, given by (20), the effective EOS parameter behaves as in radiation

Stability analysis
In the dynamical system approach, the field equations are re written with respect to the ENV, such that the solutions are fixed points. In particular, the solution space described in the previous section constitute an invariant set of fixed points. The linearization around the fixed points reveals the local stability of the theory. In fact, since all eigenvalues λ i ≥ 0, this solution set is an attractor to the past, as all trajectories to the future deviate exponentially from this solution set. Stability with and without matter source is going to be addressed, and the presence of matter is irrelevant for sufficiently big shear.

Obtaining the dynamical system
In order to describe the evolution of the correct degrees of freedom, in this section we will describe the strategy that we have adopted in order to simplify the system of equations of motion. From the E 11 , E 22 and E 33 equations in (2) we can isolate the variable related to the higher order time derivative Q 2 ; then we find a system of 3 differential equationṡ , , where f is a generic function of all the remaining ENV. Form the E 00 component of (2), we obtain a constraint equation that, as we will see, will be important to check the stability of the numerical evolution of the dynamical system. By doing a linear combinations of the field equations (26), we can obtain two additional constraints that read as So now, from the constraints (27) and (28), it is possible to write three algebraic equation for Q 2 , Σ +1 and Σ −1 , that will be function of the remaining variables and we write schematically as If now, we consider the ENV related toσ ± , defined as Σ ±2 =σ ± H , we can use its definition to derive the equationΣ ±1 = Σ ±2 − Q 1 Σ ±1 . Using (29), it is now possible to derive the equation for Σ ±2 Substituting the last equation and eq. (29) into the original dynamical system equations we finally obtain also the equation forQ 2 = f 1 (Σ ± , Q 1 , Ω, B).
Then the complete dynamical system is described by the following equationṡ where in our analysis the last equation will be integrated numerically to be compared with the algebraic relation Q 2 (Σ ± , Q 1 , Ω, B) contained in (29). Moreover we will use one of the constraints to obtain a conserved quantity to numerically check the stability of our results. The last equation is not a dynamical degree of freedom, but just an artifact to check numerically the solutions.
Looking at the above set of equations, it can be noted that there is only one additional dynamical degree of freedom compared to General Relativity, which is the first equation forQ 1 . This can be easily understood by remembering that, through a conformal transformation, this gravitational theory is equivalent to GR plus a scalar field [41].
As we have described above, the linearization of the dynamical system (31) around the solution (21) gives rise to the following eigenvalues, that will be discussed in the next subsections.

Pure geometric modes
As a first case we consider the vacuum case (without the matter modes). In this case the stability of the system is characterized by the following eigenvalues
Excluding the zero eigenvalues, it is interesting that this solution is an attractor to the past in the presence of cosmic substance for any initial values for ⌃ + and ⌃ as long as the EOS parameter 1 < w < 1/3. For 1/3 < w < 1, the solution will still be an attractor to the past for sufficiently big initial values of ⌃ + and ⌃ such that 2 is positive.

Conclusions
In the present paper we have considered the past attractor solution for the evolution of the flat anisotropic Universe in R + R 2 gravity. Our results in combination with already known solutions indicate that the properties of Universe evolution near a cosmological singularity change significantly with taking into account anisotropy and/or modifications of gravity. Indeed, the evolution of isotropic Universe is determined solely by the matter equation of state. When anisotropy is taken into account, this isotropic solution becomes future asymptotic solution, while vacuum Kasner solution becomes a past attractor (unless the stiff fluid case with Jacobs solution).
In quaratic gravity without anisotropy new vacuum isotropic solution ("false

Matter modes
Allowing perturbations in the matter sector (⌦ ⇤ , ⌦ m ) we have the following eigenvalues, Excluding the zero eigenvalues, it is interesting that this solution is an attractor to the past in the presence of cosmic substance for any initial values for ⌃ + and ⌃ as long as the EOS parameter 1 < w < 1/3. For 1/3 < w < 1, the solution will still be an attractor to the past for sufficiently big initial values of ⌃ + and ⌃ such that 2 is positive.

Conclusions
In the present paper we have considered the past attractor solution for the evolution of the flat anisotropic Universe in R + R 2 gravity. Our results in combination with already known solutions indicate that the properties of Universe evolution near a cosmological singularity change significantly with taking into account anisotropy and/or modifications of gravity. Indeed, the evolution of isotropic Universe is determined solely by the matter equation of state. When We explicitly checked that asymptotically Q 1 < −4.37. In a) it can be seen that the Ricci scalar decreases which is consistent with eq. (61) when Q 1 < −3. In b) it is explicitly shown that the numerical solution approaches eq. (19), as it should for all points of the past attractor. Even if not reported in the plot, the constraint has been checked to be satisfied numerically.
Except for the two zero eigenvalues this solution is an attractor to the past. These two zero eigenvalues appears naturally because in fact we have two-dimensional set of fixed points So, according to this theory the universe began as a generalized non-isotropic solution as shown in Figure 2. There it can be seen in panel b) that an arbitrary initial condition approaches eq. (19) which defines this solution set. Even if non reported the constraint was numerically verified.

Matter modes
Allowing perturbations in the matter sector (Ω Λ , Ω m ) we have the following eigenvalues, Excluding the zero eigenvalues, it is interesting that this solution is an attractor to the past in the presence of cosmic substance for any initial values for Σ + and Σ − as long as the EoS parameter −1 < w < 1/3. For 1/3 < w < 1, the solution will still be an attractor to the past for sufficiently big initial values of Σ + and Σ − such that λ 2 is positive. We focused on solutions which are attractors to the past, however there can be other solutions [35].

Analytic behavior
Now, starting from a remarkable property of any f (R) gravity, obtained in [14], it is possible to determine the dynamical evolution of the shear as an exact, analytical result. Considering a Lagrangian like the variation with respect to the metric gives the following field equation where f is used to denote derivative wrt to R, and represents derivative wrt to the proper time.
T ab is the energy momentum tensor of some matter source. Since we are considering a spatially homogeneous spacetime, with corresponding sources, the following combinations are zero a) T 22 −T 33 = 0, b) T 11 − T 22 /2 − T 33 /2 = 0. In the same way, taking the same combinations of the left hand side of the field equations (35), we have We can see that these equations admit as solutions Now we can connect the constants (C ± ) with the parameters p 1 , p 2 and p 3 . As in [14] we can write the scale factors given in (16) as such that assuming i g i = 1, we haveġ Taking into account (7), (8) and (9) we find thaṫ Remembering that derivatives with respect to proper time are related to derivatives with respect to dynamical time by d dτ = dt dτ d dt = H d dt , we can solve the first of (43) which, when substituted into (39) and (43), giveṡ This relation found here is the same of (10) in [14] with the constant C ± related to C 1 , C 2 and C 3 by When we consider the asymptotic solution discussed above we finally find that the relation with the coefficient p 1 , p 2 and p 3 is the following where again s = p 1 + p 2 + p 3 and (22) must be satisfied. C ± can also be expressed as It is now possible to understand the space of solutions close to the singularity, when B → 0, as long as the solution stays near the asymptotic solution given in (20), which must be fulfilled since as explained in section 3.1 the solution is an attractor to the past. First of all, by taking the trace of (2) bearing in mind (3), gives the following equation for R in absence of sources Near the singularity, when t → −∞, Q 1 = −3/s = const. and B = e 6t/s 3βH 2 0 and eq. (49) has an analytical solution given by Bessel function of the first type J a and the Bessel function of the second type Y a , which when written with respect to proper time τ , exp(3t/s) ∝ τ is given by Asymptotically, as τ → 0 and B → 0, (50) simplifies as and when s = 1 the solution is while for s = 1, which means Q 1 = −3, it is where all the Cs are constants. This asymptotic behavior of R can be substituted into (39) for the particular theory analyzed hitherto, for which we have f = 1 + 2βR, , .
When Q 1 > −3, which corresponds to s > 1, this expression gives at the singularity (t → −∞) and since 2 + Q 1 + Σ 2 − + Σ 2 + = 0 with Q 1 = −3/s = const., it results in the following asymptotic form for the Ricci scalar and since the Ricci scalar must be real then Q 1 < −2 which gives s < 3/2. We can also obtain the constant C since we know that, interchangeably when s < 1 or Q 1 < −3, the asymptotic solution set (22) must continue to be a past attractor, see section 3.1. This attractor has constant well defined values for Q 1 , Σ ± satisfying (19) and when Q 1 < −3 this will only occur if there is a particular cancellation in the denominator of (55) f → 0 giving a well defined limit for Σ ± at the singularity We have the final asymptotic form for the Ricci scalar which is valid through 0 < s < 3/2. Through (48), this last expression can be written as The numerical behavior of the Ricci scalar is shown in Figure 2 in panel a). There it can be seen that the asymptotic constant value R −1/(2β) is not reproduced exactly in the plot. The reason for that is due to the fact that, in the dynamical system described in the appendix, the following denominator 4Σ 2 + + 4Σ 2 − + 4Q 1 + B + 8 occurs in all equations and it vanishes on the attractor set. Although it is not possible to check numerically the asymptotic time evolution of the Ricci scalar it is possible to see in Figure (2) panel a) that the Ricci scalar decreases asymptotically as it is expected for Q 1 < −3.

Conclusions
In the present paper we have considered the past attractor solution for the evolution of the flat anisotropic universe in R+R 2 gravity. Our results, in combination with already known results, indicate that the properties of the universe evolution near a cosmological singularity change significantly taking into account anisotropy and/or modifications of gravity. Indeed, the evolution of isotropic universe is determined solely by the matter equation of state. When anisotropy is taken into account, this isotropic solution becomes future asymptotic solution, while generalized vacuum Kasner solution becomes a past attractor (except for stiff fluid with Jacobs solution).
In general quadratic gravity without anisotropy new vacuum isotropic solution ("false radiation" solution) is stable to the past. The anisotropic case has instead has two sub-cases, because general quadratic corrections to the gravitational action has two independent terms which can be chosen as proportional to squares of scalar curvature and the Weyl tensor. In a general situation, when these two terms are of the same order, the dynamical system describing the universe past evolution has both "false radiation" and generalized Kasner solution as attractors (the latter is, more precisely, a saddle-node fixed point). So the nature of cosmological singularity (isotropic or anisotropic) depends on initial conditions imposed.
However, since the R + R 2 inflationary model is observationally well motivated, and we have argued in Sec. II above that there exists a large range of the Riemann and Ricci curvatures where the anomalously large R 2 term dominates the Einstein term R, while a "normal-size" Weyl squared term is still small, one can expect that new solution appears. It has two parameters (so it is a two-dimensional set of solutions), and includes both isotropic "false radiation" and generalized anisotropic Kasner solution (which is a one-dimensional set) as subsets. Moreover, in some sense, it interpolates between them, because it is possible to construct line of solutions with one end being isotropic solution and the other end being a point in the generalized Kasner set.
All this intermediate points disappears when the correction proportional to Weyl square is added to the action, leaving only isotropic and generalized Kasner solutions and this represents one of the main results of our paper. However, since R 2 inflation model is observationally well motivated we can neglect the coefficient in front of the Weyl term, and we can expect that the two-dimensional set of solutions discussed in this paper could be a good approximation for realistic models in quadratic gravity.
In the present paper we have restricted the analysis to flat metrics. However positive spatial curvature could, in principle, destroy this regime and generate more complicated behavior similar to the Belinsky-Khalatnikov-Lifshitz (BKL) [42] singularity in General Relativity. We leave this problem for future analysis.