What’s the shortest path between two points - A and B? Most people instinctively say it’s a straight line between those two points. Which is true on a flat 2D plane, but what about on the surface of the earth? What’s the shortest path between two points on the surface of the earth? Let’s ask Google:

That’s clearly not a straight line. So is Google Maps bugged, or is math bugged? Have we uncovered a conspiracy by Big Globe? Actually, turns out it is a straight line! But it’s a straight line on a sphere, and straight lines on a sphere aren’t straight on a 2D projection of the sphere, which is what a map is!

In fact, the most common projection used in maps, the Mercator Projection, infamously distorts distances. This results in countries farther from the equator appearing grossly distorted (made MUCH larger).

So what gives? Can we generalize a way to always find the shortest path between two points, regardless of the surface we are working with? Turns out, yes, we can (more on this later)! This is known as Geodesics. If you have seen any videos about General Relativity on YouTube, you definitely have heard “Geodesic” at least once. In this blog, I will motivate why we need geodesics, and personally why I needed geodesics.

This blog contains math and code, but you can go through and understand the whole thing without actually doing either. You have to just trust me bro.

However, true satisfaction will only come when you actually do at least some of both, and see that it all works out. So, I will try to explain the math as well, but I will not go into too much detail. If you want to see the proofs, many good videos on YouTube explain them in detail. For the code, I will provide snippets, and a lengthier post about my path tracer will be coming soon.

Motivation

I must confess, I am not a mathematician. I am a Software Engineer. Normally, if for some reason I came across the math for geodesics, I would simply close that browser tab and move on. However, I was finally tackling a bucket-list hobby project I wanted to make for a long, long time. Unfortunately, skill and time issues were holding me back. The goal was to recreate the famous Interstellar scene of the black hole.

This is the result. If this is not motivation enough, I really don’t know what is. Anyway, now that we know the payoff, time for the math (and code)!

Before I piss off the physicists and mathematicians

In order to keep this blog somewhat reasonable, I will be glossing over the whys and only focussing on the hows. My motivation, really, is to get that pretty image of the black hole. Maybe in the future, when my ChatGPT-Chromium-VSCode wrapper finally takes off and I make a gazillion dollars and retire, I will learn and write about the mysteries of the universe in more detail. Anyway, enjoy the ride!

The Plan

The plan is to learn:

• What’s a tensor?
• What the **** is a $g_{\mu\nu}$ ?
• What’s a metric?
• The generalized distance formula
• The Bridge: The Euler-Lagrange Equation
• The Geodesic Equation
• Christoffel symbols
• A concrete example

When learning hard things, I find it helps to turn off your ape brain screaming “THIS IS TOO HARD! I DON’T UNDERSTAND! TURN AWAY NOW!”, and just follow the instructions. Often, once you go through the process once, it becomes clear that, in fact, it’s really easy, just bogged down by notation and tedium. If you are a programmer, this will probably be very easy for you. Really, you can understand the whole thing using high school calculus. Yeah you won’t be able to derive it but you can appreciate how it works, and that’s enough to get a pretty picture, so it’s good enough for me (for now).

What’s a tensor?

Let’s consult Wikipedia:

In mathematics, a tensor is an algebraic object that describes a multilinear relationship between sets of algebraic objects associated with a vector space. Tensors may map between different objects such as vectors, scalars, and even other tensors. There are many types of tensors, including scalars and vectors (which are the simplest tensors), dual vectors, multilinear maps between vector spaces, and even some operations such as the dot product. Tensors are defined independent of any basis, although they are often referred to by their components in a basis related to a particular coordinate system; those components form an array, which can be thought of as a high-dimensional matrix.

Yikes! Let me use programmer-translate on that: For our purposes, a tensor is simply a multi-dimensional array. You can define a tensor:

$$g_{\mu\nu} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$

And it’s really just a 2D array like [[1, 0], [0, 1]]. But hold on a second…

What the **** is a $g_{\mu\nu}$ ?

$g_{\mu\nu}$ says “you need two indices $\mu$ and $\nu$ , to get a scalar value from this tensor.” Yes, $\mu$ and $\nu$ are just array indices! Okay, so.. how do we index them? You can simply index them like in your favorite language, but by convention, the indices start at 1 🤢 sorry. So, the components of the tensor can be accessed like

