xPand: An algorithm for perturbing homogeneous cosmologies

xPand: An algorithm for perturbing homogeneous cosmologies

Cyril Pitrou, Xavier Roy and Obinna Umeh Institut d’Astrophysique de Paris, Université Pierre & Marie Curie—Paris VI, CNRS-UMR 7095, 98 bis, Bd Arago, 75014 Paris, France
Institut Lagrange de Paris, Sorbonne Universités, 98 bis Bd Arago, F-75014 Paris, FranceAstrophysics, Cosmology and Gravity Centre (ACGC), and Department of Mathematics and Applied Mathematics,Cape Town University, Rondebosch 7701, South Africa
September 12, 2019

In this paper, we develop in detail a fully geometrical method for deriving perturbation equations about a spatially homogeneous background. This method relies on the splitting of the background space–time and does not use any particular set of coordinates: it is implemented in terms of geometrical quantities only, using the tensor algebra package xTensor in the xAct distribution along with the extension for perturbations xPert. Our algorithm allows one to obtain the perturbation equations for all types of homogeneous cosmologies, up to any order and in all possible gauges. As applications, we recover the well-known perturbed Einstein equations for Friedmann–Lemaître–Robertson–Walker cosmologies up to second order and for Bianchi I cosmologies at first order. This work paves the way to the study of these models at higher order and to that of any other perturbed Bianchi cosmologies, by circumventing the usually too cumbersome derivation of the perturbed equations.

02.70.Wz, 98.80.Jk, 98.80.-k


Cosmological perturbation theory constitutes the cornerstone of our current understanding of the origin, evolution and formation of large-scale structures. The evolution history of perturbations is written about a fixed homogeneous and isotropic Friedmann–Lemaître–Robertson–Walker (FLRW) background space–time, and the interpretation of cosmological observations (such as WMAP, QUIET and ACT Hinshaw:2012fq ; Sievers:2013wk ; Araujo:2012yh ), within this model, converge toward a coherent and unified picture of the Universe. This picture is likely to get even clearer when the next generation of large-array cosmological observations (such as EUCLID and SKA) becomes operational. It is widely expected that these observations will generate large amount of data that will provide a percent level accuracy for the cosmological parameters.

Cosmological perturbation equations are simple and straightforward to derive at linear order, but they are inadequate for understanding the late-time evolution of the Universe, precisely when the nonlinear gravitational effects, carrying the information of the physics of current interest, come into play. Going beyond first order is a difficult task, and in some cases it becomes extremely arduous to even perform a coordinate or gauge transformation at nonlinear order by hand. To the best of our knowledge, there is no available easy-to-use software designed for cosmology and capable of deriving all equations of motion for perturbed variables. The closely related available option is the GRTensor package GRTensor , which runs on Maple or Mathematica. However, the outputs generated at linear order are already a bit complicated to understand, let alone its outputs at nonlinear orders, as it relies exclusively on a properly defined set of background coordinates each time it acts on a perturbed variable.

To fill up this gap, we have developed an algebra package for cosmological perturbation theory, called xPand xPand , which uses the tools of the tensor algebra package xTensor and an extension for perturbation, xPert xAct ; Brizuela:2008ra . The xTensor and xPert packages are part of the xAct distribution xAct that run on Mathematica and they are available under the General Public License. xPert is specifically designed to perform perturbations on arbitrary background space–times Brizuela:2006ne ; Brizuela:2008ra , but it lacks the features for specializing to a specific background space–time and a specific form for the metric perturbations, as is needed in the case of cosmology. In Brizuela:2009qd ; Brizuela:2010qu these packages were used to study perturbations about a spherically symmetric space–time, more precisely around a Schwarzschild solution of the Einstein field equations. The xPand package now allows one to derive, in a simple and user-friendly manner, all the necessary equations for cosmological perturbation variables, around any homogeneous background space–times, at any order and in any gauge. Specifically, the available type of cosmologies cover the Minkowski, FLRW (flat and curved), and Bianchi space–times, and the available predefined gauge are: general gauge (no gauge choice), comoving gauge, flat gauge, isodensity gauge, Newtonian gauge and synchronous gauge.

The paper is organized as follows. In section I, we provide a general overview of the mathematical framework on which xPert is built. In section II, we detail the splitting of the background manifold into a family of homogeneous hypersurfaces orthogonal to a fundamental observer’s velocity. In section III, we decompose the perturbed metric with respect to this foliation and define the scalar, vector and tensor perturbations. Each of these sections is supplemented by a presentation of the associated implementation in xPand, by means of several detailed examples. Finally, we summarize and discuss in section IV the features and performances of the package.

