arman1o1 blog
Astrophysics Physics WebGL GLSL

Spinning Up Gargantua: A Real-Time Kerr Black Hole in Your Browser

• ☕ 11 min read
Live Boyer–Lindquist Kerr ray tracer. Every pixel above integrates a null geodesic in Mino time through a spinning Kerr spacetime ($a/M = 0.990$ by default). Drag to orbit, open ⚙ Controls to sweep spin $a/M \in [0, 0.998]$, flip between Prograde (+a) and Retrograde (-a) disks, toggle the Ergosphere static-limit shell, turn on the D-Shadow Guide to compare Bardeen's shifted D-shadow against the dashed Schwarzschild circle, or open the simulation full-screen.

In Part 1 (Rendering Gargantua), we built a real-time ray tracer for a Schwarzschild black hole and left an honest caveat in the table at the end: the real Gargantua in Interstellar is not Schwarzschild. To keep Miller's planet surviving on a stable orbit with a 60,000-to-1 time-dilation factor, Kip Thorne had to spin Gargantua up to within $10^{-14}$ of the speed-of-light limit ($a/M \approx 1 - 10^{-14}$).

When a black hole spins, spacetime itself is dragged into a relentless whirlpool. Spherical symmetry breaks down to axial symmetry; the event horizon shrinks from $2M$ down toward $1M$; an oblate ergosphere bulges outside the horizon where nothing—not even light—can stand still against the spin; the innermost stable circular orbit (ISCO) plunges from $6M$ all the way down to the horizon for co-rotating gas while being flung out to $9M$ for counter-rotating gas; and the circular shadow flattens on its approaching edge into James Bardeen's famous D-shaped silhouette.

This post upgrades the fullscreen WebGL shader to the exact Kerr metric in Boyer–Lindquist coordinates, using Brandon Carter's fourth constant of motion and Mino time to integrate turning-point-free polynomial geodesics at 60 fps right in the canvas above.

0. What Changes When Spacetime Spins

Working in geometric units $G = c = M = 1$, a rotating black hole is completely characterized by a single dimensionless spin parameter $a \equiv J/M^2 \in [0, 1)$, where $J$ is angular momentum:

Outer Horizon ($r_+$)
Located where $\Delta(r) \equiv r^2 - 2r + a^2 = 0$, giving $r_+ = 1 + \sqrt{1 - a^2}$. As $a$ increases from $0$ to $1$, the horizon radius shrinks from $2M$ to $1M$.
Ergosphere ($r_{\text{ergo}}(\theta)$)
The static limit where $g_{tt} = 0$, given by $r_{\text{ergo}}(\theta) = 1 + \sqrt{1 - a^2\cos^2\theta}$. It touches the horizon at the poles ($\theta = 0, \pi$) and bulges out to $2M$ at the equator ($\theta = \pi/2$).
Frame Dragging ($\omega$)
Even a zero-angular-momentum observer ($L_z = 0$) is swept around the spin axis at coordinate angular velocity $\omega(r, \theta) = -g_{t\varphi}/g_{\varphi\varphi} = 2ar / A$.
Split ISCO ($r_{\text{ISCO}}$)
Co-rotating (prograde) orbits ride the frame dragging down to $r_{\text{ISCO}} \to 1M$ as $a \to 1$, while counter-rotating (retrograde) orbits are destabilized inside $r_{\text{ISCO}} \to 9M$.
Carter's Constant ($Q$)
A hidden Killing-tensor symmetry discovered by Brandon Carter (1968) that supplies the missing fourth conserved quantity, keeping 3D Kerr trajectories completely integrable without chaotic scattering.

1. The Boyer–Lindquist Metric and the Ergosphere

In Boyer–Lindquist coordinates $(t, r, \theta, \varphi)$ with $G = c = M = 1$, the Kerr line element takes the form:

Boyer–Lindquist Kerr line element ($G = c = M = 1$)
$$ds^2 = -\left(1 - \frac{2r}{\Sigma}\right)dt^2 - \frac{4ar\sin^2\theta}{\Sigma}\,dt\,d\varphi + \frac{\Sigma}{\Delta}\,dr^2 + \Sigma\,d\theta^2 + \frac{A\sin^2\theta}{\Sigma}\,d\varphi^2$$ $$\Sigma = r^2 + a^2\cos^2\theta, \qquad \Delta = r^2 - 2r + a^2, \qquad A = (r^2 + a^2)^2 - a^2\Delta\sin^2\theta$$

Setting $a = 0$ immediately gives $\Sigma = r^2$, $\Delta = r(r - 2)$, $A = r^4$, and zeroes out the off-diagonal $dt\,d\varphi$ cross term, recovering the Schwarzschild metric of Part 1.

Everything physically new in Kerr comes from two places in this metric: the off-diagonal cross term $g_{t\varphi} = -2ar\sin^2\theta/\Sigma$, which couples time to azimuthal rotation, and the separation between the roots of $g_{tt} = 0$ and $\Delta = 0$:

Static limit (ergosphere) vs. event horizon and frame-dragging rate
$$r_{\text{ergo}}(\theta) = 1 + \sqrt{1 - a^2\cos^2\theta}, \qquad r_+ = 1 + \sqrt{1 - a^2}, \qquad \omega(r, \theta) = -\frac{g_{t\varphi}}{g_{\varphi\varphi}} = \frac{2ar}{(r^2+a^2)^2 - a^2\Delta\sin^2\theta}$$

Inside the ergosphere ($r_+ < r < r_{\text{ergo}}(\theta)$), $g_{tt} > 0$: the coordinate vector $\partial_t$ becomes spacelike, forcing every timelike particle and photon to rotate with $d\varphi/dt > 0$. Toggle the gold Ergosphere pill in the simulation above to render this oblate static-limit surface directly around the shadow.

2. Carter's Constant & Polynomial Geodesics in Mino Time

In Schwarzschild spacetime, spherical symmetry meant every light ray lived in a flat 2D plane through the origin. In Kerr spacetime, frame dragging twists trajectories out of their initial plane—rays wobble across latitude $\theta$ while corkscrewing in azimuth $\varphi$.

Fortunately, Kerr spacetime possesses not just energy $E = -p_t$ and axial angular momentum $L_z = p_\varphi$, but also Carter's constant $Q$, associated with a rank-2 Killing tensor. Dividing by $E$, two conserved impact parameters label every pixel's ray from the moment it leaves the camera:

Conserved impact parameters per pixel ($\xi, \eta, K$)
$$\xi \equiv \frac{L_z}{E}, \qquad \eta \equiv \frac{Q}{E^2} = p_\theta^2 + \cos^2\theta\left(\xi^2\csc^2\theta - a^2\right), \qquad K \equiv \eta + (\xi - a)^2 \ge 0$$

Here $\xi$ measures impact distance perpendicular to the spin axis, while $\eta$ controls how high in latitude the photon can climb before turning back toward the equatorial plane.

If you write the first-order geodesic equations in affine parameter $\lambda$, both $\dot r$ and $\dot\theta$ are coupled through $\Sigma = r^2 + a^2\cos^2\theta$, and worse: $\dot\theta$ contains a $\xi^2\cot^2\theta$ pole singularity while both equations hit square-root singularities ($\pm\sqrt{R(r)}$ and $\pm\sqrt{\Theta(\theta)}$) at every radial or polar turning point. Two classic coordinate transformations eliminate both headaches at once:

  1. Mino time ($\tau$): Reparameterize the ray by $d\tau = d\lambda / \Sigma$. This absorbs the $\Sigma$ coupling factor so $r(\tau)$ and $\theta(\tau)$ evolve as completely independent 1D oscillators.
  2. Polar cosine ($\mu \equiv \cos\theta$): Since $d\mu/d\tau = -\sin\theta\,(d\theta/d\tau)$, multiplying the polar potential $\Theta(\theta)$ by $\sin^2\theta = 1 - \mu^2$ cancels the $\cot^2\theta$ singularity analytically, turning both potentials into smooth quartic polynomials $R(r)$ and $M(\mu)$!
Decoupled polynomial potentials and turning-point-free 2nd-order Mino ODEs
$$\left(\frac{dr}{d\tau}\right)^2 = R(r) = \big(r^2 + a^2 - a\xi\big)^2 - \Delta K, \qquad \left(\frac{d\mu}{d\tau}\right)^2 = M(\mu) = \eta - (\eta + \xi^2 - a^2)\mu^2 - a^2\mu^4$$ $$\frac{d^2 r}{d\tau^2} = \frac{1}{2}R'(r) = 2r\big(r^2 + a^2 - a\xi\big) - (r - 1)\,K, \qquad \frac{d^2\mu}{d\tau^2} = \frac{1}{2}M'(\mu) = -(\eta + \xi^2 - a^2)\mu - 2a^2\mu^3$$ $$\frac{d\varphi}{d\tau} = \frac{a(2r - a\xi)}{\Delta} + \frac{\xi}{1 - \mu^2}$$

Look at $\ddot r(\tau)$ and $\ddot\mu(\tau)$: there are no square roots and no divisions in either acceleration! Radial bounces and polar oscillations pass smoothly through $v_r = 0$ and $v_\mu = 0$ under a 2nd-order midpoint step, and testing whether a ray crossed the equatorial accretion disk is just checking whether $\mu = \cos\theta$ changed sign (prevMu * mu <= 0.0).

💻 GLSL: Per-Pixel Mino-Time Kerr Geodesic IntegratorGLSL 
// Conserved impact parameters (xi = L_z/E, eta = Q/E^2) from camera ray (nr, nTheta, nPhi)
float bScale = sqrt(r0 * r0 + a2);
float xi     = -bScale * sin0 * nPhi;
float dth0   = bScale * nTheta;
float cot0   = cos0 / max(sin0, 1e-4);
float eta    = dth0 * dth0 + xi * xi * cot0 * cot0 - a2 * cos0 * cos0;
float K      = max(eta + (xi - a) * (xi - a), 0.0);

// Exact initial Mino velocities satisfying vr^2 = R(r0) and vmu^2 = M(mu0)
float Delta0 = r0 * r0 - 2.0 * r0 + a2;
float R0     = (r0 * r0 + a2 - a * xi) * (r0 * r0 + a2 - a * xi) - Delta0 * K;
float vr     = sign(nr) * sqrt(max(R0, 0.0));
float vmu    = -sin0 * dth0;

for (int i = 0; i < MAX_STEPS; i++) {
  if (r <= uRp * 1.015) { captured = true; break; }
  if (r > 30.0 && vr > 0.0) break;

  float distH = max(r - uRp, 0.025);
  float ds    = clamp(0.15 * distH, 0.018, 0.85);
  float dtau  = ds / (r * r + a2);

  // Polynomial Mino accelerations: d2r/dtau2 = 0.5*R'(r), d2mu/dtau2 = 0.5*M'(mu)
  float accR1  = 2.0 * r * (r * r + a2 - a * xi) - (r - 1.0) * K;
  float accMu1 = -(eta + xi * xi - a2) * mu - 2.0 * a2 * mu * mu * mu;

  float vr_mid  = vr  + 0.5 * accR1  * dtau;
  float vmu_mid = vmu + 0.5 * accMu1 * dtau;
  float r_mid   = max(r + 0.5 * vr * dtau, uRp * 1.005);
  float mu_mid  = clamp(mu + 0.5 * vmu * dtau, -0.9995, 0.9995);

  float accR_mid  = 2.0 * r_mid * (r_mid * r_mid + a2 - a * xi) - (r_mid - 1.0) * K;
  float accMu_mid = -(eta + xi * xi - a2) * mu_mid - 2.0 * a2 * mu_mid * mu_mid * mu_mid;

  r   += vr_mid  * dtau;
  vr  += accR_mid  * dtau;
  mu  = clamp(mu + vmu_mid * dtau, -0.9998, 0.9998);
  vmu += accMu_mid * dtau;

  // Azimuthal twist including frame-dragging term (negated for backward ray integration)
  float Delta_mid = max(r_mid * r_mid - 2.0 * r_mid + a2, 0.025);
  float sin2_mid  = max(1.0 - mu_mid * mu_mid, 0.002);
  phi -= ((a * (2.0 * r_mid - a * xi)) / Delta_mid + xi / sin2_mid) * dtau;
}

3. Why Gargantua's Shadow is D-Shaped

Click the D-Shadow Guide pill in the controls above and drag the Kerr Spin · $a/M$ slider from $0.000$ up to $0.998$. At $a = 0$, the black shadow sits centered inside the dashed cyan circle of radius $\sqrt{27}M \approx 5.196M$. As you spin the hole up, the entire black silhouette slides to the right (receding side) while its left (prograde, approaching) edge flattens into a vertical gold chord.

Why? In a non-rotating hole, unstable photon orbits all live at a single sphere $r_{\text{ph}} = 3M$. In Kerr spacetime, frame dragging helps photons orbiting with the spin stay aloft much closer to the horizon, while photons fighting against the spin must stay much farther out. The single photon sphere splits into a continuous shell of spherical photon orbits bounded by the prograde ($r_{\text{ph}}^-$) and retrograde ($r_{\text{ph}}^+$) equatorial photon rings:

Prograde & retrograde photon-sphere radii and Bardeen's critical curve ($R(r_0) = R'(r_0) = 0$)
$$r_{\text{ph}}^\mp = 2M\left[1 + \cos\!\left(\frac{2}{3}\arccos(\mp |a|)\right)\right] \quad\Longrightarrow\quad r_{\text{ph}}^- \xrightarrow{a\to 1} 1M, \quad r_{\text{ph}}^+ \xrightarrow{a\to 1} 4M$$ $$\xi_c(r_0) = \frac{r_0^2(r_0 - 3) + a^2(r_0 + 1)}{a(1 - r_0)}, \qquad \eta_c(r_0) = \frac{r_0^3\big(4a^2 - r_0(r_0 - 3)^2\big)}{a^2(1 - r_0)^2}, \qquad r_0 \in [r_{\text{ph}}^-, r_{\text{ph}}^+]$$

Each spherical photon orbit radius $r_0 \in [r_{\text{ph}}^-, r_{\text{ph}}^+]$ contributes one point $(\alpha, \beta) = (-\xi_c/\sin\theta_0,\, \pm\sqrt{\Theta(\theta_0)})$ on the shadow rim. As $a \to 1$, the prograde family collapses onto $r_0 \to 1M$ with constant $\xi_c = +2M$, carving the vertical flat side of the "D".

4. Prograde vs. Retrograde Disks & Cunningham's Exact Redshift

In Part 1, our "Spin" slider was an honest stand-in that moved the inner disk edge manually inside a static Schwarzschild metric. Here, the inner edge of the accretion disk is locked to the exact general-relativistic Innermost Stable Circular Orbit ($r_{\text{ISCO}}$) derived by Bardeen, Press, and Teukolsky (1972):

Exact Kerr ISCO radius (Bardeen–Press–Teukolsky 1972)
$$Z_1 = 1 + (1 - a^2)^{1/3}\left[(1 + a)^{1/3} + (1 - a)^{1/3}\right], \qquad Z_2 = \sqrt{3a^2 + Z_1^2}$$ $$r_{\text{ISCO}}(a) = M\left(3 + Z_2 \mp \sqrt{(3 - Z_1)(3 + Z_1 + 2Z_2)}\right) \quad \begin{cases} \to 1M & \text{prograde }(-) \\ \to 9M & \text{retrograde }(+) \end{cases}$$

Click the Prograde (+a) pill in the controls to flip the disk to Retrograde (-a) at $a = 0.990$: the live telemetry HUD immediately jumps from $r_{\text{ISCO}} = 1.454M$ to $r_{\text{ISCO}} = 8.972M$, opening a cavernous dark gap around the shadow while reversing the Doppler-bright side of the disk.

Even better, computing the exact general-relativistic frequency shift $g = \nu_{\text{obs}}/\nu_{\text{emit}}$ (combining gravitational redshift, frame dragging, transverse orbital dilation, and directional Doppler beaming) is simpler in Boyer–Lindquist coordinates than in approximate Euclidean vectors. Because disk plasma on an equatorial circular orbit has four-velocity $u^\mu = u^t(1, 0, 0, \Omega_K)$ with Kerr Keplerian angular velocity $\Omega_K = \frac{\pm 1}{r^{3/2} \pm a}$, and because $E = -p_t$ and $L_z = p_\varphi$ are conserved all the way to the camera, Cunningham's (1975) exact formula evaluates $g$ directly from the crossing radius $r$ and the ray's conserved $\xi = L_z/E$:

Exact Kerr Keplerian angular velocity and Cunningham (1975) frequency shift factor $g(r, \xi)$
$$\Omega_K(r) = \frac{\pm 1}{r^{3/2} \pm a}, \qquad g(r, \xi) = \frac{-p_\mu u^\mu_{\text{obs}}}{-p_\nu u^\nu_{\text{emit}}} = \frac{1}{u^t\big(1 - \Omega_K\,\xi\big)} = \frac{\sqrt{1 - \frac{3}{r} \pm \frac{2a}{r^{3/2}}}}{1 + \frac{\pm a - \xi}{r^{3/2}}}$$

Notice that at $a = 0$, this reduces to $g = \frac{\sqrt{1 - 3/r}}{1 - \xi/r^{3/2}}$, automatically avoiding the double-counted $\sqrt{1-\beta^2}$ pitfall of naive special-relativistic decompositions.

💻 GLSL: Exact Cunningham Kerr Redshift & Disk EmissionGLSL 
// Equatorial disk plane crossing (prevMu * mu <= 0.0) at Boyer-Lindquist radius rr
if (rr > uIsco && rr < rOutDisk) {
  float x      = rr / uIsco;
  float tprof  = pow(x, -0.75) * pow(max(1.0 - inversesqrt(x), 0.0), 0.25) * 2.05;
  float Tb     = uTpeak * tprof;

  float edge = max(1.0 - sqrt(uIsco / max(rr, uIsco + 1e-4)), 0.0);
  float flux = 2.6 * pow(uIsco / rr, 2.85) * pow(edge, 0.78);

  // Signed spin (+a for prograde, -a for retrograde) and r^(3/2)
  float aSigned = a * uPrograde;
  float r32     = rr * sqrt(rr);

  // Exact Cunningham (1975) Kerr redshift + frame-dragging + Doppler factor g(r, xi)
  float numG = sqrt(max(1.0 - 3.0 / rr + 2.0 * aSigned / r32, 0.015));
  float denG = max(1.0 + (aSigned - xi * uPrograde) / r32, 0.14);
  float g    = mix(numG / max(1.0 + aSigned / r32, 0.2), numG / denG, uDoppler);

  vec3 bb = blackbody(Tb * clamp(g, 0.28, 3.4));
  bb = mix(bb, bb * vec3(1.14, 0.99, 0.75), 0.58);
  accumCol += T * bb * flux * pow(clamp(g, 0.08, 18.0), 3.6) * (0.32 + 1.32 * n) * 6.2 * env;
}

5. The Kerr Render Pipeline

Putting the pieces together, each pixel converts its view-plane ray into conserved Carter constants $(\xi, \eta)$, integrates the polynomial Mino-time equations across the ergosphere and equatorial plane, and accumulates exact Cunningham-shifted blackbody emission:

Camera ray (ξ = Lz/E, η = Q/E²) Mino Verlet integrator r̈ = ½R'(r), μ̈ = ½M'(μ) Ergosphere & φ frame-dragging twist Cunningham disk r_ISCO(a) · g(r, ξ) ACES + Bloom
Figure 1: Replacing affine parameter $\lambda$ with Mino time $\tau$ decouples the Kerr geodesic into polynomial radial and polar oscillators with exact conserved impact parameters $(\xi, \eta)$.

6. Schwarzschild vs. True Kerr: What Changed

Comparing the static spacetime of Part 1 against the spinning Kerr spacetime built here highlights how much physics hides inside that single spin parameter $a/M$:

Feature Part 1: Schwarzschild ($a = 0$) Part 3: Boyer–Lindquist Kerr ($a \to 1$)
Event horizon $r_+$ Fixed at $r_s = 2M$ Shrinks with spin: $r_+ = M(1 + \sqrt{1 - a^2}) \to 1M$
Ergosphere None ($r_{\text{ergo}} = r_+$ everywhere) Oblate static-limit shell $r_{\text{ergo}}(\theta) = M(1 + \sqrt{1 - a^2\cos^2\theta})$
Ray plane Planar great-circle motion ($\vec h = \vec r \times \dot{\vec r}$) 3D non-planar twisting governed by Carter's constant $Q$
Shadow silhouette Centered circle of radius $\sqrt{27}M \approx 5.20M$ Shifted, prograde-flattened Bardeen D-curve $(\xi_c, \eta_c)$
ISCO radius $r_{\text{ISCO}}$ Fixed at $6M$ ($3\,r_s$) Shrinks to $1M$ (prograde $+a$) or expands to $9M$ (retrograde $-a$)
Frequency shift $g$ Local special-relativistic + static redshift Exact Cunningham factor $g(r,\xi) = [u^t(1 - \Omega_K\xi)]^{-1}$

Try clicking the Gargantua (a=0.998) and Schwarzschild (a=0) preset buttons back-to-back in the live simulation at the top of the page—watching the horizon contract, the prograde disk plunge inward to $1.24M$, and the shadow snap into Bardeen's D-profile brings five decades of general-relativistic black hole optics into a single browser tab.