Getting Started with the Black Hole Perturbation Toolkit

In this tutorial, we will give an overview of the Mathematica components of the Black Hole Perturbation Toolkit.

Installation

We recommend that you install the Toolkit Mathematica packages via the Paclet server using the instructions below.

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

If[$VersionNumber≥12.1,
PacletSiteRegister["https://pacletserver.bhptoolkit.org","Black Hole Perturbation Toolkit Paclet Server"],
PacletSiteAdd["http://pacletserver.bhptoolkit.org","Black Hole Perturbation Toolkit Paclet Server"]
]

Get an updated list of packages available on the server:

If[$VersionNumber≥12.1,
PacletSiteUpdate["https://pacletserver.bhptoolkit.org"],
PacletSiteUpdate["http://pacletserver.bhptoolkit.org"]
]

Install the paclets:

PacletInstall["GeneralRelativityTensors"];
PacletInstall["KerrGeodesics"];
PacletInstall["SpinWeightedSpheroidalHarmonics"];
PacletInstall["ReggeWheeler"];
PacletInstall["Teukolsky"];
PacletInstall["PostNewtonianSelfForce"];

If a new version of a package is released and you would like to update, just run steps 2 and 3 above again and the latest available version will be installed.

Getting Help

Documentation for each of the packages is available from the Documentation Center, which can be found under the Help -> Wolfram Documentation menu item.

A full list of installed documentation is available by clicking the “Add-ons and packages” link at the bottom.

Alternatively, you can search for a documentation on any of the provided functions:

All of the packages provide Guide pages (which list the available functions), Tutorials (which provide a walkthrough of the functionality with worked examples) and Reference Pages (which provide descriptions of the functions).

