(*This notebook computes scattering of a monochromatic pulse, or of a spatially gaussian pulse off a Schwarzschild black hole*)

Monochromatic pulse in Schwarzschild

M=1;
f=1-2M/r;f1=D[f,r];

(*define tortoise coordinate*)rt=Integrate[1/f,r];

l=0;
r1=2*(1+10^(-4));

ii=0;

Nmaxx=600;
While[ii<Nmaxx+1,
w=0.002+2*10^(-3)*ii;

(*Numerical infinity is chosen such that the domain always contains a large number of wavelengths, thus ensuring we are in the wave zone*)r2=100/w;

ahor=((I (1+l+(l^2)))/(2 M (I+4 M w)));bhor=-((2 (l^3)+(l^4)+12 I M w+l (2+8 I M w)+(l^2) (3+8 I M w))/(16 (M^2) (-1+6 I M w+8 (M^2) (w^2))));

(*Expanded field equation close to the horizon to get more accurate behavior*)
Hor=(1+ahor*(r-2*M)+bhor*(r-2*M)^2)*Exp[-I*w*rt];

BC=Hor/.r->r1;
BC1=D[Hor,r]/.r->r1;
EQ=f^2*D[y[r],{r,2}]+f*f1*D[y[r],{r,1}]+(w^2-f*(l*(l+1)/r^2+f1/r))*y[r];

sol=NDSolve[{EQ==0,y[r1]==BC,y'[r1]==BC1},y,{r,r1,r2}(*,
            AccuracyGoal -> 18, PrecisionGoal -> 18,
                             WorkingPrecision -> 24*),MaxSteps->340000];

ain=-((I l (1+l))/(2 w));bin=((2 l+(l^2)-2 (l^3)-(l^4)-4 I (M) w)/(8 (w^2)));cin=((I l (12+4 l-15 (l^2)-5 (l^3)+3 (l^4)+(l^5))-12 (-2+3 l+3 (l^2)) (M) w)/(48 (w^3)));din=(1/(384 (w^4)))(-56 (l^5)-14 (l^6)+4 (l^7)+(l^8)+4 (l^3) (49+72 I M w)+(l^4) (49+144 I M w)-144 (M) w (-2 I+3 M w)-48 I l (-3 I+14 M w)-12 I (l^2) (-3 I+44 M w));
aout=((I l (1+l))/(2 w));bout=((2 l+(l^2)-2 (l^3)-(l^4)+4 I (M) w)/(8 (w^2)));cout=((15 I (l^3)+5 I (l^4)-3 I (l^5)-I (l^6)+24 (M) w-12 l (I+3 M w)-4 (l^2) (I+9 M w))/(48 (w^3)));dout=(1/(384 (w^4)))(-56 (l^5)-14 (l^6)+4 (l^7)+(l^8)+4 (l^3) (49-72 I M w)+(l^4) (49-144 I M w)-144 (M) w (2 I+3 M w)+48 I l (3 I+14 M w)+12 I (l^2) (3 I+44 M w));

yinf=y[r2]/.sol[[1]];yinf1=y'[r2]/.sol[[1]];


rtinf=r+2 (M) Log[-2 M+r];

yinfana=Cin*(1+ain/r+bin/r^2+cin/r^3+din/r^4)*Exp[-I*w*rtinf]+Cout*(1+aout/r+bout/r^2+cout/r^3+dout/r^4)*Exp[I*w*rtinf]/.r->r2;
yinfana1=D[Cin*(1+ain/r+bin/r^2)*Exp[-I*w*rtinf]+Cout*(1+aout/r+bout/r^2)*Exp[I*w*rtinf],r]/.r->r2;
solcin=Solve[{yinf==yinfana,yinf1==yinfana1},{Cin,Cout}][[1]];

answer[ii]={w,(1-Abs[Cout]^2/Abs[Cin]^2)/.solcin};ii++]

Tabe=Table[answer[ii],{ii,0,Nmaxx}];

(*Amplification Factor*)Amp=Table[{N[Tabe[[i,1]]],Tabe[[i,2]]},{i,1,Length[Tabe]}];
ListPlot[Amp,PlotRange->All]

