xAct Implementation of the Theory of Cosmological Perturbation in Bianchi I Spacetimes

This paper presents a computational algorithm to derive the theory of linear gauge invariant perturbations on anisotropic cosmological spacetimes of the Bianchi I type. Our code is based on the tensor algebra packages xTensor and xPert, within the computational infrastructure of xAct written in Mathematica. The algorithm is based on a Hamiltonian, or phase space formulation, and it provides an efficient and transparent way of isolating the gauge invariant degrees of freedom in the perturbation fields and to obtain the Hamiltonian generating their dynamics. The restriction to Friedmann--Lema\^itre--Robertson--Walker spacetimes is straightforward.


I. INTRODUCTION
Our understanding of the physics of the early universe is connected to cosmological perturbation theory. This framework describes the evolution of cosmological spacetimes, together with matter and gravitational perturbations propagating thereon. The description of perturbations is technically complicated, mainly due to the fact that general relativity is a gauge theory. Due care is needed to separate gauge artifacts from physical effects. A natural and conceptually clean strategy is to work with gauge invariant fields, which are combinations of matter and gravitational perturbations that are left invariant under diffeomorphisms, or coordinate transformations, at the desired order in the perturbative expansion. This strategy was first implemented by Bardeen [1] by expanding Einstein's equations to first order in perturbations, and it has been extensively used since then. An alternative and more geometric treatment can be obtained by working in the Hamiltonian, or phase space formulation of general relativity [2]. Here, the phase space is equipped with four constraints-the so-called scalar and vector constraints-which are the generators of gauge transformations. Working at leading order in perturbations, one defines the gauge invariant fields as those combinations that remain unchanged under the transformations generated by the linearized constraints or, in simpler words, that Poisson-commute with them. Hence, at the practical level, the task of finding gauge invariant perturbations reduces to finding combinations of linear perturbations whose Poisson brackets with the linearized constraints vanish. Furthermore, this procedure is equivalent to finding an appropriate canonical transformation, in such a way that the search for gauge invariant fields reduces to solving a Hamilton-Jacobi-like equation for the generating function of the transformation-with the added advantage that this equation becomes a set of simple algebraic equations for the unknown coefficients. This strategy was first implemented by Langlois [3] for Friedmann-Lemaître-Robertson-Walker (FLRW) cosmologies (see Refs. [4][5][6][7] for a recent discussion), and it has been extended to anisotropic Bianchi I spacetimes in [8] (see also [9,10] for a previous analysis starting from Einstein equations).
Although more systematic and geometric, the phase space derivation of linear cosmological perturbations is still extremely tedious and lengthy, in particular in the presence of anisotropies. In fact, the complexity of the calculations is one of the main barriers that researchers find when entering this field. The goal of this paper is to alleviate this issue by introducing a computational algorithm to find gauge invariant linear perturbations in Bianchi I spacetimes, and the equations of motion they satisfy in the symbolic language of Mathematica. Our algorithm is based on the package xPert written for Mathematica, which is embedded on the tensor algebra packages xTensor of the xAct distribution [11][12][13] (they are available under the General Public License). A previous, and more general implementation of perturbation theory using similar tools was published in [14]; it is more general because it applies to all orders in perturbation theory and to any spatially homogenous spacetimes in any gauge-it does not deal, however, with gauge invariant perturbations. (See also [15,16] for the application to perturbations on spherically symmetric spacetimes.) The merit of our algorithm is that we combine a series of analytical and computational techniques to make the problem of finding gauge invariant variables at leading order in perturbations tractable. We find that the phase space formulation helps enormously in these respects and, when combined with carefully thought numerical algorithms, it gives rise to a computational infrastructure of great utility to researchers in this field. This manuscript is accompanied by a Mathematica notebook, publicly available at [17], that contains a step by step implementation of the algorithms presented here. Furthermore, we have complemented the analysis with another code, based on the C programming language and publicly available in [18], to solve the equations of motions of gauge invariant perturbations and to compute observable quantities in the cosmic microwave background (CMB) starting from suitable initial data, although the details of this code are not spelled out in this paper.
This article is organized in such a way that the procedure has been separated in small steps that are, one by one, introduced theoretically and implemented computationally. Namely, section II introduces Bianchi I spacetimes and linear perturbations thereon in the Arnowit-Deser-Misner (ADM) formalism [2]. Section III introduces the generalization of the scalar-vector-tensor decomposition commonly used in FLRW backgrounds. Section IV introduces gauge invariant fields, and in section V we study their dynamics. Section VI contains a short summary and some concluding remarks.

