GR perturbations with Mathematica

Tensor calculus

Let us start with the definition of the metric and the coordinate set. From now on, unless otherwise specified, 
all indexes are lower indexes. The indexes 1,2,3 are spatial, 4 is the time index. We use  a conformal time metric.

Clear[metric];
metric:=DiagonalMatrix[
{a[t]^2,a[t]^2,a[t]^2,-a[t]^2}]

coord={x,y,z,t};

metric//TableForm

Next we define the functions that raise and lower the indexes of a tensor s using the metric. 
The Dot operation means tensor multiplication of the last index of the first tensor with the
 first index of the second tensor. So we are doing
(g^(μν))(Subscript[s, ν...])
and
(Subscript[g, μν])(s^(ν...))

raise[metric_,s_]:=Inverse[metric].s;
lower[metric_,s_]:=metric.s;

In Mathematica the operations with tensor can be defined in terms of Inner and Outer products. 
The product of two tensor A and B with index summation is Inner[Times,A,B,Plus] or
simply A.B (Dot product). The outer or tensor  product, without index summation, is Outer[Times,A,B] 
(this increases the tensor rank). This applies also to derivatives

?Inner

div[s_,var_]:=Inner[D,s,var,Plus];

vector={x,y,z,t};div[vector,coord]

gives the divergence (∂_μ)(s^μ) of the tensor s with respect to the variables var,
while

Clear[grad];grad[s_List,var_]:=Outer[D,s,var];

gives the gradient (∂_μ)(s^ν).

grad[vector,coord]//TableForm

But unfortunately this does not work for scalar functions, because Outer 
requires an indexed quantity. Then  for scalars we use the alternative definition

