Implementing General Relativity: Wormholes, spinning black holes, accretion disks, and redshift

Hiyas! This article is the third in a series on rendering general relativity - if you’re unfamiliar you may want to read these articles first: one, two. We’re going to tie up some loose ends today, and dip heavily into the science behind Interstellar (mostly by accident)¹. This is the last jumbo tutorial article I’m doing in this series - after this we’ll be moving onto numerical relativity, so its time to clear up a few straggler topics:

1. Wormholes
2. A dynamic timestep
3. Workable camera controls/consistently orienting tetrads
4. Observers with velocity
5. Redshift
6. Accretion disks
7. Spinning black holes
8. If there will be more cats (yes)

If you’re looking for a PhD student to work on numerical relativity, then I’d be very interested!

The full code for this article can be found here. First up: Wormholes

The Interstellar wormhole

The paper which describes Interstellar’s wormhole is this one. We want the configurable smooth version, which is described by equations (5a-c)

Given a coordinate system $(t, l, \theta, \phi)$, and the parameters $M$ = mass, $a$ = wormhole length and $p$ = throat radius:

Note that there’s a discontinuity in these equations at $\vert l\vert = a$ as given, so I swap (2) for a $ <= $ instead. Using the raytracer we’ve produced, we can translate this to code:

metric<valuef, 4, 4> get_metric(const tensor<valuef, 4>& position) {
    using namespace single_source;

    valuef M = 0.01;
    valuef p = 1;
    valuef a = 0.001f;

    valuef l = position[1];

    valuef x = 2 * (fabs(l) - a) / (M_PI * M);

    //the analytic differentiator I've got can't handle multiple statements, so this becomes a ternary
    valuef r = ternary(fabs(l) <= a,
                       p,
                       p + M * (x * atan(x) - 0.5f * log(1 + x*x)));


    valuef theta = position[2];

    metric<valuef, 4, 4> m;
    m[0, 0] = -1;
    m[1, 1] = 1;

    m[2, 2] = r*r;
    m[3, 3] = r*r * sin(theta)*sin(theta);

    return m;
}

With almost no changes to our code (other than updating our termination conditions, as $l$ can be $< 0$), we can immediately produce a rendering of this:

While the performance isn’t too bad at ~70ms/frame, its definitely time to fix our timestepping - particularly to cut down on spherical artifacts

Dynamic Timestepping

One of the simplest and most effective strategies for dynamic timestepping is to ensure that the distance² a ray moves is limited to some constant. The simplest version of this method looks something like this:

valuef get_timestep(v4f position, v4f velocity)
{
    v4f avelocity = fabs(velocity);
    return 0.01f/max(max(avelocity.x(), avelocity.y()), max(avelocity.z(), avelocity.w()));
}

The works reasonably well, but it still gives a fairly low universal timestep, whereas we want to concentrate our timestepping more in areas that are likely to be a bit more visually interesting, ie near the object we’re simulating in question. So something like this is a bit more interesting:

valuef get_timestep(v4f position, v4f velocity)
{
    v4f avelocity = fabs(velocity);
    valuef divisor = max(max(avelocity.x(), avelocity.y()), max(avelocity.z(), avelocity.w()));

    valuef low_precision = 0.05f/divisor;
    valuef normal_precision = 0.012f/divisor;
    valuef high_precision = 0.005f/divisor;

    return ternary(fabs(position[1]) < 10, ternary(fabs(position[1]) < 3.f, high_precision, normal_precision), low_precision);
}

This assumes that position[1] is a radial coordinate, which is not always true - and you’ll need a generic system for calculating the distance from your object in question for this to work - but its worth it for the extra performance

With this in place, we get this result, which looks pretty great:

The singularities are barely noticable, and our performance is 60ms/frame @ 1080p. Its not incredible, but interactive framerates here are workable - more advanced timestepping and some general performance improvements can improve things a lot

Watch out for your integrator!

Be aware, not all integrators work with variable timesteps. For example, in this article we’ve previously been using the leapfrog integrator:

        as_ref(velocity) = cvelocity + acceleration * dt;
        as_ref(position) = cposition + velocity.as<valuef>() * dt;

Where the velocity is updated, and then the position is updated with that new velocity. This integrator does not work with a dynamic timestep, and so you must swap to a different scheme - like euler

        as_ref(position) = cposition + cvelocity * dt;
        as_ref(velocity) = cvelocity + acceleration * dt;

Camera Controls / Orienting Tetrads

At the moment, we’re constructing a tetrad directly from the underying metric. This works great, but results in a tetrad that - in spherical metrics - tends to point directly at our object. This leads to very unintuitive camera controls as we move our camera around. In this segment we’re going to commit some crimes against coordinate systems to get them to point roughly where we want them to

What we’d like to do is have our tetrads consistently point in a specific direction - any direction would do. While this isn’t possible to do in the general case (because the tetrads and coordinates are inherently arbitrary), as it turns out - we actually can do this in a pretty wide variety of cases

The idea here is to create an overlaying cartesian coordinate system, and use that to define our new basis vectors. Then you convert these to the real coordinate system, and orthonormalise them to produce a new set of tetrad vectors

First off, lets define what we’d like our basis vectors to be in cartesian coordinates:

Now, our actual coordinate system isn’t in cartesian, so we need to convert these ideal vectors into our actual 4d coordinate system, whatever that may be. For this article everything is in spherical-like coordinates, so lets use that illustratively

