Gravito-Electromagnetic Perturbations of a Kerr-Newman Black Hole in the Slow-Rotation Limit
Reference: “Gravito-Electromagnetic Perturbations of Kerr-Newman Black Holes: Stability and Isospectrality in the Slow-Rotation Limit,” P .Pani, E. Berti, L. Gualtieir (2013).

SetDirectory[NotebookDirectory[]];

Load equations without couplings to l±1

The equations are stored in an external file. For simplicity we report them here.

(*<<"pert_eqs.mx";*)

EQ1=-2 (Q^2) (r^4) H00[r]-((((Q^2)+r (-2 M+r)))^2) (4 (Q^2)-r (4 M+(2+l+(l^2)) r)) H20[r]+(1/(l (1+l) M))2 I J m (r^2) ϵa (l (1+l) (-M (Q^2)+2 ((M^2)+(Q^2)) r-3 M (r^2)) H10[r]+I l (1+l) ((Q^2)-2 M r) ((Q^2)-2 M r+(r^2)) ω H20[r]+l (1+l) (r (I r ((Q^2)-2 M r) ω K0[r]+2 Q u2[r])+((Q^2)-2 M r) ((Q^2)-2 M r+(r^2)) H10'[r])-2 Q ((Q^2)-2 M r+(r^2)) u3'[r])+r ((-2+l+(l^2)) (r^3) ((Q^2)+r (-2 M+r)) K0[r]+2 (-2 Q r ((Q^2)+r (-2 M+r)) u1[r]+2 I Q (r^4) ω u2[r]+((Q^2)+r (-2 M+r)) (((((Q^2)+r (-2 M+r)))^2) H20'[r]+(r^2) ((-2 (Q^2)+5 M r-3 (r^2)) K0'[r]+2 Q u1'[r]-r ((Q^2)-2 M r+(r^2)) K0''[r]))));
EQ2=(1/(l (1+l) M))I J m ϵa (l (1+l) (r^2) ((Q^2)-M r) (2 (Q^2)+r (-4 M+r)) H00[r]+((Q^2)+r (-2 M+r)) (-I l (1+l) (r^3) (-(Q^2)+2 M r) ω H10[r]+l (1+l) ((Q^2)-M r) ((Q^2)-2 M r+(r^2)) H20[r]+(r^2) (4 l (1+l) Q u1[r]+4 I Q r ω u3[r]+l (1+l) r ((Q^2)-2 M r) H00'[r])))+(r^2) ((Q^2)+r (-2 M+r)) (l (1+l) r ((Q^2)+r (-2 M+r)) H10[r]-2 I ω (((((Q^2)+r (-2 M+r)))^2) H20[r]-(r^2) ((2 (Q^2)+r (-3 M+r)) K0[r]+r ((Q^2)+r (-2 M+r)) K0'[r])));
EQ3=(r^2) (-2 (Q^2)+r (4 M+l (1+l) r)) H00[r]-4 I (r^3) ((Q^2)+r (-2 M+r)) ω H10[r]-2 ((((Q^2)+r (-2 M+r)))^2) H20[r]-(r^2) ((-2+l+(l^2)) ((Q^2)+r (-2 M+r))-2 (r^4) (ω^2)) K0[r]+2 (2 Q ((Q^2)+r (-2 M+r)) u1[r]-r (2 I Q (r^2) ω u2[r]+((Q^2)+r (-2 M+r)) ((r^2) H00'[r]+(M-r) r K0'[r]+2 Q u1'[r])))-((2 I J m ϵa (-l (1+l) M ((Q^2)+r (-2 M+r)) H10[r]+I l (1+l) r (r ((Q^2)-2 M r) ω K0[r]+2 I Q u2[r])+2 Q ((Q^2)+r (-2 M+r)) u3'[r]))/(l (1+l) M));
EQ4=(1/r)(1/(l (1+l) M)2 J m ((Q^2)+r (-2 M+r)) ϵa (2 I l (1+l) r ((Q^2)-M r) (2 (Q^2)+r (-4 M+r)) H10[r]+((Q^2)-2 M r) (l (1+l) r (2 r ((Q^2)-2 M r+(r^2)) ω H20[r]+(r^3) ω K0[r]+4 I Q u2[r]-2 I r ((Q^2)-2 M r+(r^2)) H10'[r])-4 I Q ((Q^2)+r (-2 M+r)) u3'[r]))+r ((r^4) (-4 M (Q^2)+4 (M^2) r+(4+l+(l^2)) (Q^2) r-2 (2+l+(l^2)) M (r^2)+l (1+l) (r^3)) H00[r]+((Q^2)+r (-2 M+r)) (4 I (M-r) (r^5) ω H10[r]-((Q^2)+r (-2 M+r)) (4 M (Q^2)-4 (M^2) r-2 (-2+l+(l^2)) M (r^2)+r ((-4+l+(l^2)) (Q^2)+(r^2) (l+(l^2)-2 (r^2) (ω^2)))) H20[r]+2 r ((r^6) (ω^2) K0[r]-4 Q ((Q^2)-2 M r+(r^2)) u1[r]+4 I Q (r^3) ω u2[r]-2 (Q^2) (r^3) H00'[r]+3 M (r^4) H00'[r]-(r^5) H00'[r]-2 I (Q^2) (r^4) ω H10'[r]+4 I M (r^5) ω H10'[r]-2 I (r^6) ω H10'[r]+M (Q^4) H20'[r]-4 (M^2) (Q^2) r H20'[r]-(Q^4) r H20'[r]+4 (M^3) (r^2) H20'[r]+6 M (Q^2) (r^2) H20'[r]-8 (M^2) (r^3) H20'[r]-2 (Q^2) (r^3) H20'[r]+5 M (r^4) H20'[r]-(r^5) H20'[r]-2 M (Q^2) (r^2) K0'[r]+4 (M^2) (r^3) K0'[r]+2 (Q^2) (r^3) K0'[r]-6 M (r^4) K0'[r]+2 (r^5) K0'[r]+4 (Q^3) r u1'[r]-8 M Q (r^2) u1'[r]+4 Q (r^3) u1'[r]-(Q^2) (r^4) H00''[r]+2 M (r^5) H00''[r]-(r^6) H00''[r]+(r^2) ((((Q^2)-2 M r+(r^2)))^2) K0''[r]))));
EQ5=(1/(2 (r^4)))(2 (r^4) (l (1+l) r (-2 ((Q^2)-M r) H10[r]+I r ω (((Q^2)-2 M r+(r^2)) H20[r]+(r^2) K0[r])-4 Q u2[r]+r ((Q^2)+r (-2 M+r)) H10'[r])+4 Q ((Q^2)+r (-2 M+r)) u3'[r])+1/MI J m ϵa ((((-2+l+(l^2)) (r^6) ((Q^2)-2 M r) H00[r])/((Q^2)+r (-2 M+r)))+((2 (r^2) (l (1+l) (((Q^2)-2 M r) (((Q^2)+r (-2 M+r)) H20[r]-2 (r^2) K0[r])-8 Q (r^2) u1[r])-8 I Q (r^3) ω u3[r]))/(l (1+l)))+1/((((Q^2)+r (-2 M+r)))^2)(-(r^4) (8 (Q^6)+(Q^4) r (-40 M+(12+l+(l^2)) r)-2 M (r^3) (16 (M^2)-2 (5+l+(l^2)) M r+l (1+l) (r^2))+(Q^2) (r^2) (64 (M^2)-4 (8+l+(l^2)) M r+l (1+l) (r^2))) H00[r]+2 ((Q^2)+r (-2 M+r)) (2 I (M-r) (r^6) (-(Q^2)+2 M r) ω H10[r]+((Q^2)+r (-2 M+r)) (-4 (Q^6)+(Q^4) r (20 M+(-6+l+(l^2)) r)+2 M (r^3) (8 (M^2)+(-5+2 l (1+l)) M r-l (1+l) (r^2)+(r^4) (ω^2))+(Q^2) (r^2) (-32 (M^2)-4 (-4+l+(l^2)) M r+(r^2) (l+(l^2)-(r^2) (ω^2)))) H20[r]+(r^4) ((Q^2)-2 M r) (l (1+l) ((Q^2)+r (-2 M+r))-(r^4) (ω^2)) K0[r]-4 Q (r^4) ((Q^2)+r (-2 M+r)) u1[r]+4 I Q (r^7) ω u2[r]-(Q^2) (r^7) H00'[r]+M (r^8) H00'[r]+2 I (Q^4) (r^6) ω H10'[r]-8 I M (Q^2) (r^7) ω H10'[r]+8 I (M^2) (r^8) ω H10'[r]+2 I (Q^2) (r^8) ω H10'[r]-4 I M (r^9) ω H10'[r]+2 (Q^8) r H20'[r]-16 M (Q^6) (r^2) H20'[r]+48 (M^2) (Q^4) (r^3) H20'[r]+7 (Q^6) (r^3) H20'[r]-64 (M^3) (Q^2) (r^4) H20'[r]-41 M (Q^4) (r^4) H20'[r]+32 (M^4) (r^5) H20'[r]+80 (M^2) (Q^2) (r^5) H20'[r]+8 (Q^4) (r^5) H20'[r]-52 (M^3) (r^6) H20'[r]-30 M (Q^2) (r^6) H20'[r]+28 (M^2) (r^7) H20'[r]+3 (Q^2) (r^7) H20'[r]-5 M (r^8) H20'[r]-4 (Q^6) (r^3) K0'[r]+26 M (Q^4) (r^4) K0'[r]-56 (M^2) (Q^2) (r^5) K0'[r]-10 (Q^4) (r^5) K0'[r]+40 (M^3) (r^6) K0'[r]+42 M (Q^2) (r^6) K0'[r]-44 (M^2) (r^7) K0'[r]-6 (Q^2) (r^7) K0'[r]+12 M (r^8) K0'[r]+4 (Q^3) (r^5) u1'[r]-8 M Q (r^6) u1'[r]+4 Q (r^7) u1'[r]+(Q^4) (r^6) H00''[r]-4 M (Q^2) (r^7) H00''[r]+4 (M^2) (r^8) H00''[r]+(Q^2) (r^8) H00''[r]-2 M (r^9) H00''[r]+2 (r^4) (-(Q^2)+2 M r) ((((Q^2)+r (-2 M+r)))^2) K0''[r]))));
EQ6=(1/(2 r))(2 r (l (1+l) (M-r) (r^4) H00[r]-((Q^2)+r (-2 M+r)) (-I l (1+l) (r^4) ω H10[r]-l (1+l) (M-r) (-(Q^2)+(2 M-r) r) H20[r]+r (4 l (1+l) Q u1[r]+r (4 I Q ω u3[r]+l (1+l) (-(r^2) H00'[r]+((Q^2)-2 M r+(r^2)) K0'[r])))))+1/MI J m ϵa (2 I (r^4) (2 (Q^2)-3 M r) ω H00[r]+1/(l (1+l))2 ((Q^2)+r (-2 M+r)) ((l^2) (((1+l))^2) r ((Q^2)-2 M r) H10[r]-I (l (1+l) (2 (Q^2)-3 M r) ((Q^2)-2 M r+(r^2)) ω H20[r]+l (1+l) r (2 r (-2 (Q^2)+3 M r) ω K0[r]-8 I Q u2[r]+(r^2) (-(Q^2)+2 M r) ω K0'[r])+8 I Q ((Q^2)-2 M r+(r^2)) u3'[r]))));
EQ7=-(1/(r^3))(((I J m ϵa (-((I (-2+l+(l^2)) (r^5) (-(Q^2)+2 M r) ω h00[r])/((Q^2)+r (-2 M+r)))+((2 (l (1+l) (2 (Q^4)+(Q^2) r (-12 M+l (1+l) r)+M (r^2) (16 M-(6+l+(l^2)) r)) h10[r]+4 I Q (r^4) ω u4[r]+l (1+l) r (-(Q^2)+2 M r) ((Q^2)+r (-2 M+r)) h10'[r]))/(l (1+l)))))/M)-l (1+l) (r^3) (-(2 (Q^2)+r (-4 M+l (1+l) r)) h00[r]+((((Q^2)+r (-2 M+r)) (4 Q u4'[r]+l (1+l) r (I ω (2 h10[r]+r h10'[r])+r h00''[r])))/(l (1+l)))));
EQ8=(1/(l (1+l) M r))(-2 I l (1+l) (J m (6 M r+l (1+l) ((Q^2)-2 M r)) ϵa+l (1+l) M (r^4) ω) h00[r]+(l^2) (((1+l))^2) r ((-2+l+(l^2)) M ((Q^2)+r (-2 M+r))-2 J m ((Q^2)-2 M r) ϵa ω-M (r^4) (ω^2)) h10[r]+I (4 l (1+l) Q r (-J m ϵa+M (r^2) ω) u4[r]+(l^2) (((1+l))^2) r (J m ((Q^2)-2 M r) ϵa+M (r^4) ω) h00'[r]-8 J m Q ((Q^2)-2 M r+(r^2)) ϵa u4'[r]));
EQ9=(((l^2) (((1+l))^2) (r^2) ((r^4) H00[r]-((((Q^2)+r (-2 M+r)))^2) H20[r])+((2 I J m ϵa (l (1+l) r (I (r^3) ((Q^2)-2 M r) ω K0[r]-4 Q ((Q^2)+r (-2 M+r)) u2[r])+4 Q ((((Q^2)+r (-2 M+r)))^2) u3'[r]))/M))/(l (1+l) (r^2)));
EQ10=(1/r)(l (1+l) r (I (r^5) ω h00[r]+((Q^2)+r (-2 M+r)) (-2 ((Q^2)-M r) h10[r]+r ((Q^2)-2 M r+(r^2)) h10'[r]))-((I J m ϵa (-(8 (Q^4)+2 M (r^2) (12 M-(6+l+(l^2)) r)+(Q^2) r (-28 M+(8+l+(l^2)) r)) h00[r]+((2 ((Q^2)+r (-2 M+r)) (l (1+l) r (2 (Q^2)-3 M r) (I ω h10[r]+h00'[r])+2 Q ((Q^2)+r (-2 M+r)) u4'[r]))/(l (1+l)))))/M));
(*MAXWELL*)
EQM1=-2 Q (r^4) ((Q^2)-M r) H00[r]-2 Q ((Q^2)-M r) ((((Q^2)+r (-2 M+r)))^2) H20[r]-2 l (1+l) (r^4) ((Q^2)+r (-2 M+r)) u1[r]+(1/(l (1+l) M))2 I J m r ϵa (l (1+l) Q ((((Q^2)+r (-2 M+r)))^2) H10[r]+r (l (1+l) r (-I r ((Q^2)-2 M r) ω u1[r]-2 ((Q^2)-M r) u2[r])+(r^3) ((Q^2)-2 M r) (ω^2) u3[r]+2 ((Q^2)-M r) ((Q^2)+r (-2 M+r)) u3'[r]))+r (2 I (r^4) (3 (Q^2)+r (-4 M+r)) ω u2[r]+((Q^2)+r (-2 M+r)) (-2 I (r^4) ω u3[r]-Q (r^4) H00'[r]+(Q^5) H20'[r]-4 M (Q^3) r H20'[r]+4 (M^2) Q (r^2) H20'[r]+2 (Q^3) (r^2) H20'[r]-4 M Q (r^3) H20'[r]+Q (r^4) H20'[r]-2 (Q^3) (r^2) K0'[r]+4 M Q (r^3) K0'[r]-2 Q (r^4) K0'[r]+2 I (r^5) ω u2'[r]+2 (r^3) ((Q^2)-2 M r+(r^2)) u1''[r]));
EQM2=(r^2) (-Q (r^4) ω H00[r]+Q ((((Q^2)+r (-2 M+r)))^2) ω H20[r]+2 r (-Q r ((Q^2)+r (-2 M+r)) ω K0[r]-r ((Q^2)-2 M r+(r^2)) ω u1[r]-I (l (1+l) ((Q^2)-2 M r+(r^2))-(r^4) (ω^2)) u2[r]+(r^2) ((Q^2)-2 M r+(r^2)) ω u1'[r])+2 I ((((Q^2)+r (-2 M+r)))^2) u3'[r])-((J m ϵa (Q (r^2) (2 (Q^2)+r (-4 M+r)) H00[r]+Q ((((Q^2)+r (-2 M+r)))^2) H20[r]+((2 ((Q^2)-2 M r) (l (1+l) ((Q^2)+r (-2 M+r)) u1[r]+r (-2 I l (1+l) (r^2) ω u2[r]-((Q^2)-2 M r+(r^2)) (l (1+l) u1'[r]-I r ω u3'[r]))))/(l (1+l)))))/M);
EQM3=2 I J m Q (r^4) (3 (Q^2)+r (-5 M+r)) ϵa H00[r]+2 J m Q (r^5) ((Q^2)-2 M r+(r^2)) ϵa ω H10[r]+2 I J m (Q^7) ϵa H20[r]-10 I J m M (Q^5) r ϵa H20[r]+16 I J m (M^2) (Q^3) (r^2) ϵa H20[r]+6 I J m (Q^5) (r^2) ϵa H20[r]-8 I J m (M^3) Q (r^3) ϵa H20[r]-20 I J m M (Q^3) (r^3) ϵa H20[r]+16 I J m (M^2) Q (r^4) ϵa H20[r]+6 I J m (Q^3) (r^4) ϵa H20[r]-10 I J m M Q (r^5) ϵa H20[r]+2 I J m Q (r^6) ϵa H20[r]+8 I J m (Q^4) (r^2) ϵa u1[r]-28 I J m M (Q^2) (r^3) ϵa u1[r]+24 I J m (M^2) (r^4) ϵa u1[r]+8 I J m (Q^2) (r^4) ϵa u1[r]+2 I J l m (Q^2) (r^4) ϵa u1[r]+2 I J (l^2) m (Q^2) (r^4) ϵa u1[r]-12 I J m M (r^5) ϵa u1[r]-4 I J l m M (r^5) ϵa u1[r]-4 I J (l^2) m M (r^5) ϵa u1[r]+2 I l M (r^8) ω u1[r]+2 I (l^2) M (r^8) ω u1[r]-2 l M (Q^2) (r^5) u2[r]-2 (l^2) M (Q^2) (r^5) u2[r]+4 l (M^2) (r^6) u2[r]+4 (l^2) (M^2) (r^6) u2[r]-2 l M (r^7) u2[r]-2 (l^2) M (r^7) u2[r]+8 J m (Q^2) (r^5) ϵa ω u2[r]-12 J m M (r^6) ϵa ω u2[r]-2 J m (Q^2) (r^5) ϵa ω u3[r]+4 J m M (r^6) ϵa ω u3[r]-2 M (r^9) (ω^2) u3[r]-I J m (Q^3) (r^5) ϵa H00'[r]+2 I J m M Q (r^6) ϵa H00'[r]-I J m Q (r^7) ϵa H00'[r]-I J m (Q^7) r ϵa H20'[r]+6 I J m M (Q^5) (r^2) ϵa H20'[r]-12 I J m (M^2) (Q^3) (r^3) ϵa H20'[r]-3 I J m (Q^5) (r^3) ϵa H20'[r]+8 I J m (M^3) Q (r^4) ϵa H20'[r]+12 I J m M (Q^3) (r^4) ϵa H20'[r]-12 I J m (M^2) Q (r^5) ϵa H20'[r]-3 I J m (Q^3) (r^5) ϵa H20'[r]+6 I J m M Q (r^6) ϵa H20'[r]-I J m Q (r^7) ϵa H20'[r]+2 I J m (Q^5) (r^3) ϵa K0'[r]-8 I J m M (Q^3) (r^4) ϵa K0'[r]+8 I J m (M^2) Q (r^5) ϵa K0'[r]+2 I J m (Q^3) (r^5) ϵa K0'[r]-4 I J m M Q (r^6) ϵa K0'[r]-8 I J m (Q^4) (r^3) ϵa u1'[r]+28 I J m M (Q^2) (r^4) ϵa u1'[r]-24 I J m (M^2) (r^5) ϵa u1'[r]-8 I J m (Q^2) (r^5) ϵa u1'[r]+12 I J m M (r^6) ϵa u1'[r]+2 l M (Q^2) (r^6) u2'[r]+2 (l^2) M (Q^2) (r^6) u2'[r]-4 l (M^2) (r^7) u2'[r]-4 (l^2) (M^2) (r^7) u2'[r]+2 l M (r^8) u2'[r]+2 (l^2) M (r^8) u2'[r]+4 M (Q^4) (r^4) u3'[r]-12 (M^2) (Q^2) (r^5) u3'[r]+8 (M^3) (r^6) u3'[r]+4 M (Q^2) (r^6) u3'[r]-4 (M^2) (r^7) u3'[r]-2 M (r^5) ((((Q^2)-2 M r+(r^2)))^2) u3''[r];
EQM4=-2 Q (r^3) (l (1+l) M ((Q^2)+r (-2 M+r))+J m (r^2) ϵa ω) h00[r]-I Q ((Q^2)+r (-2 M+r)) (J m (10 (Q^2)+r (-16 M+(6+l+(l^2)) r)) ϵa-l (1+l) M (r^4) ω) h10[r]+r ((r^3) (-l (1+l) M ((Q^2)+r (-2 M+r))+2 J m ((Q^2)-2 M r) ϵa ω+M (r^4) (ω^2)) u4[r]+((Q^2)+r (-2 M+r)) (l (1+l) M Q (r^3) h00'[r]+2 I J m Q ((Q^2)-2 M r+(r^2)) ϵa h10'[r]+M (r^2) (-2 ((Q^2)-M r) u4'[r]+r ((Q^2)-2 M r+(r^2)) u4''[r])));
EQML=(1/(l (1+l)))(I l (1+l) Q (r^2) (J m (2 (Q^2)+r (-4 M+r)) ϵa+M (r^4) ω) H00[r]+((Q^2)+r (-2 M+r)) (2 l (1+l) M Q r ((Q^2)+r (-2 M+r)) H10[r]+I l (1+l) Q ((Q^2)+r (-2 M+r)) (-J m ϵa+M (r^2) ω) H20[r]+2 (I l (1+l) M Q (r^4) ω K0[r]+I l (1+l) (J m ((Q^2)-2 M r) ϵa+M (r^4) ω) u1[r]+r (l (1+l) M ((Q^2)+r (-2 M+r)) u2[r]+(-l (1+l) M ((Q^2)-2 M r+(r^2))+J m ((Q^2)-2 M r) ϵa ω) u3[r]+l (1+l) M r ((Q^2)-2 M r+(r^2)) (Q H10'[r]+u2'[r])))));

Axial sector

(*<<"pert_eqs.mx";*)
F[r_]:=1-2M/r+(Q^2)/(r^2);

Define the RW function:

h10[r_]:=r ((ψ[r])/(F[r]));
(*Clear[h10];*)
EQ7//Simplify;
EQ8//Simplify
EQ10//Simplify;
EQM4//Simplify;

ruleJ={};
h0pn=Solve[EQ8==0,h00'[r]][[1]]
ruleh0=Union[h0pn,D[h0pn,r],D[h0pn,{r,2}]];
Simplify[EQ7//.ruleh0];
h0n=Solve[%==0,h00[r]][[1]];
ruleh0b=Union[h0n,D[h0n,r],D[h0n,{r,2}]];
ERW1=EQ8/.ruleh0b/.ruleJ;
ψpp=Solve[ERW1==0,ψ''[r]][[1]];
ERW2=EQM4/.ruleh0b/.ψpp/.ruleJ;

ERW1b=Normal[Simplify[Series[ERW1,{J,0,1}]]]
ERW2b=Normal[Simplify[Series[ERW2,{J,0,1}]]]

Eliminating the third derivatives at first order:

u4pp=Solve[ERW2b==0,u4''[r]][[1]];
ruleu4=Union[u4pp,D[u4pp,r]];
ERW1c=Normal[Simplify[Series[ERW1b//.ruleu4,{J,0,1}]]]

Z1 and Z2

$Assumptions={M>0,Q>0,M>Q,l>1,ω>0,a>0,b>0,A>0,B>0};
Clear[Z1,Z2];
u4[r_]:=A Z1[r]+B Z2[r];
ψ[r_]:=α Z1[r]+β Z2[r];
E1=Simplify[ERW1c];
E2=Simplify[ERW2b];
systZ=Simplify[Solve[{E1==0,E2==0},{Z1''[r],Z2''[r]}]][[1]];
(*c1Z1=Simplify[Coefficient[ERW1c,Z1''[r]]]
c1Z2=Simplify[Coefficient[ERW1c,Z2''[r]]]
c2Z1=Simplify[Coefficient[ERW2b,Z1''[r]]]
c2Z2=Simplify[Coefficient[ERW2b,Z2''[r]]]
rulecoeff={};
CLZ2=Simplify[(c2Z1 ERW1-c1Z1 ERW2b)/.J->0]
CLZ1=Simplify[(c2Z2 ERW1-c1Z2 ERW2b)/.J->0]*)
(*Coefficient[CLZ2,{Z1''[r],Z1'[r],Z1[r],Z2''[r],Z2'[r],Z2[r]}]/.rulecoeff//FullSimplify
Coefficient[CLZ1,{Z1''[r],Z1'[r],Z1[r],Z2''[r],Z2'[r],Z2[r]}]/.rulecoeff//FullSimplify*)
Clear[u4,ψ];

qq1=3 M+(Sqrt[9 (M^2)+4 (-2+l+(l^2)) (Q^2)]);
qq2=3M-(Sqrt[9 (M^2)+4 (-2+l+(l^2)) (Q^2)]);
ruleqq={q1->qq1,q2->qq2};(*ruleTM={α->(I (-3 A M+√((A^2) (9 (M^2)+4 (-2+l+(l^2)) (Q^2)))) ω)/(l (1+l) (-2+l+(l^2)) Q),β->-(I B (3 A M+√((A^2) (9 (M^2)+4 (-2+l+(l^2)) (Q^2)))) ω)/(A l (1+l) (-2+l+(l^2)) Q)};*)
ruleTM={α->((4I A  ω Q )/(l (1+l) q1 )),β->((4I B  ω Q )/(l (1+l)  q2))};
c1=Series[Coefficient[Z1''[r]/.systZ,Z2[r]]/.ruleTM/.ruleqq//Simplify,{J,0,0}]
c2=Series[Coefficient[Z2''[r]/.systZ,Z1[r]]/.ruleTM/.ruleqq//Simplify,{J,0,0}]

FullSimplify[-((I B (3 M+(Sqrt[9 (M^2)+4 (-2+l+(l^2)) (Q^2)])) ω)/(l (1+l) (-2+l+(l^2)) Q))-((4I B  ω Q )/(l (1+l)  q2))/.{q1->qq1,q2->qq2}]
FullSimplify[((I A (-3 M+(Sqrt[9 (M^2)+4 (-2+l+(l^2)) (Q^2)])) ω)/(l (1+l) (-2+l+(l^2)) Q))-((4I A  ω Q )/(l (1+l) q1 ))/.{q1->qq1,q2->qq2}]

systZ

(*Z2 is Zerilli, Z1 is EM*)

DC1=DC[Z1[r]]/.Solve[(Z1''[r]/.systZ)==(1/((F[r])^2))(DC[Z1[r]]-F'[r]F[r] Z1'[r]-((ω^2))Z1[r]),DC[Z1[r]]][[1]];
DC2=DC[Z2[r]]/.Solve[(Z2''[r]/.systZ)==(1/((F[r])^2))(DC[Z2[r]]-F'[r]F[r] Z2'[r]-((ω^2))Z2[r]),DC[Z2[r]]][[1]];

Equation for Z1:

V1v=Series[Coefficient[DC1,Z1[r]],{ϵa,0,1}];
V10=FullSimplify[FullSimplify[((SeriesCoefficient[V1v,0]/.ruleTM)/(F[r] ((r^2)/(r^5))))]/.ruleqq]
V11=FullSimplify[((ϵa SeriesCoefficient[V1v,1]/.ruleTM)/(m J))]
V12=Normal[FullSimplify[((Series[Coefficient[DC1,Z1'[r]],{ϵa,0,1}])/(m J))/.ruleTM]]
W110=FullSimplify[((SeriesCoefficient[Series[Coefficient[DC1,Z2[r]],{ϵa,0,1}],0])/1)/.ruleTM/.ruleqq]
W11=FullSimplify[((SeriesCoefficient[Series[Coefficient[DC1,Z2[r]],{ϵa,0,1}],1]ϵa)/(J m))/.ruleTM]
W12=Normal[((FullSimplify[Series[Coefficient[DC1,Z2'[r]],{ϵa,0,1}]/.ruleTM])/(J m))]

Equation for Z2:

V2v=Series[Coefficient[DC2,Z2[r]],{ϵa,0,1}];
V20=FullSimplify[FullSimplify[((SeriesCoefficient[V2v,0]/.ruleTM)/((F[r] (r^2))/(r^5)))]/.ruleqq]
V21=FullSimplify[((ϵa SeriesCoefficient[V2v,1]/.ruleTM)/(J m))]
V22=Normal[FullSimplify[((Series[Coefficient[DC2,Z2'[r]],{ϵa,0,1}])/(J m))/.ruleTM]]
W210=FullSimplify[SeriesCoefficient[Series[Coefficient[DC2,Z1[r]],{ϵa,0,1}],0]/.ruleTM/.ruleqq]
W21=FullSimplify[((ϵa SeriesCoefficient[Series[Coefficient[DC2,Z1[r]],{ϵa,0,1}],1])/(J m))/.ruleTM]
W22=Normal[((FullSimplify[Series[Coefficient[DC2,Z1'[r]],{ϵa,0,1}]/.ruleTM])/(J m))]
(*

W11=Normal[FullSimplify[Series[Coefficient[DC1,Z2[r]],{ϵa,0,1}]]]
W12=Normal[FullSimplify[Series[Coefficient[DC1,Z2'[r]],{ϵa,0,1}]]]*)

FINAL EQUATIONS

Clear[Z1,Z2];
RHS1=((F[r] (r^2))/(r^5))V10 Z1[r]+ J m(V11 Z1[r]+V12 Z1'[r])+1 J m(W11 Z2[r]+W12 Z2'[r]);
eqF1=FullSimplify[DC[Z1[r]]-RHS1];
CHECK1=Simplify[(DC1-RHS1)//.Union[ruleTM,ruleqq]]
(****)
RHS2=((F[r] (r^2))/(r^5))V20 Z2[r]+1 J m(V21 Z2[r]+V22 Z2'[r])+1 J m(W21 Z1[r]+W22 Z1'[r]);
eqF2=FullSimplify[DC[Z2[r]]-RHS2];
CHECK2=Simplify[(DC2-RHS2)//.Union[ruleTM,ruleqq]]

V10m=((F[r])/(r^3))V10;
V20m=((F[r])/(r^3))V20;
V11m=ϵa J m V11;
V21m=ϵa J m V21;
V12m=ϵa J m V12;
V22m=ϵa J m V22;
W11m=ϵa J m W11;
W21m=ϵa J m W21;
W12m=ϵa J m W12;
W22m=ϵa J m W22;
(***)
Z1[r_]:=Z1m[r];
Z2[r_]:=Z2m[r];
eqZ1m=eqF1;
eqZ2m=eqF2;
Clear[Z1,Z2];

DumpSave["axial_final_KN.mx",{eqZ1m,eqZ2m,V10m,V20m,V11m,V21m,V12m,V22m,W11m,W21m,W12m,W22m}];

COEFFICIENTS for LATEX:

CHECK ZEROTH ORDER:

LV10=Simplify[((((F[r])/(r^3))V10-F[r]l(((l+1))/(r^2)))/1)-((F[r])/(r^3))(-q2+4 ((Q^2)/r))/.ruleqq]
LV20=Simplify[((((F[r])/(r^3))V20-F[r]l(((l+1))/(r^2)))/1)-((F[r])/(r^3))(-q1+4 ((Q^2)/r))/.ruleqq]

LV11=Simplify[((m at (M^2) V11)/(m at))](*factoring (OverTilde[a]) m out*)
LV21=Simplify[((m at (M^2) V21)/(m at))](*factoring (OverTilde[a]) m out*)
(*FullSimplify[LV11+LV21];*)
LV12=Simplify[((m at (M^2) V12)/(((F[r])^2)m at))];(*factoring (F^2) and (OverTilde[a]) m out*)
LV22=Simplify[((m at (M^2) V22)/(((F[r])^2)m at))];(*factoring (F^2) and (OverTilde[a]) m out*)
(*FullSimplify[LV12+LV22];*)
LW11=FullSimplify[((m at (M^2) W11)/((F[r](Q^2)m at)))((((q1/q2)))^(-1))(((-4 (-2+l+(l^2)) (Q^2))/(q2^2)))/.{B->1,A->1}](*factoring F, (OverTilde[a]) m and (Q^2) out*)
LW21=FullSimplify[((m at (M^2) W21)/((F[r](Q^2)m at)))((((q2/q1)))^(-1))(((-4 (-2+l+(l^2)) (Q^2))/(q1^2)))/.{B->1,A->1}];(*factoring F, (OverTilde[a]) m and (Q^2) out*)
LW12=FullSimplify[((m at (M^2) W12)/(((F[r])^2)(Q^2)m at))((((q1/q2)))^(-1))(((-4 (-2+l+(l^2)) (Q^2))/(q2^2)))/.{B->1,A->1}];(*factoring (F^2), (OverTilde[a]) m and (Q^2) out*)
LW22=FullSimplify[((m at (M^2) W22)/(((F[r])^2)(Q^2)m at))((((q2/q1)))^(-1))(((-4 (-2+l+(l^2)) (Q^2))/(q1^2)))/.{B->1,A->1}];(*factoring (F^2), (OverTilde[a]) m and (Q^2) out*)

Check symmetries i--> j

Simplify[LW21-(LW11/.{q1->q2,q2->q1})]
Simplify[LW22-(LW12/.{q1->q2,q2->q1})]
FullSimplify[LV11-(LV21/.{q1->q2,q2->q1})]
Simplify[LV21-(LV11/.{q1->q2,q2->q1})]
Simplify[LV22-(LV12/.{q1->q2,q2->q1})]

Asymptotic Behavior

SetDirectory[NotebookDirectory[]];
<<"axial_final_KN.mx";

Clear[ω];
F[r_]:=1-2M/r+(Q^2)/(r^2);
Z1[r_]:=z1 Exp[I ω r];
Z2[r_]:=z2 Exp[I ω r];
EE1=eqF1/.DC[Z1[r]]->((F[r])^2)Z1''[r]+F[r]F'[r]Z1'[r]+(ω^2)Z1[r];
EE2=eqF2/.DC[Z2[r]]->((F[r])^2)Z2''[r]+F[r]F'[r]Z2'[r]+(ω^2)Z2[r];
Series[{Simplify[(EE1/(Exp[I ω r]))],Simplify[(EE2/(Exp[I ω r]))]},{r,∞,0}]
Clear[Z1,Z2]
rp=M+(Sqrt[(M^2)-(Q^2)]);
bbn=-(1/2) (Sqrt[((ω (2 J m (-(Q^2)+2 M (M+(Sqrt[(M-Q) (M+Q)])))-M (8 (M^4)-8 (M^2) (Q^2)+(Q^4)+8 (M^3) (Sqrt[(M-Q) (M+Q)])-4 M (Q^2) (Sqrt[(M-Q) (M+Q)])) ω))/(M (M-Q) (M+Q)))]);
ruleH={bb->bbn};
Z1[r_]:=z1 (((r-rp))^bb);
Z2[r_]:=z2 (((r-rp))^bb);
EE1=eqF1/.DC[Z1[r]]->((F[r])^2)Z1''[r]+F[r]F'[r]Z1'[r]+(ω^2)Z1[r];
EE2=eqF2/.DC[Z2[r]]->((F[r])^2)Z2''[r]+F[r]F'[r]Z2'[r]+(ω^2)Z2[r];
Simplify[Series[{Simplify[(EE1/(((r-rp))^bb))],Simplify[(EE2/(((r-rp))^bb))]},{r,rp,0}]/.ruleH]
Clear[Z1,Z2]

solb=FullSimplify[Solve[(((4 (-1+bb) bb ((M^2)-(Q^2)))/(((M+(Sqrt[(M^2)-(Q^2)])))^4))+((4 bb (Sqrt[(M^2)-(Q^2)]) ((M^2)-(Q^2)+M (Sqrt[(M^2)-(Q^2)])))/(((M+(Sqrt[(M^2)-(Q^2)])))^5))-(1/(M (((M+(Sqrt[(M^2)-(Q^2)])))^11)))2 J m (256 (M^9)-576 (M^7) (Q^2)+432 (M^5) (Q^4)-120 (M^3) (Q^6)+9 M (Q^8)+256 (M^8) (Sqrt[(M^2)-(Q^2)])-448 (M^6) (Q^2) (Sqrt[(M^2)-(Q^2)])+240 (M^4) (Q^4) (Sqrt[(M^2)-(Q^2)])-40 (M^2) (Q^6) (Sqrt[(M^2)-(Q^2)])+(Q^8) (Sqrt[(M^2)-(Q^2)])) ω+(ω^2))==0,bb]]

$Assumptions={q>0,q<1,q>-1,ω>0,M>0};
kH0=FullSimplify[I bb F'[rp]/.ruleH/.Q->q M/.J->at (M^2),{q>0,q<1,q>-1}]
kH=FullSimplify[((- (-1+(q^2)-(Sqrt[1-(q^2)])) (Sqrt[((M ω (2 at m ((q^2)-2 (1+(Sqrt[1-(q^2)])))+M ((q^4)+8 (1+(Sqrt[1-(q^2)]))-4 (q^2) (2+(Sqrt[1-(q^2)]))) ω))/(1-(q^2)))]))/(M (((1+(Sqrt[1-(q^2)])))^3)))]

FullSimplify[Series[kH,{at,0,1}]]

ζn=2 at ((M^2)/r)-at (Q^2)(M/(r^2));
FullSimplify[((at M)/(2M rp))]

Factor[(q^4)+8 (1+(Sqrt[1-(q^2)]))-4 (q^2) (2+(Sqrt[1-(q^2)]))]

Angular velocity

ζn=2 at ((M^2)/r)-at (Q^2)(M/(r^2));
gtϕ=-ζn ((Sin[θ])^2);
gϕϕ=(r^2)((Sin[θ])^2);
gtt=-(1-2(M/r)+((Q^2)/(r^2)));
omega=((-gtϕ)/(gϕϕ));
rp=M+(Sqrt[(M^2)-(Q^2)]);
Simplify[omega]
Ωh=FullSimplify[Series[omega/.r->rp,{at,0,1}]]

Polar sector

(we follow Zerilli, http://inspirehep.net/record/93100)

Clear[F,ft01,ft02,ft12];
F[r_]:=1-2M/r+(Q^2)/(r^2);
SJ[X_]:=Normal[Series[X,{ϵa,0,1}]]
(*SetDirectory[NotebookDirectory[]];(*this load all the equations at first order in the spin*)
<<"pert_eqs.mx";*)

Check Zerilli

There are some typos in his eqs. 22, 23 and 24. In particular Eq.24 appears to be completely wrong in the charged case. Anyway this equation is never really used, but it should be satisfied automatically if the others are.
Note that Zerilli’s definition of the perturbation functions H0 and H2 is slightly different from ours. (he included background metric quantities as factors)

Using Maxwell equations to obtain a single second-order equation for fem:

u3[r_]:=0;u4[r_]:=0;(*Se set these to zero in the polar sector*)
ruleu=Solve[{((-u1[r]+((I (r^3) ω u2[r])/((Q^2)-2 M r+(r^2)))+r u1'[r])/(r^2))==ft01[r],((u1[r])/r)==ft02[r],(( r u2[r])/((Q^2)+r (-2 M+r)))==ft12[r]},{u1'[r],u1[r],u2[r]}][[1]];(*this solves u1, u2 and u1' in terms of the \tilde f quantities defined in Zerilli*)
ruleU=Union[ruleu,{D[u1'[r]->r ft01[r]+ft02[r]-I r ω ft12[r],r],D[u2[r]->((((Q^2)-2 M r+(r^2)) ft12[r])/r),r]}]

We follow Zerilli and first solve for H2:

H2n=SJ[H20[r]/.Solve[EQ9==0,H20[r]][[1]]]

Then we solve two of the Maxwell equations for \tilde f_01 and \tilde f_02:

f01n=SJ[ft01[r]/.Solve[Simplify[EQM2//.ruleU]==0,ft01[r]][[1]]]
f02n=SJ[ft02[r]/.Solve[Simplify[EQM3//.ruleU]==0,ft02[r]][[1]]]
rulef0102={ft01[r]->f01n,ft02[r]->f02n};

ft12[r_]:=((fem[r])/(F[r]));(*this is just a redefinition*)
eqfem0=SJ[Simplify[SJ[(ft02'[r]+I ω ft12[r]-ft01[r])//.Union[{ft01[r]->f01n,ft01'[r]->D[f01n,r],ft02[r]->f02n,ft02'[r]->D[f02n,r],H20[r]->H2n},ruleU]]]];
(*The equation above is a single ODEs for fem(r) and there are source terms proportional to the metric perturbations as well*)

Coefficient[eqfem0,H10[r]]//Simplify

Gravitational sector

Clear[ft01,ft02,ft12];
ft12[r_]:=((fem[r])/(F[r]));(*same redefinition as above*)
ruleu=Solve[{((-u1[r]+((I (r^3) ω u2[r])/((Q^2)-2 M r+(r^2)))+r u1'[r])/(r^2))==ft01[r],((u1[r])/r)==ft02[r],(( r u2[r])/((Q^2)+r (-2 M+r)))==ft12[r]},{u1'[r],u1[r],u2[r]}][[1]](*again, same as above*)

We insert the solution of H2 obtained before into other Einstein’s equations:

ruleH2={H20[r]->H2n,H20'[r]->D[H2n,r],H20''[r]->D[H2n,{r,2}]};
EQ2b=Simplify[SJ[EQ2//.Union[ruleH2,ruleu,rulef0102]]]
EQ5b=Simplify[SJ[EQ5//.Union[ruleH2,ruleu,rulef0102]]]
EQ6b=Simplify[SJ[EQ6//.Union[ruleH2,ruleu,rulef0102]]]

Note that now the equations contains second derivatives of H00 and K00 which are first order in the spin. We can solve the system perturbatively (a similar procedure is explained more in detail below). First, we solve two of the equations above for K0’ and H0’:

syst1a=Simplify[Solve[{EQ2b==0,EQ6b==0},{H00'[r],K0'[r]}][[1]]];
H00pn=SJ[H00'[r]/.syst1a];
K0pn=SJ[K0'[r]/.syst1a];
rulesyst1a={H00'[r]->H00pn,H00''[r]->D[H00pn,r],K0'[r]->K0pn,K0''[r]->D[K0pn,r]};

Then, at first order in J, the following are now H0’, K0’ and H1’ as only functions of H0, K0, H1 and fem, fem’ anf fem’’:

syst1=Solve[{EQ2b==0,(EQ5b//.rulesyst1a)==0,EQ6b==0},{H00'[r],K0'[r],H10'[r]}][[1]];
H00pn=Simplify[SJ[H00'[r]/.syst1]];
K0pn=Simplify[SJ[K0'[r]/.syst1]];
H10pn=Simplify[SJ[H10'[r]/.syst1]];
rulesyst1={H00'[r]->H00pn,K0'[r]->K0pn,H10'[r]->H10pn};

We can eliminate H0:

EQ3b=Simplify[EQ3//.Union[ruleH2,rulesyst1,ruleu,rulef0102]];
H0n=Simplify[SJ[H00[r]/.Solve[EQ3b==0,H00[r]][[1]]]]
ruleH0={H00[r]->H0n,H00'[r]->D[H0n,r],H00''[r]->D[H0n,{r,2}]};

Therefore, we just have three equations for H1, K and fem:

eq33=Simplify[SJ[EQ2//.Union[ruleH2,ruleH0,ruleu,rulef0102]]];
eq34=Simplify[SJ[EQ5//.Union[ruleH2,ruleH0,ruleu,rulef0102]]];
eqfem=Simplify[SJ[eqfem0//.ruleH2//.Union[rulesyst1,D[rulesyst1,r]]/.ruleH0/.ruleu]];

DumpSave["3eqs_KN_intermediate.mx",{eq33,eq34,eqfem}];(*we backup the output*)

Coupled system of eqs

SetDirectory[NotebookDirectory[]];
<<"3eqs_KN_intermediate.mx";
F[r_]:=1-2M/r+(Q^2)/(r^2);
SJ[X_]:=Normal[Series[X,{ϵa,0,1}]]

Equations 33 and 34 contains second derivatives of H1 and K and third derivatives of fem which are proportional to J. Within our perturbative scheme, we can eliminate them by taking the zeroth order of eq33, eq34 and eqfem, solve them for H10’, K’ and fem’’, then substitute the derivatives into the full form of eq33, eq34 and eqfem and finally expand at first order.

eq330=Simplify[eq33/.ϵa->0];
eq340=Simplify[eq34/.ϵa->0];
eqfem0=Simplify[eqfem/.ϵa->0];
systH1Kfem0=Simplify[Solve[{eq330==0,eq340==0,eqfem0==0},{K0'[r],H10'[r],fem''[r]}]][[1]];
eq33b=Simplify[SJ[eq33/.D[systH1Kfem0,r]]];
eq34b=Simplify[SJ[eq34/.D[systH1Kfem0,r]]];
eqfemb=Simplify[SJ[eqfem/.D[systH1Kfem0,r]]];

DumpSave["3eqs_KN_intermediate2.mx",{eq33b,eq34b,eqfemb}];

Finally, we get a system of three equations that we can solve for H10’, K’ and fem’’:

<<"3eqs_KN_intermediate2.mx";
syst3=Solve[{eq33b==0,eq34b==0,eqfemb==0},{H10'[r],K0'[r],fem''[r]}];
H10pn=Simplify[SJ[H10'[r]/.syst3[[1]]]];
K0pn=Simplify[SJ[K0'[r]/.syst3[[1]]]];
femppn=Simplify[SJ[fem''[r]/.syst3[[1]]]];

Single 2nd order:

We perform a linear transformation on K0 and H10 in order to obtain a single second order equation. The procedure is detailed in Zerilli.

eq33F=K0'[r]-K0pn;
eq34F=H10'[r]-H10pn;

K0[r_]:=Ψ1[r];
H10[r_]:=ω Ψ2[r];(*this is just a redefinition*)
E1=Simplify[eq33F];
E2=Simplify[eq34F];
systΨ=Simplify[Solve[{E1==0,E2==0},{Ψ1'[r],Ψ2'[r]}]][[1]];
Clear[K0,H10]

(*The matrix A is defined in Zerilli*)
A11=Simplify[Coefficient[(Ψ1)'[r]/.systΨ,Ψ1[r]]];
A12=Simplify[Coefficient[(Ψ1)'[r]/.systΨ,Ψ2[r]]];
A21=Simplify[Coefficient[(Ψ2)'[r]/.systΨ,Ψ1[r]]];
A22=Simplify[Coefficient[(Ψ2)'[r]/.systΨ,Ψ2[r]]];
(*Source terms:*)
S1=Simplify[((Ψ1)'[r]/.systΨ)-A11 Ψ1[r]-A12 Ψ2[r]]
S2=Simplify[((Ψ2)'[r]/.systΨ)-A21 Ψ1[r]-A22 Ψ2[r]]

Here we follow Zerilli:

Clear[V];
Clear[F11,F22,F12,F21,n];
AA={{A11,A12},{A21,A22}};
FF={{F11[r]+0ϵa F11b[r] ,F12[r]+0ϵa F12b[r]},{F21[r]+0ϵa  F21b[r],F22[r]+0ϵa  F22b[r]}};
MM={{ϵa  G1[r],1+ϵa W[r]},{-(ω^2)+V+ϵa V1[r],ϵa G2[r]}};
eqMATRIX=Simplify[n[r] Inverse[FF].(AA.FF-D[FF,r])-MM];
(*ZEROth ORDER*)
F[r_]:=1-2(M/r)+((Q^2)/(r^2));
n[r_]:=F[r];
F12[r_]:=1;
F21[r_]:=-I (1+((Q^2)-M r)/((Q^2)-2 M r+(r^2))+(8 (Q^2)-6 M r)/(-4 (Q^2)+r (6 M+(-2+l+(l^2)) r)))
V=1/((r^6) (((-4 (Q^2)+r (6 M+(-2+l+(l^2)) r)))^2))((Q^2)+r (-2 M+r)) (-32 (Q^6)+24 (Q^4) r (6 M+(-2+l+(l^2)) r)+8 (Q^2) (r^2) (-27 (M^2)-2 (-11+4 l (1+l)) M r+3 (-2+l+(l^2)) (r^2))+(r^3) (72 (M^3)+36 (-2+l+(l^2)) (M^2) r+6 (((-2+l+(l^2)))^2) M (r^2)+l (1+l) (((-2+l+(l^2)))^2) (r^3)));
F11[r_]:=(((Q^2)+r (-2 M+r)) (((-4 (Q^2)+r (6 M+(-2+l+(l^2)) r)))^2) F12[r]+2 (r^3) (-6 (M^2) r+M (2 (Q^2)-3 (-2+l+(l^2)) (r^2))+r (-2 (r^2)+l (2 (Q^2)+(r^2))+(l^2) (2 (Q^2)+(r^2)))) n[r])/(2 (r^5) (-4 (Q^2)+r (6 M+(-2+l+(l^2)) r)) n[r]);
F22[r_]:=-((I (r^3) F12[r])/((Q^2)-2 M r+(r^2)));
eqMs=Simplify[Series[eqMATRIX,{ϵa,0,0}]]

The zeroth order is solved. Let us look at the first order:

eqMs1=Simplify[Normal[Series[eqMATRIX,{ϵa,0,1}]]];

DumpSave["eqMs1.mx",eqMs1];

SetDirectory[NotebookDirectory[]];(*backup*)
<<"eqMs1.mx";

We solve for the extra coefficients defined above:

syst1=Simplify[Solve[{eqMs1[[1]][[1]]==0,eqMs1[[1]][[2]]==0,eqMs1[[2]][[1]]==0,eqMs1[[2]][[2]]==0},{V1[r],W[r],G1[r],G2[r]}]][[1]];

Simplify[eqMs1/.syst1]

Single second-order equation: R’’+Ulm R’+Vlm R=Slm [X’=dX/drs, tortoise]

Setting up the single 2nd order eq:

Clear[F];
(***two equations:*)
Ln=Solve[Kt'[r] F[r]==ϵa G1[r] Kt[r]+(1+ϵa W[r]) Lt[r]+S1t[r],Lt[r]][[1]]//Simplify;
eqKt=Simplify[Lt'[r] F[r]-(-((ω^2)-VV[r]-ϵa V1[r])Kt[r]+ϵa G2[r]Lt[r]+S2t[r])//.Union[Ln,D[Ln,r]]];
(***)
c2=Coefficient[eqKt,Kt''[r]];
c1=Normal[Simplify[Series[(((F[r])^2)/(F'[r]F[r]))((Coefficient[eqKt,Kt'[r]])/c2),{ϵa,0,1}]]]
uLM=Simplify[c1-1]
vLM=Normal[Simplify[Series[(((F[r])^2)/1)((Coefficient[eqKt,Kt[r]])/c2),{ϵa,0,1}]]]
sLM=-Normal[Simplify[Series[((((F[r])^2)eqKt)/c2)-(((F[r])^2)Kt''[r]+c1 F'[r]F[r]Kt'[r]+vLM Kt[r]),{ϵa,0,1}]]]

Source terms:

(***)
St=Simplify[SJ[n[r] Inverse[FF].{S1,S2}/.syst1/.fem''[r]->femppn]];
St1=St[[1]];
St2=St[[2]];
ruleS={S1t[r]->St1,S2t[r]->St2};

Explicit form:

ruleSingle=Union[syst1,{VV[r]->V,fem''[r]->femppn},D[syst1,r],ruleS,D[ruleS,r]];
F[r_]:=1-2(M/r)+((Q^2)/(r^2));
ULM=Simplify[uLM//.ruleSingle];
VLM=Simplify[vLM//.ruleSingle];

SLM= Simplify[sLM//.ruleSingle];

St1

EQR=SJ[((F[r])^2)R''[r]+F'[r]F[r](1+ULM)R'[r]+VLM R[r]-SLM];

K0t=(F11[r] R[r]+((F12[r])/(1+ϵa W[r]))(R'[r]F[r]-ϵa G1[r] R[r]-St1))/.syst1;
H10t=ω(F21[r] R[r]+((F22[r])/(1+ϵa W[r]))(R'[r]F[r]-ϵa G1[r] R[r]-St1))/.syst1;

EQRb=Simplify[SJ[EQR/.K0[r]->K0t]];
c2F=Coefficient[EQRb,R''[r]];
ULMF=FullSimplify[SJ[(((F[r])^2)/(F'[r]F[r]))((Coefficient[EQRb,R'[r]])/c2F)]];
VLMF=Simplify[SJ[((F[r])^2)((Coefficient[EQRb,R[r]])/c2F)]];
SLMF=Simplify[SJ[(((F[r])^2)R''[r]+F'[r]F[r]ULMF R'[r]+VLMF R[r])-(((F[r])^2)/c2F)EQRb]];

EQfem=Simplify[SJ[(fem''[r]-femppn)/.{K0[r]->K0t,H10[r]->H10t}]];

DumpSave["final_2nd_order_R.mx",{EQRb,SLMF,ULMF,VLMF,EQfem}];

Simplify[Coefficient[EQRb,R''[r]]];
Simplify[Coefficient[EQRb,R'[r]]];
Simplify[Coefficient[EQRb,R[r]]];
Simplify[Coefficient[EQRb,fem''[r]]];
Simplify[Coefficient[EQRb,fem'[r]]];
Simplify[Coefficient[EQRb,fem[r]]];
(*****)
Simplify[Coefficient[EQfem,R''[r]]]
Simplify[Coefficient[EQfem,R'[r]]]
Simplify[Coefficient[EQfem,R[r]]]
Simplify[Coefficient[EQfem,fem''[r]]]
Simplify[Coefficient[EQfem,fem'[r]]]
Simplify[Coefficient[EQfem,fem[r]]]

Final Equations:

SetDirectory[NotebookDirectory[]];
<<"final_2nd_order_R.mx";
F[r_]:=1-2M/r+(Q^2)/(r^2);
SJ[X_]:=Normal[Series[X,{ϵa,0,1}]]

The equations that we get do not seem to agree with Zerilli’s Eqs.(48) and (49). However, we show below that we can perform a linear transformation that brings our equations to those reported by Chandra in his book, which are correct.

<<"RN_intermediate_Rfem.mx";
eqR0=eqR;
eqf0=eqf;

(*Linear transformation:*)
fem[r_]:=B11[r] Z1[r]+B12[r]Z2[r];
R[r_]:=B21[r]Z1[r]+B22[r]Z2[r];

eqR=Simplify[EQRb];
eqf=Simplify[EQfem];
Clear[R,fem]

Check zeroth order: OK!

eee1=FullSimplify[eqf0-Coefficient[eqf0,Z2''[r]]/Coefficient[eqf,Z2''[r]] eqf/.ϵa->0]
eee2=FullSimplify[eqR0-Coefficient[eqR0,Z2''[r]]/Coefficient[eqR,Z2''[r]] eqR/.ϵa->0]

Continue

systZ=Solve[{eqR==0,eqf==0},{Z1''[r],Z2''[r]}];

Z1ppn=Z1''[r]/.systZ[[1]];
Z2ppn=Z2''[r]/.systZ[[1]];
(***)
CZ1Z1=Simplify[SJ[Coefficient[Z1ppn,Z1[r]]]];
CZ1Z1p=Simplify[SJ[Coefficient[Z1ppn,Z1'[r]]]];
CZ1Z2=Simplify[SJ[Coefficient[Z1ppn,Z2[r]]]];
CZ1Z2p=Simplify[SJ[Coefficient[Z1ppn,Z2'[r]]]];
(***)
CZ2Z2=Simplify[SJ[Coefficient[Z2ppn,Z2[r]]]];
CZ2Z2p=Simplify[SJ[Coefficient[Z2ppn,Z2'[r]]]];
CZ2Z1=Simplify[SJ[Coefficient[Z2ppn,Z1[r]]]];
CZ2Z1p=Simplify[SJ[Coefficient[Z2ppn,Z1'[r]]]];

systBa=DSolve[{(CZ1Z2p/.ϵa->0)==0,(CZ2Z1p/.ϵa->0)==0},{B22[r],B21[r]},r][[1]]/.{}//Simplify

(*This is just to check that C1 and C2 can be written in a more convenient and compact form:*)
q1n=3M+(Sqrt[9 (M^2)+4 (-2+l+(l^2)) (Q^2)]);
q2n=3M-(Sqrt[9 (M^2)+4 (-2+l+(l^2)) (Q^2)]);
Simplify[(((-3 I M ω+(Sqrt[-(9 (M^2)+4 (-2+l+(l^2)) (Q^2)) (ω^2)]))/(8 Q))-(-I ((q2n ω)/(8 Q))))/.{M->1,l->2,Q->(2/10)},ω>0]
Simplify[(((-3 I M ω-(Sqrt[-(9 (M^2)+4 (-2+l+(l^2)) (Q^2)) (ω^2)]))/(8 Q))-(-I ((q1n ω)/(8 Q))))/.{M->1,l->2,Q->(8/10)},ω>0]

Clear[q1,q2];
systC1=Simplify[{C[1]->- I ((q2 ω)/(8 Q)),C[2]->- I ((q1 ω)/(8 Q))}];
(*systC2=Simplify[{C[1]->- I ((q1 ω)/(8 Q)),C[2]->- I ((q2 ω)/(8 Q))}];
systCij=Simplify[{C[1]->- I ((qi ω)/(8 Q)),C[2]->- I ((qj ω)/(8 Q))}];*)
systC=systC1
systq={q1->q1n,q2->q2n};

C1pbis=Simplify[SJ[CZ1Z1p//.Union[systBa,D[systBa,r],D[systBa,{r,2}]]]/.ϵa->0]
C2pbis=Simplify[SJ[CZ2Z2p//.Union[systBa,D[systBa,r],D[systBa,{r,2}]]]/.ϵa->0]

(*eqRbis=Simplify[Z1''[r]-(CZ1Z1p Z1'[r]+CZ1Z1 Z1[r])//.Union[systBa,D[systBa,r],D[systBa,{r,2}]]];
eqfbis=Simplify[Z2''[r]-(CZ2Z2p Z2'[r]+CZ2Z2 Z2[r])//.Union[systBa,D[systBa,r],D[systBa,{r,2}]]];
C1pbis=Coefficient[eqRbis,Z1'[r]]
C2pbis=Simplify[Coefficient[eqfbis,Z2'[r]]]*)

systBb=Simplify[FullSimplify[DSolve[{(-((F[r])^2)C1pbis==F'[r]F[r]),(-((F[r])^2)C2pbis==F'[r]F[r])}/.systC,{B11[r],B12[r]},r]][[1]]/.{C[1]->α,C[2]->β}]
(*systB12=Simplify[FullSimplify[DSolve[(-((F[r])^2)C2pbis==F'[r]F[r])/.systC,B12[r],r]]][[1]]*)

systB=Union[systBb,systBa(*,D[systBa,r],D[systBa,{r,2}]*)]
systBT=Union[systB,D[systB,r],D[systB,{r,2}]];

Simplify[{B21[r],B22[r]}//.systBT/.systC]

Z1Fppn=Z1ppn//.systBT;
Z2Fppn=Z2ppn//.systBT;

systQ=systq;systQ={};
cF11p=Simplify[SJ[((-((F[r])^2))/(F'[r]F[r]))Coefficient[Z1Fppn,Z1'[r]]]/.systC/.systQ];
cF22p=Simplify[SJ[((-((F[r])^2))/(F'[r]F[r]))Coefficient[Z2Fppn,Z2'[r]]]/.systC/.systQ];
cF11=Simplify[SJ[((-((F[r])^2))/1)Coefficient[Z1Fppn,Z1[r]]]/.systC/.systQ];
cF22=Simplify[SJ[((-((F[r])^2))/1)Coefficient[Z2Fppn,Z2[r]]]/.systC/.systQ];
cF12p=Simplify[SJ[((-((F[r])^2))/1)Coefficient[Z1Fppn,Z2'[r]]]/.systC/.systQ];
cF21p=Simplify[SJ[((-((F[r])^2))/1)Coefficient[Z2Fppn,Z1'[r]]]/.systC/.systQ];
cF12=Simplify[SJ[((-((F[r])^2))/1)Coefficient[Z1Fppn,Z2[r]]]/.systC/.systQ];
cF21=Simplify[SJ[((-((F[r])^2))/1)Coefficient[Z2Fppn,Z1[r]]]/.systC/.systQ];

Simplify[SJ[(((F[r])^2)Z1Fppn+(F'[r]F[r] cF11p Z1'[r]+cF11 Z1[r]+cF12p Z2'[r]+cF12 Z2[r]))/.systC/.systQ]]
Simplify[SJ[(((F[r])^2)Z2Fppn+(F'[r]F[r] cF22p Z2'[r]+cF22 Z2[r]+cF21p Z1'[r]+cF21 Z1[r]))/.systC/.systQ]]

eqFZ1=Simplify[SJ[(((F[r])^2)Z1''[r]+(F'[r]F[r] cF11p Z1'[r]+cF11 Z1[r]+cF12p Z2'[r]+cF12 Z2[r]))/.systC/.systQ]];
eqFZ2=Simplify[SJ[(((F[r])^2)Z2''[r]+(F'[r]F[r] cF22p Z2'[r]+cF22 Z2[r]+cF21p Z1'[r]+cF21 Z1[r]))/.systC/.systQ]];

DumpSave["final_polar_Chandra_form_Z1Z2.mx",{cF11p,cF11,cF12p,cF12,cF22p,cF22,cF21p,cF21,eqFZ1,eqFZ2}];

Same form as in the notes:

Clear[Z1,Z2];
<<"final_polar_Chandra_form_Z1Z2.mx";

Eq for Z1: {\cal D} Z1=V10 Z1+ϵa V11 Z1+ϵa V12 Z1’+ϵa W11 Z2+ ϵa W12 Z2’

C11pp=Coefficient[eqFZ1,Z1''[r]];
C11p=Simplify[(((F[r])^2)/(C11pp F'[r]F[r]))Coefficient[eqFZ1,Z1'[r]]/.ϵa->0]
V10=Simplify[-((((F[r])^2)/C11pp)(Coefficient[eqFZ1,Z1[r]]/.ϵa->0)-((ω^2)-F[r]((l(l+1))/(r^2))))];
V11=Simplify[-(((F[r])^2)/C11pp)SeriesCoefficient[Series[Coefficient[eqFZ1,Z1[r]],{ϵa,0,1}],1]ϵa];
V12=Simplify[-(((F[r])^2)/C11pp)SeriesCoefficient[Series[Coefficient[eqFZ1,Z1'[r]],{ϵa,0,1}],1]ϵa];
W11=Simplify[-(((F[r])^2)/C11pp)SeriesCoefficient[Series[Coefficient[eqFZ1,Z2[r]],{ϵa,0,1}],1]ϵa];
W12=Simplify[-(((F[r])^2)/C11pp)SeriesCoefficient[Series[Coefficient[eqFZ1,Z2'[r]],{ϵa,0,1}],1]ϵa];

Eq for Z2: {\cal D} Z2=V20 Z2+ϵa V21 Z2+ϵa V22 Z2’+ϵa W21 Z1+ ϵa W22 Z1’

C22pp=Coefficient[eqFZ2,Z2''[r]];
C22p=Simplify[(((F[r])^2)/(C22pp F'[r]F[r]))Coefficient[eqFZ2,Z2'[r]]/.ϵa->0]
V20=Simplify[-((((F[r])^2)/C22pp)(Coefficient[eqFZ2,Z2[r]]/.ϵa->0)-((ω^2)-F[r]((l(l+1))/(r^2))))];
V21=Simplify[-(((F[r])^2)/C22pp)SeriesCoefficient[Series[Coefficient[eqFZ2,Z2[r]],{ϵa,0,1}],1]ϵa];
V22=Simplify[-(((F[r])^2)/C22pp)SeriesCoefficient[Series[Coefficient[eqFZ2,Z2'[r]],{ϵa,0,1}],1]ϵa];
W21=Simplify[-(((F[r])^2)/C22pp)SeriesCoefficient[Series[Coefficient[eqFZ2,Z1[r]],{ϵa,0,1}],1]ϵa];
W22=Simplify[-(((F[r])^2)/C22pp)SeriesCoefficient[Series[Coefficient[eqFZ2,Z1'[r]],{ϵa,0,1}],1]ϵa];

Z2[r]

The zeroth order potentials do not depend on ω:

D[{V10,V20},ω]//Simplify

The potentials W’s vanish either when Q=0 or when J=0:

Simplify[W21/.systq/.Q->0,M>0]
Simplify[W22/.systq/.Q->0,M>0]

Simplify[W21/.systq/.J->0,M>0]
Simplify[W22/.systq/.J->0,M>0]

Symmetries i-->j OK!

Simplify[W22-(W12/.{q1->q2,q2->q1})/.α->β]
Simplify[W11-(W21/.{q1->q2,q2->q1})/.α->β]
Simplify[V21-(V11/.{q1->q2,q2->q1})/.α->β]
Simplify[V20-(V10/.{q1->q2,q2->q1})/.α->β]
Simplify[V22-(V12/.{q1->q2,q2->q1})/.α->β]

We invert 1-->2 to agree with Chandra (this is equivalent to take the other solution of C[1] and C[2] in the system above)

Clear[Z1,Z2];
V20p=Simplify[V10+F[r]((l(l+1))/(r^2))];
V10p=Simplify[V20+F[r]((l(l+1))/(r^2))];
V21p=V11;
V11p=V21;
V22p=V12;
V12p=V22;
W21p=W11;
W11p=W21;
W22p=W12;
W12p=W22;
Z2[r_]:=Z1p[r];
Z1[r_]:=Z2p[r];
eqZ2p=eqFZ1;
eqZ1p=eqFZ2;
Clear[Z2,Z1];

DumpSave["polar_final_KN.mx",{eqZ1p,eqZ2p,V10p,V20p,V11p,V21p,V12p,V22p,W11p,W21p,W12p,W22p}];

CHECK CHANDRA 183-184 at zeroth order: OK!

Δ=(r^2)F[r];
nn=(((l-1)(l+2))/2);
ωp=nn r+3 M-2 (Q^2)/r;
W=(Δ/(r ((ωp)^2)))(2 nn r+3 M)+(1/(ωp))(nn r+M);
U=(2 nn r +3M)W+(ωp-nn r-M)-2 nn (Δ/(ωp));
V2m=FullSimplify[(Δ/(r^5))U-3 M W];
V1m=FullSimplify[(Δ/(r^5))U+3 M W];
V1p=Simplify[(Δ/(r^5))(U+(1/2)(q1-q2)W)];
V2p=Simplify[(Δ/(r^5))(U-(1/2)(q1-q2)W)];

V1=V10+F[r]((l(l+1))/(r^2));
V2=V20+F[r]((l(l+1))/(r^2));
FullSimplify[(V1p-V2)/.systq]
FullSimplify[(V2p-V1)/.systq]
Clear[W,V1,V2,U];

Explicit form of the coefficients appearing in the paper: 
[we divide by (M^2)/(J m ϵa) in order to factor out the factor \tilde{a}]

V20l=FullSimplify[Series[V20/.{q1->qi,q2->qj},{ω,0,10}]];
V21l=FullSimplify[(((M^2)V21)/(J m ϵa))/.{q1->qi,q2->qj}];
V22l=FullSimplify[((((M^2)V22)/(J m ϵa))/.{q1->qi,q2->qj})/.l (1+l)->Λ//.l (1+l)->Λ//.l (1+l) ->Λ];
W22l=FullSimplify[((((M^2)W22)/(J m ϵa (Q^2)))/.{q1->qi,q2->qj})/.l (1+l)->Λ//.l (1+l)->Λ//.l (1+l) ->Λ];
W21l=FullSimplify[((((M^2)W21)/(J m ϵa (Q^2)))/.{q1->qi,q2->qj})/.l (1+l)->Λ//.l (1+l)->Λ//.l (1+l) ->Λ];

Sometimes we divide by powers of F, in order to factor out the near-horizon behavior. We also set α=1 without loss of generality [remove comma to show the formulae]

((Subscript[V, 0])^((i,+)))/F:

Simplify[(V20l/(F[r]))/.α->1];

((Subscript[V, 1])^((i,+))): [remove comma to show the formula]

Simplify[V21l/.α->1];

((Subscript[V, 2])^((i,+)))/(F^2): [remove comma to show the formula]

FullSimplify[(V22l/((F[r])^2))]

((Subscript[W, 1])^((i,+)))/F: [remove comma to show the formula]

FullSimplify[(W21l/(F[r]))/.{α->1,β->1}]

((Subscript[W, 2])^((i,+)))/(F^2): [remove comma to show the formula]

FullSimplify[(W22l/((F[r])^2))/.{α->1,β->1}]

Near-horizon behavior:

These potentials are vanishing at the horizon:

Simplify[{V20,W21,(W22/(F[r])),(V22/(F[r]))}/.r->M+(Sqrt[(M^2)-(Q^2)])]

While this one gives the superradiant factor:

V21h=FullSimplify[V21/.r->M+(Sqrt[(M^2)-(Q^2)])]

kH=Simplify[Series[(Sqrt[(ω^2)-V21h]),{ϵa,0,1}],ω>0]
FullSimplify[Normal[kH]-(ω- ((m J/(M^2) ϵa)/(M(2+2(Sqrt[1-(Q^2)/(M^2)])-(Q^2)/(M^2))))),0<Q<M]
Simplify[Series[kH,{Q,0,2}],M>0]