Clear[gradScal];gradScal[s_,var_]:=Map[D[s,#]&, var]

gradScal[phi[x,y,z,t],coord]

Here, Map applies the function D, derivative with respect to an arbitrary variable #, to each variable var. 
We can then define the Christoffel symbols (first index up) (Subscript[Γ,    μν]^α):

Clear[christoffel2];christoffel2[metric_,var_]:=Table[Sum[Inverse[metric][[i,j]]*(D[metric[[m,j]],var[[n]]]+
			D[metric[[n,j]],var[[m]]]-D[metric[[m,n]],var[[j]]]),{j,Length[var]}],{i,Length[var]},{m,Length[var]},{n,Length[var]}]/2

A faster evaluation is as follows

Clear[christoffel];christoffel[metric_,var_]:=Inverse[metric].(grad[metric,var]+
			Transpose[grad[metric,var],{1,3,2}]-Transpose[grad[metric,var],{2,3,1}])/2

Notice that with this definition the first index is the upper one. For instance, we have

christoffel[metric,coord][[4,4,4]]

It is useful also to derive the contracted Christoffel symbol

chrContracted[metric_,var_]:=gradScal[Log[(-Det[metric])^(1/2)],var]

The D'Alambertian (Subscript[∇, μ])((∇^μ)) s  of a scalar s can be evaluated as

dalambertian[metric_,s_,var_]:=(-Det[metric])^(-1/2)*
div[(-Det[metric])^(1/2)*raise[metric,gradScal[s,var]],var]

For instance,

Simplify[dalambertian[metric,phi[t],coord]]

this is the conformal H:  confH=a*H

Simplify[dalambertian[metric,phi[t],coord]/.a'[t]->h[t]*a[t]]

We can now evaluate the Ricci tensor, mixed components

Clear[ricci2];ricci2[metric_,var_]:=
Table[Sum[D[christoffel[metric,var][[i,j,k]],var[[i]]],{i,Length[var]}],{j,Length[var]},
{k,Length[var]}]-grad[ChrContracted[metric,var],var]+chrContracted[metric,var].christoffel[metric,var]-
		Table[Sum[(christoffel[metric,var].christoffel[metric,var])[[i,j,kk,i]],
				{i,Length[var]}],{j,Length[var]},{kk,Length[var]}]

A faster expression is as follows

Clear[ricci];ricci[metric_,var_]:=
div[Transpose[christoffel[metric,var],{3,2,1}],var]-grad[chrContracted[metric,var],var]+
		chrContracted[metric,var].christoffel[metric,var]-
		Table[Sum[(christoffel[metric,var].christoffel[metric,var])[[i,j,kk,i]],
				{i,Length[var]}],{j,Length[var]},{kk,Length[var]}]

As an example, the FRW Ricci tensor is

MatrixForm[ricci[metric,coord]]

hubrep={a'[t]->h[t]*a[t],a''[t]->h'[t]*a[t]+h[t]*a'[t]};
Simplify[MatrixForm[ricci[metric,coord]]/.hubrep]
Simplify[MatrixForm[ricci[metric,coord]]//.hubrep]

and the Ricci scalar is

Simplify[Tr[raise[metric,ricci[metric,coord]]]//.hubrep]

Simplify[Sum[(raise[metric,ricci[metric,coord]])[[i,i]],{i,Length[coord]}]//.hubrep]

Let us start now with the pertubation calculation. We need first a k-space unit vector and the  argument of the plane wave modes. 
Since the problem is isotropic, we can use  for k any mode we wish, provided its modulus is k

krad=k*(x+y+z)/Sqrt[3]; (* this is (OverVector[k]).(OverVector[r]) *)

fourier=Exp[I*krad];

Then we define the perturbed metric in NEWTONIAN or LONGITUDINAL gauge in k-space as

phim[x,y,z,t]->(E^((I k (x+y+z))/(Sqrt[3]))) phim[t]

Clear[pmetric];
pmetric:=DiagonalMatrix[
{a[t]^2*(1-2*e*phim[t]*fourier),a[t]^2*(1-2*e*phim[t]*fourier),a[t]^2*(1-2*e*phim[t]*fourier),a[t]^2*(1+2*e*psim[t]*fourier)}]

coord={x,y,z,t};

pmetric//TableForm

Notice the order parameter e

The perturbed 4-velocity vector, raised index (u^μ), is

uvup={e*u[t]*fourier/a[t],e*u[t]*fourier/a[t],e*u[t]*fourier/a[t],(1-e*psim[t]*fourier)/a[t]};
uvup

or, in terms of the velocity field divergence, θ = i (Subscript[k, μ])(u^μ)

Clear[uvup];uvup={e*theta[t]*fourier/(k*I*a[t]*Sqrt[3]),e*theta[t]*fourier/(k*I*a[t]*Sqrt[3]),e*theta[t]*fourier/(k*I*a[t]*Sqrt[3]),(1-e*psim[t]*fourier)/a[t]};
uvup/.e->0
uvup

and we also need  (Subscript[u, μ]). and  (u^ν)(Subscript[u, μ])

Clear[uv];uv=lower[pmetric,uvup];

Clear[prod];prod=Normal[Series[Outer[Times,uvup,uv],
	{e,0,1}]]; (* here we start by expanding in e and taking the first order *)

Assuming p = w * rho, the perturbed EM tensor is then, 
(mixed components)
(Subscript[T, μ]^ν)=(ρ+p) (u^ν)(Subscript[u, μ]) + p (Subscript[I, μ]^ν)=(1+w)ρ (u^ν)(Subscript[u, μ])+wρ (Subscript[I, μ]^ν)
with ρ = (Subscript[ρ, 0])(1+e δ exp(i(OverVector[k]).(OverVector[r])))

rho[x,y,z,t]->rho[t](1+e*dt[t]*fourier)

tem=Normal[Series[(rho[t]*(1+w[t]+e*(1+cs2[t])*d[t]*Exp[I*krad])*prod-IdentityMatrix[4]*rho[t]*(w[t]+e*cs2[t]*d[t]*Exp[I*krad])),{e,0,1}]];

Simplify[tem/.e->0]//TableForm

and lower components

temdo=lower[pmetric,tem];

The perturbed Christoffel symbols are

chr=Simplify[Normal[Series[christoffel[pmetric,coord],{e,0,1}]]];

For instance,

chr[[1,4,1]]

Now, the background covariant conservation equation is (Subscript[T,       ;ν]^(μν)) = 0

Clear[emcb]; emcb=
Table[SeriesCoefficient[Series[
Simplify[Sum[D[tem[[j,i]],coord[[j]]],{j,4}]-
Sum[chr[[j,kk,i]]*tem[[kk,j]],{j,4},{kk,4}]+
Sum[chr[[j,j,kk]]*tem[[kk,i]],{j,4},{kk,4}]],{e,0,1}],0],{i,4}];

emcb

psol=Simplify[Solve[emcb[[4]]==0,rho'[t]]][[1,1]]

psolh=Simplify[Solve[emcb[[4]]==0,rho'[t]][[1,1]]/.hubrep]

Now, the perturbed covariant conservation equation (making use also of the background equation) is

Clear[emc]; emc=Table[
SeriesCoefficient[Series[
Simplify[Sum[D[tem[[j,i]],coord[[j]]],{j,4}]-
Sum[chr[[j,kk,i]]*tem[[kk,j]],{j,4},{kk,4}]+
Sum[chr[[j,j,kk]]*tem[[kk,i]],{j,4},{kk,4}]],{e,0,1}],1]/.psol,{i,4}];

Conservation Equations

Matter conservation equations

dust={w[t]->0,w'[t]->0,cs2[t]->0,a'[t]->h[t]*a[t],a''[t]->h'[t]*a[t]+h[t]*a'[t]};

Eq. (0)

dpsolw=Simplify[Solve[emc[[4]]==0,d'[t]]][[1,1]]

Simplify[dpsolw/.dust]

thetasolw=Simplify[Solve[emc[[4]]==0,theta[t]]][[1,1]]

Simplify[thetasolw/.dust]

Eq. (1)

thetapsolw=Simplify[Solve[emc[[1]]==0,theta'[t]]][[1,1]]

Simplify[thetapsolw/.dust]

Einstein equations

We finally perturb the Einstein equations

ei=Normal[Series[ricci[pmetric,coord],{e,0,1}]];

eis=Tr[raise[pmetric,ei]];

The total lower-index EM tensor is then (in case you have several components)

Clear[temtot];temtot=temdo;

from which we get the various components

Simplify[SeriesCoefficient[Series[a[t]^2*tem[[4,4]],{e,0,1}],1]]/fourier

Simplify[SeriesCoefficient[Series[a[t]^2*tem[[1,1]],{e,0,1}],1]]/fourier

Einstein  equations

eqmatrix=ei-pmetric*eis/2-8*Pi*temtot;

unperturbed Friedmann equation:

Simplify[Solve[SeriesCoefficient[Series[eqmatrix[[4,4]],{e,0,1}]/.dust,0]==0,h[t]]]

perturbed equations:

Eq. (0,0)

grav00=Expand[Simplify[Solve[SeriesCoefficient[Series[eqmatrix[[4,4]],{e,0,1}],1]==0,phim[t]]][[1,1]]/.dust]

Eq. (0,i)

grav0i=Simplify[Solve[SeriesCoefficient[Series[Sum[eqmatrix[[4,i]],{i,3}],{e,0,1}],1]==0,phim'[t]]/.dust][[1,1]]

Eq. (i,j)

{i,j}={1,2};gravij=Simplify[Solve[SeriesCoefficient[Series[eqmatrix[[i,j]],{e,0,1}],1]==0,phim[t]][[1,1]]]

Scalar field

Let us now evaluate the KG equation:

Unperturbed KG equation:

phisecsol=Simplify[Solve[Simplify[(dalambertian[metric,phi[t],coord]-U'[phi[t]])]==0,phi''[t]][[1,1]]/.dust]

Perturbed KG equation :

kgequation[metric_,field_,coord_]:=dalambertian[pmetric,field,coord]-U'[field]

dpsecsol=Simplify[Solve[Simplify[(SeriesCoefficient[Series[kgequation[pmetric,phi[t]+e*dp[t]*fourier,
	coord],{e,0,1}],1])]==0,dp''[t]][[1,1]]/.dust]

Before completing the calculation, we need also the phi EM tensor (lower component)

Clear[temphi];temphi[metric_,coord_,field_]:=Outer[Times,gradScal[field,coord],
		gradScal[field,coord]]-(gradScal[field,coord].raise[metric,gradScal[field,coord]])*metric/2+U[field]*metric

Clear[temphi];temphi[metric_,coord_,field_]:=Outer[Times,gradScal[phi[t]+e*dp[t]*fourier,coord],
		gradScal[phi[t]+e*dp[t]*fourier,coord]]-
(gradScal[phi[t]+e*dp[t]*fourier,coord].raise[metric,gradScal[phi[t]+e*dp[t]*fourier,coord]])*metric/2-U[phi[t]+e*dp[t]*fourier]*metric

The unperturbed tensor is

Table[Simplify[SeriesCoefficient[Series[temphi[pmetric,coord,phi[t]+e*dp[t]*fourier][[i,j]],{e,0,1}],0]],{i,4},{j,4}]//TableForm

and the perturbed term is

Simplify[SeriesCoefficient[Series[temphi[pmetric,coord,phi[t]+e*dp[t]*fourier][[4,4]],{e,0,1}],1]]/fourier

The total energy-momentum tensor is now

Clear[temtot];temtot=temdo+temphi[pmetric,coord,phi[t]+e*dp[t]*fourier];

so we have now the following equations

eqmatrix=ei-pmetric*eis/2-8*Pi*temtot;

unperturbed

Simplify[Solve[Exp[-I*krad]*SeriesCoefficient[Series[eqmatrix[[4,4]],{e,0,1}]/.hubrep,0]==0,h[t]]]

perturbed

phimsol=Expand[Simplify[Solve[SeriesCoefficient[Series[eqmatrix[[4,4]],{e,0,1}],1]==0,phim[t]]][[1,1]]]/.dust

Eq. (0,i)

phimpsol=Simplify[Solve[SeriesCoefficient[Series[Sum[eqmatrix[[4,i]],{i,3}],{e,0,1}],1]==0,phim'[t]]][[1,1]]/.dust

Equations in customized form

We can now write the equations in a more readable form, also to eventually export them in Latex

custom={w[t]->0,w'[t]->0,cs2[t]->0,a'[t]->h[t]*a[t],a''[t]->h'[t]*a[t]+h[t]*a'[t],phi[t]->ϕ,phi'[t]->ϕ',phim[t]->Φ,phim'[t]->Φ',psim[t]->Φ,
psim'[t]->Φ',xi->ζ,dp[t]->φ,dp'[t]->φ',dp''[t]->φ'',rho[t]->ρ,a[t]->a,
d[t]->δ,d'[t]->δ',theta[t]->θ,theta'[t]->θ',phi''[t]->ϕ'',h[t]->h,h'[t]->h'};

dpsolw/.custom
Simplify[thetapsolw//.custom]
Simplify[dpsecsol//.custom]
Simplify[phimsol//.custom]
phimpsol//.custom