$$g_{11} = 1, g_{12} = 0$$ $$g_{21} = 0, g_{22} = 1$$

So, really, $$g_{\mu\nu} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$ is sort of like saying g[i][j] = [[1, 0], [0, 1]].

Now, so far, $\mu$ and $\nu$ have been generic indices. You can get a value, say $g_{12}$ , but it doesn’t really mean anything. To be useful, we attach a coordinate system with a tensor! Let’s say we’re working on the 2D cartesian plane, so our coordinates are $(x, y)$ (The order is important!) This means, our first index, $\mu$ corresponds to the $x$ coordinate, and the second index, $\nu$ corresponds to the $y$ coordinate. Here’s all the components of the tensor above:

$$g_{xx} = 1, g_{xy} = 0$$ $$g_{yx} = 0, g_{yy} = 1$$

Here, $g_{\mu\nu}$ corresponds to “how much of the $\nu$ component relates to the $\mu$ component”? This probably doesn’t make sense, so I introduce…

What’s a metric?

A metric, (a.k.a. a metric tensor, the “tensor” part is often omitted to make it sound mysterious and hard), is just a tensor that helps us calculate the distance between two points in a given space and coordinate system. I copied this from Wolfram, so it must be correct. So what does it look like? Well, you’ve already seen one! The tensor I’ve been using is a metric tensor for the flat 2D plane in cartesian coordinates $(x, y)$ !

$$g_{\mu\nu} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$

Now, we can answer the question from the previous section, “how much of the $\nu$ component relates to the $\mu$ component”? Let’s take $g_{xx}$ . This tells us, “how much of the $x$ direction relates to the $x$ ” direction? They’re the same, so the value is $1$ . Now, what about $g_{xy}$ ? This tells us “how much of the $y$ direction relates to the $x$ direction?” These two directions are orthogonal (perpendicular), so the value is $0$ . Similarly, $g_{yy}$ is $1$ , and $g_{yx}$ is $0$ . Basically, moving in the $x$ direction doesn’t change your $y$ coordinate, and vice-versa.

Let’s see some other metric tensors, say for spherical metric tensor in spherical coordinates ($r, \theta, \phi$ ):

$$g_{\mu\nu} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & r^2 & 0 \\ 0 & 0 & r^2sin^2({\theta}) \end{pmatrix}$$

And the Schwarzschild metric tensor in spherical coordinates ($t, r, \theta, \phi$ ):

$$g_{\mu\nu} = \begin{pmatrix} -(1 - \frac{2GM}{r}) & 0 & 0 & 0 \\ 0 & (1 - \frac{2GM}{r})^{-1} & 0 & 0 \\ 0 & 0 & r^2 & 0 \\ 0 & 0 & 0 & r^2sin^2({\theta}) \end{pmatrix}$$

As an aside, this is the metric tensor that describes the space around a uncharged non-rotating black hole! This is the one used to create the image you saw above, although the one in Interstellar is a rotating black hole, which is described by the even more complicated Kerr Metric. However, you can’t really tell the difference. We’ll come back to this metric soon.

In any case, the metric tensor directly tells us how to calculate the distance between two points regardless of how complicated the underlying space is, as long as we are willing to put up with the tedious math. To do it we just have to integrate…

The generalized distance formula

Let me just hit you with the thing:

$$ds^2 = \Sigma{g_{\mu\nu}dx^{\mu}dx^{\nu}}$$

Let’s unpack this…

• $g_{\mu\nu}$ is the metric tensor of the space you are working with
• $x^{\mu}$ is a compact way to refer to any of the coordinates in our system. It means whatever coordinate index you are using in place of $\mu$ , replace $x^{\mu}$ with that. So in our 2D system with coordinates $(x, y)$ , $dx^x$ would be $dx$ , and $dx^y$ would be $dy$ . This will make more sense momentarily.
• What the $\Sigma$ ?! Remember, $x^{\mu}$ , $x^{\nu}$ are just placeholders. You have to replace them with all the possible combinations of your coordinate system. The $\Sigma$ says after you do that, you sum them all together.

Fun aside: Physicists are lazy, and got tired of writing $\Sigma$ over and over again. So we just don’t write the $\Sigma$ , and anytime we have a repeated index like $\mu$ or $\nu$ , assume it is summed over all possible values for those indices. This has unfortunately caught on, much to the detriment to everyone learning this for the first time, and now it’s called “Einstein notation”. If you look online, you’ll see the distance formula written as $ds^2 = {g_{\mu\nu}dx^{\mu}dx^{\nu}}$ . They are the same thing!

Let’s take our good old flat 2D plane with cartesian coordinates, $(x, y)$ . This means $x^{\mu}$ and $x^{\nu}$ have to take on the coordinates $x$ and $y$ . And we know the metric tensor for this space and coordinate system is $$g_{\mu\nu} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$

So, the distance formula in this space and coordinate system is: $$ds^2 = (g_{xx}dx^xdx^x) + (g_{xy}dx^xdx^y) + (g_{yx}dx^ydx^x) + (g_{yy}dx^ydx^y)$$

We already know how to index the metric tensor, so we find the values like this: $g_{xx} = 1, g_{xy} = 0, g_{yx} = 0, g_{yy} = 1$ . And for the confusing part, $x^x = x$ and $x^y = y$ ! We just replace the whole thing by the coordinate, it not an exponent!

Plugging everything in to the equation, we get: $$ds^2 = (1*dxdx) + (0*dxdy) + (0*dydx) + (1*dydy)$$ $$\implies ds^2 = dx^2 + dy^2$$

Does this look familiar? It should, because it’s just the Pythagorean theorem! This is saying, if we move some infinitesimal distance $dx$ in the $x$ direction, and $dy$ in the $y$ direction, the square of the total distance is the sum of squares of the individual distances! That wasn’t so hard, was it? Just a lot of symbols!

Well, that was pretty cool! However, we are still only getting the distance, not the path! To do that, we’ll need…

The Bridge: The Euler-Lagrange Equation

We have the equation for smallest distance, but what’s the shortest path corresponding to that distance? This is where the Euler-Lagrange equation comes in. This equation comes from the Calculus of Variations, and is really really cool! I won’t go into the details here, but you can watch this video to get a sense of what’s going on. All you need to know, the Euler-Lagrange equation is a way to find the function that minimizes (or maximizes) some quantity. In our case, we want to minimize the distance between two points. The Euler-Lagrange equation looks like this:

$$\frac{\partial{L}}{\partial{x^{\alpha}}} - \frac{d}{ds}\left(\frac{\partial{L}}{\partial{\dot{x^{\alpha}}}}\right) = 0$$

In our case, we want to minimize the distance, so we set $L = ds$ . Now, $\alpha$ is again a placeholder, so we have to replace it with a coordinate. And, $\dot{x^{\alpha}}$ means $\frac{dx^{\alpha}}{ds}$ , the derivative of the coordinate with respect to the distance. This is just a fancy way of saying “the rate of change of the coordinate as we move along the path”. So, we have:

$$L^2 = ds^2 = {g_{\mu\nu}dx^{\mu}dx^{\nu}}$$ $$\implies L = \sqrt{g_{\mu\nu}dx^{\mu}dx^{\nu}}$$

If you put your faith in me (and the many videos online), you can plug this into the Euler-Lagrange equation, and after a lot of algebra, you will get…

The Geodesic Equation

$$\frac{d^2x^{\alpha}}{ds^2} + \Gamma^{\alpha}_{\mu\nu}\frac{dx^{\mu}}{ds}\frac{dx^{\nu}}{ds} = 0$$

Ignoring all the indices and the $\Gamma^{\alpha}_{\mu\nu}$ , the bottom line is, if we are able to somehow solve this thing, it will give us the equations for the shortest path. However, notice the $\frac{d^2x^{\alpha}}{ds^2}$ . Here, again, $\alpha$ is a placeholder, and we have to replace it with a coordinate. Then, if we solve the equation, we will get another equation for the shortest path in that coordinate direction!.

Awesome! Now, what about the $\Gamma^{\alpha}_{\mu\nu}$ ? Well…

The Christoffel symbols

The Christoffel symbols are just functions of the metric. They are easily calculatable. But there is a slight problem. See if you can guess it.

$$\Gamma^{\alpha}_{\mu\nu} = \frac{1}{2}g^{\alpha\beta}(g_{\nu\beta,\mu} + g_{\mu\beta,\nu} - g_{\nu\mu,\beta})$$

Yes, pretty gnarly. Let’s unwrap this.

