When Two Shadows Collide: Ray-Tracing a Binary Black Hole Merger

In Part 1 (Schwarzschild), Part 2 (Morris–Thorne Wormhole), and Part 3 (Boyer–Lindquist Kerr), every spacetime we ray-traced shared a hidden luxury: integrability. A single black hole—even a spinning one—has enough continuous symmetries (energy $E$, axial angular momentum $L_z$, and Carter's constant $Q$) to reduce every photon trajectory to decoupled 1D equations.

What happens when two black holes orbit each other on the verge of merging? All continuous spatial symmetries vanish. Light rays can skim the photon sphere of Black Hole 1, cross the gravitational saddle point between the pair, wind around Black Hole 2 in a figure-eight ($\infty$), and either escape to the camera or plunge into the companion horizon. As the Simulating eXtreme Spacetimes (SXS) collaboration and the Interstellar team showed (Bohn et al., 2015), a binary black hole does not just cast two circular shadows: each primary shadow is flanked by a fractal hierarchy of thin, crescent-shaped "eyebrow" shadows bridged by a glowing saddle caustic.

This post builds a real-time WebGL ray tracer for an orbiting binary black hole system using the exact Majumdar–Papapetrou two-hole metric, Fermat's optical refractive index formulation, and a three-zone accretion model (two tidally truncated mini-disks inside a circumbinary cavity).

0. Anatomy of a Binary Black Hole Image

Measuring lengths and masses in units of the total system mass $M_{\text{tot}} = M_1 + M_2 = 1$ ($G = c = 1$), the binary geometry is governed by the separation $d = |\vec x_1 - \vec x_2|$ and mass ratio $q = M_2 / M_1 \in (0, 1]$:

Metric

Two-center Majumdar–Papapetrou spacetime in isotropic coordinates $(t, \vec x)$, where the harmonic potential $U(\vec x) = 1 + M_1/r_1 + M_2/r_2$ satisfies $\nabla^2 U = 0$ exactly.

Punctured Horizons ($\vec x_i$)

In isotropic coordinates, each point $\vec x_i$ ($r_i = |\vec x - \vec x_i| \to 0$) is not a point singularity but a regular event horizon of finite surface area $A_i = 4\pi M_i^2$ (Hartle & Hawking, 1972).

Primary Shadows

Rays captured directly by the local photon shell ($r_{\text{ph},i} \approx M_i$) of either black hole cast two facing teardrop-flattened shadows.

Eyebrow Shadows

Rays aimed just outside Primary Shadow 1 are deflected across the inter-hole saddle point and captured by Black Hole 2, carving thin crescent "eyebrows" outside the main shadows.

Three-Zone Disks

Gravitational torques clear a central cavity inside the outer circumbinary disk ($r > r_{\text{cb,in}}$) while feeding two co-moving mini-disks around $\vec x_1(t)$ and $\vec x_2(t)$.

1. The Majumdar–Papapetrou Metric and Fermat's Refractive Index

In general relativity, Einstein's field equations are nonlinear: you cannot normally superpose two black hole solutions by adding their potentials. Remarkably, Majumdar (1947) and Papapetrou (1945) discovered that for extremal black holes in Einstein–Maxwell theory, the nonlinear gravitational cross-terms cancel the Maxwell stress tensor identically, reducing the entire Einstein system to the linear 3D Laplace equation $\nabla^2 U = 0$:

For a real-time GPU ray tracer, the isotropic form of the Majumdar–Papapetrou metric offers an extraordinary simplification. Along any null geodesic ($ds^2 = 0$), the coordinate travel time $dt$ and Euclidean arc length $ds = \sqrt{dx^2 + dy^2 + dz^2}$ obey:

Applying the Euler–Lagrange equation $\frac{d}{ds}(n\hat v) = \nabla n$ to the unit tangent vector $\hat v = d\vec x/ds$ ($|\hat v| = 1$) gives the exact spatial curvature equation integrated per pixel:

// Exact Majumdar-Papapetrou potential U(x) and log-index gradient g(x) = 2 * grad(U) / U
void mpField(vec3 p, out float U, out vec3 gVec, out float r1, out float r2) {
  vec3 d1 = p - uPos1;
  vec3 d2 = p - uPos2;
  r1 = length(d1);
  r2 = length(d2);
  float r1Safe = max(r1, 0.02);
  float r2Safe = max(r2, 0.02);
  U = 1.0 + uM1 / r1Safe + uM2 / r2Safe;
  vec3 gradU = -uM1 * d1 / (r1Safe * r1Safe * r1Safe) - uM2 * d2 / (r2Safe * r2Safe * r2Safe);
  gVec = (2.0 / U) * gradU;
}

// Inside per-pixel ray loop (pos = x(s), vel = dx/ds with |vel| = 1):
float U; vec3 gVec;
mpField(pos, U, gVec, r1, r2);
if (r1 < cap1 || (r2 < cap2 && uM2 > 0.001)) { captured = true; break; }

float dMin = min(r1 - cap1, r2 - cap2);
float ds   = clamp(0.11 * max(dMin, 0.04) * (1.0 + 0.06 * rBary), 0.016, 0.75);

// Second-order midpoint step on (pos, vel) preserving |vel| = 1
vec3 acc1   = gVec - dot(vel, gVec) * vel;
vec3 posMid = pos + 0.5 * vel * ds;
vec3 velMid = normalize(vel + 0.5 * acc1 * ds);

float UMid; vec3 gMid; float r1m, r2m;
mpField(posMid, UMid, gMid, r1m, r2m);
vec3 accMid = gMid - dot(velMid, gMid) * velMid;

pos += velMid * ds;
vel  = normalize(vel + accMid * ds);

2. Figure-Eight Photon Orbits and "Eyebrow" Shadows

Click the Eyebrow Bridge (d=2.4M) preset in the simulation above and toggle Accretion Disks off for a moment so you can see the pure silhouette against the starfield. Alongside the two large black shadows, you will see thin dark arcs curving around their outer rims, linked across the center by a bright gold-cyan caustic bridge.

Where do those extra shadows come from? For a single isolated Majumdar–Papapetrou black hole ($M_2 = 0$), setting the inward optical curvature $|\vec g(r)| = \frac{2M_1}{r^2(1 + M_1/r)}$ equal to the circular orbit curvature $1/r$ gives a single photon sphere at isotropic radius $r_{\text{ph}} = M_1$, casting a single circular shadow of apparent radius $b_c = U^2(M_1)\,M_1 = 4M_1$.

When two black holes sit at $\vec x_1$ and $\vec x_2$, the potential $U(\vec x)$ develops an unstable gravitational saddle point ($\nabla U(\vec x_{\text{saddle}}) = \vec 0$) on the axis connecting them:

Because photons can cross the saddle and trade holes, the binary shadow decomposes into a self-similar hierarchy indexed by how many times a ray switches partners before plunging through a horizon:

• 0-switch rays (Primary Shadows): Rays aimed directly inside the local photon shell of Hole 1 or Hole 2 plunge straight into that hole, forming the two main shadows (flattened into facing teardrops by the companion's tidal pull).
• 1-switch rays (Primary "Eyebrows"): A narrow fan of rays aimed just outside Hole 1's shadow is bent by $\sim 180^\circ$ around Hole 1, flung across $\vec x_{\text{saddle}}$, and swallowed by Hole 2—carving a dark crescent eyebrow shadow belonging to Hole 2 right outside Hole 1's rim (and vice versa).
• $N$-switch rays (Figure-Eight Caustic Bridge): At the boundary between escaping and captured 1-switch rays sit unstable periodic figure-eight ($\infty$) photon orbits that weave back and forth around both horizons forever. The gold-cyan Photon Bridge pill highlights rays that skim within $2.4\,M_i$ of both black holes (`minR1 < 2.4 && minR2 < 2.4`), illuminating this chaotic separatrix in real time.

3. Circumbinary Cavities, Mini-Disks, and Binary Doppler Beaming

General-relativistic magnetohydrodynamic (GRMHD) simulations of supermassive black hole binaries show that gas does not form a single smooth disk all the way to the horizons. Instead, gravitational torques from the orbiting pair clear out a central cavity ($r_{\text{cb,in}} = 0.80\,d + 1.50\,M_{\text{tot}}$ in the shader, smoothly approaching the single-hole inner edge as $d \to 0$) surrounded by an outer circumbinary disk with $m=2$ spiral density waves, while streams of gas feed two tidally truncated mini-disks ($r_{\text{out},i} = 0.42\,d$) around $\vec x_1(t)$ and $\vec x_2(t)$.

When a ray crosses the equatorial plane ($y = 0$), the shader evaluates this three-zone geometry and computes the total velocity $\vec\beta(\vec x)$ of the emitting plasma by superposing the binary orbital velocity $\vec v_i = \vec\Omega_{\text{orb}} \times \vec x_i$ with the local circular orbital speed around Hole $i$. In a single-center Majumdar–Papapetrou metric ($U = 1 + M_i/r_i$), the physical circular-orbit speed measured by a local static observer is $\beta_{\text{circ}}(r_i) = \sqrt{-r_i U' / (U + r_i U')} = \sqrt{M_i / r_i}$ exactly (reaching $c$ at the photon sphere $r_{\text{ph},i} = M_i$). Combining the static lapse redshift $\sqrt{-g_{tt}} = U^{-1}(\vec x)$ with the special-relativistic boost gives:

// Mini-disk 1 around BH 1 (tidally truncated at rOut1 = 0.42 * d)
if (uSep > 0.9 && d1 > rIn1 && d1 < rOut1) {
  float x = d1 / rIn1;
  env = smoothstep(rIn1, rIn1 * 1.12, d1) * (1.0 - smoothstep(rOut1 * 0.68, rOut1, d1));
  emitFlux = 2.8 * pow(1.0 / x, 2.1) * env;
  baseTemp = 9200.0 * pow(x, -0.55);
  vec3 tang1 = normalize(vec3(-rel1.z, 0.0, rel1.x));
  vGas = uVel1 + clamp(sqrt(uM1 / max(d1, 0.25)), 0.0, 0.65) * tang1;
}
// ... (Mini-disk 2 around BH 2 and outer Circumbinary Disk with m=2 spiral arm) ...

// Exact MP lapse redshift (1 / U) * special-relativistic Doppler boost
float Uhit     = 1.0 + uM1 / max(d1, 0.15) + uM2 / max(d2, 0.15);
float grav     = 1.0 / Uhit;
float beta2    = clamp(dot(vGas, vGas), 0.0, 0.78);
float betaCosT = dot(vGas, -vel); // beta * cos(theta), since |vel| = 1 and |vGas| = beta
float dop      = sqrt(1.0 - beta2) / max(1.0 - betaCosT, 0.16);
float g        = clamp(grav * mix(1.0, dop, uDoppler), 0.18, 3.0);

vec3 bb = blackbody(baseTemp * g);
accumCol += T * bb * emitFlux * pow(g, 3.2) * (0.35 + 1.25 * turb) * 3.8;

4. The Binary Ray-Tracing Pipeline

Each frame updates the Keplerian barycentric positions $(\vec x_1(t), \vec x_2(t))$ on the CPU, hands them as uniforms to the fullscreen fragment shader, and integrates the Fermat optical curvature across the two-center potential:

5. Across the Series: Single vs. Spinning vs. Binary Black Holes

Looking back across all three black hole ray tracers shows how symmetry—and its breaking—shapes the optical appearance of curved spacetime:

Drag the Binary Separation · $d/M$ slider in the live demo slowly from $4.0M$ down to $0.0M$: watching the two mini-disks merge into the circumbinary ring while the twin shadows and their outer eyebrows coalesce into a single horizon is a front-row view of the final seconds of a binary black hole inspiral.