You can also find documentation for a specific function using Mathematica’s ? command (this will only work after loading the relevant package

?Simplify

Component-based tensor algebra

The GeneralRelativityTensors package provides a set of tools for performing component-based tensor algebra. It does not aim to reproduce the all of the features of something like xAct, but rather focuses on providing a simplified tool for doing component calculations in a specific spacetime.

Load GeneralRelativityTensors package

<<GeneralRelativityTensors`

Metrics

All calculations in GeneralRelativityTensors rely on having a spacetime metric, which is defined using the ToMetric function. This can either load a built-in metric, or define a new metric.

Built in metrics

We can load the metric for the Schwarzschild spacetime in standard coordinates

gS=ToMetric["Schw"]

TensorValues[gS]

Defining a new metric

To define a new metric, we need to have its components in the chosen basis

newMetVals={{-f[x,y],0,0},{0,g[x,y](x^2),x y g[x,y]},{0,x y g[x,y],g[x,y](x^2)}};
newMetVals//MatrixForm

Then, we can define our new metric using ToMetric

gNew=ToMetric[{"NewMetric","(g^new)"},{t,x,y},newMetVals,"Latin"]

Tensors

We can define new Tensors using ToTensor. Tensors always have associated metrics, which are specified at the time that they are defined.

T=ToTensor[{"NewTensor","T"},gNew,{q[t],r[x],s[y]},{-a}]

Accessing components of a tensor

There are two ways to access the components of a tensor. Either we can ask for a specific component:

gS[-t,-t]

or we can as for an array of all components:

TensorValues[gS]

Raising and lowering indices

Raise and lowering of indices (with the metric) is automatic based on the indices (arguments passed to the tensor):

T[a]

T[a]//TensorValues

T[x]

Special tensors

As a simple example, let’s now compute some tensors in Schwarzschild spacetime.

Christoffel symbols

First we compute the Christoffel symbols,  (Subscript[(Γ^α), βγ]) associated with the metric

chrS=ChristoffelSymbol[gS]

chrS//TensorValues

Ricci tensor

Compute the Ricci tensor:

ricTens=RicciTensor[gS]

ricTens//TensorValues

We know this should be zero, but it is not obviously so here. We need to simplify the result by using ActWith->Simplify to clean up internal expressions

RicciTensor[gS,ActWith->Simplify]//TensorValues

Ricci scalar

Compute the Ricci scalar

ricScal=RicciScalar[gS,ActWith->Simplify]

ricScal//TensorValues

Merging tensors

Sometimes we want to construct a tensor from some other tensors. For this we use MergeTensors. As an example, let’s compute the Einstein tensor.

einsExpr=ricTens[-α,-β]-gS[-α,-β] ricScal/2

This is not a single tensor yet, but we construct a tensor from it:

einS=MergeTensors[einsExpr,{"EinsteinSchwarzschild","G"},ActWith->Simplify]

As expected for Schwarzschild spacetime, this is zero

einS//TensorValues

Derivatives

Partial derivatives

D is overloaded to take derivatives with respect to indices

gSDiv=D[gS[α,β],-α]

Indices are not traced until forced

gSDiv//TensorValues//Simplify

Covariant derivatives

dgS=CovariantD[gS,-γ]

MergeTensors[dgS,ActWith->Simplify]//TensorValues

Geodesics in Kerr spacetime

The KerrGeodesics package provides a set of functions for calculating properties of bound timelike geodesic orbits about a Kerr black hole.

Load KerrGeodesics package

<<KerrGeodesics`

Geodesic equations

The geodesic equations in Kerr spacetime are

(((Σ((d r)/(dτ))))^2)== -(1-(ℰ^2))(r^4)+2 (r^3) -((a^2) (1-(ℰ^2))+(ℒ^2)+𝒬) (r^2)+2((a^2) (ℰ^2)-2 a ℰ ℒ+(ℒ^2)+𝒬) r -(a^2)𝒬
(((Σ((d z)/(dτ))))^2)==(a^2) (1-(ℰ^2)) (z^4)-((a^2) (1-(ℰ^2))+(ℒ^2)+𝒬) (z^2)+𝒬
Σ((d t)/(dτ)) == ((r (2 (a^2)  ℰ-2 a  ℒ+(a^2) ℰ r+ℰ (r^3)))/Δ)+(a^2) ℰ (z^2)
Σ((d ϕ)/(dτ)) ==((a (2  ℰ r-a ℒ))/Δ)+(ℒ/(1-(z^2)))

Introducing "Mino time" dλ ==  ((dτ)/Σ), they take the alternative form

(((((d r)/(dλ))))^2)== (-1+(ℰ^2))(r^4)+2 (r^3) +((a^2) (-1+(ℰ^2))-(ℒ^2)-𝒬) (r^2)+2((a^2) (ℰ^2)-2 a ℰ ℒ+(ℒ^2)+𝒬) r -(a^2)𝒬 ==R(r)
(((((d z)/(dλ))))^2)==-(a^2) (-1+(ℰ^2)) (z^4)+((a^2) (-1+(ℰ^2))-(ℒ^2)-𝒬) (z^2)+𝒬== Z(z)
((d t)/(dλ)) == ((r (2 (a^2)  ℰ-2 a  ℒ+(a^2) ℰ r+ℰ (r^3)))/Δ)+(a^2) ℰ (z^2)
((d ϕ)/(dλ)) ==((a (2  ℰ r-a ℒ))/Δ)+(ℒ/(1-(z^2)))

Orbital parametrization: "apex" parameters

The radial equation has four turning points

R(r)== (-1+(ℰ^2))(r-(Subscript[r, 1]))(r-(Subscript[r, 2]))(r-(Subscript[r, 3]))(r-(Subscript[r, 4]))

(for bound orbits)
(Subscript[r, 1])≥ (Subscript[r, 2])≥ (Subscript[r, 3])>(SubPlus[r])>(SubMinus[r])>(Subscript[r, 4])≥ 0

Similarly the polar equation has two turning points

Z(z)==((z^2)-(Subscript[z, 1]^2))((a^2)(-1+(ℰ^2))(z^2)-(Subscript[z, 2]^2))

(for bound orbits)
0≤ (Subscript[z, 1])≤ 1≤ (Subscript[z, 2])

The turning points (Subscript[r, 1]), (Subscript[r, 2]) and (Subscript[z, 1]) uniquely identify a bound orbit. In the KerrGeodesics package they are expressed in terms of the semi-latus rectum p, eccentricity e and inclination parameter x

(Subscript[r, 1])== (p/(1-e))
(Subscript[r, 2])== (p/(1+e))
(Subscript[z, 1])==(Sqrt[1-(x^2)])

or

p==2(((Subscript[r, 1])(Subscript[r, 2]))/((Subscript[r, 1])+(Subscript[r, 2])))
e== (((Subscript[r, 1])-(Subscript[r, 2]))/((Subscript[r, 1])+(Subscript[r, 2])))
x ==  Sign[ℒ](Sqrt[1-(Subscript[z, 1]^2)])

(circular) 0<e<1 (parabolic)
(retrograde) -1 < x < 1 (prograde)
p> (Subscript[p, LSO])(e,x)

Constants of Motion

Kerr spacetime has two Killing vector fields:

(kt^μ):=-(∂/(∂t))     and     ((kϕ)^μ):=(∂/(∂ϕ))

These are associated with two conserved quantities, the (specific) energy and angular momentum

ℰ ==  (Subscript[kt, μ])(u^μ)

ℒ ==(Subscript[kϕ, μ])(u^μ)

The spacetime also admits a Killing tensor

(Subscript[K, μν])= 2Σ (Subscript[l, (μ])(Subscript[n, ν)])+(r^2)(Subscript[g, μν])

This leads to an additional conserved quantity, the Carter constant :

𝒬 +(((ℒ-a ℰ))^2) == (u^μ)(Subscript[K, μν])(u^ν)

Analytic evaluation

In some cases, the constants may be determined analytically.

KerrGeoEnergy[a,p,e,x]

Numerical Evaluation

Even if an analytic value is not possible, the constants can be determined numerically

With[{a=0.9,p=10.,e=0.5,x=0.5},
KerrGeoConstantsOfMotion[a,p,e,x]
]

Frequencies

There are three orbital frequencies, associated with the radial, polar and azimuthal motion. Depending on parametrization, these take different values. The package currently supports three options for parametrization

With[{a=0.9,p=10.0,e=0.5,x=0.5},
{KerrGeoFrequencies[a,p,e,x,"Time"-> "Mino"],
KerrGeoFrequencies[a,p,e,x,"Time"-> "Proper"],
KerrGeoFrequencies[a,p,e,x,"Time"-> "BoyerLindquist"]}
]

Orbits

Constructing an orbit

With[
{a=0.9,p=10.0,e=0.5,x=0.5},
orbit=KerrGeoOrbit[a,p,e,x]
]

Orbit properties

Extract constants of motion and frequencies

orbit["ConstantsOfMotion"]

orbit["Frequencies"]

List all possible properties

Keys@orbit[[-1]]

Orbital trajectory

We can use an orbit object to evaluate the orbital position at a specified time

orbit[2.]

Using this, we can create a plot of the orbit

Module[
{to,ro,θo,ϕo,λ,Λr},
{to,ro,θo,ϕo}=orbit[λ];
Λr=2π/orbit["RadialFrequency"];
Show[
ParametricPlot3D[
ro{Cos[ϕo]Sin[θo],Sin[ϕo]Sin[θo],Cos[θo]},{λ,0,1Λr}
,Lighting->"Neutral",PlotPoints->20,PlotRange->20],
Graphics3D[{Black,Specularity[.5],Sphere[{0,0,0},1+(Sqrt[1-((.9)^2)])]}]
]
]

Special Orbits

Innermost stable circular (equatorial) orbit (ISCO)

(Subscript[r, 1])==(Subscript[r, 2])== (Subscript[r, 3]) && (Subscript[z, 1])==0

Analytic :

KerrGeoISCO[a,1]
KerrGeoISCO[a,-1]

Plot[{KerrGeoISCO[a,1],KerrGeoISCO[a,-1]},{a,-1,1},PlotTheme->"Detailed",ImageSize->700]

KerrGeoSeparatrix[.9,.5,.5]

Innermost stable spherical  orbit (ISSO)

(Subscript[r, 1])==(Subscript[r, 2])== (Subscript[r, 3])

KerrGeoISSO[.9,.5]

ContourPlot[
KerrGeoISSO[a,x],{x,-1,1},{a,(10^(-3)),1-(10^(-3))},PlotTheme->"Detailed",FrameLabel->{"x","a"}]

Innermost bound spherical orbit (IBSO)

e ==0 &&(Subscript[r, 1])==(Subscript[r, 2]) && ℰ==1

Analytic for polar orbits

KerrGeoIBSO[a,0]

Generic

KerrGeoIBSO[.9,.5]

Separatrix between stable and unstable orbits

(Subscript[r, 2])== (Subscript[r, 3])

With[{a=0.9,e=0.5,x=0.5},
KerrGeoSeparatrix[a,e,x]
]

Photon sphere

ℰ==∞
ℒ == ∞
e==0

Schwarzschild

KerrGeoPhotonSphereRadius[0,x]

Kerr (equatorial plane)

KerrGeoPhotonSphereRadius[a,1]
KerrGeoPhotonSphereRadius[a,-1]

Kerr (generic)

ContourPlot[
KerrGeoPhotonSphereRadius[a,x],{x,-1,1},{a,-1,1},PlotTheme->"Detailed",FrameLabel->{"x","a"}]

Parabolic or Marginally bound orbits (currently not implemented)

ℰ == 1
e == 1
(Subscript[r, 1])== ∞

Equatorial orbits

x==±1
(Subscript[z, 1])==0
𝒬 ==0

Polar orbits

x==0
(Subscript[z, 1])==1
ℒ==0

Spin-weighted spheroidal harmonics

The SpinWeightedSpheroidalHarmonics package provides a set of functions for computing spin-weighted spherical harmonics and their associated eigenvalues. The harmonics are solutions to the spin-weighed spheroidal equation

[(1/(sin θ))(d/(dθ))(sinθ (d/(dθ)))-(γ^2)(sin^2)θ-((((m+s cos θ))^2)/((sin^2)θ))-2γ s cos θ+s+2m γ+λ]S=0

Load package

<<SpinWeightedSpheroidalHarmonics`

Spin-weighted spherical harmonics

In the case γ=0, the spin-weighted spheroidal harmonics reduce to the spin-weighted spherical harmonics. We can evaluate these analytically:

SpinWeightedSphericalHarmonicY[-2,2,2,θ,ϕ]

or numerically:

SpinWeightedSphericalHarmonicY[-2,2,2,0.5,0.1]

We can also evaluate them to arbitrary precision:

SpinWeightedSphericalHarmonicY[-2,2,2,0.5`30,0.1`30]

For s=0 the spin-weighted spherical harmonics reduce to the scalar spherical harmonics

Table[SpinWeightedSphericalHarmonicY[0,l,m,θ,ϕ]-SphericalHarmonicY[l,m,θ,ϕ]//Simplify,{l,0,5},{m,-l,l}]//MatrixForm

Spin-weighted spheroidal eigenvalue

The eigenvalue is analytic in the γ=0 case

SpinWeightedSpheroidalEigenvalue[s,l,m,0]

For γ≠0, we can obtain approximate numerical values

SpinWeightedSpheroidalEigenvalue[2,3,2,0.1]

SpinWeightedSpheroidalEigenvalue can be quickly evaluated to very high precision.

Timing[SpinWeightedSpheroidalEigenvalue[2,3,2,0.1`200]]

