Basic Mathematica For Black Holes

Mathematica file for only calculating metric and tensor quantities. 

Raising and lowering indices is a pain in the a$$ but this is a good way to start learning how to do numerical GR.

It is usually better to write a consistent package that does so (I do have such a package but it is built on personal idiosyncrasies wrt choice of convention so I have decided not to share that with this course cohort.)

4D Schwarzschild (Coordinates = t,r,θ,ϕ)

Definitions

ClearAll["Global`*"]
dim=4;
LCT= LeviCivitaTensor[dim];
g = {{-(1-2M/r),0,0,0},{0,(1-2M/r)^(-1),0,0},{0, 0, r^2 ,0},{0,0,0, r^2 Sin[θ]^2}};
G = Inverse[g]//Simplify;
detg = Det[g];
(Subscript[x, 1])=t;
(Subscript[x, 2])=r;
(Subscript[x, 3
])= θ;
(Subscript[x, 4]) = ϕ;

perm[dim_,n_]:=Permutations[Sort[Flatten[(ConstantArray[a,n])/.a->Range[dim]]],{n}];
(*$Assumptions={r>0};*)

G//Simplify

g.G//Simplify

Christoffel Symbols

Γ[α_, μ_, ν_] := FullSimplify[(∑_(β=1)^dim)(1/2)G[[α,β]](D[g[[β,μ]],(Subscript[x, ν])]+D[g[[β,ν]],(Subscript[x, μ])]-D[g[[μ,ν]],(Subscript[x, β])])]