We are dealing with tangent vectors, so we need to treat $b$ as a velocity - for that reason we want the total derivative³ of the above set of equations:

Plug $b_k$ into these, and we’ll get a new set of vectors in our actual coordinate system that we’ll call $d_k$⁴. Note that we implicitly created 4-vectors from our 3-vectors here, by setting the time component to $0$

We now have 3 ($_k$) spatial vectors $d_k^\mu$, which have 4 components ($^\mu$) each. Now, calculate your tetrads as per the usual fashion, via gram-schmidt to get a valid tetrad frame $e_i^\mu$, followed by projecting your vectors $d_k^\mu$ into the local frame of reference (via the inverse tetrads). We’ll call these projected local vectors $l_k^\mu$

Note that $i$ ranges over 1-3. We can now orthonormalise (using normal euclidian geometry) the spatial part of these vectors ($\mu$ = 1-3) to produce a new set of basis vectors - which are now roughly oriented in the directions we picked initially ($b_k$) as best as is possible. To finish up, we turn these new 3-vectors back into spacelike 4-vectors by setting the time component to $0$, and then convert them into global coordinates with our old tetrads. These form the 3 new spacelike tetrad vectors

Only one vector is ‘true’

As you may have spotted, orthornormalising these vectors changes them - orthonormalising can only preserve the direction of a single vector (the one we start from). The correct vector to preserve is the ‘up’ vector that you use for your fixed mouse vertical axis - this means that the camera doesn’t roll as you move it around - which is very disorienting. The new vectors are still a bit of a compromise - we’re discarding the timelike component, which means that they’ll still inevitably be a bit ‘incorrect’ as to where they point - but in practice this scheme works very well

One other thing to note is that we only orient the initial tetrad if we’re parallel transporting tetrads⁵

Code

if(should_orient)
{
    v4f spher = GenericToSpherical(position.get());
    v3f cart = spherical_to_cartesian(spher.yzw());

    v3f bx = (v3f){1, 0, 0};
    v3f by = (v3f){0, 1, 0};
    v3f bz = (v3f){0, 0, 1};

    v3f sx = convert_velocity(cartesian_to_spherical, cart, bx);
    v3f sy = convert_velocity(cartesian_to_spherical, cart, by);
    v3f sz = convert_velocity(cartesian_to_spherical, cart, bz);

    //there's an extra step here through spherical coordinates, to more consistently orient the basis when r < 0. This is optional
    sx.x() = ternary(spher.y() < 0, -sx.x(), sx.x());
    sy.x() = ternary(spher.y() < 0, -sy.x(), sy.x());
    sz.x() = ternary(spher.y() < 0, -sz.x(), sz.x());

    v4f dx = convert_velocity(SphericalToGeneric, spher, (v4f){0.f, sx.x(), sx.y(), sx.z()});
    v4f dy = convert_velocity(SphericalToGeneric, spher, (v4f){0.f, sy.x(), sy.y(), sy.z()});
    v4f dz = convert_velocity(SphericalToGeneric, spher, (v4f){0.f, sz.x(), sz.y(), sz.z()});

    //invert the tetrads we get from the gram-schmidt process, calculated in the usual fashion
    inverse_tetrad itet = tet.invert();

    v4f lx = itet.into_frame_of_reference(dx);
    v4f ly = itet.into_frame_of_reference(dy);
    v4f lz = itet.into_frame_of_reference(dz);

    //start normalising at the y coordinate, which I define as 'up' arbitrarily
    std::array<v3f, 3> ortho = orthonormalise(ly.yzw(), lx.yzw(), lz.yzw());

    //undo the permutation above
    v4f x_basis = {0, ortho[1].x(), ortho[1].y(), ortho[1].z()};
    v4f y_basis = {0, ortho[0].x(), ortho[0].y(), ortho[0].z()};
    v4f z_basis = {0, ortho[2].x(), ortho[2].y(), ortho[2].z()};

    v4f x_out = tet.into_coordinate_space(x_basis);
    v4f y_out = tet.into_coordinate_space(y_basis);
    v4f z_out = tet.into_coordinate_space(z_basis);

    tet.v[1] = x_out;
    tet.v[2] = y_out;
    tet.v[3] = z_out;
}

Results

This results in your metric looking like this as you fly around, making it way easier to implement camera controls, as now you have a consistent usable tetrad basis

Observers with velocity / Lorentz boosts

The interstellar wormhole has neglegible gravity, so we’ll simply sit still forever if we can’t represent a moving observer. Ideally we’d like to give our observer a push in a particular direction, to represent them moving. Note that we define this velocity relative to our initial frame of reference

In a previous article, we learnt that given a set of tetrad vectors $e_i^\mu$, the 0th vector $e_0^\mu$ represents the velocity of our observer. If we want to represent an observer with a different speed, sadly we can’t just modify that component - to transform the entire tetrad is a little more involved. The standard method for this is something called a Lorentz boost (or more generally, a Lorentz transform - which may include rotations), which relates two observers moving at different speeds. A Lorentz transform in general relativity is often denoted by the symbol $\Lambda^i_{\;j}$, or $B^i_{\;j}$ for a Lorentz boost specifically here. This paper and this paper contain more information if you want to go digging

For us, we’re looking to perform a Lorentz boost in an arbitrary direction, and apply it to our basis vectors. Before we get there, we need to know what a 4-velocity is - something we’ve skimmed over a bit

4-velocities are not quite like 3-velocities

Lets imagine you have a regular good ol’ velocity in 3 dimensional space. We’ll scale our velocities so that $1$ represents the speed of light, and $0$ represents stationary relative to our observer. We’re going to examine the difference between:

1. 3-velocities parameterised by coordinate time $\frac{dx^i}{dt}$, your regular everyday concept of velocity. We will call this $v$, and its euclidian magnitude is $\vert v \vert$. $\vert v \vert = 1$ represents light
2. 4-velocities parameterised by coordinate time⁶ $\frac{dx^\mu}{dt}$
3. 4-velocities parameterised by an affine parameter $\frac{dx^\mu}{d\lambda}$
4. 4-velocities parameterised by proper time⁷ $\frac{dx^\mu}{d\tau}$
5. Timelike vs lightlike geodesics, with all of the different parameterisations

It is common to define 4-velocities as only being those velocities which are parameterised by proper time, but we’re going with the generalised understanding here. Note that this segment is dealing with line elements in minkowski as that’s where we’ll be constructing our geodesics, and as such is special relativity. That said, the language of special relativity, and general relativity are often not very similar at all, so I’m putting this all down in our terminology

Lightlike Geodesics

A lightlike geodesic is defined has having $ds^2 = 0$. If you’ve forgotten, remember that we evaluate this quantity with the metric tensor, or directly with the line element

Coordinate time parameterisation

Constructing a 4-velocity for a lightlike geodesic which is parameterised by coordinate time is easy. Definitionally for a lightlike geodesic:

We can use the line element for minkowski as such, plugging in our 3-velocity to solve for the time component of a 4-velocity

We know a ray of light moves with a speed of $\vert v \vert=1$, therefore

Giving us a final coordinate parameterised 4-velocity of $\frac{dx^\mu}{dt} = (\pm 1, v^0, v^1, v^2)$

Affine parameterisation

This is the most common parameterisation for a geodesic, and the one we will be using. Luckily, because the parameter for a lightlike geodesic has no particularly useful interpretation⁸, we simply set $\lambda = t$ at the moment of construction

In minkowski, because spacetime is trivially flat, this parameterisation will always hold. In curved spacetime, this relation will only hold at the moment of construction, and then will diverge. This is because we use different geodesic equations depending on the parameterisation we pick

Proper time parameterisation

The full line element reads:

For a lightlike geodesic, $ds^2 = 0$, so $d\tau^2 = 0$ ($d\tau$ being the change in proper time). There is therefore⁹ no proper time parameterisation of a lightlike geodesic. For this reason, it is common to state that the velocity of a lightlike geodesic is not a 4-velocity, as it can never be parameterised by proper time

Timelike Geodesics

A timelike geodesic is defined as when $ds^2 < 0$¹⁰

Coordinate time parameterisation

Lets first up construct a timelike 4-velocity $\frac{dx^\mu}{dt}$ parameterised by coordinate time. We know that the $dx^0/dt$ component must be $1$, as $dx^0 = dt$¹¹. Lets check if $(1, v^0, v^1, v^2)$ is timelike:

Therefore, if $\vert v \vert < 1$, we do get a valid timelike geodesic $\frac{dx^\mu}{dt} = (1, v^0, v^1, v^2)$ parameterised by coordinate time

Proper time parameterisation

To transform from a coordinate time paramterisation, to a proper time parameterisation, we must multiply $\frac{dx^\mu}{dt}$ by $\frac{dt}{d\tau}$. How do we calculate $d\tau$? Easy, we use the line element for minkowski, and plug in our coordinate time velocity

You may recognise this as the formula for the Lorentz factor $\gamma = \frac{dt}{d\tau}$. So to construct a 4-velocity parameterised by coordinate time from the 3-velocity $v^i$, we do:

Part of the reason why I’m spelling this out so explicitly is because all this notation is regularly used to mean the same thing, so hopefully you can come back to this in the future

Affine parameterisation

Like with lightlike geodesics, we can construct ‘an’ affine parameterisation by setting $\lambda = t$ at the moment of construction, after which the two parameters diverge. This however is very uncommon, and is only mentioned for completeness. When you do this, the parameterisation is hard to interpret physically¹²

We can also construct an affine time parameterisation by setting $\lambda = \tau$, where $d\tau = ds^2 = -1$. One very neat fact of proper time is that it is a general affine parameterisation, and so if we use a proper time parameterised geodesic and plug it through the geodesic equation specialised for the affine parameter (which is the one we use), it remains parameterised by proper time

Cat break

Calculating a Lorentz boost

We now know how to make an observer with a (timelike) 4-velocity in minkowski, by constructing it from a 3-velocity. To be very explicit, given a velocity in cartesian coordinates $d^i = (dx, dy, dz)$, where $\vert d\vert < 1$

We then want to convert $v_{local}$ to $v$ by transforming it with the tetrad, as $v^\mu = e_i^\mu v^i_{local}$, to get our new observer velocity in our curved spacetime

If we have an initial 4-velocity $u^\mu = e_0^\mu$ of our tetrad, and we want to boost the tetrad to represent an observer with a 4-velocity $v$, the formula for the Lorentz boost is this¹³:

$\delta^i_{\;j}$ is known as the kronecker delta. It is $1$ when $i==j$, and $0$ otherwise

Next up, we need to apply this Lorentz boost to our tetrad vectors. Lets say our original tetrads are $e_i^\mu$, and our boosted tetrads are $\hat{e}_i^\mu$. We already know what the new value of $\hat{e}_0^\mu$ will be, as it must be $v$, but the general transform is:

Code

//gets our new observer
v4f get_timelike_vector(v3f velocity, const tetrad& tetrads)
{
    v4f coordinate_time = {1, velocity.x(), velocity.y(), velocity.z()};

    valuef lorentz_factor = 1/sqrt(1 - (velocity.x() * velocity.x() + velocity.y() * velocity.y() + velocity.z() * velocity.z()));

    v4f proper_time = lorentz_factor * coordinate_time;

    ///put into curved spacetime
    return proper_time.x() * tetrads.v[0] + proper_time.y() * tetrads.v[1] + proper_time.z() * tetrads.v[2] + proper_time.w() * tetrads.v[3];
}

//does the actual tetrad boosting
tetrad boost_tetrad(v3f velocity, const tetrad& tetrads, const metric<valuef, 4, 4>& m)
{
    using namespace single_source;

    v4f u = tetrads.v[0];
    v4f v = get_timelike_vector(velocity, tetrads);

    v4f u_l = m.lower(u);
    v4f v_l = m.lower(v);

    valuef Y = -dot(v_l, u);

    ///https://arxiv.org/pdf/2404.05744 18
    tensor<valuef, 4, 4> B;

    for(int i=0; i < 4; i++)
    {
        for(int j=0; j < 4; j++)
        {
            valuef kronecker = (i == j) ? 1 : 0;

            B[i, j] = kronecker + ((v[i] + u[i]) * (v_l[j] + u_l[j]) / (1 + Y)) - 2 * v[i] * u_l[j];
        }
    }

    tetrad next;

    for(int a=0; a < 4; a++)
    {
        for(int i=0; i < 4; i++)
        {
            valuef sum = 0;

            for(int j=0; j < 4; j++)
            {
                sum += B[i, j] * tetrads.v[a][j];
            }

            next.v[a][i] = sum;
        }
    }

    return next;
}

Wormhole Results

With a wormhole with params M = 0.01, p = 1, a = 1, and an observer moving at 0.2c towards the wormhole:

Of course the foreground and background texture do not have to be the same, and can be distinguished by $l < 0$

M = 0.01, p = 1, a = 0.001

Redshift

The equations for redshift in general relativity are fairly straightforward. First off, we define the redshift¹⁴ $z$, as follows:

This is commonly also written:

$k^\mu$ represents our geodesic’s velocity, and $u$ is the observer’s velocity. $u^\mu_{emit}$ specifically will be where our ray terminates (and may either be the $e_0$ of a calculated tetrad, or in the case of eg an accretion disk - a known velocity), and $u^\mu_{obs} = e_0$ is our initial observer’s velocity (after boosting!). Do note that the metric tensors are evaluated at different locations

Next up, we need to work out how our light changes, from our end frame of reference (defined by $u^\mu_{emit}$), to our initial frame of reference. Light has two properties - frequency/wavelength, and intensity. The equation for wavelength looks like this¹⁵

Once we have our new wavelength, we need to calculate the intensity. We can do this by calculating the Lorentz invariant (constant in every frame of reference): $\frac{I_\nu}{\nu^3}$, where $\nu$ is your frequency, and $I_\nu$ is spectral radiance. See here 1.26 for details. Note that the quantity $I_\nu \lambda^3$ is also therefore Lorentz invariant

So, to calculate our observed intensity, we know that this quantity never changes, and can say:

This equation is linear in terms of intensity, which is often a useful property

Once we have our new wavelength(s) $\lambda_{obs}$, and $I_{obs}$, in theory we have everything we need to render our our final colour. We just have two unknowns, which are our initial intensity, and the initial wavelength(s)

Spectral Radiance

The $I_\nu$ that we’re dealing with in these equations is called spectral radiance - it is the power per unit frequency, and is not really a quantity that you often have in practice. Terminology here is confusing:

If you want the equation for transforming a radiant flux, you’re looking for:¹⁶

This is often much easier to work with, and is the quantity we get out of our accretion disk model later

Where do $I_{emit}$ and $\lambda_{emit}$ come from?

It depends what we’re simulating. For our use case - redshifting a galaxy background, you’d need frequency and intensity data across the entire sky. A good starting point is over here, luckily we live in 2024 and a significant amount of this information is simply public - unfortunately these skymaps do not come with what units their intensity data is in, making them unusable¹⁷. Still, you can go find the original surveys - though it requires significant digging which I’m not going to do in this article

If you have a blackbody radiator, it becomes fairly straightforward, as given a temperature we can redshift that directly, via the equation:

For our galaxy background, we’re instead going to implement illustrative redshift, rather than attempting to find calibration data for hundreds of surveys - which I’ve had to try very hard to stop myself from doing

Illustrative redshift

Given that accurate skymap data is hard to find, I’m going to discard physicality when it comes to the colour, and just show a measure of redshift. We can still calculate $z+1$ as per normal, and then transform our radiant flux without needing to have a specific wavelength. The only thing we actually need wavelengths for is generating a final colour, but instead we’ll essentially just interpolate between red and blue with our $z$ value. When reading from a texture, our intensity data can be defined as the $Y$ component of the $XYZ$ colour space which represents energy, see here

Where $(r, g, b)$ is a linear sRGB colour. Once we’ve calculated our new intensity via the redshift equation, its then straightforward to perform a recolouring, and I’ll show an example below. $z$ has a range of $[-1, +\infty]$, so its easiest to split into two branches

Redshift Only (z > 0)

new_colour = mix(old_colour, pure_red / 0.2126f, tanh(z));

Redshift naturally fades to black as the intensity drops. The choice of tanh to map the infinite range to $[0, 1]$ is fairly arbitrary. See below for the explanation of the constant

Blueshift Only (z < 0)

iv1pz = (1/(1 + z)) - 1;
interpolating_fraction = tanh(iv1pz);
new_colour = mix(old_colour, pure_blue / 0.0722f, interpolating_fraction);

The mapping here is more complicated to replicate the same falloff as redshift. One question you might have is where these division constants come from: notice that they’re the same constants we use to calculate $Y$. This ensures that our new colour scales in brightness better - a pure white background colour will map to an equivalently bright blue colour, and darker colours will blend towards a brighter colour as our intensity increases. Without the division, a bright background colour would first darken as we mix in the visually less bright primary, before brightening again - this U bend in brightness is visually strange

One problem specific to blueshift is that our energy is unbounded, and our pixels can become infinitely bright. Its therefore much more aesthetically pleasing to spill over the extra energy into white once we max out the blue colour - which is shown in the full code sample

Code

This is one of those things that’s a lot easier to express in code, rather than equations:

//calculate Y of XYZ
valuef energy_of(v3f v)
{
    return v.x()*0.2125f + v.y()*0.7154f + v.z()*0.0721f;
}

v3f redshift(v3f v, valuef z)
{
    using namespace single_source;

    {
        valuef iemit = energy_of(v);
        valuef iobs = iemit / pow(z+1, 4.f);

        v = (iobs / iemit) * v;

        pin(v);
    }

    valuef radiant_energy = energy_of(v);

    v3f red = {1/0.2125f, 0.f, 0.f};
    v3f green = {0, 1/0.7154, 0.f};
    v3f blue = {0.f, 0.f, 1/0.0721};

    mut_v3f result = declare_mut_e((v3f){0,0,0});

    if_e(z >= 0, [&]{
        as_ref(result) = mix(v, radiant_energy * red, tanh(z));
    });

    if_e(z < 0, [&]{
        valuef iv1pz = (1/(1 + z)) - 1;

        valuef interpolating_fraction = tanh(iv1pz);

        v3f col = mix(v, radiant_energy * blue, interpolating_fraction);

        //calculate spilling into white
        {
            valuef final_energy = energy_of(clamp(col, 0.f, 1.f));
            valuef real_energy = energy_of(col);

            valuef remaining_energy = real_energy - final_energy;

            col.x() += remaining_energy * red.x();
            col.y() += remaining_energy * green.y();
        }

        as_ref(result) = col;
    });

    as_ref(result) = clamp(result, 0.f, 1.f);

    return declare_e(result);
}

For a Schwarzschild black hole, a velocity of 0.5c boosting towards the black hole, and a starting position of $(0, 5, 0, 0)$, you end up with this

From the front and back. Notice that we blueshift in the direction of travel, and see a redshift in the opposite direction. Do be aware that you should do all of this in linear colour, and only map your colours to sRGB right at the very end of your pipeline

Spinning black holes / The Kerr Metric

Accretion Disks

No black hole is complete without an accretion disk, so today we’re going to be looking at implementing a reasonable accretion disk model. We’re going to use whats called a thin-disk model (98-100), because it has a straightforward algebraic solution. In the future, we’ll be simulating accretion disks directly

One thing we should examine first is the concept of an innermost stable orbit, or ISCO for short

ISCO

The ISCO is the radius at which orbits around a black hole become unstable. Outside of this radius we can have long lived circular orbits - inside of this radius, circular orbits do not exist. This inner region is known as the ‘plunging’ region, and anything that’s within this region rapidly enters into a black hole

For a schwarzschild black hole, the ISCO is defined as $r_{isco} = 3rs = 6M$. For a Kerr (spinning) black hole, the ISCO is defined as follows¹⁸ for a black hole with mass $M$, and spin $a$:

The retrograde ISCO is higher than the prograde ISCO, so flip the sign appropriately when calculating

Accretion Disk Regions

Because orbits within the ISCO are unstable, matter depletes from this region very quickly. For this reason, accretion disks are often modelled as having a gap between the event horizon, and the ISCO - which we will follow¹⁹. For the model we’re looking at, here equations 98-100, we have three regions for us to work with, which are called:

1. The inner region
2. The middle region
3. The outer region

Distinguishing between these three regions is done by two values

1. Gas pressure vs radiation pressure. When gas pressure < radiation pressure, we’re either in the inner or outer region
2. Whether opacity is driven by free-free interactions, or electron scattering. If free-free > electron scattering, we’re in the outer region, otherwise we’re in the inner or middle region

The details of this are interesting²⁰, but we’re going to focus on the implementation. The general idea for us is to iterate from the innermost boundary of our accretion disk ($r = r_{isco}$), and terminate at some finite radius away from the accretion disk ($r = 100M$), working out which region we’re in as we go. We know at the start, we must be in region 1 - the inner region - at the ISCO

To distinguish when we transition from region 1, to region 2, which want to calculate (gas pressure / radiation pressure), and if its > 1 move into the middle region. This quantity is labelled $\beta / (1-\beta)$ (I do not know why)

To distinguish when we transition from region 2, to region 3, we calculate the quantity $T_{ff} / T_{es}$, which is the free-free opacity / the electron scattering opacity. When this quantity is > 1, we swap to region 3²¹

Other details

To be able to implement this, we first need to calculate our spatial functions $A$ $B$ $C$ $D$ $E$ $Q_0$ $Q$, as well as $y_1$ $y_2$ $y_3$, and $y_0$. While not mentioned on the page we’re looking at, subscript $_0$ variables are evaluated at the ISCO. We do have an explicit formula for $Q_0$, which means that the only other one we need is $y_0 = (\frac{r_{isco}}{M})^{\frac{1}{2}}$

The meaning of the variables here are:

With this, we should have everything²² we need to implement this correctly

Code

The equations here - as written, are long and complicated to implement correctly. I won’t reproduce the equations algebraically in full here - it just introduces a risk of mistakes. You can find the full code for implementing this method, and producing a nice accretion disk texture over here

The core of our algorithm looks like this:

for(int steps = 0; steps < max_steps; steps++)
{
    double r = mix(isco, outer_boundary, steps / (double)max_steps);
    double x = sqrt(r/mass);

    double r_star = r / mass;

    double x_pow_m2 = mass/r;
    double x_pow_m4 = pow(mass/r, 2.);

    double A = 1 + assq * x_pow_m4 + 2 * assq * pow(x, -6.);
    double B = 1 + a_star * pow(x, -3.);
    double C = 1 - 3 * x_pow_m2 + 2 * assq * pow(x, -3.);
    double D = 1 - 2 * x_pow_m2 + assq * pow(x, -4.);
    double E = 1 + 4 * assq * x_pow_m4 - 4 * assq * pow(x,-6.) + 3 * pow(a_star, 4.) * pow(x, -8.);
    double Q = B * pow(C, -1/2.) * (1/x) * (x - x0 - (3/2.) * a_star * log(x/x0)
                                            - (3 * pow(x1 - a_star, 2.) / (x1 * (x1 - x2) * (x1 - x3))) * log((x - x1) / (x0 - x1))
                                            - (3 * pow(x2 - a_star, 2.) / (x2 * (x2 - x1) * (x2 - x3))) * log((x - x2) / (x0 - x2))
                                            - (3 * pow(x3 - a_star, 2.) / (x3 * (x3 - x1) * (x3 - x2))) * log((x - x3) / (x0 - x3)));

    if(region == region_type::INNER)
    {
        ///This is B/(1-B)
        double p_gas_p_rad = 4 * pow(10., -6.) * pow(alpha, -1/4.) * pow(m_star, -1/4.) * pow(mdot_star, -2.) * pow(r_star, 21/8.) * pow(A, -5/2.) * pow(B, 9/2.) * D * pow(E, 5/4.) * pow(Q, -2.);

        ///in the edge region, gas pressure dominates over radiation pressure
        ///in the inner region, gas pressure is less than radiation pressure
        ///in the middle region, gas pressure is greater than radiation pressure
        if(p_gas_p_rad > 1)
            region = region_type::MIDDLE;
    }

    if(region == region_type::MIDDLE)
    {
        //in the outer region opacity is free-free
        //in the middle region, opacity is electron scattering
        double Tff_Tes = (2 * pow(10., -6.)) * (pow(mdot_star, -1.)) * pow(r_star, 3/2.) * pow(A, -1.) * pow(B, 2.) * pow(D, 1/2.) * pow(E, 1/2.) * pow(Q, -1.);

        if(Tff_Tes >= 1)
            region = region_type::OUTER;
    }

    double surface_flux = 7 * pow(10., 26.) * pow(m_star, -1.) * mdot_star * pow(r_star, -3.) * pow(B, -1.) * pow(C, -1/2.) * Q;
    double T = 0;

    if(region == region_type::INNER)
        T = 5 * pow(10., 7.) * pow(alpha, -1/4.) * pow(m_star, -1/4.) * pow(r_star, -3/8.) * pow(A, -1/2.) * pow(B, 1/2.) * pow(E, 1/4.);

    ///'edge' region as well in the more developed model
    if(region == region_type::MIDDLE)
        T = 7 * pow(10., 8.) * pow(alpha, -1/5.) * pow(m_star, -1/5.) * pow(mdot_star, 2/5.) * pow(r_star, -9/10.) * pow(B, -2/5.) * pow(D, -1/5.) * pow(Q, 2/5.);

    if(region == region_type::OUTER)
        T = 2 * pow(10., 8.) * pow(alpha, -1/5.) * pow(m_star, -1/5.) * pow(mdot_star, 3/10.) * pow(r_star, -3/4.) * pow(A, -1/10.) * pow(B, -1/5.) * pow(D, -3/20.) * pow(E, 1/20.) * pow(Q, 3/10.);

    radius.push_back(r);
    brightness.push_back(surface_flux);
    temperature.push_back(T);
}

This gives pretty nice results. For $M=1$²³, $a=0.6$, $\dot{m} = 0.3$:

In practice you don’t need to use a 2d texture of an accretion disk - because its spherically symmetric, you could just take a radial slice

Colouring

One of the values we get out of our accretion disk model is a temperature $T$. If we assume that the accretion disk is a blackbody radiator, there’s a 1:1 mapping between every temperature, and the human perception of this colour. This forms a line across our colour space, and this relation from temperature to colour is known as the plankian locus. Working this out directly is nontrivial, though wikipedia also provides a useful approximation here. I’ve provided both of these methods over here, as I’d already implemented most of this in the past²⁴

If you run the accretion disk temperatures through this process, you’ll be surprised to learn the colour of a real accretion disk:

(M = 1, a=0.99)

Rendering

Rendering the accretion disk is fairly straightforward. Because its modelled as a thin disk, you can check if a geodesic crosses the equatorial plane ($\theta = n\;\pi + \pi/2$), and sample the accretion disk texture. Here, I approximate the opacity as the brightness of the texture, because we’re not physically accurately raytracing this²⁵, and do some very basic volumetrics on the disk

One key thing to note here is that we’re going to implement redshift, and to do that, we need to know the velocity of the fluid in the accretion disk. We already know that its moving approximately in circular orbits (otherwise it’d escape, or hit the black hole), but what is that as a velocity vector specifically?

Well, in the equatorial plane for circular orbits, we know that $dr = 0$, and $d\theta = 0$. This means that we’re only left to solve for $dt$, and $d\phi$.

Luckily other people have put in the legwork here, and here, and figured out that the angular velocity of a geodesic in the Kerr spacetime is

Given that we know that our geodesic must point in the $d\phi$ direction, we can work out the $d\phi$ component as being $wr$. This means that $dt = d\phi/w = r$. We now have the 4-velocity of a geodesic on our equatorial orbit

Physically accurate redshifting

The nice thing here is that with the approximation that an accretion disk is a blackbody radiator, we can directly redshift the temperature²⁶, and so end up with an accurate visual colour out of the other end:

Given that mapping a temperature to a colour is an expensive process, I simply chuck the whole thing into a big buffer indexed by temperature, and call it a day

Code

#ifdef HAS_ACCRETION_DISK
valuef period_start = floor(position.z() / pi) * pi;

valuef in_start = cposition.z() - period_start;
valuef in_end = position.z() - period_start;

valuef min_start = min(in_start, in_end);
valuef max_start = max(in_start, in_end);

///we've hit the equatorial plane
if_e(pi/2 >= min_start && pi/2 <= max_start, [&]
{
    valuef radial = position[1];

    valuef M = BH_MASS;
    valuef a = BH_SPIN;

    ///calculate the angular velocity
    valuef w = pow(M, 1.f/2.f) / (pow(radial, 3.f/2.f) + a * pow(M, 1.f/2.f));

    ///calculate the real observer velocity
    v4f observer = {radial, 0, 0, w * radial};

    valuef ds = dot_metric(observer, observer, get_metric(cposition));

    ///valid timelike circular geodesic on the +r branch of kerr
    if_e(ds < 0 && radial > 0, [&]
    {
        int buffer_size = 2048;
        valuef outer_boundary = 2 * BH_MASS * 50;

        valuef buffer_coordinate = (fabs(radial) / outer_boundary) * buffer_size;
        v3f disk = lookup(accretion_disk, buffer_coordinate, 0.f, valuef(buffer_size - 1));

        //if we've already hit another element of the disk, dim our contribution
        disk = disk * clamp(1 - declare_e(opacity), 0.f, 1.f);
        //add a fairly arbitrary amount of opacity to the tracing based on the brightness
        as_ref(opacity) = declare_e(opacity) + energy_of(disk) * 50;

        ///change the parameterisation to proper time
        observer = observer / sqrt(fabs(ds));

        pin(observer);

        #define ACCURATE_REDSHIFT
        #ifdef ACCURATE_REDSHIFT
        valuef temperature_in = lookup(temperature, buffer_coordinate, 0.f, valuef(buffer_size - 1));

        ///temperature == 0 is impossible in our disk, so indicates an invalid area
        if_e(temperature_in >= 1, [&]
        {
            valuef zp1 = get_zp1(g.position, g.velocity, initial_observer, cposition, cvelocity, observer, get_metric);

            ///https://www.jb.man.ac.uk/distance/frontiers/cmb/node7.htm
            valuef shifted_temperature = temperature_in / zp1;

            valuef old_brightness = energy_of(disk);

            ///https://arxiv.org/pdf/gr-qc/9505010 12
            valuef new_brightness = old_brightness / pow(zp1, 4.f);

            ///convert our temperature to a colour
            v3f final_colour = lookup(bbody_table, shifted_temperature, 1.f, 100000 - 1.f) * new_brightness;

            as_ref(colour_out) = declare_e(colour_out) + final_colour;
        });
        #endif

        #ifdef RAW_DISK
        as_ref(colour_out) = declare_e(colour_out) + disk;
        #endif // RAW_DISK

        if_e(opacity >= 1, [] {
            break_e();
        });
    });
});
#endif

Interstellar’s other secret wormhole / Spinning black holes

A spinning black hole is generally modelled by the Kerr (or Kerr-Newman, with charge) metric. The interior of these metrics are actually pretty interesting. In addition to a ringularity (a ring singularity), the centre of a Kerr style black hole contains a wormhole, and copious time travel - which are certainly unusal things to find²⁷. We’re going to take a trip inside!