• $\alpha$ and $\beta$ are also placeholders. You need to replace them with coordinates.
• The $g^{\mu\nu}$ . The metric tensor is always a 2D matrix. So, you can calculate the inverse of that matrix. I’m going to use the notation $(g_{\mu\nu})^{-1}$ (don’t @ me mathematicians). Note that this will also be a matrix of the same size. In fact, this is also a tensor, the inverse metric tensor! We define $g^{\mu\nu} = ((g_{\mu\nu})^{-1})_{\mu\nu}$ .
• We simply define $g_{\mu\nu,\beta} = \frac{\partial{g_{\mu\nu}}}{\partial{x^{\beta}}}$ .

Now notice, nothing here is super advanced stuff, just some linear algebra and partial derivatives. Pretty doable, we did this in high school! But the problem, as I’m sure you’ve noticed by now, is that the Christoffel symbols have $3$ indices. And, to calculate one symbol, we need $4$ indices. That’s a ridiculous amount of pure computation. Luckily, we don’t have to do it for our original use case, they’ve already been calculated, and we can just plug in the values, but if you are feeling particularly masochistic, do give it a try! If we have $n$ coordinates, there will be $n^3$ Christoffel symbols! Good luck!

Let me give you a Python snippet with a magic calculate_partial_derivative function to illustrate the point.

g = [[1, 0], [0, 1]]  # Some metric tensor
inv_g = [[1, 0], [0, 1]] # The inverse metric
coords = ['x', 'y'] # The coordinate axes

# We'll store the values like christoffel_symbols[f'{alpha}_{mu}_{nu}'] = value
christoffel_symbols = dict() 

for alpha in coords:
for mu in coords:
for nu in coords:
 value = 0.0
 symbol = f'{alpha}_{mu}_{nu}'
 alpha_idx = coords.index(alpha)

for beta in coords:
 beta_idx = coords.index(beta)

 g_nu_beta_mu = calculate_partial_derivative(nu, beta, mu)
 g_mu_beta_nu = calculate_partial_derivative(mu, beta, nu)
 g_mu_nu_beta = calculate_partial_derivative(mu, nu, beta)

 internal_sum = g_nu_beta_mu + g_mu_beta_nu - g_mu_nu_beta
 factor = 0.5 * inv_g[alpha_idx][beta_idx]

 value += factor * internal_sum

 christoffel_symbols[symbol] = value

Yes, this means it’s $O(n^4)$ . (Yes yes I know [].index() is not constant time, I don’t care). It’s extremely tedious, and very hard to do by hand. But, it’s simple. A really focussed 17 year old can do this, given enough coffee. Unfortunately I’m too old and I cannot do this, I just look it up 😀

Phew, now that we got that out of the way, let’s get back to…

A concrete example

Actually, that’s it! No new equations, no new concepts. Just plug in the Christoffel symbols you calculated (or looked up), and solve the differential equation! You will get the equations for the shortest path in your coordinate system!

Let’s take the example of the flat 2D plane in cartesian coordinates again. We already know the metric tensor is $$g_{\mu\nu} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$

Okay, so let’s take a look at the Christoffel symbol equation again.

$$\Gamma^{\alpha}_{\mu\nu} = \frac{1}{2}g^{\alpha\beta}(g_{\nu\beta,\mu} + g_{\mu\beta,\nu} - g_{\nu\mu,\beta})$$

Notice that it’s a product of $g^{\alpha\beta}$ and $(g_{\nu\beta,\mu} + g_{\mu\beta,\nu} - g_{\nu\mu,\beta})$ . This means if any one of these is zero, the full Christoffel symbol is zero!

So let’s think, when is $g^{\alpha\beta} = 0$ ? Since it’s the inverse metric tensor, and the metric tensor is a diagonal matrix, its zero only when $\alpha \neq \beta$ . This is because the inverse of a diagonal matrix has 1/(diagonal element) on the diagonal, zeros elsewhere. That is to say,

$${\begin{pmatrix} a & 0 \\ 0 & b \end{pmatrix}}^{-1} = {\begin{pmatrix} \frac{1}{a} & 0 \\ 0 & \frac{1}{b} \end{pmatrix}}$$

This is true for all sizes of inverse matrices.

This is great! Since our metric tensor is a diagonal matrix, we only need to consider the symbols where $\alpha = \beta$ ! This reduces the number of symbols to calculate significantly. However, let’s take a look at the other factor as well.

