Sunday, May 24, 2015

Geodesics in the schwarzschild metric

I'll resume making posts about biking once I get in some good rides this week. I'm planning a cool gravel ride up to the Amana colonies tomorrow. It'll be my first time using my new shoes! We've been busy here at the ranch partying. My dad Tom just got his Ph.D. so we threw a big party to celebrate. Go dad! We've been drinking and eating and generally having a good time. I even went out bowling the other night. Hanging out with James and Iju was cool, and I'm looking forward to seeing them a little more before I start out on my trip in June.

I've had a little time to play around with the geodesic equations associated with the Schwarzschild metric. I apologize in advance if I misspell Schwarzschild. If you're interested in reading a derivation of the metric and a proof that it is a solution to the Einstein field equations in vacuum, you'll have to look elsewhere.... Anyways, the line element in this metric looks like the following,
$$ds^2 = (1-\frac{r_g}{r})dt^2 - (1-\frac{r_g}{r})^{-1}dr^2 - r^2(d\theta^2 + sin^2(\theta)d\phi^2) $$
where $r_g$ is the Schwarzschild radius. If you're like me and you stop understanding things once $c = G = 1$ then
$$ r_g = \frac{2GM}{c^2} $$
Remember that the line element is related to the metric by the following relationship.
$$ ds^2 = g_{\mu\nu}dx^{\mu}dx^{\nu} $$
We could visualize the metric as a four by four matrix that we contract twice with our differential vectors.
$$ g_{\mu\nu} = \left[ \begin{array}{c} 1-\frac{r_g}{r} & 0 & 0 & 0 \\ 0 & - (1-\frac{r_g}{r})^{-1} & 0 & 0 \\ 0 & 0 & -r^2 & 0 \\ 0 & 0 & 0 & r^2sin^2(\theta)  \end{array} \right] $$
Normally the first thing we do when solving the geodesic equations for this metric is say that $\theta = \frac{\pi}{2}$, but because I'm using Python, it doesn't really matter. The job at hand is to solve the following four differential equations:
$$ \frac{d}{d\lambda}\left( \frac{dx^{\mu}}{d\lambda}  \right) + \Gamma^{\mu}_{\alpha\beta} \frac{dx^{\alpha}}{d\lambda} \frac{dx^{\beta}}{d\lambda}= 0 $$
$\lambda$ is just some parameter. The $\Gamma$ symbol is related to the metric through the following relation.
$$ \Gamma^{\gamma}_{\beta\mu} = \frac{1}{2}g^{\alpha\gamma}(g_{\alpha\beta,\mu} + g_{\alpha\mu,\beta} - g_{\beta\mu,\alpha}) $$
If you're not familiar with the upper versus lower index notation, in this situation, the upper indices on $g$ indicate that we're using the inverse metric. Calculating this looks pretty rough, but I made some Python code that calculates these "Christoffel symbols" with the help of SymPy. Check out this snippet below (I think some of it might get cut off -- just check it out on Github if you're really interested). This little chunk of code makes forming the geodesic equation much simpler. Ultimately, after plugging everything in, we get a set of reasonably nasty looking equations.
$$ \begin{array}{rcl} \ddot{t} &=& -\frac{r_g}{r(r-r_g)}\dot{r}\dot{t} \\
\ddot{r} &=& -\frac{1}{2}\frac{r_g(r-r_g)}{r^3}\dot{t}^2 + \frac{1}{2}\frac{r_g(r-r_g)}{r^3(1-\frac{r_g}{r})^2} + (r-r_g)sin^2(\theta)\dot{\theta}^2 - (r-r_g)\dot{\phi}^2 \\
\ddot{\theta} &=& -\frac{2\dot{r}\dot{\theta}}{r} - cot(\theta)\dot{\theta}^2 \\
\ddot{\phi} &=& -\frac{2\dot{r}\dot{\phi}}{r} \end{array}$$
Solving these equations numerically is cumbersome, but not super hard! We just make the same substitutions I made in my previous post and we get a system of coupled first order ODEs. I ended up just solving this using forward Euler's method, because I didn't want to manually form the Jacobian matrix for this system. That said, I do plan on making a code that will automatically form the Jacobian matrix given some system of equations and set of variables. This system would greatly benefit from backward Euler, as the step size I ended up using was really small, and it took forever to calculate the geodesics. Anyways, below is a picture of one of the geodesics I calculated, assuming that $\theta=\frac{\pi}{2}$. The red circle is a circle with Schwarzschild radius equal to one (the same radius I used in calculating the geodesic). Pretty cool!
I am by no means convinced that this solution is correct, so if you have any suggestions, make a comment!


1 comment:

  1. Nice....but check line 5 of the 2md equation - I'm not sure that's a Cot? Shldnt it be tang??

    ReplyDelete