Numerical solutions for equatorial orbits in the Kerr black hole Currently, I am trying to find timelike orbits in the Kerr metric around the equator. The problem is that no matter which parameters I choose or the method I use I can't seem to get to physically sound orbits. The solutions I got so far diverge to infinity, make some crazy nonsensical spirals or oscillate in a circular path.
Here's what I'm using currently for the differential equations:
$$r''=\frac{-1}{2r}\left(r'^2+\frac{1}{r}\left(2r(1-T^2)-r_0+\frac{r_0}{r^2}\left(aT+R\right)^2\right)\right)$$
and
$$r'^2 = -\frac{\Delta(r)}{r^2}(1+T\dot{t}+R\dot{\phi}).$$
The symbols $T$ and $R$ are constants, whereas for $\Delta$,$\dot{t}$ and $\dot{\phi}$:
$\Delta(r)=a^2+r^2-rr_0$
$\dot{\phi} = \frac{1}{\Delta}\left[\left(1-\frac{r_0}{r}\right)R-\frac{ar_0}{r}T\right]$
$\dot{t} = \frac{1}{\Delta}\left[-\left(r^2+a^2+\frac{a^2r_0}{r}\right)T-\frac{ar_0}{r}R\right]$
I chose $\phi$ as my azimuthal angle.
I tried to use these on a code I made in Python and another one in Java. Both of them seem to give the same crazy solutions.
If you need to take a look at the code:
https://github.com/icarosadero/black_holes/blob/master/geodesic.java
https://github.com/icarosadero/black_holes/blob/master/script.py
With all of that said I would like to ask whether or not those equations are right. I don't have many people near me available to help at the moment.
 A: The equatorial circular local velocity in natural units of $\rm G=M=c=1$ is
$$\rm v_{\phi}=\frac{a^2  \pm 2a \sqrt{r}+r^2}{\sqrt{a^2+(r-2)r} \ (a \mp r^{3/2})}$$
where the larger solution is retrograde, and the smaller one prograde. This can be obtained by setting
$$\rm \ddot{r}=\dot{r}=\ddot{\theta}=\dot{\theta}=\ddot{\phi}=0 \ , \ \ \theta = \pi/2  \ , \ \ v=v_{\phi}$$
and solve for $\dot{\phi}$. Then you get
$$\rm \ddot{r} = \frac{2 a r^2 \Delta \dot{t} \dot{\phi}-r^2 \left(a^2+(r-2) r\right) \dot{t}^2+\Delta \left(r^5-a^2 r^2\right) \dot{\phi}^2}{r^6} = 0$$
The total time dilation $\rm \dot{t} $ is
$$\rm \dot{t} = \frac{\varsigma}{\sqrt{1-v^2}}$$
with the gravitational time dilation
$$\varsigma = \sqrt{g^{\rm t t}} = \sqrt{\frac{\chi}{\Delta \ \Sigma}}$$
where the abbreviated terms are
$$ \rm \Sigma =a^2 \cos ^2 \theta +r^2 \ , \ \ \chi =\left(a ^2+r^2\right)^2-a ^2 \ \sin ^2 \theta  \ \Delta ,  \ \  \Delta =a^2+r^2-2 r \ $$
The conversion from $\dot{\phi}$ to $\rm v_{\phi}$ is
$$ \rm v_{\phi} = \frac{ L_{z} \sqrt{1 - v^2}}{ \rm \bar R_{\phi} } $$
so in the equatorial plane with $\rm v=v_{\phi}$
$$\rm v_{\phi} = \frac{L_z}{\sqrt{L_z^2+\bar R_{\phi}^2}}$$
where $  \rm \bar R_{\phi} $ is the axial radius of gyration
$$ {\rm \bar R_{\phi}} = \sqrt{|g_{\phi \phi}|} =  \sqrt{\frac{\chi}{\Sigma}} \ \sin \theta $$
and $\rm L_z$ is the conserved axial angular momentum
$$ {\rm L_z} = -g_{\phi \phi} \ \dot \phi- g_{\rm t \phi} \ \dot {\rm t} = {\rm \frac{\sin ^2 \theta \ (\dot{\phi} \ \Delta \ \Sigma - 2 \ a \ E \ r)}{\Sigma -2 \ r}} $$
with the total energy $\rm E$
$${\rm E} = g_{\rm t t}\ \dot {\rm t} + g_{\rm t \phi} \ \dot \phi =  \rm \surd \left(\frac{(\Sigma - 2 \ r) \left(\dot{\theta}^2 \ \Delta \ \Sigma +\dot{r}^2 \ \Sigma + \Delta \right)}{\Delta \ \Sigma }+\dot{\phi}^2 \ \Delta \ \sin ^2 \theta \right)$$
which can in the equatorial plane further be abbreviated by setting all the derivatives which we don't need to $ 0 $ and $\sin \theta=1$, like we did in the equation for $\rm \ddot{r}$. For examples and the complete equations of motion see here, here or here. If you are looking for the innermost stable orbits (ISCO) the equation is
$$\rm r_{\mathrm{isco}} =  3 + Z_2 \pm \sqrt{(3-Z_1)(3+Z_1+2Z_2)}$$
with the shorthand terms
$$\rm Z_1 = 1 + \sqrt[3]{1-a^2} \left( \sqrt[3]{1+a} + \sqrt[3]{1-a} \right)  \ , \ \ \rm Z_2 = \sqrt{3a^2 + Z_1^2}$$
see this reference.

In the upper image you have a prograde orbit (red) compared to a locally stationary ZAMO (dashed magenta) around a spinning black hole with spin parameter $\rm a=1$ at $\rm r=4$. The turquoise dot represents an object which is stationary with respect to the far away observer, for which it must locally counterrotate with the local frame dragging velocity of $\rm c/6$.