A. Summary of the theory
We start with Einstein's gravity minimally coupled to a scalar field Φ with potential energy density V (Φ) and no anisotropic stresses. We assume the spacetime manifold to be M = R × M 3 , where M 3 has R 3 topology. In the following, we will restrict ourselves to a finite volume V 0 relative to an auxiliary flat Euclidean metric δ ij defined in M 3 . 1 We will adopt a Hamiltonian formulation following ADM [2]. In this formalism, elements of the phase space are made of four real fields (Φ( x), P Φ ( x), h ij ( x), π ij ( x)) defined in M 3 , where Latin indices i, j run from 1 to 3. Here, h ij ( x) is a Riemannian metric that describes the intrinsic spatial geometry of M 3 , and π ij ( x) its conjugate momentum. The non-vanishing Poisson brackets between these fields are Dynamics is generated by the Hamiltonian which is a combination of first class constraints, and N ( x) and N i ( x) (the lapse and shift, respectively) play the role of Lagrange multipliers. Concretely, S( x) is the scalar constraint and V i ( x) are the vector or diffeomorphism constraints. In terms of the canonical variables, they have the form where κ = 8πG, and D i is the covariant derivative compatible with the spatial metric h ij , h its determinant, and (3) R its Ricci scalar curvature. Given N ( x), N i ( x), and a solution to Hamilton's equations, h ij ( x, t), the spacetime metric takes the form where t is a time variable that labels each space-like hyper-surface M 3 (t), and x i are spatial coordinates on them.
Let us now focus on the sector of the phase space of general relativity that is made of Bianchi I geometries together with small inhomogeneous perturbations. This is commonly done by considering curves γ[ǫ] in the ADM phase space that pass through the Bianchi I subspace at ǫ = 0. Expanding the phase space variables around ǫ = 0, we have: where φ, p φ ,h ij ,π ij describe a Bianchi I background geometry, and δφ (n) ( x), δp ij ( x), δπ ij (n) ( x) describe nth-order perturbations thereon. The homogeneous variables satisfy the Poisson brackets As it is common in the literature, we restrict ourselves to diagonal Bianchi I metrics, such that the phase space variablesh ij andπ ij take a diagonal form in an appropriate system of coordinates x i h ij = diag(a 2 1 , a 2 2 , a 2 3 ) ,π ij = diag π a 1 2 a 1 , π a 2 2 a 2 , π a 3 2 a 3 , where a i define the three directional scale factors; it follows from (2.7) that a i and π a j are canonically conjugate, {a i , π a j } = 1 V 0 δ ij (note that the subscripts i, j in a i and π a j are just labels, and not tensorial indices). From now on we will raise and lower all spatial indices i, j, k, ... withh ij and its inverse. The next step is to expand the constraints (2.3) and (2.4) in perturbations where the superscripts in parenthesis denote the order in our perturbative expansion. As mentioned before, we want to focus on linear perturbations. This will require, on the one hand, to keep only first order perturbations δh φ ( x)-since these will be the only perturbations in the rest of this paper, from now on we will remove the label (1)-and, on the other hand, to truncate all constraints at quadratic order in these fields.
In addition to the perturbative expansion of the constraints, we also expand the lapse function as N + δN ( x), and the shift as N i + δN i ( x), where, in the following, N is a homogeneous function, and we take N i = 0, so the background line element takes the familiar form (2.10) In the rest of this section we discuss the background sector, and leave the study of the inhomogeneous degrees of freedom for the next section. Because of homogeneity, V (0) i identically vanish (note that they are proportional to derivatives in space-like directions). Hence, the homogeneous degrees of freedom are subject to only one constraint, S (0) , which takes the form − a 1 π a 1 a 2 π a 2 − a 2 π a 2 a 3 π a 3 − a 3 π a 3 a 1 π a 1 +p 2 φ + 2hV (φ) ≈ 0, (2.11) withh = (a 1 a 2 a 3 ) 2 = a 6 the determinant ofh ij , and we have defined a ≡ (a 1 a 2 a 3 ) 1/3 as the mean scale factor. Then, the Hamiltonian (2.2) reduces to If we choose N = 1, H BI generates evolution in standard cosmic time t. Hamilton's equations of motion then readȧ These equations fully determine the dynamical evolution of the Bianchi I geometry, once suitable initial data satisfying the scalar constraint is provided. For convenience in the interpretation of the solutions, it is useful to introduce (see e.g. Ref. [8] for details), on the one hand, the average Hubble rate H =ȧ a . Its relation to the directional Hubble rates H i ≡ȧ i a i is H = 1 3 (H 1 + H 2 + H 3 ). On the other hand, the anisotropic shear σ ij is defined fromπ ij bẙ where π a is the conjugate momenta of a (its can be written in terms ofȧ as π a = − 6 κ aȧ). Equation (2.14) is equivalent to saying that the components of σ ij are related to the canonical variables π a i by Using the equations of motion (2.13), one can check that σ i = H i − H, and from this it is obvious that σ i are not all independent, but satisfy σ 1 + σ 2 + σ 3 = 0 -or in other words, σ ij is traceless. It is also convenient to define the shear squared With these definitions, the equations of motion (2.13) can be written in the more familiar form 18) and the scalar constraint (2.11) as