For s=0 the eigenvalue is related to the SpheroidalEigenvalue[ ] function build into Mathematica

With[{l=2,m=2,γ=0.5`20},
SpinWeightedSpheroidalEigenvalue[0,l,m,γ]==
SpheroidalEigenvalue[l,m,I γ]-2 m γ
]

Spin-weighted spheroidal harmonic

The SpinWeightedSpheroidalHarmonicS function can be evaluated quickly to arbitrary precision

t1=Timing[SpinWeightedSpheroidalHarmonicS[-2,2,2,0.5`40,0.5`40,0.4`40]]

We can also compute derivatives of the harmonic w.r.t. θ and ϕ

D[SpinWeightedSpheroidalHarmonicS[2,2,2,0.1,θ,0.4],{θ,2}]/.θ->0.5

The above was fast but if you want to compute the harmonic for many different values of {θ,ϕ} then there is a quicker way. If you call the function without the θ,ϕ arguments you get an object that stores some precomputed steps and that can be evaluated quickly.

S=SpinWeightedSpheroidalHarmonicS[-2,2,2,0.5`40]

t2=Timing[S[0.5`40,0.4`40]]

t1/t2

This is ~70 times faster to run. This is fast enough to plot interactively using a manipulate function

Manipulate[Block[{S,s=-2,l=2,m=2},
S=SpinWeightedSpheroidalHarmonicS[s,l,m,γ];
Show[
ListPlot[Table[{θ,S[θ,0]},{θ,0,π,0.01}],Joined->True,PlotTheme->"Detailed",PlotStyle->Red,ImageSize->600],
Plot[SpinWeightedSphericalHarmonicY[s,l,m,θ,0],{θ,0,π}],
PlotRange->All
]
],{{γ,1.},0,3}]

Let’s do the same but for all the l=2 harmonics:

plotOptions:=Sequence[PlotRange->{0,1},Frame->True,FrameLabel->{"θ","Spin-weighted spheroidal harmonic (γ=2.0)"},PlotLegends->Placed[LineLegend[Table["m = "<>ToString[m],{m,-l,l}],LegendFunction->(Framed[#,RoundingRadius->5,Background->Directive[White,Opacity[0.8]]]&)],{0.85,0.71}],ImageSize->600];

Manipulate[Block[{s=-2,l=2},Module[{sSlm},
sSlm=Table[SpinWeightedSpheroidalHarmonicS[s,l,m,γ][θ,0],{m,-l,l}];Plot[sSlm,{θ,0,π},Evaluate[plotOptions]]]],{γ,0,2}]

We can also make a nice 3D plot show all the l=2 harmonics:

plot3DOptions:=Sequence[Axes->False,Boxed->False,PlotRange->{{-1,1},{-1,1},{-1,1}},BoxRatios->{1,1,1},PlotStyle->Directive[ColorData[97,m+l],Opacity[0.8]],ViewPoint->{Front,Right},SphericalRegion->True,PlotPoints->50];

Block[{s=-2,l=2,γ=2.0},Show[
Table[
sSlm=SpinWeightedSpheroidalHarmonicS[s,l,m,γ];
SphericalPlot3D[Abs[Re[sSlm[θ,ϕ]]],{θ,0,π},{ϕ,-π,π},Evaluate[plot3DOptions]],{m,-l,l}],ImageSize->800]]

Series expansions

Eigenvalue

We can use Mathematica’s Series command to expand about γ=0 for arbitrary {s,l,m}

Series[SpinWeightedSpheroidalEigenvalue[s,l,m,γ],{γ,0,1}]

It is quicker to perform the expansion for specific values of {s,l,m}

smallγExpansion=Series[SpinWeightedSpheroidalEigenvalue[2,3,2,γ],{γ,0,5}]

We can also expand around γ =∞. This only works for specific values of {s,l,m}

largeγExpansion=Series[SpinWeightedSpheroidalEigenvalue[2,3,2,γ],{γ,∞,5}]

Let’s compare the two expansions above to the numerical results

eigenvalueTable=Table[{γ,SpinWeightedSpheroidalEigenvalue[2,3,2,γ]},{γ,0,8,0.1}];

Show[
ListPlot[eigenvalueTable,PlotTheme->"Detailed",BaseStyle->18,FrameStyle->Black,FrameLabel->{"γ","(Subscript[λ, 32])"},ImageSize->600,PlotStyle->Darker[Green]],
Plot[Evaluate[Normal[smallγExpansion]],{γ,0,10},PlotStyle->Red],
Plot[Evaluate[Normal[largeγExpansion]],{γ,1,10},PlotStyle->Blue]
]

Harmonic

We can also use Series[ ] to expand the SpinWeightedSpheroidalHarmonic about γ=0 as a series of SpinWeightedSphericalHarmonics

Series[SpinWeightedSpheroidalHarmonicS[s,l,m,γ,θ,ϕ],{γ,0,1}]

Series[SpinWeightedSpheroidalHarmonicS[2,l,2,γ,π/2,0],{γ,0,3}]/.l->3

Solutions of the Regge-Wheeler and Teukolsky equations

The ReggeWheeler package provides a set of functions for computing solutions to the Regge-Wheeler equation, given by

[(f^2)((d^2)/(dr^2))+f f' (d/dr)+((ω^2)-(Subscript[V, l]))]Ψ=(Subscript[𝒯, RW])

where f=1-2M/r.

The Teukolsky package provides a set of functions for computing solutions to the Teukolsky equation, given by

[(Δ^(-s))(d/dr)((Δ^(s+1))(d/dr))+(((K^2)-2i s(r-M)K)/Δ)+4i s ω r-λ]R=𝒯

where Δ=(r^2)-2M r+(a^2) and K=((r^2)+(a^2))ω-a m.

Load packages

Before using the functions, first load the packages

<<ReggeWheeler`
<<Teukolsky`

Homogeneous solutions

We start by considering homogeneous solutions to the Regge-Wheeler and Teukolsky equations. These are generated using the ReggeWheelerRadial and TeukolskyRadial functions, respectively. A single homogeneous solution satisfying particular boundary conditions is represented by a ReggeWheelerRadialFunction or a TeukolskyRadialFunction.

With[{s=2,l=2,ω=0.1`20},
ψ=ReggeWheelerRadial[s,l,ω]
]

