Regge-Wheeler and Zerilli equations

References:
	[1] T. Regge and  J.A. Wheeler, PR108:1063
	[2] F.J. Zerilli, PRL24:737

We derive the Regge-Wheeler [1] and Zerilli [2] equations for the linear pertubations of spherically symmetric vacuum spacetimes using xTensor and xCoba. xTensor can be used for abstract tensor calculations and xCoba for explicitly computing components. The ansatz for the metric in Regge-Wheeler coordinates is directly inserted into the Einstein equations and all components are calculated explicitly using xCoba. Notice that it would be faster to use xPert in an intermediate step for perturbations (important for more complicated situations, like higher order perturbations, gauge invariant formalisms etc).

This notebook requires the xAct package: http://www.xact.es/

Created by: Sayan Chakrabarti, Jan Steinhoff, Jorge Rocha

Basic Definitions

load xTensor and xCoba packages:

<<xAct`xCoba`;

remove dollar symbols from indices:

$PrePrint=ScreenDollarIndices;

be less verbose:

$DefInfoQ=False;
$CVVerbose=False;

fix conventions:

$CovDFormat="Prefix";
$epsilonSign=1;
$RicciSign=1;
$RiemannSign=1;
(*$TorsionSign=1;
$ExtrinsicKSign=-1;
$AccelerationSign=-1;*)

define manifold and metric:

DefManifold[H,4,{a,b,c,d,a1,a2,a3,a4,b1,b2,b3,b4}];
DefMetric[-1,g[-a,-b],CD,{";","▽"}];

check Riemann convention:

RiemannToChristoffel[RiemannCD[-d,-c,-b,a]]

xCoba and Regge Wheeler Coordinates

define Regge Wheeler coordinates/chart in xCoba: RW
indices run through {0,1,2,3} and correspond to coordinate scalars {t[],r[],θ[],ϕ[]}.

DefChart[RW,H,{0,1,2,3},{t[],r[],θ[],ϕ[]}];

define metric components of spherically symmetric background:

DefScalarFunction/@{f,bb};
metricBG=DiagonalMatrix[{-f[r[]],bb[r[]],r[]^2,r[]^2Sin[θ[]]^2}];
metricBG//MatrixForm

a rule inserting the functions f and bb explicitly:

DefConstantSymbol[M];
backg={bb->Function[r,1/(1-2M/r)],f->Function[r,1-2M/r]};

define ϵ for counting the perturbation order and a function “linearize” which Taylor-expands to first order:

DefConstantSymbol[ϵ];
linearize[expr_]:=Collect[Series[expr,{ϵ,0,1}]//Normal,ϵ,Together];

define frequency ω, angular momentum quantum numbers, Legendre polynomial PL, and spherical harmonics YLM. we use the shorthand notation ll:=l(l+1).

DefConstantSymbol/@{ω,ll,m};
DefScalarFunction[PL];
PL''[θ_]:=-(ll-m^2 Csc[θ]^2) PL[θ]-Cot[θ]PL'[θ];
YLM:=E^(-I m ϕ[]-I ω t[])PL[θ[]];

even metric perturbation before gauge fixing, see (13) in [1].

DefScalarFunction/@{H0,H1,H2,KK,h0,h1,GG};metricPertEvenGeneric=ϵ {{f[r[]]H0[r[]]YLM,H1[r[]]YLM,h0[r[]] D[YLM,θ[]],h0[r[]] D[YLM,ϕ[]]},{H1[r[]]YLM,bb[r[]]H2[r[]]YLM,h1[r[]] D[YLM,θ[]],h1[r[]] D[YLM,ϕ[]]},{h0[r[]] D[YLM,θ[]],h1[r[]] D[YLM,θ[]],r[]^2(KK[r[]]YLM+GG[r[]]D[YLM,{θ[],2}]),r[]^2 GG[r[]](D[D[YLM,ϕ[]],θ[]]-((Cos[θ[]])/(Sin[θ[]]))D[YLM,ϕ[]])},{h0[r[]] D[YLM,ϕ[]],h1[r[]] D[YLM,ϕ[]],r[]^2 GG[r[]](D[D[YLM,ϕ[]],θ[]]-((Cos[θ[]])/(Sin[θ[]]))D[YLM,ϕ[]]),r[]^2(KK[r[]]Sin[θ[]]^2 YLM+GG[r[]](D[YLM,{ϕ[],2}]+Sin[θ[]]Cos[θ[]]D[YLM,θ[]]))}};
metricPertEvenGeneric//MatrixForm

odd metric perturbation before gauge fixing, see (12) in [1].

DefScalarFunction/@{hh0,hh1,hh2};metricPertOddGeneric=ϵ {{0,0,((-hh0[r[]])/(Sin[θ[]])) D[YLM,ϕ[]],hh0[r[]]Sin[θ[]] D[YLM,θ[]]},{0,0,((-hh1[r[]])/(Sin[θ[]])) D[YLM,ϕ[]],hh1[r[]]Sin[θ[]] D[YLM,θ[]]},{((-hh0[r[]])/(Sin[θ[]])) D[YLM,ϕ[]],((-hh1[r[]])/(Sin[θ[]])) D[YLM,ϕ[]],hh2[r[]]((1/(Sin[θ[]]))D[D[YLM,ϕ[]],θ[]]-((Cos[θ[]])/((Sin[θ[]])^2))D[YLM,ϕ[]]),(1/2)hh2[r[]]((1/(Sin[θ[]]))D[YLM,{ϕ[],2}]+Cos[θ[]]D[YLM,θ[]]-Sin[θ[]]D[YLM,{θ[],2}])},{hh0[r[]]Sin[θ[]] D[YLM,θ[]],hh1[r[]]Sin[θ[]] D[YLM,θ[]],(1/2)hh2[r[]]((1/(Sin[θ[]]))D[YLM,{ϕ[],2}]+Cos[θ[]]D[YLM,θ[]]-Sin[θ[]]D[YLM,{θ[],2}]),-hh2[r[]](Sin[θ[]]D[D[YLM,ϕ[]],θ[]]-Cos[θ[]]D[YLM,ϕ[]])}};
metricPertOddGeneric//MatrixForm

tell xCoba that the metric in RW coordinates is given by the matrix metricBG.
the minus in front of RW means we define lower-index components.

?MetricCompute

MetricInBasis[g,-RW,metricBG];

calculate all curvature tensors etc:

MetricCompute[g,RW,"Christoffel"[1, -1, -1],CVSimplify->linearize];

Gauge fixing the even sector

now we will gauge fix so that the components (Subscript[h, 02]), (Subscript[h, 03]), (Subscript[h, 12]), (Subscript[h, 13]), (Subscript[h, 23]) of the even metric perturbation all vanish (Regge-Wheeler gauge).
under an infinitesimal coordinate transformation (x^a)->(x^a)+(ξ^a) the metric perturbation changes as (Subscript[h, ab])-> (Subscript[h, ab])+δ (Subscript[h, ab]), δ (Subscript[h, ab]):=(Subscript[∇, a])(Subscript[ξ, b])+(Subscript[∇, b])(Subscript[ξ, a]).

DefTensor[ξ[-a],H];
(CD[-a][ξ[-b]]+CD[-b][ξ[-a]])

choose the form of (ξ^a) appropriate to simplify the even metric perturbation, see (19) in [1].

DefScalarFunction/@{M0,M1,M2};
ComponentValue[ComponentArray[ξ[{a,RW}]],{M0[r[]]YLM,M1[r[]]YLM,M2[r[]]D[YLM,θ[]],((M2[r[]])/(Sin[θ[]]^2))D[YLM,ϕ[]]}]

ChangeComponents[ξ[-{a,RW}],ξ[{a,RW}]]

δh=(CD[-a][ξ[-b]]+CD[-b][ξ[-a]])//ToBasis[RW]//ToBasis[RW]//CovDToChristoffel//TraceBasisDummy//ComponentArray//ToValues//ToValues;
δh//MatrixForm

we get three independent equations for M0, M1 and M2 from the requirement that the components (Subscript[h, 02]), (Subscript[h, 03]), (Subscript[h, 12]), (Subscript[h, 13]), (Subscript[h, 23]) of the new metric perturbation vanish:

metricPertEvenNew=metricPertEvenGeneric+ϵ δh;Simplify[((E^(+I ω t+I m ϕ))/(ϵ PL'[θ]))metricPertEvenNew[[1,3]]]==0
Simplify[I((E^(+I ω t+I m ϕ))/(ϵ m PL[θ]))metricPertEvenNew[[1,4]]]==0
Simplify[((E^(+I ω t+I m ϕ))/(ϵ PL'[θ]))metricPertEvenNew[[2,3]]]==0
Simplify[I((E^(+I ω t+I m ϕ))/(ϵ m PL[θ]))metricPertEvenNew[[2,4]]]==0
Simplify[-I((E^(+I ω t+I m ϕ))/(ϵ m ((r[])^2)(Cot[θ] PL[θ]-PL'[θ])))metricPertEvenNew[[3,4]]]==0

determine the function M2:

Solve[metricPertEvenNew[[3,4]]==0,M2[r[]]][[1,1]]

a rule inserting the function M2 explicitly:

replaceM2=M2->Function[r,-((GG[r])/2)];

determine functions M0 and M1:

Solve[(metricPertEvenNew[[1,3]]/.replaceM2)==0,M0[r[]]][[1,1]]
Solve[(metricPertEvenNew[[2,3]]/.replaceM2)==0,M1[r[]]][[1,1]]

rule inserting functions M0, M1 and M2 explicitly:

replaceM={M0->Function[r,((2 h0[r]+I ω GG[r] (r^2))/(2 f[r]))],M1->Function[r,((-2 h1[r]+(r^2) GG'[r])/(2 bb[r]))],replaceM2};

insert the functions M0, M1 and M2 just found in the new metric perturbation to obtain the gauge-fixed form for the even metric perturbation.

Simplify[metricPertEvenNew/.replaceM]//MatrixForm

redefine functions H0, H1, H2 and KK to simplify the presentation.

DefScalarFunction/@{H0new,H1new,H2new,KKnew};Simplify[Solve[-(( (2 bb[r] (-2 I ω h0[r]-f[r] H0+(ω^2) GG[r] (r^2))+f'[r] (-2 h1[r]+(r^2) GG'[r])))/(2 bb[r]))==f[r]H0new[r],H0]][[1,1]]
Simplify[Solve[((((2 h0[r]+I ω GG[r] (r^2)) f'[r]+2 f[r] (I ω h1[r]+H1-I ω GG[r] r-I ω (r^2) GG'[r]-h0'[r])))/(2 f[r]))==H1new[r],H1]][[1,1]]
Simplify[Solve[(1/(2 bb[r])) (2 ((bb[r])^2) H2+bb'[r] (2 h1[r]-(r^2) GG'[r])+2 bb[r] (2 r GG'[r]-2 h1'[r]+(r^2) GG''[r]))==bb[r]H2new[r],H2]][[1,1]]
Simplify[Solve[(( r (-2 h1[r]+r (bb[r] KK+r GG'[r])))/(bb[r]))==(r^2)KKnew[r],KK]][[1,1]]

replaceHK={H0[r[]]->(((-2 I ω h0[r[]]+f[r[]] H0new[r[]]+(ω^2) GG[r[]] ((r[])^2)))/(f[r[]]))+((f'[r[]](-2 h1[r[]]+((r[])^2) GG'[r[]]))/(2 bb[r[]]f[r[]])), H1[r[]]->-(2 h0[r[]]+I ω GG[r[]] ((r[])^2)) ((f'[r[]])/(2 f[r[]]))+ (-I ω h1[r[]]+H1new[r[]]+I ω GG[r[]]r[]+I ω ((r[])^2) GG'[r[]]+h0'[r[]]),H2[r[]]->H2new[r[]]+((bb'[r[]] (-2 h1[r[]]+((r[])^2) GG'[r[]])-2 bb[r[]] (2 r[] GG'[r[]]-2 h1'[r[]]+((r[])^2) GG''[r[]]))/(2 ((bb[r[]])^2))),KK[r[]]->((((r[])^2)bb[r[]] KKnew[r[]]+2 h1[r[]] r[]-((r[])^3) GG'[r[]])/(bb[r[]] ((r[])^2)))};

Simplify[metricPertEvenNew/.replaceM/.replaceHK]//MatrixForm

finally we obtain the even metric perturbation in Regge-Wheeler gauge, see (20) in [1]:

%/.H0new->H0/.H1new->H1/.H2new->H2/.KKnew->KK//MatrixForm

Gauge fixing the odd sector

now we will gauge fix so that the components (Subscript[h, 22]), (Subscript[h, 23]), (Subscript[h, 33]) of the odd metric perturbation all vanish (Regge-Wheeler gauge).
under an infinitesimal coordinate transformation (x^a)->(x^a)+(η^a) the metric perturbation changes as (Subscript[h, ab])-> (Subscript[h, ab])+δ (Subscript[h, ab]), δ (Subscript[h, ab]):=(Subscript[∇, a])(Subscript[η, b])+(Subscript[∇, b])(Subscript[η, a]).

DefTensor[η[-a],H];
(CD[-a][η[-b]]+CD[-b][η[-a]])

choose the form of (η^a) appropriate to simplify the odd metric perturbation, see (17) in [1].

DefScalarFunction/@{Λ};
ComponentValue[ComponentArray[η[{a,RW}]],{0,0,-((Λ[r[]])/(((r[])^2)Sin[θ[]]))D[YLM,ϕ[]],((Λ[r[]])/(((r[])^2)Sin[θ[]]))D[YLM,θ[]]}]

ChangeComponents[η[-{a,RW}],η[{a,RW}]]

δhodd=(CD[-a][η[-b]]+CD[-b][η[-a]])//ToBasis[RW]//ToBasis[RW]//CovDToChristoffel//TraceBasisDummy//ComponentArray//ToValues//ToValues//ToValues;
δhodd//MatrixForm

we get a single equation for Λ from the requirement that the components (Subscript[h, 22]), (Subscript[h, 23]), (Subscript[h, 33]) of the new metric perturbation vanish:

metricPertOddNew=metricPertOddGeneric+ϵ δhodd;
Simplify[((E^(I (ω t+m ϕ)))/(I m ϵ))metricPertOddNew[[3,3]]]==0
Simplify[-2((E^(I (ω t+m ϕ)))/ϵ)metricPertOddNew[[3,4]]]==0
Simplify[-((E^(I (ω t+m ϕ)))/(I m ϵ))metricPertOddNew[[4,4]]]==0

determine the function Λ:

Solve[metricPertOddNew[[3,3]]==0,Λ[r[]]][[1,1]]

a rule inserting the function Λ explicitly:

replaceΛ=Λ->Function[r,((hh2[r])/2)];

insert the function Λ just found in the new metric perturbation to obtain the gauge-fixed form for the odd metric perturbation.

Simplify[metricPertOddNew/.replaceΛ]//MatrixForm

redefine functions hh0 and hh1 to simplify the presentation.

DefScalarFunction/@{hh0new,hh1new};Simplify[Solve[(hh0+(ω/(2I)) hh2[r])==  hh0new[r],hh0]][[1,1]]
Simplify[Solve[(-(1/r) hh2[r]+(hh1+(1/2)hh2'[r]))==hh1new[r],hh1]][[1,1]]

replacehh={hh0[r[]]->hh0new[r[]]+(1/2) I ω hh2[r[]], hh1[r[]]->hh1new[r[]]+((hh2[r[]])/(r[]))-((hh2'[r[]])/2)};

Simplify[metricPertOddNew/.replaceΛ/.replacehh]//MatrixForm

finally we obtain the odd metric perturbation in Regge-Wheeler gauge, see (18) in [1]:

%/.hh0new->h0/.hh1new->h1//MatrixForm

Combining the background with the even and odd sectors of the metric perturbation

Clear[H0,H1,H2,KK,h0,h1]

even metric perturbation in Regge-Wheeler gauge, see (20) in [1].
notice that this results from a gauge fixing from the generic ansatz in (13).

DefScalarFunction/@{H0,H1,H2,KK};metricPertEven=ϵ YLM{{f[r[]]H0[r[]],H1[r[]],0,0},{H1[r[]],bb[r[]]H2[r[]],0,0},{0,0,KK[r[]]r[]^2,0},{0,0,0,KK[r[]]r[]^2Sin[θ[]]^2}};
metricPertEven//MatrixForm

odd metric perturbation in Regge Wheeler gauge, see (18) in [1].
notice that we do not set m=0, so we have more components than in [1] here.
further this results from a gauge fixing from the generic ansatz in (12).

DefScalarFunction/@{h0,h1};
DefConstantSymbol[odd];metricPertOdd=ϵ odd{{0,0,-D[YLM,ϕ[]]/Sin[θ[]]h0[r[]],Sin[θ[]]D[YLM,θ[]]h0[r[]]},{0,0,-D[YLM,ϕ[]]/Sin[θ[]]h1[r[]],Sin[θ[]]D[YLM,θ[]]h1[r[]]},{-D[YLM,ϕ[]]/Sin[θ[]]h0[r[]],-D[YLM,ϕ[]]/Sin[θ[]]h1[r[]],0,0},{Sin[θ[]]D[YLM,θ[]]h0[r[]],Sin[θ[]]D[YLM,θ[]]h1[r[]],0,0}};
metricPertOdd//MatrixForm

put it all together:

metricFull=metricBG+metricPertEven+metricPertOdd;

tell xCoba that the metric in RW coordinates is given by the matrix metricFull.
the minus in front of RW means we define lower-index components.

MetricInBasis[g,-RW,metricFull];

calculate all curvature tensors etc:

MetricCompute[g,RW,All,CVSimplify->linearize];

components in RW coordintes are indicated by colored indices:

RicciCD[{a,RW},{b,-RW}]
RicciCD[{0,RW},{1,-RW}]

not all index positions are calculated by default.
here is an example which calculates Ricci with mixed index positions:

ChangeComponents[RicciCD[{a,RW},-{b,RW}],RicciCD[-{a,RW},-{b,RW}]];

show result for mixed Ricci by calling ToValues:

RicciCD[{0,RW},{1,-RW}]
RicciCD[{0,RW},{1,-RW}]//ToValues

we need to apply ToValues repeatedly until all components are inserted:

RicciCD[{0,RW},{1,-RW}]//ToValues//ToValues

Field Equations

Einstein field equations. we use a mixed index version (which is better for extending the calculation to neutron stars).

EFE=RicciCD[c,-b]-1/2delta[c,-b]RicciScalarCD[];
EFE==0

change all indices to Regge Wheeler (RW) basis.

EFEbasis=EFE//ToBasis[RW]

convert to matrix:

EFEmatrix=EFEbasis//ComponentArray

footnote: if necessary the sum convention can be performed using TraceBasisDummy:

DefTensor[u[a],H];
u[a]u[-a]//ToBasis[RW]//ComponentArray
u[a]u[-a]//ToBasis[RW]//ComponentArray//TraceBasisDummy

insert component values:

EFEvalues=EFEmatrix//ToValues//ToValues

remove brackets from coordinates and linearize:

EFERW=EFEvalues/.{t[]->t,r[]->r,θ[]->θ,ϕ[]->ϕ}//linearize

Odd Part: Regge-Wheeler Equation

take the odd part of the field equations:

EFERWodd=Coefficient[EFERW,odd,1]//Together

solve the θ-ϕ-component for the metric function h0[r]. Notice that lists/arrays in Mathematica always start with 1, so the coordinate numbers are shifted up by +1 in EFERWodd.

h0sol=h0[r]/.Solve[EFERWodd[[3,4]]==0,h0[r]]//FullSimplify//First

insert h0[r] back into the field equations:

EFERWoddWithh0=EFERWodd/.h0-> Function[r,h0sol//Evaluate]

solve the r-θ-component for h1’’[r]:

h1ppsol=h1''[r]/.Solve[0==EFERWoddWithh0[[2,4]],h1''[r]]//Simplify//First

define Regge-Wheeler Q[r] function, eq. (23) in [1]:

Qdef=h1-> Function[r,r/Sqrt[f[r]/bb[r]] Q[r]];

derive the Regge-Wheeler equation in terms of the ordinary radial coordinate r:

RWeq=h1''[r]-h1ppsol/.Qdef//Simplify

define tortoise coordinate rT, see between (24) and (25) in [1]:

DefScalarFunction[rT,PrintAs->"(SubStar[r])"];
rTdef={rT'[r]-> 1/Sqrt[f[r]/bb[r]],rT''[r]-> D[1/Sqrt[f[r]/bb[r]],r]};

notice the effect of the PrintAs option:

rT

change radial coordinate to the tortoise one:

RWeqTort=RWeq/.Q->Function[r,Q[rT[r]]]//Simplify

replace derivatives of the tortoise coordinate w.r.t. r-coordinate:

RWeqTort2=RWeqTort/.rTdef//Simplify

normalize the coefficient in front of Q''[(SubStar[r])[r]] to 1 and insert Schwarzschild metric for the background functions f[r] and bb[r]:

RWeqSchwarz=RWeqTort2 Sqrt[f[r]/bb[r]]^3/r/.backg

simplify the result:

RWE=Collect[RWeqSchwarz,{Q[_],Q''[_]},FullSimplify]

Even Part: Zerilli Equation

take the even part of the field equations:

EFERWeven=Coefficient[EFERW/.odd-> 0,ϵ,1]//Simplify

solve for H2[r]:

H2sol=H2[r]/.Solve[EFERWeven[[3,4]]==0,H2[r]]//First

derive the starting equations in the Zerilli paper, see (F3) in [2]:

zerillieq0=Solve[{EFERWeven[[3,1]]==0,EFERWeven[[1,2]]==0,EFERWeven[[2,3]]==0},{H1'[r],KK'[r],H0'[r]}]//Simplify//First

solve the algebraic equation (r-r-component) for H0[r]:

H0sol=H0[r]/.Solve[0==EFERWeven[[2,2]]/.zerillieq0,H0[r]]//Simplify//First

absorb a factor of ω in R[r] compared to H1[r]:

Rdef=H1->Function[r,ω R[r]];

obtain eq. (1) in [2]:

Zeq1rhs={KK'[r],H1'[r]/ω}/.zerillieq0/.H2[r]-> H2sol/.H0[r]-> H0sol/.Rdef;
Zeq1={KK'[r],R'[r]}-Zeq1rhs

replace R and K by “hatted” functions Rh and Kh:

Zeq1replaced=Zeq1/.{KK-> Function[r,f1[r] Kh[r]+ f2[r] Rh[r]],R-> Function[r,f3[r] Kh[r]+f4[r] Rh[r]]}

change to tortoise coordinate:

Zeq1tor=Zeq1replaced/.{Kh-> Function[r,Kh[rT[r]]],Rh-> Function[r,Rh[rT[r]]]}/.rTdef

solve for “hatted” functions Kh and Rh:

Zeqsol=Solve[Zeq1tor=={0,0},{Kh'[rT[r]],Rh'[rT[r]]}]//Simplify//First

assume f2[r] is 1 for simplicity. we need to find just one solution that works ;)

Zeqsol=Zeqsol/.f2->Function[r,1];

formulate requirements on the transformtion:

requirement1=Kh'[rT[r]]-Rh[rT[r]]/.Zeqsol;
requirement2=Coefficient[Rh'[rT[r]]/.Zeqsol,Rh[rT[r]]];

the requirements lead to algebraic and differential equations for f1, f3, and f4. it is enough to pick the algebraic equations and check the differential ones afterwards.
we start by solving for f4[r]:

f4sol=f4[r]/.Solve[0==Coefficient[Coefficient[requirement1,Rh[rT[r]]],ω^2],f4[r]]//First

solve for f3[r]:

f3sol=f3[r]/.Solve[0==Coefficient[requirement1/.f4->Function[r,f4sol//Evaluate],Rh[rT[r]]],f3[r]]//First//Simplify

solve for f1[r]:

f1sol=f1[r]/.Solve[0==requirement2/.{f4-> Function[r,f4sol//Evaluate],f3-> Function[r,f3sol//Evaluate]},f1[r]]//First//Simplify

plug in solutions for f1, f3, and f4:

Zeqsol1=Zeqsol/.f4-> Function[r,f4sol//Evaluate]/.f3-> Function[r,f3sol//Evaluate]/.f1-> Function[r,f1sol//Evaluate]//Simplify

simplify the result. use λ instead of ll to allow direct comparison with [2].

Zeqs=Collect[Zeqsol1/.backg/.ll->2(λ+1),{Kh[_],Rh[_],ω},FullSimplify]

transform the two first-order equations to one second order equation :

ZE=Kh''[r]-(Rh'[rT[r]]/.Zeqs)

Results

result for the Regge-Wheeler equation. notice that it has the form of a Schoedinger equation, but the radial coordinate in the argument of Q and in the potential are not the same. the same holds for the Zerilli equation.

RWE==0

isolate the potential:

V=-Coefficient[RWE,Q[rT[r]]]+ω^2

obtain rT in terms of r:

rT==Integrate[rT'[r]/.rTdef/.backg,r]//PowerExpand

plot the potential for M=1, l=2, or ll:=l(l+1)=6, over the tortoise coordinate. the horizon is located at negative infinity in tortoise coordinate.

table=Table[({r+2 M Log[-2 M+r],V}/.r->(2M+10^i))/.{M->1,ll->6},{i,-4,2,1/100}];
ListLinePlot[table,AxesOrigin->{-20,0}]

result for the Zerilli equation. notice that it is possible to transform it to the Regge-Wheeler equation (borrowing results from “supersymmetric quantum mechanics”). so the odd and even spectra are actually the same! in practice one would just solve the simpler Regge-Wheeler equation, e.g., numerically.

ZE==0