B. Implementation in Mathematica
We will begin here the description of the main steps carried out in the Mathematica notebook [17].

Preliminaries
The perturbative expansions are carried out employing the package xPert [13]. Hence, the first step is to load this package with the following command

In[1] := <<xAct'xPert'
We define the three-dimensional manifold M3 with abstract indices {i, j, k, l, m, n}: The spatial slices will be parameterized by a time variable t. We define it with the command

In[3] := DefParameter[t,PrintAs->"t"];
We also define the gravitational coupling constant: The (Riemannian) spatial metric h and its covariant derivative CD are introduced by means of Note that we have allowed the spatial metric h to depend on the time parameter t. We now introduce perturbations of the metric: We define the scalar field: and its potential: We incorporate perturbations of the scalar field with We define next canonical momenta. The momentum conjugate to the spatial metric is introduced as: In [10]

Scalar and vector constraints
We start introducing the diffeomorphism constraints defined in (2.4): In [14] As we will see below, only the linear term in the perturbative expansion of these constraints, called V (1) i ( x) above, will be relevant in the description of linearized perturbations (recall also that V (0) i are identically zero). They are defined in the notebook by: In this expression, we have imposed that the partial spatial derivatives of the background degrees of freedom vanish because of homogeneity. Let us focus now on the scalar constraint. We are going to compute each of its terms, written in (2.3), separately. First of all, the three-dimensional Ricci curvature is We now expand this term in perturbations by using In [17] where we have imposed again homogeneity of the background metric, ∂ i h jk = 0, and we have put to zero the second order perturbations, δh ij = 0. On the other hand, the first term in (2.3) (the "kinetic" term of the gravitational sector) is In [18] where we have imposed δh (2) ij = 0 and δπ (2) ij = 0 in the last line. The terms in (2.3) that depend on the scalar field are In [19]  In addition, the conjugate variable π a to the average scale factor a is defined as Expression (2.14) above is implemented as With this, we express the first-order diffeomorphism constraints in terms of the shear as and in terms of shear: In[30] := S1b = S1a/.bgmomrule//ToCanonical; The part of the scalar constraint that is quadratic in perturbations will be discussed in Sec. V.