$(g_{\nu\beta,\mu} + g_{\mu\beta,\nu} - g_{\nu\mu,\beta})$ , where $g_{\mu\nu,\beta} = \frac{\partial{g_{\mu\nu}}}{\partial{x^{\beta}}}$ . So each term is a partial derivative of some term of the metric tensor. But wait.. our metric tensor is constant! For any $\mu, \nu$ , $g_{\mu\nu}$ is either $0$ or $1$ . It doesn’t matter what you differentiate a constant by, the result is always $0$ ! This means, all our Christoffel symbols are $0$ .

So, then, what does our geodesic equation look like?

$$\frac{d^2x^{\alpha}}{ds^2} + \underbrace{\Gamma^{\alpha}_{\mu\nu}\frac{dx^{\mu}}{ds}\frac{dx^{\nu}}{ds}}_{0} = 0$$ $$\implies \frac{d^2x^{\alpha}}{ds^2} = 0$$

Wow! A massive simplification! All that’s left is to plug our coordinates in, and solve the equation.

$$\frac{d^2x}{ds^2} = 0$$ $$\frac{d^2y}{ds^2} = 0$$

Solving these (use Wolfram Alpha if you want), we get:

$$\begin{equation}x = as + b\end{equation}$$ $$\begin{equation}y = cs + d\end{equation}$$

Where $a, b, c$ and $d$ are constants that depend on our initial conditions (where we start and which direction we’re going). Turns out, this is just the equation of a straight line, but in the parametric form, with $s$ as the parameter. You can get rid of the parameter, by re-arranging equation $(1)$ :

$$s = \frac{x - b}{a}$$

and plugging it into equation $(2)$ :

$$y = c\left(\frac{x - b}{a}\right) + d$$ $$\begin{equation}\implies y = \left(\frac{c}{a}\right)x + \left(d - \frac{bc}{a}\right)\end{equation}$$

Now we see it! This is the classic, $y = mx + c$ equation of a straight line!

You can, if you want, plug in $(x_1, y_1)$ and $(x_2, y_2)$ into equations $(1)$ and $(2)$ to find the constants, and then plug them into equation $(3)$ , and you will get

$$y - y_1 = \left(\frac{y_2 - y_1}{x_2 - x_1}\right)(x - x_1)$$

the two point form of a line going from $(x_1, y_1)$ to $(x_2, y_2)$ .

So, after all that heavy machinery, we have finally proved that the shortest path between two points in a flat 2D plane is a straight line! However, now we have a method that works for any space and coordinate system! Just find the metric tensor, calculate the Christoffel symbols, and solve the geodesic equation! Easier said than done, but it’s possible.

Back to reality

We’ve gone through a lot of math, but now it’s time to tackle the real world! Let’s take a look at the geodesic equation again…

$$\frac{d^2x^{\alpha}}{ds^2} + \Gamma^{\alpha}_{\mu\nu}\frac{dx^{\mu}}{ds}\frac{dx^{\nu}}{ds} = 0$$

Often, it’s either impractical or straight up impossible to calculate the general solution to this equation. I mean, were you able to solve the equations above for the sphere? 99.9% people can’t. But most importantly, a computer also can’t! So we must rely on numerical methods, like the Euler method or the more standard, Runge-Kutta methods, specifically the fourth order Runge-Kutta method, known as RK4.

The problem is, RK4 is made specifically for first order differential equations, and the geodesic equation is a second order differential equation. However, we can convert it to a system of first order differential equations by introducing a new variable, $v^{\alpha} = \frac{dx^{\alpha}}{ds}$ . This is just the velocity in the $\alpha$ direction. So, we have a system of two first order differential equations:

$$\begin{equation} \frac{dx^{\alpha}}{ds} = v^{\alpha}\end{equation}$$ $$\frac{dv^{\alpha}}{ds} + \Gamma^{\alpha}_{\mu\nu}v^{\mu}v^{\nu} = 0$$ $$\begin{equation}\implies \frac{dv^{\alpha}}{ds} = -\Gamma^{\alpha}_{\mu\nu}v^{\mu}v^{\nu}\end{equation}$$

So, we can just track the velocities of our path in each direction, and update them as we go along the path. Then, we select some small time step, $\Delta t$ , and update our position and velocity like this for each coordinate per timestep. We have 2 simultaneous equations to solve for each coordinate, and RK4 can do that for us!

Aside: The Euler method

Before we do that, let me just quickly explain the Euler method, since it’s simpler to understand. The Euler method is a way to solve first order differential equations of the form $\frac{dy}{dt} = f(t, y)$ . The idea is simple, if you know the value of $y$ at some point $t$ , and you know the derivative at that point, you can estimate the value of $y$ at some small distance $\Delta t$ away from $t$ by assuming the derivative is constant over that small distance. So, you would have:

$$\begin{equation}y(t + \Delta t) \approx y(t) + \Delta t f(t, y(t))\end{equation}$$

This is a simple and intuitive method, but it’s not very accurate, especially for larger values of $h$ . The error can accumulate quickly, leading to inaccurate results. That’s why we often use more advanced methods like RK4, which provide better accuracy by taking into account the curvature of the solution over each timestep. Which leads us to…

The RK4 method

Runge-Kutta methods help us numerically solve ordinary differential equations (ODEs) of the form $\frac{dy}{dt} = f(t, y)$ . Here, $y$ is the function we want to find, and $f(t, y)$ is a given function that describes how $y$ changes with respect to $t$ .

Runge-Kutta is actually a family of methods, but the most commonly used one is the fourth order method, RK4. The idea is to take multiple “samples” of the derivative at different points within the timestep, and then combine these samples to get a better estimate of the value of $y$ at the end of the timestep.

The real math nerds have already figured out exactly what kind of steps to take and the weights to give them to get the best accuracy. With all that, we first calucate 4 samples $k_1, k_2, k_3, k_4$ like so:

$$\begin{equation}k_1 = f\left(t_n, y_n\right)\end{equation}$$ $$\begin{equation}k_2 = f\left(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} \cdot k_1\right)\end{equation}$$ $$\begin{equation}k_3 = f\left(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} \cdot k_2\right)\end{equation}$$ $$\begin{equation}k_4 = f\left(t_n + \Delta t, y_n + \Delta t \cdot k_3\right)\end{equation}$$

And then, we combine them to get the new value of $y$ at $t + \Delta t$ :

$$\begin{equation}t_{n+1} = t_n + \Delta t\end{equation}$$ $$\begin{equation}y_{n+1} = y_n + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4)\end{equation}$$

And now, the point $(t_{n+1}, y_{n+1})$ is our new point, and we can repeat the process for the next timestep. Note that to calculate a step $n+1$ , we only need the values from step $n$ . This means we can do this iteratively, without needing to store all the previous steps. Since we are doing it for a second order differential equation, we have to do this for both equations, for each coordinate. We call these changing values our “state”. So, if we have $n$ coordinates, our state will be a vector $$\text{state}_i = [x_i^1, x_i^2, ..., x_i^n, v_i^1, v^2, ..., v_i^n]$$ And at the next timestep, $i+1$ , we will have $$\text{state}_{i+1} = [x_{i+1}^1, x_{i+1}^2, ..., x_{i+1}^n, v_{i+1}^1, v^2, ..., v_{i+1}^n]$$ This continues until we reach our destination, or we decide we have had enough.

This method is incredibly accurate, especially when compared to Euler’s method, and is widely used in scientific computing. Stare at this figure and see for yourself how overpowered RK4 is compared to Euler’s method:

If you want more math, check out this incredible video on YouTube:

A little demo

Let’s now look at a demo, where we calculate the shortest path between two points on a sphere with radius $R$ . The metric tensor for a sphere in spherical coordinates $(\theta, \phi)$ is:

$$g_{\mu\nu} = \begin{pmatrix} R^2 & 0 \\ 0 & R^2sin^2({\theta}) \end{pmatrix}$$

And the non-zero Christoffel symbols are:

$$\Gamma^\theta_{\phi\phi} = -sin(\theta)cos(\theta)$$ $$\Gamma^\phi_{\theta\phi} = \Gamma^\phi_{\phi\theta} = cot(\theta)$$

Putting these into the geodesic equation, we get:

$$\frac{d^2\theta}{ds^2} - \sin(\theta)\cos(\theta)\left(\frac{d\phi}{ds}\right)^2 = 0$$ $$\frac{d^2\phi}{ds^2} + 2\cot(\theta)\frac{d\theta}{ds}\frac{d\phi}{ds} = 0$$

