Getting Started with the Black Hole Perturbation Toolkit

In this tutorial, we will give an overview of  some 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["KerrGeodesics"];
PacletInstall["SpinWeightedSpheroidalHarmonics"];
PacletInstall["ReggeWheeler"];
PacletInstall["Teukolsky"];

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

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 for bound orbits: "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,10Λ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  hyperbolic bound orbits

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

KerrGeoOrbit[0.5,10,1.5,1]["ConstantsOfMotion"]

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"]

ψ["In"]["Amplitudes"]

Methods

There are multiple supported methods for computing homogeneous solutions including 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. 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"]]

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).

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

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 generic orbit is supported directly, and the framework exists for supporting other inhomogeneous solutions of interest.

s=-2, generic 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 at a specific radius:

mode[11]

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 solutions

We can obtain solutions as post-Newtonian expansions (currently only for circular orbits):

orbit=KerrGeoOrbit[a,p,0,1];

ψPN[2,2]=TeukolskyPointParticleModePN[0,2,2,orbit,{η,3}]

ψPN[2,2][r]

Self-force

p=10.`32;
lmax=10;

orbit=KerrGeoOrbit[0,p,0,1];

Monitor[Do[ψ[l,m]=TeukolskyPointParticleMode[0,l,m,orbit],{l,0,lmax},{m,-l,l}];,{l,m}]

Monitor[Do[Fr[l]=Re[Sum[ψ[l,m]["ExtendedHomogeneous" -> "ℐ"]'[p]SphericalHarmonicY[l,m,π/2,0],{m,-l,l,2}]];,{l,0,lmax}],{l,m}]

With[{E0=orbit["Energy"],r0=orbit["p"],L=orbit["AngularMomentum"],M=1},
Flr[-1]=-((E0 r0 (1+2 l) )/(2 (r0-2M) ((L^2)+(r0^2))));
Flr[0]=((2 (E0^2) (r0^3)+(2 M-r0) ((L^2)+(r0^2))) EllipticE[((L^2)/((L^2)+(r0^2)))]+(-(E0^2) (r0^3)+(2 M-r0) ((L^2)+(r0^2))) EllipticK[((L^2)/((L^2)+(r0^2)))])/(π r0 (-2 M+r0) ((((L^2)+(r0^2)))^(3/2)));
Flr[2]=((-8 (E0^4) (L-r0) (r0^10) (L+r0)+4 (E0^2) (r0^3) ((L^2)+(r0^2)) (9 (L^6) M+26 (L^4) M (r0^2)+23 (L^2) M (r0^4)+(L^2) (r0^5)+14 M (r0^6)-3 (r0^7))+(2 M-r0) ((((L^2)+(r0^2)))^2) (28 (L^6) M+82 (L^4) M (r0^2)+82 (L^2) M (r0^4)-(L^2) (r0^5)+32 M (r0^6)-3 (r0^7))) EllipticE[((L^2)/((L^2)+(r0^2)))]+((E0^4) (r0^10) (3 (L^2)-5 (r0^2))-(E0^2) (r0^5) ((L^2)+(r0^2)) (18 (L^4) M+34 (L^2) M (r0^2)+(L^2) (r0^3)+32 M (r0^4)-7 (r0^5))-(2 M-r0) (r0^2) ((((L^2)+(r0^2)))^2) (14 (L^4) M+28 (L^2) M (r0^2)+16 M (r0^4)-(r0^5))) EllipticK[((L^2)/((L^2)+(r0^2)))])/(2 (-1+2 l) (3+2 l) π (r0^6) (-2 M+r0) ((((L^2)+(r0^2)))^(7/2)));
]

Show[ListLogLogPlot[{Abs@Table[{l,Fr[l]},{l,0,lmax}],Abs@Table[{l,Fr[l]-Flr[-1]},{l,0,lmax}],Abs@Table[{l,Fr[l]-Flr[-1]-Flr[0]},{l,0,lmax}],Abs@Table[{l,Fr[l]-Flr[-1]-Flr[0]-Flr[2]},{l,0,lmax}]},Joined->True,PlotMarkers->Automatic]]

Show[ListPlot[Accumulate[Table[Fr[l]-A-B-Flr[2],{l,0,lmax}]],DataRange->{0,lmax}],Plot[1.378448 (10^(-5)),{l,0,lmax}]]