I Perturbations around a general space–time

In this section we briefly review the algorithm of xPert Brizuela:2008ra , which constitutes the basis of our method. For more details about perturbation theory in the context of cosmology, we refer the reader to, e.g., Bruni1997 ; Nakamura2007 .

i.1 General framework

Let us consider in what follows a background manifold along with its perturbed manifold . Both are related by means of a diffeomorphism . Tensorial quantities are thus transported from one manifold to the other with the help of the pull-back , the push-forward , and their respective inverses. The metric of the perturbed manifold relates to that of the background as


Here and in the sequel, we use boldface symbols for tensorial quantities, an over-bar for background quantities111 This convention differs from the one adopted in Brizuela:2008ra . We however opt for this choice as it reflects the standard usage in cosmological perturbation theory. , and the notation (resp. ) for the total (resp.  order) perturbation of a tensor . One may prefer to write instead of , as the definition of the perturbations depends on the diffeomorphism , that is on the choice of the gauge. We however choose to omit this reference for the sake of clarity, and in order not to burden the notation unnecessarily, we shall moreover use the short-hand: , for any perturbed quantity.

Unless otherwise specified, when we write down the components of a tensor, these should be understood as expressed in a general basis222 We use Greek letters (, , , , , …) for space–time indices. (this holds equally for the background and perturbed tensors, and for the perturbations). Since all perturbation orders live on the background manifold, as they are the result of the pull-back of a tensorial quantity living on the perturbed manifold, we shall raise and lower indices using the background metric. We have for instance


for the order metric perturbations .

i.2 Expansion of the curvature tensors

i.2.1 Mathematical framework

The inverse of the metric tensor is given by the relation


Expanding it into


and making use of the definition , we obtain the order perturbation of :


where the sum runs over the sorted partitions of for positive integers, such that . Note that (e.g., we have at first order: ).

By means of relation (5), we can express the perturbation of the connection components as Brizuela:2008ra


where the last term of the right-hand side is defined by


The perturbation of the Riemann tensor is given in all generality by


and for a connection compatible with the metric, we have, from equation (6),


The symbol denotes the repetition of the preceding expression with indices and exchanged. The perturbation of the Ricci tensor is simply obtained by contracting the second and fourth indices of in the previous expression, and the perturbation of the Ricci scalar, , is written:


At last, the perturbation of the Einstein tensor is expressed according to


i.2.2 Implementation in xPert

All the perturbative expansions expounded above are implemented in the package xPert Brizuela:2008ra . For completeness, we here briefly review its main commands. The package can be loaded by evaluating

In[1] :=  <<xAct‘xPert‘

(Version and copyright messages)

We first define the four-dimensional manifold M with abstract indices :

In[2] :=  DefManifold[ M, 4, ];

and then we define the ambient metric g of negative signature, along with its associated covariant derivative CD:

In[3] :=  DefMetric[ -1, g[-,-], CD, {";",""}, PrintAs->"" ];

where and g respectively correspond to and . Several tensors related to this metric are automatically defined at the same time (e.g., all the curvature tensors). Note that in xTensor, the covariant indices of a tensor are represented with a minus sign (g[-, -] means ), while the latter is omitted for contravariant indices (g[, ] means ).

Upon defining the perturbations dg of the metric g with the command:

In[4] :=  DefMetricPerturbation[ g, dg, ];

where is the perturbative parameter to be used in the expansions, it becomes feasible to evaluate the perturbation of any tensor associated with the metric. For instance, the perturbation at first order of the Ricci scalar is simply obtained by evaluating

In[5] :=  ExpandPerturbation@Perturbed[ RicciScalarCD[], 1 ] // ContractMetric // ToCanonical

Out[5] := 

The functions ExpandPerturbation and Perturbed are used to evaluate the perturbation of any expression up to any order (here at first order), the function ContractMetric removes the background metric tensor through contraction on dummy indices (i.e. repeated indices), and the function ToCanonical simplifies the result, gathering together the terms which are equal up to symmetries. Further details can be found in Brizuela:2008ra .

i.3 Conformal transformation