Now, if somehow you manage to solve this, you get the equation of a Great Circle. A great circle is just the biggest circle you can draw on the surface of a sphere. Spoiler alert: It’s the circle whose center is at the sphere’s center. Its radius is equal to the sphere’s radius. It splits the sphere into two equal parts. You get it!

So the geodesic equation tells us, the shortest path from A to B on the surface of a sphere is along the great circle containing the two points. Let’s check out a interactive demo demonstrating this fact:

Pretty cool, right? Now that you’ve seen how to do this, you KNOW how to do it for any space and coordinate system! Just find the metric tensor, calculate the Christoffel symbols, and solve the geodesic equation, OR, use RK4 to get the path, and plug in the initial conditions! But wait a second…

The punchline

How do we find the initial conditions? We need to know the initial position, and the initial velocity (direction). The initial position is easy, it’s just the coordinates of point A. But what about the initial velocity? We want to go from A to B, but how do we know which direction to start going in? Well okay, we could find the initial conditions for the flat plane and sphere, but what about a black hole? Or some weird curved space? How do we find the initial velocity then?

Turns out, we can’t! This is a problem known as the Two-Point Boundary Value Problem. It’s really really hard to answer, “If our path starts at A and ends at point B, what should our initial velocity be?"" Because for example, it might be a very long complex path, that may me impractical to compute. There may be infinite shortest paths between two points, for example, if two points are on the opposite ends of a sphere, any great circle passing through both points is a shortest path. So how do we find the right initial velocity? There even may be points that are unreachable from point A, no matter what initial velocity we choose! For example, if point A is inside the event horizon of a black hole, and point B is outside, there is no initial velocity that will get us from A to B, because nothing can escape a black hole.

Well that sucks. But there is a way out. We can use a method called the Shooting Method. The idea is simple, we just guess an initial velocity, and see where we end up. If we end up at point B, great! If not, we adjust our initial velocity based on how far off we were, and try again. We keep doing this until we get close enough to point B. This is a numerical method, and it may not always converge to the right solution, but it’s often good enough for practical purposes.

So, in a cruel game by mother nature, in general, we cannot find the shortest path between two points analytically. We have to use numerical methods, and even then, we have to guess the initial conditions. Anyway, this isn’t enough to stop us from recreating the black hole from Interstellar! This is because we care about the path any light ray takes, not necessarily the shortest path between two points. And we can trace a path just fine! So, without further ado…

The Payoff

Now equipped with all this knowledge, we can finally recreate the black hole from Interstellar! Well sort of. Instead of a rotating black hole, we’ll just do a non-rotating black hole, which is described by the Schwarzschild metric. It’s MUCH simpler than the Kerr metric, and you can’t really tell the difference in the final image anyway. Also, to make an image, we don’t find the shortest path between two points (doesn’t make sense, what would those points even be?), we find the path of light rays from the camera in a random direction, and see where they end up. If they end up at the accretion disk, we color that pixel accordingly, if they end up at the black hole, we color it black. If they end up somewhere else, we color it the background color. This is called ray tracing, and it’s a common technique in computer graphics.

So one final time, let’s go through the motions.

1. The metric tensor for a Schwarzschild black hole in spherical coordinates $(t, r, \theta, \phi)$ is:

$$g_{\mu\nu} = \begin{pmatrix} -(1 - \frac{2GM}{r}) & 0 & 0 & 0 \\ 0 & (1 - \frac{2GM}{r})^{-1} & 0 & 0 \\ 0 & 0 & r^2 & 0 \\ 0 & 0 & 0 & r^2sin^2({\theta}) \end{pmatrix}$$

2. The non-zero Christoffel symbols are:

• Time components:

$$\Gamma^t_{tr} = \Gamma^t_{rt} = \frac{GM}{r^2(1 - \frac{2GM}{r})}$$

• Radial components:

$$\Gamma^r_{tt} = \frac{GM(1 - \frac{2GM}{r})}{r^2}$$ $$\Gamma^r_{rr} = -\frac{GM}{r^2(1 - \frac{2GM}{r})}$$ $$\Gamma^r_{\theta\theta} = -(1 - \frac{2GM}{r})r$$ $$\Gamma^r_{\phi\phi} = -(1 - \frac{2GM}{r})r\sin^2(\theta)$$

