GR calculations

L.A.  2024

Definitions

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 now the static radial  metric.

(* r, theta, phi, t *)

Clear[metric];

metric:=DiagonalMatrix[{Exp[2*lam[r]],r^2,r^2*Sin[theta]^2,-Exp[2*phi[r]]}]

coord={r,theta,phi,t};

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.

Raise[metric_,s_]:=Inverse[metric].s;
Lower[metric_,s_]:=metric.s;

Ordinary derivative and divergence can be obtained with the commands Grad and Div, respectively. E.g.

vector={f1[r,t],f2[r,t],f3[r,t],f4[r,t]};

(* (Subscript[V,  ,α]^α)  *) 
Div[vector,coord]

(* (Subscript[V,  ,β]^α)  *) 
Grad[vector,coord]

Since the multiplication is always between the last index of the first tensor and the first index of the second tensor, if we need to use another index we have to Transpose the tensor first. We can then define the Christoffel symbols (first index up):

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

For instance, we have

Simplify[Christoffel[metric,coord]]//MatrixForm

It is useful also to derive the contracted Christoffel symbol

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

Simplify[ChrContracted[metric,coord]]

We can now evaluate the Ricci tensor (index down)

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  Ricci tensor is

Simplify[Ricci[metric,coord]]

and the Trace is

ricciscalar=Simplify[Tr[Raise[metric,Ricci[metric,coord]]]]

Components of Einstein tensor

(* g00 component *)

g00=Simplify[Ricci[metric,coord][[4,4]]-metric[[4,4]]*ricciscalar/2]

(* grr component *)

grr=Simplify[Ricci[metric,coord][[1,1]]-metric[[1,1]]*ricciscalar/2]

(* gthetatheta component *)

gtt=Simplify[Ricci[metric,coord][[2,2]]-metric[[2,2]]*ricciscalar/2]

(* gphiphi component *)

gpp=Simplify[Ricci[metric,coord][[3,3]]-metric[[3,3]]*ricciscalar/2]

Schwarzschild metric. Kretschmann scalar.

The Kretschmann scalar is defined as

In vacuum, the Ricci tensor for the Schwarzschild metric vanished, so the space - time is Ricci flat . But since the Riemann tensor does not vanish, the space - time is not globally flat . Let’s evaluate then the Riemann tensor and the Kretschmann scalar.

We can now adopt the Schwarszschild metric in vacuum

(* r, theta, phi, t *)

lam[r_]=-1/2*Log[(1-2*M/r)];
phi[r_]=1/2*Log[(1-2*M/r)];

and obtain the Riemann tensor

chrd[metric_,var_]=Simplify[Grad[Christoffel[metric,var],var]];

(* Riemann all indexes low *)

riemann[metric_,var_]:=Lower[metric,Transpose[chrd[metric,var],3<->4]-chrd[metric,var]-
Table[Sum[Christoffel[metric,var][[alpha,nu,sigma]]*Christoffel[metric,var][[sigma,beta,mu]]-Christoffel[metric,var][[alpha,sigma,mu]]*Christoffel[metric,var][[sigma,nu,beta]],{sigma,4}],{alpha,4},{beta,4},{mu,4},{nu,4}]]

riedown=Simplify[riemann[metric,coord]];

(* Inverse metric *)
metricup=Simplify[Inverse[metric]]

(* Riemann all indexes up *)

riemannup[metric_,var_]:=Table[Sum[metricup[[iup,i]]*metricup[[jup,j]]*metricup[[kup,k]]*metricup[[mup,m]]*riedown[[i,j,k,m]],{i,4},{j,4},{k,4},{m,4}],{iup,4},{jup,4},{kup,4},{mup,4}]

rieup=Simplify[riemannup[metric,coord]];

Finally :

Simplify[Sum[riedown[[i,j,k,m]]*rieup[[i,j,k,m]],{i,4},{j,4},{k,4},{m,4}]]

which is indeed not zero and show the singularity at the origin .

FLRW metric (flat space)

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

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

Components of Einstein tensor

(* g00 component *)

g00=Simplify[Ricci[metric,coord][[4,4]]-metric[[4,4]]*ricciscalar/2]

(* grr component *)

grr=Simplify[Ricci[metric,coord][[1,1]]-metric[[1,1]]*ricciscalar/2]

(* gthetatheta component *)

grr=Simplify[Ricci[metric,coord][[2,2]]-metric[[2,2]]*ricciscalar/2]

(* gphiphi component *)

grr=Simplify[Ricci[metric,coord][[3,3]]-metric[[3,3]]*ricciscalar/2]

4 - velocity and EMT (mixed form)

Clear[uvup];

uvup={0,0,0,1/a[t]}

Clear[uv];uv=Lower[metric,uvup];

Clear[prod];prod=Outer[Times,uvup,uv]

tem=Simplify[(rho[t]*(1+w[t])*prod+IdentityMatrix[4]*rho[t]*(w[t]))];
tem//TableForm

temdo=Lower[metric,tem];

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

chr=Christoffel[metric,coord];

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
Solve[emcb[[4]]==0,rho'[t]][[1,1]]

Einstein equations (index down)

ei=Ricci[metric,coord];

eis=Tr[Raise[metric,ei]];


eqmatrix=ei-metric*eis/2-8*Pi*temdo;

ah={a'[t]->h[t]*a[t],a''[t]->h'[t]*a[t]+h[t]*a'[t]};

g00=Simplify[eqmatrix[[4,4]]/.ah]

g11=Simplify[eqmatrix[[1,1]]//.ah]

Simplify[Solve[g00==0,h[t]]]

Simplify[Solve[g11==0,h'[t]]]