In xTensor, the first metric defined on the manifold is the one that is used to raise and lower any tensor indices. For our purpose, we have chosen it to be the conformal metric (cf input In[3]), which is different from, but conformally related to, the background metric of the physical space–time. This choice ensures that the conventional way of moving the indices of perturbed fields (equation (2)) is well satisfied within our algorithm. We however need, now, to relate the tensorial quantities one considers in perturbation theory, namely those living on the background manifold of the physical space–time, to those we have defined or shall define on the conformal background manifold M.

We detail this important point in the rest of this section.

i.3.1 Mathematical framework

Let us denote by the metric of the physical space–time, and by its background value. The metrics and are related by the conformal transformation


with being the scale factor of the background space–time333We stress again that the first metric defined on M is the one that is used to raise and lower any tensor indices. Defining another metric, say , on the same manifold from the function DefMetric actually creates two objects: (i) the tensor , with internal notation f, and (ii) the tensor , with internal notation Inv[f]. These have the following properties:

which explains the notation we use in equation (12). . Substituting into the first expression of (12) the perturbative expansion (1) for and its counterpart for , we extend the conformal transformation to the background:


and to the perturbed level:


The associated Levi-Civita connections and , on the one hand, and and , on the other, are related by444By definition the scale factor used in the transformation (12) is not to be perturbed. Hence we have: , and thus: .


for any 1-form . Using equations (12) and (13), we can write the quantities and as


with and , and where the parentheses indicate symmetrization over the indices enclosed. We can then formulate the order perturbation of as555The background and perturbed connections, on the one hand, and the background and conformal connections, on the other hand, are respectively related by

for any 1-form . Performing a conformal transformation on the former expression and perturbing the latter respectively yields
These relations can only be compatible if
It can be checked directly from equations (6), (14) and (17) that this is indeed the case. This shows the equivalence between the transformations and . The latter approach actually proves itself to be faster within xPand. It is therefore the one that we have coded in the function ToxPand (see section IV.2).


With the help of the two previous relations, we are now able to provide the correspondence we seek.

For instance, the Riemann tensor associated with the metric is given by (see appendix D of 1984ucp..book…..W )


where the brackets indicate anti-symmetrization over the indices enclosed. Perturbing this expression and using equation (17), we can finally relate to and recover the usual quantities studied in perturbation theory.

i.3.2 Implementation in xPand

The xTensor package provides the tools to define a metric conformally related to another, thanks to the option ConformalTo of the function DefMetric. We have encapsulated this in xPand within the function DefConformalMetric, which furthermore ensures the transitivity of several conformal transformations.

Let us load the package xPand:

In[6] :=  <<xAct‘xPand‘

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Package xAct‘xPand‘ version 0.4.0, {2013,02,08}

CopyRight (C) 2012-2013, Cyril Pitrou, Xavier Roy and Obinna Umeh under the GPL.

By evaluating the command

In[7] :=  DefConformalMetric[ g, a ];

we define the scalar factor a[] and the metric ga2 conformally related to g through a (ga2 thus corresponds to the background metric of the physical space–time). To obtain the expression of any tensorial quantities living on the manifold described by ga2 in terms of those defined on M, one then simply has to use the xPand function Conformal.

For instance, we obtain for the Riemann tensor associated with ga2 the expression

In[8] :=  Conformal[ g, ga2 ][ RiemannCD[-,-,-,] ]

Out[8] := 
Out[8] :=    

which coincides with equation (18).

The conformal transformation may as well be performed on quantities that are not related to a metric. For a general tensor, it is defined as


where is the conformal weight of the tensor . The default value of is chosen to be zero in order to leave the norm invariant under a conformal transformation:

In[9] :=  DefTensor[ W[-], M ];

In[10] :=  Conformal[ g, ga2 ][ W[-] ]

Out[10] := 

In[11] :=  Conformal[ g, ga2 ][ W[] ]

Out[11] := 

The conformal weight can however be modified for each tensor with the xPand function ConformalWeight. This can be of use for instance to preserve the geodesic character of light-like vectors666We will make use of this prescription in a future version of xPand to implement the derivation of the (perturbed) null geodesic equation. :

In[12] :=  DefTensor[ k[-], M ];

In[13] :=  ConformalWeight[ k ] ^ = -1;

In[14] :=  ConformalWeight[ k[-] ]

Out[14] := 

In[15] :=  ConformalWeight[ k[] ]

Out[15] := 

In[16] :=  Conformal[ g, ga2 ][ k[] CD[-]@k[-] ]

Out[16] := 