With[{s=-2,l=2,m=2,a=0.3`20,ω=0.1`20},
R=TeukolskyRadial[s,l,m,a,ω]
]

As we see from the result above, the default is to produce a pair of functions representing  “in” (purely ingoing from past infinity) and “up” (purely upgoing from the past horizon) boundary conditions.

Our normalisation convention is that the homogeneous solutions have transmission coefficient equal to 1.

These functions can be treated as a normal Mathematica functions that can be evaluated at any value of the radius. Let’s try with r=10

R["In"][10]

ψ["In"][10]

We can also compute derivatives in the normal way for functions

R["In"]'[10]

Higher derivatives are also supported

R["In"]''[10]

R["In"]'''''[10]

We can even evaluate on a grid of points

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

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

Finally, we can query a TeukolskyRadialFunction about its properties. For example, we can obtain the values of the asymptotic amplitudes (as expected from our normalisation convention, the transmission coefficient is 1)

R["In"]["Amplitudes"]

X["In"]["Amplitudes"]

Methods

There are currently three supported methods for computing homogeneous solutions: MST, NumericalIntegration and HeunC. The MST method is the default and can produce highly accurate results, but often needs to work beyond machine precision. The NumericalIntegration method can produce fast results on a grid of points, but is not always as accurate as the MST method. The HeunC method uses the new HeunC function included in Mathematica 12.1, but only currently supports computing solutions satisfying “in” boundary conditions. It is often fast, but fails to produce accurate (or even correct) results in some cases.

In most cases, all three methods produce consistent numerical results.

MST method

The MST method represents solutions of the Teukolsky equation in terms of a convergent series of hypergeometric functions, (see M. Sasaki and H. Tagoshi, Living Rev. Relativity, 6, 6 (2003) for details). The MST method is the default, but can also be specified using the Method -> “MST” option. The method has a single sub-option “RenormalizedAngularMomentum”, which is passed as the Method option to RenormalizedAngularMomentum. The default options are equivalent to the following:

R==With[{s=-2,l=2,m=2,a=0.3`20,ω=0.1`20},TeukolskyRadial[s,l,m,a,ω,Method->{"MST","RenormalizedAngularMomentum"->"Monodromy"}]]

Renormalized angular momentum

The renormalized angular momentum, ν,  is central to the MST method.

With[{s=-2,l=2,m=2,a=0.3`20,ω=0.1`20},
RenormalizedAngularMomentum[s,l,m,a,ω]
]

The renormalized angular momentum is in general complex. If we plot it, its real and imaginary parts appear to have an interesting structure.

With[{s=-2,l=2,m=2,a=0.3`20},
Plot[Evaluate[ReIm[RenormalizedAngularMomentum[s,l,m,a,ω]]],{ω,0,6}]//Quiet
]

This is a consequence of the interesting property that cos(2πν) is real and sometimes goes outside the range [-1,1]

With[{s=-2,l=2,m=2,a=0.3`20},
Plot[Cos[2π RenormalizedAngularMomentum[s,l,m,a,ω]],{ω,0,0.65}]//Quiet
]

The RenormalizedAngularMomentum function has a single option, Method. The default is “Monodromy” and is generally the most reliable choice provided it is given enough digits of working precision.

NumericalIntegration method

The NumericalIntegration method uses numerical differential equation solving (NDSolve) to compute solutions of the radial Teukolsky equation. It uses hyperboloidal slicing to maintain high accuracy and avoid typical problems associated with the long-range potential in the Teukolsky equation. The NumericalIntegration method can be specified using an option of the form Method -> {“NumericalIntegration”, “Domain” -> ...}.  There should be one numerical domain specified for each solution (satisfying a given boundary condition). For example,

With[{s=-2,l=2,m=2,a=0.3`20,ω=0.1`20},
Rin=TeukolskyRadial[s,l,m,a,ω,"BoundaryConditions"->"In",Method->{"NumericalIntegration","Domain"->{3,20}}]
]

Once a TeukolskyRadialFunction is constructed using the NumericalIntegration method, it can be efficiently evaluated on a grid of points.

Table[Rin[r],{r,3,20}]

HeunC method

The HeunC method uses Mathematica’s built-in HeunC method (available since version 12.1) to compute solutions of the radial Teukolsky equation. It is reasonably fast and accurate in many instances and is as reliable as the HeunC function itself (which works in many, but not all cases typically of interest). Currently, only “In” boundary conditions are supported as no version of HeunC directly maps onto a solution satisfying “Up” boundary conditions.

With[{s=-2,l=2,m=2,a=0.3`20,ω=0.1`20},
RHeun=TeukolskyRadial[s,l,m,a,ω,Method->"HeunC"]
]

All three methods produce the same result

RHeun["In"][10]

Rin[10]

R["In"][10]

NumericalIntegration:   faster (particularly with machine precision)
MST: high-accuracy
HeunC: quite fast, quite accurate but only for “in” solutions

Inhomogeneous solutions

Inhomogeneous solutions can be obtained by variation of parameters, integrating the homogeneous solutions against a source. In practice, this is rarely done by direct integration. Instead, efficient methods are developed on a case-by-case basis.

Inhomogeneous solutions are represented by a TeukolskyMode (or a ReggeWheelerMode), which stores information about the radial functions, angular functions and normalisation coefficients. So far, the cases of a point mass (s=-2), electric (s=-1) or scalar (s=0) charge on a circular, equatorial orbit is supported directly, although the framework exists for supporting any inhomogeneous solution of interest.

s=-2, circular, equatorial orbit

First, we construct an orbit using the KerrGeodesics package. We consider a circular orbit of radius (Subscript[r, 0])=10M in Kerr spacetime with a=0.5 and compute the s=2, l=2, m=2 mode

With[{a=0.5`32,r0=10.`32,e=0,x=1, s=-2,l=2,m=2,n=0,k=0},
orbit=KerrGeoOrbit[a,r0,e,x];
mode=TeukolskyPointParticleMode[s,l,m,n,k,orbit]
]

This has computed the asymptotic amplitudes at the horizon and infinity. We can access these from the TeukolskyMode object

mode["Amplitudes"]

We can also evaluate the radial and angular functions associated with the mode

mode["RadialFunctions"]

mode["AngularFunction"]

Fluxes of energy and angular momentum

It is often of physical interest to know the fluxes of energy an angular momentum in a given mode. These can be accessed as properties of the mode

mode["Fluxes"]

mode["EnergyFlux"]

mode["AngularMomentumFlux"]

Post-Newtonian expansions

The PostNewtonialSelfForce package provides analytic approximations of Black Hole Perturbation Theory quantities as post-Newtonian expansions.

Load PostNewtonialSelfForce package

<<PostNewtonianSelfForce`

Conventions and common expansion variables

({{Variable, Description, Relation to other variables}, {p, Semi-latus rectum, p=((2 r_min r_max)/(M (r_min + r_max)))}, {v, Characteristic velocity, v=(Sqrt[1/p])}, {y, Inverse separation, y=(((M (Subscript[Ω, ϕ])))^(2/3))}, {e, Eccentricity (Darwin), e=(((Subscript[r, max]) - (Subscript[r, min]))/((Subscript[r, max]) + (Subscript[r, min])))}, {q, Spin of the large (Kerr) black hole, q=a/M}})

Accessing post-Newtonian expansions

The PostNewtonianExpansion function is the main way to access any series available in the package. It takes a path-like list describing the desired quantity. For example, we can load the energy flux through future null infinity for an eccentric orbit around a Schwarzschild black hole. This returns a PostNewtonianSeries object

ℱ=PostNewtonianExpansion[{"Schwarzschild","Eccentric","Flux","EInfinity"}]

Properties of PostNewtonianSeries

ℱ["Name"]

ℱ["Description"]

ℱ["Authors"]

ℱ["References"]

Working with series

Loading the full series

ℱ["Series"]

ℱ["Series"]+O[y]^(17/2)

Accessing subsets of the series

PostNewtonianCoefficient[energy, order, logpower] functionality to retrieve individual PN terms.

petersMathews=PostNewtonianCoefficient[ℱ,5,0]

Resummed series

Some series have been resummed to pull out singular factors such as (1/(1-(e^2)))

resummedSeries=PostNewtonianCoefficient[ℱ,23/2,1]

We can separately extract the singular factor and the regular series

resummedSeries["SingularFactor"]

resummedSeries["SeriesData"]

We can also re-expand the singular factor

resummedSeries+O[e]^11