Plunging photons out of stable orbits

Following work from a previous blog post concerning circular orbits in the Kerr metric, this post will investigate the region within the innermost stable circular orbit (ISCO). I'd also like to re-derive results for general static, axis-symmetric metrics, such that the energy EE and angular momentum LzL_z from the previous blog post may be more compactly expressed. This will differ from the SageMath method by introducing the Keplerian frequency, and manipulating the equations by hand.

I aim to calculate velocity vectors of massive geodesics in the so-called plunging region, i.e. the region within the ISCO, and use this to calculate redshifts in our class of spacetimes. In the literature, these results rely on having explicit first-order equations of motion. Since LzL_z, Carter's constant QQ, and EE are known for orbits at the ISCO, the plunging region trajectories are calculated by substituting these into the first-order equations, thereby lifting the condition on the radial component of the velocity.

In practice, this is equivalent to pushing the photon an infinitesimal amount towards the black hole. Consequently, if these first-order equations are unknown, we can instead create numerical approximations for the analytic solutions by tracing this slightly offset trajectory with Gradus.jl (I finally settled on a name).

Stable orbits in static, axis-symmetric metrics

Following the previous blog post, and elaborating on the method of Johannsen (2013)[1], we start by restricting ourselves to θ=π2\theta = \frac{\pi}{2}, with θ˙=0\dot{\theta} = 0, and require, for stable orbits, that

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

Here, as in the previous post, the dot denotes differentiation with respect to the affine parameter λ\lambda.

We write the geodesic equation

ddλu˙σ+Γσνεσu˙νu˙ε=0, \frac{\text{d}}{\text{d} \lambda} \dot{u}^\sigma + \Gamma^{ \sigma}_{\phantom{ \sigma} \nu \varepsilon} \dot{u}^\nu \dot{u}^\varepsilon = 0,

and use the definition of the Christoffel symbols,

Γσνεσ=12gσα(uνgαε+uεgανuαgνε). \Gamma^{ \sigma}_{\phantom{ \sigma} \nu \varepsilon} = \frac{1}{2} g^{\sigma \alpha} \left( \frac{\partial}{\partial u^\nu} g_{\alpha \varepsilon} + \frac{\partial}{\partial u^\varepsilon} g_{\alpha \nu} - \frac{\partial}{\partial u^\alpha} g_{\nu \varepsilon} \right).

to expand the geodesic equation as

gσαddλu˙σ+12(uνgαε+uεgανuαgνε)u˙νu˙ε=0. g_{\sigma \alpha} \frac{\text{d}}{\text{d} \lambda} \dot{u}^\sigma + \frac{1}{2} \left( \frac{\partial}{\partial u^\nu} g_{\alpha \varepsilon} + \frac{\partial}{\partial u^\varepsilon} g_{\alpha \nu} - \frac{\partial}{\partial u^\alpha} g_{\nu \varepsilon} \right) \dot{u}^\nu \dot{u}^\varepsilon = 0.

Notice the indices of the terms in bracket may be interchanged in the sum over the velocity vectors given the symmetry of the tensors involved, permitting us to combine terms and rearrange

gσαddλu˙σ+u˙νu˙εuνgαε=12u˙μu˙βuαgμβ. g_{\sigma \alpha} \frac{\text{d}}{\text{d} \lambda} \dot{u}^\sigma + \dot{u}^\nu \dot{u}^\varepsilon \frac{\partial}{\partial u^\nu} g_{\alpha \varepsilon} = \frac{1}{2} \dot{u}^\mu \dot{u}^\beta \frac{\partial}{\partial u^\alpha} g_{\mu \beta}.

Next, we (ab)use the chain rule to express

u˙εddλgαε=u˙εddλuνuνgαε=u˙εu˙νuνgαε, \dot{u}^\varepsilon \frac{\text{d}}{\text{d} \lambda} g_{\alpha \varepsilon} = \dot{u}^\varepsilon \frac{\text{d}}{\text{d} \lambda} u^\nu \frac{\partial}{\partial u^\nu} g_{\alpha \varepsilon} = \dot{u}^\varepsilon \dot{u}^\nu \frac{\partial}{\partial u^\nu} g_{\alpha \varepsilon},

and therefore

gσαddλu˙σ+u˙εddλgαε=12u˙μu˙βuαgμβ, g_{\sigma \alpha} \frac{\text{d}}{\text{d} \lambda} \dot{u}^\sigma + \dot{u}^\varepsilon \frac{\text{d}}{\text{d} \lambda} g_{\alpha \varepsilon} = \frac{1}{2} \dot{u}^\mu \dot{u}^\beta \frac{\partial}{\partial u^\alpha} g_{\mu \beta},

or, using the product rule,

ddλ(gσαu˙σ)=12u˙μu˙βuαgμβ. \frac{\text{d}}{\text{d} \lambda} \left( g_{\sigma\alpha} \dot{u}^\sigma \right) = \frac{1}{2} \dot{u}^\mu \dot{u}^\beta \frac{\partial}{\partial u^\alpha} g_{\mu \beta}.

Examining just the radial derivative components (α=r\alpha = r), along with the condition that u˙r=0\dot{u}^r = 0, forces the LHS to disappear for our class of metrics, with the only non-zero terms of the RHS for stable orbits being

0=rgtt(u˙t)2+2rgtϕ u˙t u˙ϕ+rgϕϕ(u˙ϕ)2. 0 = \frac{\partial}{\partial r} g_{tt} (\dot{u}^t)^2 + 2 \frac{\partial}{\partial r} g_{t\phi} \, \dot{u}^t \dot{u}^\phi +\frac{\partial}{\partial r} g_{\phi \phi} (\dot{u}^\phi)^2.

This is just an alternative method to achieve a result from the previous blog post, which I wanted to derive motivated by the Johannsen paper.

Keplerian frequency

We introduce the Keplerian frequency, derived from Kepler's third law, which is the rate at which the angular coordinate changes with coordinate time. As an aside, the Keplerian frequency also allows us to calculate the rotation limits of various astrophysical bodies, as a star rotating faster than this frequency would promote particles into higher orbits, effectively ripping itself apart (something for another post).

The Keplerian frequency for our purposes may be defined

Ωϕ:=u˙ϕu˙t, \Omega_\phi := \frac{\dot{u}^\phi}{\dot{u}^t},

and therefore Eq. (9) may be expressed

Ωϕ2 rgϕϕ+2Ωϕ rgtϕ+rgtt=0, \Omega_\phi^2 \, \frac{\partial}{\partial r} g_{\phi \phi} + 2 \Omega_\phi \, \frac{\partial}{\partial r} g_{t \phi} + \frac{\partial}{\partial r} g_{t t} = 0,

inviting us to solve for the Keplerian frequency

Ωϕ=(rgtϕ±(rgtϕ)2rgttrgϕϕ )(rgϕϕ)1. \Omega_\phi = \left( - \frac{\partial}{\partial r} g_{t \phi} \pm \sqrt{ \left( \frac{\partial}{\partial r} g_{t \phi} \right)^2 - \frac{\partial}{\partial r} g_{t t} \frac{\partial}{\partial r} g_{\phi \phi} } \,\right) \left( \frac{\partial}{\partial r} g_{\phi \phi} \right)^{-1}.

As an exercise, substituting in the Kerr metric terms yields

Ωϕ=1Ma±r3/2. \Omega_\phi = \frac{1}{ M a \pm r^{3/2} }.

Energy and angular momentum

To solve for the energy and angular momentum of a given stable circular orbit, we use that

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

expanding components to find

(u˙rr)2=1grr(gttttE22gtϕtϕELz+gϕϕϕϕLz2+μ2), \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),

where we have identified covariant momenta u˙t=E\dot{u}_t = -E and u˙ϕ=Lz\dot{u}_\phi = L_z, and used to the freedom to raise and lower indices under contraction. We will write an effective potential (u˙r)2=Veff(r)(\dot{u}^r)^2 = V_\text{eff}(r), and solve

Veff(r)=0,andddrVeff(r)=0, V_\text{eff}(r) = 0, \quad \text{and} \quad \frac{\text{d}}{\text{d} r} V_\text{eff}(r) = 0,

