Slowly rotating solution for a BH with Schwarzschild background

Exit[]

$Assumptions={r≥0,M>0,Element[L,Integers],0<=θ≤π,0≤ϕ≤2π};

(* This set up the working directory of the notebook. *)
path=SetDirectory[NotebookDirectory[]];
SetDirectory[path];

(* This load the package to compute the needed tensor quantities *)
<<GR`

Hartle - Thorne metric [See ApJ 150:1005 (1967) for details]

ds=-F[r](1+2h)(dt^2)+(((1+2m))/(B[r]))(dr^2)+(r^2)(1+2k)(((dθ)^2)+((Sin[θ])^2)(((dϕ-Ω dt))^2));

h=(β^2)(h0[r]+h2[r] p2[θ]);
m=(β^2)(m0[r]+m2[r] p2[θ]);
k=(β^2)k2[r] p2[θ];
Ω=β ω[r];

S=Collect[ds,β,FullSimplify];

var={(dt^2),(dr^2),((dθ)^2),((dϕ)^2),β}

Coefficient[Collect[Series[S,{β,0,3}]/2,var,FullSimplify],dϕ,1]

Coefficient[Collect[Series[S,{β,0,3}]/2,var,FullSimplify],dϕ,1]

Metric variables

(* Here I define all the metric variables that will be used afterwards *)

(* Legendre polynomial for l=2 *)
p2[θ_]=LegendreP[2,Cos[θ]];

(* The metric tensor is defined in terms of a bookeeping parameter β which tracks the BH's spin order. Remember to make series expansion everytime you compute something. *)
βn=2;

g00=-F[r]-(β^2) (2 F[r] (h0[r]+h2[r] p2[θ])-(r^2) ((Sin[θ])^2) ((ω[r])^2));
grr=(1/(B[r]))+((2 (β^2) (m0[r]+m2[r] p2[θ]))/(B[r]));
gθθ=(r^2)+2 (r^2) (β^2) k2[r] p2[θ];
gϕϕ= (r^2) ((Sin[θ])^2)+2 (r^2) (β^2) k2[r] p2[θ] ((Sin[θ])^2);
gtϕ=-(r^2) β ((Sin[θ])^2) ω[r];

g[r_,θ_]={{g00,0,0,gtϕ},{0,grr,0,0},{0,0,gθθ,0},{gtϕ,0,0,gϕϕ}};

Ig[r_,θ_]=Collect[Series[Inverse[g[r,θ]],{β,0,βn}],β,Simplify];

y={t,r,θ,ϕ};

(* Christoffel Symbols *)
Γ=Collect[Series[Chr[g[r,θ],y],{β,0,βn}],β];

(* Riemann tensor fully covariant  *)
Riem=Collect[Series[Table[Sum[g[r,θ][[l1,l2]]Riemann[g[r,θ],y][[l1,l3,l4,l5]],{l1,1,4}],{l2,1,4},{l3,1,4},{l4,1,4},{l5,1,4}],{β,0,βn}],β];

(* Riemann tensor fully contracovariant  *)
Riemu=Collect[Series[Table[Sum[Ig[r,θ][[l1,l5]]Ig[r,θ][[l2,l6]]Ig[r,θ][[l3,l7]]Ig[r,θ][[l4,l8]]Riem[[l5,l6,l7,l8]],{l5,1,4},{l6,1,4},{l7,1,4},{l8,1,4}],{l1,1,4},{l2,1,4},{l3,1,4},{l4,1,4}],{β,0,βn}],β];

(* Kretchmann scalar  *)

𝒦=Series[Sum[Riemu[[l1,l2,l3,l4]]Riem[[l1,l2,l3,l4]],{l1,1,4},{l2,1,4},{l3,1,4},{l4,1,4}],{β,0,βn}]//Normal;

Field’s equations

(* Einstein equations expanded at the desired order in the spin parameter *)

(* Ricci tensor (Subscript[R, αβ]) *)
Ric=Series[Ricci[g[r,θ],y],{β,0,βn}]//Normal;

(* Ricci scalar *)
RicS=Series[RicciScalar[g[r,θ],y],{β,0,βn}]//Normal;

ℰdd=Ric-(1/2)Collect[Series[g[r,θ] RicS,{β,0,βn}],β]//Normal;

(* I prefer to work with the mixed components (Subscript[(G^α), β]) *)
ℰ=Collect[Series[Table[Sum[Ig[r,θ][[l1,l3]]ℰdd[[l3,l2]],{l3,1,4}],{l1,1,4},{l2,1,4}],{β,0,βn}],β]//Normal;

Background solution

DSolve[Coefficient[ℰ[[1,1]],β,0]==0,B[r],r]
B[r_]=%[[1,1,2]]//.C[1]->-2M

(* The costant is fixed by requring asymtotic flatness *)

DSolve[Coefficient[ℰ[[2,2]],β,0]==0,F[r],r]//Simplify
F[r_]=%[[1,1,2]]//.C[1]->1

𝒪(χ) equations

(* Find the spin linear contribution in the Einstein's equations. Note that (Subscript[ℰ, 14]) is the only non zero compoenent at this order *)

ℰ141=Series[Coefficient[ℰ[[1,4]],β,1],{ϵ,0,1}]//Simplify

(* We integrate the equation and fix the constants of integration by asking the metric to be asymptotically flat, i.e. C[2] -> 0,and to satisfy the far field limit of an isolated source with (Subscript[g, tϕ]) = -2J/r *)

DSolve[Coefficient[ℰ141,ϵ,0]==0,ω[r],r]

ω[r_]=DSolve[Coefficient[ℰ141,ϵ,0]==0,ω[r],r][[1,1,2]]//.{C[1]->-6J,C[2]->0}

𝒪((χ^2)) equations

Monopole component

(* Here I integrate over Sin[θ] dθ using the ortogonality properties of Legendre polinomials to separate the contributions from l=0, l=2. For example, top get only the l=0 I integrate ∫(Subscript[G, tt]) Sin[θ] dt.  For the l=2 I integrate ∫(Subscript[G, tt]) (Subscript[P, 2])[Sin[θ]]Sin[θ] *)

var={h0[r],m0[r],h2[r],m2[r],k2[r],h0'[r],m0'[r],h2'[r],m2'[r],k2'[r],h0''[r],m0''[r],h2''[r],m2''[r],k2''[r]};

eq011=Integrate[Collect[Coefficient[ℰ[[1,1]],β,2],var,Simplify]Sin[θ],{θ,0,π}]

eq022=Integrate[Collect[Coefficient[ℰ[[2,2]],β,2],var,Simplify]Sin[θ],{θ,0,π}]

(* From the Series at infinity you see that C[1] -> 0 *)
DSolve[eq011==0,m0[r],r][[1,1,2]];
Series[%,{r,∞,1}]
m0[r_]=DSolve[eq011==0,m0[r],r][[1,1,2]]//.C[1]->0

DSolve[Coefficient[eq022,ϵ,0]==0,h0[r],r][[1,1,2]];
Series[%,{r,∞,1}]
h0[r_]=DSolve[Coefficient[eq022,ϵ,0]==0,h0[r],r][[1,1,2]]//.C[1]->0

Quadrupole component

(* First isolate the quadrupolar component by integrating over the Legendre polynomial *)

eq211=Integrate[Collect[Coefficient[ℰ[[1,1]],β,2],var,Simplify]LegendreP[2,Cos[θ]]Sin[θ],{θ,0,π}]

eq222=Integrate[Collect[Coefficient[ℰ[[2,2]],β,2],var,Simplify]LegendreP[2,Cos[θ]]Sin[θ],{θ,0,π}]

eq23344=Integrate[(Collect[Coefficient[ℰ[[3,3]],β,2],var,Simplify]+Collect[Coefficient[ℰ[[4,4]],β,2],var,Simplify])LegendreP[2,Cos[θ]]Sin[θ],{θ,0,π}]

(* Herafter I "massage" the solution to find the easiest form to compute (h2, m2, k2) *)

sol=Solve[{eq211==0,eq222==0,eq23344==0},{k2''[r],h2''[r],m2[r]}][[1]]//Simplify

m2[r_]=sol[[3,2]]

sol2=Solve[{eq211==0,eq23344==0},{k2'[r],h2''[r]}][[1]]//Simplify

Eq1=eq211//.{sol2[[1]],k2''[r]->D[sol2[[1,2]],r]}//FullSimplify

Eq2=eq222//.{sol2[[1]],k2''[r]->D[sol2[[1,2]],r]}//FullSimplify

Eq3=eq23344//.{sol2[[1]],k2''[r]->D[sol2[[1,2]],r]}//FullSimplify

solh2k2=DSolve[{Simplify[Eq1]==0,k2'[r]==sol2[[1,2]]},{k2[r],h2[r]},r][[1]]

(* I fix the constant C[2] by requiring no logaritmic divergence *)
Collect[solh2k2[[1,2]],{Log[-2 M+r],Log[r]},Simplify]
Solve[((3 (2 M-r) r (4 (J^2) (M^2)-C[2]) Log[r])/(32 (M^8)))==0,C[2]]

(* C[1] again set by asymptotic flatness *)
Series[solh2k2[[1,2]],{r,∞,1}]

h2[r_]=solh2k2[[1,2]]//.{C[1]->0,C[2]->4 (J^2) (M^2)}//Simplify

k2[r_]=solh2k2[[2,2]]//.{C[1]->0,C[2]->4 (J^2) (M^2)}//Simplify

m2[r]//Simplify

Kretchmann scalar

Collect[Series[𝒦,{β,0,2}],β,Simplify]

Total metric

g00=-F[r]-(β^2) (2 F[r] (h0[r]+h2[r] p2[θ])-(r^2) ((Sin[θ])^2) ((ω[r])^2));
grr=(1/(B[r]))+((2 (β^2) (m0[r]+m2[r] p2[θ]))/(B[r]));
gθθ=(r^2)+2 (r^2) (β^2) k2[r] p2[θ];
gϕϕ= (r^2) ((Sin[θ])^2)+2 (r^2) (β^2) k2[r] p2[θ] ((Sin[θ])^2);
gtϕ=-(r^2) β ((Sin[θ])^2) ω[r];

g00

gϕϕ

g00

(* Export the final metric at the second order in the spin *)
Export["sloly_metric.m",Collect[Normal[Series[{g00,grr,gθθ,gϕϕ,gtϕ}//.{J->(M^2)χ,β->1},{χ,0,2}]],{χ,r},Simplify]]

Geodesic properties at the second order

(* Import the metric and define some auxiliary function to extrac the spin components order by order *)

p2[θ_]=LegendreP[2,Cos[θ]];

{g00,grr,gθθ,gϕϕ,gtϕ}=Import["sloly_metric.m"];

ℊ={{g00,0,0,gtϕ},{0,grr,0,0},{0,0,gθθ,0},{gtϕ,0,0,gϕϕ}};

CC[F_,n1_]:=Coefficient[F,χ,n1](χ^n1)

(* Now I will solve all the relevant equations order by order *)

Horizon

HorEq=Series[g00 gϕϕ-((gtϕ)^2),{χ,0,2}]//Normal;

rH0=Solve[CC[HorEq,0]==0,r][[2,1,2]]

rH2=Series[Solve[CC[HorEq,0]+(CC[HorEq,2]//.r->rH0)==0,r][[2,1,2]],{χ,0,2}]//Normal

(* Check that the equation is solved by the expression we found *)

Series[HorEq//.r->rH2,{χ,0,2}]//Normal

ISCO

(* Here I define the effective potential, the energy and angular momentum, and the equatorial angular velocity. 
Note that the potential is given as a function of the parameter s:  
s=0 for massless particles and s=1 for massive particles. *)

{𝒱[s_],ωϕ}=With[{gtt=ℊ[[1,1]],grr=ℊ[[2,2]],gtϕ=ℊ[[1,4]],gϕϕ=ℊ[[4,4]]},{(1/grr)((((e^2)gϕϕ+2e ℒ gtϕ+(ℒ^2)gtt)/(((gtϕ)^2)-gtt gϕϕ))-s),((-D[gtϕ,r]+Sqrt[((D[gtϕ,r])^2)-D[gtt,r]D[gϕϕ,r]])/(D[gϕϕ,r]))}]//.θ->π/2;

{ℯ,ℓ}=With[{gtt=ℊ[[1,1]],grr=ℊ[[2,2]],gtϕ=ℊ[[1,4]],gϕϕ=ℊ[[4,4]]},{Normal[Series[-((gtt+gtϕ ωϕ)/(Sqrt[-gtt-2gtϕ  ωϕ-gϕϕ ((ωϕ)^2)])),{χ,0,2}]],Normal[Series[((gtϕ+gϕϕ ωϕ)/(Sqrt[-gtt-2gtϕ  ωϕ-gϕϕ ((ωϕ)^2)])),{χ,0,2}]]}]//.θ->π/2;

dd𝒱=Series[D[𝒱[1]//.θ->π/2,{r,2}]//.{e->ℯ,ℒ->ℓ,θ->π/2},{χ,0,2}]//Normal;

rI0=Solve[CC[dd𝒱,0]==0,r][[1,1,2]]

fac[r_]=((Coefficient[Simplify[CC[dd𝒱,0,0]],r-6M])^(-1));

rI1=Solve[fac[r]CC[dd𝒱,0]+(fac[r]CC[dd𝒱,1]//.r->rI0)==0,r][[1,1,2]]//Simplify

rI2=Series[Solve[fac[r]CC[dd𝒱,0]+(fac[r](CC[dd𝒱,1]+CC[dd𝒱,2])//.r->rI1)==0,r][[1,1,2]],{χ,0,2}]//Normal

(* Check that the equation is solved by the expression we found *)

Series[fac[r]dd𝒱//.r->rI2,{χ,0,2}]//Normal

Ergo-region

rE0=Solve[CC[g00,0]==0,r][[1,1,2]]

rE2=Collect[Series[Solve[CC[g00,0]+(Simplify[CC[g00,1]+CC[g00,2]]//.r->rE0)==0,r][[1,1,2]],{χ,0,2}],χ,Simplify]

(* Check that the equation is solved by the expression we found *)

Series[g00//.r->rE2,{ϵ,0,1},{χ,0,3}]//Normal//FullSimplify

Angular velocity of the horizon and moment of Inertia

rH2

gtϕ

gϕϕ

Ωh=Series[-((gtϕ)/(gϕϕ))//.{r->rH2},{χ,0,2}];

Ωh

ℐ=Series[(((M^2)χ)/(Assuming[M>0,Ωh (M^3)])),{χ,0,2}];

Collect[Ωh,χ,Simplify]//.θ->π/2

ℐ