BHPToolkit Workshop

Leanne Durkan, PhD Student, UCD
Supervisors: Dr. Niels Warburton and Prof. Adrian Ottewill

Package Contributors: B. Wardell, N. Warburton, M. Casals, A. Ottewill, C. Kavanagh, D. Aruquipa, L. Durkan

MST and Numerical Integration Methods for the Regge-Wheeler Package

Installation

Add the Black Hole Perturbation Toolkit paclet server to Mathematica’s list of paclet servers:

PacletSiteAdd["http://pacletserver.bhptoolkit.org","Black Hole Perturbation Toolkit Paclet Server"]

Install the ReggeWheeler package:

PacletInstall["ReggeWheeler"]

Get an updated list of packages available on the server:

PacletSiteUpdate["http://pacletserver.bhptoolkit.org"]

To see which paclets are available:

PacletFindRemote["*","Location"->"http://pacletserver.bhptoolkit.org"]

Finish installing the ReggeWheeler package :

<<ReggeWheeler`

Regge-Wheeler Package

See documentation for further details and tutorials.

Function Usage

?ReggeWheelerRadial

Defining relevant parameters

(*Convention: M=1*)
s0=2;
l0=2;
ω0=0.1`32;

Standard output (default MST):

R=ReggeWheelerRadial[s0,l0,ω0]

Ingoing solution (MST):

R["In"]

Ingoing solution at a specific radius (MST):

r0=10;

R["In"][r0]

We can also compute derivatives in the normal way for functions (derivative of the MST series, not a numerical derivative of previous result):

R["In"]'[r0]

We can also evaluate on a grid of points:

R["In"][{10,11,12}]

Numerical integration option for ingoing solution over a given domain (options should be strings):

ReggeWheelerRadial[s0,l0,ω0, Method-> {"NumericalIntegration", "Domain"-> {10,100}}, "BoundaryConditions"-> "In"]

Numerical integration option at a specific radius:

ReggeWheelerRadial[s0,l0,ω0, Method-> {"NumericalIntegration", "Domain"-> {10,100}}, "BoundaryConditions"-> "In"][r0]

Numerical integration option for ingoing and upgoing solution over a given domain:

ReggeWheelerRadial[s0,l0,ω0, Method-> {"NumericalIntegration", "Domain"-> {"Up"->{10,100},"In"-> {3,10}}}]

Numerical method:   faster, lower precision
MST method: slower, higher precision

Full set of properties:

R["In"]["s"]

R["In"]["l"]

R["In"]["ω"]

R["In"]["Eigenvalue"]

R["In"]["Method"]

R["In"]["BoundaryConditions"]

R["In"]["Amplitudes"]

R["In"]["Domain"]

Renormalized Angular Momentum:

?RenormalizedAngularMomentum

m0=1;
(*Ω=((((M/(r0^3))))^(1/2));
ω0=N[m0 Ω,22];*)
a0=0; (*Regge-Wheeler -> Schwarzschild*)
λ0=l0(l0+1)-s0(s0+1);
ν0=RenormalizedAngularMomentum[s0,l0,m0,a0,ω0,λ0]

R["In"]["Method"]

Further Options

Options[ReggeWheelerRadial]

Example:

ReggeWheelerRadial[s0,l0,ω0,PrecisionGoal->8,AccuracyGoal->Infinity]["In"][10]

Comparison with Regge Wheeler Numerical Integration Method with MST: Boundary Conditions

See documentation: Boundary Conditions

Asymptotic Boundary Conditions

MST expression for (Subscript[ψ, in]) valid over the domain      (SubPlus[r]) ≤ r < ∞     →     Can obtain  (Subscript[ψ, in]) on ℋ
MST expression for (Subscript[ψ, in]) on ℐ pending, but is defined for Teukolsky case.
MST expression for (Subscript[ψ, up]) valid over the domain    (SubPlus[r]) < r ≤ ∞      →     Can obtain (Subscript[ψ, up]) on ℐ
MST expression for (Subscript[ψ, up]) on ℋ pending, but is defined for Teukolsky case.

Upgoing Case: Near ℐ the behaviour should be (Subscript[ψ, up])≃(e^(i ω (SubStar[r])))

Can leave in terms of r as similar behaviour to (SuperStar[r]) near ℐ

{ψIn,ψUp}=Values[ReggeWheelerRadial[s0,l0,ω0]];

ψUpℐ=Table[Flatten[{r,ReIm[ψUp[r]]}],{r,10000,10100,5}];

ψUpℐN=ReggeWheelerRadial[s0,l0,ω0,"BoundaryConditions"->"Up",Method->{"NumericalIntegration","Domain"->{10000,10100}}];

MST and Numerical methods agree near ℐ:

Show[ListPlot[{ψUpℐ[[;;,{1,2}]],ψUpℐ[[;;,{1,3}]]},PlotTheme->"Detailed"],Plot[Evaluate[ReIm[ψUpℐN[r]]],{r,10000,10100}]]

MST and Numerical methods agree with asymptotic behaviour on ℐ:

M=1;
rstar[r_]:=r+2M Log[(r/(2M))-1];

Show[ListPlot[{ψUpℐ[[;;,{1,2}]],ψUpℐ[[;;,{1,3}]]},PlotTheme->"Detailed"],Plot[Evaluate[ReIm[ψUpℐN[r]]],{r,10000,10100}],Plot[Evaluate[ReIm[Exp[I ω0 rstar[r]]]],{r,10000,10100},PlotStyle->{Directive[Darker[ColorData[97,1]],Dotted],Directive[Darker[ColorData[97,2]],Dotted]}]]

Ingoing Case: Near ℋ the behaviour should be (Subscript[ψ, in])≃(e^(-i ω (SubStar[r]))) (Exercise)

Need to use (SuperStar[r]) near ℋ. Obtain r in terms of tortoise coordinate:

Solve[rstar==r+2M Log[(r/(2M))-1],r]

rt[rstar_]:=2M (1+ProductLog[(Sqrt[(E^(-2M+rstar))])]);

1.) Do the MST and Numerical methods agree near ℋ?

ψInℋ=Table[Flatten[{rstar,ReIm[ψIn[rt[rstar]]]}],{rstar,-100,0,5}];

ψInℋN=ReggeWheelerRadial[s0,l0,ω0,"BoundaryConditions"->"In",Method->{"NumericalIntegration","Domain"->{rt[-50],rt[0]}}];

MST and Numerical methods agree near ℋ somewhat:

Show[ListPlot[{ψInℋ[[;;,{1,2}]],ψInℋ[[;;,{1,3}]]},PlotTheme->"Detailed"],Plot[Evaluate[ReIm[ψInℋN[rt[rstar]]]],{rstar,-50,0}],PlotRange->All]

2.) Do the MST and Numerical methods agree with asymptotic behaviour on ℋ?

Near ℋ the behaviour should be (Subscript[ψ, in])≃(e^(-i ω (SubStar[r]))). This is satisfied by the MST evaluation, but not currently by numerical integration:

Show[ListPlot[{ψInℋ[[;;,{1,2}]],ψInℋ[[;;,{1,3}]]},PlotTheme->"Detailed"],Plot[Evaluate[ReIm[ψInℋN[rt[rstar]]]],{rstar,-50,0}],Plot[Evaluate[ReIm[Exp[-I ω0 rstar]]],{rstar,-100,0},PlotStyle->{Directive[Darker[ColorData[97,1]],Dotted],Directive[Darker[ColorData[97,2]],Dotted]}],PlotRange->All]

Comparison with Regge Wheeler Numerical Integration Method with MST:  Field Equations

See documentation: Field Equations

Regge-Wheeler equation

RW[ψ_,s_,l_,r_]:=(1-2M/r)ψ''[r]+((2M)/(r^2))ψ'[r]+((((1-2M/r))^(-1))((ω0)^2)-(((l(l+1))/(r^2))+((2M(1-(s^2)))/(r^3))))ψ[r];

{ψInMST,ψUpMST}=Values[ReggeWheelerRadial[s0,l0,ω0]];

Ingoing MST

Showing that the MST solution to (Subscript[ψ, in]) satisfies the Regge-Wheeler equation:

