# A Equations for computing null geodesics in Schwarzschild spacetime

#### LightAintSoFast

My project for obtaining my master's degree in computer science involved ray tracing in Schwarzschild spacetime in order to render images of black holes. These light rays had to be computed numerically using the geodesic equation. However, I ran into a problem. The geodesic equation is given as:
$$\frac{d^2 x^\mu}{d \lambda^2} + \Gamma^{\mu}_{\alpha_\beta} \frac{dx^\alpha}{d \lambda} \frac{dx^\beta}{d \lambda} = 0$$
We can find the Christoffel symbols and thus obtain an expression for $\frac{d^2 r}{d \lambda^2}$ and $\frac{d^2 \varphi}{d \lambda^2}$. We can then numerically integrate these two second order ODEs in order to find $r(\lambda)$ and $\varphi(\lambda)$. This suffices to trace the light ray, because we assume that $\theta = \pi/2$ and we are only interested in the shape of the curve that the light ray produces, so $t(\lambda)$ is unimportant. The issue with this is that the parameter $\lambda$ carries no physical meaning as far as I'm aware. I don't want to go into the details of my ray tracing project too much, but for the purpose of my code it was important that I could assume that the light ray traveled at a constant velocity. In other words, when performing a numerical integration on these equations I needed a parameter step of $\Delta$ to always correspond to a coordinate distance step of $\Delta$ as well. This is equivalent to demanding that the parameterization of the curve be chosen such that the curve always has a constant spacial velocity of 1. I set out to produce a new parameterization, $r(\ell)$ and $\varphi(\ell)$ that would satisfy this condition, which is given as:
$$\Big( \frac{\partial r}{d \ell} \Big)^2 + r^2 \Big( \frac{\partial r}{d \ell} \Big)^2 = 1$$
$\frac{\partial r}{d \ell}$ is the radial velocity of the curve while $r \frac{\partial \varphi}{d \ell}$ is the tangential velocity. The squared magnitude of the total velocity is simply the sum of their squares, which we are requiring to always be 1. Now we come to the parameter $\ell$, which is distance as measured by a stationary observer at infinity. This is very intuitive if you think about it for a little bit: the functions $r(\ell)$ and $\varphi(\ell)$ tell you the position of the light ray after it has traveled a distance of $\ell$ through space (as measured by our stationary observer at infinity). Since we only care about the shape of the curve, we can suppress the $t$ coordinate of Schwarzschild spacetime and simply think about the light ray as a curve through 3 dimensional space. Thus parameterizing the curve according to spacial distance makes sense, and results in the curve having a constant velocity of 1. The question is: how do we numerically compute $r(\ell)$ and $\varphi(\ell)$? To start off we can write the geodesic equation in terms of coordinate time $t$. (This is all on wikipedia: https://en.wikipedia.org/wiki/Geodesics_in_general_relativity)

$$\frac{d^2 x^\mu}{d t^2} = - \Gamma^{\mu}_{\alpha_\beta} \frac{dx^\alpha}{d t} \frac{dx^\beta}{d t} + \Gamma^{0}_{\alpha_\beta} \frac{dx^\alpha}{d t} \frac{dx^\beta}{d t} \frac{dx^\mu}{d t}$$
Computing the Christoffel symbols is a bit of a chore but we end up with the two equations:

$$\rddot = -\frac{c^2 r_s}{2r^2} \rrs + \frac{3r_s}{2r^2} \rrs^{-1} \rdot^2 + (r - r_s) \phidot^2$$
$$\phiddot = -\frac{2}{r} \rdot \phidot + \frac{r_s}{r^2} \rrs^{-1} \rdot \phidot$$
Where the dot of course refers to derivative with respect to $t$. Now, we can find $\partialdd{r}{\ell}$ by using the chain rule.
$$\partialdd{r}{\ell} = \partialdd{r}{t} \partialdt{t}{\ell} + \partiald{r}{t} \partialdd{t}{\ell} = \rddot \partialdt{t}{\ell} + \rdot \partialdd{t}{\ell}$$
We simply need to know $\partiald{t}{\ell}$ and $\partialdd{t}{\ell}$. Notice that
$$\partiald{t}{\ell} = \frac{1}{\partiald{\ell}{t}} = \frac{1}{\ldot}$$
Let's think about what $\ldot$ is. It is observed change in spacial position of the light ray versus observed changed in time. This is the delayed velocity and can be given in terms of $r, \rdot,$ and $\phidot$ as $\ldot = \sqrt{\rdot^2 + r^2 \phidot^2}$. We can also use the second derivative rule for inverse functions to find
$$\partialdd{t}{\ell} = - \frac{\lddot}{\ldot^3}$$
The value $\lddot$ can be obtained by deriving $\ldot$ with respect to time, resulting in an expression containing $r, \rdot, \phidot, \rddot$ and $\phiddot$. But as we have seen, $\rddot$ and $\phiddot$ are functions of $r, \rdot$ and $\phidot$. Thus the whole expression for $\partialdd{r}{\ell}$ can be given entirely in terms of $r, \rdot$ and $\phidot$. Plugging everything in and simplifying is messy, but in the end (using the fact that $\partiald{r}{\ell} = \frac{\rdot}{\ldot}$ and $\partiald{\vp}{\ell} = \frac{\phidot}{\ldot}$) we can obtain the result we need:
$$\partialdd{r}{\ell} = -\frac{3r_s}{2} r^2 \Big( \partiald{\vp}{\ell} \Big)^4 + r \partialdt{\vp}{\ell}$$
Where $r_s$ is the Schwarzschild radius of the gravitating mass. There are two terms in this equation: the first is an attractive term that is proportional to both the current radius and the angular velocity. The second term is repulsive and is just the classic centripital acceleration. As it turns out, we don't even need to find $\partialdd{\vp}{\ell}$. Using the fact that the velocity of the curve is constant we trivially get an equation for $\partiald{\vp}{\ell}$:
$$\partiald{\vp}{\ell} = \frac{1}{r} \sqrt{1 - \partialdt{r}{\ell}}$$
These two coupled differential equations are enough to numerically compute $r(\ell)$ and $\vp(\ell)$. I used mathematica to try some concrete examples (setting $r_s = 1$ for simplicity). Let's start with the initial conditions for a light ray in orbit at the photon sphere. We can always set $\vp(0) = 0$ because it doesn't matter where our $\vp$ coordinate starts. We also set $r(0) = 1.5$ and $r'(0) = 0$.

The light ray just goes in a circle, as expected. Orbits about the photon sphere are unstable, so slight perturbations will result in the light ray either being captured by the black hole or going off to infinity. Let's set $r'(0) = 0.01$. In this case the light ray should escape to infinity because it is at the photon sphere with a positive radial velocity.

The equations seem to work properly. Now, one of the reasons I am making this post is because, interestingly enough, these equations that I've derived can apparently be used to trace light rays past the event horizon. Both of the differential equations are perectly well defined at $r = r_s$. If I set $r'(0) = -0.01$ we would expect the light ray to be captured by the black hole. Here is the result of that computation:

The light ray spirals in to the origin and hits the singularity. Mathematica gives an error if you try to trace the light ray past the origin, which is to be expected.

I think this is super interesting! I have tried numerically computing the light rays using both $\lambda$ and $t$ as parameters and, as I expected, light rays will not cross the event horizon at a finite parameter value. Here is a plot of $r(\lambda)$ for a light ray emitted at $r(0) = 5$ traveling directly towards the origin:

As you can see, it approaches the event horizon but never crosses it. The plot for $r(t)$ looks similar. But here is the plot for $r(\ell)$:

It just travels straight past the event horizon and hits the origin, no problem. I think this is cool because I have been reading that different coordinate systems such as Kruskal-Szekeres coordinates are used to make it so that light rays will fall past the event horizon in finite parameter distance. But my approach just uses the traditional Schwarzschild coordinates and it still works. I wanted to share my finding with other people, mainly because I'm not a physicist and am curious as to whether this kind of thing has already been done. I apologize in advance if I made mistakes in my post or if I've not been clear enough in my writing. And if anyone is interested, here is an example of an image I rendered using this technique to trace light rays (I didn't embed it because it's very large):

Last edited:
Related Special and General Relativity News on Phys.org

#### PeterDonis

Mentor
while $r \frac{\partial r}{d \ell}$ is the tangential velocit
I assume you mean $r \frac{\partial \varphi}{\partial \ell}$, correct?

#### PeterDonis

Mentor
these equations that I've derived can apparently be used to trace light rays past the event horizon
This can't work the way you are assuming, because at and inside the horizon, the spacetime is no longer static, so there is no such thing as "space" independent of time. So your technique of eliminating the time coordinate and just looking at the "path in space" of the light ray can't be correct.

#### LightAintSoFast

I assume you mean $r \frac{\partial \varphi}{\partial \ell}$, correct?
Yes, I just fixed it.
This can't work the way you are assuming, because at and inside the horizon, the spacetime is no longer static, so there is no such thing as "space" independent of time. So your technique of eliminating the time coordinate and just looking at the "path in space" of the light ray can't be correct.
Hmm that could be the case. Anyways every time I've tried tracing the light ray past the event horizon it always spirals in to the origin. But maybe the actual path it's computing is not correct? The Schwarzschild metric in terms of Schwarzschild coordinates is still defined below the event horizon, so in principle measures of spacial distance should still carry meaning shouldn't they?

#### PeterDonis

Mentor
maybe the actual path it's computing is not correct?
"Path" as you've defined it doesn't even have a well-defined meaning below the horizon, so whatever you're computing can't be correct as an actual description of something.

The Schwarzschild metric in terms of Schwarzschild coordinates is still defined below the event horizon
Yes, but...

in principle measures of spacial distance should still carry meaning shouldn't they?
Not as you're defining "spatial distance", because the $t$ coordinate, which is the one the metric does not depend on, is not timelike at or below the horizon. That means that the 3-surfaces orthogonal to $t$, which are the ones you are calling "space" (and then you're considering the $\theta = \pi / 2$ "equatorial plane" within them), are not spacelike below the horizon. And if they're not spacelike, you can't interpret them as "space".

