Hunting circular orbits of the Kerr metric

I'm going to try and re-derive 50 year old results in a day, programmatically, with modern computers, and see how close I can get.

We'll look at results for orbits in the Kerr metric, using a mathematical approach that can be scripted for other static, axis-symmetric spacetimes (though generic analytic results exist, to be explored in a future blog post), which is the underlying motivation; can we get a good enough of an approximation of tailored, analytic solutions to merit writing generic code for alternative metrics?

Of particular importance when studying metrics are the permitted circular orbits, especially for creating observational tests, as the orbits are used to model processes in e.g. accreting matter. Such trajectories also allow for the calculation of redshift values for discs when ray tracing, which in turn allows us to visualize and simulate other processes, such as flux. The little exploration here will hopefully lead to automated stable orbit discovery with our integrator, and (hopefully) allow us to eventually simulate arbitrary disc profiles for different metrics.


The main reference for this analysis is Bardeen, Press and Teukolsky (1972)[1], which explores orbits in rotating black holes, and various emergent processes. The paper covers an awful lot, but we are only interested (for now) in Section II.

The standard, static, axis-symmetric, and asymptotically flat spacetime line element is

ds2=e2νdt2+e2ψ(dϕωdt)2+e2μ1dr2+e2μ2dθ2, \text{d} s^2 = - \text{e}^{2\nu} \text{d} t^2 + \text{e}^{2\psi} \left( \text{d}\phi - \omega \text{d} t \right)^2 + \text{e}^{2\mu_1} \text{d} r^2 + \text{e}^{2\mu_2} \text{d} \theta^2,

where ω=gϕtgϕϕ\omega = -\frac{ g^{ }_{\phantom{ } \phi t}}{ g^{ }_{\phantom{ } \phi \phi}}.

We differentiate stationary and static static, where a stationary means the components of the metric do not change with time (tgμν=0\frac{\partial}{\partial t} g_{\mu\nu} = 0), and static adds symmetry under time parity, or simply that ttt \rightarrow -t does not change the metric.

Then axis-symmetric implies invariance under translations in ϕ\phi (things look the same when you rotate), and asymptotically flat means the metric reduces to the Minkowski metric, η\eta, for rr \rightarrow \infty.

These metrics govern the structure of spacetime in the vicinity of some central singularity, and are all we need to study a black hole and its orbits.

This metric becomes Kerr when

e2ν=ΣΔA,e2ψ=AΣsin2θ,e2μ1=ΣΔ,e2μ2=Σ \text{e}^{2\nu} = \frac{\Sigma \Delta}{A}, \quad \text{e}^{2\psi} = \frac{A}{\Sigma} \sin^2 \theta, \quad \text{e}^{2\mu_1} = \frac{\Sigma}{\Delta}, \quad \text{e}^{2\mu_2} = \Sigma

and ω=2MarA\omega = \frac{2Mar}{A}. The symbols have their usual meaning, e.g. from Cunningham (1975)[2], namely

Σ=r2+a2cos2θ,Δ=r22Mr+a2,A=(r2+a2)2a2Δsin2θ. \Sigma = r^2 + a^2 \cos^2 \theta, \quad \Delta = r^2 - 2 M r + a^2, \quad A = (r^2+a^2)^2 - a^2 \Delta \sin^2 \theta.

There's a future blog post in my take on the physical significance of all of these symbols, but that's not for today.

Circular orbits in the Kerr metric

It is a known result that circular orbits exist in the equatorial plane of static axis-symmetric spacetimes, i.e. where θ=π2\theta = \frac{\pi}{2} and θ˙=0\dot{\theta} = 0, as shown by Carter (1968)[3]. We can further constrain our mathematical picture by noting that circular orbits have no motion in the radial components,

u˙r=0,andru˙r=0. \dot{u}^r = 0, \quad \text{and} \quad \frac{\partial}{\partial r}\dot{u}^r = 0.

We can use the invariance of four velocities to express these constraints in terms of the metric,

