I am pleased to announce a major update to the KerrGeodesics module of the Black Hole Perturbation Toolkit. This module calculates the properties of bound, timelike geodesics in Kerr spacetime and this major update brings the ability to work with fully generic (non-equatorial) geodesics. Key to the release of this update was the donation of code by Maarten van de Meent (who implemented the analytic formula of Fujita and Hikida) and Charles Evans, Zac Nasipak and Thomas Osburn (who implemented their own algorithm). This mixing of the best codes out there is exactly how I hoped the Toolkit would grow so a big thank you for these contributions. In addition to a slew of bug fixes, the documentation and website have also been updated.

It is now extremely easy to work with generic geodesics in Kerr spacetime. For instance, the geodesic plotted at the top of this post is generated using:

orbit = KerrGeoOrbit[0.998, 3, 0.6, Cos[π/4]];
{t, r, θ, φ} = orbit["Trajectory"];

Once you have the trajectory it is easily plotted via, e.g.,

   {r[λ] Sin[θ[λ]] Cos[φ[λ]], r[λ] Sin[θ[λ]] Sin[φ[λ]], r[λ] Cos[θ[λ]]}, {λ, 0, 20}, ImageSize -> 700, Boxed -> False, Axes -> False, 
   PlotStyle -> Red, PlotRange -> All],
 Graphics3D[{Black, Sphere[{0, 0, 0}, 1 + Sqrt[1 - 0.998^2]]}ß

In addition to computing the orbital trajectory there are many quantities that can be calculated. A complete list of the currently available functions is

KerrGeoEnergy[a, p, e, x]
KerrGeoAngularMomentum[a, p, e, x]
KerrGeoCarterConstant[a, p, e ,x]
KerrGeoConstantsOfMotion[a, p, e, x]
KerrGeoFrequencies[a, p, e, x]
KerrGeoOrbit[a, p, e, x]
KerrGeoPhotonSphereRadius[a, x]
KerrGeoISCO[a, x]
KerrGeoIBSO[a, x]
KerrGeoISSO[a, x]
KerrGeoSeparatrix[a, e, x]
KerrGeoBoundOrbitQ[a, p, e, x]

Many of these have additional options that can be set. For instance, the orbital frequencies can be computed w.r.t. Boyer-Lindquist or Mino time (see the documentation for more detail). Each orbit is parametrized by the black hole spin ‘a’, the semi-latus rectum ‘p’, the orbital eccentricity ‘e’ and the inclination angle ‘x_inc’. This is defined via

\( x_{inc} = Cos(\theta_{inc})\)    where     \( \theta_{inc} = \pi/2 – Sign(L_z)\theta_{min}  \)

where \( \theta_{min} \) is the minimum theta angle obtained during the orbital motion. Prograde orbits occur for \(x_{inc} > 0\), retrograde ones for \( x_{inc} <0 \).

So please download and use the KerrGeodesics package and give us feedback. There is a tutorial built into the Mathematica documentation and you can also find example notes books in the Mathematica Toolkit Examples section of the Toolkit. If you find any bugs or there’s a feature missing you would like to see implemented, please report them to the issue tracker.