Axial and Polar perturbations with a source term given by an infalling particle

Exit[]

$Assumptions={r≥0,M>0,Element[L,Integers]};

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

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

Field’s equations

Tensor harmonic decomposition [For the tensor harmonics I use the definition of Sago et al. PRD 67, 10 (2003)]

(* Note that I write the metric perturbation by explicitng the frequency dependent part, so I will not have any time-deroivative in the resulting equations *)

This is the general form for the set of tensorial harmonics

a0={{Y[θ,ϕ],0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
a1=(I/(Sqrt[2])){{0,Y[θ,ϕ],0,0},{Y[θ,ϕ],0,0,0},{0,0,0,0},{0,0,0,0}};
a={{0,0,0,0},{0,Y[θ,ϕ],0,0},{0,0,0,0},{0,0,0,0}};
b0=((I r)/(((2l(l+1)))^(1/2))){{0,0,D[Y[θ,ϕ],θ],D[Y[θ,ϕ],ϕ]},{0,0,0,0},{D[Y[θ,ϕ],θ],0,0,0},{D[Y[θ,ϕ],ϕ],0,0,0}};
b=(r/(((2l(l+1)))^(1/2))){{0,0,0,0},{0,0,D[Y[θ,ϕ],θ],D[Y[θ,ϕ],ϕ]},{0,D[Y[θ,ϕ],θ],0,0},{0,D[Y[θ,ϕ],ϕ],0,0}};
c0=(r/(((2l(l+1)))^(1/2))){{0,0,Sθ[θ,ϕ],Sϕ[θ,ϕ]},{0,0,0,0},{Sθ[θ,ϕ],0,0,0},{Sϕ[θ,ϕ],0,0,0}};
c=((I r)/(((2l(l+1)))^(1/2))){{0,0,0,0},{0,0,Sθ[θ,ϕ],Sϕ[θ,ϕ]},{0,Sθ[θ,ϕ],0,0},{0,Sϕ[θ,ϕ],0,0}};
d=-((I (r^2))/(((2l(l+1)(l-1)(l+2)))^(1/2))){{0,0,0,0},{0,0,0,0},{0,0,-((X[θ,ϕ])/(Sin[θ])),Sin[θ]W[θ,ϕ]},
{0,0,Sin[θ]W[θ,ϕ],Sin[θ]X[θ,ϕ]}};
g=((r^2)/(Sqrt[2])){{0,0,0,0},{0,0,0,0},{0,0,Y[θ,ϕ],0},{0,0,0,((Sin[θ])^2)Y[θ,ϕ]}};
f=((r^2)/(((2l(l+1)(l-1)(l+2)))^(1/2))){{0,0,0,0},{0,0,0,0},{0,0,W[θ,ϕ],X[θ,ϕ]},{0,0,X[θ,ϕ],-((Sin[θ])^2)W[θ,ϕ]}};

a0//MatrixForm

Metric perturbations

h2[r_]=0;

η0[r_]=η1[r_]=G[r_]=0;

δgA=(-((((2l(l+1)))^(1/2))/r)h0[r]c0+I((((2l(l+1)))^(1/2))/r)h1[r]c+((((2l(l+1)(l-1)(l+2)))^(1/2))/(2(r^2)))h2[r]d)Exp[-I ω t];

δgP=(F[r]H0[r]a0-I Sqrt[2]H1[r]a1+((H2[r])/(F[r]))a-(I/r)(((2l(l+1)))^(1/2))η0[r]b0+(1/r)(((2l(l+1)))^(1/2))η1[r]b+
((((1/2)l(l+1)(l-1)(l+2)))^(1/2))G[r]f+(Sqrt[2]K[r]-((l(l+1))/(Sqrt[2]))G[r])g)Exp[-I ω t];

Compute the perturbed field’s equations

(* My bookeeping parameter for the perturbation is ϵ *)

Sϵ[F_]:=Series[F,{ϵ,0,1}]//Normal

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

ℊ0={{-(1-2M/r),0,0,0},{0,(1/((1-2M/r))),0,0},{0,0,(r^2),0},{0,0,0,(r^2)((Sin[θ])^2)}};

ℊ=ℊ0+ϵ δgP+ϵ δgA;

{Sθ[θ_,ϕ_],Sϕ[θ_,ϕ_]}:={((D[Y[θ,ϕ],ϕ])/(Sin[θ])),-Sin[θ]D[Y[θ,ϕ],θ]};

(* I compute the axial and polar components of the field's euqations. Due the the background spherical symmetry they decouple, so I treat them separately *)
ℰ=Sϵ[Ricci[ℊ,y]]-(1/2)Sϵ[ℊ RicciScalar[ℊ,y]];

(* Export the Files so I don't need to compute them every time *)
Export["field_equations_GR.m",ℰ]

Separate Axial and Polar components

Simplifications rules Y(θ,ϕ)

(* These are "a lot" of spherycal harmonics relations which deeply simplify the euqations. *)

rulesY={Derivative[2,0][Y][θ,ϕ]->-l(l+1)Y[θ,ϕ]-(1/((Sin[θ])^2))Derivative[0,2][Y][θ,ϕ]-Cot[θ]Derivative[1,0][Y][θ,ϕ],Derivative[3,0][Y][θ,ϕ]->l (1+l) Cot[θ] Y[θ,ϕ]-(1+l+(l^2)) Derivative[1,0][Y][θ,ϕ]+((Csc[θ])^2) (3 Cot[θ] Derivative[0,2][Y][θ,ϕ]+2 Derivative[1,0][Y][θ,ϕ]-Derivative[1,2][Y][θ,ϕ]),Derivative[4,0][Y][θ,ϕ]->l (1+l) (1+l+(l^2)-3 ((Csc[θ])^2)) Y[θ,ϕ]+Csc[θ] (Csc[θ] (7+l+(l^2)-11 ((Csc[θ])^2)) Derivative[0,2][Y][θ,ϕ]+((1+2 l (1+l)) Cos[θ]-6 Cot[θ] Csc[θ]) Derivative[1,0][Y][θ,ϕ]+Csc[θ] (5 Cot[θ] Derivative[1,2][Y][θ,ϕ]-Derivative[2,2][Y][θ,ϕ])),Derivative[2,1][Y][θ,ϕ]->-l (1+l) Derivative[0,1][Y][θ,ϕ]-((Csc[θ])^2) Derivative[0,3][Y][θ,ϕ]-Cot[θ] Derivative[1,1][Y][θ,ϕ],Derivative[2,2][Y][θ,ϕ]->-l (1+l) Derivative[0,2][Y][θ,ϕ]-((Csc[θ])^2) Derivative[0,4][Y][θ,ϕ]-Cot[θ] Derivative[1,2][Y][θ,ϕ]};

rulesY2={Derivative[1,1][Y][θ,ϕ]->(1/2) X[θ,ϕ]+Cot[θ] Derivative[0,1][Y][θ,ϕ],Derivative[0,2][Y][θ,ϕ]->-(1/2) Sin[θ] (Sin[θ] (W[θ,ϕ]+l (1+l) Y[θ,ϕ])+2 Cos[θ] Derivative[1,0][Y][θ,ϕ]),Derivative[2,0][Y][θ,ϕ]->(1/2) (W[θ,ϕ]-l (1+l) Y[θ,ϕ]),Derivative[3,0][Y][θ,ϕ]->(1/2) (-3 Cot[θ] W[θ,ϕ]-l (1+l) Cot[θ] Y[θ,ϕ]-2 (-2+l+(l^2)+((Csc[θ])^2)) Derivative[1,0][Y][θ,ϕ]-2 ((Csc[θ])^2) Derivative[1,2][Y][θ,ϕ]),Derivative[4,0][Y][θ,ϕ]->(1/2) (-(7+2 l (1+l)-11 ((Csc[θ])^2)) W[θ,ϕ]+2 ((Csc[θ])^4) Derivative[0,4][Y][θ,ϕ]+Cot[θ] (5 l (1+l) Cot[θ] Y[θ,ϕ]+2 (-6+5 ((Csc[θ])^2)) Derivative[1,0][Y][θ,ϕ]+12 ((Csc[θ])^2) Derivative[1,2][Y][θ,ϕ])),Derivative[2,1][Y][θ,ϕ]->-(1/2) Cot[θ] X[θ,ϕ]-(l+(l^2)+((Cot[θ])^2)) Derivative[0,1][Y][θ,ϕ]-((Csc[θ])^2) Derivative[0,3][Y][θ,ϕ],Derivative[2,2][Y][θ,ϕ]->-((Csc[θ])^2) Derivative[0,4][Y][θ,ϕ]+(1/2) l (1+l) Sin[θ] (Sin[θ] (W[θ,ϕ]+l (1+l) Y[θ,ϕ])+2 Cos[θ] Derivative[1,0][Y][θ,ϕ])-Cot[θ] Derivative[1,2][Y][θ,ϕ],Derivative[3,1][Y][θ,ϕ]->(1/4) (-2 (1+l+(l^2)-2 ((Csc[θ])^2)) X[θ,ϕ]+((Csc[θ])^3) ((7 Cos[θ]+Cos[3 θ]) Derivative[0,1][Y][θ,ϕ]+12 Cos[θ] Derivative[0,3][Y][θ,ϕ]-4 Sin[θ] Derivative[1,3][Y][θ,ϕ]))};

rulesXW={X[θ,ϕ]->2D[D[Y[θ,ϕ],θ],ϕ]-2D[Cot[θ]Y[θ,ϕ],ϕ],Derivative[1,0][X][θ,ϕ]->D[2D[D[Y[θ,ϕ],θ],ϕ]-2D[Cot[θ]Y[θ,ϕ],ϕ],θ],Derivative[0,1][X][θ,ϕ]->D[2D[D[Y[θ,ϕ],θ],ϕ]-2D[Cot[θ]Y[θ,ϕ],ϕ],ϕ],Derivative[1,1][X][θ,ϕ]->D[D[2D[D[Y[θ,ϕ],θ],ϕ]-2D[Cot[θ]Y[θ,ϕ],ϕ],ϕ],θ],
W[θ,ϕ]->D[Y[θ,ϕ],{θ,2}]-Cot[θ]D[Y[θ,ϕ],θ]-(1/((Sin[θ])^2))D[Y[θ,ϕ],{ϕ,2}],
Derivative[1,0][W][θ,ϕ]->D[D[Y[θ,ϕ],{θ,2}]-Cot[θ]D[Y[θ,ϕ],θ]-(1/((Sin[θ])^2))D[Y[θ,ϕ],{ϕ,2}],θ],Derivative[1,1][W][θ,ϕ]->D[D[D[Y[θ,ϕ],{θ,2}]-Cot[θ]D[Y[θ,ϕ],θ]-(1/((Sin[θ])^2))D[Y[θ,ϕ],{ϕ,2}],θ],ϕ],Derivative[0,1][W][θ,ϕ]->D[D[Y[θ,ϕ],{θ,2}]-Cot[θ]D[Y[θ,ϕ],θ]-(1/((Sin[θ])^2))D[Y[θ,ϕ],{ϕ,2}],ϕ]};

Einstein’s equations in the RW Gauge

(* In this section I isolate each component of the field's equations. Note that I work in the Regge-Wheeler Gauge, i.e. (Subscript[h, 2]) = (Subscript[η, 1]) = (Subscript[η, 0]) = G = 0. Moreover, I only keep the Einstein's equations components needed to derive the master equation for the axial modes, namely the Regge-Wheeler equation *)

(* If you compute everything correclty polar and axial compoenent do not speak each other, i.e. they are decoupled *)

h2[r_]=0;

ℰ13=Collect[-2((Coefficient[ℰ[[1,3]],ϵ,1])/((E^(-I t ω))Csc[θ]))//.Join[rulesY,rulesXW],{Y[θ,ϕ],Derivative[1,0][Y][θ,ϕ],Derivative[0,1][Y][θ,ϕ]},Simplify]

ℰ23=Collect[((Coefficient[ℰ[[2,3]],ϵ,1])/((E^(-I t ω))Csc[θ]))//.Join[rulesY,rulesXW],{Y[θ,ϕ],Derivative[1,0][Y][θ,ϕ],Derivative[0,1][Y][θ,ϕ],h0''[r],h0'[r],h0[r],h1''[r],h1'[r],h1[r],h2[r],ω},Simplify]

ℰ33=Collect[((Coefficient[ℰ[[3,3]],ϵ,1])/(E^(-I t ω)))//.rulesY2,{Y[θ,ϕ],Derivative[1,0][Y][θ,ϕ],Derivative[0,1][Y][θ,ϕ]},Simplify]

Axial perturbations ℓ ≥ 2

eq1A=Assuming[(l^2)+l==λ,Collect[Coefficient[-((F[r])^(-1))ℰ13,Derivative[0,1][Y][θ,ϕ]],{h0''[r],h1'[r],h1[r],𝒬0[r]},FullSimplify]]

eq2A=Assuming[(l^2)+l==λ,Collect[2F[r] Coefficient[-ℰ23,Derivative[0,1][Y][θ,ϕ]],{h0'[r],h0[r],h1[r],𝒬[r]},FullSimplify]]

eq3A=Assuming[(l^2)+l==λ,Collect[Coefficient[-2((ℰ33)/(Csc[θ])),X[θ,ϕ]],{h0[r],h1[r],h1'[r],𝒟[r]},Simplify]]

Rearrange to find the Regge-Wheeler equation

(* Use the Regge-Wheeler function R = h1(1-2M/r)/r *)
rulesRW={h1[r]->r((R[r])/(F[r])),h1'[r]->D[r((R[r])/(F[r])),r],h1''[r]->D[r((R[r])/(F[r])),{r,2}],F[r]->1-((2M)/r),F'[r]->D[1-((2M)/r),r]}//Simplify;

(* Derive h0 *)
sh0=Collect[Solve[(eq3A//.rulesRW)==0,h0[r]][[1,1,2]],{R'[r],R[r],𝒟[r]},Simplify]

(* Second order equation in R. Note that this is in the standard 'r' *)
RW=Collect[eq2A//.Join[rulesRW,{h0'[r]->D[sh0,r],h0[r]->sh0}],{R''[r],R'[r],R[r]},Simplify]

(* Finally change to the tortoise coordinate to find the canonical Schroedinger-like form of the equation *)

rulesT={R''[r]->-((F'[r])/((F[r])^2))R'[𝓇]+((R''[𝓇])/((F[r])^2)),R'[r]->((R'[𝓇])/(F[r])),F[r]->1-((2M)/r),F'[r]->D[1-((2M)/r),r]};

RWr=Collect[-F[r]/r RW//.rulesT,{R''[𝓇],R[r],𝒬[r],𝒟'[r],𝒟[r]},FullSimplify]