So far, by applying Conformal then ExpandPerturbation@Perturbed on a given expression defined on M, one obtains the perturbation of its conformal transformation in terms of the tensors defined on M, the metric g, its perturbations dg, the connection , and the scale factor a. To end up with the usual expressions of perturbation theory, one needs to perform a splitting of the background manifold (section II), then decompose each perturbed fields into its spatial and temporal parts and finally parameterize the perturbations of the metric (section III).

Ii splitting of the background manifold

ii.1 Induced metric

The assumption that the background space–time possesses a set of (three-dimensional) homogeneous surfaces provides a natural choice for the slicing. We foliate the background manifold by means of this family, and we denote by the unit time-like vector () normal to it. The metric of is decomposed as


where represents the induced metric of the spatial hypersufaces777For the sake of clarity, let us note that are not the perturbations of . From the definition of together with relation (20), we can actually relate them as . .

The acceleration of the so-called Eulerian observers satisfies in all generality Gourgoulhon:2012ue


with being the lapse function. stands for the connection of the three-surfaces associated with (), and it is related to the four-covariant derivative as


for any spatial tensor field888We recall that the operator loses its character of derivative when it is applied to non-spatial tensors. More precisely, we are not allowed to use the Leibniz rule anymore, as one can realize upon writing for instance

for any scalar field . . Since the lapse is homogeneous in the configuration at stake, the acceleration vanishes and the observers are in geodesic motion. We can therefore label each hypersurface by their proper time and write: . In addition, being hypersurface-forming by construction, its vorticity vanishes; this property yields


where the equivalence stems from the null acceleration. For comprehensive reviews on the formalism, we refer the reader to, e.g., Smarr ; Gourgoulhon:2012ue .

ii.2 Extrinsic curvature

Another tensor we shall make use of is the symmetric extrinsic curvature tensor, which characterizes the way the three-surfaces are embedded into the background manifold. It satisfies the relation


where we have chosen a positive sign for the right-hand side999This convention does not affect the Einstein equations as written in terms of the kinematical quantities of the Eulerian observers. . From the decomposition (20) along with the vanishing of the acceleration and the unitary of , we can reformulate expression (24) as


Since the volume expansion of the background space–time is entirely contained in the scale factor , and owing to the conformal transformation (13), the trace of the extrinsic curvature vanishes: . As a result, we have for general Bianchi cosmologies: , with being the shear of the Eulerian observers; and for FLRW cosmologies: .

ii.3 Curvature tensors

The splitting of the four-Riemann tensor can be constructed from its different projections onto the spatial slices and the congruence of the observers. It is written as


where stands for the three-Riemann curvature of the hypersurfaces. The over-dot indicates the covariant derivative along the world-lines of the observers (for any tensor field , we have: ). The purely spatial projection of (26) only calls upon the first two terms, and it drives the Gauss relation


The three-space and one-time projection gives, from the next two terms, the Codazzi relation


and the last non-null projection (two-space and two-time) provides, from the last two terms, an evolution equation for the extrinsic curvature. These decompositions can be performed in xAct with the function GaussCodazzi.

For FLRW space–times, the curvature tensors of the hypersurfaces cast the form


with being the curvature parameter (equal to zero for flat FLRW cosmologies). The corresponding expressions for general Bianchi space–times are more involved, as they require the introduction of the constants of structures. We detail their derivation in appendix B.

ii.4 Derivatives

In order to achieve the splitting, we are left with the decomposition of the covariant derivative in terms of the induced derivative . For general spatial tensors (namely, for spatial tensors defined within or defined within then mapped onto ), the relation between the two derivatives reads


Even though our formalism is purely geometrical, we aim at eventually providing, for the perturbations, partial differential equations with respect to the proper time of the Eulerian observers. When considering the four-dimensional basis built to address the Bianchi classification (refer to appendix B), the Lie derivative along the direction of precisely comes down to . It is accordingly more appropriate for our purpose to use the Lie derivative rather than the dot derivative.

The relation between and is written as101010Note that for a spatial tensor , the quantity is also spatial.


and provides us with the following reformulation of (30):


Finally, the last expression we shall need is the commutation rule between the derivatives and . For general spatial tensors, it is given by


where we have made use of relations (22) and (28) for its derivation.

ii.5 Implementation in xPand

The splitting of the background manifold is performed by the xPand function SetSlicing. It can be applied to the following spatially homogeneous cosmologies: "Minkowski", "FLFlat", "FLCurved", "BianchiI", "BianchiA", "BianchiB" and "Anisotropic".

