How to Integrate the geodesic equations numerically?

[Matter_Matters]

Hello there, I've been considering the geodesic equations of motion for a test particle in Schwarzschild geometry for some time now. Similar to what we can do with the Kepler problem I would like to be able to numerically integrate the equations of motion. I'm quite interested to see how coordinate time and proper time change over a simulation and also the trajectory itself. But the problem is.. I have no idea how to do that. Any suggestions?

 

[Ibix]

See https://www.physicsforums.com/threads/null-geodesics-in-schwarzschild-spacetime.895174/

Note that there's a mistake early on in the thread that gets corrected further down. That's for null geodesics, but it's fairly easy to generalise for time-like geodesics. I haven't tried space-like geodesics but I suspect they're fairly straightforward too.

 

[pervect]

Hello there, I've been considering the geodesic equations of motion for a test particle in Schwarzschild geometry for some time now. Similar to what we can do with the Kepler problem I would like to be able to numerically integrate the equations of motion. I'm quite interested to see how coordinate time and proper time change over a simulation and also the trajectory itself. But the problem is.. I have no idea how to do that. Any suggestions?

The Schwarzschild geodesic has a couple of conserved quantites that mean you don't have to numerically integrate the equations, you can solve them analytically by taking advantage of the existence of these constants of motion.

Sean Caroll discusses this in his lecture notes on GR, though it's not an easy read. The applicable concept is called "Killing Vectors", a quote from the relevant section (found at https://www.preposterousuniverse.com/grnotes/) is:

There does not seem to be much hope for simply solving this set of coupled equations by inspection. Fortunately our task is greatly simplified by the high degree of symmetry of the Schwarzschild metric. We know that there are four Killing vectors: three for the spherical symmetry, and one for time translations. Each of these will lead to a constant of the motion

For another online reference, http://www.fourmilab.ch/gravitation/orbits/ might or might not be easier to follow - it is taken pretty much from MTW's treatment in "Gravitation", which would be a better source than the above (being more detailed and of better provenance) - if you can get a hold of it. This approach is a bit more physical, they give the conserved quantites names, like "energy at infinity" and "angular momentum". The resulting equations are rather similar to how we might solve for motion in the Newtonian case by taking advantage of energy and angular momentum conservation.

It might help to review the "effective potential" method for the discussion in the second link. Orbits occur in a plane, so they're basically two dimensional. So one can introduce unit vectors $\hat{r}$ and $\hat{\phi}$ in the Newtonian case to describe the orbital plane, the former being in the radial direction, the other being at right angles to the radial vector in the orbital plane.

By the conservation of angular momentum, one can write a relationship between the the rate of change of $\phi$ and the radial distance r. This gives the non-radial component of the velocity. Knowing the radial and transverse-to-radial components of the velocity, we can get the total kinetic energy, and using the fact that energy is conserved, we can solve for the radial and transverse components of the velocity at any point.

The "effective potential" method is basically a matter of interpreting this solution as the motion of a 1 dimensional particle in a force field that gives the proper solution for the radial component of the motion. Many texts should discuss the effective potential method, including Goldstein's "Classical Mechanics".

 

[m4r35n357]

[To the OP]
Yes, this is most definitely the way to go. There are some, er, "fascinating" technical issues you will need to deal with, but it's far better than just blasting away at the second order ODEs in the geodesic equation (particularly if you like to cross horizons!).

Oh, couldn't resist. Here is my take on the fourmilab animations, with a spinning black hole (it has one solar mass, and is orbited by a test particle similar to Mercury). Also does Newtonian orbit for "comparison".

 

[Matter_Matters]

Absolutely fantastic responses both very helpful and also a lot I didn't think about. Being a "trained mathematician" I was always satisfied to leave the geodesic equations as they are and have no other concerns i.e. I was completely content with just having a set of DE's and that was that. However, now that my needs require more involvement I find myself almost lost in the technicalities of things. Still though, absolutely fascinating to study. I do find it interesting that such topics are not heavily discussed in the literature, after all, nearly every student would know a member of the Runge-Kutta family. So extending this should be a topic of interest. Don't you agree?

 

[m4r35n357]

Absolutely fantastic responses both very helpful and also a lot I didn't think about. Being a "trained mathematician" I was always satisfied to leave the geodesic equations as they are and have no other concerns i.e. I was completely content with just having a set of DE's and that was that. However, now that my needs require more involvement I find myself almost lost in the technicalities of things. Still though, absolutely fascinating to study. I do find it interesting that such topics are not heavily discussed in the literature, after all, nearly every student would know a member of the Runge-Kutta family. So extending this should be a topic of interest. Don't you agree?

Since you are a trained mathematician (I most definitely am NOT), I will try to be brief ;) The seminal paper on geodesics IMO is this one by Wilkins. The first order equations for the Kerr metric are all there in plain sight in section II.

One other topic I would bring up is the Symplectic Integrator, which is frequently used over RK4 for long term stability.

And finally, a plug for my proudest "contribution" to the field, in case you like to mess with working code.

 

[pervect]

The end result is that if you consider the motion in an equatorial orbital plane (t,r, $\phi$), the solution for writing the geodesic equation is three functions $t(\tau), r(\tau), \phi(\tau)$.

We can define $E = (1-2m/r) dt / d\tau$ and $L = r^2 (d\phi/d\tau)$ and re-write two of these equations as $dE/d\tau=0, dL/d\tau=0$, using the chain rule to expand them. The last equation comes from the Schwarzschild metric itself, which has various forms. One is:

$$c^2 \,{d \tau}^{2} = \left(1 - 2m/r \right) c^2 \,dt^2 - \left(1- 2m/r \right)^{-1} \,dr^2 - r^2 \left(d\theta^2 + \sin^2\theta \, d\phi^2\right)$$

but the fourmilab reference assumes c=1, for instance, and by assuming an equatorial orbit we have previously set $\theta=0$. Making these assumptions and dividing both sides by $d\tau^2$ should basically give the last equation.

I think I posted something more detailed along these lines a long time ago, but I'm not finding it in search.