The metric for a Kerr-Newman black hole in ingoing²⁸ coordinates looks like this:

where

$Q$ is our charge, $M$ is our mass parameter, and $a$ is the black hole’s spin

Kerr Metric Code

#define BH_MASS 1
#define BH_SPIN 0.99

metric<valuef, 4, 4> get_metric(const tensor<valuef, 4>& position) {
    using namespace single_source;

    valuef r = position[1];
    valuef theta = position[2];

    metric<valuef, 4, 4> m;

    float M = BH_MASS;
    float a = BH_SPIN;

    valuef rs = 2 * M;

    valuef ct = cos(theta);
    valuef st = sin(theta);

    valuef R2 = r*r + a * a * ct * ct;
    valuef D = r*r + a * a - rs * r;

    valuef dv = -(1 - (rs * r) / R2);
    valuef dv_dr = 2;
    valuef dv_dphi = -(2 * a * st * st / R2) * (rs * r);
    valuef dr_dphi = -2 * a * st * st;
    valuef dtheta = R2;
    valuef dphi = -(st * st / R2) * (D * a * a * st * st - pow(a * a + r*r, 2.f));

    ///v, r, theta, phi
    m[0, 0] = dv;
    m[1, 0] = 0.5f * dv_dr;
    m[0, 1] = 0.5f * dv_dr;

    m[3, 0] = 0.5f * dv_dphi;
    m[0, 3] = 0.5f * dv_dphi;

    m[1, 3] = 0.5f * dr_dphi;
    m[3, 1] = 0.5f * dr_dphi;

    m[2, 2] = dtheta;
    m[3, 3] = dphi;

    return m;
}

This gives us a pretty nice looking black hole, which looks like this when rendered with an accretion disk (including redshift):

Parallel transport numerical accuracy

One thing to note is that the interior of Kerr is very numerically unstable, due to very high accelerations on the equatorial plane - it behaves very poorly in general. To make this article work, I had to upgrade the parallel transport code to be second order, to increase the accuracy sufficiently. This is a very basic second order integrator with nothing fancy otherwise:

v4f transport2(v4f what, v4f position, v4f next_position, v4f velocity, v4f next_velocity, valuef dt, auto&& get_metric)
{
    using namespace single_source;

    tensor<valuef, 4, 4, 4> christoff2 = calculate_christoff2(position, get_metric);

    pin(christoff2);

    v4f f_x = parallel_transport_get_change(what, velocity, christoff2);

    v4f intermediate_next = what + f_x * dt;

    tensor<valuef, 4, 4, 4> nchristoff2 = calculate_christoff2(next_position, get_metric);

    pin(nchristoff2);

    return what + 0.5f * dt * (f_x + parallel_transport_get_change(intermediate_next, next_velocity, nchristoff2));
}

Taking a trip into a spinning black hole

The ringularity in Kerr is extremely cool, and is one of my favourite things to simulate, so lets have a look at how this looks with an accretion disk!

Purrfect!

While I partly built this for the memes, its worth noting that the interior of a black hole and the wormhole itself is generally considered to be a mathematical artifact, and that it wouldn’t show up in reality. With that in mind, nobody really has a clue whats inside a black hole²⁹!

The end

That’s the end of this small series of rendering arbitrary analytic metric tensors - there’s a lot more we could do here and I’ve been really crunching these articles down to the bones (possibly a bit too much). There’ll be articles in the future on this topic, but they’ll likely be on specific elements of this topic rather than the more broad approach I’ve taken here. The focus after this is largely going to be on numerical relativity - simulating the evolution of spacetime itself, and how to do that on a GPU

With that in mind, here’s the current list of topics that I’m going to get to in the future:

1. Examining the BSSN equations and the ADM formalism in numerical relativity, to evolve some simple numerical spacetimes
2. Binary black hole mergers
3. Fast laplacian solving
4. Gravitational wave extraction
5. Approximate numerical raytracing
6. Full numerical raytracing
7. Stability modifications, and constraint damping
8. Integrators for the EFEs
9. Examining and calculating mass in GR
10. Boundary conditions
11. Particle system dynamics, including massless and massive accretion disks
12. Investigating the hoop conjecture, mainly because its cool
13. Negative energy in general relativity
14. Relativistic eularian hydrodynamics
15. Neutron star initial conditions, and equations of state
16. Relativistic smooth particle hydrodynamics, which is a fun phrase to say
17. Relativistic magnetohydrodynamics
18. CCZ4/Z4cc/constraint damped formalisms, and friends
19. Solving elliptic equations

On top of this, there are a number of misc topics to cover in analytic metrics:

1. 4d Triangle rasterisation
2. Time travel, and CTCs
3. More interesting metric tensors, like the misner spacetime, cosmic strings, or analytic binary black hole metrics
4. Near-massless (dust) particle systems
5. Near-massless (dust) accretion disks
6. Wave propagation in curved spacetimes
7. 256x anisotropic texture filtering
8. Relativistic volumetrics

This is going to take a hot minute to write up in a tutorially format - though its all largely work I’ve done previously (except smooth particle hydrodynamics, and magnetohydrodynamics). This kind of information doesn’t really exist outside of papers - and there are no gpu implementations - so it seems good to write it up

Until then, good luck!