III. SCALAR-VECTOR-TENSOR DECOMPOSITION
A. Summary of the theory Linear perturbations satisfy, via (2.1) and (2.7), the canonical Poisson brackets For convenience, we Fourier expand these perturbations and their conjugate momenta 2 where k · x = k i x i and such that k i is time independent (the comoving wavevector). The Poisson brackets (3.1) become We now perform a generalization of the scalar-vector-tensor decomposition of δh ij ( k) and δπ ij ( k) that is commonly used in FLRW spacetimes. Although this decomposition is adapted to the rotational invariance of FLRW geometries, it will also be useful in Bianchi I, since it will allow us to work with variables that become the familiar scalar, vector, and tensor modes when the background geometry isotropizes (as it quickly happens if there is a phase of inflation). We define now a basis of 3 × 3 symmetric matrices as Here,k is the unit vector (with respect toh ij ) in the direction of k, andê 1 ,ê 2 are two unit vectors orthogonal among themselves and tok. 3 We now define γ n ( k) and π n ( k) as the components of δh ij ( k) and δπ ij ( k), respectively, in this basis In FLRW spacetimes γ n and π n are called scalar modes for n = 1, 2, vector modes for n = 3, 4, and tensor modes for n = 5, 6, due to their properties under rotations around the directionk. We will keep using these names along this paper. The non-zero Poisson brackets of these modes are Furthermore, we define and we denote all the degrees of freedom in perturbations as γ α ( k) and π α ( k) with α = 0, · · · , 6. It will be useful in the next section to define the products of the shear tensor σ ij and A ij (n) as for n = 2, · · · , 6 (σ (1) vanishes, because it is proportional to the trace of σ ij ). Note however that σ (n) (k) is not the Fourier transform of any of the components of σ ij .
We can now write the linear constraints S (1) ( x) and V (1) i ( x) in terms of the variables γ α ( k) and π α ( k). For this purpose, we first expand the constraints in Fourier modes S (1) ( x) = kS (1) ( k)e i k· x and V (1) (1) i ( k)e i k· x , and then we replace (3.6). Explicit expressions are provided in Appendix B of Ref. [8].

B. Implementation in Mathematica
We begin by defining the vectors k,ê 1 , andê 2 as follows We implement the traces and orthogonality properties of these matrices as automatic rules using the command AutomaticRules. However, for convenience, we express the matrices A and similarly for the matrices corresponding to vector and tensor modes. Finally we define γ 0 ( k) and π 0 ( k): Below we use this function to verify that the Poisson algebra of linear constraints is closed. It is important to keep in mind that the conjugate variable to γ α ( k) is π α (− k). We will take this into account, although we will not implement it explicitly in the notebook for the sake of simplicity. We now define the components σ (n) (k) of the shear tensor following Eq. (3.9) In Consequently, these variables are not gauge invariant. 4