ψInMSTi=Interpolation[Table[{r,ψInMST[r]},{r,8,22,0.1`32}],InterpolationOrder->8];

RWψInMSTi=Monitor[Table[Flatten[{r,ReIm[RW[ψInMSTi,s0,l0,r]]}],{r,10,20,0.01}],r];

RWψInMST=Monitor[Table[Flatten[{r,ReIm[RW[ψInMST,s0,l0,r]]}],{r,10,20,0.1}],r];

ListPlot[{RWψInMST[[;;,{1,2}]],RWψInMST[[;;,{1,3}]]},PlotTheme->"Detailed",PlotRange-> All]

ListPlot[{RWψInMSTi[[;;,{1,2}]],RWψInMSTi[[;;,{1,3}]]},PlotTheme->"Detailed",PlotRange-> All]

Ingoing Numerical (Exercise)

Does the numerical integration solution to (Subscript[ψ, in]) satisfy the Regge-Wheeler equation?

ψInNI=ReggeWheelerRadial[s0,l0,ω0,"BoundaryConditions"->"In",Method->{"NumericalIntegration","Domain"->{10,20}}];

RWψInNI=Monitor[Table[Flatten[{r,ReIm[RW[ψInNI,s0,l0,r]]}],{r,10,20,0.01}],r];

ListPlot[{RWψInNI[[;;,{1,2}]],RWψInNI[[;;,{1,3}]]},PlotTheme->"Detailed",PlotRange-> All]

Upgoing MST (Exercise)

Does the MST solution to (Subscript[ψ, up]) satisfy the Regge-Wheeler equation?

ψUpMSTi=Interpolation[Table[{r,ψUpMST[r]},{r,8,22,0.1`32}],InterpolationOrder->8];

RWψUpMSTi=Monitor[Table[Flatten[{r,ReIm[RW[ψUpMSTi,s0,l0,r]]}],{r,10,20,0.01}],r];

RWψUpMST=Monitor[Table[Flatten[{r,ReIm[RW[ψUpMST,s0,l0,r]]}],{r,10,20,0.1}],r];

ListPlot[{RWψUpMSTi[[;;,{1,2}]],RWψUpMSTi[[;;,{1,3}]]},PlotTheme->"Detailed",PlotRange-> All]

ListPlot[{RWψUpMST[[;;,{1,2}]],RWψUpMST[[;;,{1,3}]]},PlotTheme->"Detailed",PlotRange-> All]

Upgoing Numerical (Exercise)

Does the numerical integration solution to (Subscript[ψ, up]) satisfy the Regge-Wheeler equation?

ψUpNI=ReggeWheelerRadial[s0,l0,ω0,"BoundaryConditions"->"Up",Method->{"NumericalIntegration","Domain"->{10,20}}];

RWψUpNI=Monitor[Table[Flatten[{r,ReIm[RW[ψUpNI,s0,l0,r]]}],{r,10,20,0.01}],r];

ListPlot[{RWψUpNI[[;;,{1,2}]],RWψUpNI[[;;,{1,3}]]},PlotTheme->"Detailed",PlotRange-> All]

Exercise: Calculating the Retarded Solution for a Circular Orbit

M=1;
r0=10;
Ω=((((M/(r0^3))))^(1/2));
s0=2;
l0=2;
m0=2;
q0=1;
ω0=N[Ω,22];

{ψInMST,ψUpMST}=Values[ReggeWheelerRadial[s0,l0,ω0]];

Define Wronskian

Clear[W]
W[l_,r1_]:=W[l,r1]=((1-2M/r1)ψInMST'[r0]ψUpMST[r1]-(1-2M/r1)ψUpMST'[r0] ψInMST[r1])/.r->r1;

Define matching

Clear[cup]
cup[l_,m_,r0_]:=cup[l,m,r0]=-ψInMST[r0]4 π q0 Conjugate[SphericalHarmonicY[l,m,π/2,0]]/(W[l,r0]r0^2Sqrt[1-3M/r0])

Clear[cin]
cin[l_,m_,r0_]:=cin[l,m,r0]=-ψUpMST[r0]4 π q0 Conjugate[SphericalHarmonicY[l,m,π/2,0]]/(W[l,r0]r0^2Sqrt[1-3M/r0])

Define retarded solution

ψ0ret[l_,m_,r_,r0_]:=If[r≤r0,cin[l,m,r0]ψInMST[r],cup[l,m,r0]ψUpMST[r]]

4.) Plot retarded solution:

ψ0retdata=Monitor[Table[Flatten[{r,ReIm[ψ0ret[l0,m0,r,r0]]}],{r,3,200,1}],r];

ListLogLinearPlot[{ψ0retdata[[;;,{1,2}]],ψ0retdata[[;;,{1,3}]]},PlotRange->All,Joined->True,PlotTheme->"Detailed",ImageSize->800,BaseStyle->20,FrameLabel->{"r/M","(Subscript[ψ, 0])"}]

Further Exercises

1). See ReggeWheeler documentation for examples on calculating fluxes of energy and angular momentum.
2). Try the above examples using the TeukolskyRadialFunction instead from the Teukolsky package.