gσνu˙σσu˙νν=μ2, g^{ }_{\phantom{ } \sigma \nu} \dot{u}^{ \sigma}_{\phantom{ \sigma} } \dot{u}^{ \nu}_{\phantom{ \nu} } = - \mu^2,

where, for time-like geodesics, the effective mass parameter μ=1\mu = 1. In order to rearrange for the radial component, this equation may be expanded,

grru˙rru˙rr+gttu˙ttu˙tt+2gtϕu˙ttu˙ϕϕ+gϕϕu˙ttu˙tt=μ2, g^{ }_{\phantom{ } rr} \dot{u}^{ r}_{\phantom{ r} } \dot{u}^{ r}_{\phantom{ r} } + g^{ }_{\phantom{ } tt} \dot{u}^{ t}_{\phantom{ t} } \dot{u}^{ t}_{\phantom{ t} } + 2 g^{ }_{\phantom{ } t\phi} \dot{u}^{ t}_{\phantom{ t} } \dot{u}^{ \phi}_{\phantom{ \phi} } + g^{ }_{\phantom{ } \phi\phi} \dot{u}^{ t}_{\phantom{ t} } \dot{u}^{ t}_{\phantom{ t} } = - \mu^2,

and, using the freedom to raise and lower indices under contraction, rewritten

(u˙rr)2=1grr(gtttt(u˙t)2+2gtϕtϕu˙tu˙ϕ+gϕϕϕϕ(u˙ϕ)2+μ2). \left( \dot{u}^{ r}_{\phantom{ r} } \right)^2 = \frac{-1}{ g^{ }_{\phantom{ } rr} } \left( g^{ tt}_{\phantom{ tt} } \left( \dot{u}^{ }_{\phantom{ } t} \right)^2 + 2 g^{ t\phi}_{\phantom{ t\phi} } \dot{u}^{ }_{\phantom{ } t} \dot{u}^{ }_{\phantom{ } \phi} + g^{ \phi\phi}_{\phantom{ \phi\phi} } \left( \dot{u}^{ }_{\phantom{ } \phi} \right)^2 + \mu^2 \right).

In this form, the usual identification is permitted for energy u˙t=E \dot{u}^{ }_{\phantom{ } t} = - E and angular momentum u˙ϕ=Lz \dot{u}^{ }_{\phantom{ } \phi} = L_z from the covariant momenta,

(u˙rr)2=1grr(gttttE22gtϕtϕELz+gϕϕϕϕLz2+μ2). \therefore \quad \left( \dot{u}^{ r}_{\phantom{ r} } \right)^2 = \frac{-1}{ g^{ }_{\phantom{ } rr} } \left( g^{ tt}_{\phantom{ tt} } E^2 - 2 g^{ t\phi}_{\phantom{ t\phi} } E L_z + g^{ \phi\phi}_{\phantom{ \phi\phi} } L_z^2 + \mu^2 \right).

The first constraint from eq. (4) may now be calculated from this expression, namely

u˙rr=±1eμ1e2ν(E+ωLz)2e2ΨLz2μ2=0. \dot{u}^{ r}_{\phantom{ r} } = \frac{\pm 1}{ \text{e}^{\mu_1} } \sqrt{ \text{e}^{-2\nu} \left( E + \omega L_z \right)^2 - \text{e}^{-2\Psi} L_z^2 - \mu^2 } = 0.

If we rephrase the above constraint in terms of effective potentials, we find the expression cleans up to

u˙rr=±R e2μ1=0. \dot{u}^{ r}_{\phantom{ r} } = \pm \sqrt{ \mathcal{R} \, \text{e}^{-2\mu_1} } = 0.

Our constraints may be modified now to R=0\mathcal{R} = 0 and rR=0\frac{\partial}{ \partial r}\mathcal{R} = 0, the latter of which is directly

(E+ωLz)2re2ν+2(E+ωLz)e2νLzrωLz2re2Ψ=0, \left( E + \omega L_z \right)^2 \frac{\partial}{\partial r} \text{e}^{-2\nu} + 2 \left( E + \omega L_z \right) \text{e}^{-2\nu} L_z \frac{\partial}{\partial r} \omega - L_z^2 \frac{\partial}{\partial r} \text{e}^{-2 \Psi} = 0,

