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

Quit[];

SetDirectory[NotebookDirectory[]];

Teukolsky’s equation

Here we write Teukolsky’s equation with in wave-like form. In principle any kind of perturbation in different theories can be set in this form, just changing the potential.

Black hole parameter

$Assumptions={M>0,a> 0,r>0,a<M,ω>0};
Δ[r_]:=(r^2)+(a^2)-2 M r;
rplus=M+(Sqrt[(M^2)-(a^2)]);
rminus=M-(Sqrt[(M^2)-(a^2)]);
Ωh=(a/((rplus^2)+(a^2)));
f[r_]:=((Δ[r])/((r^2)+(a^2)));

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

Wave equation

V=(((KK^2)-2I s(r-M)KK+Δ[r](4I s ω r-λ))/((((r^2)+(a^2)))^2))-((G[r])^2)-f[r]G'[r];
G[r_]:=((s(r-M))/((r^2)+(a^2)))+((r Δ[r])/((((r^2)+(a^2)))^2));
KK=((r^2)+(a^2))ω-a m;
eqn=((f[r])^2)R''[r]+f'[r]f[r]R'[r] +V R[r]//Simplify;

VV=V/.λ->A[s,l,m,ω]+(a^2)(ω^2)-2a m ω

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)

R[r_]:=(((r-rplus))^AA)(((r-rplus))^(-s/2))(((r-rminus))^(-s/2)) h[r];
eqn1=(eqn/(((r-rplus))^AA))(((r-rplus))^(s/2))(((r-rminus))^(s/2))//Simplify;
Clear[R]

ORDH=3;
h[r_]:=Sum[hh[i] (((r-rplus))^i),{i,0,ORDH}]; 
ruleH={AA->((I (a m-2 M (M+(Sqrt[-(a^2)+(M^2)])) ω))/(2 (Sqrt[-(a^2)+(M^2)])))};
ss=Series[eqn1//.ruleH,{r,rplus,ORDH}]//Chop//Simplify;
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 terms is needed for s=-2 to achieve numerical stability.

R[r_]:=Exp[kk r](r^AA)H[r];
eq1b=(eqn/((E^(kk  r))(r^AA) ))//Simplify;
Clear[R]

ORDINF=4;
H[r_]:=Sum[hhinf[i](r^(-i)),{i,0,ORDINF+10}];
ruleINF={(kk^2)->-(ω^2),AA->((-I s ω-2 M (ω^2))/kk),(kk^4)->(ω^4),(kk^6)->-(ω^6)};
ss1=Series[{eq1b}//.ruleINF,{r,∞,ORDINF}]//Chop;

eqs1=Table[SeriesCoefficient[ss1[[1]],i]==0,{i,0,ORDINF}];
yinf=Table[hhinf[i],{i,1,ORDINF-1}];
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.

Qinf=(((E^(kk  r))(r^((-I s ω-2 M (ω^2))/kk))Sum[((hhinf[i])/(r^i)),{i,0,ORDINF-1}])//.Union[systINF1,{hhinf[0]->BB1,kk-> -I ω}])+ (((E^(kk  r))(r^((-I s ω-2 M (ω^2))/kk))Sum[((hhinf[i])/(r^i)),{i,0,ORDINF-1}])//.Union[systINF1,{hhinf[0]->CC1,kk->I ω}]);
Qpinf=D[Qinf,r];

solBC=Solve[{(Qinf/.r->INF1)==Y1,(Qpinf/.r->INF1)==Y1p},{BB1,CC1}][[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.m"]
(*DumpSave["series_super.mx",{systINF1,ORDINF,seriesH,ORDH,eqn,solBC}];*)

Integrator

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

INT[ω0_,a0_,s0_,l0_,m0_,EPS_,INF0_,NITMAX_]:=(
ruleP={AccuracyGoal->80,PrecisionGoal->12,MaxSteps->(10^6)};

k0n=I ω;

param={ω->ω0,l->l0,m-> m0,a-> a0,s-> s0,hh[0]->1,kk->k0n,k0->k0n};
M=1;

(*calculating the eigenvalues (Subscript[A, lmj])*)
km=1/2*Abs[m-s];kp=1/2*Abs[m+s];Alm=l*(l+1)-s*(s+1)-((2m (s^2))/(l(l+1)))a ω//.param;
(*NITMAX=150;*)
γ=Function[n,2*a*ω*(n+km+kp+s)//.param];
β=Function[n,n*(n-1)+2*n*(km+kp+1-2*a*ω)-(2*a*ω*(2*km+s+1)-(km+kp)*(km+kp+1))-((a^2)*(ω^2)+s*(s+1)+Sep)//.param];
α=Function[n,-2*(n+1)*(n+2*km+1)//.param];
Leaver31[Sep_]:=Module[{Rp},
For[{n=NITMAX;Rp=-1.0;},n>0,{Rp=((γ[n])/(β[n]-α[n]*Rp));n--;}
];Rp];Leaver33[Sep_]:=((β[0])/(α[0]))-Leaver31[Sep];

A=FindRoot[Leaver33[Sep]==0,{Sep,Alm}][[1,2]];
λ=A+(a^2)(ω^2)-2a m ω//.param;


(*Integration*)
rplus=M+(Sqrt[(M^2)-(a^2)]);
rh=rplus//.param;
rhn=rh(1+EPS);
INF=Abs[INF0]//.param;

Rg=(((r-rplus))^((I (a m-2 M (M+(Sqrt[-(a^2)+(M^2)])) ω))/(2 (Sqrt[-(a^2)+(M^2)]))))(((r-rplus))^(-s/2))(((r-rminus))^(-s/2)) Sum[hh[i](((r-rplus))^i),{i,0,ORDH}]//.Union[param,seriesH];
BCs={R[rhn]==Rg/.r->rhn,R'[rhn]==D[Rg,r]/.r->rhn}//.param;
EQs={eqn==0}//.param;
sol=NDSolve[Union[EQs,BCs],R[r],{r,rhn,INF},ruleP];
Rn=R[r]/.sol[[1]];
Rnp=D[Rn,r];

(*Match the result with the expansion at infinity*)
B=(((λ+s(s+1)))^2)+4m a ω-4(a^2)(ω^2); 
Cc=B(((((λ+s(s+1))-2))^2)+36a ω m-36(a^2)(ω^2))+(2(λ+s(s+1))-1)(96(a^2)(ω^2)-48a ω m)+144(ω^2)((M^2)-(a^2));

CCn=CC1//.Union[solBC,param,{r-> INF,INF1-> INF,Y1-> Rn,Y1p-> Rnp,kk-> k0n}];
BBn=BB1//.Union[solBC,param,{r-> INF,INF1-> INF,Y1-> Rn,Y1p-> Rnp,kk->-k0n}];

(*Amplification Factor for different s*)
If[s0== 0,Z=(((((Abs[CCn])/(Abs[BBn]))))^2)//.param ];
If[s0==-1,Z=(B/(16(ω^4)))(((((Abs[CCn])/(Abs[BBn]))))^2)//.param];
If[s0==-2,Z=(Cc/(256(ω^8)))(((((Abs[CCn])/(Abs[BBn]))))^2)//.param] ;
If[s0==1,Z=((((B/(16(ω^4)))))^(-1))(((((Abs[CCn])/(Abs[BBn]))))^2)//.param];
If[s0==2,Z=((((Cc/(256(ω^8)))))^(-1))(((((Abs[CCn])/(Abs[BBn]))))^2)//.param] ;
Z-1);

ww=0.5;
Timing[INT[0.5,0.99,-2,2,2,1 (10^(-3)),(100/ww),200]]