Amplification factors for charged bosonic waves in a BH Background
Reference: Richard Brito, Vitor Cardoso and Paolo Pani, “Black-hole superradiance: from fundamental physics to astrophysics” (2015)

Quit[];

SetDirectory[NotebookDirectory[]];

Equation

Here we write the Schroedinger-like equation. In principle any kind of perturbation in different theories can be set in this form, just changing the potential and the metric coefficients. In this example we specialize to D=4 Reissner-Nordstrom BHs and a massless charged scalar.

Parameters

$Assumptions={M>0,a> 0,r>0,a<M,ω>0,Q<M,q>0,Q>0,r>M+(Sqrt[(M^2)-(Q^2)])};

f[r_]:=1-2 M/r+Q^2/r^2;
B[r_]:=f[r];
om[r_]:=0;
U=((l (l+1))/(r^2))+((f'[r])/r);
rp=M+(Sqrt[(M^2)-(Q^2)]);
rm=M-(Sqrt[(M^2)-(Q^2)]);
V=((ω^2)-((2 m om[r] ω)/(r^2))-f[r] U)- ((2q Q ω)/r)+(((q^2)(Q^2))/(r^2));
(***)

To avoid complications, we prescribe kH=(Sqrt[V(rp)]), kinf=(Sqrt[V(infinity)]) and hte tortoise coordinate:

kH=ω-((q Q)/(M+(Sqrt[(M^2)-(Q^2)])));
kinf=ω;
rs[r_]:=r+M Log[(r-rp)(r-rm)]+((2 (M^2)-(Q^2))/(2(Sqrt[(M^2)-(Q^2)]))) Log[((r-rp)/(r-rm))];

Equations and Series    [no need to run this if the external file has been already generated]

Wave equation

Clear[rs];
eqn=f[r] B[r] ψ''[r]+(1/2) (B'[r] f[r]+f'[r] B[r]) ψ'[r]+V ψ[r];

Transformation:

ψ[r_]:=Exp[I k rs[r]]h[r];
eqS=Simplify[(eqn/(E^(I k rs[r])))/.rs'[r]->(1/(Sqrt[f[r]B[r]]))/.D[rs'[r]->(1/(Sqrt[f[r]B[r]])),r]];
Clear[ψ];

Series Horizon (this can take some time)

Here we do the standard series expansion at the horizon. (The order of the expansion, ORDH, can be set freely)

ORDH=3;
eqSH=eqS/.k->-kH;
h[r_]:=Sum[hh[i] (((r-rp))^i),{i,0,ORDH}]; 
ruleH={};
ss=Series[eqSH//.ruleH,{r,rp,ORDH}];
eqsH=Table[SeriesCoefficient[ss,i]==0,{i,1,ORDH}];
yh=Table[hh[i],{i,1,ORDH}];
Length[eqsH]
Length[yh]
seriesH=Solve[eqsH,yh][[1]]//Simplify;
seriesH=Union[seriesH,ruleH];
Clear[h]

Series infinity (this can take some time)

Here we do the standard series expansion at infinity. Again the order is arbitrary, but a high number of termmight be needed  to achieve numerical stability.

ORDINF=3;
eqSinf=eqS;
h[r_]:=Sum[hhinf[i](r^(-i)),{i,0,ORDINF+10}];
ruleINF={(k^2)->(kinf^2),(k^4)->(kinf^4),(k^6)->(kinf^6)};
ss1=Series[{eqSinf}//.ruleINF,{r,∞,ORDINF}];

eqs1=Table[SeriesCoefficient[ss1[[1]],i]==0,{i,2,ORDINF}];
yinf=Table[hhinf[i],{i,1,ORDINF-1}];
Length[eqs1]
Length[yinf]
systINF1=Solve[eqs1,yinf][[1]]//Simplify;
Clear[h,R]

This is needed to extract the coefficients of the asymptotic expansion at infinity with enough precision.

rs[r_]:=r+M Log[(r-rp)(r-rm)]+((2 (M^2)-(Q^2))/(2(Sqrt[(M^2)-(Q^2)]))) Log[((r-rp)/(r-rm))];Qinf=((Exp[I k rs[r]]Sum[((hhinf[i])/(r^i)),{i,0,ORDINF-1}])//.Union[systINF1,{hhinf[0]->Ain,k-> -kinf}])+ ((Exp[I k rs[r]]Sum[((hhinf[i])/(r^i)),{i,0,ORDINF-1}])//.Union[systINF1,{hhinf[0]->Aout,k->kinf}]);
Qpinf=D[Qinf,r];

solBC=Solve[{(Qinf/.r->INF1)==Y1,(Qpinf/.r->INF1)==Y1p},{Ain,Aout}][[1]];

We store the equations to save time in future :

SetDirectory[NotebookDirectory[]];
CreateDirectory["auxiliary_files"];
Put[{systINF1,ORDINF,seriesH,ORDH,eqn,solBC},"auxiliary_files/series_SR_charged.m"]
(*DumpSave["series_super.mx",{systINF1,ORDINF,seriesH,ORDH,eqn,solBC}];*)

Integrator

SetDirectory[NotebookDirectory[]];
ext=Get["auxiliary_files/series_SR_charged.m"];
systINF1=ext[[1]];
ORDINF=ext[[2]];
seriesH=ext[[3]];
ORDH=ext[[4]];
eqn=ext[[5]];
solBC=ext[[6]];

INT[ω0_,l0_,m0_,Q0_,q0_,a0_,s0_,EPS_,INF0_]:=(
ruleP={AccuracyGoal->80,PrecisionGoal->12,MaxSteps->(10^6)};
param={ω->ω0,l->l0,m-> m0,a-> a0,s-> s0,hh[0]->1,q->q0,Q->Q0,M->1};
(*Integration*)
rh=rp//.param;
rhn=rh(1+EPS);
INF=(INF0/ω)//.param;
ψg=Exp[-I kH rs[r]]Sum[hh[i](((r-rp))^i),{i,0,ORDH}]//.Union[param,seriesH];
BCs={ψ[rhn]==ψg/.r->rhn,ψ'[rhn]==D[ψg,r]/.r->rhn}//.param;
EQs={eqn==0}//.param;
sol=NDSolve[Union[EQs,BCs],ψ[r],{r,rhn,INF},ruleP];
ψn=ψ[r]/.sol[[1]];
ψpn=D[ψn,r];
(****)
Ainn=Ain//.Union[solBC,param,{r-> INF,INF1-> INF,Y1-> ψn,Y1p-> ψpn}];
Aoutn=Aout//.Union[solBC,param,{r-> INF,INF1-> INF,Y1-> ψn,Y1p-> ψpn}];

(*Amplification Factors*)
Z=(((((Abs[Aoutn])/(Abs[Ainn]))))^2)//.param;
Z-1);

INT[(10^(-1)),0,0,0.5,1,0,0,1 (10^(-4)),3 (10^3)]