calculated from eq. (9).

The term Rre2μ1\mathcal{R} \frac{\partial}{\partial r} \text{e}^{-2\mu_1} is omitted, as it is manifestly zero, since R=0\mathcal{R} = 0. Implicitly it is assumed that EE and LzL_z are constant with respect to rr over an orbit. Indeed, this assumption is valid as EE and LzL_z are constants of motion for any given geodesic, whereas the E(r)E(r) and Lz(r)L_z(r) which we later solve for are functions in the phase space. This distinction is subtle but important to note.

Here we will specialize to the Kerr metric, with subsequent manipulations performed using SageMath, and the occasional small bit of hand manipulation. Substituting in the metric components from eq. (2), we obtain

R=2Mr(μ2r2+(EaLz)2)+((E2μ2)(a2+r2)Lz2)Σ=0, \mathcal{R} = 2 M r {\left(\mu^{2} r^{2} + {\left(E a - L_z\right)}^{2}\right)} + {\left({\left(E^{2} - \mu^{2}\right)} {\left(a^{2} + r^{2}\right)} - L_z^{2}\right)} \Sigma = 0,


rR=6 Mμ2r2+2 (EaLz)2M+2 (E2μ2)rΣ+((E2μ2)(a2+r2)Lz2)rΣ=0. \frac{\partial}{\partial r} \mathcal{R} = 6 \, M \mu^{2} r^{2} + 2 \, {\left(E a - L_z\right)}^{2} M + 2 \, {\left(E^{2} - \mu^{2}\right)} r \Sigma + {\left({\left(E^{2} - \mu^{2}\right)} {\left(a^{2} + r^{2}\right)} - L_z^{2}\right)} \frac{\partial}{\partial r}\Sigma = 0.

It is worth noting that μ\mu has reappeared in the derivative term. This is due to the complexity of eq. (11) when taking the derivative in SageMath – instead, SageMath prefers to work with eq. (12) to produce more visually friendly results. Both will eventually solve the system identically, so the choice is largely stylistic.

In the equatorial plane, Σ=r2\Sigma = r^2; therefore

R=2Mr(μ2r2+(EaLz)2)+r2((E2μ2)(a2+r2)Lz2)=0 \mathcal{R} = 2 M r {\left(\mu^{2} r^{2} + {\left(E a - L_z\right)}^{2}\right)} + r^{2} {\left({\left(E^{2} - \mu^{2}\right)} {\left(a^{2} + r^{2}\right)} - L_z^{2}\right)} = 0


rR=6Mμ2r2+2M(EaLz)2+2r3(E2μ2)+2r((E2μ2)(a2+r2)Lz2)=0. \frac{\partial}{\partial r} \mathcal{R} = 6 M \mu^{2} r^{2} + 2 M {\left(E a - L_z\right)}^{2} + 2 r^3 \left( E^2 - \mu^2 \right) + 2r \left( {\left(E^{2} - \mu^{2}\right)} {\left(a^{2} + r^{2}\right)} - L_z^{2} \right) = 0.

Noticing a similar term in the two equations, we can rearrange eq. (14) for

2r((E2μ2)(a2+r2)Lz2)=4M(μ2r2+(EaLz)2), 2 r {\left({\left(E^{2} - \mu^{2}\right)} {\left(a^{2} + r^{2}\right)} - L_z^{2}\right)} = -4 M {\left(\mu^{2} r^{2} + {\left(E a - L_z\right)}^{2}\right)},

and substitute into eq. (15) to obtain

rR=6Mμ2r2+2M(EaLz)2+2r3(E2μ2)4M(μ2r2+(EaLz)2)=0, \frac{\partial}{\partial r} \mathcal{R} = 6 M \mu^{2} r^{2} + 2 M {\left(E a - L_z\right)}^{2} + 2 r^3 \left( E^2 - \mu^2 \right) - 4 M {\left(\mu^{2} r^{2} + {\left(E a - L_z\right)}^{2}\right)} = 0,

which may be simplified further:

rR=2Mμ2r2+2r3(E2μ2)2M(EaLz)2=0. \therefore \quad \frac{\partial}{\partial r} \mathcal{R} = 2 M \mu^{2} r^{2} + 2 r^{3} {\left(E^{2} - \mu^{2}\right)} - 2 M {\left(E a - L_z\right)}^{2} = 0.

Our constrains in eq. (4) are now eq. (14) and (18) respectively, which may be solved simultaneously to find EE and LzL_z satisfying circular orbits. In SageMath, this is fiddly, and sometimes SageMath will fail after a minute, or run for seemingly hours without terminating. But one soon develops a sense for what SageMath can and cannot do, and after this, it is generally quite straight forward to solve, even for other metrics. We find there are four solutions in total, representing the intersection of the positive and negative energy regimes with prograde and retrograde orbits,

Eμ=±r57Mr4+16M2r33(4M3+Ma2)r2+5M2a2r±2MaΔMrr5Mr2(2r2Mr+4Δ), \frac{E}{\mu} = \pm \sqrt{ \frac{ r^5 - 7Mr^4 + 16M^2r^3- 3(4M^3 + Ma^2)r^2 + 5M^2a^2r \pm 2 M a \Delta \sqrt{M r} } { r^5 - M r^2 \left(2 r^2 - Mr + 4 \Delta \right) } },

where the first ±\pm denotes the positive and negative energies, and the second ±\pm the prograde and retrograde respectively. The angular momentum can be rearranged as

Lzμ=Eμa±rM1r+r(Eμ)2. \frac{Lz}{\mu} = \frac{E}{\mu} a \pm \frac{r}{\sqrt{M}} \sqrt{ 1 - r + r \left(\frac{E}{\mu} \right)^2 }.

At the time of writing, this is the most simplified form I have been able to manipulate these equations into, which seems pretty okay since we haven't had to even touch any advanced techniques.

