How to render a black hole

The image in the header is not a texture, a video, or a 3D model. It's the output of a fragment shader that takes every pixel on the screen, shoots an imaginary photon backwards into the scene, and asks: where would this photon have come from if light were bending around a one solar mass of curvature here?

That's a lot of work for one pixel. We do it sixty times a second, for two million pixels at once, on a laptop GPU, with no precomputed lookup tables.

You can poke every step. Drag inside any sandbox to look around.

Starting from nothing: straight rays

Before we add a black hole, here's the baseline. Each pixel shoots a ray into space. The ray goes in a straight line. It samples a procedural starfield: bright sparse stars, a denser field underneath, and a faint nebula on top of that. No physics, just trigonometry and a hash function.

Now we curve those rays. From here on, everything we add is just a question of what the ray hits along the way.

Gravity, by itself

Massive objects curve spacetime. Light follows the shortest path through that curvature, called a geodesic. Far from a mass, geodesics look like straight lines. Close to a mass, they bend. Close enough to a black hole, they wrap.

Before any of that becomes a fragment shader, it helps to see geodesics in a flat top-down view. Each line in the playground is a photon launched from the left edge, aimed straight to the right. Without gravity, they'd all be parallel. With gravity, the ones that pass closer to the central mass curve inward. The closest ones are captured.

The bend angle is set almost entirely by the impact parameter (how close the ray came to the mass), not by where the ray started. And the transition between escape and capture is sharp: rays just outside the photon sphere swerve and escape, rays just inside it spiral in. That cliff is what creates the hard black silhouette around a real black hole.

The equation we'll integrate

For a Schwarzschild (non-rotating, uncharged) black hole, a photon's spatial trajectory obeys

where x is the photon's position relative to the hole, r = |x|, RS is the Schwarzschild radius (the event horizon), and L = |x × v| is the angular momentum. L is conserved along the geodesic, so we compute it once.

In code, we integrate this with velocity-Verlet: a symplectic, energy-conserving integrator that lets us take big steps far from the hole and small ones near it without the trajectory blowing up.

1// Schwarzschild geodesic step (in fragment shader)
2// d²x/dλ² = -1.5 · RS · L² / r⁵ · x   where L = |x × v|
3vec3 Lvec = cross(pos, vel);
4float L2 = dot(Lvec, Lvec);
5float gravCoeff = -1.5 * RS * L2;
6 
7for (int i = 0; i < 200; i++) {
8  float r = length(pos);
9  float h = 0.16 * clamp(r - 0.4 * RS, 0.06, 3.5);  // adaptive step
10 
11  // a(x) = gravCoeff / r⁵ · x
12  float invR5 = 1.0 / (r*r*r*r*r);
13  vec3 acc  = (gravCoeff * invR5) * pos;
14 
15  // Velocity-Verlet — symplectic, conserves energy
16  vec3 p1 = pos + vel * h + 0.5 * acc * h * h;
17  float invR15 = 1.0 / pow(length(p1), 5.0);
18  vec3 acc1 = (gravCoeff * invR15) * p1;
19  vec3 v1 = vel + 0.5 * (acc + acc1) * h;
20 
21  if (length(p1) < RS * 0.35) { absorbed = true; break; }
22  pos = p1; vel = v1;
23}

Replace the straight ray from Step 01 with this loop and you get gravitational lensing for free. No disk yet, just stars seen through curved spacetime.

Pull the step count down to 30 and you'll see the integration break: light begins to wobble and stars near the rim shimmer as the integrator skips important geometry. 200 steps is what we use in production. Adaptive step size means it only pays that cost near the hole, and most rays terminate in 20 or 30 iterations anyway.

A disk of glowing plasma

Real black holes aren't lit. They're black. What we see in famous images, and what we'll render here, is the accretion disk: a thin layer of plasma spiralling inward, heated by friction, glowing in the X-ray and visible band before it falls past the event horizon.

The temperature isn't uniform. Plasma orbiting closer to the hole is moving faster, has lost more potential energy, and radiates more brightly. Novikov and Thorne worked out the radial profile in 1973:

The r_isco is the innermost stable circular orbit. Inside it, plasma can't hold an orbit and plunges straight into the hole. For a non-spinning Schwarzschild hole, r_isco = 3·RS. Plasma at exactly the ISCO has zero temperature (the second factor goes to zero); the disk peaks slightly outside it and falls off as r^(−3/4) going out.

Temperature is then mapped to color through a blackbody approximation. Hot things glow blue, moderate things glow white, cold things glow red. The shader doesn't compute Planck's law from first principles. That would be expensive and wrong-looking under tone mapping anyway, so instead it uses a hand-fit four-color gradient that lines up with the real curve in the visible range.

The full disk shader combines the temperature profile with a few extra layers: FBM turbulence oriented along Keplerian streamlines (so it flows correctly), concentric ring modulation for visible structure, and a soft inner/outer fade. Here's the temperature half on its own, without any gravitational lensing yet:

1// Novikov-Thorne temperature profile (1973)
2// T(r) ∝ r^(-3/4) · (1 − √(r_isco / r))^(1/4)
3float xr = ISCO / r;
4float tProfile = pow(ISCO / r, 0.75)
5               * pow(max(0.001, 1.0 - sqrt(xr)), 0.25);
6 
7// Gravitational redshift: light from deep in the well loses energy
8float gRedshift = sqrt(max(0.01, 1.0 - RS / r));
9tProfile *= gRedshift;
10 
11vec3 col = blackbodyColor(tProfile);

Wrapping the disk around the hole

Now we put the two pieces together. Same geodesic loop as Step 02, but each time a ray crosses the disk's plane (the y=0 equator) we sample its color and accumulate it. Light from the far side of the disk gets bent upward by the hole's gravity, around the silhouette, and into our eye as a halo arching over the top of what should be empty space. That halo is what makes a black hole look like a black hole instead of a dark circle.

Lower the tilt to near zero and the disk becomes a flat band crossing the silhouette. Increase it and the top arc thickens. The disk isn't getting wider; you're seeing more of the back of the disk being lensed up into view.

Why one side is brighter

Compare the previous sandbox to the hero shader at the top of this article. They look almost the same. Almost. The hero is brighter on one side, and that asymmetry is the most photogenic bit of physics in the whole thing.

The plasma in the disk isn't sitting still. It's orbiting at a meaningful fraction of the speed of light. The side moving toward us is Doppler-boosted: its photons arrive more frequently, blue-shifted, and the relativistic beaming concentrates them in our direction. The side moving away is the inverse: red-shifted, dimmed, defocused.

The brightness boost goes as the Doppler factor cubed. Small velocity differences become large brightness differences.

1// Doppler beaming from Keplerian orbital motion.
2// orbDir is tangent to the orbit at this point on the disk.
3float orbSpeed = sqrt(0.5 * RS / max(r, DISK_IN));
4vec3  orbDir   = normalize(vec3(-hit.z, 0.0, hit.x));
5 
6float dop = 1.0 + 2.0 * dot(normalize(vel), orbDir) * orbSpeed;
7dop = max(0.15, dop);
8 
9// Brightness scales as Doppler^3 (relativistic beaming).
10// Frequency shifts blue-ward on approach, red on retreat.
11float boost     = dop * dop * dop;
12float colorTemp = tProfile * pow(dop, 1.8) * 1.2;
13vec3  col       = blackbodyColor(colorTemp) * intensity * boost;

The photon ring

If you stare at Step 05 with the camera tilted sharply, you'll notice a thin bright line right on the edge of the silhouette. That's the photon ring: light that orbited the hole once (or twice, or more) before escaping. Photons that pass right at the photon sphere can make a half-loop, three-quarters, even a full loop, picking up successive views of the disk on each pass.

In code this falls out naturally from one detail: a ray can cross the disk plane more than once. The first crossing is the disk in front of the hole. The second is the back of the disk, wrapped over the top. The third, fourth, fifth get progressively concentrated near r ≈ 1.5·RS, forming a stack of ever-thinner rings hugging the silhouette.

1// A ray can cross the y=0 disk plane MULTIPLE times after lensing.
2// Crossing 1: the disk in front of the hole.
3// Crossing 2: the back of the disk, wrapped over the top.
4// Crossing 3+: the photon ring — light orbiting near r=1.5·RS.
5if (pos.y * p1.y < 0.0 && diskAccum.a < 0.97) {
6  float t = pos.y / (pos.y - p1.y);
7  vec3 hit = mix(pos, p1, t);
8  vec4 dc = shadeDisk(hit, vel, time);
9 
10  // Crisp 1px photon ring would alias horribly; attenuate later passes.
11  if (diskCrossings >= 2) { dc.rgb *= 0.15; dc.a *= 0.15; }
12 
13  // Alpha-over compositing (front-to-back).
14  diskAccum.rgb += dc.rgb * dc.a * (1.0 - diskAccum.a);
15  diskAccum.a   += dc.a * (1.0 - diskAccum.a);
16  diskCrossings++;
17}

Making it 60fps on a laptop

Every pixel of a 1080p screen runs the whole ray-marching loop, sixty times a second. The worst-lit pixels (the rays near the hole) do up to 200 Verlet steps each, which works out to around 124 million steps per second for the dense pixels alone.

The biggest savings come from variable step size. Far from the hole, geodesics are nearly straight and we coast in strides as large as h = 3.5 in scene units. Close in, curvature is steep and we drop to h ≈ 0.06. The expression h = 0.16 · clamp(r − 0.4·RS, 0.06, 3.5) picks the size each step. Almost all rays spend almost all of their iterations far from the hole, so the average step is large even though the worst-case step is small.

Loop invariants are pulled out: angular momentum L² and the constant factor in the acceleration are computed once per ray. Inside the loop the only divisions are 1/r² and 1/r⁵, and the latter is reused as (invR²)² / r to save a multiply. Small savings, two million pixels, a hundred iterations each.

Termination is aggressive. A ray outside r = 25 and still moving outward isn't coming back, so we break. A ray past r = 55 is gone regardless. A ray inside 0.35·RS has crossed the event horizon; we mark it absorbed and break. Most pixels exit in fewer than 50 of the 200 nominal iterations.

Putting it all together

The full shader stitches every step together: star field for the sky, geodesic integration for the bending, Novikov-Thorne plus blackbody for the disk's color, alpha-over compositing for multiple crossings, Doppler beaming for the asymmetry, filmic tone curve at the end. About 560 lines of HTML total, half of which is the GLSL fragment shader. Drop it in an <iframe> anywhere and it just runs.

Everything is turned on here, every parameter exposed. Try the chromatic dispersion slider: it's a spectral hue mapping on the disk we didn't cover in the article, but it's the same blackbody mapping with an extra dimension.

Where to go from here

Two papers and one implementation, in escalating commitment. rantonels' Starless is the more rigorous offline cousin of this shader, with full Kerr support and beautiful explanatory diagrams. The Wikipedia article on Schwarzschild geodesics derives the equation we used from the metric tensor in three pages of pleasant tensor calculus. And the Interstellar look is documented in Kip Thorne's Gravitational Lensing by Spinning Black Holes in Astrophysics: roughly this article, plus a year of work and a Nobel laureate.

The source is under 600 lines of HTML. View it, copy it, save it. The whole file works as a standalone background, hero, or wallpaper without any build step.