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, in WebGPU compute shaders or, where there is no WebGPU, through WebGL 2: fourth-order Runge–Kutta with a step size proportional to the distance from the center: 5% of $r$ near the hole, up to 30% far away. On average a ray needs 25–45 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. Often this is clear earlier: below the disk’s inner edge a ray is checked against the radial potential $R(r) = E^2 r^4 + (2a^2E^2 - 2aLE - K)\,r^2 + 2MKr - a^2Q$, with $K = Q + (L - aE)^2$. If $R$ does not vanish between the ray and the horizon, there is no turning point left, and the tracing can stop at once. Carter’s constant $Q$ comes straight from the Cartesian components of the momentum. The minima of $R$ come from the cubic $R'(r) = 0$. A test checks that stopping this way never changes where a ray ends.

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, the mass and the spin: $T_\text{peak} \propto (\dot M\, F_\text{max}(a))^{1/4}$. A fast spin brings the disk’s edge in toward the horizon and heats it up. At the same rate, a non-spinning hole’s peak is 0.45 times that at spin 0.9, and at spin 0.998 it is 1.94 times higher. Real disks are hotter — they shine in ultraviolet and X-rays — and in visible light they would look uniformly bluish white. We take a rate at which a hole of spin 0.9 peaks at about 7500 K, so that the colors land in the visible part of the spectrum; you can change the rate in the settings. The disk’s brightness is then evened out halfway (in log), so it doesn’t go black or white: the difference stays visible, just softened. The sky is left alone.

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”.

No switch swaps one picture for another. Each is a dial from 0 to 1 that slides to its new setting. The frequency shifts enter as $g_\text{grav}^{\,k_1}\, g_\text{dopp}^{\,k_2}$. For the bending of light the dial is the mass of the metric itself: it goes down to zero together with the spin ($M \to kM$, $a \to ka$). So every state along the way is an exact Kerr solution, a lighter and lighter hole, and at $k = 0$ space is flat. Meanwhile the hole itself is drawn as a black ball the size of its horizon. The spin, the temperature and everything else change just as continuously when you go from one hole to another.

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 hot spot shows this best. It is a clump of gas about $1\,M$ across on a Kepler orbit at $7\,M$, noticeably hotter than the disk around it. For every pixel its position is taken at the moment the light left the disk. So the spot has several images arriving at different times: the direct one and the one that went round the hole. Passing behind the hole, the spot is stretched by the lens into an arc over the shadow; on the side coming towards us it flares up with Doppler beaming. In 2018 the GRAVITY interferometer saw flares like this around Sagittarius A*: the spot went round in about 45 minutes, and an orbit at $7\,M$ for a mass of 4.3 million suns takes 43.

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 shaders of both engines mirror 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, WebGPU and WebGL 2) 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. GRAVITY Collaboration. Detection of orbital motions near the last stable circular orbit of the massive black hole SgrA*. A&A 618, L10 (2018).
  11. Event Horizon Telescope Collaboration. First M87 Event Horizon Telescope results (2019); First Sagittarius A* results (2022).
  12. NASA Scientific Visualization Studio. Deep Star Maps 2020. Yale Bright Star Catalogue, 5th ed.