Tracing rays…

Black Hole RU

Black Hole

This isn’t a drawing. For every pixel, your computer is solving Einstein’s equations right now: where the light came from as it swung around the spinning hole, and what became of it along the way.

Drag to fly around it. The article with the formulas is below.

How it works

The picture above is neither a drawing nor a video. For every pixel, your browser solves the equations of motion of light in the spacetime of a spinning black hole, and works out the color from the laws of radiation. Below: which formulas are at work, what has been simplified, and how we checked that it all adds up.

Kerr spacetime

A spinning, uncharged black hole is described by Kerr’s solution (1963). In units where the gravitational constant, the speed of light and the mass of the hole are all equal to one ($G = c = M = 1$), it has a single parameter — the spin $a = J/M$, which runs from 0 to 1. From here on, all distances are measured in gravitational radii $GM/c^2$: for Sagittarius A* that’s 6.3 million km, for a hole of 10 solar masses, 15 km.

We compute in Cartesian Kerr–Schild coordinates, in which the metric is flat space plus a “correction” along a null vector $l_\mu$:

$$g_{\mu\nu} = \eta_{\mu\nu} + f\, l_\mu l_\nu,\qquad f = \frac{2 r^3}{r^4 + a^2 z^2},$$ $$l_\mu = \Big(1,\ \frac{r x + a y}{r^2 + a^2},\ \frac{r y - a x}{r^2 + a^2},\ \frac{z}{r}\Big),$$

where $r$ is defined implicitly by $r^4 - (x^2 + y^2 + z^2 - a^2)\, r^2 - a^2 z^2 = 0$. The great virtue of these coordinates is that they have no singularity at the horizon. So the very same program follows light outside the hole, through the horizon, and inside it when the camera falls in.

The horizons lie where $\Delta = r^2 - 2r + a^2$ goes to zero: $r_\pm = 1 \pm \sqrt{1 - a^2}$. The outer one, $r_+$, is the event horizon. Between it and the surface $r = 1 + \sqrt{1 - a^2\cos^2\theta}$ lies the ergosphere: there you can’t stay in one place, because space itself drags everything around the axis. That’s why a stationary camera isn’t allowed any closer to the hole than the ergosphere.

How light travels

A ray of light is a null geodesic. The most convenient way to write it down is as motion in Hamiltonian mechanics, with momentum $p_\mu$ and the Hamiltonian

$$H = \tfrac12\, g^{\mu\nu} p_\mu p_\nu = 0,\qquad \frac{dx^\mu}{d\lambda} = \frac{\partial H}{\partial p_\mu},\qquad \frac{dp_\mu}{d\lambda} = -\frac{\partial H}{\partial x^\mu}.$$

In Kerr–Schild form the inverse metric is just as simple: $g^{\mu\nu} = \eta^{\mu\nu} - f\, l^\mu l^\nu$. We take the derivatives with respect to the coordinates analytically, with no numerical differentiation. The integration runs on the graphics card: fourth-order Runge–Kutta with a step size proportional to the distance from the center. On average a ray needs 40–60 steps.

Rays are traced backward in time — from the camera to the source. A pixel’s direction is a direction in the camera’s own frame of reference (an orthonormal tetrad built around its 4-velocity). As a result, aberration and frequency shift for a moving camera come out automatically. A ray ends up in one of three places:

  • on the disk — it crossed the equatorial plane between the inner and outer edges; the crossing point is refined by the secant method within the step;
  • on the sky — it got farther out than 15,000 M; its direction of travel is exactly the point on the sky the light came from;
  • in the hole — it ended up inside the innermost circular photon orbit and is moving inward. Light doesn’t turn back from there, so the pixel is black.

Along the ray, the energy $E = -p_t$, the angular momentum $L = p_\phi$ and the Carter constant $Q$ are conserved. The reference integrator (Dormand–Prince with error control, in double precision) holds $H$ at zero to within $10^{-10}$, and $Q$ to a relative error of $10^{-8}$.

The shadow and the photon ring

A Kerr black hole has spherical photon orbits: light on them winds around the hole forever. An orbit of radius $r$ carries the constants (Bardeen, 1973; Teo, 2003)

$$\xi = \frac{L}{E} = -\frac{r^3 - 3r^2 + a^2 r + a^2}{a\,(r - 1)},\qquad \eta = \frac{Q}{E^2} = -\frac{r^3\,(r^3 - 6r^2 + 9r - 4a^2)}{a^2\,(r - 1)^2},$$

and $r$ ranges over the interval between the prograde and retrograde circular photon orbits, $r_{\text{ph}} = 2\big(1 + \cos\big(\tfrac23 \arccos(\mp a)\big)\big)$. Rays with these constants that reach the camera trace out the edge of the shadow. The “Analytic shadow” dashed line is built exactly this way: the constants are converted into a momentum in Boyer–Lindquist coordinates, then into Kerr–Schild coordinates, then into the camera’s frame. The renderer plays no part in it.

For a non-spinning hole, the shadow is a circle of radius $\sqrt{27}\,M \approx 5.2\,M$ (the horizon is $2M$). With spin, the shadow shifts sideways and flattens on the side rotating toward us: rays traveling with the rotation can get closer in. The thin bright line along the edge of the shadow is the photon ring — images of the disk stacked on top of one another after wrapping around the hole once, twice or more.

The disk

The disk follows the Novikov–Thorne model (1973): thin and opaque, with gas on circular Keplerian orbits at angular velocity $\Omega = 1/(r^{3/2} + a)$. Its inner edge is the innermost stable circular orbit (ISCO; Bardeen, Press and Teukolsky, 1972):

$$r_{\text{ISCO}} = 3 + Z_2 - \sqrt{(3 - Z_1)(3 + Z_1 + 2 Z_2)},$$ $$Z_1 = 1 + \sqrt[3]{1 - a^2}\,\big(\sqrt[3]{1 + a} + \sqrt[3]{1 - a}\big),\quad Z_2 = \sqrt{3a^2 + Z_1^2}.$$

That’s 6 M with no spin and 1.24 M at $a = 0.998$. The energy flux per unit area is given by the Page–Thorne formula (1974):

$$F(r) = -\frac{\dot M}{4\pi \sqrt{-g}}\,\frac{\Omega_{,r}}{(E - \Omega L)^2} \int_{r_{\text{ISCO}}}^{r} (E - \Omega L)\, L_{,r}\, dr.$$

Every patch of the disk glows as a blackbody: $\sigma T^4 = F$. So at the very inner edge the temperature drops to zero (there’s no friction there), it peaks a little farther out, and far away it falls off as $T \propto r^{-3/4}$. The absolute temperature depends on the accretion rate and the mass. Real disks are hotter — they shine in ultraviolet and X-rays — and in visible light they would look uniformly bluish white. We take a peak of about 6500 K so that the colors land in the visible part of the spectrum; you can change it in the settings.

Color and brightness

Light arriving from the disk is shifted in frequency by a factor of $g$:

$$g = \frac{\nu_{\text{obs}}}{\nu_{\text{em}}} = \frac{-p_\mu u^\mu_{\text{camera}}}{-p_\mu u^\mu_{\text{gas}}} = \frac{1}{u^t\,(-p_t - \Omega\, p_\phi)},$$

if the momentum is normalized to unit energy in the camera’s frame. Everything is already folded into $g$: the gas moving toward us or away from us (Doppler), its clocks running slow on the orbit, and the light climbing out of the gravitational well.

From here, Liouville’s theorem does the work: $I_\nu/\nu^3$ stays constant along a ray. So a blackbody at temperature $T$, seen with a shift $g$, is exactly a blackbody at temperature $gT$. The color comes from a “temperature → color” table, which we compute by integrating Planck’s law against the CIE 1931 color matching functions (in the approximation by Wyman, Sloan and Shirley, 2013). The brightness growing as $g^4$ falls out of this on its own.

The switches in the settings pull the effects apart. With “Doppler effect” off, the gas is replaced by a zero-angular-momentum observer (ZAMO): it is at rest relative to the local space, and only frame dragging carries it around the hole. With “Gravitational shift” off, we divide out the shift for such an observer. Turn off both and you get $g = 1$, as in “Interstellar”.

Time $t$ is integrated along the ray too, so the pattern on the disk is taken at the moment the light was emitted, not the moment it was received. Light from the far side of the disk takes longer to arrive, so its image “lags” slightly behind.

The sky

The background is Earth’s real sky: the Milky Way from NASA’s “Deep Star Maps 2020” panorama and the 9,096 stars of the Yale Bright Star Catalogue. The stars are drawn not as a texture but as points. Each pixel knows, from its neighboring rays, which patch of sky lands in it — and if a star falls in that patch, the star’s flux is multiplied by the magnification of the gravitational lens

$$\mu = \frac{\Omega_{\text{pixel}}}{\Omega_{\text{its footprint on the sky}}}.$$

That’s why a star near the edge of the shadow flares up and splits in two, and its images move toward each other. Nobody draws the multiple images on purpose: different pixels simply receive light from the same star. A star’s color is that of a blackbody at its effective temperature, also shifted by $g$.

The fall

The “Fall in” button releases the camera with a small sideways push, chosen by computing the whole world line in advance: the strongest push with which the camera still falls into the hole and never passes through the disk. Its 4-velocity obeys the same Hamilton’s equations, only with $H = -\tfrac12$. The integrator is Dormand–Prince in double precision, stepping in the camera’s proper time $\tau$. Playback is sped up far from the hole and slowed down close to it, so that every “octave” of distance takes about the same amount of time. The panel shows the true $r$ and $\tau$.

A falling observer doesn’t notice the horizon: locally there’s nothing special about it. But past it, the radius turns into time — it shrinks as inevitably as a clock ticks. In a non-spinning hole, the path ends at the singularity. From the horizon to it takes no more than $\pi M$ of proper time: about a minute for Sagittarius A*, 0.15 ms for a hole of 10 solar masses. In a spinning hole, you first run into the inner horizon $r_-$. There, radiation falling in behind you is blueshifted without limit, and the classical picture itself (Kerr’s solution) stops being reliable. That’s where we stop.

Tidal forces near the horizon scale as $M/r^3$. A stellar-mass hole would stretch a person apart long before the horizon. A supermassive one can be crossed without feeling a thing.

What’s simplified

  • The disk is infinitely thin. There’s no corona, no jet and no hot, thick flow — which is what actually dominates in M87* and Sgr A*: the EHT images were taken at radio wavelengths, and it’s exactly that kind of flow that glows there.
  • The clumps in the disk are an illustration: fractal noise that rotates with the gas and gets stretched by differential rotation. The temperature profile, velocities and shifts don’t depend on it; the “Turbulence” switch shows a smooth disk.
  • The disk temperature is chosen so that the colors are visible, and the sky is brightened so that stars can be seen next to the disk.
  • There’s no emission from inside the ISCO, and no polarization.
  • Under strong shifts, the color of the Milky Way panorama is recalculated approximately, as if it were a 6500 K blackbody; the stars are handled exactly.
  • The graphics card computes with 32-bit numbers. Rays that keep winding around near the horizon for more than 700 steps are counted as captured — the reference integrator gives the same answer for them.

How it was checked

The physics is written once, in JavaScript, in double precision; the shader mirrors it line for line. Automated tests compare it against known answers:

TestExpectedResult
Kerr–Schild metric transformed into Boyer–Lindquist coordinatesBoyer–Lindquist form±10⁻¹⁰
Analytic gradient of the Hamiltonian vs. finite differencesagreement±10⁻⁷
ISCO at a = 0 and a = 0.9986 M; 1.2369 Mexact
Schwarzschild critical impact parameter√27 M = 5.19615 M±2·10⁻⁶
Deflection of light far from the hole4M/b±1%
Shadow edge at a = 0.5, 0.9, 0.998 vs. Bardeen’s curveagreement±2·10⁻⁴ of the width
Frequency shift of a Schwarzschild disk seen from above√(1 − 3M/r)±10⁻⁷
Page–Thorne flux: the integral vs. the closed-form formulaagreement±2·10⁻⁴
Blackbody color at 3000, 6500, 10,000 K vs. the CIE Planckian locus(x, y)±0.003
Free-fall time to the horizon vs. the cycloidagreement±0.01 M
Fast RK4 (as in the shader) vs. the referencethe same frame98.5% of pixels; r ±0.2%; sky ±0.05°
The GPU’s ray map (32-bit, in a real browser) vs. the reference, 4 views and spinsthe same frame≥ 99.9% of pixels; r ±0.4%; sky ±0.1°

Sources

  1. R. P. Kerr. Gravitational field of a spinning mass as an example of algebraically special metrics. Phys. Rev. Lett. 11, 237 (1963).
  2. J. M. Bardeen, W. H. Press, S. A. Teukolsky. Rotating black holes: locally nonrotating frames, energy extraction, and scalar synchrotron radiation. ApJ 178, 347 (1972).
  3. J. M. Bardeen. Timelike and null geodesics in the Kerr metric. In Black Holes (Les Astres Occlus), 1973.
  4. I. D. Novikov, K. S. Thorne. Astrophysics of black holes. In Black Holes (Les Astres Occlus), 1973.
  5. D. N. Page, K. S. Thorne. Disk-accretion onto a black hole. ApJ 191, 499 (1974).
  6. J.-P. Luminet. Image of a spherical black hole with thin accretion disk. A&A 75, 228 (1979).
  7. E. Teo. Spherical photon orbits around a Kerr black hole. Gen. Rel. Grav. 35, 1909 (2003).
  8. O. James, E. von Tunzelmann, P. Franklin, K. S. Thorne. Gravitational lensing by spinning black holes in astrophysics, and in the movie Interstellar. Class. Quantum Grav. 32, 065001 (2015).
  9. C. Wyman, P.-P. Sloan, P. Shirley. Simple analytic approximations to the CIE XYZ color matching functions. JCGT 2(2), 2013.
  10. Event Horizon Telescope Collaboration. First M87 Event Horizon Telescope results (2019); First Sagittarius A* results (2022).
  11. NASA Scientific Visualization Studio. Deep Star Maps 2020. Yale Bright Star Catalogue, 5th ed.