Kerr geodesics differential equations in equatorial plane

With friend, we are writing an interactive educational simulation of particle falling into a black hole.

Currently we use Schwarzschild geodesics. However, we want to generalize it to the case of rotating (and perhaps rotating and charged) black hole. We are mostly interested in the equatorial plane, as then we can plot it on a 2D tablet.

So, what are the differential equations for a particle (with given initial position and velocity) falling in the Kerr (or Kerr-Newman) metrics in the equatorial plane?

I'm interested in an explicit form (plug & play - should work after insertion of the black hole parameters (i.e. $M, L, Q$) and the initial contitions (i.e. $\vec{x}, \vec{v}, q$); $Q$ and $q$ are optional, as Kerr metrics is nice by itself).

Side notes:

Yes, I know the general procedure. Just I'm short of time (so now I'm even no longer coding it). So I may self-answer, but rather later than sooner.

It's almost in Chapter 20 of something: Geodesic motion in Kerr spacetime (i.e. (20.25) and (20.31) for the equations of motion; (20.18) and (20.19) for energy and angular momentum). However, some parameters are not introduced (perhaps there are in the previous chapters...).

-
I don't have the book handy at the moment, but you should (at least, that's where I would) check Chapter 4 of B. O'Neill's The Geometry of Kerr Black Holes –  Willie Wong Nov 7 '12 at 15:25

I'll follow Gravitation by Misner, Thorne, and Wheeler (hereafter MTW), which is the standard reference textbook encyclopedic tome for the field despite its age.

Let $\lambda$ parametrize the path such that the derivative with respect to it gives the 4-momentum. Using Boyer-Lindquist coordinates, MTW Box 33.5 gives $$\left(\frac{\mathrm{d}r}{\mathrm{d}\lambda}\right)^2 = \frac{1}{r^4} \left(\alpha E^2 - 2\beta E + \gamma_0\right),$$ where $$E = \frac{1}{\alpha} \left(\beta + \sqrt{\beta^2 - \alpha\gamma_0 + \alpha r^4(p^r)^2}\right)$$ is a constant of the motion (energy at infinity) and we have $$\alpha = \left(r^2 + a^2\right)^2 - \Delta a^2 \\ \beta = \left(L_z a + qQr\right) \left(r^2 + a^2\right) - L_z a\Delta \\ \gamma_0 = \left(L_z a + qQr\right)^2 - \Delta L_z^2 - m^2r^2\Delta \\ \Delta = r^2 - 2Mr + a^2 + Q^2.$$ Here $m$ is the test particle's rest mass and $L_z$ is its (conserved) angular momentum at infinity. $a = L/M$ is the black holes's angular momentum per unit mass. ($L$ not related to $L_z$ - sorry about that.)

For the azimuthal motion, I turn to MTW Eq. 33.32c, which states (after setting $\theta = \pi/2$) $$\frac{\mathrm{d}\phi}{\mathrm{d}\lambda} = -\frac{1}{r^2} \left(\frac{aP}{\Delta} - aE + L_z\right).$$ Here we define $$P = E \left(r^2 + a^2\right) - L_z a - qQr.$$

The final step is finding the relation between time $t$ (of the Boyer-Lindquist variety, which means it's not crazy) and $\lambda$. MTW Eq. 33.32d tells us (again after setting $\theta = \pi/2$) $$\frac{\mathrm{d}t}{\mathrm{d}\lambda} = \frac{1}{r^2} \left(\frac{P}{\Delta} \left(r^2 + a^2\right) - a^2E + aL_z\right).$$

Hope this helps. I remember coding something similar (alright, plugging the ODEs into Mathematica) once upon a time. It seemed to work reasonably well without needing any fancy numerical techniques to ensure stability... at least for a few orbits, after which I couldn't tell what it was supposed to be doing.

-
If anyone reading this wants a character-building experience, I recommend Exercise 33.7 of MTW, which derives these equations in all their gory detail. I should mention though: I'm reasonably sure there's a sign error or two buried in some of the intermediate equations given in the exercise - don't loose as much sleep as I did when you come across it. –  Chris White Nov 8 '12 at 5:00
Thanks! However, what does "MTW" stand for? (A book I guess.) –  Piotr Migdal Nov 8 '12 at 9:18
Yes, Misner, Thorne and Wheeler. text book "Gravitation" –  Jim Graber Nov 8 '12 at 20:22
@Jim Thanks. Edited. I sometimes forget the book has a title other than those initials. –  Chris White Nov 8 '12 at 22:28
@PiotrMigdal: I would also advise reading a bit about the behavior of the Kerr orbits--in particular, the co-rotating orbits and the anti-co-rotating orbits will behave differently and have different points of instability. This should be a decent check on your code. –  Jerry Schirmer Nov 8 '12 at 22:36