Next Article in Journal
A-Statistical Convergence Properties of Kantorovich Type λ-Bernstein Operators Via (p, q)-Calculus
Next Article in Special Issue
From Fractional Quantum Mechanics to Quantum Cosmology: An Overture
Previous Article in Journal
Existence of Bounded Solutions to a Modified Version of the Bagley–Torvik Equation
Previous Article in Special Issue
The Effect of a Positive Cosmological Constant on the Bounce of Loop Quantum Cosmology
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Physics and Astronomy, Louisiana State University, Baton Rouge, LA 70803, USA
2
Department of Physics, National Institute of Technology Karnataka, Surathkal, Mangaluru 575025, India
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(2), 290; https://doi.org/10.3390/math8020290
Submission received: 31 December 2019 / Revised: 4 February 2020 / Accepted: 15 February 2020 / Published: 20 February 2020
(This article belongs to the Special Issue Mathematical and Computational Cosmology)

Abstract

:
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ître–Robertson–Walker spacetimes is straightforward.

1. 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 [4,5,6,7] for a recent discussion), and it was extended to anisotropic Bianchi I spacetimes in [8] (see also [9,10] for a previous analysis starting from Einstein’s 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). For the application to perturbations on spherically symmetric spacetimes see also [14,15]. The merit of our algorithm is that we combine a series of analytical and computational techniques to make the problem tractable. We find that the phase space formulation enormously helps 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 [16], 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 [17], 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 into small steps that are, one by one, introduced theoretically and implemented computationally. Namely, Section 2 introduces Bianchi I spacetimes and linear perturbations thereon in the Arnowitt–Deser–Misner (ADM) formalism [2]. Section 3 introduces the generalization of the scalar-vector-tensor decomposition commonly used in FLRW backgrounds. Section 4 introduces gauge invariant fields, and in Section 5, we study their dynamics. Section 6 contains a short summary and some concluding remarks.

2. Perturbed Bianchi I Spacetimes in the ADM Formalism

