Massive spin-2 fields on black hole spacetimes 
Reference: “Massive spin-2 fields on black hole spacetimes: instability of the Schwarzschild and Kerr solutions and bounds on graviton mass,” R. Brito, V.Cardoso, P .Pani (2013).

SetDirectory[NotebookDirectory[]];

We use a slightly different notation than in the paper:

1) Axial functions:
h0 -> hh0; h1 -> hh1; h2 -> hh2

2) Polar functions:
η0 -> h0; η1 -> h1;

3) ã -> a with ã=((ϵ J)/(M^2)) and ϵ << 1;

Perturbation equations

Coefficients at first order in the rotation

(*POLAR FUNCTIONS*)
A0=(1/(r^4))((-2 (M^2)+2 M r (Λ+(r^2) (μ^2))-(r^2) (Λ+(r^2) ((μ^2)-(ω^2)))) H0[r]-4 I M (r^2) ω H1[r]+6 (M^2) H2[r]-4 M r H2[r]-8 (M^2) KK[r]+4 M r KK[r]+4 (M^2) r H0'[r]-6 M (r^2) H0'[r]+2 (r^3) H0'[r]+4 (M^2) (r^2) H0''[r]-4 M (r^3) H0''[r]+(r^4) H0''[r]);
A1=(1/((2 M-r) (r^3)))((4 M Λ-2 r Λ) h0[r]+2 I M (r^2) ω H0[r]+4 (M^2) H1[r]-4 M r H1[r]+2 (r^2) H1[r]-2 M r Λ H1[r]+(r^2) Λ H1[r]-2 M (r^3) (μ^2) H1[r]+(r^4) (μ^2) H1[r]-(r^4) (ω^2) H1[r]+2 I M (r^2) ω H2[r]-4 (M^2) r H1'[r]+6 M (r^2) H1'[r]-2 (r^3) H1'[r]-4 (M^2) (r^2) H1''[r]+4 M (r^3) H1''[r]-(r^4) H1''[r]);
A2=(1/((r^3) (((-2 M+r))^2)))(2 M (3 M-2 r) r H0[r]+4 (((-2 M+r))^2) Λ h1[r]+r (-4 I M (r^2) ω H1[r]+(-18 (M^2)+2 M r (8+Λ+(r^2) (μ^2))-(r^2) (4+Λ+(r^2) ((μ^2)-(ω^2)))) H2[r]+(2 M-r) (4 (3 M-r) KK[r]+r (2 (M-r) H2'[r]+(2 M-r) r H2''[r]))));
A3=-(1/((2 M-r) (r^2)))2 (2 M r (-2 M+r) H0[r]-2 (((-2 M+r))^2) Λ h1[r]+r (2 (6 (M^2)-5 M r+(r^2)) H2[r]+(-16 (M^2)+2 M r (6+Λ+(r^2) (μ^2))-(r^2) (2+Λ+(r^2) ((μ^2)-(ω^2)))) KK[r]+(2 M-r) r (2 (M-r) KK'[r]+(2 M-r) r KK''[r])));
A4=((-Λ h0[r]+I (r^2) ω H0[r]-2 M H1[r]+2 r H1[r]-2 M r H1'[r]+(r^2) H1'[r])/(r^2));
A5=((M r H0[r]+(2 M-r) Λ h1[r]+r (I (r^2) ω H1[r]+(-3 M+2 r) H2[r]+(2 M-r) (2 KK[r]-r H2'[r])))/((r^2) (-2 M+r)));
AL=-H0[r]+H2[r]+2 KK[r];
C0=((4 a (M^2) (M (-4 M+r) h0[r]+(2 M-r) (r^2) (I r ω H0[r]+H1[r]-h0'[r])))/((2 M-r) (r^6)));
C1=((2 a (M^2) (I (4 M-3 r) (r^2) ω h0[r]+(2 M-r) ((4 M-3 r) r H0[r]+(-6 M+4 r) h1[r]+r (-2 I (r^2) ω H1[r]-r H2[r]+(-2 M+r) h1'[r]))))/((r^5) (((-2 M+r))^2)));
C2=((4 a (M^2) ((12 (M^2)-17 M r+6 (r^2)) h0[r]+I (2 M-r) r ((4 M-3 r) ω h1[r]+I (4 M-3 r) H1[r]+(r^2) ω H2[r])))/((((2 M-r))^3) (r^4)));
C3=-((4 a (M^2) (4 M h0[r]+r (-3 I (2 M-r) ω h1[r]+2 I (r^2) ω KK[r]+3 (-2 M+r) h0'[r])))/((2 M-r) (r^3)));
C4=((2 a (M^2) (I r ω h0[r]+(2 M-r) H0[r]))/((r^3) (-2 M+r)));
C5=-((2 a (M^2) ((4 M-3 r) h0[r]+I (2 M-r) r (ω h1[r]+I H1[r])))/((r^3) (((-2 M+r))^2)));
CL=((4 a (M^2) h0[r])/((2 M-r) (r^2)));
α0=(((8 (M^2)-2 M r (2+Λ+(r^2) (μ^2))+(r^2) (Λ+(r^2) ((μ^2)-(ω^2)))) h0[r]-(2 M-r) r (2 I M ω h1[r]+(2 M-r) (2 H1[r]+r h0''[r])))/((2 M-r) (r^3)));
α1=(1/((r^3) (((-2 M+r))^2)))(2 (r^2) (((-2 M+r))^2) (-2+Λ) GG[r]-2 I M (r^3) ω h0[r]+(2 M-r) ((16 (M^2)-2 M r (8+Λ+(r^2) (μ^2))+(r^2) (4+Λ+(r^2) ((μ^2)-(ω^2)))) h1[r]-(2 M-r) r (-2 r H2[r]+2 r KK[r]-4 M h1'[r]+2 M r h1''[r]-(r^2) h1''[r])));
α2=(((2 M-r) (r^2) (-2+Λ) GG[r]+I (r^3) ω h0[r]+(2 M-r) (2 (M-r) h1[r]-r (r KK[r]+(-2 M+r) h1'[r])))/((r^2) (-2 M+r)));
βt0=((2 a (M^2) (2 (2 M-r) (r^2) (-2+Λ) GG[r]+2 I (r^3) ω h0[r]+(2 M-r) (20 M h1[r]-11 r h1[r]+2 (r^2) KK[r])))/((2 M-r) (r^5)));
βt1=-((2 a (M^2) ((-16 M+9 r) h0[r]-2 I (2 M-r) r ω h1[r]))/((r^3) (((-2 M+r))^2)));
βt2=((2 a (M^2) h0[r])/((r^2) (-2 M+r)));
ζ0=-((a (M^2) (-2 I (r^2) Λ ω h0[r]+(2 M-r) (-2 M H0[r]+Λ h1[r]-6 I (r^2) ω H1[r]+6 M H2[r]-8 r H2[r]-8 M KK[r]+8 r KK[r]+12 M r H0'[r]-6 (r^2) H0'[r]+2 (r^2) KK'[r])))/((2 M-r) (r^4)));
ζ1=((a (M^2) ((4 M-3 r) Λ h0[r]+2 (I (2 M-r) r Λ ω h1[r]+6 (2 (M^2)-3 M r+(r^2)) H1[r]+r (-3 I (2 M-r) r ω H2[r]+I (4 M-3 r) r ω KK[r]+3 (((-2 M+r))^2) H1'[r]))))/((r^3) (((-2 M+r))^2)));
ζ2=-((a (M^2) (Λ h0[r]+2 I (r^2) ω KK[r]))/((2 M-r) (r^2)));
ξ0=-((a (M^2) (4 r (-2 M+r) GG[r]+2 I (r^2) ω h0[r]+(2 M-r) (-h1[r]+2 (r^2) GG'[r])))/((2 M-r) (r^4)));
ξ1=((a (M^2) (2 I (4 M-3 r) (r^2) ω GG[r]+(-4 M+3 r) h0[r]-2 I (2 M-r) r ω h1[r]))/((r^3) (((-2 M+r))^2)));
ξ2=((a (M^2) (-2 I (r^2) ω GG[r]+h0[r]))/((2 M-r) (r^2)));
Ks=-(((r^2) ((2 M-r) (-2+Λ+(r^2) (μ^2))+(r^3) (ω^2)) GG[r]+(2 M-r) ((4 M-2 r) h1[r]+(r^2) (2 (M-r) GG'[r]+(2 M-r) r GG''[r])))/((2 M-r) (r^2)));
Ksh=-((4 I a (M^2) ω GG[r])/(2 M-r));
Kf=((6 a (M^2) (-2 h0[r]+r (I ω h1[r]+h0'[r])))/(r^3));

(*AXIAL FUNCTIONS*)
At0=((8 a (M^2) Λ hh0[r])/(r^5));
At1=((4 a (M^2) Λ hh1[r])/(r^5));
At2=0;
At3=0;
At4=0;
At5=0;
AtL=0;
B0=((4 a (M^2) (M (4 M-r) hh0[r]+(2 M-r) (r^2) hh0'[r]))/((2 M-r) (r^6)));
B1=((2 a (M^2) (-I (4 M-3 r) (r^2) ω hh0[r]+(2 M-r) ((6 M-4 r) hh1[r]+(2 M-r) r hh1'[r])))/((r^5) (((-2 M+r))^2)));
B2=-((4 a (M^2) (4 M-3 r) ((3 M-2 r) hh0[r]+I (2 M-r) r ω hh1[r]))/((((2 M-r))^3) (r^4)));
B3=((4 a (M^2) (4 M hh0[r]-3 I (2 M-r) r (ω hh1[r]-I hh0'[r])))/((2 M-r) (r^3)));
B4=((2 I a (M^2) ω hh0[r])/((2 M-r) (r^2)));
B5=((2 a (M^2) ((4 M-3 r) hh0[r]+I (2 M-r) r ω hh1[r]))/((r^3) (((-2 M+r))^2)));
BL=-((4 a (M^2) hh0[r])/((2 M-r) (r^2)));
αt0=((2 a (M^2) (2 I (r^3) ω hh0[r]+(2 M-r) ((20 M-11 r) hh1[r]-2 (-2+Λ) hh2[r])))/((2 M-r) (r^5)));
αt1=((2 a (M^2) ((16 M-9 r) hh0[r]+2 I (2 M-r) r ω hh1[r]))/((r^3) (((-2 M+r))^2)));
αt2=((2 a (M^2) hh0[r])/((r^2) (-2 M+r)));
β0=-(((8 (M^2)-2 M r (2+Λ+(r^2) (μ^2))+(r^2) (Λ+(r^2) ((μ^2)-(ω^2)))) hh0[r]-(2 M-r) r (2 I M ω hh1[r]+(2 M-r) r hh0''[r]))/((2 M-r) (r^3)));
β1=-((-2 I M (r^3) ω hh0[r]+(2 M-r) ((16 (M^2)-2 M r (8+Λ+(r^2) (μ^2))+(r^2) (4+Λ+(r^2) ((μ^2)-(ω^2)))) hh1[r]-(2 M-r) (2 (-2+Λ) hh2[r]-r (4 M hh1'[r]+r (-2 M+r) hh1''[r]))))/((r^3) (((-2 M+r))^2)));
β2=-((I (r^3) ω hh0[r]+(2 M-r) (2 (M-r) hh1[r]-(-2+Λ) hh2[r]+(2 M-r) r hh1'[r]))/((r^2) (-2 M+r)));
η0=((a (M^2) Λ (2 I (r^2) ω hh0[r]+(-2 M+r) hh1[r]))/((2 M-r) (r^4)));
η1=((a (M^2) Λ ((4 M-3 r) hh0[r]+2 I (2 M-r) r ω hh1[r]))/((r^3) (((-2 M+r))^2)));
η2=((a (M^2) Λ hh0[r])/((r^2) (-2 M+r)));
χ0=((a (M^2) (2 I (r^3) ω hh0[r]-(2 M-r) (r hh1[r]-8 hh2[r]+2 r hh2'[r])))/((2 M-r) (r^5)));
χ1=((a (M^2) ((4 M-3 r) hh0[r]+2 I ω ((2 M-r) r hh1[r]+(4 M-3 r) hh2[r])))/((r^3) (((-2 M+r))^2)));
χ2=((a (M^2) (hh0[r]+2 I ω hh2[r]))/((r^2) (-2 M+r)));
Kt=((2 r (((-2 M+r))^2) hh1[r]+(-(2 M-r) (8 M+r (-4+Λ+(r^2) (μ^2)))-(r^4) (ω^2)) hh2[r]+r (2 (6 (M^2)-5 M r+(r^2)) hh2'[r]-r (((-2 M+r))^2) hh2''[r]))/((2 M-r) (r^3)));
Kth=-((4 I a (M^2) ω hh2[r])/((2 M-r) (r^2)));
Kg=((6 a (M^2) (-2 hh0[r]+r (I ω hh1[r]+hh0'[r])))/(r^3));

(*DumpSave["coefficients_Richard/coefficients.mx",{A0,A1,A2,A3,A4,A5,AL,B0,B1,B2,B3,B4,B5,BL,At0,At1,At2,At3,At4,At5,AtL,C0,C1,C2,C3,C4,C5,CL,α0,α1,α2,αt0,αt1,αt2,β0,β1,β2,βt0,βt1,βt2,η0,η1,η2,ξ0,ξ1,ξ2,χ0,χ1,χ2,ζ0,ζ1,ζ2,Ks,Ksh,Kt,Kth,Kf,Kg}];*)

Equations without the couplings l±1

The coefficients are written in an external file but we state them in the previous subsection. For ã=0 we recover the Schwarzschild perturbations.

(*<<"coefficients_Richard/coefficients.mx"*)

ruleΛ=Λ->l(l+1);

(*FIRST GROUP*)
(*Polar*)
EQ1=FullSimplify[(A0+I m C0)/.ruleΛ];
EQ2=FullSimplify[(A1+I m C1)/.ruleΛ];
EQ3=FullSimplify[(A2+I m C2)/.ruleΛ];
EQ4=FullSimplify[(A3+I m C3)/.ruleΛ];

(*SECOND GROUP*)
(*Polar*)
EQ5=FullSimplify[(l(l+1)α0+I m((l-1)(l+2)ξ0-βt0-ζ0))/.ruleΛ];
EQ6=FullSimplify[(l(l+1)α1+I m((l-1)(l+2)ξ1-βt1-ζ1))/.ruleΛ];
(*Axial*)
EQ7=FullSimplify[(l(l+1)β0+I m((l-1)(l+2)χ0+αt0+η0))/.ruleΛ];
EQ8=FullSimplify[(l(l+1)β1+I m((l-1)(l+2)χ1+αt1+η1))/.ruleΛ];

(*THIRD GROUP*)
(*Polar*)
EQ9=FullSimplify[(l(l+1)(Ks+I m Ksh)-I m Kf)/.ruleΛ];
(*Axial*)
EQ10=FullSimplify[(l(l+1)(Kt+I m Kth)+I m Kg)/.ruleΛ];

(*CONSTRAINTS*)
(*Trace: Polar*)
EQT=FullSimplify[(AL+I m CL)/.ruleΛ];
(*Tranverse constraint: Polar*)
EQD1=FullSimplify[(A4+I m C4)/.ruleΛ];
EQD2=FullSimplify[(A5+I m C5)/.ruleΛ];
EQD3=FullSimplify[(l(l+1)α2+I m((l-1)(l+2)ξ2-βt2-ζ2))/.ruleΛ];
(*Tranverse constraint: Axial*)
EQD4=FullSimplify[(l(l+1)β2+I m((l-1)(l+2)χ2+αt2+η2))/.ruleΛ];

Axial Sector

Simplify the equations

Eliminating h0 and defining functions Q and Z

Clear[hh0,hh1,hh2]

f[r_]:=1-((2M)/r);
hh1[r_]:=((Q[r])/( f[r]))(1-((a(l(l+1)+2)m (M^2))/(l(l+1) (r^3)ω)));hh2[r_]:=Z[r]r(1-((a(l(l+1)-2)m (M^2))/(l(l+1) (r^3) ω)));
hh0p=FullSimplify[Solve[EQD4== 0,hh0[r]]][[1]]
rulehh0=Union[hh0p,D[hh0p,r],D[hh0p,{r,2}]];

Equations for Q and Z

EEQ=Series[(EQ8)//.rulehh0,{a,0,1}]//Normal//Simplify;
EEZ=Series[(EQ10)//.rulehh0,{a,0,1}]//Normal//Simplify;
EQs=Solve[{EEQ==0,EEZ==0},{Q''[r],Z''[r]}][[1]]//Simplify;
EEQ1=Series[((f[r])^2)EQs[[1]][[2]],{a,0,1}]//Normal//FullSimplify;
EEZ1=Series[((f[r])^2)EQs[[2]][[2]],{a,0,1}]//Normal//FullSimplify;

Potentials and source terms

DQ=-Coefficient[EEQ1,Q'[r]]//FullSimplify
DZ=-Coefficient[EEZ1,Z'[r]]//FullSimplify
VQ=-(Coefficient[EEQ1,Q[r]])//FullSimplify
VZ=-(Coefficient[EEZ1,Z[r]])//FullSimplify
SQ=Coefficient[EEQ1,Z[r]]//FullSimplify
SZ=Coefficient[EEZ1,Q[r]]//FullSimplify

Check at first order: OK!

ruleTT=Union[rulehh0,EQs,D[EQs,r],D[EQs,{r,2}],{}];

Series[EQ7//.ruleTT,{a,0,1}]//Simplify
Series[EQ8//.ruleTT,{a,0,1}]//Simplify
Series[EQ10//.ruleTT,{a,0,1}]//Simplify
Series[EQD4//.ruleTT,{a,0,1}]//Simplify

Axial equations at first order in ã

EQfinal=((f[r])^2)Q''+f[r]f'[r]Q'[r]+VQ Q[r]-SQ Z[r]
EZfinal=((f[r])^2)Z''+f[r]f'[r]Z'[r]+VZ Z[r]-SZ Q[r]

Axial dipole (l=1)

EQdipole=EQfinal/.l-> 1

(*DumpSave["eqs_Richard/axial_final.mx",{EQfinal,EZfinal,EQdipole}];*)

Massless limit (ã=0)

Clear[hh0,hh1,hh2]

Berndtson-Zerilli transformations

λ=(1/2)(l-1)(l+2);
hh2[r_]:=(1/(((I ω))^2))(((-6M+(3+2λ)r)/(3 r))ψ2[r]+(r-2M)ψ2'[r]+ψ1[r]);
hh1[r_]:=(1/(((I ω))^2))(-((2λ)/(3r))ψ2[r]-((2λ)/3)ψ2'[r]+(2/r)ψ1[r]-ψ1'[r]);
hh0[r_]:=(1/(I ω))(((2λ)/3)ψ2[r]+ψ1[r]);

ψ2 and ψ1 satisfy Regge-Wheeler’s equation [((d^2)/(SubStar[(dr^2)]))+(ω^2)-f(r)(((l(l+1))/(r^2))+(1-(s^2))((2M)/(r^3)))](Subscript[ψ, s])=0 for s=2 and s=1 respectively

R2=ψ2''[r]->-(((2M)/(r^2))(1-((2M)/r))ψ2'[r]+((ω^2)-(1-((2M)/r))(((l(l+1))/(r^2))-((6M)/(r^3))))ψ2[r])(((1-((2M)/r)))^(-2));
R1=ψ1''[r]->-(((2M)/(r^2))(1-((2M)/r))ψ1'[r]+((ω^2)-(1-((2M)/r))(((l(l+1))/(r^2))))ψ1[r])(((1-((2M)/r)))^(-2));
rulenomass=Union[{R2,R1}];
ruleT=Union[rulenomass,D[rulenomass,r],D[rulenomass,{r,2}],{}];

Check:OK!

EQ7//.Union[{a-> 0,μ-> 0},ruleT]//Simplify
EQ8//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQ10//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQD4//.Union[{a->0,μ-> 0},ruleT]//Simplify

Polar Sector

Simplify the equations

These computations can take some minutes.

Eliminating H0, H1, H2, and η0:

H2p=Solve[{EQT==0},{H2[r]}][[1]]//Simplify
rule1=Union[H2p,D[H2p,r],D[H2p,{r,2}],D[H2p,{r,3}],{}];
eeD2=Series[EQD2//.rule1,{a,0,1}]//Normal//Simplify;
H1p=Solve[eeD2==0,H1[r]][[1]]//Simplify
rule2=Union[rule1,H1p,D[H1p,r],D[H1p,{r,2}],D[H1p,{r,3}],{}];
eeD3=Series[EQD3//.rule2,{a,0,1}]//Normal//Simplify;
h0p=Solve[eeD3==0,h0[r]][[1]]//Simplify
rule3=Union[rule2,h0p,D[h0p,r],D[h0p,{r,2}],D[h0p,{r,3}],{}];
ee4=Series[EQ4//.rule3,{a,0,1}]//Normal//Simplify;
H0p=Solve[ee4==0,H0[r]][[1]]//Simplify
rule4=Union[rule3,H0p,D[H0p,r],D[H0p,{r,2}],D[H0p,{r,3}],{}];

Eliminating the third derivatives of K and η1:

ee5=Series[EQ5//.rule4,{a,0,1}]//Normal//Simplify;
KKp3=Solve[ee5==0,KK'''[r]][[1]]//Simplify;
rule5=Union[rule4,KKp3,D[KKp3,r],D[KKp3,{r,2}],{}];
ee1=Series[EQ1//.rule5,{a,0,1}]//Normal//Simplify;
h1p4=Solve[ee1==0,h1''''[r]][[1]]//Simplify;
rule6=Union[rule5,h1p4,D[h1p4,r],D[h1p4,{r,2}],{}];
eeD1=Series[EQD1//.rule6,{a,0,1}]//Normal//Simplify;
h1p3=Solve[eeD1==0,h1'''[r]][[1]]//Simplify;

Second-order ODE’s for K and η1

rulePo={H0[r]-> Series[(H0p[[1]][[2]]),{a,0,1}]//Normal//Simplify,H1[r]-> Series[(H1p[[1]][[2]]//.Union[rule4,KKp3,h1p3]),{a,0,1}]//Normal//Simplify,H2[r]-> Series[(H2p[[1]][[2]]//.rule4),{a,0,1}]//Normal//Simplify,h0[r]-> Series[h0p[[1]][[2]],{a,0,1}]//Normal//Simplify,KK'''[r]-> Series[(KKp3[[1]][[2]]//.Union[h1p3,D[h1p3,r]]),{a,0,1}]//Normal//Simplify,h1'''[r]-> Series[h1p3[[1]][[2]],{a,0,1}]//Normal//Simplify};
ruleT=Union[rulePo,D[rulePo,r],D[rulePo,{r,2}],D[rulePo,{r,3}],{}];

ee3=Series[EQ3//.ruleT,{a,0,1}]//Normal//Simplify;
ee6=Series[EQ6//.ruleT,{a,0,1}]//Normal//Simplify;

These equations are quite formidable. We’ll state them in the next subsection.

EQs=Solve[{ee3==0,ee6==0},{KK''[r],h1''[r]}][[1]]//Simplify;
rule7=Union[ruleT,EQs,{}];

Second-order ODE for G

ee9=Series[EQ9//.rule7,{a,0,1}]//Normal//Simplify;
GGp2=Solve[ee9==0,GG''[r]][[1]]//Simplify;

Check at first order: OK!

rulePofinal={H0[r]-> Series[(rulePo[[1]][[2]]//.EQs),{a,0,1}]//Normal//Simplify,H1[r]-> Series[(rulePo[[2]][[2]]//.EQs),{a,0,1}]//Normal//Simplify,H2[r]-> Series[(rulePo[[3]][[2]]//.EQs),{a,0,1}]//Normal//Simplify,h0[r]-> Series[(rulePo[[4]][[2]]//.EQs),{a,0,1}]//Normal//Simplify,KK''[r]-> Series[EQs[[1]][[2]],{a,0,1}]//Normal//Simplify,h1''[r]-> Series[EQs[[2]][[2]],{a,0,1}]//Normal//Simplify,GGp2[[1]]};
ruleTT=Union[rulePofinal,D[rulePofinal,r],D[rulePofinal,{r,2}],D[rulePofinal,{r,3}],{}];

Series[EQ1//.ruleTT,{a,0,1}]//Simplify
Series[EQ2//.ruleTT,{a,0,1}]//Simplify
Series[EQ3//.ruleTT,{a,0,1}]//Simplify
Series[EQ4//.ruleTT,{a,0,1}]//Simplify
Series[EQ5//.ruleTT,{a,0,1}]//Simplify
Series[EQ6//.ruleTT,{a,0,1}]//Simplify
Series[EQ9//.ruleTT,{a,0,1}]//Simplify
Series[EQD1//.ruleTT,{a,0,1}]//Simplify
Series[EQD2//.ruleTT,{a,0,1}]//Simplify
Series[EQD3//.ruleTT,{a,0,1}]//Simplify
Series[EQT//.ruleTT,{a,0,1}]//Simplify

Polar equations at first order in ã

f[r_]:=1-((2M)/r);
EEK=Series[((f[r])^2)(KK''[r]-EQs[[1]][[2]]),{a,0,1}]//Normal//Simplify;
EEh1=Series[((f[r])^2)(h1''[r]-EQs[[2]][[2]]),{a,0,1}]//Normal//Simplify;
EEGG=Series[((f[r])^2)(GG''[r]-GGp2[[1]][[2]]),{a,0,1}]//Normal//Simplify;

Equation for K

αα1=Coefficient[EEK,KK'[r]]//Simplify;
ββ1=Coefficient[EEK,KK[r]]//Simplify;
γγ1=-(((Coefficient[EEK,h1'[r]]))/(l(l+1)))//Simplify;
δδ1=-(((Coefficient[EEK,h1[r]]))/(l(l+1)))//Simplify;
σσ1=-(((Coefficient[EEK,GG'[r]]))/(l(l+1)(l(l+1)-2)))//Simplify;
ρρ1=-(((Coefficient[EEK,GG[r]]))/(l(l+1)(l(l+1)-2)))//Simplify;
EKfinal=((f[r])^2)KK''[r]+αα1 KK'[r]+ββ1 KK[r]-γγ1 l(l+1) h1'[r]-δδ1 l(l+1) h1[r]-σσ1 l(l+1)(l(l+1)-2) GG'[r]-l(l+1)(l(l+1)-2)ρρ1 GG[r];

EEK

Equation for η1

αα2=Coefficient[EEh1,h1'[r]]//Simplify;
ββ2=Coefficient[EEh1,h1[r]]//Simplify;
γγ2=-Coefficient[EEh1,KK'[r]]//Simplify;
δδ2=-Coefficient[EEh1,KK[r]]//Simplify;
σσ2=-(((Coefficient[EEh1,GG'[r]]))/(l(l+1)(l(l+1)-2)))//Simplify;
ρρ2=-(((Coefficient[EEh1,GG[r]]))/(l(l+1)(l(l+1)-2)))//Simplify;
Eh1final=((f[r])^2)h1''[r]+αα2 h1'[r]+ββ2 h1[r]-γγ2  KK'[r]-δδ2 KK[r]-σσ2 l(l+1)(l(l+1)-2) GG'[r]-l(l+1)(l(l+1)-2)ρρ2 GG[r];

EEh1

Equation for G

αα3=Coefficient[EEGG,GG'[r]]//Simplify;
ββ3=Coefficient[EEGG,GG[r]]//Simplify;
γγ3=-Coefficient[EEGG,KK'[r]]//Simplify;
δδ3=-Coefficient[EEGG,KK[r]]//Simplify;
σσ3=-Coefficient[EEGG,h1'[r]]//Simplify;
ρρ3=-Coefficient[EEGG,h1[r]]//Simplify;
EGfinal=((f[r])^2)GG''[r]+αα3 GG'[r]+ββ3 GG[r]-γγ3  KK'[r]-δδ3 KK[r]-σσ3 h1'[r]-ρρ3 h1[r];

EEGG

(*DumpSave["eqs_Richard/polar_final.mx", {EKfinal, Eh1final, EGfinal}];*)

Asymptotic behavior

(*<<"eqs_Richard/polar_final.mx";*)

Define (Subscript[Z, 1]), (Subscript[Z, 2]), (Subscript[Z, 3]) such that they have the usual boundary counditions and the source terms vanish at the horizon and at infinity.

KK[r_]:=((Z1[r])/r);h1[r_]:= Z2[r]((f[r])^(-1));GG[r_]:=((Z3[r])/r);
EZ1=Simplify[EKfinal r];
EZ2=Simplify[Eh1final f[r]];
EZ3=Simplify[EGfinal r];

Asymptotic behavior of Z1

Limit[Coefficient[EZ1,Z1[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ1,Z1[r]],r->2M]//Simplify
Limit[Coefficient[EZ1,Z2'[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ1,Z2'[r]],r->2M]//Simplify
Limit[Coefficient[EZ1,Z2[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ1,Z2[r]],r->2M]//Simplify
Limit[Coefficient[EZ1,Z3'[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ1,Z3'[r]],r->2M]//Simplify
Limit[Coefficient[EZ1,Z3[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ1,Z3[r]],r->2M]//Simplify

Asymptotic behavior of Z2

Limit[Coefficient[EZ2,Z2[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ2,Z2[r]],r->2M]//Simplify
Limit[Coefficient[EZ2,Z1'[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ2,Z1'[r]],r->2M]//Simplify
Limit[Coefficient[EZ2,Z1[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ2,Z1[r]],r->2M]//Simplify
Limit[Coefficient[EZ2,Z3'[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ2,Z3'[r]],r->2M]//Simplify
Limit[Coefficient[EZ2,Z3[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ2,Z3[r]],r->2M]//Simplify

Asymptotic behavior of Z3

Limit[Coefficient[EZ3,Z3[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ3,Z3[r]],r->2M]//Simplify
Limit[Coefficient[EZ3,Z1'[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ3,Z1'[r]],r->2M]//Simplify
Limit[Coefficient[EZ3,Z1[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ3,Z1[r]],r->2M]//Simplify
Limit[Coefficient[EZ3,Z2'[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ3,Z2'[r]],r->2M]//Simplify
Limit[Coefficient[EZ3,Z2[r]],r->Infinity]//Simplify
Limit[Coefficient[EZ3,Z2[r]],r->2M]//Simplify

Polar Dipole (l=1)

Clear[GG,KK,h1]

(*<<"eqs_Richard/polar_final.mx";*)

EKdipole=EKfinal/.l->1//Simplify
Eh1dipole=Eh1final/.l->1//Simplify

The function G is not defined for l=1

Coefficient[EKdipole,GG[r]]//Simplify
Coefficient[EKdipole,GG'[r]]//Simplify
Coefficient[Eh1dipole,GG[r]]//Simplify
Coefficient[Eh1dipole,GG'[r]]//Simplify

Monopole (l=0)

Only one equation for the monopole.

EKmono=EKfinal//.{l->0,m-> 0}//Simplify

Clear[KK,H1]

Monopole equation can be written in simplified form.

eeD1=EQD1//.{l-> 0,m->0}//Simplify;
eeT=EQT//.{l-> 0,m->0}//Simplify;
Hp=Solve[{eeT==0,eeD1==0},{H2[r],H0[r]}][[1]]//Simplify
rule1=Union[Hp,D[Hp,r],D[Hp,{r,2}],{}];

K and H1 can be written as a linear comination of a single function ψ0.

KK[r_]:=(((2 r-4M)/(r^4))-(((μ^2) (-10 M+6 r+(r^3) (μ^2)))/(2 M r+(r^4) (μ^2))))ψ0[r]+((-2r+4M)/(r^3))ψ0'[r];
H1[r_]:=(((2I ω(r-M))/((2M-r)(r^2)))+((6 I r (μ^2) ω)/(2 M+(r^3) (μ^2))))ψ0[r]+((2I ω)/r)ψ0'[r];

ee2=EQ2//.Union[{l-> 0,m->0},rule1]//Simplify;
ee4=EQ4//.Union[{l-> 0,m->0},rule1]//Simplify;
ψp=Solve[{ee2==0,ee4==0},{ψ0''[r],ψ0'''[r]}][[1]]//Simplify
rule2=Union[rule1,{ψp[[1]]},{D[ψp[[1]],r]},{D[ψp[[1]],{r,2}]},{}];

Check: OK!

ee1=EQ1//.Union[{l-> 0,m->0},rule2]//Simplify
ee2=EQ2//.Union[{l-> 0,m->0},rule2]//Simplify
ee3=EQ3//.Union[{l-> 0,m->0},rule2]//Simplify
ee4=EQ4//.Union[{l-> 0,m->0},rule2]//Simplify
eD1=EQD1//.Union[{l-> 0,m->0},rule2]//Simplify
eD2=EQD2//.Union[{l-> 0,m->0},rule2]//Simplify
eT=EQT//.Union[{l-> 0,m->0},rule2]//Simplify

Monopole Equation

EQmono=(((1-((2M)/r)))^2)(ψ0''[r]-ψp[[1]][[2]])//Simplify

Potential

Dψ0=Coefficient[EQmono,ψ0'[r]]//FullSimplify
Vψ0=Coefficient[EQmono,ψ0[r]]-((ω^2)-(1-((2M)/r))(((2M)/(r^3))+(μ^2)))//FullSimplify

It reduces to the scalar equation with l=0 in the massless limit.

EQmono/.μ-> 0//Simplify

Massless limit (ã=0)

Clear[H1,KK,GG,h1]

Berndtson-Zerilli transformations

λ=(1/2)(l-1)(l+2);
H0[r_]:=-λ(1+λ)*(18(M^4)+3(-3+4λ)(M^3)r+(((I ω))^2)(λ^2)(r^6)-3λ M (r^3)(1+λ-2(((I ω))^2)(r^2))+(M^2)(6(λ^2)(r^2)+9(((I ω))^2)(r^4)))/(3 (((I ω))^2) (2M-r)(r^4)(((3M+λ r))^2))ψ2[r]-λ(1+λ)*((M(-3M+(3+λ)r))/(3 (((I ω))^2) (r^3)(3M+λ r)))ψ2'[r]+((4(1+λ) I ω )/(2M-r))ψ1[r]+((4(1+λ) M)/(I ω (r^3)))ψ1'[r]+((2(-2(M^2)+M r+(((I ω))^2)(r^4)))/((2M-r)(r^4)))ψ0[r]+((2M)/(r^3))ψ0'[r];
H1[r_]:=-(((3(M^2)+3M λ r-λ (r^2))λ(1+λ))/(3 I ω (2M-r)(r^2)(3M+λ r)))ψ2[r]-((λ(1+λ))/(3 I ω r))ψ2'[r]+(((2+2λ))/(2M r-(r^2)))ψ1[r]+((4(1+λ))/r)ψ1'[r]+((2I ω(r-M))/((2M-r)(r^2)))ψ0[r]+((2I ω)/r)ψ0'[r];
H2[r_]:=-λ(1+λ)*(-54(M^4)+3(9-16λ)(M^3)r+(λ^2)(r^4)(2+2λ+(((I ω))^2)(r^2))+9(M^2)(r^2)(2λ-2(λ^2)+(((I ω))^2)(r^2))+λ M (r^3)(3+5λ-4(λ^2)+6(((I ω))^2)(r^2)))/(3 (((I ω))^2) (2M-r)(r^4)(((3M+λ r))^2))ψ2[r]+λ(1+λ)*(((-9(M^2)+(3-5λ)M r+2λ (r^2)))/(3 (((I ω))^2) (r^3)(3M+λ r)))ψ2'[r]+((4(1+λ) (-4(1+λ) M+r(2+2λ+(((I ω))^2)(r^2))) )/((2M-r)(r^3)I ω))ψ1[r]-((4(1+λ)(M-r))/(I ω (r^3)))ψ1'[r]+((2(6(M^2)-(11+4λ)M r+(r^2)(4+2λ+(((I ω))^2)(r^2))))/((2M-r)(r^4)))ψ0[r]+((-6M+4r)/(r^3))ψ0'[r];
KK[r_]:=-λ(1+λ)*((6(M^2)+3M λ r+λ(1+λ)(r^2))/(3  (((I ω))^2)(r^4)(3M+λ r)))ψ2[r]+λ(λ+1)*((2M-r)/(3 (((I ω))^2) (r^3)))ψ2'[r]+((4(((1+λ))^2))/(I ω (r^3)))ψ1[r]+((2(1+λ)(2M-r))/(I ω (r^3)))ψ1'[r]+((2(2+λ )r-4M)/(r^4))ψ0[r]+((-2r+4M)/(r^3))ψ0'[r];
h0[r_]:=-λ((6(M^2)+3M λ r+λ(1+λ)(r^2))/(3  I ω (r^2)(3M+λ r)))ψ2[r]+λ((2M-r)/(3  I ω r))ψ2'[r]+((4(1+λ))/r)ψ1[r]+(1-((2M)/r))ψ1'[r]+((2I ω)/r)ψ0[r];
h1[r_]:=-λ(9(((I ω))^2)(M^2)(r^3)+2(λ^3)(r^2)(r-2M)+(λ^2)r(-12(M^2)+2M r+2(r^2)+(((I ω))^2)(r^4))+3λ M(-4(M^2)+(r^2)+2(((I ω))^2)(r^4)))/(3  (((I ω ))^2)(2M-r)(r^2)(((3M+λ r))^2))ψ2[r]+λ(((3+λ)M+λ(2+λ)r)/(3 (((I ω))^2) r(3M+λ r)))ψ2'[r]+((8M+8λ M-4r-4λ r+(((I ω))^2)(r^3))/(I ω (r^2)(2 M - r)))ψ1[r]-((2(2M+r+2λ r))/(I ω (r^2)))ψ1'[r]+(4/(r^2))ψ0[r]-(2/r)ψ0'[r];
GG[r_]:=(4(λ^3)(r^2)+(λ^4)(r^2)+27(((I ω))^2)(M^2)(r^2)+9λ M(M+2(((I ω))^2)(r^3))+3(λ^2)((M^2)+M r+(r^2)+(((I ω))^2)(r^4)))/(6 (((I ω))^2) (r^3)(((3M+λ r))^2))ψ2[r]-λ(((3+2λ)(2M-r))/(6 (((I ω))^2) (r^2)(3M+λ r)))ψ2'[r]-((2(1+λ))/(I ω (r^3)))ψ1[r]+((2M-r)/(I ω (r^3)))ψ1'[r]-(1/(r^3))ψ0[r];

ψ0 and ψ1 satisfy Regge-Wheeler’s equation [((d^2)/(SubStar[(dr^2)]))+(ω^2)-f(r)(((l(l+1))/(r^2))+(1-(s^2))((2M)/(r^3)))](Subscript[ψ, s])=0, for s=0 and s=1, respectively, and ψ2 satisfy Zerilli’s equation  [((d^2)/(SubStar[(dr^2)]))+(ω^2)-f(r)(((18(M^3)+18λ (M^2)r+6(λ^2)M (r^2)+2(λ^2)(1+λ)(r^3)))/((r^3)(((3M+λ r))^2))))](Subscript[ψ, 2])=0.

LZ=ψ2''[r]->-(((2M)/(r^2))(1-((2M)/r))ψ2'[r]+((ω^2)+((2(2M-r)(9(M^3)+9λ (M^2)r+3(λ^2)M (r^2)+(λ^2)(1+λ)(r^3)))/((r^4)(((3M+λ r))^2))))ψ2[r])(((1-((2M)/r)))^(-2));
R1=ψ1''[r]->-(((2M)/(r^2))(1-((2M)/r))ψ1'[r]+((ω^2)-(1-((2M)/r))(((l(l+1))/(r^2))))ψ1[r])(((1-((2M)/r)))^(-2));
R0=ψ0''[r]->-(((2M)/(r^2))(1-((2M)/r))ψ0'[r]+((ω^2)-(1-((2M)/r))(((l(l+1))/(r^2))+((2M)/(r^3))))ψ0[r])(((1-((2M)/r)))^(-2));
rulenomass=Union[{LZ,R1,R0}];
ruleT=Union[rulenomass,D[rulenomass,r],D[rulenomass,{r,2}],{}];

Check:OK!

EQ1//.Union[{a-> 0,μ-> 0},ruleT]//Simplify
EQ2//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQ3//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQ4//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQ5//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQ6//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQ9//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQT//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQD1//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQD2//.Union[{a->0,μ-> 0},ruleT]//Simplify
EQD3//.Union[{a->0,μ-> 0},ruleT]//Simplify

Continued fractions:

M=1;
q=sign(Sqrt[(μ^2)-(ω^2)])
sign=1 for bound states
sign=-1 for QNM

Axial Dipole (l=1)

α[n_]:=(n+1)(n+1-4I ω)+I m a(n+1);
β[n_]:= 2(1-(n^2)-n)+(q^(-1))(1+2 n-4 I ω) (ω^2)-3 q(1+2 n-4 I ω)+(2n +1)4 I ω-4 (q^2)++12 (ω^2)+ m a (-I (1+2n)-3I q+(3/ω)(1/4)-2 ω+ ((I (ω^2))/q));
γ[n_]:=(q^(-2))(n q+(q^2)-3q-2I q ω-(ω^2))(n q+(q^2)+3q-2I q ω-(ω^2))+m (I (n+q)-(9/ω)+2 ω-((I (ω^2))/q)) a;
δ[n_]:=((45 m a)/(2 ω));
ρ[n_]:=-((21 m a)/ω);
σ[n_]:=((27 m a)/(4 ω));

Axial (l≥2)

α11[n_]:=(n+1)(n+1-4I ω)+I m a(n+1);
α12[n_]:=0;
α21[n_]:=0;
α22[n_]:=α11[n];
(***)
β11[n_]:=2-Λ-4 (q^2)+(q^(-1))(1+2 n-4 I ω) (ω^2)-3 q(1+2 n-4 I ω)+ 2(1-(n^2)-n)+(2n +1)4 I ω+12 (ω^2)+ m a (-I (1+2n)-3I q+((3  (2+Λ))/(8 Λ ω))-2 ω+ ((I (ω^2))/q));
β12[n_]:=(Λ-2)-(((m (-2+Λ) (-1+6 (ω^2))) a)/(2 (Λ ω)));
β21[n_]:=-2-(((m (-10+12 (q^2)+3 Λ+12 (ω^2))) a)/(4 (Λ ω)));
β22[n_]:=β11[n]-3(1+((a m)/(2 Λ ω)));
(***)
γ11[n_]:=(q^(-2))((q^4)+(q^3) (2n-4 I ω)-q(2n-4 I ω) (ω^2)+(ω^4)-(q^2)(9-(n^2)+4n I ω+6(ω^2)))+m (I (n+q)-((9 (2+Λ))/( 2Λ ω))+2 ω-((I (ω^2))/q)) a;
γ12[n_]:=(6-3 Λ)-((3 (m (-2+Λ) (6+Λ-4 (ω^2))) a)/(4 (Λ ω)));
γ21[n_]:=((3 m (-10+4 (q^2)+3 Λ+4 (ω^2)) a)/(4 Λ ω));
γ22[n_]:=γ11[n]+9(1-((a m (6+Λ))/(4 Λ ω)));
(***)
δ11[n_]:=((45 m (2+Λ) a)/(4 Λ ω));
δ12[n_]:=((3 m (-2+Λ) (14+3 Λ) a)/(4 Λ ω));
δ21[n_]:=((m (30-9 Λ) a)/(4 Λ ω));
δ22[n_]:=((9 m (-2+Λ) a)/(2 Λ ω));
(***)
ρ11[n_]:=-((21 (m (2+Λ)) a)/(2 (Λ ω)));
ρ12[n_]:=-(((m (-2+Λ) (38+9 Λ)) a)/(4 (Λ ω)));
ρ21[n_]:=((m (-10+3 Λ) a)/(4 Λ ω));
ρ22[n_]:=-((15 (m (-2+Λ)) a)/(4 (Λ ω)));
(***)
σ11[n_]:=((27 m (2+Λ) a)/(8 Λ ω));
σ12[n_]:=((3 m (-2+Λ) (4+Λ) a)/(4 Λ ω));
σ21[n_]:=0;
σ22[n_]:=((9 m (-2+Λ) a)/(8 Λ ω));
(******************************)
α[n_]:={{α11[n],α12[n]},{α21[n],α22[n]}};
β[n_]:={{β11[n],β12[n]},{β21[n],β22[n]}};
γ[n_]:={{γ11[n],γ12[n]},{γ21[n],γ22[n]}};
δ[n_]:={{δ11[n],δ12[n]},{δ21[n],δ22[n]}};
ρ[n_]:={{ρ11[n],ρ12[n]},{ρ21[n],ρ22[n]}};
σ[n_]:={{σ11[n],σ12[n]},{σ21[n],σ22[n]}};