As an additional note, with SageMath, many of the preceding steps in preparing the constraints are not actually necessary, but would provide results that are unwieldy for such a blog post. I was also partially curious to see if simplifying the constraint equations dramatically would make any difference to the solutions (it didn't really).

A comparison to Bardeen et al. (1972)

The Bardeen et al. paper quotes a more beautiful analytic solution

Eμ=r322Mr12±aM12r3/4(r323Mr12±2aM12)12, \frac{E}{\mu} = \frac{ r^\frac{3}{2} - 2M r^\frac{1}{2} \pm a M^\frac{1}{2} }{ r^{3/4}\left( r^\frac{3}{2} - 3Mr^\frac{1}{2} \pm 2 a M^\frac{1}{2} \right)^\frac{1}{2} },


Lzμ=±M12(r22aM12r12+a2)r3/4(r323Mr12±2aM12)12, \frac{L_z}{\mu} = \frac{ \pm M^\frac{1}{2} \left( r^2 \mp 2 a M^\frac{1}{2} r^\frac{1}{2} + a^2 \right) }{ r^{3/4}\left( r^\frac{3}{2} - 3Mr^\frac{1}{2} \pm 2 a M^\frac{1}{2} \right)^\frac{1}{2} },

both of which are now ubiquitous in the literature.

They are also substantially different from eq. (19) and eq. (20). We can try to find the degree of similarity by plotting the equations with r>0.1r > 0.1 to avoid the central singularity:


The plots for energy (left panel) are in good agreement up until the energy becomes negative in the Bardeen case (though such energies are allowed in e.g. the Penrose process). The SageMath derivation fails here, presumably due to an ambiguity when squaring roots, which loses sign information. Indeed, we do have a negative energy solution from SageMath, but we would require a priori knowledge of which solution to apply in which domain. This wouldn't be ideal, and without a more refined solution, the SageMath approach would not even suggest negative energies are possible.

The consequence of the energy sign is also pronounced in the angular momentum (right panel). For the a=0.998a=0.998 case, the plots are still in good agreement, but the divergence for a=0.0a=0.0 is quite apparent, and worse still for a=0.988a=-0.988 (not depicted, as the diagram became difficult to annotate). Even when using the Bardeen et al. energy solution in the SageMath expression for LzL_z, some of the features of the curves do not disappear, such as the sharp cusp at r2r\sim2, which, again, is problematic.

Fortunately, however, the differences in the two solutions is for our purposes pedagogical, since the divergence lies within the photon orbit, at which radii stable circular massive orbits are long since prohibited.

Special orbit boundaries

There are many special boundaries in the Kerr metric. For purely circular orbits, we consider three of these boundaries and the respective possible orbits within them, as they are arguably most relevant to our study – these are the innermost stable circular orbits (ISCO), the marginally bound orbits, and the photon orbit.

There is of course also the event horizon, but this is not an orbital boundary. We will draw reference to this boundary, but only with the need of understanding that we cease to consider geodesics after they would pass this horizon.


The innermost stable circular orbit, as the name suggests, are the set of orbits which are both circular and energetically parabolic (stable). Stability requires the energy be a minimum; this may be rephrased as

2r2R0, \frac{\partial^2}{\partial r^2} \mathcal{R} \geq 0,

which solves for an inequality rrISCOr \geq r_\text{ISCO}, which is in the Bardeen paper. Physically, these are the orbits where perturbations to the potential will not cause the orbit to diverge; it is the energetically preferred state, if you will.

Marginally bound orbits

In direct contrast to the ISCO, we consider unstable circular orbits at r<rISCOr < r_\text{ISCO}. All such unbound circular orbits have hyperbolic energies (will diverge from circular orbits when perturbed), with a specific radius hosting the marginally bound orbits, defined as rr where Eμ=1\frac{E}{\mu} = 1. Any small perturbation of such a marginally bound orbit outwards will escape the black hole, but with asymptotically zero velocity.

The marginally bound radius is may be solved for the Kerr metric as

rmb=2Ma+2M2aM, r_{\text{mb}} = 2M \mp a + 2 \sqrt{M^2 \mp aM },

corresponding to first (infalling) radius at which the condition Eμ=1\frac{E}{\mu} = 1 is satisfied.

Photon orbit

Photons are null geodesics with μ=0\mu=0, and, due to this, cannot form circular orbits anywhere except at a one particular radius, known as the photon radius, denoted rphr_\text{ph}. This radius is where Eμ\frac{E}{\mu} \rightarrow \infty, or, equivalently, where

r23Mr±2aMr=0, r^2 - 3Mr \pm 2 a \sqrt{M r} = 0,

i.e. the denominator from eq. (21) vanishes. This equation is solved by Bardeen, with the solution

rph=2M[1+cos(23arccos(aM))]. r_\text{ph} = 2M \left[ 1 + \cos\left( \frac{2}{3} \arccos \left( \frac{ \mp a}{M} \right) \right) \right].

The photon orbit represents, if you will, the innermost possible (unstable) circular orbit, which requires velocity equal to the speed of light. Orbits between with rphr<rISCOr_\text{ph} \leq r < r_\text{ISCO} can only form unstable circular orbits, which may quickly and catastrophically collapse.

Analogous to the marginally bound orbits, the unstable orbits at the photon orbit, if perturbed outwards, will arrive at infinity still traveling at the speed of light.

Tracing circular orbits

With solutions for the energy and angular momentum, and a brief overview of relevant orbital boundaries, we can begin to simulate some of these orbits.

Tracing stable circular orbits when EE and LzL_z is known is relatively simple. We only need to calculate u˙ϕ\dot{u}^\phi, since our integrator automatically constrains u˙t\dot{u}^t via eq. (5). We find

u˙ϕ=gϕϕLzgttE, \dot{u}^\phi = g^{\phi\phi} L_z - g^{tt}E,

using our eq. (19) and eq. (20) with μ=1.0\mu=1.0 for time-like geodesics. Setting up the integrator is then quite straight forward:

using GeodesicTracer
using ComputedGeodesicEquations
using GeodesicBase

# define `energy` and `angmom` functions
Δ(M, a, r) = r^2 + a^2 - 2 * M * r
function energy(M, a, r)
    var = 2 * M * a * Δ(M, a, r) * sqrt(M * r)
    denom = r^5 - M * r^2 * (4 * Δ(M, a, r) + 2r^2 - M * r)
    en = (
        5 * M^2 * a^2 * r + 16 * M^2 * r^3 - 7 * M * r^4 + r^5 -
        3 * (4 * M^3 + M * a^2) * r^2 - var
    sqrt(abs(en / denom))

angmom(E, M, a, r) = E * a + (r / √M) * sqrt(abs(1 - r + r * E^2))
angmom(M, a, r) = angmom(energy(M, a, r), M, a, r)

# function for calculating the initial radial velocity
function phi_vel(m, u)
    metric = GeodesicBase.metric(m, u)
    inv_metric = inv(metric)

    E = my_energy(m.M, m.a, u[2])
    Lz = my_angmom(E, m.M, m.a, u[2])

    inv_metric[4, 4] * Lz - inv_metric[4, 1] * E

This code snippet defines utility functions for calculating the initial radial velocity component for a given metric and position. The actual tracing code is a single function call:

using StaticArrays

m = BoyerLindquist(M=1.0, a=0.0)
u = @SVector [0.0, 6.0, π/2, 0.0]
v = @SVector [0.0, 0.0, 0.0, phi_vel(m, u)]

sol = tracegeodesics(
    (0.0, 1200.0),
    μ = 1.0,
    abstol = 1e-10,
    reltol = 1e-10,

Tracing a couple of configurations and adding a boundary for the event horizon:


Note: for performance reasons, the integrator terminates at 1.1 rs1.1 \, r_\text{s} to avoid wasting time on doomed orbits.

Above are two orbits corresponding to circular EE and LzL_z for M=1.0M=1.0 and a=0.0a=0.0. The dashed line has an orbit at r=5.6r = 5.6, whereas the solid line r=6.0r = 6.0, i.e. an orbit just within the ISCO, and one exactly on the ISCO.

The integration is run for a total affine time of 1200 s1200\,\text{s}, as up until around 1000 s1000\,\text{s} both orbits are circular. In theory, the inner orbit is still above rmbr_\text{mb} and should be able to remain a circular orbit, but due to the long integration times, we see the machine error contribute enough that the orbit becomes unstable, collapses and falls into the black hole. Other values of r<rISCOr < r_\text{ISCO} produce radically different unstable orbits, which all diverge in one way or another if the integrator is run for long enough.

Orbit tracing regimes

We can investigate the stability of the integrator by setting the integration time to 10,000 s10,000\,\text{s}, and tracing all possible circular orbits for a specific MM and aa. A measure of deviation we employ is designed to tell us if the orbit collapsed inwards or outwards, and is calculated

Qd=1N rinit2i=1ri2, \mathcal{Q}_\text{d} = \frac{1}{N\,r_\text{init}^2} \sum_{i=1} r_i^2,

such that we can distinguish three cases:

Note that the final two cases are not sufficient conditions for the orbit to fall into or completely escape the black hole respectively. This measure instead tells us the majority of the orbit was either within or beyond the initial radius, but some orbits may have crossed back and forth.


This figure is for M=1.0M=1.0 and a=0.6a=-0.6, since the effect of aa is only to translate the distribution left and right, and this configuration illustrates the regimes clearly. The radii considered are in steps of Δr=0.01\Delta r = 0.01, and the pattern of Qd\mathcal{Q}_\text{d} is highly dependent on this, indicating that whether the orbit collapses or expands is due to machine precision. A few interesting things to note; the integrator is able to hold stable orbits all the way until the ISCO for long integration times. Between the marginally bound radius and the ISCO, deviations from the orbit are relatively small, which would indicate that the orbit remains circular for most of the integration time, before diverging – which is what we would expect from hyperbolic energies when dealing with numeric simulations. Between the photon orbit and the marginally bound orbit is where we see orbits that either collapse or expand extremely far away.

Within the photon horizon, the integration terminates after one or two steps towards the event horizon, which is why we gradually see Qd\mathcal{Q}_\text{d} approach unity, as only a single point remains in the integration result.

Particularly pronounced are the different boundaries, causing clearly different general behavior for the integrator. To investigate this further, we'll vary aa and introduce a measure of stability, coined as simply the normalized squared sum of the residuals in rr, expressed

Qs=1Ni=1(ririnit1)2, \mathcal{Q}_\text{s} = \sqrt{ \frac{1}{N} \sum_{i=1} \left( \frac{r_i}{r_\text{init}} - 1 \right)^2},

chosen to not penalize greater radii orbits for small fluctuations – hence divide by initial radius rinitr_\text{init} – and to not penalize orbits with many steps by dividing by total integrated points NN. The motivation for using Qs\mathcal{Q}_\text{s} over Qd\mathcal{Q}_\text{d} is to be able to clearly differentiate regions where circular orbits are possible without penalizing eventual divergence in the same way Qd\mathcal{Q}_\text{d} does.

Applying this measure for a range of aa and rr:


Here, Qs<109\mathcal{Q}_\text{s} < 10^{-9} have been coloured transparent.

Overlayed are the different orbit boundaries. These boundaries align well with the different integration regimes, which is always nice to see. There are however a few anomalies, namely those orbits close to the event horizon for highly retrograde black holes, a0.998a \rightarrow -0.998, where we see very poor stability, as the geodesics are launched out of the system from small energetic perturbations.

The small region of white within the event horizon in the bottom left corner is due to the integration being able to take a single step before it realizes any violation and terminates; it is, of course, a pure pathology.

Up until this point we have been able to use SageMath to derive results for EE and LzL_z for circular orbits which hold all the way until rphr_\text{ph}. We have also confirmed that these values lead to circular orbits, and that the stability (and presence) of these orbits aligns with the different special boundaries. But what if we don't have good solutions for EE and LzL_z?

Automated circular orbit discovery

Generalizing the last section; can we reliably and quickly find circular orbits without any a priori knowledge of the velocity vectors, other than that we are searching in the equatorial plane?

The temptation here is to brute-force search over the u˙ϕ\dot{u}^\phi parameter space until some stability condition is met, but we can be more elegant with this by using Optim.jl, which optimizes (i.e. finds the minima of) univariate functions. Optim.jl searches over a bounded range of any single argument function using the golden-section search algorithm.

For our purposes, the implementation needs a function which describes the degree of circularity of our orbits, which is effectively just Qs\mathcal{Q}_\text{s} from eq. (29).

The boilerplate is a few utility functions:

function trace_single_geodesic(m, u, vϕ)
    v = @SVector [0.0, 0.0, 0.0, vϕ]
        m, u, v, 
        # use a fairly long, but not excessively long integration time
        (0.0, 1000.0), 

Qs(rs) = sqrt(sum((rs ./ rs[1] .- 1.0).^2) / length(rs))

# utility function to trace and evaluate some trial `vϕ`
function estimate_stability(m, u, vϕ)
    sol = trace_single_geodesic(m, u, vϕ)
    rs = selectdim(sol, 1, 6)

Now we come to the crux, which is wrapping estimate_stability into a univariate function for Optim.jl. This is achieved through

using Optim

function find_vϕ_for_orbit(m, r_init; lower_bound=0.0, upper_bound=0.1)
    u = @SVector [0.0, r_init, deg2rad(90.0), 0.0]
    res = optimize(
        vϕ -> estimate_stability(m, u, vϕ),

The limitation of the univariate optimizer is the bounding window of possible values; estimating the bounding box too low risks excluding the true optimum, and estimating the bounding box too large risks local minima and bad high iteration count, which in turns leads to poor convergence. To try and mitigate this, we'll use a rolling window that continuously rescales these boundaries after each discovered orbit.

We can bake an assumption into this rolling window, which is that stable orbits closer to the singularity are more likely to have a higher initial u˙ϕ\dot{u}^\phi than those further out, and solve the system backwards, starting with a conservative flat-space estimate for the velocity, and progressively work inwards towards the black hole. Half of this assumption is pretty robust, since the metrics we are dealing with are asymptotically flat, but it may not always be the case that orbits closer to the singularity rotate faster.

In code, we express this rolling window as

function find_vϕ_for_orbit_range(m, rs; lower=0.0, upper=0.1)    
    map(reverse(rs)) do r
        vϕ = find_vϕ_for_orbit(m, r; lower_bound=lower, upper_bound=upper)
        # continuously adjust based on heuristic
        lower = vϕ * 0.9
        upper = vϕ * 2.0# call reverse twice so the order is maintained 
    end |> reverse

We can then solve for circular orbits over a range of rr with

m = BoyerLindquist(M=1.0, a=-0.4)
vϕ = find_vϕ_for_orbit(m, 10.0)

r_range = 1.1:0.1:10.0
vϕs = @time find_vϕ_for_orbit_range(m, r_range)

The solver takes on average 35 iterations (36 function calls) to converge to a minima of Qs\mathcal{Q}_\text{s} per orbit. Shorter integration times will also reliably converge in approximately the same number of function calls, with the routine being able to estimate stability after as little as a quarter of a full orbit reliably, up until approximately the ISCO. Within the ISCO, higher integration times are needed to accurately find solutions which illustrate the different orbital boundaries.

Calculating the energy for each orbit and plotting these again eq. (19) shows very good agreement to within 10910^{-9} for numerical accuracy:


Of course, the region within the photon orbit is impossible to probe with this technique, since the integrator cannot handle this domain. The negative energy regimes are therefore also impossible to discover with this method. Nevertheless, it proves reliable for finding stable orbits, up to the initial domain of the bounding window.

I'd also note here that this prototype code is able to discover those 100 stable orbits in just under a second, with approximately O(nn) scaling. Some aspects of this method could be optimized, but that is beyond the intention of this blog post.


We can get good approximations of analytic solutions using both SageMath and integrator-driven automated circular orbit discovery. Both methods are faced with limitations, however, in how they (fail to) probe within the photon radius.

In the case of the SageMath, the negative energy description of the Kerr metric is entirely lost, and the angular momenta within the photon radius are inaccurate at best. When using the integrator, we cannot even probe within this radius reliably. It also has a practical limitation; SageMath is an unreliable solver at the best of times. Often the expressions it produces when invoking its solve function are a mess, if it is able to solve them at all. I do not know how well this work will apply to metrics more complex than Kerr (which is the intended purpose after all). Though fortunately it seems that some fairly generic operations can go some distance to remedying this.

SageMath merits further investigation, and I will make the relevant SageMath notebooks for recreating the analysis in this blog available at some point.

The integrator suffers from an initial guess problem for the rolling window bounds on the possible range of u˙ϕ\dot{u}^\phi, but other than that is able to work without issue directly from the input metric. It is also able to describe enough of the space of circular orbits that the different orbital boundaries are visible, and will no doubt prove to be useful when studying other metrics.

Overall, I am pleased with how well the integrator performs for such orbits, and am confident these methods should be applicable to other spacetimes too.

Future work

The next step is to apply this work to general static, axis-symmetric spacetimes and modified Kerr metrics, and to assemble a full pipeline for either deriving and wrapping symbolic expressions into Julia functions, or a pipeline for automatic orbit discovery in the integrator ecosystem. This will then be used to calculate the redshift values for geometrically thin discs in the equatorial plane, to move the integrator in the direction of observational results purely from metrics.

Another avenue to later explore are the class of spherical orbits, which are not constrained to the equatorial plane, but nevertheless are able to maintain a constant radius.

A future post may also investigate ways of creating or even machine learning angular velocity tables for circular orbits, and how these could be shared and use in integration models.


[1] J M Bardeen, W H Press, S A Teukolsky (1972). Rotating black holes: locally nonrotating frames, energy extraction, and scalar synchrotron radiation. APJ, 178:347-169.
[2] C T Cunningham (1975). The effects of redshifts and focusing on the spectrum of an accretion disk around a Kerr black hole. APJ, 202:788-802.
[3] B Carter (1968). Global structure of the Kerr family of gravitational fields. PHYS REV 174:5.