Column[Select[({Subsuperscript[Γ, #[[2]]*10 + #[[
         3]], #[[1]]], Apply[Γ, #]} & /@ 
    perm[dim, 3]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor (First index up, next three down)

R[α_, β_, γ_, δ_] := R[α,β,γ,δ]=FullSimplify[
 D[Γ[α, β, δ], Subscript[
   x, γ]] - 
  D[Γ[α, γ, β], Subscript[
   x, δ]] + (∑_(μ=1)^dim)(Γ[μ,β,δ]Γ[α,μ,γ]-Γ[μ,β,γ]Γ[α,μ,δ])]
Column[Select[({Subsuperscript[R,#[[2]]*100+#[[3]]*10+#[[4]],#[[1]]],Apply[R,#]}&/@perm[dim,4]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

Ricci Tensor

Ricci[α_, β_] := FullSimplify[(∑_(γ=1)^dim)R[γ,α,γ,β]]

Column[Select[({Subscript[Ricci, #[[1]]*10 + #[[2]]], 
      Apply[Ricci, #]} & /@ 
    perm[dim, 2]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor (all lower indices)

lowerRiemann[a_,b_,c_,d_] := lowerRiemann[a,b,c,d] = FullSimplify[(∑_(x=1)^dim)R[x,b,c,d] g[[a,x]]]

Column[Select[({Subscript[R, #[[1]]*1000+#[[2]]*100+#[[3]]*10+#[[4]]], 
      Apply[lowerRiemann, #]} & /@ 
    perm[dim, 4]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor (All upper)

ClearAll[upperRiemann]

upperRiemann[a_,b_,c_,d_] := upperRiemann[a,b,c,d] = FullSimplify[(∑_(k=1)^dim)(∑_(l=1)^dim)(∑_(m=1)^dim)(∑_(n=1)^dim)lowerRiemann[k,l,m,n] G[[k,a]] G[[l,b]] G[[m,c]] G[[n,d]]]

Column[Select[({Superscript[R, #[[1]]*1000+#[[2]]*100+#[[3]]*10+#[[4]]], 
      Apply[upperRiemann, #]} & /@ 
    perm[dim, 4]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Ricci Scalar (= 0 for vacuum solution)

RicciS = FullSimplify[(∑_(α=1)^dim)(∑_(β=1)^dim)G[[α,β]]Ricci[β,α]]

Kretschmann Scalar

ClearAll[KretschmanS]

KretschmanS = Simplify[Sum[lowerRiemann[a,b,c,d]*upperRiemann[a,b,c,d],{a,dim},{b,dim},{c,dim},{d,dim}], M≥ 0]

Einstein equations

Λ = 0; (* Cosmological constant is zero *)
Ein[α_,β_]:=FullSimplify[Ricci[α,β]-1/2 g[[α,β]]*RicciS+g[[α,β]] Λ];
Column[({Subscript[Ein,#[[1]]*10+#[[2]]],Apply[Ein,#]}&/@perm[dim,2])]

4D Kerr (Coordinates = t,r,θ,ϕ)

Definitions

ClearAll["Global`*"]
dim=4;
LCT= LeviCivitaTensor[dim];
xmu={t,r,θ,ϕ};
(Subscript[x, 1])=t;
(Subscript[x, 2])=r;
(Subscript[x, 3
])= θ;
(Subscript[x, 4]) = ϕ;
perm[dim_,n_]:=Permutations[Sort[Flatten[(ConstantArray[a,n])/.a->Range[dim]]],{n}];
(*$Assumptions={r>0};*)

a = J/M;
e = 0;
Σ = r^2 + a^2 Cos[θ]^2;
Δ = r^2 - 2 M r + a^2 + e^2;

met =-(((Δ - (a^(2 ))((Sin[θ])^2))/Σ)) d[t]^2 + (Σ/Δ) d[r]^2- 2 a Sin[θ]^2((((r^2) + (a^2)- Δ)/Σ)) d[t]d[ϕ] + Σ d[θ]^2 + (((((((r^(2 ))+ (a^2)))^2) - Δ (a^(2 ))((Sin[θ])^2))/Σ))Sin[θ]^2  d[ϕ]^2;
g=Table[If[a==b,1,1/2]Coefficient[met,d[xmu[[a]]]d[xmu[[b]]]],{a,1,dim},{b,1,dim}];
G= Inverse[g]//Simplify;
sqrtdetmet = Normal[Series[Sqrt[-Det[met]],{ϵ,0,1}]]//Simplify;

lct  =sqrtdetmet* LeviCivitaTensor[dim];
lctup = ((-1)/sqrtdetmet)LeviCivitaTensor[dim];
detg = Det[g];

g//Simplify

g.G//Simplify

Christoffel Symbols

Γ[α_, μ_, ν_] := FullSimplify[(∑_(β=1)^dim)(1/2)G[[α,β]](D[g[[β,μ]],(Subscript[x, ν])]+D[g[[β,ν]],(Subscript[x, μ])]-D[g[[μ,ν]],(Subscript[x, β])])]

Column[Select[({Subsuperscript[Γ, #[[2]]*10 + #[[
         3]], #[[1]]], Apply[Γ, #]} & /@ 
    perm[dim, 3]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor (First index up, next three down)

R[α_, β_, γ_, δ_] := R[α,β,γ,δ]=FullSimplify[
 D[Γ[α, β, δ], Subscript[
   x, γ]] - 
  D[Γ[α, γ, β], Subscript[
   x, δ]] + (∑_(μ=1)^dim)(Γ[μ,β,δ]Γ[α,μ,γ]-Γ[μ,β,γ]Γ[α,μ,δ])]
Column[Select[({Subsuperscript[R,#[[2]]*100+#[[3]]*10+#[[4]],#[[1]]],Apply[R,#]}&/@perm[dim,4]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

Ricci Tensor

Ricci[α_, β_] := FullSimplify[(∑_(γ=1)^dim)R[γ,α,γ,β]]

Column[Select[({Subscript[Ricci, #[[1]]*10 + #[[2]]], 
      Apply[Ricci, #]} & /@ 
    perm[dim, 2]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor (all lower indices)

lowerRiemann[a_,b_,c_,d_] := lowerRiemann[a,b,c,d] = FullSimplify[(∑_(x=1)^dim)R[x,b,c,d] g[[a,x]]]

Column[Select[({Subscript[R, #[[1]]*1000+#[[2]]*100+#[[3]]*10+#[[4]]], 
      Apply[lowerRiemann, #]} & /@ 
    perm[dim, 4]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor (All upper)

ClearAll[upperRiemann]

upperRiemann[a_,b_,c_,d_] := upperRiemann[a,b,c,d] = FullSimplify[(∑_(k=1)^dim)(∑_(l=1)^dim)(∑_(m=1)^dim)(∑_(n=1)^dim)lowerRiemann[k,l,m,n] G[[k,a]] G[[l,b]] G[[m,c]] G[[n,d]]]

Column[Select[({Superscript[R, #[[1]]*1000+#[[2]]*100+#[[3]]*10+#[[4]]], 
      Apply[upperRiemann, #]} & /@ 
    perm[dim, 4]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Ricci Scalar (= 0 for vacuum solution)

RicciS = FullSimplify[(∑_(α=1)^dim)(∑_(β=1)^dim)G[[α,β]]Ricci[β,α]]

Kretschmann Scalar

ClearAll[KretschmanS]

KretschmanS = Sum[lowerRiemann[a,b,c,d]*upperRiemann[a,b,c,d],{a,dim},{b,dim},{c,dim},{d,dim}]

FullSimplify[KretschmanS, M≥ 0]

Einstein equations

Λ = 0; (* Cosmological constant is zero *)
Ein[α_,β_]:=FullSimplify[Ricci[α,β]-1/2 g[[α,β]]*RicciS+g[[α,β]] Λ];
Column[({Subscript[Ein,#[[1]]*10+#[[2]]],Apply[Ein,#]}&/@perm[dim,2])]

3D BTZ (Schwarzschild coordinates = t,r,θ)

Definitions

ClearAll["Global`*"]
dim=3;
$Assumptions={r>0, t>0, M>0};
LCT= LeviCivitaTensor[dim];
J = 0;
Np = Sqrt[-M + r^2/l^2 + J^2/(4*r^2)];
Nf = - J/(2r^2);
g = {{-Np^2 + Nf^2,0,r^2* Nf},{0,1/Np^2,0},{r^2 * Nf,0 , r^2}};
G = Inverse[g];
detg = Det[g];
(Subscript[x, 1])=t;
(Subscript[x, 2])=r;
(Subscript[x, 3
])= ϕ;
perm[dim_,n_]:=Permutations[Sort[Flatten[(ConstantArray[a,n])/.a->Range[dim]]],{n}];
$Assumptions={z>0};

Christoffel Symbols

Γ[α_, μ_, ν_] := FullSimplify[(∑_(β=1)^dim)(1/2)G[[α,β]](D[g[[β,μ]],(Subscript[x, ν])]+D[g[[β,ν]],(Subscript[x, μ])]-D[g[[μ,ν]],(Subscript[x, β])])]

Column[Select[({Subsuperscript[Γ, #[[2]]*10 + #[[
         3]], #[[1]]], Apply[Γ, #]} & /@ 
    perm[dim, 3]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor

R[α_, β_, γ_, δ_] :=FullSimplify[
 D[Γ[α, β, δ], Subscript[
   x, γ]] - 
  D[Γ[α, γ, β], Subscript[
   x, δ]] + (∑_(μ=1)^dim)(Γ[μ,β,δ]Γ[α,μ,γ]-Γ[μ,β,γ]Γ[α,μ,δ])]
Column[Select[({Subsuperscript[R,#[[2]]*100+#[[3]]*10+#[[4]],#[[1]]],Apply[R,#]}&/@perm[dim,4]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

Ricci Tensor

Ricci[α_, β_] := FullSimplify[(∑_(γ=1)^dim)R[γ,α,γ,β]]

Column[Select[({Subscript[Ricci, #[[1]]*10 + #[[2]]], 
      Apply[Ricci, #]} & /@ 
    perm[dim, 2]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Ricci Scalar

RicciS = FullSimplify[(∑_(α=1)^dim)(∑_(β=1)^dim)G[[α,β]]Ricci[β,α]]

Λ = -1/l^2;
Ein[α_,β_]:=FullSimplify[Ricci[α,β]-1/2 g[[α,β]]*RicciS+g[[α,β]] Λ];
Column[Select[({Subscript[Ein,#[[1]]*10+#[[2]]],Apply[Ein,#]}&/@perm[dim,2]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

Riccimod[σ_,ν_]:=FullSimplify[Ricci[σ,ν]-1/4*RicciS*g[[σ,ν]]]

Cotton[μ_, ν_] := FullSimplify[(∑_(α=1)^dim)(∑_(κ=1)^dim)(∑_(σ=1)^dim)(1/detg*g[[μ,α]]*LCT[[α,κ,σ]]*(D[Riccimod[σ,ν],(Subscript[x, κ])]-((∑_(δ=1)^dim)(Riccimod[δ,ν]*Γ[δ,σ,κ]))-((∑_(δ=1)^dim)(Riccimod[σ,δ]*Γ[δ,ν,κ]))))]
Column[Select[({Subscript[Cotton,#[[1]]*10+#[[2]]],Apply[Cotton,#]}&/@perm[dim,2]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

FullSimplify[(∑_(α=1)^dim)(∑_(β=1)^dim)(G[[α,β]]*Cotton[α,β])](* Sanity Check! = 0*)

TMG[α_,β_]:=FullSimplify[Ein[α,β]+1/m*Cotton[α,β]];(*sign error in Cotton tensor?*)Column[Select[({Subscript[TMG,#[[1]]*10+#[[2]]],Apply[TMG,#]}&/@perm[dim,2]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

3D BTZ (EF coordinates = w, r, ϕ)

Definitions

ClearAll["Global`*"]
dim=3;
LCT= LeviCivitaTensor[dim];
g = {{-r^2/l^2 + m, 1, -j/2},{1,0,0},{-j/2, 0, r^2}};
G = Inverse[g];
detg = Det[g];
(Subscript[x, 1])=w;
(Subscript[x, 2])=r;
(Subscript[x, 3
])= ϕ ;
perm[dim_,n_]:=Permutations[Sort[Flatten[(ConstantArray[a,n])/.a->Range[dim]]],{n}];
$Assumptions={z>0};

Christoffel Symbols

Γ[α_, μ_, ν_] := FullSimplify[(∑_(β=1)^dim)(1/2)G[[α,β]](D[g[[β,μ]],(Subscript[x, ν])]+D[g[[β,ν]],(Subscript[x, μ])]-D[g[[μ,ν]],(Subscript[x, β])])]

Column[Select[({Subsuperscript[Γ, #[[2]]*10 + #[[
         3]], #[[1]]], Apply[Γ, #]} & /@ 
    perm[dim, 3]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Riemann Tensor

R[α_, β_, γ_, δ_] :=FullSimplify[
 D[Γ[α, β, δ], Subscript[
   x, γ]] - 
  D[Γ[α, γ, β], Subscript[
   x, δ]] + (∑_(μ=1)^dim)(Γ[μ,β,δ]Γ[α,μ,γ]-Γ[μ,β,γ]Γ[α,μ,δ])]
Column[Select[({Subsuperscript[R,#[[2]]*100+#[[3]]*10+#[[4]],#[[1]]],Apply[R,#]}&/@perm[dim,4]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

Ricci Tensor

Ricci[α_, β_] := FullSimplify[(∑_(γ=1)^dim)R[γ,α,γ,β]]

Column[Select[({Subscript[Ricci, #[[1]]*10 + #[[2]]], 
      Apply[Ricci, #]} & /@ 
    perm[dim, 2]), (! NumericQ[#[[2]]] || #[[2]] != 0) &]]

Ricci Scalar

RicciS = FullSimplify[(∑_(α=1)^dim)(∑_(β=1)^dim)G[[α,β]]Ricci[β,α]]

Λ = -1/l^2;
Ein[α_,β_]:=FullSimplify[Ricci[α,β]-1/2 g[[α,β]]*RicciS+g[[α,β]] Λ];
Column[Select[({Subscript[Ein,#[[1]]*10+#[[2]]],Apply[Ein,#]}&/@perm[dim,2]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

Riccimod[σ_,ν_]:=FullSimplify[Ricci[σ,ν]-1/4*RicciS*g[[σ,ν]]]

Cotton[μ_, ν_] := FullSimplify[(∑_(α=1)^dim)(∑_(κ=1)^dim)(∑_(σ=1)^dim)(1/detg*g[[μ,α]]*LCT[[α,κ,σ]]*(D[Riccimod[σ,ν],(Subscript[x, κ])]-((∑_(δ=1)^dim)(Riccimod[δ,ν]*Γ[δ,σ,κ]))-((∑_(δ=1)^dim)(Riccimod[σ,δ]*Γ[δ,ν,κ]))))]
Column[Select[({Subscript[Cotton,#[[1]]*10+#[[2]]],Apply[Cotton,#]}&/@perm[dim,2]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]

FullSimplify[(∑_(α=1)^dim)(∑_(β=1)^dim)(G[[α,β]]*Cotton[α,β])](* Sanity Check! = 0*)

TMG[α_,β_]:=FullSimplify[Ein[α,β]+1/m*Cotton[α,β]];(*sign error in Cotton tensor?*)Column[Select[({Subscript[TMG,#[[1]]*10+#[[2]]],Apply[TMG,#]}&/@perm[dim,2]),(!NumericQ[#[[2]]]||#[[2]]≠0)&]]