A. Theory
We have seven degrees of freedom (per Fourier mode) in configuration variables in the perturbations, γ α ( k). They are subject to four first class constraints,S (1) ( k) ≈ 0 andṼ (1) i ( k) ≈ 0, that are the generators of gauge transformations for each mode. Each constraint reduces the number of independent configuration variables by one. In consequence, we are left with 7 − 4 = 3 physical configuration fields, together with their conjugate momenta. We will isolate these degrees of freedom by identifying gauge invariant fields. In the Lagrangian formalism, this was accomplished in Ref. [9]. In the Hamiltonian framework, the gauge invariant variables are defined as fields that are left invariant by the gauge flow generated by the linear constraints, or equivalently, that Poisson-commute with them.
A conceptually simple and elegant procedure to find gauge invariant variables can be obtained by using the ideas of [19], which were in part applied to FLRW cosmologies in [3]. The details of this procedure in Bianchi I cosmologies can be found in [8]. In summary, the idea is to find a canonical transformation from γ α ( k), π α ( k) to new variables Γ α ( k), Π α ( k) such that the new momenta Π α ( k) for α = 3, 4, 5, 6 are proportional to the four constraintsS (1) ( k) andṼ (1) i ( k), respectively. More concretely, we demand (1) (1) j ( k) , (4.1) where the factors 1/| k|, with | k| ≡ k i k i = 0, have been introduced for dimensional reasons, and the imaginary unit for convenience. This choice is possible because the constraints are first class, and it can be done globally in the perturbed phase space because of the linearity of the system. If equation (4.1) is satisfied, then the canonical commutation relations guarantee that Γ α ( k) and Π α ( k) for α = 0, 1, 2 Poisson-commute with the constraints, and hence they are gauge invariant. This procedure also guarantees that gauge invariant fields and pure gauge ones are decoupled in the Hamiltonian (as we will see explicitly below), and hence dynamics does not mix them. One can then consistently focus attention on gauge invariant perturbations. Interestingly, the task of finding such a canonical transformation reduces to solving a Hamilton-Jacobi-like equation for a generating function, and furthermore, by working in Fourier space, this equation reduces to algebraic equations that are easy to solve in Mathematica. More concretely, we start with a generating function in Fourier space that we choose to be of type 2-i.e. it depends on old variables γ α and new momenta Π α -and from which the rest of variables are given by where B αβ and A αβ are matrices whose components depend on background variables, but not on perturbations, and furthermore A αβ is symmetric. Equations (4.3) provide then a set of algebraic relations for the 77 unknown coefficients B αβ and A αβ , although only 38 of them are independent. Hence, there is freedom in choosing gauge invariant variables. As mentioned above, we want to choose gauge invariant fields that in the isotropic limit reduce to the familiar comoving curvature perturbations and the two tensor modes. Indeed, gauge invariant fields Γ α ( k) satisfying this property can be identified by inspection, just by looking at the Poisson brackets of γ α and the linear constraints. They are This will be our choice of gauge invariant fields. Other choices are possible, and all of them can be derived by using the companion notebook [17]. In the isotropic limit σ (n) → 0, Γ 1 and Γ 2 reduce to the familiar two polarizations of transverse and traceless tensor modes, and Γ 0 becomes proportional to the comoving curvature perturbation Identifying the form of the new gauge invariant variables Γ α in terms of γ α will simplify the computation of G( k) (for instance, this determines the value of some of the coefficients B αβ ), but it is important to emphasize that the method, and therefore its implementation in the algorithm reported in our manuscript, allows us to work with gauge invariant variables and pure gauge ones in a systematic way, without having to identify them a priori. To our knowledge, numerical tools meeting all these requirements (i.e. handle efficiently complicated phase space functions and keep the construction as general as possible) are not publicly available, at least in the context of cosmological perturbation theory or similar settings.
In addition, we will demand (see next section) the Hamiltonian to be "diagonal" in the new configuration variables and momenta, i.e. we will eliminate "cross terms" of the type Γ α Π β . This aesthetic condition will impose further restriction in the coefficients A αβ 's. The rest of free coefficients can be equated to zero for simplicity. Once all the coefficients A αβ and B αβ are specified, the form of the conjugate momenta Π α for α = 0, 1, 2 are obtained from (4.3).

B. Implementation in Mathematica
We summarize here the main steps of the code [17]. We start defining the new variables Γ α ( k), Π α ( k): and similarly for α = 1, . . . , 6. As explained above, rather than solving all the equations that constrain the coefficients A αβ and B αβ , we simplify the calculation by identifying suitable gauge invariant variables (see equation (4.6) above). These variables are named in the code as Γ0new, Γ1new , Γ2new: We now follow the strategy described above. Namely, we first obtain an expression for the old momenta π α in terms of γ α and Π β by taking derivative of the generating function with respect to γ α . These expressions can then be substituted in the scalar and vector constraints, to express them in terms of old configuration variables γ α and new momenta Π β . By equating these constraints to Π α for α = 3, 4, 5, 6 as indicated in (4.1), we obtain 44 algebraic equations for the coefficients A αβ , out of which 38 are independent. Some of the remaining coefficients can be solved by demanding that the new Hamiltonian does not contain cross terms between new momenta and new configuration variables. Finally, the remaining coefficients are set to zero for the sake of simplicity. See the notebook [17] for details.