simultaneously for EE and LzL_z. This approach is very similar to last time, except now we will use the Keplerian frequency in Eq. (10), with the observation that u˙ϕ=gϕϕLzgtϕE\dot{u}^\phi = g^{\phi \phi} L_z - g^{t\phi} E and u˙t=gtϕLzgttE\dot{u}^t = g^{t \phi} L_z - g^{t t} E, inferred from the operation of the metric. Consequently, the Keplerian frequency is

Ωϕ=gϕϕLzgtϕEgtϕLzgttE. \Omega_\phi = \frac{g^{\phi \phi} L_z - g^{t\phi} E}{g^{t \phi} L_z - g^{t t} E}.

Johannsen uses the definition of Ωϕ\Omega_\phi in Eq. (12) along with Eq. (16) and Eq. (17) to derive the energy and angular momentum, but I was unable to recreate this. Instead, I came up with method to derive these using only the latter two and the metric condition

gμνgμρ=δνρν, g^{\mu \nu} g_{\mu \rho} = \delta^{ \nu}_{\phantom{ \nu} \rho},

manifest for any metric tensor. In particular, for the class of static, axis-symmetric metrics, with ν=ρ\nu = \rho one obtains

gtϕ gtϕ+gϕϕ gϕϕ=gtϕ gtϕ+gtt gtt=1, g^{t \phi} \, g_{t \phi} + g^{\phi \phi} \, g_{\phi \phi} = g^{t \phi} \, g_{t \phi} + g^{t t} \, g_{t t} = 1,

and with ν = ρ\nu \ \cancel{=} \ \rho,

gtϕ gtt+gϕϕ gtϕ=gtϕ gϕϕ+gtt gtϕ=0. g^{t \phi}\, g_{tt} + g^{\phi\phi} \, g_{t \phi} = g^{t \phi}\, g_{\phi \phi} + g^{tt}\, g_{t \phi} = 0.

We begin by re-arranging Eq. (17) into

Lz=gtϕgttΩϕgϕϕgtϕΩϕE. L_z = \frac{g^{t \phi} - g^{t t} \Omega_\phi}{g^{\phi \phi} - g^{t\phi} \Omega_\phi} E.

Multiplying in the fraction by gtϕg_{t\phi}

Lz=gtϕ gtϕgtϕ gttΩϕgtϕ gϕϕgtϕ gtϕΩϕE, L_z = \frac{g_{t\phi} \, g^{t \phi} - g_{t\phi} \, g^{t t} \Omega_\phi}{g_{t\phi} \, g^{\phi \phi} - g_{t\phi} \, g^{t\phi} \Omega_\phi} E,

and using Eq. (20) to substitute relevant symbols, we find

Lz=(gtϕ+gϕϕΩϕgtt+gtϕΩϕ)E= AB E, L_z = - \left( \frac{g_{t \phi} + g_{\phi \phi} \Omega_\phi}{g_{tt} + g_{t \phi} \Omega_\phi} \right) E = - \, \frac{\mathcal{A}}{\mathcal{B}}\, E,

where we use A\mathcal{A} and B\mathcal{B} for convenience to denote the numerator and denominator respectively. Substituting this expression for LzL_z into Eq. (16) yields

E2[gtt+2gtϕ(AB)+gϕϕ(AB)2]+μ2=0, E^2 \left[ g^{tt} + 2 g^{t\phi} \left( \frac{\mathcal{A}}{\mathcal{B}} \right) + g^{\phi \phi}\left( \frac{\mathcal{A}}{\mathcal{B}} \right)^2 \right] + \mu^2 = 0,

or equivalently

E2μ2=B2gttB22gtϕABgϕϕA2, \frac{E^2}{\mu^2} = \frac{ \mathcal{B}^2 }{ - g^{tt} \mathcal{B}^2 - 2 g^{t\phi} \mathcal{A} \mathcal{B} - g^{\phi \phi} \mathcal{A}^2 },

and, by symmetry of Eq. (23), one can immediately write

Lz2μ2=A2(gttB2+2gtϕAB+gϕϕA2). \frac{L_z^2}{\mu^2} = \frac{ \mathcal{A}^2 }{ - \left( g^{tt} \mathcal{B}^2 + 2 g^{t\phi} \mathcal{A} \mathcal{B} + g^{\phi \phi} \mathcal{A}^2 \right) }.

