How to simulate black hole accretion

Every black-hole image you have seen, including the one behind this page, is made the same way. Pick a spacetime, follow light through it, and add up what that light picks up on its way to the camera. This page walks through those steps in the simplest setting: a non-rotating black hole and a thin disk of glowing gas.

Throughout, $G = c = 1$, and distances are measured in units of the black hole's mass $M$. For M87*, $M$ is about 64 times the Earth–Sun distance.

1Spacetime is the stage

A non-rotating black hole is described by the Schwarzschild metric. It tells us how to measure time and distance at every point around a mass $M$:

$$ds^2 = -\left(1-\frac{2M}{r}\right)dt^2 + \left(1-\frac{2M}{r}\right)^{-1}dr^2 + r^2\left(d\theta^2 + \sin^2\theta\,d\phi^2\right)$$

At $r = 2M$ the first term vanishes: this is the event horizon. Nothing that crosses it, light included, comes back out.

2Light follows geodesics

Light travels along the straightest possible paths in this curved spacetime, called null geodesics. They obey the geodesic equation, together with the condition that light has zero interval:

$$\frac{d^2x^\mu}{d\lambda^2} + \Gamma^\mu_{\alpha\beta}\,\frac{dx^\alpha}{d\lambda}\,\frac{dx^\beta}{d\lambda} = 0, \qquad\qquad g_{\mu\nu}\,\frac{dx^\mu}{d\lambda}\,\frac{dx^\nu}{d\lambda} = 0$$

The Christoffel symbols $\Gamma^\mu_{\alpha\beta}$ encode how the metric bends the path. The spacetime does not change in time and is the same in every direction, so each photon keeps a constant energy $E$ and angular momentum $L$. Their ratio is the impact parameter, $b = L/E$: how far from the black hole the ray would pass if gravity were switched off.

With $u = 1/r$, the whole path of a photon in its orbital plane reduces to one equation:

$$\frac{d^2u}{d\phi^2} + u = 3Mu^2$$

Without the right-hand side, the solutions are straight lines. The term $3Mu^2$ is general relativity: it is small far away and dominant close in.

Impact parameter b5.60 M

Rays arrive from the left, each with its own impact parameter. The black disk is the horizon, $r = 2M$; the dashed circle is the photon sphere, $r = 3M$. Move the slider to follow a single ray. Near $b \approx 5.2M$ it can circle the black hole before escaping or falling in.

The diagram shows the key number in black-hole imaging. Rays with

$$b < b_c = 3\sqrt{3}\,M \approx 5.2\,M$$

fall into the horizon; rays with larger $b$ escape, deflected. Rays very close to $b_c$ wind around the photon sphere at $r = 3M$, possibly several times, before they leave. Seen from far away, the captured rays form the dark shadow, and the winding ones pile up into a thin, bright photon ring at its edge.

3Ray tracing: run the light backwards

Most rays leaving the disk never reach us, so we do not trace light forwards from the gas. We start at the camera instead. Each pixel of the image defines one ray, with coordinates $(\alpha, \beta)$ on the sky and an impact parameter

$$b = \sqrt{\alpha^2 + \beta^2}$$

We integrate the geodesic equation backwards in time from the camera, and stop at the first thing the ray meets:

The horizon. The pixel is black.
The disk. Record the radius $r$ where the ray crosses it, and how fast the gas there is moving.
Infinity. The pixel shows the background sky.

In code, the whole method is a loop over pixels:

for each pixel (alpha, beta):
    b = sqrt(alpha**2 + beta**2)          # impact parameter
    start the ray at the camera, heading back towards the black hole
    repeat:
        take a small step along   d²u/dφ² = 3 M u² − u
        if r < 2M:                 pixel = black;           stop
        if the ray crosses the disk: pixel = intensity(r);  stop
        if r > r_far:              pixel = sky;             stop

4From a hit to a brightness

The light that reaches the camera is not the light the gas emitted. Its frequency is shifted by the ratio

$$g = \frac{\nu_{\rm obs}}{\nu_{\rm em}}$$

Along a ray, the quantity $I_\nu/\nu^3$ is conserved, so the observed intensity follows directly from the emitted one:

$$I_{\rm obs} = g^3\, I_{\rm em}$$

Two effects set $g$. Climbing out of the black hole's gravity costs energy: for gas at rest, $g = \sqrt{1 - 2M/r}$, which falls to zero at the horizon. Gas in the disk also orbits, at angular velocity $\Omega = \sqrt{M/r^3}$, so light sent towards us from the side moving our way is boosted. For a circular orbit,

$$g = \frac{\sqrt{1 - 3M/r}}{1 - \Omega\,\lambda}$$

where $\lambda$ is the photon's angular momentum per unit energy about the black hole's axis, positive when the light leaves in the direction the gas is moving. With $g^3$ applied, one side of the disk outshines the other, the inner edge reddens and dims, and the image takes on its familiar lopsided ring.

5What this leaves out

This is the skeleton of a real calculation. Research codes add the black hole's spin (the Kerr metric), let the gas be thick and partly transparent so intensity builds up along the whole ray (radiative transfer), take the gas from GRMHD simulations, and finally pass the image through a model of the telescope array. Each of those is a page of its own, to come.