A. Theory
The dynamics of perturbations is generated by the second order scalar constraint d 3 x N S (2) ( x). Let us notice that the second order vector constraints V (2) i do not contribute since the homogeneous part of the shift N i is zero and the next contribution would come from δN i ( x) V (2) i , which is third order in perturbations. If we expand the fields inside S (2) ( x) in Fourier modes, the expression for d 3 x N S (2) ( x) can be written as k NS (2) ( k), whereS (2) ( k) is quadratic in perturbations. In each term, one perturbation is evaluated at k and the other at − k.
From the expression forS (2) ( k) in terms of δh ij and δπ ij , we obtain a Hamiltonian for the new variables Γ α ( k) and Π α ( k), by first implementing the change from δh ij ( k), δπ ij ( k) to γ( k), π( k) and then from the later to Γ α ( k), Π α ( k). One has to keep in mind that these transformations involve coefficients that depend on functions of the Bianchi I background geometry, and therefore they have to be understood as time-dependent quantities. This means that the final Hamiltonian is equal to the original one in new variables, plus the time derivative of the generating function of the canonical transformation, where the time derivative only affects the background functions. The first canonical transformation can be implemented by the following generating function ij ( k) depend on time. The second canonical transformation from γ α ( k), π α ( k) to Γ α ( k), Π α ( k) is defined by the generating function (4.2).
After implementing these transformations, one can check that gauge invariant fields decouple from pure gauge ones, and one obtains the following Hamiltonian for the former (see [8] for further details) 2) where δ µ,µ ′ is the Kronecker delta, and k 2 (t) ≡ a 2 (t) k i k j = a 2 (t) The (time-dependent) effective potentials U µµ ′ are symmetric in µ and µ ′ and they become diagonal (i.e. proportional to δ µµ ′ ) in the isotropic limit. But in presence of anisotropies, they couple gauge invariant perturbations among themselves. They are explicitly written in Appendix A.
The equations of motion are (we use cosmic time) Combining these equations into second order differential equations we obtain This is a set of three coupled, second order, ordinary differential equations for each wavevector k, and they reduce to the familiar (decoupled) equations for scalar and tensor perturbations in the isotropic FLRW limit

B. Implementation in Mathematica
We start with the expression k NS (2) ( k) in terms of δh ij ( k) and δπ ij ( k). As mentioned above, each termS (2) ( k) is quadratic in perturbations, with one field evaluated at k and the other at − k. In our Mathematica notebook this will remain implicit, since the code will be significantly simpler in this way.
We first obtain an expression forS (2) ( k) in terms of δh ij ( k) and δπ ij ( k): where, as before, we have replaced spatial derivatives by i k. The second step is to move from δφ( k), δp φ ( k), δh ij ( k), δπ ij ( k) to γ α ( k), π β ( k) : In We have verified that there is no coupling between gauge invariant and pure gauge variables, once the background constraint is also imposed. Therefore, we can focus on the part of the Hamiltonian that contains only gauge-invariant fields. Nevertheless, this part still contains cross terms between Π µ and Γ µ , with µ = 0, 1, 2. Fortunately, and as discussed above, the generating function G( k) still has some free parameters that can be fixed by requiring that these terms vanish. Concretely, we require the coefficients multiplying the cross terms ( Finally, in order to focus on gauge invariant fields, we set pure gauge variables to zero in the code. This yields an expression for the Hamiltonian that we denote by H2f, which is further simplified by means of the rule

VI. DISCUSSION
We have described in this paper the main steps of a computer code written in the symbolic language of Mathematica, and made publicly available in [17], that derives gauge invariant linear perturbations in Bianchi type I cosmological spacetimes in a Hamiltonian or phase space approach. In this formulation, gauge invariant linear perturbations are defined by fields that Poisson-commute with the linear constraints of the theory, and they can be found systematically by solving a canonical transformation that identifies some of the new momenta with the constraints. We have described in detail the implementation of this procedure in Mathematica. Our code provides an efficient tool to explore different choices of linear gauge invariant fields, and to derive the equations of motion they satisfy. It can also be used to work with gauge dependent variables, after choosing a gauge, and to relate physical observables written in different gauges. Furthermore, we have complemented this analysis with a computer code, based on the C programing language, available in [18], that solves the equations of motions and computes observables that can be compared with current and future data from the cosmic microwave background (CMB). Our computer codes should be of great utility for researchers interested in cosmological perturbations in Bianchi spacetimes and their consequences for the CMB, as well as in isotropic FLRW, which is obtained by simply putting the anisotropies to zero. Our theoretical and numerical analysis can also be of interest for pedagogical purposes, since it provides a step by step, guided way of implementing the theory of gauge invariant cosmological perturbations in a computer.