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

• If you're trying to integrate stable orbits with Runge-Kutta, you're gonna have a bad time. Jun 3 '19 at 0:44
• Are you sure that the orbits starting outside the event horizon don't diverge to infinity, make crazy nonsensical spirals, or oscillate in a circular path? The first two cases are essentially what happens with the Schwarzschild metric. Jun 3 '19 at 0:44
• @HiddenBabel Oh I know the pain, it did work for the Schwarzchild case so I'm assuming it can handle this one as well. Jun 3 '19 at 0:47
• @PeterShor Huh, I was really expecting at least a few of them to be stable. I haven't found many articles talking abou the stability of the orbits so I'm kind of lost. Jun 3 '19 at 0:48
• @ÍcaroLorran How did you get the equation system? Jun 3 '19 at 4:08

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

• ISO-31/XI: According to the ISO regulations and the IPU recommendations, italic symbols should be used only to denote those mathematical and physical entities that may assume different values, typically those symbols that play the role of physical variables.... Any other symbol that was not dealt with in the preceding subsection must be set in roman font Jun 7 '19 at 12:07