The physics behind a black hole shader

Back in August 2023 I had to submit a project for an elective WebGL programming course in university. I chose to create an interactive and animated black hole demo. The primary idea was to show how the gravitational pull of black holes “bends” light and how the result looks from different points of view.

Then in the December of that same year, I sat down and started writing a blog post (as you may notice from the “Created” date), explaining how some of the more complicated parts were made. Funnily enough, I started finding errors… a lot of errors.

As it turned out, my implementation was completely wrong, but somehow I fudged the formulas around in just the right way to make it look correct. So the post got delayed and sparsely rejuvenated over the following months.

For this reason I chose to take up two related electives this year: one on general relativity and one on ray tracing, with the hopes of gaining enough understanding of the material to rework the original solution into something proper.

This article represents the first of a two-part series about this very shader. Here, we’ll talk about physics and how to generally model them. In the second part, we’ll explore ray tracing and computer graphics with WebGL.

I’ll assume you don’t know much about physics, but are familiar with mathematics and at the very least some linear algebra.

What do the physics say?

The “big words” version of what we’re going to explore is:

Using Einstein’s field equations, we can precisely calculate the curvatures in space(time) at every point around a black hole. Obviously we’re going to use Schwarzschild metric on a (Schwarzschild) black hole. To further simplify our calculations, we can make approximations using celestial mechanics.

This is not necessary to understand, but it will make our goal much clearer.

Now, let’s try and unravel everything.

What is gravity?

Everything with mass applies a “pulling” force to everything else with mass. It’s easy to see this is true, but we get into trouble when trying to find a formula which predicts/calculates that force.

Newton gave it a good shot, with his 3 laws of motion, which proved to be correct for about a century and a half. They even helped the discovery of Neptune, since measurements showed that the orbit of Uranus didn’t meet expectations if there wasn’t another planet meddling with things. The real crux of the issue came when the same problem occurred with Mercury’s orbit, but no new planet was found.

As it turned out, Newton’s formulae had only been approximations, they worked when light (and therefore causality, or the speed at which changes propagate) had infinite speed.

Albert Einstein’s general theory of relativity showed that gravity wasn’t a “force which pulled”, it was actually an effect which bent space. And that only happened because the (constant) speed of light wasn’t infinite.

What does it mean to bend space?

Think of a flat surface, like a table. What is the shortest distance between two points? Obviously, that’s a straight line segment. If you then extend that segment, so it has an infinite length, it becomes a line¹.

However, what if our surface isn’t flat? Imagine two points on a sphere, what is the shortest distance between them? Obviously, it’s an arc of a circle and extending it to it’s maximum possible length, we get a great circle.

Image by CheCheDaWaff on Wikipedia

Now, what does it mean to go “forwards”? On a flat surface, we can define it as “following a (straight) line segment”, on our sphere we can say “following an arc of a great circle”.

Better yet, in geometry there are geodesics, the general term for a curve on any surface which is the shortest² path between any two points. So now we can say “going forwards means going along a geodesic”.

It’s important to note that an arbitrary geodesic doesn’t have most properties you might be used to. For example, on a flat surface, we can have parallel lines, but on our sphere, great circles will always intersect³.

This is what gravity does, it doesn’t pull or push, it bends, and more precisely, it bends geodesics. Furthermore, we live in three spacial dimensions, so we can move in three different directions.

The resultant effect is that if you travel perfectly forwards in a spaceship and pass by an object with mass (like a planet) from a close enough distance, your path will be curved, without you steering your ship.

Although not yet mentioned, you need to know that photons (light) don’t have mass! In effect, in a universe where there is only light bouncing around, Newton’s equations still work.

On that note, the point that light doesn’t travel at an infinite speed is significant. When the opposite is true, General Theory of Relativity collapses into the Newtonian equations!

Going back on track, Einstein gives us the necessary equations for those geodesics, the Einstein Field Equations.

How to handle the EFEs?

Although they’re called equations the result/output we get from solving them isn’t a number! It’s a function, which more or less takes the parameters of objects in the universe and returns how they will all behave. To be more precise, the output isn’t really a function, it’s a “metric”, which is a distinctly different mathematical object.

In the 1910s Schwarzschild solved the equations for a universe where there is only one non-rotating⁴ particle with mass, then in the 1960s the equations were solved when that particle could rotate, that’s it. There are plenty of other mathematical solutions, but we have been able to relate the “real” world to only those two.

This may sound unsatisfactory, but in practice, it’s quite enough. In the solar system, >99.86% of mass is in the Sun, so assuming all other objects are massless, we get enough precision to explain Mercury’s orbit. We can apply the same idea to our simulation, leaving all objects which aren’t our black hole without mass.

Schwarzschild’s metric and accompanying black hole properties are the simplest possible solution, so we’ll use that.

What is a black hole?

Let’s imagine atoms as a set of glass marbles, scattered on a table:

Image from By Okkisafire on Wikipedia

Now, start pushing from all sides of the table to it’s center. In the beginning, you’ll mostly roll them, until they all get in contact with each other. At that point (assuming they all stay flat on the table), the glass will keep the marbles together, it’s strength opposing your force.

By increasing the force with which you push on them, after a certain threshold, the glass will break, turning the marbles into shards. Those shards will still have some fight in them, so after increasing your force by another level, you’ll turn all of those shards into tiny glass granules.

A similar occurrence happens with atoms. Normally atoms are spread apart, by pushing them together, at a certain point (ignoring the orbiting electrons) all nuclei will be “touching” and the nuclear force will keep them together. By further increasing your pushing force you can overwhelm the nuclear force and “shatter” the nuclei into their subatomic particles.

Assuming an enormous explosion doesn’t happen, and you can keep increasing your force, it’s also possible to overwhelm the “subatomic forces” that keep those subatomic particles “together”⁵. However, those subatomic particles are more commonly referred to as elementary particles, meaning there are no “granules” into which the elementary particles can break down.

It’s not possible for those particle to be completely unbreakable, since that would mean they could theoretically sustain infinite amounts of force. Such an assumption is unthinkable.

Then that leaves us with only one solution: laws of matter “break”. If nothing can keep the elementary particles “together”, so keep matter “together”, then the result is a sort of “collapse”, all of that matter just “deflates” into a single, infinitely small point (but the amount of mass isn’t infinite, it’s the sum of mass from the particles that collapsed). That point is called a black hole.

This object is very bizarre, generally it doesn’t have too much properties, but the ones it has are really unique. In our case, the only important part is it’s effect on gravity, and more precisely, light.

Since a black hole is a single infinitely small point, then all of it’s mass is concentrated on that point, making it alter the geodesics in a relatively small area of space a lot. So much so, if you get too close, you will not be able to overwhelm gravity and “merge” into that infinitesimal point:

The white lines are our geodesics
Image by OpenStax University Physics on Wikipedia

A fascinating realization is that light always travels along geodesics (“forwards”), therefore this extreme bending allows us to view otherwise obstructed objects. We call that effect gravitational lensing. Though not unique to black holes, that’s when the effect is strongest:

Notice how the gray arrows show the path light travels from a galaxy that is behind (obstructed by) the black hole

I’ll abstain from explaining the specifics of light, but you may think of it as a small particle, “similar” to an electron, travelling through space. In fact, Newton understood light to be just that, and called the underlying particle a “corpuscle”. This approximation will serve us well.

What about speed and time?

Two of the most significant things Einstein’s theory shows, which I’ve left out thus far, is that time is just another point in space (a fourth “spacial” direction in a way) and everything is relative. Which also signifies that our geodesics are actually four dimensional, or in other words, gravity also bends the observed passage of time!

This is relevant for light, since it’s speed is constant and it always travels along a geodesic, but to an outside observer it doesn’t look that way. However, in our simulation, we only care about how thing look, so it will make our live much simpler if we try to ignore reality wherever we can.

Therefore, to calculate by how much light will move and change, we’ll start with the EFEs and then simplify and approximate until we get something useable.

Finding a formula

Let’s get down to business, our starting point, the Schwarzschild metric, in Schwarzschild coordinates, looks like this:

YOUCH!

The variables themselves aren’t that bad: c is the speed of light, τ is time, Ω, t and r are Schwarzschild coordinates.

However, as you can see, there are two big problems: we have to calculate derivatives (that’s what the d signifies) and convert from Schwarzschild coordinates to normal (euclidean) coordinates, that the computer understands “natively”. This would be too slow to render in real time.

Thankfully Riccardo Antonelli was able to do all of our coordinate conversions and use Newton’s ideas to simplify the living crap out of that monster, turning it into something digestible:

Explaining his steps and the formula is way beyond the scope of my understanding and this post. He explains it in more mathematical detail here.

But now, finally, we can get a plan for our simulation:

• represent light as particles
• move each particle in a direction and amount specified by the EFE-based formula
• rinse and repeat until the light hits something

Question now is: what do all of those letters represent? This turned out to be a much harder problem than it should have, I actually had to refer Ricardo Antonelli’s and Dmitry Brant’s implementation in their projects, alongside copious articles on Wikipedia to fully grasp how to put it to work.

Applying the formula to “bend” light’s path

The direction, alongside the amount by which the particle moves, can be represented with a single vector, which we’ll call light’s velocity.

Referring to the formula again:

F→ returns a force vector (F for “force” and ^→ for “vector”): any force acting on an object has a strength (magnitude) and a direction. Imagine you are pushing a box: you are applying some force to it in a “forward” direction. If you pull the box, you’ll apply some force in a “backward” direction.

This force vector conveys how much and in what direction exactly the black hole is “pulling” on the corpuscle.

Although the function takes in only one argument, r, that doesn’t mean that h or r^ aren’t also arguments:

• r is the distance from the center of the black hole (the singularity) to the current position of the corpuscle.

• h is the magnitude of the “specific angular momentum” of the light “particle”.

  “Normal” linear momentum specifies how much “force” at what speed it takes to stop (or change) the speed and direction of an object in linear motion (that’s why it’s in kg*m/s). Imagine a big heavy train on a perfectly straight rails going at a fast speed. How do you forcefully stop it?

  You can put a really heavy barrier on the tracks, in front of the train. Since the train is moving at a really great speed, this barrier will have to be much much heavier (with much more mass) than the train itself.

  Alternatively, you can “launch” something against the train, meaning you can put a second train, moving in the opposite direction of the first train, so they hit head-on. Since the second train would also be moving, if it’s travelling at the same speed and has the same mass as the first train, upon collision they would both stop⁶. If the second travels at twice the speed, then it can have half of the mass of the first to produce the same effect.

  The property with which we can “prove” why both a larger mass at no speed, and a smaller mass at a faster speed, can stop the train, is the momentum.

  Angular momentum specifies how much “force” at what speed it takes to stop (or change) the speed and direction of a rotating object. Imagine the train scenario but on a circular track.

  The specific angular momentum is the angular momentum of an object, divided by it’s mass. It’s magnitude kind of tells us which property is more prevalent in the momentum, the mass or the speed. For the same angular momentum, a long (big magnitude) specific angular momentum means our object is small and moves fast, while for a short (small magnitude) specific angular momentum out object is large and moves slowly.

  Another simplification we need to is make our black hole stationary! Then, rather than finding the angular momentum⁷ and dividing by the mass⁸, we’ll use another formula:

  

  where in this specific formula only r is the vector from the singularity to the corpuscle and v is light’s velocity.

  The cross product returns a vector, so we’ll need to take the magnitude of the result and plug that into the formula. It’s also important that we calculate it only once, in the beginning. So, we shouldn’t be calculating the magnitude of the specific angular momentum every time the light “particle” is moved, rather we calculate it once, after it gets close enough to be affected by the black hole.

• r^ is the normalized (unit, hence the ^ on top) vector from the center of the black hole to the light particle:

  

Finally, putting all needed variables together in one proper function, we get:

As you may notice, this can be simplified, if we substitute r^:

In effect, we pass the distance from the black hole to the corpuscle and the function returns how light’s velocity changes.

Finally, we know what to give the formula, and what it will return. The task at hand is making an algorithm which utilises it. For that we’ll need to get even more familiar with black holes.

Properties of a black hole

To summarize, we’re working with a stationary Schwarzschild black hole. Meaning our black hole has only one property: mass⁹. Everything else can be derived from it.

There are three radii, which will be the backbone of our black hole: the event horizon, the photon sphere and ISCO.

The event horizon, also called the Schwarzschild radius, is the distance from the singularity, which once passed, nothing can return and escape. The formula for it is:

where G is the gravitational constant and M is the mass of the black hole (in kg).

The photon sphere is the distance from the singularity, in which photons start orbiting the black hole (and can either escape or pass the event horizon). Formula:

where G and M are the same as before, c is the speed of light in vacuum and rs is the aforementioned Schwarzschild radius.

Finally, the innermost stable circular orbit (ISCO), formally it’s the distance, before which you can orbit the black hole, but once passed, you’ll get sucked in and either get shot out or be forever lost past the event horizon. For our purposes, it’s fine to approximate it as the maximum radius after which space-time isn’t bent¹⁰ and objects aren’t affected. Formula:

Putting them all together, we should expect light to travel like this: Animation by Hugo Spinelli on Wikipedia

Light-bending algorithm

Finally, now we have everything we need to devise an algorithm for our simulation.

• For light we’ll store two vectors: one for it’s position and one for it’s velocity.
• For the black hole we’ll store a vector and a quantity: the position of the black hole’s center and it’s mass. For convenience, we could just store the Schwarzschild radius, since all of our radii can be expressed only by it.

Once the distance between a corpuscle and the center is less than or equal to the ISCO, we can begin augmenting light’s velocity. On small steps, we move the light forward a bit with the velocity and then update the velocity with the black hole’s pull.

If light enters the event horizon, then we know it cannot escape, and stop moving it. If it passes the ISCO then it has left the black hole and will never approach it again.

If light enters the photon sphere, then we’ll have to “track” how long it’s been there, if it decides to orbit in a stable manner, then it will stay there indefinitely. The simplest way to “track” this is just to limit how many times we can move the light forward.

An algorithm doing this would look a bit like this (in pseudocode, inspired by Python):

as a reminder, here is the definition of F→:

That’s it! We have all of our theoretical background, formulae and algorithm to start work on our actual program.

To Be Continued!

Appendix

LaTeX source for equations

The images for my custom LaTeX equations were created on latexeditor.lagrida.com, by using Firefox’s built-in ability to screenshot HTML nodes, in particular the MathJax node with the equation.

Antonelli’s formula - explained

\begin{align*}
    & \overrightarrow{F}(r_l,h) = -\frac{3}{2} \; h^2 \; \frac{\hat{r_l}}{\|r_l\|^5} \quad \text{, where}\\
    & h = \|r_l \cross v_l\|\\
    & \hat{r_l} = \frac{r_l}{\|r_l\|}\\
    & \|x\| = \text{magnitude of }x\\
    & r_l = \text{vector from black hole center to light "particle"}\\
    & v_l = \text{velocity of light "particle"}
\end{align*}

Antonelli’s formula - simplified

\begin{align*}
    & \overrightarrow{F}(r_l,h) = -\frac{3}{2} \; h^2 \; \frac{r_l}{\|r_l\|^6} \quad \text{, where}\\
    & h = \|r_l \cross v_l\| \text{ (calculated right before being affected by black hole)}\\
    & \|x\| = \text{magnitude of }x\\
    & r_l = \text{vector from black hole center to light "particle"}\\
    & v_l = \text{velocity of light "particle"}
\end{align*}

Algorithm

\begin{align*}
    & p_{cbh} = \text{position of center of black hole}\\
    & r_s = \text{Schwarzchild distance from } p_{cbn} \\
    & r_{ps} = \frac{2}{3}r_s \quad \text{(photon sphere)}\\
    & r_{ISCO} = 3 r_s \\
    \\
    & p_l = \text{position of light "particle"}\\
    & v_l = \text{velocity of light "particle"}\\
    \\
    & r_l = \| p_l-p_{cbh} \| \quad \text{(distance between black hole and light)}\\
    & \text{If } r_l \leq r_{ISCO} : \\
    & \quad h = \| r_l \cross v_l \|\\
    & \quad \text{LightMovesLimit} = \infty \\
    & \quad \text{Repeat maximum LightMovesLimit times}:\\
    & \qquad \text{If } r_l \leq r_s : \\
    & \qquad \quad \text{return "Light never escapes"} \\
    & \qquad \text{If } r_l \leq r_{ps} : \\
    & \qquad \quad \text{LightMovesLimit} = 100 \\
    & \qquad \text{If } r_l > r_{ISCO} : \\
    & \qquad \quad \text{return "Light escaped black hole"} \\
    & \qquad p_l = p_l + v_l \quad \text{(move light)}\\
    & \qquad r_l = \| p_l-p_{cbh} \| \quad \text{(update distance to $p_{cbh}$)}\\
    & \qquad v_l = v_l + \overrightarrow{F}(r_l,h) \quad \text{(update light velocity)} \\
    & \quad \text{return "Light stuck in photon sphere, never escapes"}
\end{align*}