Quasinormal modes of Schwarzschild BHs via direct integration (a la Chandrasekhar-Detweiler)

[Chandrasekhar-Detweiler, Proc. R. Soc. Lond. A344 441-452 (1975); notebook written by Paolo Pani]

We solve the Regge-Wheeler equation (due to isospectrality, this describes the full spectrum)

f[r_]:=1-2M/r;
V=f[r](((l(l+1))/(r^2))-6(M/(r^3)));
EQ=((f[r])^2)ψ''[r]+f'[r]f[r]ψ'[r]+((ω^2)-V)ψ[r]

Series at the horizon

We solve the field equation in the near-horizon limit requiring purely ingoing wave boundary condition.
We use the ansatz Ψ=(((r-2M))^(-2I ω M))(∑_(n=0)^∞)(((Subscript[h, n])(r-2M))^n)
Higher truncation order (>~10) is needed for stability.

ORDH=13;
(*ψ[r_]:=(((r-2M))^(-2I ω M))h[r];
f[r_]:=1-2M/r;
V=f[r](((l(l+1))/(r^2))-6(M/(r^3)));
EQ=((f[r])^2)ψ''[r]+f'[r]f[r]ψ'[r]+((ω^2)-V)ψ[r];
Simplify[EQ]*)
h[r_]:=Sum[hh[i] (((r-2M))^i),{i,0,ORDH+5}];
EQH=(-r (l+(l^2)-(r^2) (ω^2))+2 M (3+I r ω+(r^2) (ω^2))) h[r]+r (2 (M-2 I M r ω) h'[r]+r (-2 M+r) h''[r]);
ss=Simplify[Series[EQH,{r,2M,ORDH}]];
eqsH=Table[SeriesCoefficient[ss,i]==0,{i,0,ORDH}];
yh=Table[hh[i],{i,1,ORDH+1}];
seriesH=Simplify[Solve[eqsH,yh]][[1]];
Length[yh]
Length[eqsH]
Clear[ψ,h]

Series (infinity)

We solve the field equation in the large distance limit requiring no incoming wave boundary condition
We use the ansatz Ψ=(e^(I ω (SubStar[r])))(∑_(n=0)^∞)((Subscript[g, n])/(r^n)), where (SubStar[r]) is the tortoise coordinate.
Higher truncation order (>~10) is needed for stability.

ORDINF=14;
(*ψ[r_]:=Exp[I ω r]g[r];

EQ=((f[r])^2)ψ''[r]+f'[r]f[r]ψ'[r]+((ω^2)-V)ψ[r];
Simplify[EQ]*)
g[r_]:=Sum[((ginf[i])/(r^i)),{i,0,ORDINF+5}];
f[r_]:=1-2M/r;
V=f[r](((l(l+1))/(r^2))-6(M/(r^3)));
EQINF=g''[r] ((f[r])^2)+g'[r] f'[r]f[r]+2 I ω f[r] g'[r]-V g[r];
ss=Simplify[Series[EQINF,{r,∞,ORDINF}]];
eqsINF=Table[SeriesCoefficient[ss,i]==0,{i,2,ORDINF}];
yinf=Table[ginf[i],{i,1,ORDINF-1}];
seriesINF=Simplify[Solve[eqsINF,yinf]][[1]];
Length[yinf]
Length[eqsINF]
Clear[ψ,g]

Integration

We perform two integrations: (I) from the horizon to a matching point r=rm imposing the purely ingoing wave boundary condition, (II) from infinity to the matching point, imposing no incoming waves at infinity. Then, we construct a linear combination of the two solutions and fix the ratio of the coefficients in order for the linear combination to be (C^0) at the matching point. However, for a generic frequency the linear combination will have discontinous derivatives. Finally, we find the eigenvalue by imposing that the linear combination is (C^1)

INT[ω0_,rm0_,l0_,EPS0_,INF0_]:=(
ruleP={AccuracyGoal->12,PrecisionGoal->12};
param={M->1,l->l0,ω->ω0,hh[0]->1,ginf[0]->1};
rh=2M/.param;
EPS=EPS0;
rinf=INF0 M/.param;
rhn=(1+EPS)rh;
(****INTEGRATION FROM THE HORIZON***)
ψgH=(((r-2M))^(-2 I ω M))Sum[hh[i] (((r-2M))^i),{i,0,ORDH+1}]//.Union[param,seriesH];
BCs={ψ[rhn]==ψgH/.r->rhn,ψ'[rhn]==D[ψgH,r]/.r->rhn};
EQs={((-(-((6 M)/(r^3))+((l (1+l))/(r^2))) (1-((2 M)/r))+(ω^2)) ψ[r]+((2 M (1-((2 M)/r)) ψ'[r])/(r^2))+(((1-((2 M)/r)))^2) ψ''[r])==0}/.param;
rm=rm0;
solH=NDSolve[Union[EQs,BCs],ψ,{r,rhn,rm},ruleP];
ψnH=ψ[r]/.solH[[1]];
ψnHp=D[ψnH,r];
(*******)
(***INTEGRATION FROM INFINITY****)
rs=r+2M Log[(r/(2M))-1]/.param;(*tortoise coordinate*)
ψgINF=Exp[I ω rs]Sum[((ginf[i])/(r^i)),{i,0,ORDINF-1}]//.Union[param,seriesINF];
BCs={ψ[rinf]==ψgINF/.r->rinf,ψ'[rinf]==D[ψgINF,r]/.r->rinf};
EQs={((-(-((6 M)/(r^3))+((l (1+l))/(r^2))) (1-((2 M)/r))+(ω^2)) ψ[r]+((2 M (1-((2 M)/r)) ψ'[r])/(r^2))+(((1-((2 M)/r)))^2) ψ''[r])==0}/.param;
rm=rm0;
solINF=NDSolve[Union[EQs,BCs],ψ,{r,rinf,rm},ruleP];
ψnINF=ψ[r]/.solINF[[1]];
ψnINFp=D[ψnINF,r];
(*******)
b=((ψnH)/(ψnINF))/.r->rm;(*we fix the coefficients of the linear combination in order for ψn to be continuous at r=rm*)
ψn=Piecewise[{{ψnH,r<rm},{b ψnINF,r≥rm}}];
Δψp=((D[ψnH,r]-D[b ψnINF,r])/(D[ψnH,r]))/.r->rm;
{Δψp}(*the function INT returns the jump of the first derivative at the matching point for a given ω*)
);

Example: for a generic ω, the function ψ is continuous but its derivative is discontinuous at the matching point

INT[0.1-0.1I,6,2,(10^(-4)),(10^2)]
Plot[Abs[ψn],{r,0.5 rm,2 rm}]

Root finder: we search for ω such that the first derivative is continuous

Note that the value of the numerical infinity is not very large. This must be compensated with a large number of terms in the series expansion above.

Fundamental mode for l=2

ωguess=0.37367 -0.08896I;
ωguess=0.4 -0.1 I;
f4[y_?NumericQ]:=INT[y,6,2,(10^(-4)),3 (10^1)][[1]];
y0n=y/.FindRoot[f4[y]==0,{y,ωguess}][[1]];
check=Abs[f4[y0n]];
{y0n,check}

The solution does not depend on the matching point:

ωguess=0.4 -0.1 I;
f4[y_?NumericQ]:=INT[y,5,2,(10^(-4)),3 (10^1)][[1]];
y0n=y/.FindRoot[f4[y]==0,{y,ωguess}][[1]];
check=Abs[f4[y0n]];
{y0n,check}
(****)
f4[y_?NumericQ]:=INT[y,6,2,(10^(-4)),3 (10^1)][[1]];
y0n=y/.FindRoot[f4[y]==0,{y,ωguess}][[1]];
check=Abs[f4[y0n]];
{y0n,check}
(****)
f4[y_?NumericQ]:=INT[y,7,2,(10^(-4)),3 (10^1)][[1]];
y0n=y/.FindRoot[f4[y]==0,{y,ωguess}][[1]];
check=Abs[f4[y0n]];
{y0n,check}

We also find the first overtone for l=2:

ωguess=0.35 -0.2 I;
f4[y_?NumericQ]:=INT[y,6,2,(10^(-4)),2 (10^1)][[1]];
y0n=y/.FindRoot[f4[y]==0,{y,ωguess}][[1]];
check=Abs[f4[y0n]];
{y0n,check}

References

[1] S. Chandrasekhar and S. Detweiler,"The quasi-normal modes of the Schwarzschild black hole.",Proc.Roy.Soc.Lond.A344 (1975) 441-452.

[2] C. Molina, P. Pani, V.Cardoso and L. Gualtieri,"Gravitational signature of Schwarzschild black holes in dynamical Chern-Simons gravity,"
Phys.Rev.D81 (2010) 124021; e-Print:arXiv:1004.4007[gr-qc].