Kerr geodesics differential equations in equatorial plane

Question

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...).

Answer 1

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.