Examining this common denominator by expanding all A\mathcal{A} and B\mathcal{B} and collecting powers of Ωϕ\Omega_\phi gives us

Ωϕ2(gtt gtϕ gtϕ+2gtϕ gtϕ gϕϕ+gϕϕ gϕϕ gϕϕ)+ 2 Ωϕ(gϕϕ gtϕ gtϕ+gtϕ gϕϕ gtt+gtt gtt gtϕ+gϕϕ gϕϕ gtϕ)+ 2 Ωϕ(gtt gtt gtt+2gtϕ gtϕ gtt+gϕϕ gϕϕ gϕϕ),\begin{aligned} \Omega_\phi^2 &\left( g^{tt} \, g_{t\phi} \, g_{t\phi} + 2 g^{t\phi} \, g_{t\phi} \, g_{\phi\phi} + g^{\phi\phi} \, g_{\phi\phi} \, g_{\phi\phi} \right)\\ +\ 2 \,\Omega_\phi &\left( g^{\phi\phi} \, g_{t\phi} \, g_{t\phi} + g^{t \phi} \, g_{\phi\phi} \, g_{tt} + g^{tt} \, g_{tt} \, g_{t\phi} + g^{\phi \phi} \, g_{\phi\phi} \, g_{t\phi} \right)\\ +\ \phantom{2 \,\Omega_\phi} &\left( g^{tt} \, g_{tt} \, g_{tt} + 2g^{t\phi} \, g_{t\phi} \, g_{tt} + g^{\phi\phi} \, g_{\phi\phi} \, g_{\phi\phi} \right), \end{aligned}

dropping the negative sign temporarily.

The next steps are pretty methodological, so I will demonstrate the mechanics for Ωϕ2\Omega_\phi^2 terms, and then state the results for the other powers in the interest of brevity. First we use Eq. (19) to simplify the first term

Ωϕ2(gtϕ gtϕ gϕϕ+2gtϕ gtϕ gϕϕ+gϕϕ gϕϕ gϕϕ), \Omega_\phi^2 \left( -g^{t\phi} \, g_{t\phi} \, g_{\phi\phi} + 2 g^{t\phi} \, g_{t\phi} \, g_{\phi\phi} + g^{\phi\phi} \, g_{\phi\phi} \, g_{\phi\phi} \right),

which may then be combined to write

Ωϕ2 gϕϕ(gtϕ gtϕ+gϕϕ gϕϕ), \Omega_\phi^2 \, g_{\phi\phi} \left( g^{t\phi} \, g_{t\phi} + g^{\phi\phi} \, g_{\phi\phi} \right),

and, via Eq. (20),

Ωϕ2 gϕϕ. \implies \Omega_\phi^2 \, g_{\phi\phi}.

Similarly for the other powers. After working our way through this piece of algebra, we obtain

Eμ=±gtt+gtϕ Ωϕgϕϕ Ωϕ22gtϕ Ωϕgtt, \frac{E}{\mu} = \pm \frac{g_{tt} + g_{t\phi} \, \Omega_\phi}{ \sqrt{ -g_{\phi \phi} \, \Omega_\phi^2 - 2 g_{t\phi} \, \Omega_\phi - g_{tt} } },


Lzμ=±gtϕ+gϕϕ Ωϕgϕϕ Ωϕ22gtϕ Ωϕgtt. \frac{L_z}{\mu} = \pm \frac{g_{t\phi} + g_{\phi\phi} \, \Omega_\phi}{ \sqrt{ -g_{\phi \phi} \, \Omega_\phi^2 - 2 g_{t\phi} \, \Omega_\phi - g_{tt} } }.

Tracing orbits with the analytic constants

To verify these results, we can attempt to trace orbits around a Kerr black hole. We first define the relevant formulae

using Gradus
using StaticArrays

# implementation for energy
# fix ourselves to the negative energy solution arbitrarily
function energy(g, Ωϕ)
    -(g[1] + g[5] * Ωϕ) / √(-g[1] - 2g[5]*Ωϕ - g[4]*Ωϕ^2)

# implementation for angular momentum
# include a flag for whether we want the prograde or retrograde orbit
function angmom(g, Ωϕ, prograde::Bool)
    res = (g[5] + g[4] * Ωϕ) / √(-g[1] - 2g[5]*Ωϕ - g[4]*Ωϕ^2)
    if prograde

# keplerian frequency
# requires choice of positive or negative domain
function Ω(m, rθ, positive::Bool)
    # require derivatives with respect to radial coordinate
    jacs = Gradus.GeodesicTracer.metric_jacobian(m, rθ)
    ∂rg = jacs[:, 1]

    Δ = √(∂rg[5]^2 - ∂rg[1] * ∂rg[4])
    if positive 
        -(∂rg[5] + Δ) / ∂rg[4]
        -(∂rg[5] - Δ) / ∂rg[4]

# calculate all relevant metric quantities for stable circular orbit system
function stable_system(m, r; contra_rotating=false, prograde=true)
    let rθ = @SVector([r, π/2])
        g = Gradus.GeodesicTracer.metric_components(m, rθ)
        ginv = Gradus.GeodesicTracer.inverse_metric_components(g)

        Ωϕ = Ω(m, rθ, contra_rotating)
        E = energy(g, Ωϕ)
        L = angmom(g, Ωϕ, prograde)


This code is all relatively straight forward, using Gradus.jl to calculate the metric components for us, and just using them where needed. We can then write a utility function for creating stable 4-velocities:

function fourvelocity(m, r; kwargs...)
    g, ginv, E, L = stable_system(m, r; kwargs...) 
    vt = -E * ginv[1] + L * ginv[5]
    vϕ = -E * ginv[5] + L * ginv[4]

    @SVector[vt, 0.0, 0.0, vϕ]

and trace as usual

# choice of metric
m = BoyerLindquistAD(M=1.0, a=0.0)

# range of radii
rs = range(Gradus.isco(m), 10.0, 7)

# generate position and init velocities
us = [@SVector([0.0, r, deg2rad(90), 0.0]) for r in rs]
vs = map(r -> fourvelocity(m, r), rs)