ii.5.1 Construction of the spatial hypersurfaces

From the ambient metric g, SetSlicing first defines the unit normal vector n, the induced metric h of the hypersurfaces, the associated covariant derivative cd, the associated scale factor a[h] and the conformal metric gah2. It then specifies the expressions of the intrinsic and extrinsic curvature tensors according to the type of cosmologies chosen by the user (hence, the geometry of the model is fully described).

Let us illustrate with the cosmologies "FLCurved" and "BianchiA" (for the latter most sophisticated case, refer to appendix B). With the evaluation

In[17] :=  SetSlicing[ g, n, h, cd, {"|",""}, "FLCurved" ];

the extrinsic curvature K[h][-,-] is set to zero, and the expressions of the three-curvature tensors are implemented following relations (29). For instance, we have for the three-Ricci tensor

In[18] :=  Riccicd[-,-]

Out[18] := 

For "BianchiA" cosmologies,

In[19] :=  SetSlicing[ g, nA, hA, cdA, {"|",""}, "BianchiA" ];

SetSlicing defines the spatial constants of structure (see appendix B) and then constructs a function allowing one to express the intrinsic curvature tensors in terms of them. For the three-Ricci tensor, we then have

In[20] :=  RiccicdA[-b,-c] // ToConstantsOfStructure[hA]

Out[20] := 

which is equivalent to equation (60) 111111Equation (60) is recovered by making use of the Jacobi identity (44) on Out[20]. This latter relation is not implemented, as it implies a symmetry among several terms that xTensor does not yet handle. However, once the constants of structure are expanded according to the parameterization (47), the Jacobi identity reduces to the constraint (48), which is automatically applied in xPand. . The constants of structure can be further expanded, following the usual Schücking, Kundt and Behr (SKB) decomposition 1969CMaPh..12..108E , by means of the xPand function ToBianchiType[hA].

We summarize in table 1 the different evaluations performed by SetSlicing with respect to the type of cosmologies.

Space–time Extrinsic curvature Three-curvature tensors Constants of structure
Minkowski Null Null Null
FLFlat Null Null Null
FLCurved Null Eqs. (29)
BianchiI Null Null
BianchiA   Eqs. (59), (60), (61)
BianchiB Eqs. (59), (60), (61)
Anisotropic Eqs. (59), (60), (61)
Table 1: Evaluations performed by SetSlicing for the extrinsic curvature, the three-curvature tensors and the constants of structure, according to the type of homogeneous cosmologies. In terms of g (or, equivalently, ), the dynamics of "Minkowski" and "FLFlat" space–times are identical. The difference lies in the conformal transformation (13): for the former models, the scale factor is set to 1. For "BianchiA", "BianchiB" and "Anisotropic" space–times, the three-curvature tensors can be formulated in terms of the constants of structure by means of the function ToConstantsOfStructure[] and, for the two first models, can be further expanded with the help of ToBianchiType[]. The "Anisotropic" space–time is used for hypersurfaces whose dimension differs from three (for these models the SKB decomposition does not apply).

ii.5.2 Derivatives

For anisotropic cosmologies, SetSlicing creates a set of automatic rules to handle the action of the derivatives on the three-curvature tensors. The aim is to obtain a formulation in terms of the extrinsic curvature tensor and/or the constants of structure, which fully identify the geometry of the model. Note that such rules are not necessary for isotropic models: the expressions of the intrinsic curvature tensors are simple enough to let xTensor perform the evaluation.

The action of the derivative is implemented following equation (32), that of following equation (53) together with (55), and the application of on a three-curvature tensor automatically calls the function ToConstantsOfStructure[] so as to use the property (58). Regarding the three-covariant derivative, we have for instance

In[21] :=  cdA[-]@RiccicdA[-,-]

Out[21] := 
Out[21] :=    

Lastly, SetSlicing constructs an automatic rule to perform the commutation (33) for any expression. This serves to make sure that the Lie derivative will first act on a given tensor, so as to recover the usual formulation of the perturbed equations121212In line of this comment, let us note that SetSlicing also creates internal rules for the commutation of several . These will enforce the appearance of Laplacians, and move a covariant derivative closer to a tensor when a divergence is present. . For "FLCurved" cosmologies, we have

In[22] :=  DefTensor[ V[-], M, OrthogonalTo->{n[]}, ProjectedWith->{h[,-]} ];

In[23] :=  LieD[n[]]@cd[-]@V[-]

Out[23] :=  </