Numerics in frequency - domain for Schwarzschild: scattering of a Gaussian

M=1;
f=1-2M/r;f1=D[f,r];

(*define tortoise coordinate*)rt=Integrate[1/f,r];

l=1;

(*Initial data I, which I consider to be only psidot*)Initdata=Exp[-(r+2 Log[-2+r]-r0)^2/sigma^2];

(*r0 is the location, in tortoise coordinates, of the peak of the Gaussian initial data*)r0=10;

(*sigma is the width of the Gaussian wavepacket*)sigma=6;

r1=2*(1+10^(-4));
(*[r1,r2] define the integration interval over the source and avoids integrating over too large and unnecessary a domain:*)
rsource1=r1;
rsource2=r0+4*sigma;

ii=0;

Nmaxx=600;
While[ii<Nmaxx+1,
w=0.002+2*10^(-3)*ii;

(*Numerical infinity is chosen such that the domain always contains a large number of wavelengths, thus ensuring we are in the wave zone*)r2=100/w;

ahor=((I (1+l+(l^2)))/(2 M (I+4 M w)));bhor=-((2 (l^3)+(l^4)+12 I M w+l (2+8 I M w)+(l^2) (3+8 I M w))/(16 (M^2) (-1+6 I M w+8 (M^2) (w^2))));

(*Expanded field equation close to the horizon to get more accurate behavior*)Hor=(1+ahor*(r-2*M)+bhor*(r-2*M)^2)*Exp[-I*w*rt];

BC=Hor/.r->r1;
BC1=D[Hor,r]/.r->r1;
EQ=f^2*D[y[r],{r,2}]+f*f1*D[y[r],{r,1}]+(w^2-f*(l*(l+1)/r^2+f1/r))*y[r];

sol=NDSolve[{EQ==0,y[r1]==BC,y'[r1]==BC1},y,{r,r1,r2}(*,
            AccuracyGoal -> 18, PrecisionGoal -> 18,
                             WorkingPrecision -> 24*),MaxSteps->340000];

ain=-((I l (1+l))/(2 w));bin=((2 l+(l^2)-2 (l^3)-(l^4)-4 I (M) w)/(8 (w^2)));cin=((I l (12+4 l-15 (l^2)-5 (l^3)+3 (l^4)+(l^5))-12 (-2+3 l+3 (l^2)) (M) w)/(48 (w^3)));din=(1/(384 (w^4)))(-56 (l^5)-14 (l^6)+4 (l^7)+(l^8)+4 (l^3) (49+72 I M w)+(l^4) (49+144 I M w)-144 (M) w (-2 I+3 M w)-48 I l (-3 I+14 M w)-12 I (l^2) (-3 I+44 M w));
aout=((I l (1+l))/(2 w));bout=((2 l+(l^2)-2 (l^3)-(l^4)+4 I (M) w)/(8 (w^2)));cout=((15 I (l^3)+5 I (l^4)-3 I (l^5)-I (l^6)+24 (M) w-12 l (I+3 M w)-4 (l^2) (I+9 M w))/(48 (w^3)));dout=(1/(384 (w^4)))(-56 (l^5)-14 (l^6)+4 (l^7)+(l^8)+4 (l^3) (49-72 I M w)+(l^4) (49-144 I M w)-144 (M) w (2 I+3 M w)+48 I l (3 I+14 M w)+12 I (l^2) (3 I+44 M w));

yinf=y[r2]/.sol[[1]];yinf1=y'[r2]/.sol[[1]];


rtinf=r+2 (M) Log[-2 M+r];

yinfana=Cin*(1+ain/r+bin/r^2+cin/r^3+din/r^4)*Exp[-I*w*rtinf]+Cout*(1+aout/r+bout/r^2+cout/r^3+dout/r^4)*Exp[I*w*rtinf]/.r->r2;
yinfana1=D[Cin*(1+ain/r+bin/r^2)*Exp[-I*w*rtinf]+Cout*(1+aout/r+bout/r^2)*Exp[I*w*rtinf],r]/.r->r2;
solcin=Solve[{yinf==yinfana,yinf1==yinfana1},{Cin,Cout}][[1]];

Aint=NIntegrate[Initdata*y[r]/.sol//Release,{r,rsource1,rsource2}(*,WorkingPrecision->8,PrecisionGoal->8,MaxRecursion->6*)][[1]];

answer[ii]={w,Cin/.solcin[[1]],Psi=Aint/(2*I*w*Cin)/.solcin[[1]],P=w^2*Abs[Psi]^2};ii++]

Tabe=Table[answer[ii],{ii,0,Nmaxx}];

(*Complex, fourier-domain waveform*)Wave=Table[{N[Tabe[[i,1]]],Tabe[[i,3]],Tabe[[i,4]]},{i,1,Length[Tabe]}];

(*Flux*)Tabepower1=Table[{Tabe[[i,1]],Tabe[[i,4]]},{i,1,Length[Tabe]}];ListPlot[Tabepower1,PlotRange->All]

(*Bromwich integral to compute inverse*)Wavetime2[t_]:=1/2*Pi*Sum[Re[Wave[[i,2]]*Exp[-I*Wave[[i,1]]*(t)]]*(Wave[[i+1,1]]-Wave[[i,1]]),{i,1,Length[Wave]-1}];

Plot[Wavetime2[t],{t,-100,100},PlotRange->All]

(*selection of the ringdown phase to estimate the QNM frequency*)ringdownpart=Table[{t-30,Wavetime2[t]},{t,30,60}];

ListPlot[ringdownpart]

par=FindFit[ringdownpart,a0*Exp[-wi*t]*Cos[wr*t+a1]+a2,{a0,a1,a2,wr,wi},t]

Leaver method yields (see BCW, https://arxiv.org/abs/0905.2975):w= 0.2929361329967782` -0.09765998891161318` I

Kerr scalars

M=1;
ar=0.9;
rplus=M+Sqrt[M^2-ar^2];
rminus=M-Sqrt[M^2-ar^2];
Omega=ar/(2*M*rplus);
Delta=r^2+ar^2-2*M*r;Delta1=D[Delta,r];kabs=Abs[m];

l=1;
m=1;

r1=rplus*(1+10^(-6));
r2=500/w;

i=0;

While[i<100,
w=0.01+0.01*i;
NITMAX=100;γang=Function[n,2*w*ar*(n+kabs)];
βang=Function[n,n*(n-1)+2*n*(kabs+1-2*ar*w)-(2*ar*w*(kabs+1)-kabs*(kabs+1))-ar^2*w^2-Akjm];
αang=Function[n,-2*(n+1)*(n+Abs[m]+1)];Leaver31Ang[Sep_]:=Module[{Rn},
For[{n=NITMAX;Rn=-1.0;},n>0,{Rn=((γang[n])/(βang[n]-αang[n]*Rn));n--;}];Rn];Leaver33Ang[Sep_]:=((βang[0])/(αang[0]))-Leaver31Ang[Sep];FR=FindRoot[Leaver33Ang[Sep]==0,{Akjm,l*(l+1)}];

Alm=FR[[1]][[2]];
sigma=(rplus^2+ar^2)*(w-m*Omega)/(rplus-rminus);

BC=(r-rplus)^(-I*sigma)/.r->r1;
BC1=D[(r-rplus)^(-I*sigma),r]/.r->r1;

sol=NDSolve[{Delta*D[y[r],{r,2}]+Delta1*D[y[r],{r,1}]+((w*(r^2+ar^2)-m*ar)^2/Delta-Alm+2*m*ar*w-w^2*ar^2)*y[r]==0,y[r1]==BC,y'[r1]==BC1},y,{r,r1,r2},
            (*AccuracyGoal -> 14, PrecisionGoal -> 14,
                             WorkingPrecision -> 20,*)MaxSteps->1000000];
yr2=y[r2]/.sol[[1]];y1r2=y'[r2]/.sol[[1]];
rt=((r (rminus-rplus)+((ar^2)+(rminus^2)) Log[r-rminus]-((ar^2)+(rplus^2)) Log[r-rplus])/(rminus-rplus));
aout=((I (Alm+(ar^2) (w^2)))/(2 w));
bout=-(((-2+Alm) Alm-4 I M w+2 ar (ar+Alm ar-4 I m M) (w^2)+(ar^4) (w^4))/(8 (w^2)));
ain=-((I (Alm+(ar^2) (w^2)))/(2 w));
bin=-(((-2+Alm) Alm+4 I M w+2 ar (ar+Alm ar+4 I m M) (w^2)+(ar^4) (w^4))/(8 (w^2)));
fim=Solve[{yr2==Aout(1+aout/r+bout/r^2)*Exp[I*w*rt]*r^(-1)+Ain(1+ain/r+bin/r^2)*Exp[-I*w*rt]*r^(-1)/.r->r2,y1r2==D[Aout(1+aout/r+bout/r^2)*Exp[I*w*rt]*r^(-1)+Ain(1+ain/r+bin/r^2)*Exp[-I*w*rt]*r^(-1),r]/.r->r2},{Ain,Aout}];Ain0=fim[[1]][[1]][[2]];Aout0=fim[[1]][[2]][[2]];
Absor[i]={w,(1-Abs[Aout0]^2/Abs[Ain0]^2)};
i=i+1];

(*This gives 1-RR^2, which is negative since we are in the superradiant regime*)Absor2=Table[{Absor[i][[1]],Absor[i][[2]]},{i,0,50}]

ListPlot[Absor2]

Expansion at infinity

Integrate[(r^2+ar^2)/((r-rplus)*(r-rminus)),r]

y=Aout(1+aout/r+bout/r^2)*Exp[I*w*rt]*r^(-1)*0+Ain(1+ain/r+bin/r^2)*Exp[-I*w*rt]*r^(-1);

aout=((I (Alm+(ar^2) (w^2)))/(2 w));
bout=-(((-2+Alm) Alm-4 I M w+2 ar (ar+Alm ar-4 I m M) (w^2)+(ar^4) (w^4))/(8 (w^2)));
ain=-((I (Alm+(ar^2) (w^2)))/(2 w));
bin->-(((-2+Alm) Alm+4 I M w+2 ar (ar+Alm ar+4 I m M) (w^2)+(ar^4) (w^4))/(8 (w^2)));

rplus=M+Sqrt[M^2-ar^2];
rminus=M-Sqrt[M^2-ar^2];rt=((r (rminus-rplus)+((ar^2)+(rminus^2)) Log[r-rminus]-((ar^2)+(rplus^2)) Log[r-rplus])/(rminus-rplus));Delta=r^2+ar^2-2*M*r;Delta1=D[Delta,r];
WaveEq=Simplify[Exp[I*w*rt]*(Delta*D[y,{r,2}]+Delta1*D[y,{r,1}]+((w*(r^2+ar^2)-m*ar)^2/Delta-Alm+2*m*ar*w-w^2*ar^2)*y)];coef=Simplify[Normal[Series[WaveEq,{r,Infinity,2}]]]

((I Ain ((Alm^2)+2 Alm (-1+(ar^2) (w^2))+w (4 I M (1+2 ar m w)+w (2 (ar^2)+8 bin+(ar^4) (w^2)))))/(2 (r^2) w))

FullSimplify[Solve[coef==0,bin]]

References

[1]E. Barausse, V. Cardoso and P. Pani, "Can environmental effects spoil precision gravitational-wave astrophysics?", Phys.Rev.D89 (2014) no.10,104059

[2]R. Vicente, V. Cardoso and J. Lopes, "The Penrose process, superradiance and ergoregion instabilities", Phys.Rev.D89 (2014) no.10,104059; e-Print:arXiv:1803.08060[gr-qc]

[3]R. Brito, V. Cardoso and P. Pani, "Superradiance", Lect.Notes Phys.906 (2015) pp.1-237

[4]https://centra.tecnico.ulisboa.pt/network/grit/files/