2.1. 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 the R 3 topology. In the following, we will restrict ourselves to a finite volume V 0 relative to an auxiliary flat Euclidean metric δ i j defined in M 3 (This is equivalent to putting the universe in a box of arbitrarily large, but finite volume V 0 , with periodic boundary conditions. We do this for convenience in the expressions below. The volume V 0 will not impact predictions, and it can be taken to infinity at the end of the calculations.). 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 i j ( x ) , π i j ( x ) ) defined in M 3 , where Latin indices i , j run from one to three. Here, h i j ( x ) is a Riemannian metric that describes the intrinsic spatial geometry of M 3 and π i j ( x ) its conjugate momentum. The non-vanishing Poisson brackets between these fields are:
{ Φ ( x ) , P Φ ( x ) } = δ ( 3 ) ( x x ) , { h i j ( x ) , π k l ( x ) } = δ ( i k δ j ) l δ ( 3 ) ( x x ) .
where δ ( i k δ j ) l 1 2 ( δ i k δ j l + δ j k δ i l ) . Dynamics is generated by the Hamiltonian:
H = d 3 x N ( x ) S ( x ) + N i ( x ) V i ( x ) ,
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:
S ( x ) = 2 κ h π i j π i j 1 2 π 2 h 2 κ ( 3 ) R + 1 2 h P Φ 2 + h V ( Φ ) + h 2 D i Φ D i Φ 0 ,
V i ( x ) = 2 h h i j D k ( h 1 / 2 π k j ) + P Φ D i Φ 0 ,
where κ = 8 π G , and D i is the covariant derivative compatible with the spatial metric h i j , h its determinant, and ( 3 ) R its Ricci scalar curvature. Given N ( x ) , N i ( x ) , and a solution to Hamilton’s equations, h i j ( x , t ) , the spacetime metric takes the form:
d s 2 = ( N 2 N i N i ) d t 2 + 2 N i d x i d t + h i j d x i d x j ,
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:
h i j ( x , ϵ ) = h ˚ i j + ϵ δ h i j ( 1 ) ( x ) + + ϵ n n ! δ h i j ( n ) ( x ) + , π i j ( x , ϵ ) = π ˚ i j + ϵ δ π i j ( 1 ) ( x ) + + ϵ n n ! δ π i j ( n ) ( x ) + , Φ ( x , ϵ ) = ϕ + ϵ δ ϕ ( 1 ) ( x ) + + ϵ n n ! δ ϕ ( n ) ( x ) + , P Φ ( x , ϵ ) = p ϕ + ϵ δ p ϕ ( 1 ) ( x ) + + ϵ n n ! δ p ϕ ( n ) ( x ) + ,
where ϕ , p ϕ , h ˚ i j , π ˚ i j describe a Bianchi I background geometry, and δ ϕ ( n ) ( x ) , δ p ϕ ( n ) ( x ) , δ h i j ( n ) ( x ) , δ π i j ( n ) ( x ) describe the n th -order perturbations thereon. The homogeneous variables satisfy the Poisson brackets:
{ ϕ , p ϕ } = 1 V 0 , { h ˚ i j , π ˚ k l } = 1 V 0 δ ( i k δ j ) l .
As is common in the literature, we restrict ourselves to diagonal Bianchi I metrics, such that the phase space variables h ˚ i j and π ˚ i j take a diagonal form in an appropriate system of coordinates x i :
h ˚ i j = diag ( a 1 2 , a 2 2 , a 3 2 ) , π ˚ i j = 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 (7) that a i and π a j are canonically conjugate, { a i , π a j } = 1 V 0 δ i j (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 , with h ˚ i j and its inverse.
The next step is to expand the constraints (3) and (4) in perturbations:
S ( x ) = S ( 0 ) + S ( 1 ) ( x ) + S ( 2 ) ( x ) + S ( 3 ) ( x ) + , V i ( x ) = V i ( 0 ) + V i ( 1 ) ( x ) + V i ( 2 ) ( x ) + V i ( 3 ) ( x ) + ,
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 i j ( 1 ) ( x ) , δ π i j ( 1 ) ( x ) , δ ϕ ( 1 ) ( x ) , and δ p ϕ ( 1 ) ( 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:
d s 2 = N 2 d t 2 + h ˚ i j d x i d x j .
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 i ( 0 ) 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:
S ( 0 ) = 1 2 h ˚ [ κ a 1 2 π a 1 2 2 + a 2 2 π a 2 2 2 + a 3 2 π a 3 2 2 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 + 2 h ˚ V ( ϕ ¯ ) ] 0 ,
with h ˚ = ( a 1 a 2 a 3 ) 2 = a 6 the determinant of h ˚ i j , and we have defined a ( a 1 a 2 a 3 ) 1 / 3 as the mean scale factor. Then, the Hamiltonian (2) reduces to:
H BI = M 3 d 3 x N S ( 0 ) = V 0 N S ( 0 ) .
If we choose N = 1 , H BI generates evolution in standard cosmic time t. Hamilton’s equations of motion then read:
a ˙ i = { a i , H BI } , π ˙ a i = { π a i , H BI } , ϕ ˙ = { ϕ , H BI } , p ϕ ˙ = { p ϕ , H BI } .
These equations fully determine the dynamical evolution of the Bianchi I geometry, once suitable initial data satisfying the scalar constraint are provided. For convenience in the interpretation of the solutions, it is useful to introduce (see, e.g., [8] for details), on the one hand, the average Hubble rate H = a ˙ a . Its relation to the directional Hubble rates H i a ˙ i a i is H = 1 3 H 1 + H 2 + H 3 . On the other hand, the anisotropic shear σ i j is defined from π ˚ i j by:
π ˚ i j = h 1 / 6 6 π a h ˚ i j + h 1 / 2 2 κ σ i j ,
where π a is the conjugate momenta of a (it can be written in terms of a ˙ as π a = 6 κ a a ˙ ). Equation (14) is equivalent to saying that the components of σ i j :
σ i j = diag ( a 1 2 σ 1 , a 2 2 σ 2 , a 3 2 σ 3 ) ,
are related to the canonical variables π a i by:
π a i = 1 κ a 3 a i ( σ i 2 H ) .
Using the equations of motion (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, σ i j is traceless. It is also convenient to define the shear squared:
σ 2 = σ i j σ i j = σ 1 2 + σ 2 2 + σ 3 2 = ( H 1 H ) 2 + ( H 2 H ) 2 + ( H 3 H ) 2 .
With these definitions, the equations of motion (13) can be written in the more familiar form:
a ¨ a = κ 6 [ ρ + 3 P ] σ 2 3 , ϕ ¨ + 3 a ˙ a ϕ ˙ + d V ( ϕ ) d ϕ = 0 ,
and the scalar constraint (11) as:
H 2 = κ 3 ρ + σ 2 6 .
Here, ρ 1 2 ϕ ˙ 2 + V ( ϕ ) and P 1 2 ϕ ˙ 2 V ( ϕ ) are the energy and pressure densities of ϕ , respectively. The expressions above are equivalent to the diagonal components of Einstein’s equations. For  σ 2 = 0 , they reduce to the Friedmann–Lemaître theory of isotropic cosmologies. On the other hand, the equations of motion for the shear are:
σ ˙ j i = 3 H σ j i ,
whose solutions are simply σ i = Σ i / a 3 , where Σ i are constants constrained by Σ 1 + Σ 2 + Σ 3 = 0 . This implies σ 2 = Σ 2 a 6 , with Σ 2 Σ 1 2 + Σ 2 2 + Σ 3 2 .

2.2. Implementation in Mathematica

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

2.2.1. 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 } :
In[2] := DefManifold[M3,3, { 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:
In[4] := DefConstantSymbol[ κ ];
The (Riemannian) spatial metric h and its covariant derivative CD are introduced by means of:
In[5] := DefMetric[1,h[- i ,- j ],CD,{";"," D "}, Otherdependencies-> { t } ,WeightedWithBasis->AIndex];
Note that we have allowed the spatial metric h to depend on the time parameter t . We now introduce perturbations of the metric:
In[6] := DefMetricPerturbation[h, δ h, ϵ ];
We define the scalar field:
In[7] := DefTensor[ ϕ [],{M3,t}];
and its potential:
In[8] := DefScalarFunction[V];
We incorporate perturbations of the scalar field with:
In[9] := DefTensorPerturbation[ δ ϕ [LI[1]], ϕ [],{M3,t}];
We define next canonical momenta. The momentum conjugate to the spatial metric is introduced as:
In[10] := DefTensor[P[i,j],{M3,t}];
and its perturbations are defined by:
In[11] := DefTensorPerturbation[ δ P[LI[1],i,j],P[i,j],{M3,t}];
Finally, we define the momentum conjugate to the scalar field:
In[12] := DefTensor[P ϕ [],{M3,t}];
and its corresponding perturbation:
In[13] := DefTensorPerturbation[ δ P ϕ [LI[1]],P ϕ [],{M3,t}];

2.2.2. Scalar and Vector Constraints

We start introducing the diffeomorphism constraints defined in (4):
In[14] := diffeo=-2PD[-k]@(h[-i,-j]P[j,k])+P[k,j]PD[-i]@h[-k,-j]+P ϕ []PD[-i]@ ϕ [];
As we will see below, only the linear term in the perturbative expansion of these constraints, called V i ( 1 ) ( x ) above, will be relevant in the description of linearized perturbations (recall also that V i ( 0 ) are identically zero). They are defined in the notebook by:
In[15] := (Perturbed[diffeo,1]/ ϵ );
   diffeo1 = %/.MakeRule[PD[-i]@ ϕ [],0]/.MakeRule[PD[-i]@h[-j,-k],0]
   /.MakeRule[PD[-i]@P[j,k],0];
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 (3), separately. First of all, the three-dimensional Ricci curvature is:
In[16] := ricci=(Deth[])] (1/2)h[j,k]RiemannCD[-j,-i,-k,i]//RiemannToChristoffel
   //ChristoffelToMetric//Simplification//NoScalar;
We now expand this term in perturbations by using:
In[17] := Perturbed[ricci,2]/.MakeRule[{PD[-i]@h[-j,-k],0}]//ExpandPerturbation;
   r2 = %/.MakeRule[h[LI[2],-i,-j],0];
where we have imposed again homogeneity of the background metric, i h j k = 0 , and we have put to zero the second order perturbations, δ h i j ( 2 ) = 0 .
On the other hand, the first term in (3) (the “kinetic” term of the gravitational sector) is:
In[18] := (Deth[]) (-1/2)P[i,j]P[k,l](h[-i,-k]h[-j,-l]-h[-i,-j]h[-k,-l]/2);
   ExpandPerturbation[Perturbed[%,2]];
   pipi = %/.MakeRule[{ δ P[LI[2],-i,-j],0}]/.MakeRule[{ δ h[LI[2],-i,-j],0}];
where we have imposed δ h i j ( 2 ) = 0 and δ π i j ( 2 ) = 0 in the last line.
The terms in (3) that depend on the scalar field are:
In[19] := 1/2Deth[] (-1/2)P ϕ [] 2+Deth[] (1/2)(1/2PD[-i]@ ϕ []PD[-j]@ ϕ []h[i,j]
   +V[ ϕ []]);
   ExpandPerturbation[Perturbed[%,2]];
   matter = %/.MakeRule[{ δ h[LI[2],-i,-j],0}]/.MakeRule[{ δ ϕ [LI[2]],0}]
   /.MakeRule[{ δ P ϕ [LI[2]],0}]/.MakeRule[{PD[-i]@ ϕ [],0}];
Putting everything together, the scalar constraint, up to second order in perturbations, is:
In[20] := S = Series[(2 κ )pipi-1/(2 κ )r2+matter,{ ϵ ,0,2}];
The scalar constraint contributes with 0th, 1st and 2nd order terms in the perturbative expansion. Let us identify each one. We first introduce the shear tensor as:
In[21] := DefTensor[ σ [i,j],{M3,t},Symmetric[{i,j}]];
and the shear σ σ 2
In[22] := DefTensor[ σ b[],{M3,t}];
We impose that the shear is traceless and symmetric, and its relation with σ b[] 2, with the following automatic rules:
In[23] := AutomaticRules[ σ ,MakeRule[{ σ [i,j]h[-i,-j],0}]];
   AutomaticRules[ σ ,MakeRule[{ σ [i,-i],0}]];
   AutomaticRules[ σ ,MakeRule[{ σ [i,j] σ [k,l]h[-i,-k]h[-j,-l], σ b[] 2}]];
In addition, the conjugate variable π a to the average scale factor a is defined as:
In[24] := DefTensor[ π a[],{M3,t}];
Expression (14) above is implemented as:
In[25] := bgmomrule=MakeRule[{P[i,j], π a[]/6 Deth[] (1/6) h[i,j]+Deth[] (1/2)/(2 κ )
    σ [i,j]}];
With this, we express the first-order diffeomorphism constraints in terms of the shear as:
In[26] := diffeo1/.bgmomrule//org//ChristoffelToMetric//Simplification//NoScalar;
   diffeoa=%/.MakeRule[{PD[-i]@h[-j,-k],0}];
Similarly, the zeroth order scalar constraint is:
In[27] := S0 = SeriesCoefficient[S,0];
We write it in terms of σ 2 and π a by:
In[28] := S0a=S0/.bgmomrule//ToCanonical;
In a similar way, we define the part of the scalar constraint that is linear in perturbations as:
In[29] := S1a = SeriesCoefficient[S,1];
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 Section 5.

3. Scalar-Vector-Tensor Decomposition

3.1. Summary of the Theory

Linear perturbations satisfy, via (1) and (7), the canonical Poisson brackets:
{ δ ϕ ( x ) , δ p ϕ ( x ) } = δ ( 3 ) ( x x ) 1 V 0 ; { δ h i j ( x ) , δ π k l ( x ) } = δ ( i k δ j ) l δ ( 3 ) ( x x ) 1 V 0 .
For convenience, we Fourier expand these perturbations and their conjugate momenta (The Fourier expansion of fields is adapted to the fiducial cell of volume V 0 , so the wavenumbers k will take values on a discrete lattice k 2 π / ( V 0 ) 1 / 3 Z 3 . In the limit V 0 , one recovers k R 3 ).
δ ϕ ( x ) = k 0 δ ϕ ˜ ( k ) e i k · x ; δ p ϕ ( x ) = k 0 δ p ˜ ϕ ( k ) e i k · x ,
δ h i j ( x ) = k 0 δ h ˜ i j ( k ) e i k · x ; δ π i j ( x ) = k 0 δ π ˜ i j ( k ) e i k · x ,
where k · x = k i x i and such that k i is time independent (the comoving wavevector).
The Poisson brackets (21) become:
{ δ ϕ ˜ ( k ) , δ p ˜ ϕ ( k ) } = V 0 1 δ k , k ; { δ h ˜ i j ( k ) , δ π ˜ k l ( k ) } = V 0 1 δ ( i k δ j ) l δ k , k .
We now perform a generalization of the scalar-vector-tensor decomposition of δ h ˜ i j ( k ) and δ π ˜ i j ( 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:
A i j ( 1 ) = h ˚ i j 3 , A i j ( 4 ) = 1 2 k ^ i e ^ 2 j + k ^ j e ^ 2 i , A i j ( 2 ) = 3 2 k ^ i k ^ j h ˚ i j 3 , A i j ( 5 ) = 1 2 e ^ 1 i e ^ 1 j e ^ 2 i e ^ 2 j , A i j ( 3 ) = 1 2 k ^ i e ^ 1 j + k ^ j e ^ 1 i , A i j ( 6 ) = 1 2 e ^ 1 i e ^ 2 j + e ^ 1 j e ^ 2 i .
Here, k ^ is the unit vector (with respect to h ˚ i j ) in the direction of k , and e ^ 1 , e ^ 2 are two unit vectors orthogonal among themselves and to k ^ (Note that the three unit vectors k ^ , e ^ 1 , e ^ 2 are time dependent. This is because, on the one hand, the norm of k i is time dependent and, on the other hand, the unit vectors e ^ 1 , e ^ 2 need to rotate in time to remain orthogonal to k ^ , unless k ^ points in one of the principal directions. For the details on how to compute the time dependence of the unit vectors k ^ , e ^ 1 , and e ^ 2 , see Appendix A in [8]). We now define γ n ( k ) and π n ( k ) as the components of δ h ˜ i j ( k ) and δ π ˜ i j ( k ) , respectively, in this basis:
δ h ˜ i j ( k ) = n = 1 6 γ n ( k ) A i j ( n ) ( k ^ ) ; δ π ˜ i j ( k ) = n = 1 6 π n ( k ) A ( n ) i j ( k ^ ) .
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 direction k ^ . We will keep using these names throughout this paper. The non-zero Poisson brackets of these modes are:
{ γ n ( k ) , π m ( k ) } = V 0 1 δ n m δ k , k .
Furthermore, we define:
γ 0 4 κ δ ϕ ˜ ( k ) , π 0 1 / 4 κ δ p ˜ ϕ ( k ) ,
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 σ i j and A ( n ) i j as:
σ ( n ) ( k ^ ) σ i j A ( n ) i j ( k ^ ) ,
for n = 2 , , 6 ( σ ( 1 ) vanishes, because it is proportional to the trace of σ i j ). Note however that σ ( n ) ( k ^ ) is not the Fourier transform of any of the components of σ i j .
We can now write the linear constraints S ( 1 ) ( x ) and V i ( 1 ) ( x ) in terms of the variables γ α ( k ) and π α ( k ) . For this purpose, we first expand the constraints in Fourier modes S ( 1 ) ( x ) = k S ˜ ( 1 ) ( k ) e i k · x and V i ( 1 ) ( x ) = k V ˜ i ( 1 ) ( k ) e i k · x , and then, we replace (26). Explicit expressions are provided in Appendix B of [8].

3.2. Implementation in Mathematica

We begin by defining the vectors k , e ^ 1 , and e ^ 2 as follows:
In[31] := DefTensor[kv[-i],{M3,t}];
In[32] := DefTensor[e1[-i],{M3,t},OrthogonalTo->kv[i]];
In[33] := DefTensor[e2[-i],{M3,t},OrthogonalTo->{kv[i],e1[i]}];
We define the norm of k as:
In[34] := DefTensor[k[],{M3,t}];
In[35] := AutomaticRules[kv,MakeRule[{kv[-i]kv[-b]h[i,j],k[] 2}]];
and we add automatic rules to indicate that e ^ 1 and e ^ 2 are unit vectors:
In[36] := AutomaticRules[e1,MakeRule[{e1[-i]e1[-b]h[i,j],1}]];
In[37] := AutomaticRules[e2,MakeRule[{e2[-i]e2[-b]h[i,j],1}]];
The scalar, vector, and tensor modes are defined as follows. First, the symmetric matrix A i j ( 1 ) is introduced as:
In[38] := DefTensor[A1[-i,-j],{M3,t},Symmetric[{-i,-j}]];
and γ 1 ( k ) and π 1 ( k ) as:
In[39] := DefTensor[ γ 1[],{M3,t}];
In[40] := DefTensor[ π 1[],{M3,t}];
In the same way, we define all other tensors A2[-i,-j], , A6[-i,-j] and the modes γ 2[], , γ 6[], and π 2[], , π 6[].
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 i j ( n ) in terms of the background metric h and the orthogonal vectors kv[-i], e1[-i], and e2[-i] with the command MakeRule. For instance, for the matrices associated with scalar modes, we define:
In[41] := A1rule=MakeRule[{A1[-i,-j],h[-i,-j]/Sqrt[3]}];
In[42] := A2rule=MakeRule[{A2[-i,-j],Sqrt[3/2](kv[-i]kv[-j]/k[] 2-h[-i,-j]/3)}];
and similarly for the matrices corresponding to vector and tensor modes. Finally, we define γ 0 ( k ) and π 0 ( k ) :
In[43] := DefTensor[ γ 0[],{M3,t}];
In[44] := DefTensor[ π 0[],{M3,t}];
and their relation with the perturbations of the scalar field and its momentum:
In[45] := γ 0rule=MakeRule[{ δ ϕ [LI[1]], γ 0[]/Sqrt[4 κ ]}];
    π 0rule=MakeRule[{ δ P ϕ [LI[1]], π 0[]Sqrt[4 κ ]}];
We can now introduce the rules that implement the decomposition (26):
In[46] := moderule1=MakeRule[{ δ h[LI[1],-i,-j], γ 1[]A1[-i,-j]+ γ 2[]A2[-i,-j]+ γ 3[]A3[-i,-j]
   + γ 4[]A4[-i,-j]+ γ 5[]A5[-i,-j]+ γ 6[]A6[-i,-j]}];
   moderule2=MakeRule[ δ P[LI[1],i,j], π 1[]A1[i,j]+ π 2[]A2[i,j]+ π 3[]A3[i,j]
   + π 4[]A4[i,j]+ π 5[]A5[i,j]+ π 6[]A6[i,j]];
Next, we implement the canonical Poisson brackets between γ α ( k ) and π α ( k ) by means of the following function:
In[47] := PoissonBracket[f_,g_,q_List,p_List]/;Length[q]==Length[p]:=D[f,{q}].D[g,{p}]
   -D[f,{p}].D[g,{q}];
   PoissonBracket[most__,q:Except[_List],p_]:=PoissonBracket[most,{q},p];
   PoissonBracket[most__,p:Except[_List]]:=PoissonBracket[most,{p}];
If we introduce the following arrays:
In[48] := Q={ γ 0[], γ 1[], γ 2[], γ 3[], γ 4[], γ 5[], γ 6[]};
   PQ={ π 0[], π 1[], π 2[], π 3[], π 4[], π 5[], π 6[]};
one can compute the Poisson brackets of any two phase space functions of perturbations. For instance:
In[49] := PoissonBracket[ γ 0[], π 0[],Q,PQ];
Out[49] := 1
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 Equation (29):
In[50] := DefTensor[ σ 2[],{M3,t}];
and similarly for σ 3[], , σ 6[]. The relation of these quantities and σ i j is implemented via the rule:
In[51] := sheardecomposition=MakeRule[{ σ [-i,-j], σ 2[]A2[-i,-j]+ σ 3[]A3[-i,-j]
   + σ 4[]A4[-i,-j]+ σ 5[]A5[-i,-j]+ σ 6[]A6[-i,-j]}];
It will be useful in the next sections to define the following rule:
In[52] := σ [-i,-j] σ [i,j]/.sheardecomposition//org;
    σ brule=MakeRule[{ σ b[],% (1/2)}];
We finish this section by writing the linear constraints in Fourier space S ˜ ( 1 ) ( k ) and V ˜ i ( 1 ) ( k ) and in terms of the new variables γ α ( k ) and π α ( k ) . For the diffeomorphism constraints V ˜ i ( 1 ) ( k ) , this is implemented by applying the rules:
In[53] := diffeoa/.MakeRule[{PD[-b]@ δ h[LI[1],-i,-k],I*kv[-j] δ h[LI[1],-i,-k]}]
   /.MakeRule[{PD[-i]@ δ P[LI[1],j,k],I*kv[-d] δ P[LI[1],j,k]}]
   /.MakeRule[{PD[-i]@ δ ϕ [LI[1]],I*kv[-i] δ ϕ [LI[1]]}]/.moderule1
   /.moderule2/. γ 0rule/. π 0rule//ContractMetric;
We further simplify the final expression:
In[54] := diffeob = %/.sheardecomposition/.A1rule/.A2rule/.A3rule/.A4rule
   /.A5rule/.A6rule//org;
Now, we decompose these constraints in their projections in the directions k ^ , e ^ 1 , and e ^ 2 :
In[55] := diffeob1=kv[i]diffeob//ToCanonical;
   diffeob2=e1[i]diffeob//ToCanonical;
   diffeob3=e2[i]diffeob//ToCanonical;
We proceed in a similar way with S ˜ ( 1 ) ( k ) . We do it in several steps. In the first one, we substitute i by i k i :
In[56] := S1b/.MakeRule[{PD[-i]@ δ h[LI[1],-j,-k],I*kv[-j] δ h[LI[1],-i,-k]}]
   /.MakeRule[{PD[-j]@ δ h[LI[1],-i,-k],I*kv[-j] δ h[LI[1],-i,-k]}]
   /.MakeRule[{PD[-i]@kv[i],0}]//ContractMetric;
We then apply the SVT decomposition by means of:
In[57] := %/.moderule1/.moderule2/. γ 0rule/. π 0rule//ToCanonical;
Finally, we simplify the result:
In[58] := S1c=%/.sheardecomposition/.A2rule/.A3rule/.A4rule/.A5rule/.A6rule
   /. σ brule//org;
One can now check that these constraints have vanishing Poisson brackets (modulo the background constraint). For instance,
In[59] := PoissonBracket[S1c,diffeob1,Q,PQ]/.bgconstraintrule//org
Out[59] := 0
where we have evaluated the background constraint on-shell. We have checked that the Poisson brackets with the remaining constraints all vanish, in these cases identically (see the Mathematica notebook [16]). Hence, they are first class constraints, as they must be from the view point of general relativity.
One can also check that none of the modes γ α or their conjugate momenta π α commute with these constraints. We show here a couple of examples:
In[60] := PoissonBracket[ γ 0[],S1c,Q,PQ]//org
Out[60] := 2 κ P ϕ h ˜ ˜
In[61] := PoissonBracket[ γ 1[],S1c,Q,PQ]//org
Out[61] := κ π a 3 h ˜ ˜ 1 / 3
Consequently, these variables are not gauge invariant (this contrasts with the isotropic case, where tensor modes γ 5 , γ 6 , and their momenta are gauge invariant).

4. Gauge Invariant Variables

4.1. 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 V ˜ i ( 1 ) ( k ) 0 , that are the generators of gauge transformations for each mode. Each constraint reduces the number of independent configuration variables by one. As a 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 [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 [18], 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 constraints S ˜ ( 1 ) ( k ) and V ˜ i ( 1 ) ( k ) , respectively. More concretely, we demand:
Π 3 ( k ) = 1 | k | S ˜ ( 1 ) ( k ) , Π 4 ( k ) = 1 i | k | k ^ j V ˜ j ( 1 ) ( k ) , Π 5 ( k ) = 1 i | k | e ^ 1 j V ˜ j ( 1 ) ( k ) , Π 6 ( k ) = 1 i | k | e ^ 2 j V ˜ j ( 1 ) ( k )
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 (30) 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, the 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:
G ( k ) = B α β ( k ) Π α ( k ) γ β ( k ) + A α β ( k ) γ α ( k ) γ β ( k ) ,
that we choose to be of type 2—i.e., it depends on old variables γ α and new momenta Π α —and from which the rest of the variables are given by:
π α ( k ) = G ( γ β , Π β ) γ α ( k ) , Γ α ( k ) = G ( γ β , Π β ) Π α ( k ) .
where B α β and A α β are matrices whose components depend on background variables, but not on perturbations, and furthermore, A α β is symmetric. Equations (32) 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:
Γ 0 ( k ) = γ 0 + κ p ϕ 1 / 6 κ a π a + a 3 σ ( 2 ) 2 γ 1 γ 2 ,
Γ 1 ( k ) = γ 5 + a 2 σ ( 5 ) 1 / 6 κ π a + a 2 σ ( 2 ) 2 γ 1 γ 2 ,
Γ 2 ( k ) = γ 6 + a 2 σ ( 6 ) 1 / 6 κ π a + a 2 σ ( 2 ) 2 γ 1 γ 2 ,
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 [16]. 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 R ( k ) 1 4 κ a z Γ 0 , where z = 6 κ p ϕ π a = ϕ ˙ H a .
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 the 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 (32).

4.2. Implementation in Mathematica

We summarize here the main steps of the code [16]. We start defining the new variables Γ α ( k ) , Π α ( k ) :
In[62] := DefTensor[ Γ 0[],{M3,t}];
   DefTensor[ Π 0[],{M3,t}];
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 (35) above). These variables are named in the code as Γ 0new, Γ 1new, Γ 2new:
In[63] := Γ 0new= γ 0[]+(6 Sqrt[2 κ ] P ϕ γ 1[])/((Deth[] (1/6)) (Sqrt[6] κ π a[]
   +6 (Deth[] (1/6)) 2 σ 2[]))-(6 Sqrt[ κ ] P ϕ []  γ 2[])/((Deth[] (1/6)) (Sqrt[6]  κ   π a[]
   +6 (Deth[] (1/6)) 2 σ 2[]));
and similarly for Γ 1new, Γ 2new. With this choice, we can determine some of the coefficients B α β . The remaining coefficients in B α β , namely those with α , β = 3 , 4 , 5 , 6 , are written as unknowns CI[], EI[], FI[], and JI[], with I= { 0 , , 6 } . For example, we define:
In[64] := DefTensor[C0[],M3,t];
   DefTensor[E0[],M3,t];
   DefTensor[F0[],M3,t];
   DefTensor[J0[],M3,t];
and similar definitions for the remaining unknowns. We use these coefficients to define:
In[65] := Γ 3new= γ 0[]C0[]+…+ γ 6[]C6[];
    Γ 4new= γ 0[]E0[]+…+ γ 6[]E6[];
    Γ 5new= γ 0[]F0[]+…+ γ 6[]F6[];
    Γ 6new= γ 0[]J0[]+…+ γ 6[]J6[];
The coefficients A α β are denoted by AIJ[], with I,J= { 0 , , 6 } . For instance, we define:
In[66] := DefTensor[A00[],M3,t];
and similarly for the remaining ones.
Using these coefficients and the definitions of new variables in terms of them, we introduce the definition of the generating function as:
In[67] := G= Γ 0new Π 0[]+ Γ 1new Π 1[]++ Γ 5new Π 5[]+ Γ 6new Π 6[]+A00[]  γ 0[]  γ 0[]
   +2A01[] γ 0[] γ 1[]++2A56[] γ 5[] γ 6[]+A66[] γ 6[] γ 6[]//NoScalar;
We now follow the strategy described above. Namely, we first obtain an expression for the old momenta π α in terms of γ α and Π β by taking the 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 (30), 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 [16] for details.

5. Dynamics of Gauge Invariant Perturbations

5.1. 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 i ( 2 ) 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 i ( 2 ) , 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 N S ˜ ( 2 ) ( k ) , where S ˜ ( 2 ) ( k ) is quadratic in perturbations. In each term, one perturbation is evaluated at k and the other at k .
From the expression for S ˜ ( 2 ) ( k ) in terms of δ h ˜ i j and δ π ˜ i j , we obtain a Hamiltonian for the new variables Γ α ( k ) and Π α ( k ) , by first implementing the change from δ h ˜ i j ( k ) , δ π ˜ i j ( k ) to γ ( k ) , π ( k ) and then from the latter 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:
G γ ( k ) = δ π ˜ i j ( k ) n = 1 6 A i j ( n ) ( k ) γ n ( k ) ,
where we have chosen it to be of the third type (i.e., it depends on new configuration variables and old momenta). Recall that the matrices A i j ( n ) ( k ) depend on time. The second canonical transformation from γ α ( k ) , π α ( k ) to Γ α ( k ) , Π α ( k ) is defined by the generating function (31).
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):
H pert = N ( t ) V 0 2 a ( t ) k μ , μ = 0 2 4 κ a 2 ( t ) δ μ , μ | Π μ ( k ) | 2 + a 2 ( t ) 4 κ δ μ , μ k 2 ( t ) + U μ μ ( t , k ) Γ μ ( k ) Γ ¯ μ ( k ) ,
where δ μ , μ is the Kronecker delta and k 2 ( t ) a 2 ( t ) k i k j = a 2 ( t ) k 1 2 a 1 2 ( t ) + k 2 2 a 2 2 ( t ) + k 3 2 a 3 2 ( t ) . If we choose N = 1 , this Hamiltonian generates evolution in proper time t and in conformal time if N = a . The (time-dependent) effective potentials U μ μ are symmetric in μ and μ , and they become diagonal (i.e., proportional to δ μ μ ) in the isotropic limit. However, in the 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):
Γ ˙ μ ( k ) = { Γ μ ( k ) , H pert } = 4 κ a 3 Π μ ( k ) , Π ˙ μ ( k ) = { Π μ ( k ) , H pert } = a 4 κ μ = 0 2 ( δ μ μ k 2 + U μ μ ) Γ μ ( k ) .
Combining these equations into second order differential equations, we obtain:
Γ ¨ μ + 3 H Γ ˙ μ + k 2 a 2 Γ μ + 1 a 2 μ = 0 2 U μ μ Γ μ = 0 .
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

5.2. Implementation in Mathematica

We start with the expression k N S ˜ ( 2 ) ( k ) in terms of δ h ˜ i j ( k ) and δ π ˜ i j ( k ) . As mentioned above, each term S ˜ ( 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 for S ˜ ( 2 ) ( k ) in terms of δ h ˜ i j ( k ) and δ π ˜ i j ( k ) :
In[68] := S2a=SeriesCoefficient[S,2];
In[69] := S2a/.MakeRule[PD[-i]@PD[-j]@ δ h[LI[1],-k,-l],-kv[-i]kv[-j] δ h[LI[1],-k,-l]]
   /.MakeRule[PD[-i]@ δ ϕ [LI[1]]PD[-j]@ δ ϕ [LI[1]],kv[-i]kv[-j] δ ϕ [LI[1]] δ ϕ [LI[1]]]
   /.MakeRule[PD[-i]@ δ h[LI[1],-j,-k] PD[-d]@ δ h[LI[1],-l,-m],kv[-i]kv[-l]
    δ h[LI[1],-j,-k] δ h[LI[1],-l,-m]]/.MakeRule[h[i,j]kv[-i]kv[-j],k[] 2]
   /.MakeRule[PD[-j]@PD[-l]@ δ h[LI[2],-i,-k],0]
   /.MakeRule[PD[-k]@PD[-l]@ δ h[LI[2],-i,-j],0]/.bgmomrule//org;
where, as before, we have replaced spatial derivatives by i k . The second step is to move from δ ϕ ˜ ( k ) , δ p ˜ ϕ ( k ) , δ h ˜ i j ( k ) , δ π ˜ i j ( k ) to γ α ( k ) , π β ( k ) :
In[70] := %/.moderule1/.moderule2/. γ 0rule/. π 0rule//org;
and we also decompose the shear in its components σ ( n ) :
In[71] := %/.sheardecomposition//org;
We further simply our final result with:
In[72] := S2b=%/.A2rule/.A3rule/.A4rule/.A5rule/.A6rule/. σ brule//org;
Finally, we want to write the Hamiltonian in terms of the new momenta Π α and configuration variables Γ α . In order to do so, we introduce in the notebook the expressions for old configuration variables as functions of new ones as γ 0old1, , γ 6old1. Similarly, we introduce expressions for old momenta in terms of new momenta and old configuration variables, which are denoted by π 0old1, , π 6old1. One can easily combine these two sets of expressions to replace old configuration variables and momenta in any expressions by new ones. Thus, we obtain:
In[73] := S2b/.MakeRule[ π 0[], π 0old1]/.MakeRule[ π 1[], π 1old1]/.MakeRule[ π 2[], π 2old1]
   /.MakeRule[ π 3[], π 3old1]/.MakeRule[ π 4[], π 4old1]/.MakeRule[ π 5[], π 5old1]
   /.MakeRule[ π 6[], π 6old1];
   S2c=%/.MakeRule[ γ 0[], γ 0old1]/.MakeRule[ γ 1[], γ 1old1]/.MakeRule[ γ 2[], γ 2old1]
   /.MakeRule[ γ 3[], γ 3old1]/.MakeRule[ γ 4[], γ 4old1]/.MakeRule[ γ 5[], γ 5old1]
   /.MakeRule[ γ 6[], γ 6old1];
As explained above, the final Hamiltonian is obtained from this expression by adding the time derivative of the generating functions (36) and (31). We denote by dG γ dt1 and dGdt2 their time derivatives in the companion notebook, respectively. The addition of these terms requires some additional simplifications. Concretely, we apply the following rule in several intermediate steps:
In[74] := bgconstraintrule=MakeRule[V[ ϕ []],-1/(2 κ )( σ 2[] 2+ σ 3[] 2+ σ 4[] 2+ σ 5[] 2+ σ 6[] 2)
   -1/Deth[] (1/2)(P ϕ [] 2/(2Sqrt[Deth[]])-( κ π a[] 2)/(12Deth[] (1/6)))];
It uses the background constraint to replace the potential V ( ϕ ) of the scalar field by the other background variables. We get:
In[75] := S2d = S2c+dGdt1+dG γ dt2/.bgconstraintrule;
This expression can be further simplified. On the one hand, it contains terms proportional to the new momenta Π 3 , , Π 6 , which, by construction, are constrained to vanish. Therefore, we set them equal to zero with the rule:
In[76] := S2e=S2d/.MakeRule[ Π 3[],0]/.MakeRule[ Π 4[],0]/.MakeRule[ Π 5[],0]/.MakeRule[ Π 6[],0]
   //org;
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 ( Γ 0[] Π 0[]), ( Γ 1[] Π 0[]), ( Γ 2[] Π 0[]), ( Γ 1[] Π 1[]), ( Γ 1[] Π 2[]), and ( Γ 2[] Π 2[]), to vanish. These conditions fix the coefficients A01[], A05[], A06[], A56[], A12[], A26[] in (31). As a bonus, this automatically guarantees that the remaining cross terms ( Γ 0[] Π 1[]), ( Γ 1[] Π 0[]), ( Γ 0[] Π 2[]), ( Γ 2[] Π 0[]), ( Γ 1[] Π 2[]), and ( Γ 2[] Π 1[]) all vanish. As mentioned before, we set the remaining free parameters to zero, in order to have a Hamiltonian as simple as possible.
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:
In[77] := H2g=H2f/.bgconstraintrule//org;
This final Hamiltonian can be simplified and written in the form given in (37). Equations of motion for the gauge invariant variables can be derived straightforwardly from it.

6. 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 [16] 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 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 complemented this analysis with a computer code, based on the C programing language, available in [17], 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.   

Author Contributions

All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by the NSF CAREER grant PHY-1552603 and the Hearne Institute for Theoretical Physics. We acknowledge the use of high performance computing resources provided by Louisiana State University (http://www.hpc.lsu.edu), Baton Rouge, USA.

Acknowledgments

We have benefited from discussions with Abhay Ashtekar, Mar Bastero-Gil, Brajesh Gupt, Guillermo A. Mena Marugán, Jorge Pullin, Parampreet Singh, and Edward Wilson-Ewing.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Hamiltonian of Gauge Invariant Perturbations

The potentials U μ μ in the Hamiltonian (37) that generates the dynamics of gauge invariant perturbations are given by the following expressions:
U 00 = a 2 V ϕ ϕ 2 κ p ϕ 2 F 2 a 3 + 2 κ F 1 κ p ϕ 2 p a 3 a 5 + 2 V ϕ p ϕ , U 01 = U 10 = 2 κ a 2 a 2 p ϕ σ ( 5 ) F 2 + a 5 V ϕ σ ( 5 ) F 1 a 2 p ϕ G 5 F 1 + κ 6 p ϕ p a σ ( 5 ) F 1 , U 02 = U 20 = 2 κ a 2 a 2 p ϕ σ ( 6 ) F 2 + a 5 V ϕ σ ( 6 ) F 1 a 2 p ϕ G 6 F 1 + κ 6 p ϕ p a σ ( 6 ) F 1 , U 12 = U 21 = 2 σ ( 5 ) σ ( 6 ) a 2 a 3 F 2 + 2 3 κ a p a F 1 2 a 3 σ ( 6 ) G 5 + 2 a 3 σ ( 5 ) G 6 F 1 U 22 = 2 a 2 σ ( 5 ) 2 + κ p a σ ( 2 ) 6 a 2 2 3 G 2 + 4 3 κ a p a σ ( 6 ) 2 F 1 4 a 3 σ ( 6 ) F 1 G 6 2 a 3 σ ( 6 ) 2 F 2 ,
with V ϕ d V / d ϕ , V ϕ ϕ d 2 V / d ϕ 2 , and:
F 1 = κ p a 2 a 3 + 3 2 σ ( 2 ) a 2 κ ρ + σ ( 3 ) 2 + σ ( 4 ) 2 + σ ( 5 ) 2 + σ ( 6 ) 2 , F 2 = 3 κ V a κ 2 p a 2 3 a 5 + κ p a σ ( 2 ) 2 6 a 3 + 3 2 G 2 a F 1 κ 2 p ϕ 2 p a a 8 + 2 σ ( 3 ) G 3 + 2 σ ( 4 ) G 4 + 2 σ ( 5 ) G 5 + 2 σ ( 6 ) G 6 ) 2 κ ρ + σ ( 3 ) 2 + σ ( 4 ) 2 + σ ( 5 ) 2 + σ ( 6 ) 2 , G 2 = κ p a σ ( 2 ) 2 a 2 3 2 σ ( 3 ) 2 + σ ( 4 ) 2 , G 3 = κ p a σ ( 3 ) 2 a 2 + 1 2 3 σ ( 2 ) σ ( 3 ) σ ( 3 ) σ ( 5 ) σ ( 4 ) σ ( 6 ) , G 4 = κ p a σ ( 4 ) 2 a 2 + 1 2 3 σ ( 2 ) σ ( 4 ) + σ ( 4 ) σ ( 5 ) σ ( 3 ) σ ( 6 ) , G 5 = κ p a σ ( 5 ) 2 a 2 + 1 2 ( σ ( 3 ) 2 σ ( 4 ) 2 ) , G 6 = κ p a σ ( 6 ) 2 a 2 + 2 σ ( 3 ) σ ( 4 ) .
Note that the expressions above have an implicit dependence on k coming from σ ( n ) ( k ) .

References

  1. Bardeen, J.M. Gauge-invariant cosmological perturbations. Phys. Rev. D 1980, 22, 1882. [Google Scholar] [CrossRef]
  2. Arnowitt, R.L.; Deser, S.; Misner, C.W. The Dynamics of general relativity. arXiv 2004, arXiv:gr-qc/0405109. [Google Scholar]
  3. Langlois, D. Hamiltonian formalism and gauge invariance for linear perturbations in inflation. Class. Quantum Grav. 1994, 11, 389. [Google Scholar] [CrossRef]
  4. Agullo, I.; Ashtekar, A.; Nelson, W. Extension of the quantum theory of cosmological perturbations to the Planck era. Phys. Rev. D 2013, 87, 043507. [Google Scholar] [CrossRef] [Green Version]
  5. Agullo, I.; Ashtekar, A.; Nelson, W. The pre-inflationary dynamics of Loop Quantum Cosmology: Confronting quantum gravity with observations. Class. Quantum Grav. 2013, 30, 085014. [Google Scholar] [CrossRef] [Green Version]
  6. Castelló Gomar, L.; Fernández-Méndez, M.; Mena Marugán, G.A.; Olmedo, J. Cosmological perturbations in Hybrid Loop Quantum Cosmology: Mukhanov-Sasaki variables. Phys. Rev. D 2014, 90, 064015. [Google Scholar] [CrossRef] [Green Version]
  7. Castelló Gomar, L.; Martín Benito, M.; Mena Marugán, G.A. Gauge-Invariant Perturbations in Hybrid Quantum Cosmology. JCAP 2015, 6, 045. [Google Scholar] [CrossRef] [Green Version]
  8. Agulló, I.; Olmedo, J.; Sreenath, V. Hamiltonian theory of classical and quantum perturbations in Bianchi I spacetimes. 2020; in preparation. [Google Scholar]
  9. Pereira, T.S.; Pitrou, C.; Uzan, J.-P. Theory of cosmological perturbations in an anisotropic universe. JCAP 2007, 0709, 006. [Google Scholar] [CrossRef] [Green Version]
  10. Pitrou, C.; Pereira, T.S.; Uzan, J.-P. Predictions from an anisotropic inflationary era. JCAP 2008, 0804, 004. [Google Scholar] [CrossRef]
  11. Martín-García, J.M. xAct, Efficient Tensor Computer Algebra for Mathematica. 2004. Available online: http://www.xact.es (accessed on 19 February 2020).
  12. Brizuela, D.; Martín-García, J.M.; Marugan, G.A.M. Second-and higher-order perturbations of a spherical spacetime. Phys. Rev. D 2006, 74, 044039. [Google Scholar] [CrossRef] [Green Version]
  13. Brizuela, D.; Martín-García, J.M.; Marugan, G.A.M. xPert: Computer algebra for metric perturbation theory. Gen. Rel. Grav. 2009, 41, 2415. Available online: http://www.xact.es/xPert (accessed on 19 February 2020). [CrossRef] [Green Version]
  14. Brizuela, D.; Martín-García, J.M.; Tiglio, M. Complete gauge-invariant formalism for arbitrary second-order perturbations of a Schwarzschild black hole. Phys. Rev. D 2009, 80, 024021. [Google Scholar] [CrossRef] [Green Version]
  15. Brizuela, D.; Martín-García, J.M.; Sperhake, U.; Kokkotas, K.D. High-order perturbations of a spherical collapsing star. Phys. Rev. D 2010, 82, 104039. [Google Scholar] [CrossRef] [Green Version]
  16. Agullo, I.; Olmedo, J.; Sreenath, V. Available online: http://bitbucket.org/jolmedo/bianchii-perts/src/master/ (accessed on 19 February 2020).
  17. Olmedo, J.; Agullo, I.; Sreenath, V. Available online: http://bitbucket.org/jolmedo/cosmo-perts/src/master/ (accessed on 19 February 2020).
  18. Goldberg, J.; Newman, E.T.; Roveli, C. On Hamiltonian systems with first-class constraints. J. Math. Phys. 1991, 32, 2739. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Agullo, I.; Olmedo, J.; Sreenath, V. xAct Implementation of the Theory of Cosmological Perturbation in Bianchi I Spacetimes. Mathematics 2020, 8, 290. https://doi.org/10.3390/math8020290

AMA Style

Agullo I, Olmedo J, Sreenath V. xAct Implementation of the Theory of Cosmological Perturbation in Bianchi I Spacetimes. Mathematics. 2020; 8(2):290. https://doi.org/10.3390/math8020290

Chicago/Turabian Style

Agullo, Ivan, Javier Olmedo, and Vijayakumar Sreenath. 2020. "xAct Implementation of the Theory of Cosmological Perturbation in Bianchi I Spacetimes" Mathematics 8, no. 2: 290. https://doi.org/10.3390/math8020290

APA Style

Agullo, I., Olmedo, J., & Sreenath, V. (2020). xAct Implementation of the Theory of Cosmological Perturbation in Bianchi I Spacetimes. Mathematics, 8(2), 290. https://doi.org/10.3390/math8020290

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop