path: root / src / blog / the-physics-behind-a-black-hole-shader / index.du

'def title The physics behind a black hole shader
'def  desc A deep dive into some of the physics (and mathematics) required for implementing a black hole shader

# {title}

{pagenav 22 min;12 December 2023;06 October 2024}

:contents
:numbered
1. {uhref What do the physics say?}
   .numbered
   1. {uhref What is gravity?}
   2. {uhref What does it mean to bend space?}
   3. {uhref How to handle the EFEs?}
   4. {uhref What is a black hole?}
   5. {uhref What about speed and time?}

2. {uhref Finding a formula}
3. {uhref Applying the formula to “bend” light’s path}
   .numbered
   1. {uhref Properties of a black hole}

4. {uhref Light-bending algorithm}
5. {uhref Appendix}



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 {i 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 {i 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:

'def EFE https://en.wikipedia.org/wiki/Einstein_field_equations

Using {url Einstein's field equations;{EFE}}, we can precisely calculate the curvatures in space(time) at every point around a black hole.
Obviously we're going to use {url Schwarzschild metric;https://en.wikipedia.org/wiki/Schwarzschild_metric} on a (Schwarzschild) black hole.
To further simplify our calculations, we can make approximations using {url celestial mechanics;https://en.wikipedia.org/wiki/Celestial_mechanics}.

This is not {i necessary} to understand, but it will make our goal {b 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 {url 3 laws of motion;https://en.wikipedia.org/wiki/Newton%27s_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 {url Mercury's orbit;https://en.wikipedia.org/wiki/Mercury_%28planet%29#Advance_of_perihelion}, but no new planet was found.

As it turned out, Newton's formulae had only been approximations, they worked when light (and therefore {url causality, or the speed at which changes propagate;https://www.youtube.com/watch?v=msVuCEs8Ydo}) 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 {b bent} space.
And that only happened because the ({b 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{fn This is the mathematical definition of a line, obviously in conversation we refer to a line segment as a line. But this distinction will be useful.}.

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 {url great circle;https://en.wikipedia.org/wiki/Great_circle}.

{img ./img/Illustration_of_great-circle_distance.svg.png}{br}
{i Image by CheCheDaWaff on {url Wikipedia;https://commons.wikimedia.org/w/index.php?curid=48187293}}

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 {url geodesics;https://en.wikipedia.org/wiki/Geodesic}, the general term for a curve on any surface which is the shortest{fn This actually isn't completely true, a geodesic is the "locally shortest distance", I've skipped the details since they aren't really that important for our purposes. More info can be found on {url Wikipedia;https://en.wikipedia.org/wiki/Geodesic#Introduction}.} 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{fn You might think about the latitudinal lines on a globe, they are indeed parallel (equidistant), but if you take any two points on the same such line, the distance between them won't be the smallest.}.

This is what gravity does, it doesn't pull or push, it {i bends}, and more precisely, it {b bends geodesics}.
Furthermore, we live in three spacial dimensions, so we can move in three {i 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 {b photons} (light) {b 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 {b 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 {url Einstein Field Equations;{EFE}}.

### 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 {i more or less} takes the parameters of objects in the universe and returns {i how} they will all behave.
To be more precise, the output isn't {i 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{fn He also assumed the body didn't have an electric charge and the {url cosmological constant;https://en.wikipedia.org/wiki/Cosmological_constant} is zero.} 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, {url >99.86% of mass is in the Sun;https://en.wikipedia.org/wiki/Solar_System}, 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:

{img https://upload.wikimedia.org/wikipedia/commons/thumb/0/0e/Marbles_from_Indonesia.jpg/1280px-Marbles_from_Indonesia.jpg}{br}
{i Image from By Okkisafire on {url Wikipedia;https://commons.wikimedia.org/w/index.php?curid=30338734}}

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 {url nuclear force;https://en.wikipedia.org/wiki/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"{fn Technically, it's more accurate to say "matter apart", since after the collapse matter is as close together as physically possible, and the nuclear forces keep everything in order but also separated. However, considering the intuitive approach I'm striving for, "matter together" is simpler to understand.}.
However, those subatomic particles are more commonly referred to as {url elementary particles;https://en.wikipedia.org/wiki/Elementary_particle}, 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 {b 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 {b a lot}.
So much so, if you get too close, you will not be able to overwhelm gravity and "merge" into that infinitesimal point:

{img https://upload.wikimedia.org/wikipedia/commons/9/91/CNX_UPhysics_13_07_spacetime.png?utm_source=commons.wikimedia.org&utm_campaign=index&utm_content=original}{br}
{i The white lines are our geodesics}{br}
{i Image by OpenStax University Physics on {url Wikipedia;https://commons.wikimedia.org/w/index.php?curid=64307901}}

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 {url gravitational lensing;https://en.wikipedia.org/wiki/Gravitational_lens}.
Though not unique to black holes, that's when the effect is strongest:

{img https://upload.wikimedia.org/wikipedia/commons/0/02/Gravitational_lens-full.jpg}{br}
{i 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 {b 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 {url Schwarzschild metric;https://en.wikipedia.org/wiki/Schwarzschild_metric}, in Schwarzschild coordinates, looks like this:

{img ./img/Schwarzschild-metric-in-Schwarzschild-coordinates-wikipedia.png}

{b YOUCH!}

The variables themselves aren't that bad: {` c} is the speed of light, {` &tau\;} is time, {` &Omega\;}, {` 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 {b too slow} to render in real time.

Thankfully {url Riccardo Antonelli;https://github.com/rantonels} 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:

{img ./img/riccardo-antonelli-schwarzschild-black-hole-force-field-formula.png}

Explaining his steps and the formula is way beyond the scope of my understanding and this post.
He explains it in more mathematical detail {url here;http://spiro.fisica.unipd.it/~antonell/schwarzschild/}.

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

.bulleted
- 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 {url Ricardo Antonelli's;https://github.com/rantonels/starless/blob/7a0c3c4efdc0666f53dcfa48d310fb2164be2d1a/tracer.py} and {url Dmitry Brant's;https://dmitrybrant.com/2018/12/11/ray-tracing-black-holes} 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:

{img ./img/riccardo-antonelli-schwarzschild-black-hole-force-field-formula.png}

{` F{^ &rarr\;}} returns a force vector ({` F} for "force" and {` {^ &rarr\;}} 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:

:unordered
- {` r} is the distance from the center of the black hole (the singularity) to the {b current} position of the corpuscle.

- {` h} is the magnitude of the {url "specific angular momentum";https://en.wikipedia.org/wiki/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 {i 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{fn Of course, I'm assuming the trains are strong enough they won't break upon the collision.}.
  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) {i specific} angular momentum means our object is small and moves fast, while for a short (small magnitude) {i specific} angular momentum out object is large and moves slowly.

  Another simplification we need to is {b make our black hole stationary}!
  Then, rather than finding the angular momentum{fn Angular momentum is actually {url defined using the cross product;https://en.wikipedia.org/wiki/Specific_angular_momentum#Definition}} and dividing by the mass{fn As mentioned, light has no mass, so this version of the formula wouldn't work.}, we'll use another formula:

  {img ./img/specific-angular-momentum-wikipedia.png}

  where {b 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 {b 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.

  You might think that makes no sense, since we cannot divide by zero, so the specific angular momentum shouldn't exist, but then with the cross product we find a value for it.
  You would be correct: Antonelli does some funky approximations, combining Newtonian physics with General theory of relativity, to make the formula that simple.

- {` r{^ ^}} is the normalized ({url unit;https://en.wikipedia.org/wiki/Unit_vector}, hence the {` ^} on top) vector from the center of the black hole to the light particle:

  {img ./img/normalize-vector-wikipedia.png}


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

{img Refer to "Antonelli’s formula - explained" in the Appendix;./img/variables-expanded.png}

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

{img Refer to "Antonelli’s formula - simplified" in the Appendix;./img/antonelli-formula-expanded.png}

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{fn The {url no-hair theorem;https://en.wikipedia.org/wiki/No-hair_theorem} states that stationary black holes are characterized by only their mass, angular momentum and electric charge. By definition of a Schwarzschild black hole, the last two are null.}.
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 {url ISCO;https://en.wikipedia.org/wiki/Innermost_stable_circular_orbit}.

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:

{img r_s = (2GM)/(c^2);./img/Schwarzschild-radius-wikipedia.png}

where {` G} is the {url gravitational constant;https://en.wikipedia.org/wiki/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:

{img r = r_s * 3/2;./img/photon-sphere-wikipedia.png}

where {` G} and {` M} are the same as before, {` c} is the speed of light in vacuum and {` r{^ s}} 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{fn It's not really true that gravity stopps affecting spacetime after the ISCO, even though things can orbit Earth that doesn't mean that Earth doesn't bend spacetime. However, for the purposes of a simple shader, this assumption produces perfectly useable results.} and objects aren't affected.
Formula:

{img ./img/isco-wikipedia.png}

Putting them all together, we should expect light to travel like this:
{img https://upload.wikimedia.org/wikipedia/commons/2/27/Black_Hole_Shadow.gif}{br}
{i Animation by Hugo Spinelli on {url Wikipedia;https://commons.wikimedia.org/w/index.php?curid=138168679}}

## Light-bending algorithm

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

.bulleted
- 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):

{img Refer to "Algorithm" section of the Appendix;./img/algorithm.png}

as a reminder, here is the definition of {` F{^ &rarr\;}}:

{img ./img/antonelli-formula-expanded.png}

{br}

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 {url latexeditor.lagrida.com;https://latexeditor.lagrida.com/}, by using Firefox's built-in ability to screenshot HTML nodes, in particular the {url MathJax;https://www.mathjax.org/} node with the equation.

#### Antonelli's formula - explained

```tex
'snippet 1 $ ./antonelli_formula.tex

#### Antonelli's formula - simplified

```tex
'snippet 1 $ ./antonelli_formula_simplified.tex

#### Algorithm

```tex
'snippet 1 $ ./algorithm.tex

{footnotes}