DEV Community

Cover image for Tracing light around a spinning black hole in the browser (WebGL2, no libraries)
sishiya (しし屋)
sishiya (しし屋)

Posted on

Tracing light around a spinning black hole in the browser (WebGL2, no libraries)

#physics #simulation #javascript #webgl

The black hole in the movie Interstellar has a flat disk that seems to wrap over the top and under the bottom of the hole. That is not an artistic choice. It is what gravity does to light.

To see it for myself, I made a page that traces light, ray by ray, around a spinning (Kerr) black hole and draws the result. It is plain WebGL2 with no libraries.

https://sishiya.github.io/manabi-lab/black-hole/?lang=en

Click any point on the image and the page shows where that light came from, as a side view and a top view of its path. This post covers how it works and where I got stuck.

What it does

Like ordinary ray tracing, it follows light backwards from the camera. The difference is that light does not travel in straight lines. For every pixel, the page numerically integrates the geodesic equations of Kerr spacetime.

Each ray ends in one of three ways:

  • It falls into the hole → black (the shadow)
  • It crosses the disk plane inside the disk → disk color
  • It escapes far away → look up the star map in its final direction

The equations

Units are G = c = M = 1 (length and time measured in black-hole masses), in Boyer–Lindquist coordinates (r, θ, φ). With the photon energy fixed at p_t = −1, there are two conserved quantities: the angular momentum b = L/E and Carter's constant q = Q/E².

Σ = r² + a²cos²θ,  Δ = r² − 2r + a²
P = r² + a² − ab,  K = (b − a)² + q,  R = P² − ΔK

dr/dλ   = Δ p_r / Σ
dθ/dλ   = p_θ / Σ
dφ/dλ   = (aP/Δ − a + b/sin²θ) / Σ
dp_r/dλ = (R′/Δ − RΔ′/Δ² − Δ′p_r²) / (2Σ)     (R′ = 4rP − Δ′K, Δ′ = 2r − 2)
dp_θ/dλ = (b²cosθ/sin³θ − a²sinθcosθ) / Σ
Enter fullscreen mode Exit fullscreen mode

For light the Hamiltonian is H = 0, which removes the derivative of Σ. I derived these myself and checked them against the appendix of the paper by the team that made the movie's images (James et al. 2015). The sign of φ is flipped compared with the paper, because the paper writes the equations for tracing backwards in time.

In JavaScript (y = [r, θ, φ, p_r, p_θ]):

function deriv(y, b, q, a, out) {
  const r = y[0], th = y[1], pr = y[3], pth = y[4];
  let s = Math.sin(th); const c = Math.cos(th);
  if (Math.abs(s) < 1e-6) s = s < 0 ? -1e-6 : 1e-6;
  const S = r * r + a * a * c * c, D = r * r - 2 * r + a * a;
  const P = r * r + a * a - a * b, K = (b - a) * (b - a) + q, R = P * P - D * K;
  const Dp = 2 * r - 2, Rp = 4 * r * P - Dp * K;
  out[0] = D * pr / S;
  out[1] = pth / S;
  out[2] = (a * P / D - a + b / (s * s)) / S;
  out[3] = (Rp / D - R * Dp / (D * D) - Dp * pr * pr) / (2 * S);
  out[4] = (b * b * c / (s * s * s) - a * a * s * c) / S;
  return out;
}
Enter fullscreen mode Exit fullscreen mode

This is stepped with 4th-order Runge–Kutta, with a negative step because we go backwards.

The camera: an observer swept along with space

Near a spinning black hole, space itself is dragged around in the direction of the spin. So the camera is a ZAMO (zero-angular-momentum observer, Bardeen, Press and Teukolsky 1972), which rotates along with that dragging.

From the pixel direction n, the page computes the photon energy E_F that this observer measures, and the conserved quantities b and q:

const EF = 1 / (z.alpha + z.omega * z.varpi * nph);   // photon energy seen by the camera
const pr = EF * z.rho / Math.sqrt(z.D) * nr;
const pth = EF * z.rho * nth;
const b = EF * z.varpi * nph;                         // angular momentum L/E
const q = pth * pth + c * c * (b * b / (s * s) - a * a); // Carter's constant
Enter fullscreen mode Exit fullscreen mode

alpha is how fast clocks tick there, omega is the angular velocity of frame dragging, and varpi is the circumferential radius.

Step size

Paths bend more sharply close to the hole, so the step is proportional to the distance from the horizon:

const STEP_K = 0.06, STEP_MIN = 0.004;
function stepLen(r, rh) { return Math.max(STEP_K * (r - rh), STEP_MIN); }
Enter fullscreen mode Exit fullscreen mode

A ray stops when:

  • r < r₊ + 0.01 → the hole (r₊ is the event horizon)
  • it crosses the equatorial plane between the disk's inner and outer edges → the disk
  • r > max(2 r_cam, 60) and moving outward → the sky
  • 900 steps pass without any of these → the hole (light circling the photon orbit)

Where I got stuck: a vertical line down the middle

The first render had a flickering vertical line down the center of the image.

Rays in the central column have b ≈ 0, so they pass right over or under the spin axis. There sinθ is small, and the b/sin²θ term in dφ/dλ becomes very large. With a step based only on r, φ jumped a lot in a single step and the integration broke down.

The fix was to shrink the step so that one step never turns more than 0.04 rad:

function rk4(y, h, b, q, a) {
  deriv(y, b, q, a, _k1);
  const lim = 0.04 / (Math.abs(_k1[2]) + Math.abs(_k1[1]) + 1e-6);
  if (Math.abs(h) > lim) h = Math.sign(h) * lim;
  // ...then an ordinary RK4 step
}
Enter fullscreen mode Exit fullscreen mode

The line disappeared.

A related problem: rays with b exactly 0 cross the axis, and θ goes negative (or beyond π). Computing the disk hit point with θ = π/2 then swapped the near and far sides. The shader adds π to φ when sinθ < 0, and the CPU side interpolates between path points to find the position.

Disk color: the frequency shift

The disk gas follows Keplerian orbits. The inner edge is the ISCO (innermost stable circular orbit, 3.83M for a = 0.6), and the outer edge is 18.7M, the same as the figure in the paper.

The frequency shift g of the light is:

g = E_F / (u^t (1 − Ω b))
Ω = 1/(r^{3/2} + a)
u^t = (r^{3/2} + a) / (r^{3/4} √(r^{3/2} − 3r^{1/2} + 2a))
Enter fullscreen mode Exit fullscreen mode

At the movie's viewing angle (a = 0.6), g is about 1.36 on the left side, where the gas comes toward us, and about 0.76 on the band above the hole.

For clarity, the movie left out this shift in color and brightness (the paper explains why). So the page lets you switch between "movie" and "realistic". In realistic mode the color temperature is multiplied by g and the brightness by g⁴: the left side gets bright and the right side gets dim.

Speed: two levels

With up to 900 integration steps per pixel, the full screen cannot be traced every frame.

  • Coarse level (about 160,000 pixels): retraced in one go whenever the camera moves
  • Fine level (2× screen resolution, up to 5 million pixels): starts 140 ms after the camera stops, in 192 px tiles from the center outward, as many tiles per frame as the frame time allows

The tracing pass only writes its result (kind, r and φ where it hit, g) to an RGBA32F texture. The disk pattern (which rotates), color and bloom are painted every frame by separate shaders.

Also, the camera's φ (dragging left and right) is a rotation around the spin axis, so no retracing is needed: the page just rotates the sky and the disk pattern. Rays are retraced only when θ, distance or spin change.

On my PC (NVIDIA), the whole fine level takes about 20 ms.

Writing the same equations twice, on GPU and CPU

The ray diagram for a clicked point is retraced on the CPU (double precision, up to 4000 steps). Since the same equations live in both the shader and JavaScript, I keep a check that compares them at the same point: the value read back from the GPU texture against the CPU result.

__bh.renderNow();
__bh.readGeo(0.04, 0.5);       // GPU: [1, ~6.8, ~3.0, ~0.76]  (disk, r, φ, g)
__bh.pickAt(0.04, 0.5).rHit;   // CPU: ~6.8
Enter fullscreen mode Exit fullscreen mode

What the bands in the movie really are

At the movie's angle (a = 0.6, r = 28, θ = 85°), tracing up the central vertical line gives:

Height on screen What you see How much the light bent
Upper part Top face of the disk behind the hole (the band above the hole) ~41°
Just above the shadow edge Bottom face of the disk in front (light loops under and comes back) ~210°
Just below the shadow edge Top face of the disk in front ~270°
Lower part Bottom face of the disk behind the hole (the band below the hole) ~100°

The band above the hole is the top of the disk behind the hole. Its light bends over the top of the hole to reach us. The band below is the underside of the same disk, and its light reaches us by passing under the hole.

Still missing

  • The thin rings at the shadow edge (second and higher images) are only 1–3 pixels wide and can look dotted
  • Aberration from camera motion, disk thickness, and the disk lighting itself are not included

Main sources

  • O. James, E. von Tunzelmann, P. Franklin, K. S. Thorne (2015) "Gravitational lensing by spinning black holes in astrophysics, and in the movie Interstellar", Classical and Quantum Gravity 32, 065001 (arXiv:1502.03808)
  • J. M. Bardeen, W. H. Press, S. A. Teukolsky (1972) "Rotating black holes: locally nonrotating frames, energy extraction, and scalar synchrotron radiation", The Astrophysical Journal 178
  • N. I. Shakura, R. A. Sunyaev (1973) "Black holes in binary systems. Observational appearance", Astronomy and Astrophysics 24
  • K. S. Thorne (2014) The Science of Interstellar, W. W. Norton

This page is part of Manabi Lab, a series of small learning apps that run in the browser (English and Japanese). The code and this post were written together with AI (Claude).

https://sishiya.github.io/manabi-lab/?lang=en

Top comments (0)