• Theta components:

$$\Gamma^\theta_{r\theta} = \Gamma^\theta_{\theta r} = \frac{1}{r}$$ $$\Gamma^\theta_{\phi\phi} = -\sin(\theta)\cos(\theta)$$

• Phi components:

$$\Gamma^\phi_{r\phi} = \Gamma^\phi_{\phi r} = \frac{1}{r}$$ $$\Gamma^\phi_{\theta\phi} = \Gamma^\phi_{\phi\theta} = \cot(\theta)$$

3. Plugging these into the geodesic equation, we get 4 coupled second order differential equations, one for each coordinate. We convert them to 8 coupled first order differential equations using the velocity substitution we did earlier. So we have 8 equations of the form $\frac{dx^{\alpha}}{ds} = v^{\alpha}$ and $\frac{dv^{\alpha}}{ds} = -\Gamma^{\alpha}_{\mu\nu}v^{\mu}v^{\nu}$ , where $\alpha, \mu, \nu$ take on the values $t, r, \theta, \phi$ . This is our system of equations to solve. Since we have 4 coordinates, our state will be a vector of size 8, $$\text{state}_i = [t_i, r_i, \theta_i, \phi_i, v_i^t, v_i^r, v_i^\theta, v_i^\phi]$$

4. We will now solve the equations analytically. Just kidding, not even Math Jesus can do that. So we use RK4 to numerically solve the equations. We start at point A, with some guessed initial velocity, and see where we end up. We keep track of where each timestep leads us, and we draw a continuous line from the camera to the end point. If we end up at the accretion disk, we color that pixel accordingly, if we end up at the black hole, we color it black. If we end up somewhere else, we color it the background color. We repeat this for every pixel in the image, and we get a beautiful image of a black hole!

Here’s a demo of a “path tracer” (literally traces the path of light rays), that shoots 100 rays in a 10x10 grid of “pixels” from an origin (the “Camera”). You can see the curved paths of the rays as they bend around the black hole. You can also see how the light “blue shifts” as it gets closer to the black hole, and “red shifts” as it moves away from the black hole. This is because of the gravitational time dilation caused by the black hole. The closer you are to the black hole, the slower time moves. So, light that is closer to the black hole has a higher frequency (bluer), and light that is farther away has a lower frequency (redder). Pretty cool, right?

In a real path tracer, we shoot many rays per pixel, and collect the final color of the rays by calculating where they end up. If they collide with something, we sample the surface’s color, if they end up at the black hole, we color it black, if they end up somewhere else, we color it the background color. We then average the colors of all the rays to get the final color of the pixel. Over time, shooting millions of rays, thousands of samples per pixel, we get a beautiful, high quality image of a black hole! It looks something like this:

But..but.. where’s the glowing disk? The glowing disk is not a part of the black hole, it’s matter swirling around the black hole, called the accretion disk. To model that accurately, we will need a few more things, which we are definitely not going to do at the moment. However, I have cheated a bit and added a fake accretion disk to the scene so that it looks cool. Adding a real accretion disk would require volumetrics at the least, and that’s a lot of work. I plan on adding this in the future, but for now, it’s a flat, procedurally textured flat disk around the black hole. So, with everything in place, we get this:

You can clearly see the gravitational lensing effect of the black hole, bending the light from the distant stars. And of course, the accretion disk is actually flat, but it is warped by the black hole’s gravity, making it look like a warped ring around the black hole, as made famous by Interstellar. Also the photon sphere, the thin ring of light around the black hole, is visible clearly. This is the path of light that is just barely able to escape the black hole’s gravity, and is bent around the black hole multiple times before escaping. I think that’s a great place to stop.

Conclusion

Phew! That was a long journey, but I hope you enjoyed it. We started with the simple problem of finding the shortest path between two points in a flat 2D plane, and ended up recreating the black hole from Interstellar! Along the way, we learned about the metric tensor, the Christoffel symbols, the geodesic equation, and numerical methods like RK4. We also saw how these concepts are used in general relativity to describe the curvature of spacetime around massive objects like black holes.

I hope your attention span hasn’t been fried by TikTok and reels and shorts and you actually managed to reach the end of this post. If you did, congratulations! Hope you are now eager to implement this yourself! Please do reach out to me if you have any questions, or if you want to share your implementation! I would love to see it.