- #1
LightAintSoFast
- 3
- 0
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 spatial 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 spatial 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:
##
\newcommand{\partiald}[2]{\frac{d #1}{d #2}}
\newcommand{\partialdt}[2]{\Big( \partiald{#1}{#2} \Big)^2}
\newcommand{\partialdd}[2]{\frac{d^2 #1}{d {#2}^2}}
\newcommand{\vp}{\varphi}
\newcommand{\rrs}{\Big(1 - \frac{r_s}{r} \Big)}
\newcommand{\rdot}{\dot{r}}
\newcommand{\rddot}{\ddot{r}}
\newcommand{\phiddot}{\ddot{\vp}}
\newcommand{\phidot}{\dot{\vp}}
\newcommand{\ldot}{\dot{\ell}}
\newcommand{\lddot}{\ddot{\ell}}##
$$ \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 spatial 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 into 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):
$$\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 spatial 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 spatial 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:
##
\newcommand{\partiald}[2]{\frac{d #1}{d #2}}
\newcommand{\partialdt}[2]{\Big( \partiald{#1}{#2} \Big)^2}
\newcommand{\partialdd}[2]{\frac{d^2 #1}{d {#2}^2}}
\newcommand{\vp}{\varphi}
\newcommand{\rrs}{\Big(1 - \frac{r_s}{r} \Big)}
\newcommand{\rdot}{\dot{r}}
\newcommand{\rddot}{\ddot{r}}
\newcommand{\phiddot}{\ddot{\vp}}
\newcommand{\phidot}{\dot{\vp}}
\newcommand{\ldot}{\dot{\ell}}
\newcommand{\lddot}{\ddot{\ell}}##
$$ \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 spatial 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 into 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: