Written by Mauricio Richartz [mauricio.richartz@ufabc.edu.br]

Quit[];

Schwarzschild metric

This section calculates the geometric tensors for the Schwarzschild metric  and separates the Klein-Gordon equation. It is based on notebooks by J. Hartle’s (http://web.physics.ucsb.edu/~gravitybook/mathematica.html) and by E. Berti  (https://www.icts.res.in/event/page/3071). 

Comment: Each tensor name is followed by a series of u and d symbols - they indicate which indices are contravariant and which are covariant. For example, Riemannuddd is the Riemann tensor with one up index and 3 down indices.

Comment2: the indices run from μ=1 to μ=n, where n is the dimension of the spacetime.

ClearAll["Global`*"]

Metric

(* Coordinates *) x={t,r,χ,ϕ};

$Assumptions={-1≤χ≤1,0≤ ϕ≤ 2 π,r≥0};

(* The Schwarzschild metric in "Rational Polynomial" coordinates, i.e. χ=Cos[θ]: *)

gdd={{g00,0,0,0},
{0,g11,0,0},
{0,0,g22,0},
{0,0,0,g33}};
g00=-(1-((2M)/r));
g11=(((1-((2M)/r)))^(-1));
g22=((r^2)/(1-(χ^2)));
g33=(r^2)*(1-(χ^2)) ;

gdd//MatrixForm

guu=Inverse[gdd]//FullSimplify; (*inverse metric (g^(μν))*)
dim=Length[gdd]; (* number of spacetime dimensions *)
det=Det[gdd]//FullSimplify; (*metric determinant: g*) 
sqdetg=Sqrt[-det]//FullSimplify; (*square root of metric determinant: (Sqrt[-g])*)

Geometry

Comment: you don’t need to run this section if you are only interested in separating the KG equation.  
Below we calculate, in the following order, the: Christoffel symbols Γ^{i}_{jk} - > (Subscript[(Γ^i), jk]) , Riemann tensor R^{k}_{jia} -> (Subscript[(R^k), jia]), the Ricci tensor Ricci_{ja} = R^{k}_{jka} -> (Subscript[R, ja])=(Subscript[(R^k), jka]) , the Ricci scalar RicSca = g^{ja} R_{ja} -> (Subscript[R, ja])(g^ja) and the Einstein tensor Eins_{ja} = Ricci_{ja} - (1/2)*g_{ja}RicSca -> (Subscript[G, ja])=(Subscript[R, ja])-R (Subscript[g, ja]).

(* (Subscript[(Γ^i), jk]) *)
Γ=ParallelTable[Simplify[Sum[(1/2) guu[[k,α]] (D[gdd[[α,i]],x[[j]]]+D[gdd[[α,j]],x[[i]]]-D[gdd[[j,i]],x[[α]]]), 
{α,1,dim}]],
{k,1,dim},
{j,1,dim},
{i,1,dim}];
(* (Subscript[(R^k), jil]) *)
 Riemannuddd=ParallelTable[
Simplify[
D[Γ[[l,i,k]],x[[j]]]-D[Γ[[l,j,i]],x[[k]]]+
Sum[
Γ[[l,j,β]]Γ[[β,k,i]]-Γ[[l,k,β]]Γ[[β,j,i]],
{β,1,dim}]],
{l,1,dim},
{i,1,dim},
{j,1,dim},
{k,1,dim}];
(* (Subscript[R, ji]) *)
Riccidd=ParallelTable[
Simplify[
Sum[Riemannuddd[[β,j,β,i]],{β,1,dim}]],
{j,1,dim},
{i,1,dim}
];
(* (R^ja) *)
Ricciuu=ParallelTable[Simplify[
Sum[guu[[i,α]] guu[[j,β]] Riccidd[[α,β]],{α,1,dim},{β,1,dim}]],
{i,1,dim},
{j,1,dim}];
(* (Subscript[(R^i), i]) *)
RicSca=Simplify[Sum[guu[[α,β]] Riccidd[[α,β]],
{α,1,dim},
{β,1,dim}]];
(* (Subscript[G, ji]) *)
Einstein=ParallelTable[
Simplify[Riccidd[[j,i]]-(1/2)RicSca gdd[[j,i]]],
{j,1,dim},
{i,1,dim}];
(* (Subscript[R, ijkl]) *)
Riemanndddd=ParallelTable[
Simplify[Sum[gdd[[i,α]] Riemannuddd[[α,j,k,l]],{α,1,dim}]],
{i,1,dim},
{j,1,dim},
{k,1,dim},
{l,1,dim}
];
(* (R^ijkl) *)
Riemannuuuu=ParallelTable[Simplify[
Sum[guu[[α,i]]guu[[β,j]]guu[[γ,k]]guu[[δ,l]] Riemanndddd[[α,β,γ,δ]],{α,1,dim},{β,1,dim},{γ,1,dim},{δ,1,dim}]],
{i,1,dim},
{j,1,dim},
{k,1,dim},
{l,1,dim}];
(* (R^ijkl)(Subscript[R, ijkl]) *)
Kretschmann=Simplify[Sum[
Riemanndddd[[α,β,γ,δ]] Riemannuuuu[[α,β,γ,δ]],
{α,1,dim},
{β,1,dim},
{γ,1,dim},
{δ,1,dim}]];

Geometry display

Comment: you don’t need to run this section if you are only interested in separating the KG equation.  
Below we display the nonzero components of the tensors

(* generates the list of non-zero Christoffel symbols *)

listChristoffel:=Table[If[UnsameQ[Γ[[i,j,k]],0],{StringJoin[ToString[Subsuperscript["Γ",StringJoin[ToString[j],ToString[k]],i],TraditionalForm]," ="],Γ[[i,j,k]]}] ,{i,1,dim},{j,1,dim},{k,1,dim}]

(* displays the non-zero Christoffel symbols - symmetries are not taken into account *)

TableForm[Partition[DeleteCases[Flatten[listChristoffel],Null],2],TableSpacing->{3,1}]//FullSimplify

(* displays the non-zero Christoffel symbols -  taking into account the symmetry (Subscript[(Γ^i), jk]) = (Subscript[(Γ^i), kj]) *)

listChristoffel2:=Table[If[UnsameQ[Γ[[i,j,k]],0],{StringJoin[ToString[Subsuperscript["Γ",StringJoin[ToString[j],ToString[k]],i],TraditionalForm]," ="],Γ[[i,j,k]]}] ,{i,1,dim},{j,1,dim},{k,1,j}]

TableForm[Partition[DeleteCases[Flatten[listChristoffel2],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero (Subscript[(R^k), jia]) taking into account the symmetry (Subscript[(R^k), jia]) = (Subscript[(R^k), jai])  *)

listRiemannuddd:=Table[If[UnsameQ[Riemannuddd[[k,j,i,l]],0],{StringJoin[ToString[Subscript[Superscript["R",k],StringJoin[ToString[j],ToString[i],ToString[l]]],TraditionalForm]," ="],Riemannuddd[[k,j,i,l]]}] ,{k,1,dim},{j,1,dim},{i,1,dim},{l,1,i}]

(* displays the non-zero (Subscript[(R^k), jia]) components taking into account the symmetry (Subscript[(R^k), jia]) = (Subscript[(R^k), jai]) *)

TableForm[Partition[DeleteCases[Flatten[listRiemannuddd],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero (Subscript[R, kjia]) taking into account the symmetries (Subscript[R, kjia]) = (Subscript[R, kjai]) and (Subscript[R, kjia]) = (Subscript[R, jkia]) - the symmetry (Subscript[R, jkia]) = (Subscript[R, iajk ])is still in there   *)

listRiemanndddd:=Table[If[UnsameQ[Riemanndddd[[k,j,i,l]],0],{StringJoin[ToString[Subscript["R",StringJoin[ToString[k],ToString[j],ToString[i],ToString[l]]],TraditionalForm]," ="],Riemanndddd[[k,j,i,l]]}] ,{k,1,dim},{j,1,k},{i,1,dim},{l,1,i}]

(* displays the non-zero Riemanndddd (Subscript[R, kjia]) taking into account the symmetries (Subscript[R, kjia]) = (Subscript[R, kjai]) and (Subscript[R, kjia]) = (Subscript[R, jkia]) - the symmetry (Subscript[R, jkia]) = (Subscript[R, iajk])is still in there  *)

TableForm[Partition[DeleteCases[Flatten[listRiemanndddd],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero (Subscript[R, ij]) components - symmetries are not taken into account  *)

listRiccidd:=Table[If[UnsameQ[Riccidd[[i,j]],0],{StringJoin[ToString[Subscript["R",StringJoin[ToString[i],ToString[j]]],TraditionalForm]," ="],Riccidd[[i,j]]}] ,{i,1,dim},{j,1,dim}]

(* displays the non-zero (Subscript[R, ij]) components - symmetries are not taken into account *)

TableForm[Partition[DeleteCases[Flatten[listRiccidd],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero Einstein (Subscript[G, ij]) components  - symmetries are not taken into account *)

listEinstein:=Table[If[UnsameQ[Einstein[[i,j]],0],{StringJoin[ToString[Subscript["G",StringJoin[ToString[i],ToString[j]]],TraditionalForm]," ="],Einstein[[i,j]]}] ,{i,1,dim},{j,1,dim}]

(* displays the non-zero Einstein components - symmetries are not taken into account *)

TableForm[Partition[DeleteCases[Flatten[listEinstein],Null],2],TableSpacing->{3,1}]//FullSimplify

(*displays the Ricci Scalar *)

RicSca

(* displays the Kretschmann scalar *)

Kretschmann

Separation of the Klein-Gordon equation

(* You only need metric to run this - Geometry and Geometry Display are not needed *)

(* Taking into account the axisymmetry and the stationarity of the metric, we intrduce the following ansatz: *)

Phi[t,r,χ,ϕ]=Exp[-I ω t]*Exp[I m ϕ]*R[r]*S[χ];

(* Below we calculate the box operator: *)

BoxPhi=Sum[(1/sqdetg)D[sqdetg guu[[α,β]]D[Phi[t,r,χ,ϕ],x[[α]]],x[[β]]],{α,1,dim},{β,1,dim}];

(* Below we calculate the left hand side of the Klein-Gordon equation: *)

KG=BoxPhi-(μ^2)*Phi[t,r,χ,ϕ];

KG//FullSimplify

(* We now try to separate r and χ variables *)

KG01=Collect[(KG/(Phi[t,r,χ,ϕ])),{R''[r],R'[r],S''[χ],S'[χ]},FullSimplify]

(* Multiplying by (r^2), we see that the equation separates: *)

KG02=KG01*(r^2)//FullSimplify

(* Therefore, if we define λ as a separation constant, we can write: *)

(* (r^2) (μ^2) + (((r^3) (ω^2))/(-2 M+r)) + ((2 (-M+r) R'[r]+r (-2 M+r) R''[r])/(R[r])) = λ

(* and *)

(* ((m^2)/(-1+(χ^2)))-((2 χ S'[χ]+(-1+(χ^2)) S''[χ])/(S[χ])) = -λ *)

(* The separated equations are, therefore: *)

eq1a=(-(r^2) (μ^2) + (((r^3) (ω^2))/(r-2 M)))*R[r] + 2 (r-M) R'[r]+r (r-2 M) R''[r] == λ*R[r];

eq2a=-((m^2)/(1-(χ^2)))*S[χ]-2 χ S'[χ]+(1-(χ^2)) S''[χ] == -λ*S[χ];

(* Note that the second equation is the associated Legendre equation. If we require that its solutions are well-behaved (i.e. non-singular) at χ=-1 and χ=1 (equivalently, at θ=0 and θ=π), the separation constant must be λ=L*(L+1), where L is an integer. *)

(* We can, therefore, rewrite the equations as: *)

eq1b=(-(r^2) (μ^2)*(r-2 M) + (r^3) (ω^2))*R[r] + 2*(r-M)*(r-2 M)*R'[r]+r*(((r-2 M))^2)*R''[r] == λ*(r-2 M)*R[r]

eq2b=(1-(χ^2)) S''[χ]-2 χ S'[χ]+(λ-((m^2)/(1-(χ^2))))*S[χ];

(* Transformation of the radial equation to a Schroedinger-like equation: *)

(*Changing R[r] to ((ψ[r])/r) *) eq1c=eq1b/.{R->Function[r,((ψ[r])/r)]}//FullSimplify

(* Introducing the tortoise coordinate y, so that dy/dr = (r^2)/((r^2)-2Mr) :=g'[r]. Therefore dψ/dr = (dψ/dy)*(dy/dr) and (d^2)ψ/(dr^2) = (dψ/dy)*((d^2)y/(dr^2))+((d^2)ψ/(dy^2))*(((dy/dr))^2)      *)

eq1c/.ψ->(ψ[g[#]]&)/.{g'[r]->((r^2)/((r^2)-2*M*r)),g''[r]-> D[((r^2)/((r^2)-2*M*r)),r]}/.{g[r]->y}//FullSimplify

λ=L*(L+1);

(* Comparing the expression above with (((d^2)ψ)/(dy^2))+((ω^2)-V(r))*ψ = 0, we find: *)

V[r_]=-(((2 M-r) (2 M+r λ+(r^3) (μ^2)))/(r^4));

(* Explicit expression for the tortoise coordinate: *)

Integrate[(r^2)/((r^2)-2*M*r),r,GeneratedParameters->C]//Simplify

(* Setting the constant (C[1]) = -2*M*Log[2M], we have: *)

y[r_]=r+2 M Log[((-2 M+r)/(2 M))];

Kerr metric

This section calculates the geometric tensors for the Kerr metric and separates the Klein-Gordon equation. It is based on notebooks by J. Hartle’s (http://web.physics.ucsb.edu/~gravitybook/mathematica.html) and by E. Berti  (https://www.icts.res.in/event/page/3071). 

Comment: Each tensor name is followed by a series of u and d symbols - they indicate which indices are contravariant and which are covariant. For example, Riemannuddd is the Riemann tensor with one up index and 3 down indices.

Comment2: the indices run from μ=1 to μ=n, where n is the dimension of the spacetime.

ClearAll["Global`*"]

Metric

(* Coordinates *) x={t,r,χ,ϕ};

$Assumptions={-1≤χ≤1,0≤ ϕ≤ 2 π,r≥0};

(* The Kerr metric in Boyer-Lindquist "Rational Polynomial" coordinates, i.e. χ=Cos[θ]: *)

Δ=(r^2)-2*M*r+(a^2);Σ=(r^2)+(a^2)*(χ^2);

gdd={{g00,0,0,g03},
{0,g11,0,0},
{0,0,g22,0},
{g30,0,0,g33}};
g00=-(1-((2 M r)/Σ));
g01=-((2 a M r (1-(χ^2)))/Σ);
g11=+(Σ/Δ);
g22=+(Σ/(1-(χ^2)));
g03=-((2 a M r (1-(χ^2)))/Σ);g30=g03;
g33=+((r^2)+(a^2)+((2 M r (a^2)(1-(χ^2)))/Σ))*(1-(χ^2));

gdd//MatrixForm

guu=Inverse[gdd]//FullSimplify; (*inverse metric (g^(μν))*)
dim=Length[gdd]; (* number of spacetime dimensions *)
det=Det[gdd]//FullSimplify; (*metric determinant: g*) 
sqdetg=Sqrt[-det]//FullSimplify; (*square root of metric determinant: (Sqrt[-g])*)

Geometry

Comment: you don’t need to run this section if you are only interested in separating the KG equation.  
Below we calculate, in the following order, the: Christoffel symbols Γ^{i}_{jk} - > (Subscript[(Γ^i), jk]) , Riemann tensor R^{k}_{jia} -> (Subscript[(R^k), jia]), the Ricci tensor Ricci_{ja} = R^{k}_{jka} -> (Subscript[R, ja])=(Subscript[(R^k), jka]) , the Ricci scalar RicSca = g^{ja} R_{ja} -> (Subscript[R, ja])(g^ja) and the Einstein tensor Eins_{ja} = Ricci_{ja} - (1/2)*g_{ja}RicSca -> (Subscript[G, ja])=(Subscript[R, ja])-R (Subscript[g, ja]).

(* (Subscript[(Γ^i), jk]) *)
Γ=ParallelTable[Simplify[Sum[(1/2) guu[[k,α]] (D[gdd[[α,i]],x[[j]]]+D[gdd[[α,j]],x[[i]]]-D[gdd[[j,i]],x[[α]]]), 
{α,1,dim}]],
{k,1,dim},
{j,1,dim},
{i,1,dim}];
(* (Subscript[(R^k), jil]) *)
 Riemannuddd=ParallelTable[
Simplify[
D[Γ[[l,i,k]],x[[j]]]-D[Γ[[l,j,i]],x[[k]]]+
Sum[
Γ[[l,j,β]]Γ[[β,k,i]]-Γ[[l,k,β]]Γ[[β,j,i]],
{β,1,dim}]],
{l,1,dim},
{i,1,dim},
{j,1,dim},
{k,1,dim}];
(* (Subscript[R, ji]) *)
Riccidd=ParallelTable[
Simplify[
Sum[Riemannuddd[[β,j,β,i]],{β,1,dim}]],
{j,1,dim},
{i,1,dim}
];
(* (R^ja) *)
Ricciuu=ParallelTable[Simplify[
Sum[guu[[i,α]] guu[[j,β]] Riccidd[[α,β]],{α,1,dim},{β,1,dim}]],
{i,1,dim},
{j,1,dim}];
(* (Subscript[(R^i), i]) *)
RicSca=Simplify[Sum[guu[[α,β]] Riccidd[[α,β]],
{α,1,dim},
{β,1,dim}]];
(* (Subscript[G, ji]) *)
Einstein=ParallelTable[
Simplify[Riccidd[[j,i]]-(1/2)RicSca gdd[[j,i]]],
{j,1,dim},
{i,1,dim}];
(* (Subscript[R, ijkl]) *)
Riemanndddd=ParallelTable[
Simplify[Sum[gdd[[i,α]] Riemannuddd[[α,j,k,l]],{α,1,dim}]],
{i,1,dim},
{j,1,dim},
{k,1,dim},
{l,1,dim}
];
(* (R^ijkl) *)
Riemannuuuu=ParallelTable[Simplify[
Sum[guu[[α,i]]guu[[β,j]]guu[[γ,k]]guu[[δ,l]] Riemanndddd[[α,β,γ,δ]],{α,1,dim},{β,1,dim},{γ,1,dim},{δ,1,dim}]],
{i,1,dim},
{j,1,dim},
{k,1,dim},
{l,1,dim}];
(* (R^ijkl)(Subscript[R, ijkl]) *)
Kretschmann=Simplify[Sum[
Riemanndddd[[α,β,γ,δ]] Riemannuuuu[[α,β,γ,δ]],
{α,1,dim},
{β,1,dim},
{γ,1,dim},
{δ,1,dim}]];

Geometry display

Comment: you don’t need to run this section if you are only interested in separating the KG equation.  
Below we display the nonzero components of the tensors

(* generates the list of non-zero Christoffel symbols *)

listChristoffel:=Table[If[UnsameQ[Γ[[i,j,k]],0],{StringJoin[ToString[Subsuperscript["Γ",StringJoin[ToString[j],ToString[k]],i],TraditionalForm]," ="],Γ[[i,j,k]]}] ,{i,1,dim},{j,1,dim},{k,1,dim}]

(* displays the non-zero Christoffel symbols - symmetries are not taken into account *)

TableForm[Partition[DeleteCases[Flatten[listChristoffel],Null],2],TableSpacing->{3,1}]//FullSimplify

(* displays the non-zero Christoffel symbols -  taking into account the symmetry (Subscript[(Γ^i), jk]) = (Subscript[(Γ^i), kj]) *)

listChristoffel2:=Table[If[UnsameQ[Γ[[i,j,k]],0],{StringJoin[ToString[Subsuperscript["Γ",StringJoin[ToString[j],ToString[k]],i],TraditionalForm]," ="],Γ[[i,j,k]]}] ,{i,1,dim},{j,1,dim},{k,1,j}]

TableForm[Partition[DeleteCases[Flatten[listChristoffel2],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero (Subscript[(R^k), jia]) taking into account the symmetry (Subscript[(R^k), jia]) = (Subscript[(R^k), jai])  *)

listRiemannuddd:=Table[If[UnsameQ[Riemannuddd[[k,j,i,l]],0],{StringJoin[ToString[Subscript[Superscript["R",k],StringJoin[ToString[j],ToString[i],ToString[l]]],TraditionalForm]," ="],Riemannuddd[[k,j,i,l]]}] ,{k,1,dim},{j,1,dim},{i,1,dim},{l,1,i}]

(* displays the non-zero (Subscript[(R^k), jia]) components taking into account the symmetry (Subscript[(R^k), jia]) = (Subscript[(R^k), jai]) *)

TableForm[Partition[DeleteCases[Flatten[listRiemannuddd],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero (Subscript[R, kjia]) taking into account the symmetries (Subscript[R, kjia]) = (Subscript[R, kjai]) and (Subscript[R, kjia]) = (Subscript[R, jkia]) - the symmetry (Subscript[R, jkia]) = (Subscript[R, iajk ])is still in there   *)

listRiemanndddd:=Table[If[UnsameQ[Riemanndddd[[k,j,i,l]],0],{StringJoin[ToString[Subscript["R",StringJoin[ToString[k],ToString[j],ToString[i],ToString[l]]],TraditionalForm]," ="],Riemanndddd[[k,j,i,l]]}] ,{k,1,dim},{j,1,k},{i,1,dim},{l,1,i}]

(* displays the non-zero Riemanndddd (Subscript[R, kjia]) taking into account the symmetries (Subscript[R, kjia]) = (Subscript[R, kjai]) and (Subscript[R, kjia]) = (Subscript[R, jkia]) - the symmetry (Subscript[R, jkia]) = (Subscript[R, iajk])is still in there  *)

TableForm[Partition[DeleteCases[Flatten[listRiemanndddd],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero (Subscript[R, ij]) components - symmetries are not taken into account  *)

listRiccidd:=Table[If[UnsameQ[Riccidd[[i,j]],0],{StringJoin[ToString[Subscript["R",StringJoin[ToString[i],ToString[j]]],TraditionalForm]," ="],Riccidd[[i,j]]}] ,{i,1,dim},{j,1,dim}]

(* displays the non-zero (Subscript[R, ij]) components - symmetries are not taken into account *)

TableForm[Partition[DeleteCases[Flatten[listRiccidd],Null],2],TableSpacing->{3,1}]//FullSimplify

(* generates the list of non-zero Einstein (Subscript[G, ij]) components  - symmetries are not taken into account *)

listEinstein:=Table[If[UnsameQ[Einstein[[i,j]],0],{StringJoin[ToString[Subscript["G",StringJoin[ToString[i],ToString[j]]],TraditionalForm]," ="],Einstein[[i,j]]}] ,{i,1,dim},{j,1,dim}]

(* displays the non-zero Einstein components - symmetries are not taken into account *)

TableForm[Partition[DeleteCases[Flatten[listEinstein],Null],2],TableSpacing->{3,1}]//FullSimplify

(*displays the Ricci Scalar *)

RicSca

(* displays the Kretschmann scalar *)

Kretschmann

Separation of the Klein-Gordon equation

(* You only need metric to run this - Geometry and Geometry Display are not needed *)

(* Taking into account the axisymmetry and the stationarity of the metric, we introduce the following ansatz: *)

Phi[t,r,χ,ϕ]=Exp[-I ω t]*Exp[I m ϕ]*R[r]*S[χ];

(* Below we calculate the box operator: *)

BoxPhi=Sum[(1/sqdetg)D[sqdetg guu[[α,β]]D[Phi[t,r,χ,ϕ],x[[α]]],x[[β]]],{α,1,dim},{β,1,dim}];

(* Below we calculate the left hand side of the Klein-Gordon equation: *)

KG=BoxPhi-(μ^2)*Phi[t,r,χ,ϕ];

KG//FullSimplify

(* We now try to separate r and χ variables *)

KG01=Collect[(KG/(Phi[t,r,χ,ϕ])),{R''[r],R'[r],S''[χ],S'[χ]},FullSimplify]

(* Multiplying by (r^2)+(a^2) (χ^2), we obtain: *)

KG02=KG01*((r^2)+(a^2) (χ^2))//FullSimplify

(*If the following term: 

(1/(((a^2)+r (-2 M+r)) (-1+(χ^2))))((m^2) (-2 M r+(r^2)+(a^2) (χ^2))-4 a m M r (-1+(χ^2)) ω+(-1+χ) (1+χ) (-((a^2)+r (-2 M+r)) (μ^2) ((r^2)+(a^2) (χ^2))+((r^4)+(a^4) (χ^2)+(a^2) r (2 M+r+(-2 M+r) (χ^2))) (ω^2)))

can be written as f(r) + g(χ),then the equation separates. Let us check that it is indeed the case. *)

H[r_,χ_]=(1/(((a^2)+r (-2 M+r)) (-1+(χ^2))))((m^2) (-2 M r+(r^2)+(a^2) (χ^2))-4 a m M r (-1+(χ^2)) ω+(-1+χ) (1+χ) (-((a^2)+r (-2 M+r)) (μ^2) ((r^2)+(a^2) (χ^2))+((r^4)+(a^4) (χ^2)+(a^2) r (2 M+r+(-2 M+r) (χ^2))) (ω^2)));

D[H[r,χ],r,χ]//Simplify

(* Therefore: *)

g[χ]=Integrate[D[H[r,χ],χ]//Simplify,χ]

f[r_]=H[r,χ]-g[χ]//FullSimplify

(* Note that f does is independent of χ, as expected *)

(* Therefore, if we define λ as a separation constant, we can write: *)

(* ((2 (-M+r) R'[r]+((a^2)+r (-2 M+r)) R''[r])/(R[r])) + f[r] = λ

(* and *)

(* -((2 χ S'[χ]+(-1+(χ^2)) S''[χ])/(S[χ]))+g[χ] = -λ *)

(* The separated equations are, therefore: *)

eq1a=(((-4 a m M r ω+(r^3) ((2 M-r) (μ^2)+r (ω^2))+(a^2) ((m^2)+r (-r (μ^2)+(2 M+r) (ω^2))))/(Δ[r])))*R[r] + Δ'[r]*R'[r]+Δ[r]*R''[r] == λ*R[r];

eq2a=(-((m^2)/(1-(χ^2)))-(a^2) (χ^2) ((μ^2)-(ω^2)))*S[χ]-2 χ S'[χ]+(1-(χ^2)) S''[χ] == -λ*S[χ];

(* We can, therefore, rewrite the equations as: *)

eq1b=(Δ^2)*R''[r]+Δ*D[Δ,r]*R'[r]+(h[r]-λ*Δ)*R[r] ==0

h[r_]=-4 a m M r ω+(r^3) ((2 M-r) (μ^2)+r (ω^2))+(a^2) ((m^2)+r (-r (μ^2)+(2 M+r) (ω^2)));

eq2b=(1-(χ^2)) S''[χ]-2 χ S'[χ]+(λ-((m^2)/(1-(χ^2))))*S[χ];

(* Transformation of the radial equation to a Schroedinger-like equation: *)

(*Changing R[r] to ((ψ[r])/r) *) eq1c=eq1b/.{R->Function[r,((ψ[r])/(Sqrt[(r^2)+(a^2)]))]}//FullSimplify

(* Introducing the tortoise coordinate y, so that dy/dr = (r^2)/Δ[r] :=g'[r]. Therefore dψ/dr = (dψ/dy)*(dy/dr) and (d^2)ψ/(dr^2) = (dψ/dy)*((d^2)y/(dr^2))+((d^2)ψ/(dy^2))*(((dy/dr))^2)      *)

eq1c/.ψ->(ψ[g[#]]&)/.{g'[r]->(((r^2)+(a^2))/Δ),g''[r]-> D[(((r^2)+(a^2))/Δ),r]}/.{g[r]->y}//FullSimplify

(* Comparing the expression above with (((d^2)ψ)/(dy^2))+((ω^2)-V(r,ω))*ψ = 0, we find: *)

V[r_,ω_]=(ω^2)-(1/((((r^2)+(a^2)))^4))(-(r^4) (-2 M+r) (2 M+r λ+(r^3) (μ^2))-4 (a^5) m M r ω-8 (a^3) m M (r^3) ω-4 a m M (r^5) ω+(r^8) (ω^2)+(a^6) (-1+(m^2)-λ-(r^2) (μ^2)+r (2 M+r) (ω^2))+(a^2) (r^2) (-8 (M^2)+(r^2) (-1+(m^2)-3 λ+3 (r^2) (-(μ^2)+(ω^2)))+2 M r (2+2 λ+(r^2) (2 (μ^2)+(ω^2))))+(a^4) r (2 (-1+(m^2)) r-3 r (λ+(r^2) (μ-ω) (μ+ω))+2 M (3+λ+(r^2) ((μ^2)+2 (ω^2)))));

V[r,ω]//FullSimplify

(* No limite a->0, recuperamos a solucao de Schwarzschild *)

V[r,ω]/.{a->0}//Simplify

Angular eigenvalues of the spin-0 spheroidal harmonics

This section calculates the angular eigenvalues of the spheroidal equation associated with massless scalar fields around a Kerr black hole. The numerical method used is the continued fraction method.

ClearAll["Global`*"]

Angular eigenvalues - CFM

Behaviour near the singular points

eq=D[(1-(χ^2))*D[S[χ],χ],χ]+(A+(γ^2)*(χ^2)-((m^2)/(1-(χ^2))))*S[χ]

(* Solutions near χ=-1 *)

Clear[S,C1];

S[χ_]=(((1+χ))^C1);

Series[(eq/(((1+χ))^C1))//Simplify,{χ,-1,-1}]//FullSimplify

(* Therefore, we have C1=m/2 or C1=-m/2;

(* Solution 2 *)

Clear[S,C2];

S[χ_]=(((1-χ))^C2);

Series[(eq/(((1-χ))^C2))//Simplify,{χ,1,-1}]//FullSimplify

(* Therefore, we have C2=m/2 or C2=-m/2;

Algebra of the angular equation

Clear[eq,eq2,S,A,γ,m,χ,z,f];

$Assumptions={m∈Reals};

(* Angular equation of the angular part S of a scalar field around Kerr: ψ=R[r]*S[θ]*(e^(+imϕ))*(e^(-iωt)): *)

eq[χ_]=D[(1-(χ^2))*D[S[χ],χ],χ]+(A+(γ^2)*(χ^2)-((m^2)/(1-(χ^2))))*S[χ]

S[χ_]:=Exp[γ*χ]*(((1+χ))^J1)*(((1-χ))^J2)*f[χ] (* We want S[χ] to be well behaved - i.e. regular - when χ->-1 and χ-> +1  *)

(*The regular solutions have J1 = J2 = Abs[m/2]: *)

J1=Abs[(m/2)];J2=Abs[(m/2)];

(* We want to find a series solution around χ=-1. Therefore, we change to a new variable z = χ + 1 *)

z[χ_]:=1+χ;χ[w_]:=InverseFunction[z][w];

(* The change of variable is performed below: *)

eq2[z_]=eq[χ]/.{f[χ]->f[z],f'[χ]->f'[z]*z'[χ],f''[χ]->f'[z]*z''[χ]+f''[z]*(((z'[χ]))^2)} /.χ->χ[z]//FullSimplify

(* Eliminating the overall factor that appears above, we find: *)

eq3[z_]:=(-A+(m^2)+(-2+2 z-γ) γ+Abs[m]+2 (-1+z) γ Abs[m]) f[z]+2 (-1+z+(-2+z) z γ+(-1+z) Abs[m]) f'[z]+(-2+z) z f''[z];

Let us try to find a solution for eq3[z] that has the for f[z]=∑ (Subscript[b, j])(z^j). Since eq3[z] is a linear ordinary differential equation, instead of substituting f[z]=∑(Subscript[b, j])(z^j) into eq3[z], we can substitute (Subscript[f, j])[z]=(Subscript[b, j])(z^j) into eq3[z] and only the sum over j:

fj[z_]:=b[j]*(z^j);

eq3[z_,j_]=eq3[z]//.{f->Function[z,fj[z]]}//FullSimplify

(* A "trick" to obtain the recurrence relations in Mathematica is the following: *)

eqfinal[z_]=(z^(-n))*Sum[eq3[z,j],{j,n-3,n+3}]//Simplify; (* Instead of summing from j=0 to j=∞, we sum from j = n-3 up to j=n+3. We later multiply the result by (x^(-n)) to facilitate the algebraic manipulation in Mathematica *)

F[j_]=SeriesCoefficient[eqfinal[z],{z,0,0}]

(* Collecting the terms: *)

Coefficient[F[n],b[n-1]]//FullSimplify

Coefficient[F[n],b[n]]//FullSimplify

Coefficient[F[n],b[n+1]]//FullSimplify

(* We therefore have: *)

αang[n_]:=-2 (1+n) (1+n+Abs[m]);

βang[n_]:=-A+(m^2)+n+(n^2)-4 n γ-γ (2+γ)+(1+2 n-2 γ) Abs[m];

γang[n_]:=2 γ (n+Abs[m]);

Angular eigenvalues and comparison with the Mathematica native function

Clear[αang,βang,γang,Nol,n,m,γ,n,l,eigenv,μ,ω]

(* The coefficients of the recursion relation *)

αang[n_]:=-2 (1+n) (1+n+Abs[m]);

βang[n_]:=-A+(m^2)+n+(n^2)-4 n γ-γ (2+γ)+(1+2 n-2 γ) Abs[m];

γang[n_]:=2 γ (n+Abs[m]);

(*Nollert's improvement below - setting it to zero corresponds to not using the improvement *)

Nol[n_]=0;

EQI[Ninv_?NumericQ,Nterms_?NumericQ]:=Block[{Rn1,Rn2},Rn1=Nol[Nterms];Rn2=0;Do[Rn1=((-γang[n])/(βang[n]+αang[n]*Rn1)),{n,Nterms,1+Ninv,-1}];Do[Rn2=((-γang[n])/(βang[n-1]+αang[n-2]*Rn2)),{n,1,Ninv}];βang[Ninv]+Rn2*αang[Ninv-1]+Rn1*αang[Ninv]];

eigenv[Ninv_?NumericQ,Nterms_?NumericQ,γγ_?NumericQ,ll_?NumericQ,mm_?NumericQ]:=Block[{$MinPrecision=20,γ=γγ,l=ll,m=mm,Ainit,eqn},Ainit=l*(l+1);eqn[A_?NumericQ]=EQI[Ninv,Nterms];A/.FindRoot[{eqn[A]==0},{A,Ainit},WorkingPrecision->20]]

eigenv[0,50,(1/2),1,1]

(*Calculation using the native function in Mathematica: *)

mateigen[γ_,l_,m_]:=-(γ^2)+SpheroidalEigenvalue[l,m,I*γ];

N[mateigen[(1/2),1,1],20]

Scattering problem in the frequency domain

This notebook calculates the reflection coefficient associated with scattering processes of massless scalar fields around a Kerr black hole

ClearAll["Global`*"]

Scattering problem in the frequency domain - massless scalar field around Kerr

(* We assume M=1. The scalar field is assumed massless (μ=0). x is the radial coordinate so that the location of the horizon is xp. ϵ is a small parameter and x1 is far away from the horizon. l is relatd to the angular eigenvalue and m is the azimuthal number. *)

ϵ=(10^(-5));xp:=1+(Sqrt[1-(a^2)]);x0:=xp+ϵ;x1=500;l=1;m=1;

(* The differential equation is: *)

eq[x_]:=((Δ[x])^2)*D[R[x],x,x]+Δ[x]*(2*x-2)*D[R[x],x]+(((K[x])^2)-λ*Δ[x])*R[x];

Δ[x_]:=(x^2)-2*x+(a^2);

K[x_]:=ω*((x^2)+(a^2))-a*m;  λ:=-2*a*m*ω+(ω^2)*(a^2)+A;

A:=SpheroidalEigenvalue[l,m,I*a*ω]-(((a*ω))^2); (* This is how the separation constant relates to the Mathematica-defined spheroidal eigenvalues*)

(* The solution near r = rplus is R = (((x-xp))^b). Therefore, we have the following initial conditions:  *)

b:=((I(a m-2 xp*ω))/(2 (-1+xp)));

iv1:=(((ϵ))^b); (*wave function at x=xp+ϵ*)
iv2:=(((ϵ))^(b-1))*b; (* derivative of the wave function at x=xp+ϵ*)

(*The function below solves the differential equation for a given frequency ω and a given spin a *)

sol[ω0_,a0_]:=Block[{ω=ω0,a=a0},NDSolve[{eq[x]==0,R[x0]==iv1,R'[x0]==iv2},R,{x,x0,x1}]]

(* The function below extracts the wave function R1 and its derivative R1p at the endpoint x=x1*)

extract[ω_,a_]:=(R1=(R[x1]/.sol[ω,a])[[1]];R1p=(R'[x1]/.sol[ω,a])[[1]];);

(* We calculate the reflection coefficient from R1 and R1p using the asymptotic solutions - see below*)

f1[x_]:=Yin*Exp[-I*ω*x]*(x^(-1-2*I*ω))+Yout*Exp[I*ω*x]*(x^(-1+2*I*ω));

FullSimplify[Solve[{f1[xx1]==R1,(f1'[x]/.x->xx1)==R1p},{Yin,Yout}]]

(*Calculating the reflection coefficient for a given frequency ω and a given spin a*) ref[ω_,a_]:=(extract[ω,a];Yin=(((E^(I x1 ω)) (x1^(1+2 I ω)) (I (R1+R1p x1)+R1 (2+x1) ω))/(2 (2+x1) ω));Yout=(((E^(-I x1 ω)) (x1^(1-2 I ω)) (-I (R1+R1p x1)+R1 (2+x1) ω))/(2 (2+x1) ω));Abs[((((Yout/Yin)))^2)]);

(* Example: reflection coefficient for ω = 0.4 and a = 0.99 *)

ref[4/10,99/100]

(* Generates a table of reflection coefficients for a given spin a *)

tab[a_]:=Monitor[Table[{ω,ref[ω,a]},{ω,0.001,1.0*((a*m)/(2*(1+(Sqrt[1-(a^2)])))),((a*m)/(100*(1+(Sqrt[1-(a^2)]))))}],{ω,a}];

t0=tab[0.999];t1=tab[0.99];t2=tab[0.95];t3=tab[0.9];t4=tab[0.8];

(* Plot of the reflection coefficients for a = 0.8, a=0.9, a=0.95, a=0.99, and a=0.999 *)

f1=ListPlot[{t0,t1,t2,t3,t4},PlotLegends->{0.999,0.99,0.95,0.9,0.8}]

(* The same plots as before, using now the MaTeX package for LaTeX commands - https://github.com/szhorvat/MaTeX *)

<<MaTeX`

texStyle={FontFamily->"Latin Modern Roman",FontSize->18,GrayLevel[0]};

fig=ListLinePlot[{t0,t1,t2,t3,t4},PlotStyle->Thickness[0.005],PlotRange->Automatic,PlotTheme->{"Detailed"},AspectRatio->1/2,ImageSize->800,FrameLabel->{MaTeX["\\omega/m\\Omega",FontSize->20],MaTeX["\\mathrm{R}_{\\omega}",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["a=0.999",FontSize->20],MaTeX["a=0.99",FontSize->20],MaTeX["a=0.95",FontSize->20],MaTeX["a=0.9",FontSize->20],MaTeX["a=0.8",FontSize->20]},{{0.1,0.7}}]]

Radial eigenvalues - finding the recurrence relation for Leaver’s method

This section determined the 3-term recurrence relation associated with the continued fraction method for the radial part of a massive scalar field around a Kerr black hole. The calculations for the angular part are performed in the section “Angular eigenvalues....”

ClearAll["Global`*"]

Massive scalar field in Kerr

(* We work here with the massless scalar field in a Kerr background. It is straightforward to generalize the calculations here to the general massless spin-s fields using the Teukolsky equations *)

(* Instead of a and M, I will use rp=rplus= M + (Sqrt[(M^2)-(a^2)])and rm=rminus= M - (Sqrt[(M^2)-(a^2)]). To convert back and find (M,a) in terms of (rp,rm), we use the relations 2 M=rp+rm; (a^2)=rp*rm *)

Δ[r_]=(r-rp)*(r-rm);

K[r_]=ω*((r^2)+rp*rm)-a*m;

k=(Sqrt[(μ^2)-(ω^2)]);(* note that the definition of k is slightly different than in the lectures *)

(* Following the notation from appendix A of 1312.2326.pdf or from 1212.1477.pdf for λ and Alm, the equations are: *)

eqrad:=((Δ[r])^2)*D[R[r],r,r]+Δ[r]*D[Δ[r],r]*D[R[r],r]+(((K[r])^2)-(λ+(μ^2)*(r^2))*(r-rp)*(r-rm))*R[r](* = 0 *) ;

ang=D[S[θ],θ,θ]+((Cos[θ])/(Sin[θ]))*D[S[θ],θ]+(-(k^2)*(a^2)*((Cos[θ])^2)-((m^2)/((Sin[θ])^2))+Alm)*S[θ](* = 0 *)

(* λ = Alm + (a^2)*(ω^2) - 2*a*m*ω; *) (* Convention in 1312.2326.pdf *)

(* Alm=(k^2)*(a^2)+SpheroidalEigenvalue[l,m,k*a] *) (* Alm in terms of the function implemented in Mathematica  - Leaver's method can also be used to find it. *)

(* S[θ]=SpheroidalPS[l,m,k*a,Cos[θ]] *)  (* The angular function using the function implemented in Mathematica - Leaver's method can also be used to find it. *)

Solutions at infinity

Clear[R,AI,BI,fI,a1,a2];

R[r_]:=(r^AI)*Exp[BI*r]*(1+(a1/r)+(a2/(r^2))(* + ... *)); (* Let's try this ansatz to find the asymptotic solution. We will only need AI and BI in the following, but a1, a2, ... can also be determined *)

Series[(eqrad/((r^AI)*Exp[BI*r])),{r,∞,-1}]//FullSimplify

(* From the (r^4) term we find BI, from the (r^3) term we find AI:  *)

Solve[{((BI^2)-(μ^2)+(ω^2))==0,2 (1+AI) BI+(BI^2) (a1-2 (rm+rp))+(-a1+rm+rp) (μ^2)+a1 (ω^2)==0},{BI,AI}]//FullSimplify

(* The two pairs found above can be written as: *)

BI1=-k;AI1=-1-(((rm+rp)*((μ^2)-2 (ω^2)))/(2*k));
(*and*)
BI2=k;AI2=-1+(((rm+rp)*((μ^2)-2 (ω^2)))/(2*k));

Solutions near the event horizon

Clear[R,AH,a1,a2];

R[r_]=(((r-rp))^AH)*(1+a1*(r-rp)+a2*(((r-rp))^2)(* + ... *));(* Let's try this ansatz to find the asymptotic solution. We will only need AH in the following, but a1, a2, ... can also be determined *)

Series[(eqrad/(((r-rp))^AH)),{r,rp,2}]//FullSimplify

Solve[(AH^2) (((rm-rp))^2)+(((a m-rp (rm+rp) ω))^2)==0,AH](* From the zeroth order term, we find AH *)

AH1=((I (-a m+rp (rm+rp) ω))/(rm-rp));(* Corresponds to the ingoing solution at the horizon *)

AH2=((I (-a m+rp (rm+rp) ω))/(-rm+rp));(* Corresponds to the outgoing solution at the horizon *)

Leaver’s method algebra

ClearAll["Global`*"]

(* Transformations between a,M and rp,rm:  2 M=rp+rm; (a^2)=rp*rm; rp = M + (Sqrt[(M^2)-(a^2)]); rm = M - (Sqrt[(M^2)-(a^2)]); *)

$Assumptions={rp>rm>0,M>a≥0};

Δ[r_]=(r-rp)*(r-rm);

K[r_]=ω*((r^2)+rp*rm)-a*m;

k=(Sqrt[(μ^2)-(ω^2)]);

(* Following the notation from appendix A of 1312.2326.pdf for λ and Alm, the Teukolsky eqs are: *)

eqrad=((Δ[r])^2)*D[R[r],r,r]+Δ[r]*D[Δ[r],r]*D[R[r],r]+(((K[r])^2)-(λ+(μ^2)*(r^2))*(r-rp)*(r-rm))*R[r];

R[r_]=(((r-rm))^(AI-AH))*(((r-rp))^AH)*Exp[-k*r]*f[r];(* k, AI and AH were defined/determined in the previous section; f[r] is a power series in ((r-rp)/(r-rm))*)

AH=((I (-a m+rp (rm+rp) ω))/(rm-rp));AI=-1-(((rm+rp)* ((μ^2)-2 (ω^2)))/(2*k));

(*change of coordinates: x = g[r] = ((r-rp)/(r-rm)) ;*)

g[r_]=((r-rp)/(r-rm));

gin[x_]=InverseFunction[g][x]; (*The inverse transformation r=gin[x], giving r in terms of x *);

(* Rewrites the radial Teukolsky equation in terms of x = ((r-rp)/(r-rm)) instad of r:  *)

FullSimplify[eqrad/.{f[r]->f[x],f'[r]->f'[x]*g'[r],f''[r]->f'[x]*g''[r]+f''[x]*((g'[r])^2)}/.r->gin[x]]

(*Since the expression above must be equal to zero, we simplify it by eliminating the overall factor that multiplies it: *)

eqrad2=(((rm-rp) (-4 (rp^2)+(((rm+rp))^2) x) (μ^4)-4 (rm-rp) (-1+x-λ) (ω^2)-8 I rp (rm+rp) (-1+x) (ω^3)-8 rp (rm+rp) (rp (-2+x)+rm x) (ω^4)+2 (rm^2) (μ^2) (Sqrt[(μ-ω) (μ+ω)])-8 rm rp (μ^2) (Sqrt[(μ-ω) (μ+ω)])+6 (rp^2) (μ^2) (Sqrt[(μ-ω) (μ+ω)])+4 (rm^2) x (μ^2) (Sqrt[(μ-ω) (μ+ω)])-4 (rp^2) x (μ^2) (Sqrt[(μ-ω) (μ+ω)])+4 I (rm^2) rp (μ^2) ω (Sqrt[(μ-ω) (μ+ω)])-8 I rm (rp^2) (μ^2) ω (Sqrt[(μ-ω) (μ+ω)])-12 I (rp^3) (μ^2) ω (Sqrt[(μ-ω) (μ+ω)])+4 I (rm^2) rp x (μ^2) ω (Sqrt[(μ-ω) (μ+ω)])+8 I rm (rp^2) x (μ^2) ω (Sqrt[(μ-ω) (μ+ω)])+4 I (rp^3) x (μ^2) ω (Sqrt[(μ-ω) (μ+ω)])+8 rm rp (ω^2) (Sqrt[(μ-ω) (μ+ω)])-8 (rp^2) (ω^2) (Sqrt[(μ-ω) (μ+ω)])-8 (rm^2) x (ω^2) (Sqrt[(μ-ω) (μ+ω)])+8 (rp^2) x (ω^2) (Sqrt[(μ-ω) (μ+ω)])+16 I rm (rp^2) (ω^3) (Sqrt[(μ-ω) (μ+ω)])+16 I (rp^3) (ω^3) (Sqrt[(μ-ω) (μ+ω)])-8 I (rm^2) rp x (ω^3) (Sqrt[(μ-ω) (μ+ω)])-16 I rm (rp^2) x (ω^3) (Sqrt[(μ-ω) (μ+ω)])-8 I (rp^3) x (ω^3) (Sqrt[(μ-ω) (μ+ω)])+4 (μ^2) ((rm-rp) (-1+x-λ)+2 I rp (rm+rp) (-1+x) ω+rp (-rp (3 rm+5 rp)+2 (((rm+rp))^2) x) (ω^2))+4 a m (2 (ω^2) (-I-2 rp ω-2 I rp (Sqrt[(μ-ω) (μ+ω)])+x (I+(rm+rp) (ω+I (Sqrt[(μ-ω) (μ+ω)]))))+(μ^2) (2 I+4 rp ω-I ((rm-3 rp) (Sqrt[(μ-ω) (μ+ω)])+x (2+(rm+rp) (-2 I ω+(Sqrt[(μ-ω) (μ+ω)]))))))) f[x]+4 ((-2 I a m (((-1+x))^2) (μ-ω) (μ+ω)+(rm^2) x (Sqrt[(μ-ω) (μ+ω)]) ((1+x) (μ^2)-2 x (ω^2))+rm (μ-ω) (μ+ω) (1-4 x+3 (x^2)+2 I rp (((-1+x))^2) ω-4 rp x (Sqrt[(μ-ω) (μ+ω)]))+rp (-(-1+x) (-1+3 x) (μ-ω) (μ+ω)+2 I rp (((-1+x))^2) (μ-ω) ω (μ+ω)+3 rp x (μ^2) (Sqrt[(μ-ω) (μ+ω)])-rp (x^2) (μ^2) (Sqrt[(μ-ω) (μ+ω)])-4 rp x (ω^2) (Sqrt[(μ-ω) (μ+ω)])+2 rp (x^2) (ω^2) (Sqrt[(μ-ω) (μ+ω)]))) f'[x]+(rm-rp) (((-1+x))^2) x (μ-ω) (μ+ω) f''[x]));

f[x_]=b[i]*(x^i); (* each term in the power series expansion has this form *)

eq[x_,i_]=Collect[eqrad2,{b[i],x},FullSimplify] (* We substitute the expression above into eqrad2 *)

tab=Table[(x^i),{i,n-10,n+10}];

series=Collect[Sum[eq[x,i],{i,n-3,n+3}],tab] (* Because of the linearity of the equation, instead of summing the power series terms directly in the expression above for f[x], we can sum them now: *)

(* From the results above we find out that the the recurrence relation includes three coefficients - we can extract them from the (x^(n ))power of the expression above *)

rec=SeriesCoefficient[(x^(-n))*series,{x,0,0}]

col=Collect[rec,{b[n-1],b[n],b[n+1]},FullSimplify] (*Here we collect the three terms in the recurrence relation*)

(* Denoting the reccurence relation by α[n]b[n+1]+β[n]b[n]+γ[n]b[n-1]== 0, we have (from the expression above): *)

γ[n_]=col/.{b[n-1]->1,b[n]->0,b[n+1]->0}

β[n_]=col/.{b[n-1]->0,b[n]->1,b[n+1]->0}

α[n_]=col/.{b[n-1]->0,b[n]->0,b[n+1]->1}

Nollert’s improvement

$Assumptions={};

rest[n_]:=((γ[n])/(M[n-1]))+β[n]+α[n]*M[n]; (*M[n] = a[n+1]/a[n]  = rest in cont fraction *)

M[n_]:=F0+(F1/(Sqrt[n]))+(F2/n) (* + ... *);

Series[rest[n],{n,∞,-2}]//Simplify

(* Equating the first coefficient to zero, we find F0=1 *)

F0=1;

Series[rest[n],{n,∞,-1}]//Simplify

(* We now find F1: *)

Solve[ ((F1^2)+2 (rm-rp) (Sqrt[(μ^2)-(ω^2)])) ==0,F1]//FullSimplify

(* There are two possibilites, but only one of them leads to a convergent series solution *)

QNMs for massless scalar fields in a Kerr black hole

This section calculates the quasinormal modes of massive scalar fields around a Kerr black hole.
The recurrence relations were determined in previous sections. Here we use the coefficients in the recurrence relation to study quasinormal modes.

Main references: Proc.Roy.Soc.Lond.A402,285 (1985);  Phys.Rev.D 47,5253;  https://arxiv.org/pdf/gr-qc/0511111.pdf;  https://arxiv.org/pdf/1501.06570.pdf

ClearAll["Global`*"]

QNMs for massless scalar fields in a Kerr black hole

Continued fraction equation for the angular eigenvalues

(* We write below the coefficients in the recurrence relation α[n]a[n+1]+β[n]a[n]+γ[n]a[n-1] *)

α1[n_]:=-4 (1+n) (1+n+Abs[m]);

β1[n_]:=-2 A+(m^2)+2 n (1+n)-4 a ω-4 a (2 n) ω-2 (a^2) (ω^2)+(1+2 n-4 a ω) Abs[m]+(1+2 n) Abs[m]+Abs[(m^2)];

γ1[n_]:=2 a ω (2 (n)+2*Abs[m]);

Nol1[n_]=0;(* Nollert's improvement term - setting it to zero means that the improvement is not used *)

(* Continued fraction equation: *)

EQAng[Ninv_?NumericQ,Nterms_?NumericQ]:=Block[{Rn1,Rn2},Rn1=Nol1[Nterms];Rn2=0;Do[Rn1=((-γ1[n])/(β1[n]+α1[n]*Rn1)),{n,Nterms,1+Ninv,-1}];Do[Rn2=((-γ1[n])/(β1[n-1]+α1[n-2]*Rn2)),{n,1,Ninv}];β1[Ninv]+Rn2*α1[Ninv-1]+Rn1*α1[Ninv]];

Continues fraction equation for the radial eigenvalues

rp=M+(Sqrt[(M^2)-(a^2)])(*location of the event horizon*);rm=M-(Sqrt[(M^2)-(a^2)])(* location of the Cauchy horizon *);

(* We write below the coefficients in the recurrence relation α[n]a[n+1]+β[n]a[n]+γ[n]a[n-1] *)

α[n_?NumericQ]:=(1+n) (-2 I a m+(rm-rp) (1+n)+2 I rp (rm+rp) ω);

β[n_?NumericQ]:=-4 I (1+2 n) (rp^2) ω-8 (rp^2) (rm+rp) (ω^2)+2 a m (I+2 I n+4 rp ω)+(-rm+rp) (1+A+2 n (1+n)+a ω (-2 m+a ω));

γ[n_?NumericQ]:=(n-2 I (rm+rp) ω) (-2 I a m+(rm-rp) (n)+2 I rp (rm+rp) ω);

(* Nollert's improvement term - setting it to zero means that the improvement is not used *)

Nol[n_?NumericQ]=1-(((Sqrt[2]) (Sqrt[I (rm-rp) ω]))/(Sqrt[n]))+((-(3/4)-2 I rp ω)/n);

(* Continued fraction equation: *)

EQRad[Ninv_?NumericQ,Nterms_?NumericQ]:=Block[{Rn1,Rn2},Rn1=Nol[Nterms];Rn2=0;Do[Rn1=((-γ[n])/(β[n]+α[n]*Rn1)),{n,Nterms,1+Ninv,-1}];Do[Rn2=((-γ[n])/(β[n-1]+α[n-2]*Rn2)),{n,1,Ninv}];β[Ninv]+Rn2*α[Ninv-1]+Rn1*α[Ninv]];

QNMs for Kerr

(* The procedure below determines the eigenvalues A and ω for a set of parameters a, M, m, where Nterms is the number of terms in the continued fraction, Ninv is the number of inversions, and Ainit is the initial guess for the eigenvalue. Remember that, when a*ω=0, we have A = l*(l+1) - therefore, it is usually a good idea to use l*(l+1) as the initial guess. *)

eigenv[Ainit_,ωinit_,Ninv_,Nterms_,MM_,aa_,mm_]:=Block[{$MinPrecision=20,a=aa,M=MM,m=mm,eqnang,eqnrad,A,ω},eqnang[A_?NumericQ,ω_?NumericQ]=EQAng[Ninv,Nterms];eqnrad[A_?NumericQ,ω_?NumericQ]=EQRad[Ninv,Nterms];{A,ω}/.FindRoot[{eqnang[A,ω],eqnrad[A,ω]},{{A,Ainit},{ω,ωinit}},WorkingPrecision->20]]

(*An example of calculation - the result is of the form {A,ω} *)

eigenv[2(*Ainit*),(1/10)-(1/10)*I(*ωinit*),0(*Ninv*),100(*Nterms*),1(*M*),(1/10)(*a*),1(*m*)]

(* We can write a procedure to calculate the quasinormal modes for M=1 and varying a: *)

qnmtable:=Block[{Ninv=0,Nterms=100,M=1,m=1,a=(1/1000),δa=(1/100),Ainit=2,ωinit=(1/10)-(1/10)*I,sol,listvals={}},
Monitor[While[a<M,
sol=eigenv[Ainit,ωinit,Ninv,Nterms,M,a,m];(* finds the eigenvalues *)
AppendTo[listvals,{a,sol[[1]],Re[sol[[2]]],Im[sol[[2]]]}]; (* adds the eigenvectors to the solution vector *)
{Ainit,ωinit}=sol;(* redefines the initial guesses based on the previous solutions *)
a+=δa(* Increments the value of a *)
],a];
listvals (* Returns the list of values in the form {a,A,Re[ω],Im[ω]} *)
];

tab=qnmtable;

(* We plot below the real part of the quasinormal mode as a function of the spin a *)

ListLinePlot[tab[[All,{1,3}]],PlotTheme->"Detailed",FrameLabel->{a,Re[ω]}]

(* We plot below the imaginary part of the quasinormal mode as a function of the spin a *)

ListLinePlot[tab[[All,{1,4}]],PlotTheme->"Detailed",FrameLabel->{a,Im[ω]}]

(* We plot below the imaginary part of the quasinormal mode as a function of the real part *)

ListLinePlot[tab[[All,{3,4}]],PlotTheme->"Detailed",FrameLabel->{Re[ω],Im[ω]},PlotRange->{{0.28,0.5},{-0.1,-0.03}}]

(* The same plots as before, using now the MaTeX package for LaTeX commands - https://github.com/szhorvat/MaTeX *)

<<MaTeX`

texStyle={FontFamily->"Latin Modern Roman",FontSize->18,GrayLevel[0]};

fig1=ListLinePlot[{tab[[All,{1,3}]]},PlotStyle->Thickness[0.0075],PlotRange->Automatic,PlotTheme->{"Detailed"},AspectRatio->1,ImageSize->400,FrameLabel->{MaTeX["a / M",FontSize->20],MaTeX["\\mathrm{Re}\\left(M \\omega \\right)",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted]]

fig2=ListLinePlot[{tab[[All,{1,4}]]},PlotStyle->Thickness[0.0075],PlotRange->Automatic,PlotTheme->{"Detailed"},AspectRatio->1,ImageSize->400,FrameLabel->{MaTeX["a / M",FontSize->20],MaTeX["\\mathrm{Im}\\left(M \\omega \\right)",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted]]

fig3=ListLinePlot[{tab[[All,{3,4}]]},PlotStyle->Thickness[0.0075],PlotRange->Automatic,PlotTheme->{"Detailed"},AspectRatio->1,ImageSize->400,FrameLabel->{MaTeX["\\mathrm{Re}\\left(M \\omega \\right)",FontSize->20],MaTeX["\\mathrm{Im}\\left(M \\omega \\right)",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted]]

fig=GraphicsRow[{fig1,fig2},ImageSize->800]

QBSs for scalar fields in a Kerr black hole

This section calculates the quasibound states of massive scalar fields around a Kerr black hole.
The recurrence relations were determined in previous sections. Here we use the coefficients in the recurrence relation to study quasibound states.

ClearAll["Global`*"]

QBSs for scalar fields in a Kerr black hole

Leaver’s method implementation

k=(Sqrt[(μ^2)-(ω^2)]);

Alm=(k^2)*(a^2)+SpheroidalEigenvalue[l,m,k*a];

rp=M+(Sqrt[(M^2)-(a^2)]);rm=M-(Sqrt[(M^2)-(a^2)]);

Ω=(a/(2*M*rp));

λ=Alm+(a^2)*(ω^2)-2*a*m*ω;  (* Convention in 1312.2326.pdf *)

(* Denoting the reccurence relation by α[n]a[n+1]+β[n]a[n]+γ[n]a[n-1]== 0, we have: *)

γ[n_?NumericQ]=(4 (n^2) (rm-rp) (μ-ω) (μ+ω)+4 n (rm+rp) (rm (Sqrt[(μ-ω) (μ+ω)]) ((μ^2)-2 (ω^2))+rp (2 (ω^2) (-I ω+(Sqrt[(μ-ω) (μ+ω)]))-(μ^2) (-2 I ω+(Sqrt[(μ-ω) (μ+ω)]))))-4 I a m (2 n (μ-ω) (μ+ω)+(rm+rp) (-2 (ω^2) (-I ω+(Sqrt[(μ-ω) (μ+ω)]))+(μ^2) (-2 I ω+(Sqrt[(μ-ω) (μ+ω)]))))+(((rm+rp))^2) (rm (μ^4)-rp ((μ^4)+8 (ω^3) (ω+I (Sqrt[(μ-ω) (μ+ω)]))-4 I (μ^2) ω (-2 I ω+(Sqrt[(μ-ω) (μ+ω)])))))  ;

β[n_?NumericQ]=2 ((rm^2) (μ^2) (Sqrt[(μ-ω) (μ+ω)]) (1+2 n+2 I rp ω)+2 a m (-2 (ω^2) (I+2 I n+2 rp ω+2 I rp (Sqrt[(μ-ω) (μ+ω)]))+(μ^2) (2 I+4 I n+4 rp ω-I rm (Sqrt[(μ-ω) (μ+ω)])+3 I rp (Sqrt[(μ-ω) (μ+ω)])))+2 rm (-2 rp (μ^2) (Sqrt[(μ-ω) (μ+ω)])-4 n rp (μ^2) (Sqrt[(μ-ω) (μ+ω)])-2 I (rp^2) (μ^2) ω (Sqrt[(μ-ω) (μ+ω)])+2 rp (ω^2) (Sqrt[(μ-ω) (μ+ω)])+4 n rp (ω^2) (Sqrt[(μ-ω) (μ+ω)])+4 I (rp^2) (ω^3) (Sqrt[(μ-ω) (μ+ω)])-(μ-ω) (μ+ω) (1+2 (n^2)+λ+n (2+4 I rp ω)+rp (rp (μ^2)+2 I ω+4 rp (ω^2))))+rp (2 (rp^2) (μ^4)+2 (ω^2) (-1-2 (n^2)-λ+2 rp (I+2 rp ω) (ω+I (Sqrt[(μ-ω) (μ+ω)]))+n (-2+4 I rp ω-4 rp (Sqrt[(μ-ω) (μ+ω)])))+(μ^2) (2+4 (n^2)+2 λ+n (4-8 I rp ω+6 rp (Sqrt[(μ-ω) (μ+ω)]))+rp (3 (Sqrt[(μ-ω) (μ+ω)])-2 ω (2 I+5 rp ω+3 I rp (Sqrt[(μ-ω) (μ+ω)]))))));

α[n_?NumericQ]=4 (1+n) (μ-ω) (μ+ω) (-2 I a m+(1+n) (rm-rp)+2 I rp (rm+rp) ω);

(* Below we have (Subscript[a, n+1])/(Subscript[a, n]) for Nollert's improvement: *)

Nol[n_?NumericQ]=1-(((Sqrt[2]) (Sqrt[(rm-rp) ((((μ-ω) (μ+ω)))^(3/2))]))/((Sqrt[n]) (Sqrt[-(μ^2)+(ω^2)])));

(* Below we have the continued fraction equation - Eq. 14 in Leaver, Proc.Roy.Soc.Lond.A402,285 (1985) with Nollerts improvement: Nollert
Phys.Rev.D 47,5253  *)
(* Ninv is the number of inversions and Nterms is the total number of terms *)

EQI[Ninv_?NumericQ,Nterms_?NumericQ]:=Block[{Rn1,Rn2},Rn1=Nol[Nterms];Rn2=0;Do[Rn1=((-γ[n])/(β[n]+α[n]*Rn1)),{n,Nterms,1+Ninv,-1}];Do[Rn2=((-γ[n])/(β[n-1]+α[n-2]*Rn2)),{n,1,Ninv}];β[Ninv]+Rn2*α[Ninv-1]+Rn1*α[Ninv]];

(* The procedure below finds a root of the equation above using ωinit as the initial guess; MM is the mass of th bh, aa is the angular momentum of the bh, (ll,mm) are the orbital and azimuthal numbers of the scalar field, and μμ is the mass of the scalar field  *)

insteq[ωinit_?NumericQ,Ninv_?NumericQ,Nterms_?NumericQ,MM_?NumericQ,aa_?NumericQ,ll_?NumericQ,mm_?NumericQ,μμ_?NumericQ]:=Block[{$MinPrecision=50,M=MM,a=aa,l=ll,m=mm,μ=μμ,eqn},eqn[ω_?NumericQ]=EQI[Ninv,Nterms];ω/.FindRoot[{eqn[ω]==0},{ω,ωinit},WorkingPrecision->15]]

(* Let's test the procedure defined above. For that, we shall use Dolan's results in https://arxiv.org/pdf/0705.2880.pdf *)

(* First test: table I on page 18 *)

(*a = 0.7, μ = 0.187, Im[ω]=3.33*(10^(-10)) *) μt=(187/1000);insteq[(186/1000)+2*(10^(-8))*I,0,100,1,7/10,1,1,μt]

(*a = 0.95, μ = 0.343, Im[ω]=4.88*(10^(-8))*) μt=(343/1000);insteq[(337/1000)+5*(10^(-6))*I,0,100,1,95/100,1,1,μt]

(* The procedure seems to be very sensitive to the real part of the initial guess (at least in comparison with the imaginary part). For instance, if we change sligthly the real part, we get an overtone instead of the fundamental unstable mode: *)

(*a = 0.95, μ = 0.343, Im[ω]=4.88*(10^(-8))*) μt=(343/1000);insteq[(340/1000)+5*(10^(-6))*I,0,100,1,95/100,1,1,μt]

Finding the maximum instability (i.e. the value of μ for which the instability is maximized)

(* The procedure below finds the instability as a function of μ for M=1 and chosen a, l, m. The parameter δμ is the resolution in μ space, Nin is the number of inversions in the continued fraction and Nte is the number of terms. Data is stored in a list of the form {μ,Re[ω],Im[ω]}  *)

maxinst[aa_?NumericQ,ll_?NumericQ,mm_?NumericQ]:=(Off[Interpolation::inhr,InterpolatingFunction::dmval];Block[{ωt,μv,δμ=(1/200),Nin=0,Nte=200,μfin,int=0,numint,tim=AbsoluteTime[],Ωt=Ω/.{M->1,a->aa},sol,list={},gre,gim},μv=(25/100)(* initial scalar mass *);μfin=(45/100)(* final scalar mass *);ωt=(98/100)*μv+2*(10^(-8))*I (* initial guess - it is very important - by changing this we may find different overtones or the method may fail (if the inital guess is far off) *);numint=((μfin-μv)/(δμ))(* estimate for number of iterations *);PrintTemporary["Expected iterations = ",numint];PrintTemporary["See below {time elapsed (sec),current mass,current iteration}:"];(* While loop that uses the previous mode found as an initial guess for the next "FindRoot". Actually, it interpolates the previous modes found to predict the next one *)Monitor[While[μv<μfin,sol=insteq[ωt,Nin,Nte,1,aa,ll,mm,μv];AppendTo[list,{μv,Re[sol],Im[sol]}];gre=Interpolation[list[[All,{1,2}]]];gim=Interpolation[list[[All,{1,3}]]];
μv=μv+δμ;
ωt=gre[μv]+I*gim[μv];
int=int+1],{AbsoluteTime[]-tim,μv,int+1}];Export["data1.dat",list](*name of the file where you want to save the data*);Print["Total Time = ",AbsoluteTime[]-tim," sec, Number of iterations = ",int]];On[Interpolation::inhr,InterpolatingFunction::dmval]);

maxinst[95/100,1,1]

data1=Import["data1.dat"]//ToExpression;

ListPlot[data1[[All,{1,3}]],AxesLabel->{"μ","Im(Mω)"}]

(* close-up look at the plot above *)

ListPlot[data1[[All,{1,3}]],PlotRange->{-(10^(-7)),(10^(-7))},AxesLabel->{"μ","Re(Mω)"}]

ListPlot[data1[[All,{1,2}]],AxesLabel->{"μ","Re(Mω)"}]

(* The same plots as before, using now the MaTeX package for LaTeX commands - https://github.com/szhorvat/MaTeX *)

<<MaTeX`

texStyle={FontFamily->"Latin Modern Roman",FontSize->18,GrayLevel[0]};

fig1=ListLinePlot[data1[[All,{1,3}]],PlotStyle->Thickness[0.007],PlotRange->{{0.25,0.4},{-(10^(-7)),(10^(-7))}},PlotTheme->{"Detailed"},AspectRatio->1,ImageSize->400,FrameLabel->{MaTeX["\\mu M",FontSize->20],MaTeX["\\mathrm{Im}\\left( M \\omega \\right)",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["(a,\\ell,m)=(0.95,1,1)",FontSize->16]},{{0.4,0.89}}]]

fig2=ListLinePlot[data1[[All,{1,2}]],PlotStyle->Thickness[0.007],PlotRange->{{0.25,0.4},{0.25,0.39}},PlotTheme->{"Detailed"},AspectRatio->1,ImageSize->400,FrameLabel->{MaTeX["\\mu M",FontSize->20],MaTeX["\\mathrm{Re}\\left( M \\omega \\right)",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["(a,\\ell,m)=(0.95,1,1)",FontSize->16]},{{0.4,0.89}}]]

fig=GraphicsRow[{fig1,fig2},ImageSize->850]

Ω=(0.95/(2*1*(1+(Sqrt[1-(0.95^2)]))));f[x_,y_,z_]:={x,(y/x),(y/(1*Ω))};(* function to produce new data set of the form {μ,Re(Mω)/μ,Re(Mω)/mΩ} *)

data2=f@@@data1; (*applies the function to the data *)

figb1=ListLinePlot[data2[[All,{1,2}]],PlotStyle->Thickness[0.007],PlotRange->Automatic,PlotTheme->{"Detailed"},AspectRatio->1,ImageSize->400,FrameLabel->{MaTeX["\\mu M",FontSize->20],MaTeX["\\mathrm{Re}\\left(\\omega \\right) / \\mu ",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["(a,\\ell,m)=(0.95,1,1)",FontSize->16]},{{0.7,0.82}}]]

figb2=ListLinePlot[data2[[All,{1,3}]],PlotStyle->Thickness[0.007],PlotRange->Automatic,PlotTheme->{"Detailed"},AspectRatio->1,ImageSize->400,FrameLabel->{MaTeX["\\mu M",FontSize->20],MaTeX["\\mathrm{Re}\\left(\\omega \\right) / (m \\Omega) ",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["(a,\\ell,m)=(0.95,1,1)",FontSize->16]},{{0.35,0.82}}]]

figb=GraphicsRow[{figb1,figb2},ImageSize->850]

Scattering problem in the time domain

Quit[]

Scattering problem in the time domain - massless scalar field around Schwarzschild

In this section we solve the scattering problem of a massless scalar field around a Schwarzschild black hole in the time domain. We follow the method described in Phys.Rev. D49 (1994) 883-889 [arXiv:gr-qc/9307009], by Gundlach, Price and Pullin.

(* Following Gundlach-Price-Pullin, we denote u = t - (SubStar[r]) as the retarded time and v = t + (Subscript[r,  *]) as the advanced time. We create an u-v grid, such that v[i]=i*dv e u[j]=j*du. Therefore, we identify ψ[v,u]=ψ[i*dv,j*du] as ψ[i,j]. *)

(* We first define the parameters: *)

M=1/2(*mass of the black hole*);rh=2*M(*location of the event horizon*);l=1(*orbital quantum number*);umax=200(*max value of the u coordinate*);umin=0(*min value of the u coordinate*);vmax=200(*max value of the v coordinate*);vmin=0(*min value of the v coordinate*);
nu=999;nv=999; (* number of points in the grid *);
du:=((umax-umin)/(nu+1))(*step size for u *); dv:=((umax-umin)/(nv+1)) (*step size for v *);i0=0 (* initial i value *);j0=0 (* initial j value *);

Inverting r* in terms of r:

(* Definition of the tortoise coordinate: *)

rs[r_]:=r+2*M*Log[(r-2*M)/(2*M)];

Plot[rs[r],{r,rh,5*rh}]

(*Inversion of the tortoise coordinate: y = r(rs) *)

y=InverseFunction[rs];

y[rs]

Plot[y[x],{x,-10,10}]

(* Effective potential for the scalar field (in terms of the usual coordinate r): *)

V0[r_]:=(1-((2*M)/r))*(((l*(l+1))/(r^2))+((2*M)/(r^3)));

Plot[V0[r],{r,rh,5*rh}]

(* Effective potential for the scalar field (in terms of the tortoise coordinate rs): *)

V[rs_]:=V0[y[rs]];

Plot[V[rs],{rs,-10,10}]

Discretization of the wave equation -4(Subscript[Ψ, ,uv])=  (Subscript[V, l])(r)*Ψ  and solution of the problem assuming a gaussian profile for the initial wave.

(*To compute the tortoise coordinate as a function of u and v, we use r*=(1/4)*((Subscript[v, N])+(Subscript[v, W])-(Subscript[u, N])-(Subscript[u, E]))=(1/4)*((2*i+1)*dv-(2*j+1)*du) - see the paper  *)

disc[i_,j_]:=ψ[i+1,j]+ψ[i,j+1]-ψ[i,j]-((du*dv)/8)*V[(1/4)*((2*i+1)*dv-(2*j+1)*du)]*(ψ[i,j+1]+ψ[i+1,j]);

(* The wave equation is ψ[i+1,j+1]=disc[i,j] *)

(* Gaussian to be used as the initial condition *) gaus[x_]:=a*Exp[((-(((x-b))^2))/(2*(c^2)))];a=1/100;b=30;c=3;

Plot[gaus[v],{v,vmin,vmax},PlotRange->{{0,100},All}]

(*initial condition: gaussian at u=u0 and constant at v=v0 *)( Do[ψ[i,j0]=gaus[i*dv],{i,i0,nv}];Do[ψ[i0,j]=gaus[0],{j,j0,nu}])

Monitor[Do[Do[ψ[i+1,j+1]=N[disc[i,j],1000],{i,i0,nv-1}],{j,j0,nu-1}],{i,j}](* This is where we actually sole the equation. It can take a while if nv and nu are large numbers. *)

Extraction of the wave function ψ[t,r*] for fixed r*.

(* We would like to: (1) transform ψ[v,u] into ψ[t,r*], (2) fix r*=rext (extraction radius), and (3) obtain ψ[t,rext]. After that, we can plot ψ[t,rext], from tmin to tmax. Since the grid has a finite size, we need to determine what values we should use for tmin and tmax *)

rext=30;tmin=rext;tmax=Min[vmax-rext,umax+rext];

tab=Flatten[Table[{i*dv,j*du,ψ[i,j]},{i,i0,nv},{j,j0,nu}],1] ;(* produces a table of values *)

g=Interpolation[tab];(* interpolates the data to produce a function psi in terms of v and u *)

psi[t_,rs_]:=g[t+rs,t-rs];(* psi in terms of t and r* *)

Plot[psi[t,rext],{t,tmin+1,tmax-1},PlotRange->All]

(* The oscillations that decay with time are the quadinormal modes. The "tail" that appears after the oscillations is the so-called "late-time tail": *)

LogPlot[{Abs[psi[t,rext]]},{t,tmin,tmax},PlotRange->All]

Animate[Plot[psi[t,r],{r,Max[-t,-199+t],Min[t,199-t]},PlotRange->{{-99,99},{-0.01,0.01}},PlotLabel->"t/2M="<>ToString[t]],{t,1,98,1}]

<<MaTeX`

texStyle={FontFamily->"Latin Modern Roman",FontSize->18,GrayLevel[0]};

t=20;fig1=Plot[psi[t,r],{r,Max[-t,-199+t],Min[t,199-t]},PlotRange->{{-70,70},{-0.01,0.01}},PlotStyle->Thickness[0.01],PlotTheme->"Detailed",AspectRatio->1,FrameLabel->{MaTeX["r_* / 2M",FontSize->20],MaTeX[" \\psi(t,r_*) ",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["t/2M=20",FontSize->20]},{{0.5,0.1}}]]

t=70;fig2=Plot[psi[t,r],{r,Max[-t,-199+t],Min[t,199-t]},PlotRange->{{-70,70},{-0.01,0.01}},PlotStyle->Thickness[0.01],PlotTheme->"Detailed",AspectRatio->1,FrameLabel->{MaTeX["r_* / 2M",FontSize->20],MaTeX[" \\psi(t,r_*) ",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["t/2M=70",FontSize->20]},{{0.5,0.1}}]]

fig=GraphicsRow[{fig1,fig2},ImageSize->800]

figb=LogPlot[{Abs[psi[t,rext]]},{t,tmin,tmax},PlotRange->All,PlotStyle->Thickness[0.007],PlotTheme->"Detailed",AspectRatio->1/2,ImageSize->600,FrameLabel->{MaTeX["t / 2M",FontSize->20],MaTeX[" \\left| \\psi(t,r_*) \\right| ",FontSize->20]},LabelStyle->texStyle,Mesh->None,MeshStyle->Transparent,GridLines->Automatic,GridLinesStyle->Directive[GrayLevel[0.5], Dotted],PlotLegends->Placed[{MaTeX["r_*/2M=30",FontSize->20]},{{0.5,0.1}}]]

SetDirectory[NotebookDirectory[]];

Export["scat_time.pdf",fig]

Export["scat_time2.pdf",figb]