# trace the geodesics for a fairly long duration to ensure stability
simsols = @time tracegeodesics(
    us, vs,
    (0.0, 2000.0);
    abstol=1e-10, reltol=1e-10,


Excellent, the orbits are indeed perfectly stable!

Plunging orbits

The plunging region is the region of unstable circular orbits, where the energies are hyperbolic – any small deviation from the trajectory will cause the particle to spiral inwards towards the black hole. As stated in the introduction, these trajectories are traditionally calculated using analytic first-order equations of motion, by substituting EE and LzL_z at the ISCO, along with the condition of equatorial orbits Q=0Q=0, and calculating u˙ν\dot{u}^\nu.

Four-velocity calculations

Instead, we will attempt to calculate u˙ν\dot{u}^\nu from metric terms at the ISCO.

We need minimally u˙r\dot{u}^r and u˙ϕ\dot{u}^\phi, since u˙θ=0\dot{u}^\theta = 0 and u˙t\dot{u}^t is constrained by the integrator. We already have u˙ϕ\dot{u}^\phi from the previous section, and u˙r\dot{u}^r may be found from Eq. (15).

A possible implementation for our plunging four-velocity may then be

function plunging_fourvelocity(m, r; kwargs...) where {T}
    g, ginv, E, L = stable_system(m, r; kwargs...) 
    vt = -E * ginv[1] + L * ginv[5]
    vϕ = -E * ginv[5] + L * ginv[4]
    nom = ginv[1] * E^2 - 2ginv[5] * E * L + ginv[4] * L^2
    denom = -g[2]

    @SVector[vt, -sqrt(abs(nom/denom)), 0.0, vϕ]

Note that we have hardcoded the negative choice of radial component, since we are infalling. Tracing this with

u_plunging = @SVector[0.0, Gradus.isco(m), deg2rad(90), 0.0]
v_plunging = CircularOrbits.plunging_fourvelocity(m, Gradus.isco(m))

sol = @time tracegeodesics(
    u_plunging, v_plunging,
    # need a very large integration time, else we never leave the ISCO
    (0.0, 18_000.0);
    abstol=1e-11, reltol=1e-11,

and plotting:


This is good, but the long integration times are perhaps not ideal, since most of the time is just spent tracing the circular orbit.

Dropping photons from sub-ISCO radii

An alternative method to the above is to drop the photon into the black hole, by moving the orbiting radius a tiny fraction closer. This only requires the u˙ϕ\dot{u}^\phi component of the velocity at the ISCO, and introduces δr\delta r as a free parameter which we can tweak. In code

u_dropped = @SVector[0.0, Gradus.isco(m) - δr, deg2rad(90), 0.0]
v_dropped = fourvelocity(m, Gradus.isco(m))

sol = @time tracegeodesics(
    u_dropped, v_dropped,
    (0.0, 500.0);
    abstol=1e-11, reltol=1e-11,

Tracing for different δr\delta r:


Here, the solid line is δr=1011\delta r = 10^{-11}, whereas the dashed line has δr=0.001\delta r = 0.001. The latter has significantly lower integration time, as one might expect, but falls nearly identically, modulo some initial rotation uϕu^\phi. Plotting the components of velocity as a function of rr:


The above plot requires the keyword argument closet_approach=1.0001 to be passed to tracegeodesics, as default termination is 1% beyond the inner horizon radius. Also, I used a=0.2a=0.2 in the above plot to include non-zero gtϕg_{t\phi} terms.

The Cunningham values are from Cunningham (1975)[2], calculated analytically from first-order equations of motion, and evenly sampled for the sake of illustration.

Note that all three traced curves agree: the one directly from plunging_fourvelocity, and the two dropped from a slightly closer radius. Setting δr\delta r too high causes the paths to diverge from expectation.

Interpolating plunging velocities

We can attempt to check the deviation from Cunningham's calculations of the above methods – to do this, we need to interpolate the velocity component calculated by the tracer as a function of radius. This may be done with Interpolations.jl:

using Interpolations

# need to sort our values
# drop the first index as it is duplicated
I = sortperm(pl_sol[2,:])[2:end]
r_sol = sol[2,:][I]
ut_sol = sol[5,:][I]

interp_ut = LinearInterpolation(r_sol, ut_sol)

Comparing the values with

using Statistics

isco = Gradus.isco(m)
# range of points to compute error over
# need to shrink the domain a little to avoid 
# stepping outside of the interpolated domain
r_range = range(
    Gradus.inner_radius(m) + 0.1, 
    isco - 0.1, 

# Cunningham analytic calculation
ut_true = map(
    r -> Gradus.AccretionFormulae.RedshiftFunctions.uᵗ(m.M, isco, r, m.a), 

# find interpolated values
test_t = map(interp_ut, r_range)

# absolute error
err_t = abs.(test_t .- ut_true)

The same is done for uru^r and uϕu^\phi. These components are then used to create a measure of error distance between the true and interpolated velocities, using the l2l_2 norm,

err=(Δut)2+(Δur)2+(Δuϕ)2, \text{err} = \sqrt{(\Delta u^t)^2 + (\Delta u^r)^2 + (\Delta u^\phi)^2},

where Δuν\Delta u^\nu denotes the difference between the velocity component of the two methods.

As an example, for various offsets, this error by radius is:


The above figure shows the rolling mean, with window size 10310^3 over a total 10510^5 points. Interesting to note is that the error close to the event horizon is left relatively unchanged for different δr\delta r. Instead, perhaps unsurprisingly, the offset primarily impacts the initial error at r=riscor = r_\text{isco}. To me, this implies there may be some error in the integration over this region that I will need to check, but investigating this is beyond the scope of this post.

To examine how this error changes with offset and integrator tolerance, we will plot the mean and standard deviation over the lifetime of the plunging trajectory. For δr\delta r in some range between 101110^{-11} and 10110^{-1}, we find these scale:


Both measures scale relatively linearly up until some cutoff δr\delta r, at which point we find a minimum deviation in error, before a plateau. This plateau and the location of the cutoff scales with the reltol of the integrator, implying the mean error is indeed a numerical error, and not directly a methodological error.

Note that the overall best mean error is very approximately reltol\sqrt{\texttt{reltol}}, which is useful for anticipating expected errors, and that this best-mean plateau also occurs at approximately δr=reltol10\delta r = \frac{\sqrt{\texttt{reltol}}}{10}, giving us some way to tune the offset to the integrator tolerances.

We can, however, overall expect the relative error to be no worse than 1 part in 100 using this interpolation method, which should be sufficient for our purposes (though I will note this seems quite high to me). We can test this by recreating a known rendering result, namely the relative redshift of an equatorial accretion disc.

Plunging redshift

A standard result in relativity is that the energy of a particle with four-momentum pνp^\nu, measured by an observer moving with u˙μ\dot{u}^\mu is

E=gμνpμu˙ν>0. E = - g_{\mu \nu} p^\mu \dot{u}^\nu > 0.

If we consider the redshift of the accretion disc as viewed by a distant observer as the ratio of energies measured at the disc and at the observer,

z=EobsEdisc=gσρ pσ u˙obsρgμν pμ u˙discν, z = \frac{E_\text{obs}}{E_\text{disc}} = \frac{ g_{\sigma \rho} \, p^\sigma \, \dot{u}_\text{obs}^\rho }{ g_{\mu \nu} \, p^\mu \, \dot{u}_\text{disc}^\nu },

then we can colour the redshift for each pixel in our traced image. Using Tullio.jl for the Einsum notation, we can easily express a PointFunction which will do this for us:

using Tullio

# need the ISCO radius so we can differentiate plunging and regular orbits
isco = Gradus.isco(m)
# fixed observer velocity
v_obs = @SVector [1.0, 0.0, 0.0, 0.0]
# fixed observer position
# have to put ourselves further away than usual to ensure the metric-flat space
# approximation holds
u_obs = @SVector [0.0, 10_000.0, deg2rad(85), 0.0]

# an interpolation which returns the plunging velocities for us
# (this is part of the standard Gradus.jl exports)
pl_int = interpolate_plunging_velocities(m)

# metric matrix at observer
g_obs = metric(m, u_obs)

redshift = PointFunction((m, gp, mt) -> begin
    let r = gp.u2[2]
        v_disc = if r < isco
            # plunging region
            vtemp = pl_int(r)
            # we have to reverse radial velocity due to backwards tracing convention
            # see
            @SVector [vtemp[1], -vtemp[2], vtemp[3], vtemp[4]]
            # regular circular orbit
            CircularOrbits.fourvelocity(m, r)

        # get metric matrix at position on disc
        g = metric(m, gp.u2)
        @tullio E_disc := - g[i,j] * gp.v2[i] * v_disc[j]
        @tullio E_obs := - g_obs[i,j] * gp.v1[i] * v_obs[j]
        E_obs / E_disc

Note that we position the observer at 10410^4 instead of 103 rg10^3 \, r_\text{g}. We do this to ensure the metric is closer to being flat space, and therefore that Eobs1E_\text{obs} \approx 1 for all geodesics. At 103 rg10^3 \, r_\text{g}, the deviation in EE is in the order of 11 part in 100100 for geodesics directed towards the edges of the accretion disc, but at this alternate distance, the deviation is approximately 1 part in 1000. In the literature, the common practice is to assume Eobs=1E_\text{obs} = 1, and hard-code this result into the redshift calculation, as is done with common Cunningham implementations, but we will instead use the full definition of redshift, precisely as in Eq. (35), and suffer the consequences.

We can then render an image in the usual way:

# set inner radius within the ISCO
d = GeometricThinDisc(
    1.0, 50.0, deg2rad(90)

# small utility function
render_redshift(m, rspf) = rendergeodesics(
    m, u_obs, d, 20_000.0;
    pf = rspf ∘ ConstPointFunctions.filter_intersected,
    abstol=1e-11, reltol=1e-11,

# comparison metric
m_cmp = BoyerLindquistAD(M=m.M, a=m.a)

# render the redshift images
img_ad = render_redshift(m, redshift)
img_cmp = render_redshift(m_cmp, ConstPointFunctions.redshift)

Subtracting these images to compare the difference between our interpolated redshift and the Cunningham analytic values for a=0a=0:


We see the error is in the order of 10410^{-4}, though fortunately this error is relatively robust to integrator tolerances, retaining the order of magnitude even at tolerances of 10810^{-8}.

For first-order method comparison, where m_cmp = BoyerLindquistFO(M=m.M, a=m.a), the situation is a little worse, but not dire:


This is all pretty good news. Sort of. A difference of 10410^{-4} is still about 5 orders of magnitude larger than I would desire. The main difference is clearly within the ISCO, and judging by the features in the noise, is caused by some geometrical anomaly. I will investigate this more at a later date.

I have a strong suspicion the differences in the first-order method outside of the ISCO, which look considerably more like noise, are correlated with numerical differences in how the β\beta impact parameter is mapped, but I would need to test to confirm – which, again, is beyond the scope of this post.

But, with a proof of concept, we just have one missing ingredient before we can calculate this redshift for other metrics!

Determining the ISCO

If we want this interpolation method to work on general static, axis-symmetric methods, the only thing left to define is a function to find the ISCO.

The ISCO is the radius corresponding to

rE=0, \frac{\partial}{\partial r} E = 0,

which may be found relatively easily via Eq. (31). The strategy here is to use automatic differentiation to compute the gradient, and use a root finding tool to optimize finding the minima. With the help of ForwardDiff.jl and Roots.jl, this is pretty easy and fast:

using ForwardDiff
using Roots

function isco(m, lower_bound upper_bound)
    dE(r) = ForwardDiff.derivative(
        x ->, x), 
    d2E(r) = ForwardDiff.derivative(dE, r)

    find_zero((dE, d2E), (lower_bound, upper_bound))

We have to use the bounded solvers, since our domain is undefined below a certain rr, and throws a DomainError. The upper bound may be quite confidently set to something large, e.g. 100 rg100 \, r_\text{g}, however there is generally no safe lower bound we can use that works for all aa. Instead, an approach is determine the lower bound coarsely as a radius for which E/μ>1E / \mu > 1. This may be programmatically found with a simple brute-force:

function find_lower_isco_bound(m; upper_bound = 100.0, step=-0.005)
    # iterate in reverse with a negative step
    for r in upper_bound:step:1.0
        if, r) > 1.0
            return r
    # for type stability
    return 0.0

This function, despite its crudeness, performs fine, and seems to be able to find a suitable solution for all currently implemented metrics:

julia> @btime find_lower_isco_bound($m)
  92.847 μs (0 allocations: 0 bytes)

Adding a dispatch to make our isco API a little friendlier:

isco(m; upper_bound = 100.0) = isco(m, find_lower_isco_bound(m; upper_bound), upper_bound)

Testing this for a few different metric setups:


Seems to work well! We should now be able to visualise redshift patterns for any of our static, axis-symmetric metrics.

Redshifts in a non-GR spacetime

Just for a laugh, let's see what the redshift of an equitorial accretion disc for an Einstein-Maxwell-Dilaton-Axion singularity would look like. The metric components for this I have taken from García et al. (1995)[3], allowing reuse of all of our above code, with just one alteration:

m = DilatonAxionAD(M=1.0, a=0.998, β=1.0, b=1.0)

This metric is a little fiddly to render at the moment, because I probably have some mathematical typos in it, and, since the Schwarzschild equivalent radius is currently not implemented, the cutoff on the inner horizon is underestimated. To avoid uneccessarily long integrations, we can just fix the minimum timestep to something an order of magnitude smaller than our tolerances:

img = @time rendergeodesics(
    m, u_obs, d, 2000.0;
    pf = redshift ∘ ConstPointFunctions.filter_intersected,
    abstol=1e-10, reltol=1e-10,

We set verbose=false to avoid excessive warnings concerning iteration counts, and produce the following render:


Gorgeous! But no-doubt highly unphysical, and probably litered with mathematical typos. There are also clearly numerical anomalies and instabilities visible in the render, which appears to be related to close-to-ISCO redshifts – plently of bug hunting and optimizing to be done.

At some point I will endevour to properly follow the EMDA derivation and motivations, and maybe there's a blog post in that, but this brief tourist's curiousity outside of GR is sufficient for the spectacle, and shows the method of calculating redshift does indeed "work" for static, axis-symmetric metrics.


[1] T Johannsen (2013). Regular black hole metric with three constants of motion. PHYS REV D 88:4.
[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] A García, D Galtsov, and O Kechkin (1995). Class of stationary axisymmetric solutions of the Einstein-Maxwell-Dilaton-Axion field equations. PHYS REV L 74:1276.