#### PeterDonis

Mentor
Kruskal-Szekeres coordinates
You might want to read this Insights article on Schwarzschild spacetime (and the next one, part 3 in the series), which uses Kruskal coordinates to show how the coordinate singularity at $r = 2m$ in Schwarzschild coordinates actually works:

Note in particular what the surfaces of constant $t$ look like.

#### LightAintSoFast

You might want to read this Insights article on Schwarzschild spacetime (and the next one, part 3 in the series), which uses Kruskal coordinates to show how the coordinate singularity at $r = 2m$ in Schwarzschild coordinates actually works:

Note in particular what the surfaces of constant $t$ look like.
Thanks I'll have a look. I think I understand your point about how spacial distance doesn't work the same below the event horizon. I'm going to think about this more later because it's late where I am.

#### m4r35n357

Thanks I'll have a look. I think I understand your point about how spacial distance doesn't work the same below the event horizon. I'm going to think about this more later because it's late where I am.
I've spent a lot of time writing various simulations like this, and my advice is to always find a way of checking your work. In the early days I used GRorbits. It is often a bit fiddly (to say the least) to arrange comparing like with like but it is the only way to be sure, as they say ;)

[EDIT] This is a reference for Kerr spacetime, but you might find it useful if you don't already know of it.

Last edited:

#### Ibix

You can actually solve the equation for radial infall and get a closed form solution for $\lambda(r)$. The solution is valid above the horizon and below, because one can define a set of coordinates with the same functional form as Schwarzschild's exterior coordinates inside the horizon. However the $r$ coordinate is, as Peter says, timelike below the horizon. So your plot switches from being something like "angular position as a function of distance from the horizon" to "angular position as a function of time to impact on the singularity" ($r$ isn't directly either time or distance, but you get the idea).

"Equations for computing null geodesics in Schwarzschild spacetime"

### Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving