Astrophysically relevant bound trajectories around a Kerr black hole

TL;DR

This paper presents new analytic solutions for non-equatorial eccentric bound trajectories around Kerr black holes, enabling faster computations and detailed parameter space analysis relevant for gravitational wave studies.

Contribution

The authors derive novel closed-form solutions and translation relations for bound trajectories around Kerr black holes, including specialized formulas for spherical and separatrix orbits, advancing the modeling of gravitational wave sources.

Findings

Derived new closed-form solutions for non-equatorial bound trajectories.

Established translation relations between orbital parameters and constants of motion.

Provided exact expressions for frequencies and orbit sketches with astrophysical applications.

Abstract

We derive alternate and new closed-form analytic solutions for the non-equatorial eccentric bound trajectories, $\{ \phi \left( r, \theta \right)$, $\ t \left( r, \theta \right),\ r \left( \theta \right) \}$, around a Kerr black hole by using the transformation $1/r=\mu \left(1+ e \cos \chi \right)$. The application of the solutions is straightforward and numerically fast. We obtain and implement translation relations between energy and angular momentum of the particle, ($E$, $L$), and eccentricity and inverse-latus rectum, ($e$, $\mu$), for a given spin, $a$, and Carter's constant, $Q$, to write the trajectory completely in the ($e$, $\mu$, $a$, $Q$) parameter space. The bound orbit conditions are obtained and implemented to select the allowed combination of parameters ($e$, $\mu$, $a$, $Q$). We also derive specialized formulae for spherical and separatrix orbits. A study of the non-equatorial analog of the previously studied equatorial separatrix orbits is carried out where a homoclinic orbit asymptotes to an energetically bound spherical orbit. Such orbits simultaneously represent an eccentric orbit and an unstable spherical orbit, both of which share the same $E$ and $L$ values. We present exact expressions for $e$ and $\mu$ as functions of the radius of the corresponding unstable spherical orbit, $r_s$, $a$, and $Q$, and their trajectories, for ($Q\neq0$) separatrix orbits; they are shown to reduce to the equatorial case. These formulae have applications to study the gravitational waveforms from EMRIs besides relativistic precession and phase space explorations. We obtain closed-form expressions of the fundamental frequencies of non-equatorial eccentric trajectories that are equivalent to the previously obtained quadrature forms and also numerically match with the equivalent formulae previously derived. We sketch several orbits and discuss their astrophysical applications.

Tables

Table 1. Table 1: Glossary of symbols used.

Boyer Lindquist coordinates
$t$	Time coordinate	$r$	Radial distance from the black hole
$\theta$	Polar angle	$\varphi$	Azimuthal angle
${\rho }^2$	$r^2+a^{2}co{s}^2\theta$	$a$	Spin of the black hole
Common physical parameters
$u$	$\frac{1}{r}$	$\tau$	Proper time
$r_{+}$	Horizon radius	$Q$	Carter’s constant
$E$	Energy per unit rest mass of the test particle	$L$	z component of Angular momentum per unit
			rest mass of the test particle
$p_t$	Generalized momentum for $t$ coordinate	$p_{\varphi }$	Generalized momentum for $\varphi$ coordinate
$p_r$	Generalized momentum for $r$ coordinate	$m_0$	=0 for photon orbits and $=1$ for particle orbits
$V_{eff}$	Radial effective potential for an eccentric test	$\mathscr{H}$	Relativistic Hamiltonian for the geodesic motion
	particle trajectory
$r_1$	apastron distance ($=r_a$)	$r_2$	periastron distance ($=r_p$)
$r_3$	Third turning point of the test particle	$r_4$	Innermost turning point of the test particle
$e$	eccentricity parameter	$\mu$	inverse latus-rectum parameter
Integrals of motion
$\chi$	defined by $u=\mu \left(1+e\cos \chi \right)$	$\psi$	$\frac{\chi }{2}-\frac{\pi }{2}$
$y$	$1+e\cos \chi$	$I$	Terminology used for radial integrals
$H$	Terminology used for $\theta$ integrals
Spherical and separatrix orbits
$r_s$	radius of spherical orbit	$r_c$	radius of circular orbit
$e_s$	eccentricity of the separatrix orbits	${\mu }_s$	inverse latus-rectum of the separatrix orbits
$Z$	ISCO radius	$X$	Light radius
Fundamental frequencies
${\nu }_{\varphi }$	Azimuthal frequency	${\nu }_r$	Radial frequency
${\nu }_{\theta }$	Vertical oscillation frequency

Table 2. Table 2: This table summarizes all the integrals solved in § 2.2 , 2.4 to calculate the integrals of motion in the Kerr geometry, where all the constants are defined in the text.

	Analytic solution of $\varphi$ and $t$ for $Q\ne 0$
	$\varphi -{\varphi }_0=\frac{1}{2}\left(C_3{I}_3+C_4{I}_4-H_1\right)$; $t-t_0=\frac{1}{2}\left(C_5{I}_5+C_6{I}_6+C_7{I}_3+C_8{I}_4+H_2\right)$
$I_1$	${\int }_{r_a}^r\frac{1}{\Delta \sqrt{R}}\frac{\partial R}{\partial L}d{r}^{{}^{\prime }}=-\left[C_3{I}_3+C_4{I}_4\right]$
$I_2$	${\int }_{r_a}^r\frac{1}{\Delta \sqrt{R}}\frac{\partial R}{\partial E}d{r}^{{}^{\prime }}=\left[C_5{I}_5+C_6{I}_6+C_7{I}_3+C_8{I}_4\right],$
$I_3$	${\int }_0^{\psi }\frac{d\psi }{\left(1+p_2^2{\sin }^2\psi \right)\sqrt{1-m^2{\sin }^2\psi }\sqrt{1-n^2{\sin }^2\psi }}=\frac{1}{\sqrt{1-m^2}\left(m^2+p_2^2\right)}\left[m^{2}F(\alpha ,k^2)+p_2^2\Pi (\frac{-p_2^2-m^2}{1-m^2},\alpha ,k^2)\right]$
	where $\sin \alpha =\frac{\sqrt{1-m^2}\sin \psi }{\sqrt{1-m^2{\sin }^2\psi }}$, $\psi =\frac{\chi }{2}-\frac{\pi }{2}$ and $1/r=\mu \left(1+e\cos \chi \right)$
$I_4$	${\int }_0^{\psi }\frac{d\psi }{\left(1+p_3^2{\sin }^2\psi \right)\sqrt{1-m^2{\sin }^2\psi }\sqrt{1-n^2{\sin }^2\psi }}=\frac{1}{\sqrt{1-m^2}\left(m^2+p_3^2\right)}\left[m^{2}F(\alpha ,k^2)+p_3^2\Pi (\frac{-p_3^2-m^2}{1-m^2},\alpha ,k^2)\right]$
$I_5$	${\int }_0^{\psi }\frac{d\psi }{{\left(1+p_1^2{\sin }^2\psi \right)}^2\sqrt{1-m^2{\sin }^2\psi }\sqrt{1-n^2{\sin }^2\psi }}$
	$=\frac{1}{\sqrt{1-m^2}{\left(m^2+p_1^2\right)}^2}\left[m^{4}F(\alpha ,k^2)+2p_1^2{m}^2\Pi (s^2,\alpha ,k^2)+p_1^4{I}_7(\alpha ,e,\mu ,a,Q)\right],$
$I_6$	${\int }_0^{\psi }\frac{d\psi }{\left(1+p_1^2{\sin }^2\psi \right)\sqrt{1-m^2{\sin }^2\psi }\sqrt{1-n^2{\sin }^2\psi }}=\frac{1}{\sqrt{1-m^2}\left(m^2+p_1^2\right)}\left[m^{2}F(\alpha ,k^2)+p_1^2\Pi (s^2,\alpha ,k^2)\right]$
$I_7$	${\int }_0^{\psi }\frac{d\psi }{{\left(1+s^2{\sin }^2\psi \right)}^2\sqrt{1-k^2{\sin }^2\psi }}=\frac{s^4\sin \alpha \cos \alpha \sqrt{1-k^2{\sin }^2\alpha }}{2\left(1-s^2\right)\left(k^2-s^2\right)\left(1-s^2{\sin }^2\alpha \right)}+\frac{\left[s^4-2s^2\left(1+k^2\right)+3k^2\right]}{2\left(1-s^2\right)\left(k^2-s^2\right)}\Pi (s^2,\alpha ,k^2)$
	$-\frac{1}{2\left(1-s^2\right)}F(\alpha ,k^2)-\frac{s^2}{2\left(1-s^2\right)\left(k^2-s^2\right)}K(\alpha ,k^2)$
$H_1$	${\int }_{{\theta }_0}^{\theta }\frac{1}{\sqrt{\Theta }}\frac{\partial \Theta }{\partial L}d{\theta }^{{}^{\prime }}=\frac{2L}{z_{+}a\sqrt{1-E^2}}\{\Pi (z_{-}^2,\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})-\Pi (z_{-}^2,\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})-$
	$F(\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})+F(\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})\}$
$H_2$	${\int }_{{\theta }_0}^{\theta }\frac{1}{\sqrt{\Theta }}\frac{\partial \Theta }{\partial E}d{\theta }^{{}^{\prime }}=\frac{2Ea{z}_{+}}{\sqrt{1-E^2}}\{K(\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})-K(\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})-$
	$F(\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})+F(\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2})\}.$
	Analytic solution of $\varphi$ and $t$ for $Q=0$
	$\varphi -{\varphi }_0=a_1\Pi (-p_2^2,\psi ,m^2)+b_1\Pi (-p_3^2,\psi ,m^2)$
	$t-t_0=a_2\left[\frac{p_1^4\sin \psi \cos \psi \sqrt{1-m^2{\sin }^2\psi }}{2\left(1+p_1^2\right)\left(m^2+p_1^2\right)\left(1+p_1^2{\sin }^2\psi \right)}-\frac{F(\psi ,m^2)}{2\left(1+p_1^2\right)}+\frac{p_1^{2}K(\psi ,m^2)}{2\left(1+p_1^2\right)\left(m^2+p_1^2\right)}\right]+d_2\Pi (-p_3^2,\psi ,m^2)$
	$+c_2\Pi (-p_2^2,\psi ,m^2)+\Pi (-p_1^2,\psi ,m^2)\left\{a_2\frac{\left[p_1^4+2p_1^2\left(1+m^2\right)+3m^2\right]}{2\left(1+p_1^2\right)\left(m^2+p_1^2\right)}+b_2\right\}$

Table 3. Table 3: This table summarizes the trajectory solution derived in § 3.3 for the non-equatorial separatrix orbits.

	Analytic solutions
$S_3$	$=\frac{1}{\left(1+p_2^2\right)}\left[\frac{p_2^2}{\sqrt{-\left(p_2^2+m^2\right)}}\ln \sqrt{\frac{\sqrt{1-m^2}+\sqrt{-\left(p_2^2+m^2\right)}\sin \alpha }{\sqrt{1-m^2}-\sqrt{-\left(p_2^2+m^2\right)}\sin \alpha }}+\frac{\ln \left(\tan \alpha +\sec \alpha \right)}{\sqrt{1-m^2}}\right]$
$S_4$	$=\frac{1}{\left(1+p_3^2\right)}\left[\frac{p_3^2}{\sqrt{-\left(p_3^2+m^2\right)}}\ln \sqrt{\frac{\sqrt{1-m^2}+\sqrt{-\left(p_3^2+m^2\right)}\sin \alpha }{\sqrt{1-m^2}-\sqrt{-\left(p_3^2+m^2\right)}\sin \alpha }}+\frac{\ln \left(\tan \alpha +\sec \alpha \right)}{\sqrt{1-m^2}}\right]$
$S_5$	$=\frac{1}{\sqrt{1-m^2}{\left(m^2+p_1^2\right)}^2}\left[\frac{m^2\left(m^2-p_1^2{m}^2+2p_1^2\right)}{\left(1+p_1^2\right)}\ln \left(\tan \alpha +\sec \alpha \right)+p_1^4{S}_7+\frac{2p_1^2{m}^2\left(1-m^2\right)}{\left(1+p_1^2\right)}\mid s\mid {\tan }^{-1}\left[\mid s\mid \sin \alpha \right]\right]$
$S_6$	$=\frac{\ln \left(\tan \alpha +\sec \alpha \right)}{\sqrt{1-m^2}\left(1+p_1^2\right)}+\frac{p_1^2\sqrt{1-m^2}}{\left(m^2+p_1^2\right)\left(1+p_1^2\right)}\mid s\mid {\tan }^{-1}\left[\mid s\mid \sin \alpha \right]$
$S_7$	$=\frac{1}{2{\left(1-s^2\right)}^2}\left[\frac{s^4\sin \alpha {\cos }^2\alpha }{\left(1-s^2{\sin }^2\alpha \right)}+2\ln \left(\tan \alpha +\sec \alpha \right)-s^2\sin \alpha +\left(3-s^2\right)\mid s\mid {\tan }^{-1}\left(\mid s\mid \sin \alpha \right)\right]$

Table 4. Table 4: This following table summarizes the values of conic parameters ( e 𝑒 e , μ 𝜇 \mu ) chosen in the listed orbit simulations to study eccentric, homoclinic and spherical orbits for different ( a 𝑎 a , Q 𝑄 Q ) combinations for both prograde and retrograde cases constructed using Eqs. ( 12 ).

Type of	Orbit $\#$	Inverse latus-	Eccentricity	Spin of	Carter’s constant	Varying
orbit		rectum of the orbit	of the orbit	the black hole		parameter
		$\mu$	$e$	$a$	$Q$
Eccentric orbits	E1	0.1	0.6	0.2	3	$a$
	E2	0.1	0.6	0.8	3
	E3	0.1	0.6	0.2	8	$a$
	E4	0.1	0.6	0.8	8
	E5	0.1	0.6	-0.2	3	$a$
	E6	0.1	0.6	-0.8	3
	E7	0.1	0.6	-0.2	8	$a$
	E8	0.1	0.6	-0.8	8
Homoclinic orbits	H1	0.153	0.6	0.2	3	$a$ and $e$
	H2	0.208	0.2	0.5	3
	H3	0.153	0.5	0.2	8	$a$
	H4	0.172	0.5	0.5	8
	H5	0.127	0.5	-0.2	3	$a$ and $e$
	H6	0.127	0.2	-0.5	3
	H7	0.134	0.5	-0.2	8	$a$
	H8	0.123	0.5	-0.5	8
Spherical orbits	S1	0.222	0	0.5	3	$a$ and $Q$
	S2	0.144	0	-0.5	8
Zoom-whirl	Z1	0.155	0.5	0.2	5	$a$
	Z2	0.226	0.5	0.8	5
	Z3	0.142	0.8	0.2	5	$a$
	Z4	0.212	0.8	0.8	5
	Z5	0.162	0.5	0.5	10	$a$
	Z6	0.179	0.5	0.8	10

Table 5. Table 5: This table summarizes the fundamental frequency formulae derived using the long time average method in the Kerr geometry. The explicit expressions for integrals I 1 subscript 𝐼 1 I_{1} , I 2 subscript 𝐼 2 I_{2} and I 8 subscript 𝐼 8 I_{8} are summarized in table 2 .

${\nu }_r(e,\mu ,a,Q)$	$\frac{F(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})}{M\left\{\left[I_2(\frac{\pi }{2},e,\mu ,a,Q)+2E{I}_8(\frac{\pi }{2},e,\mu ,a,Q)\right]F(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})-2a^2{z}_{+}^{2}E{I}_8(\frac{\pi }{2},e,\mu ,a,Q)K(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})\right\}}$
${\nu }_{\theta }(e,\mu ,a,Q)$	$\frac{a\sqrt{1-E^2}{z}_{+}{I}_8(\frac{\pi }{2},e,\mu ,a,Q)}{2M\left\{\left[I_2(\frac{\pi }{2},e,\mu ,a,Q)+2E{I}_8(\frac{\pi }{2},e,\mu ,a,Q)\right]F(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})-2a^2{z}_{+}^{2}E{I}_8(\frac{\pi }{2},e,\mu ,a,Q)K(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})\right\}}$
${\nu }_{\varphi }(e,\mu ,a,Q)$	$\frac{\left\{\left[-I_1(\frac{\pi }{2},e,\mu ,a,Q)-2L{I}_8(\frac{\pi }{2},e,\mu ,a,Q)\right]F(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})+2L{I}_8(\frac{\pi }{2},e,\mu ,a,Q)\Pi (z_{-}^2,\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})\right\}}{2\pi M\left\{\left[I_2(\frac{\pi }{2},e,\mu ,a,Q)+2E{I}_8(\frac{\pi }{2},e,\mu ,a,Q)\right]F(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})-2a^2{z}_{+}^{2}E{I}_8(\frac{\pi }{2},e,\mu ,a,Q)K(\frac{\pi }{2},\frac{z_{-}^2}{z_{+}^2})\right\}}$

Equations

$$\mathrm{d}{s}^2=-\left(1-\frac{2r}{{\rho }^2}\right)\mathrm{d}{t}^2-\frac{4ar{\sin }^2\theta }{{\rho }^2}\mathrm{d}\varphi \mathrm{d}t+\frac{{\rho }^2}{r^2-2r+a^2}\mathrm{d}{r}^2+{\rho }^2\mathrm{d}{\theta }^2$$

$$+\left(r^2+a^2+\frac{2r{a}^2{\sin }^2\theta }{{\rho }^2}\right){\sin }^2\theta \mathrm{d}{\varphi }^2,$$

$$\mathcal{H}=$$

$$=$$

$$m_0{\rho }^2\frac{\mathrm{d}t}{\mathrm{d}\tau }=$$

$$m_0{\rho }^2\frac{\mathrm{d}r}{\mathrm{d}\tau }=$$

$$m_0{\rho }^2\frac{\mathrm{d}\theta }{\mathrm{d}\tau }=$$

$$m_0{\rho }^2\frac{\mathrm{d}\varphi }{\mathrm{d}\tau }=$$

$$R=$$

$$\Theta =$$

$$\tau -{\tau }_0=$$

$$\varphi -{\varphi }_0=$$

$$t-t_0=$$

$${\int }_{r_0}^r\frac{\mathrm{d}{r}^{{}^{\prime }}}{\sqrt{R}}=$$

$$\frac{\left(E^2-m_0^2\right)}{2}=\frac{m_0^2{\rho }^4}{2r^4}{\left(\frac{\mathrm{d}r}{\mathrm{d}\tau }\right)}^2+V_{eff}\left(r,a,E,L,Q\right),$$

$$V_{eff}\left(r,a,E,L,Q\right)=-\frac{m_0^2}{r}+\frac{L^2-a^2\left(E^2-m_0^2\right)+Q}{2r^2}-\frac{{\left(L-aE\right)}^2+Q}{r^3}+\frac{a^{2}Q}{2r^4}.$$

$$E\left(e,\mu ,a,Q\right)={\left[1-{\mu }^3{\left(1-e^2\right)}^2\left(\mu a^{2}Q-Q-x^2\right)-\mu \left(1-e^2\right)\right]}^{1/2},$$

$$x^2\left(e,\mu ,a,Q\right)=\frac{-S-\sqrt{S^2-4PR}}{2P},$$

$$P\left(e,\mu ,a,Q\right)=$$

$$S\left(e,\mu ,a,Q\right)=$$

$$R\left(e,\mu ,a,Q\right)=$$

$$I_1\left(\alpha ,e,\mu ,a,Q\right)=-\left[C_3{I}_3\left(\alpha ,e,\mu ,a,Q\right)+C_4{I}_4\left(\alpha ,e,\mu ,a,Q\right)\right],$$

$$I_2\left(\alpha ,e,\mu ,a,Q\right)$$

$$I_3\left(\alpha ,e,\mu ,a,Q\right)=$$

$$I_4\left(\alpha ,e,\mu ,a,Q\right)=$$

$$I_5\left(\alpha ,e,\mu ,a,Q\right)$$

$$I_6\left(\alpha ,e,\mu ,a,Q\right)$$

$$I_7\left(\alpha ,e,\mu ,a,Q\right)$$

$$I_8\left(\alpha ,e,\mu ,a,Q\right)=$$

$$C_3=$$

$$C_4=$$

$$C_5=$$

$$C_7=$$

$$C_8=$$

$$A=$$

$$B=$$

$$C=$$

$$n^2=$$

$$k^2=$$

$$p_1^2=$$

$$x_{1,2}=$$

$$F\left(\alpha ,k^2\right)=$$

$$K\left(\alpha ,k^2\right)=$$

$$\Pi \left(s^2,\alpha ,k^2\right)=$$

$$H_1\left(\theta ,{\theta }_0,e,\mu ,a,Q\right)=\frac{2L}{z_{+}a\sqrt{1-E^2}}\{F\left(\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)-F\left(\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)+$$

$$\Pi \left(z_{-}^2,\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)-\Pi \left(z_{-}^2,\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)\},$$

$$H_2\left(\theta ,{\theta }_0,e,\mu ,a,Q\right)=\frac{2Ea{z}_{+}}{\sqrt{1-E^2}}\{K\left(\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)-F\left(\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)-$$

$$K\left(\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)+F\left(\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)\},$$

$$H_3\left(\theta ,{\theta }_0,e,\mu ,a,Q\right)=\frac{1}{a\sqrt{1-E^2}{z}_{+}}\left\{F\left(\arcsin \left(\frac{\cos {\theta }_0}{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)-F\left(\arcsin \left(\frac{\cos \theta }{z_{-}}\right),\frac{z_{-}^2}{z_{+}^2}\right)\right\},$$

Full text

††: Class. Quantum Grav.

1 Introduction

It has now been established with observational evidence that black holes with masses ranging from $4$-$20M_{\odot}$ in X-ray binaries, to $10^{6}$-$10^{9}M_{\odot}$ in galactic nuclei, are ubiquitous. One among such important evidences is the recent detection of gravitational waves from the black hole binary merger [1] and more such events are awaited to be detected by the planned LISA mission [2]. One of the major objectives of the LISA mission is the detection of gravitational waves from EMRIs, most probably to be sourced from the compact objects spiralling and finally plunging onto the super-massive black hole (SMBH) in galactic nuclei. The dynamics of EMRIs is widely accepted as representative of test-particle motion, evolving adiabatically, in the spacetime of a rotating black hole. Understanding of such strong gravity regimes involves them using the Kerr metric [3], which is a vacuum solution of Einstein’s equation for a rotating black hole. The study of trajectories around the black holes is critical to our understanding of the physical processes and their observational consequences [4].

The trajectories in Kerr and Schwarzschild [5] geometries have been studied extensively. Some of these results are covered in a pioneering work by [6] in an elegant manner. The key idea that the general trajectory in Kerr background can be expressed in terms of quadratures, was first given in [7]. In [8], the energy, $E$, and angular momentum, $L$, were expressed in terms of the circular orbit radius and the spin parameter $a$; the specific solution for the radius of the innermost stable circular orbit (ISCO) was also derived. The necessary conditions for bound geodesics for spherical orbits and the dragging of nodes along the direction of spin of a black hole was discussed [9]. The formulae have proved to be extremely useful in predicting observables in astrophysical applications like accretion around the black holes. For example, a general solution for a star in orbit around a rotating black hole was expressed in terms of quadratures [10] using the formulation given by [7]; the resulting integrals have been calculated numerically. The general expression in terms of quadratures for fundamental orbital frequencies $\nu_{\theta}$, $\nu_{\phi}$ and $\nu_{r}$, for a general eccentric orbit, were first derived by [11], where different cases for circular and equatorial orbits are also discussed but complete analytic trajectories were not calculated. An exact solution for non-spherical polar trajectories in Kerr geometry and an exact analytic expression for $t(r)$ for eccentric orbits in the equatorial plane were derived [12]. These were used to obtain the expressions for the periapsis advance and Lense-Thirring frequencies. The time-like geodesics were expressed in terms of quadratures involving hyper-elliptic, elliptic and Abelian integrals for Kerr and Kerr-(anti) de Sitter spacetimes including cosmological constant [13] and applied in a semi-analytic treatment of Lense Thirring effect.

A time-like orbital parameter $\lambda$ called Mino time [14] was introduced to decouple the $r$ and $\theta$ equations, which was then used to express a wider class of trajectory functions in terms of the orbital frequencies $\nu_{\theta}$, $\nu_{\phi}$ and $\nu_{r}$ [15]. These methods are applied to calculate closed-form solutions of the trajectories and their orbital frequencies [16], using the roots of the effective potential. However, the solutions are expressed in terms of Mino time. The commensurability of radial and longitudinal frequencies, their resonance conditions for orbits in Kerr geometry, and their location in terms of spin and orbital parameter values were studied using numerical implementations of Carlson’s elliptic integrals [17]. Considering the problem of the precession of a test gyroscope in the equatorial plane of Kerr geometry, the analytic expressions to transform energy, angular momentum of the orbiting test particle, and spin of the black hole ($E$, $L$, $a$) to eccentricity, inverse-latus rectum of the bound orbit ($e$, $\mu$, $a$) parameters were obtained [18]. The expressions for radial and orbital frequencies are obtained to the order $e^{2}$ for bound orbits and analytically for the marginally bound homoclinic orbits [19]. The dynamical studies of an important family of Kerr orbits called separatrix or homoclinic orbits are important for computing the transition of spiralling to plunge in EMRIs emitting gravitational waves [20, 21]. The test particles (compact objects in this case) transit through an eccentric separatrix orbit in EMRIs, while progressing adiabatically, before they merge with the massive black hole.

This paper is an expanded version of the published article [23]. In this paper, we study the generic bound trajectories, which are eccentric and inclined, around a Kerr black hole, and then we investigate the non- equatorial separatrix orbits as a special case. We have solved the equations of motion and produce alternate and new closed-form solutions for the trajectories in terms of elliptic integrals without using Mino time, \{\phi\left(r,\theta\right),\ t\left(r,\theta\right),r\left(\theta\right)\ $\mathrm{or}\ \theta\left(r\right)\}$, which makes them numerically faster. We also implement the essential bound orbit conditions to choose the parameters ($e$, $\mu$, $a$, $Q$) of an allowed bound orbit, derived from the essential conditions on the parameters for the elliptic integrals involved in the trajectory solutions. We find that the $e-\mu$ space is more convenient for integrating the equations of motion as the turning points of the integrands are naturally specified in terms of the bound orbit conditions. The exact solutions for the trajectories are found in terms of not overly long expressions involving elliptic functions. We implement the translation formulae between ($E$, $L$) and ($e$, $\mu$) parameters that help us to express the trajectory solutions completely in the ($e$, $\mu$, $a$, $Q$) parameter space which we call the conic parameter space. We then study the case of non-equatorial separatrix trajectories in the conic parameter space. First, we write the essential equations for the important radii like innermost stable spherical orbit ($ISSO$), marginally bound spherical orbit ($MBSO$), and spherical light radius. Similar to the equatorial separatrix orbits, the non-equatorial separatrix or homoclinic trajectories asymptote to an energetically bound unstable spherical orbit, where the spherical orbit radius can vary between $MBSO$ and $ISSO$. We show that the formulae for ($e$, $\mu$) for the non-equatorial separatrix orbits can be expressed as functions of the radius of the corresponding spherical orbit, $r_{s}$, $a$, and $Q$, which also reduce to their equatorial counterpart [20] by implementing the limit $Q\rightarrow 0$. These formulae are obtained by using the expressions of $E$ and $L$ for the spherical orbits. Next, we derive the exact solutions for the non-equatorial separatrix trajectories by reducing our general eccentric trajectory solutions to this case. These solutions are important for investigating the behaviour of gravitational waveforms emitted by inspiralling and inclined test objects near non-equatorial separatrix trajectories in the case of EMRIs.

The ab-initio specification of the allowed geometry of bound orbits in the parameter space is crucial for the calculation of the orbital trajectories and its frequencies. These criteria are used in building, studying and sketching different types of trajectories around a Kerr black hole: for instance, spherical, non-equatorial eccentric, non-equatorial separatrix and zoom-whirl orbits, using our closed-form expressions for trajectories are constructed. We also derive closed-form analytic expressions for the fundamental frequencies of the general non-equatorial trajectories as functions of elliptic integrals around the Kerr black holes. We use a time-averaging method on the first-order equations of motion to derive these frequencies and show that our closed-form analytic expressions of frequencies match with the formulae given by [11] which were left in quadrature forms. We also reduce the general forms to the equatorial case, which is also a new form that is easier to implement and faster by a factor of $\sim 20$.

This paper is organized as follows (see Fig. 1): in §2, we review the basic equations describing $\{r,\ \theta,\ \phi,\ t\}$ motion around the Kerr black hole using Hamiltonian dynamics. In §2.1, we write the translation formulae from $(e,\ \mu,\,a,\ Q)$ to $(E,\ L,\ a,\ Q)$ parameter space. In §2.2, we derive the exact closed-form solutions for the trajectories by solving all involving integrals and writing them in terms of elliptic integrals, \{\phi\left(r,\theta\right),\ t\left(r,\theta\right),r\left(\theta\right)\ $\mathrm{or}\ \theta\left(r\right)\}$. In §2.3, we give the essential bound orbit conditions on $(e,\ \mu,\,a,\ Q)$ parameters applicable to the astrophysical situations. In §2.4, we reduce the analytic solutions to the case of equatorial plane. In §3.1, we derive the formulae for $E$ and $L$ for spherical orbits as a function of radius $r_{s}$, $a$, and $Q$. In §3.2, we write the equations for the radii $ISSO$, $MBSO$, and spherical light radius. We then derive the exact expressions for $e$ and $\mu$ for the non-equatorial separatrix trajectories. In §3.3, we derive the trajectory solutions for the non-equatorial separatrix orbits. In §4, we sketch and discuss various bound trajectories around the Kerr black hole. In §5, we derive the closed-form expressions of the fundamental frequencies in terms of elliptic integrals by the long time averaging method without using Mino time. In §6, we conduct consistency checks by reducing the separatrix trajectories to the equatorial case, and also match the azimuthal to polar frequency ratio, $\nu_{\phi}/\nu_{\theta}$, with the spherical orbits case derived by [9]. We discuss possible applications of our trajectory solutions and frequency formulae in §7. We summarize and conclude our results in §8 and §9 respectively. In Table 1 a glossary of symbols is provided.

2 Integrals of motion and bound orbits around Kerr black hole

In this section, we first set up the basic equations defining the integrals of motion of the general eccentric orbit with $Q\neq 0$ around a Kerr black hole. We then write the formulae defining the transformation from conic parameters ($e$, $\mu$) to dynamical parameters ($E$, $L$) for bound orbits and also provide the conditions for the selection of the parameters ($e$, $\mu$, $a$, $Q$) for the bound orbits. These results are essential for expressing the integrals of motion in ($e$, $\mu$, $a$, $Q$) space for bound orbits. Finally, we derive and present an alternate and simple form of the analytic solutions for the integrals of motion in terms of standard elliptic integrals using the transformation $1/r=\mu\left(1+e\cos\chi\right)$. Such transformations lead to a compact and useful trajectory solution for the non-equatorial and eccentric orbits around a rotating black hole. Later, we reduce these results to a simpler form for equatorial eccentric trajectories.

Considering the Kerr metric [3] for a rotating black hole of mass $M$ in the Boyer-Lindquist coordinates, $x^{\alpha}=(t,r,\theta,\phi)$, in geometrical units $G=c=1$

[TABLE]

where $a=J/M^{2}$ is the specific angular momentum of the black hole and $\rho^{2}=r^{2}+a^{2}\cos^{2}\theta$. We have written the variables $\rho,r$, and $t$ are in units of $M$. The relativistic and conservative Hamiltonian for the geodesic motion of a test particle in Kerr spacetime [7, 24, 25]:

[TABLE]

where $m_{0}$ is the particle’s rest mass, $p_{\beta}$ are the conjugate momenta associated with particle’s coordinates. To derive the complete set of constants of motion, a canonical transformation, ($q^{\alpha}$, $p_{\beta}$)$\rightarrow$($Q^{\alpha}$, $P_{\beta}$), can be found such that the Hamiltonian is cyclic and the new set of momenta, $P_{\beta}$, are conserved along the world-line of the particle. A characteristic function is obtained to generate such transformation and Hamilton’s equations are used to obtain the first-order equations of motion [7, 24, 25]:

[TABLE]

We have written the variables $\rho,r$, and $t$ in the units of $M$. The integrals of motion have also been derived to be [7, 11]

[TABLE]

The equation for the radial motion around the Kerr black hole, Eq. (3b), can be expressed in the form

[TABLE]

where the term on the LHS represents the total energy, the first term on the RHS represents the radial kinetic energy and the second term on the RHS represents the radial effective potential given by

[TABLE]

2.1 Translation relations between ($E$, $L$) and ($e$, $\mu$)

We present the transformation of energy, angular momentum, and Carter’s constant ($E$, $L$, $Q$) space of the test particle to the eccentricity, inverse-latus rectum ($e$, $\mu$, $Q$) space of its corresponding bound orbit. These relations can be derived if $R(r)$ is factorized and the periastron $r_{p}$ and apastron $r_{a}$ of the orbit are substituted as $1/\mu\left(1+e\right)$ and $1/\mu\left(1-e\right)$ respectively. Hence, the formulae connecting ($E$, $L$) and ($e$, $\mu$) parameters for bound orbits (a derivation of these formulae is given in A) are

[TABLE]

These expressions are used to derive analytic results for the integrals of motion, given in §2.2, completely in the ($e$, $\mu$, $a$, $Q$) parameter space.

2.2 Analytic solutions of integrals of motion

Next, we solve for the integrals of motion, i.e. Eqs. (4b-4d) and reduce them to a simple form involving elliptic integrals. We first derive the expressions for the radial integrals $I_{1}$ and $I_{2}$. We assume the starting point of the radial motion to be apastron point of the bound orbit, $r_{0}=r_{a}$. The steps taken to obtain the reduced form of the radial integrals are as follows:

We make the substitution $1/r^{{}^{\prime}}=\mu\left(1+e\cos\chi\right)$ and implement the method of partial fractions. 2. 2.

Then make the substitutions, $\cos\chi=2\cos^{2}\frac{\chi}{2}-1$ and $\psi=\frac{\chi}{2}-\frac{\pi}{2}$. 3. 3.

Implement the variable transformation given by $\sin\alpha=\frac{\sqrt{1-m^{2}}\sin\psi}{\sqrt{1-m^{2}\sin^{2}\psi}}$, where $m^{2}$ is defined by Eq. (9i).

As a result the integrals of motion are expressed as functions of standard elliptic integrals, given by

[TABLE]

where

[TABLE]

and where the variables $E$, $L$ and $x$ can be written as functions of ($e$, $\mu$, $a$, $Q$) using Eqs. (7a-7e), which makes all the integrals to be only functions of ($e$, $\mu$, $a$, $Q$). The definition of the elliptic integrals involved, is given below [26]:

[TABLE]

A complete derivation of these integrals is given in B. Next, to solve the integrals $H_{1}$, $H_{2}$, and $H_{3}$ of Eqs. (4b- 4d), we make the substitutions $z=\cos\theta^{{}^{\prime}}$ and $z=z_{-}\sin\beta$ [16] which reduces these integrals to

[TABLE]

Hence, the equations of motion can now be written in short as

[TABLE]

where $I_{3}$, $I_{4}$, $I_{5}$, $I_{6}$, $I_{8}$, $H_{1}$, $H_{2}$ and $H_{3}$ are given by Eqs. (8-8h; 11c) respectively. Hence, all the integrals are written explicitly as functions of parameters ($e$, $\mu$, $a$, $Q$) through variables $\alpha\left(e,\mu,a,Q;\chi\right)$ and $\beta\left(e,\mu,a,Q;\theta\right)$ which are directly used to calculate ($r$, $\theta$, $t$) through Eqs. (12a-12c). The radial motion, which varies with the $\alpha$, is assumed to have the starting point at the apastron distance, $r_{a}$ or $\alpha=0$, of the orbit and the starting point of the polar motion, $\beta_{0}$ or $\theta_{0}$, is an extra variable which can be chosen in the range $\{\theta_{-},\pi-\theta_{-}\}$. This is tantamount to shifting the starting point of the motion in time or adjusting the initial value of the observed time, $t_{0}$.

Once the initial points are fixed ($\alpha=0$, $\theta=\theta_{0}$), Eqs. (12a) and (12b) are used to calculate $\phi\left(r,\theta\right)$ and $t\left(r,\theta\right)$ respectively, whereas Eq. (12c) gives $r\left(\theta\right)$ or $\theta\left(r\right)$, which can be used to obtain $t\left(r\right)$ or $t\left(\theta\right)$ and $\phi\left(r\right)$ or $\phi\left(\theta\right)$.

The elegant alternate forms presented here help us to write useful and simpler expressions of ($\phi$, $t$) for the equatorial eccentric trajectories, as shown later in §2.4. Also, these results can be used to reduce the radial integrals for non-equatorial separatrix trajectories in the form of logarithmic and trigonometric functions, presented in §3.3, which are useful in the study of gravitational waves from EMRIs.

2.3 Bound orbit conditions in conic parameter space

The bound orbit regions have been studied and divided in the ($E$, $L$) space according to the different types of possible $r$ motion [27]. The most relevant astrophysical bound orbit region corresponds to the case where $E<1$ and there are four real roots of $R(r)$, $r_{1}>r_{2}>r_{3}>r_{4}>0$, such that the bound orbit either exists between $r_{1}$ and $r_{2}$ or $r_{3}$ and $r_{4}$, this has been defined as region III in the ($E$, $L$) plane by [27]. Since $r_{1}$ and $r_{2}$ are the outer most turning points of the effective potential, the bound orbit should exist between these two in the astrophysical situations. We can implement this constraint in the ($e$, $\mu$, $a$, $Q$) space by imposing the condition $k^{2}<1$ on the parameter $k$ used in the radial integrals in §2.2, where we have assumed that a bound orbit exists between $r_{1}$ and $r_{2}$, which requires $k^{2}<1$ as an essential condition for the elliptic integrals to have real values, Eqs. (10a-10c). This further implies

[TABLE]

Hence, the necessary and independent conditions for this region can be collectively given as

[TABLE]

There exists unstable bound orbits for $E>1$ specified as region IV in the ($E$, $L$) plane by [27], where the bound orbit exists between $r_{2}$ and $r_{3}$. Such a situation is not important from the astrophysical point of view, because the particle will follow a bound trajectory between the outermost turning points, i.e. $r_{1}$ and $r_{2}$, and hence the above conditions, Eq. (14), together represent a necessary and sufficient condition for the existence of bound orbits.

2.4 Equatorial bound orbits

In this section, we apply the integrals of motion, Eqs. (12a, 12b), to the eccentric equatorial trajectories, where $Q=0$ ($\theta=\pi/2$). We show that the forms derived in §2.2 reduce to very compact expressions of ($\phi$, $t$) involving trigonometric functions and elliptic integrals for the equatorial eccentric orbits. We implement the limit, $Q\rightarrow 0$ which leads to $A\rightarrow 0$, Eq. (9f), and reduces the factors $\left(1+x_{1}\right)$, $A\left(1+x_{2}\right)$, using Eq. (9l), to

[TABLE]

The substitution of these reduced expressions of $m^{2}$ and $n^{2}$ further simplifies the integrals $I_{3}$, $I_{4}$, $I_{5}$, and $I_{6}$, as shown in C, which finally yields the expressions for azimuthal angle and time coordinate for equatorial trajectories to be given by

[TABLE]

The brackets of the factor $\left[1-\mu^{2}x^{2}\left(3-e^{2}-2e\right)\right]$ in the expressions of $c_{2}$ and $d_{2}$ above were missing in the journal version [23], and has been corrected here. The corresponding fundamental frequency formulae for the equatorial trajectories are

[TABLE]

These compact expressions, Eqs. (16, 17), for the equatorial eccentric trajectories have their importance in various astrophysical studies, in addition to, gyroscope precession and phase space studies.

3 Non-equatorial separatrix trajectories

The separatrix orbits have been studied for the equatorial plane around a Kerr black hole [28, 20]. They have been shown as homoclinic orbits which asymptote to an energetically bound and unstable circular orbit. Here, we discuss the non-equatorial counterpart of these separatrix trajectories where these orbits asymptote to an energetically bound, unstable spherical orbit. These non-equatorial homoclinic trajectories are critical in calculating the evolution of test objects transiting from inspiral to plunge, which is not always confined to the equatorial plane, as in EMRIs emitting gravitational radiation.

In this section, we first deduce the expressions of $E$ and $L$ for the spherical orbits as functions of the radius $r_{s}$, and ($a$, $Q$). We then derive the exact expressions for the conic parameters ($e$, $\mu$) for non-equatorial separatrix orbits as a function of the radius of the corresponding spherical orbit, $r_{s}$, and ($a$, $Q$). We also show that these formulae reduce to the equatorial case, previously derived in [20], when $Q\rightarrow 0$ is applied. Next, we derive the exact analytic expressions for the non-equatorial separatrix trajectories by reducing it from the general trajectory formulae, Eqs. (12a-12c). We find that in this case, the radial part of the solutions can be reduced to a form that involves only logarithmic and trigonometric functions.

3.1 Energy and angular momentum of spherical orbits

Spherical orbits are the non-equatorial counterparts of circular orbits and set a crucial signpost in the dynamical study of non-equatorial and separatrix trajectories. The exact expressions for energy and angular momentum for the spherical orbits can be derived by substituting $e=0$ and $\mu=1/r_{s}$, where $r_{s}$ is the radius of the orbit, in the expressions for $E$, $L$, and $x$ given by Eqs. (7a-7e), which yields

[TABLE]

Similar formulae were derived in terms of inclination angle using an approximation in [29], whereas we have written the exact form in terms of the fundamental parameters and constant of motion $Q$. These formulae reduce to the energy and angular momentum formulae for circular orbits when $Q=0$ is substituted [8]:

[TABLE]

3.2 Exact expressions of conic variables for non-equatorial separatrix orbits

Similar to the case of equatorial plane, the non-equatorial separatrix trajectories can be parametrized by the radius of unstable spherical orbits, $r_{s}$, for a given combination of $a$ and $Q$, where $r_{s}$ varies from MBSO to ISSO. The energy and angular momentum of the separatrix orbits can be determined by Eqs. (18a-18c) by varying $r_{s}$ between the extrema MBSO and ISSO radii. In the ($e$, $\mu$) plane, these homoclinic orbits forms the boundary (other than $e=0$ and $e=1$ curves) of the allowed bound orbit region defined by Eq. (14) for a fixed $a$ and $Q$; see red curve in Fig. 2. The locus of this boundary in the ($e$, $\mu$) plane is obtained when equality is applied to the inequality Eq. (14a), which results in

[TABLE]

ISSO is a homoclinic orbit with $e=0$ and MBSO is a homoclinic orbit with $e=1$; hence the endpoints of the separatrix curve (red curve in Fig. 2) represents the ISSO and MBSO radii. At these endpoints, the parameter $\mu$ takes values as described below:

[TABLE]

The equations for ISSO and MBSO radii can be obtained using the equation of separatrix curve, Eq. (20), by plugging in ($e=0$, $\mu=1/r_{s}$) and ($e=1$, $\mu=1/2r_{s}$) to derive ISSO and MBSO respectively (as shown in D). Hence, the equations for these radii are given by

[TABLE]

for ISSO and

[TABLE]

for MBSO. The light radius for the spherical orbits can be obtained by equating the denominator of Eq. (18a) to zero, so that $E\rightarrow\infty$, which has the well known form for the equatorial light radius [8] given by

[TABLE]

Fig. 3 shows the contours of these radii in the ($r_{s}$, $a$) plane for various $Q$ values.

The effective potential diagram for the non-equatorial separatrix orbits shows double roots ($r_{2}=r_{3}$) of $R(r)$ at the periastron of the eccentric orbit and it also represents the spherical orbit radius, $r_{s}$ (see Fig 2). One of the remaining two roots of $R(r)$ represents the apastron ($=r_{1}>r_{s}$) of the eccentric orbit and the other inner root ($=r_{4}<r_{s}$) is not the part of bound trajectory.

Now, following a similar method used in [20], we derive the expressions for $e$ and $\mu$ for separatrix orbits with $Q\neq 0$. We write $R(r)=0$ in the form

[TABLE]

For the separatrix orbits, Eq. (25a) can be written as

[TABLE]

where $u_{s}=1/r_{s}$, $u_{1}=1/r_{1}$ apastron of the orbit, and $u_{4}=1/r_{4}$ corresponds to the inner most root of $R(r)$. The comparison of $u^{3}$ and constant term of the above equation with those of Eq. (25a) further gives the expression

[TABLE]

The conic parameters for such an orbit are given by

[TABLE]

where the substitution of $u_{1}$ and $u_{s}=1/r_{s}$ yields

[TABLE]

These expressions reduce to the ($e$, $\mu$) formulae for the equatorial separatrix orbits (see E for the details) when the limit $Q\rightarrow 0$ is implemented, to the forms previously derived by [20]:

[TABLE]

3.3 Exact forms for the non-equatorial separatrix trajectories

In this section, we show the reduction of our general trajectory solutions, Eqs. 12, for the case of separatrix orbits with $Q\neq 0$ to simple expressions. The separatrix or homoclinic orbits represent a curve in the ($e$, $\mu$) plane for a fixed $a$ and $Q$ combination, Fig. 2, which is also the boundary of the bound orbit region defined by Eqs. (14). This separatrix curve is defined by Eq. (20), which gives us the relation

[TABLE]

this further reduces the expressions of $A$, $B$, $C$ (Eqs. (9f-9h)) and correspondingly the expressions of $n^{2}$ and $m^{2}$ to

[TABLE]

The integrals governing the vertical motion ($\theta$ integrals) given by Eqs. (11a, 11b, 11c) retain their same form as they do not involve $k^{2}=1$, whereas, the radial integrals given by Eqs. (8a-8h) reduce further, when $k^{2}=1$ is substituted. The elliptic integrals reduce to forms involving trigonometric and logarithmic functions using the following identities given here

[TABLE]

The final and simple expressions for the azimuthal angle, $\left(\phi-\phi_{0}\right)$, $\left(t-t_{0}\right)$, and the equation relating $r-\theta$ motion for the non-equatorial separatrix trajectories (see F for the derivation) are given by

[TABLE]

where integrals $S_{3}-S_{7}$ are summarized in Table 3, and $H_{1}\left(\theta,\theta_{0},e,\mu,a,Q\right)$, $H_{2}\left(\theta,\theta_{0},e,\mu,a,Q\right)$ are given by Eq. (11a), (11b) respectively.

These expressions have their utility in evaluating the trajectory evolution of inspiralling objects near the separatrix, and just before plunging, for extreme mass ratio inspirals (EMRIs) in gravitational wave astronomy [30, 31, 32].

4 Trajectories

The analytic solution of the integrals of motion presented in this paper in §2.2 provides a direct and exact recipe to study bound trajectories without involving numerical integrations. These expressions have their utility in calculating extreme mass ratio inspirals (EMRIs) in gravitational wave astronomy, where numerical models consider an adiabatic progression through series of geodesics around a Kerr black hole [30, 31, 32]. We now discuss various kinds of bound geodesics around Kerr black hole using our analytic solution for the integrals of motion. We use the translation formulae, Eqs. (7a-7e), to obtain the integrals of motion only in terms of ($e$, $\mu$, $a$, $Q$) parameters. To sketch the trajectories, we have chosen the starting point for the trajectories to be ($\beta_{0}=\pi/2$, $\alpha=0$) as it follows from Eq. (4d). We use Eq. (12c) to calculate corresponding small change in $\theta$ or $\beta$ with the small change in $r$ or $\alpha$ and substitute corresponding ($r$, $\theta$) or ($\alpha$, $\beta$) values in Eqs. (12a) and (12b) to calculate ($\phi$, $t$).

There are various possible kinds of bound orbits. Here, we take up the each case and sketch these trajectories for different combinations of ($a$, $Q$), where the parameters values are tabulated in the Table 4. We take up slow rotating ($a=0.2$) and fast rotating black hole situations ($a=0.5$ or $a=0.8$), with both prograde and retrograde cases, for various $Q$ values. The various features of these orbits are enumerated below:

Eccentric orbits: Figs. 5 and 6 represent eccentric bound prograde and retrograde trajectories respectively, where the parameter values are depicted in the Table 4. The particle periodically oscillates between the periastron and the apastron, and is also bound between $\theta=\arccos\left(z_{-}\right)$ and $\theta=\arccos\left(-z_{-}\right)$ as shown in ($t$-$r$) and ($t$-$\theta$) plots in Figs. 5 and 6, whereas ($t$-$\phi$) plots depict that $\phi$ varies between 0 to $2\pi$. We have fixed ($e$, $\mu$) of the plotted trajectories and show the variation with change in $a$ and $Q$ parameters. The motion of the trajectory increases in the vertical direction with increase in $Q$ parameter.

Homoclinic/Separatrix orbits: Homoclinic orbits are the separatrices between eccentric bound and plunge orbits, where the particle asymptotically approaches the unstable spherical/circular orbit in both the distant past and the distant future. The energy and angular momentum of the orbiting particle simultaneously correspond to a stable eccentric bound orbit and an unstable spherical/circular orbit. Separatrix orbits in the equatorial plane of a Kerr black hole are well studied, [20, 28, 33]. The homoclinic orbits form an important group in Kerr dynamics as they represent the transition between inspiral and plunge orbits and hence, have their significance in the study of gravitational wave spectrum under the adiabatic approximation. The homoclinic or separatrix orbits correspond to the boundary of the region in ($e,\mu,a,Q$) space, defined by Eq. (20). Separatrix orbits with $Q\neq 0$ also have similar features as the equatorial separatrix orbits, where the particle asymptotically approaches the unstable spherical orbit. Figs. 7 and 8 show prograde and retrograde non-equatorial homoclinic/separatrix orbits respectively (see Table 4 for parameter values). We see from H3 and H4 trajectories that an increase in spin parameter, $a$, increases the range of $\theta$. The orbit initially follows an eccentric path and asymptotically approaches the periastron radius which also corresponds to the unstable spherical orbit radius as shown in ($t$-$r$) plots of Figs. 7 and 8. 3. 3.

Spherical orbits: Fig. 9 shows prograde and retrograde innermost stable spherical orbits (ISSO), which are also the homoclinic orbits with $e=0$. All the spherical stable orbits exist outside ISSO, whereas unstable spherical orbits are found between ISSO and MBSO. 4. 4.

Zoom-whirl orbits: Zoom whirl orbits are orbits where the particle takes a finite number of revolutions at the periastron before going back to the apastron, which is an extreme form of the periastron precession. Their significance in gravitational astronomy has been studied for the case of equatorial Kerr orbits [21]. Here, we discuss zoom-whirl orbits with $Q\neq 0$ as shown in Fig. 10, where the particle takes finite revolutions with varying $\theta$ at the periastron before turning back to the apastron. We have chosen the value of $\mu$ very near to the separatrix, where usually the zoom whirl behavior is seen, for different values of ($e$, $a$, $Q$) combinations. As expected, the particle spends more time at the periastron, compared to the time taken at apastron, to take a finite number of revolutions which is making the $t-r$ plots appear flatter near the periastron, (see Fig. 10). We again see that the increase in $a$ increases the range of vertical motion of the orbit like for the eccentric orbits case. Homoclinic/Separatrix orbit family is the limiting case of the zoom-whirl orbit family where the particle takes infinite revolutions as it asymptotes to the unstable spherical orbit.

Now, we discuss how different kinds of orbits are distributed in the bound orbit region in the ($e$, $\mu$) plane defined by the Eq. (14) for a fixed combination of ($a$, $Q$). We fix $a=0.5$ and $Q=5$ and show the shaded bound orbit region in Fig. 4, that represents the eccentric orbits allowed. The black curve which is the boundary of the shaded region represents homoclinic or separatrix orbits. The curve defined by $e=0$ represents all the spherical orbits with its end point at ISSO, which intersects with the separatrix line. We fix $e=0.5$ depicted by the red curve in Fig. 4 and take different values of $\mu$, as depicted by the black dots on the red curve, and plot the corresponding trajectories and study their corresponding behavior.

We see from Figs. 11 and 12, that for a fixed $e=0.5$, as $\mu$ is increased, the trajectory shows zoom-whirl behavior as it gets closer to the separatrix or homoclinic orbit for the corresponding $e$ value. It can be seen in the $t$-$r$ plot of Fig. 12 that the particle spends some time at the periastron which clearly depicts the zoom-whirl behavior. Hence, it can be said that zoom-whirl behavior is a near separatrix phenomenon and can occur at any eccentricity.

5 Fundamental frequencies

In this section, we derive the expressions for fundamental frequencies ($\nu_{\phi}$, $\nu_{r}$, $\nu_{\theta}$) in terms of the integrals derived analytically in §2.2. We take a long time average of Eq. (4d) on both the sides so that

[TABLE]

The similar expression can also be derived using formulae given in [16]. Again, we take a long time average of Eq. (4c), so that

[TABLE]

The limits of integral $I_{2}$ are $\alpha=\{0,\pi/2\}$, and that of $H_{2}$ are $\beta_{0}=\{\pi/2,-\pi/2\}$. The substitution of $H_{2}\left(-\frac{\pi}{2},\frac{\pi}{2},e,\mu,a,Q\right)$ and $\displaystyle{\frac{\nu_{\theta}}{\nu_{r}}}$ from Eqs. (11b) and (35c) in the above equations give

[TABLE]

Similarly, taking the long time average of Eq. (4b) and the substitution of $H_{1}$ and $H_{2}$ from Eqs. (11a) and (11b) yields

[TABLE]

where $I_{1}$, $I_{2}$, and $I_{8}$ are given by Eqs. (8a)-(8g) and (8h). Hence, the fundamental frequencies are explicit functions of input parameters ($e$, $\mu$, $a$, $Q$), which can be chosen using the bound orbit conditions presented in §2.3. These frequency formulae also match with the quadrature formulae derived in [11]; but here we have explicitly solved the integrals $I_{1}$, $I_{2}$ and $I_{8}$ in §2.2.

6 Consistency check with previous results

In order to verify our results, we have reduced our formulae for the non-equatorial separatrix trajectories, Eqs. (34c), to the case of equatorial separatrix orbits and found that they are consistent with earlier results derived in [20]. We also found that the frequency ratio, $\nu_{\phi}/\nu_{\theta}$, from Eqs. (37b, 37c), reduce to the case of maximally rotating black hole, $a=1$, for spherical orbits previously derived in [9]. See G for these derivations.

7 Applications

There are various important applications of our analytic solutions of the general non-equatorial trajectories and the fundamental frequencies for astrophysical studies as discussed below:

Gravitational waves: One of the crucial applications of our trajectory solution is the case of gravitational waves from the extreme-mass ratio inspirals (EMRIs). Our analytic formulae are directly applicable for the frequency domain calculation of the gravitational waves using the Teukolsky formalism, [34] or Kludge scheme [35], and the orbits can be computed more accurately than the numerical calculations [36]. Also, the homoclinic orbits, which are the separatrix between plunge and bound geodesics [22, 28], have their importance to study the zoom-whirl behavior of inspirals near separatrix [21, 37]. In this paper, we provide the analytic formulae for eccentricity and inverse-latus rectum, ($e$, $\mu$), for non-equatorial separatrix orbits which are crucial for the selection of these orbits for the study of gravitational waveforms in the Kerr geometry. 2. 2.

Relativistic precession: The exact analytic formula for azimuthal angle, $\phi-\phi_{0}$, is useful to find the precession of the orbits in the astrophysical systems like planets, black hole, and double pulsar systems. PSR J0737-3039 is one example of a double pulsar system having two pulsars, PSR J0737-3039A and PSR J0737-3039B having 23 ms [38] and 2.8 s [39] period respectively, which is useful to study the relativistic precession phenomenon valid in a strong gravitational field. The periastron advance was estimated in this source using the first PK parameter, $\dot{\omega}$ [40]. Our exact analytic results can be used to make a more accurate estimation of the relativistic advance of the periastron in pulsar systems where one component is having a major spin contribution. 3. 3.

Quasi-periodic oscillations (QPOs): QPOs are broad peaks seen in the Fourier power spectrum of the Neutron star X-ray binaries (NSXRB) and black hole X-ray binaries (BHXRB). The relativistic precession (RP) model was introduced [41] to explain the kHz QPOs in NSXRB and later applied to BHXRB [42]. The RP model can be used to calculate the black hole parameters assuming a circular or eccentric orbit is giving rise to a pair of observed high-frequency QPOs and a singular and nearly simultaneous corresponding Type-C QPO, [43, 44], where our exact formulae for the fundamental frequencies are applicable. 4. 4.

Gyroscope precession: The calculation of precession of spin of a test gyroscope is another application for the test of general relativity. In previous studies, approximate expressions were used for the fundamental frequencies as a series expansion in terms of eccentricity up to order $e^{2}$ around a Kerr black hole for the stable bound orbits in the equatorial plane [18]. Our exact analytic results are useful to estimate more accurate results which are useful to explain the reported results of geodetic drift rate and frame-dragging drift rate by the Gravity Probe B (GP-B) [45]. 5. 5.

Phase space study: Study of dynamics of Kerr orbits by Poincaré maps is also well discussed [22]. Our closed-form solutions are directly applicable to the study of the extreme chaotic behavior of orbits like Zoom-whirl orbits, which are extreme forms of perihelion precession [21], and their phase space structures.

8 Summary

The summary of this paper is given below:

We first translate the parameters ($E$, $L$, $a$, $Q$) to ($e$, $\mu$, $a$, $Q$) using the translation formulae, Eqs. (7) to completely describe the trajectory solution in the ($e$, $\mu$, $a$, $Q$) space. We then select the allowed bound orbit by choosing the parameters ($e$, $\mu$, $a$, $Q$) using the bound orbit conditions, Eqs. (14). 2. 2.

We have derived the closed-form analytic solutions of the general eccentric trajectory in the Kerr geometry as function of elliptic integrals, $\{\phi\left(r,\theta\right),t\left(r,\theta\right),r\left(\theta\right)\}$, Eqs. (12a12c). These trajectories around a Kerr black hole were previously derived in terms of Mino time [16], $\lambda$, subject to the initial conditions on $\mathrm{d}r\left(0\right)/\mathrm{d}\lambda$ and $\mathrm{d}\theta\left(0\right)/\mathrm{d}\lambda$. The application of our trajectory solution to the various possible studies is numerically faster and does not require any selection of initial conditions. We choose the starting point of the trajectory as the apastron of the orbit, $r_{a}$ or $\alpha=0$, and the initial polar angle, $\theta_{0}$ or $\beta_{0}$, is an extra parameter which can be arbitrarily chosen between maximum and minimum allowed $\theta$ range for a given $Q$. The input variables for plotting the trajectories are $\alpha$ and $\beta$ which define the range of $r$ and $\theta$ for a fixed combination of ($e$, $\mu$, $a$, $Q$). These results are summarized in Table 2. 3. 3.

We have derived the formulae for $E$ and $L$ for the spherical orbits as functions of radius $r_{s}$, $a$, and $Q$, given by Eqs. (18c). 4. 4.

We have derived the equations for ISSO, MBSO, and spherical light radius, Eqs. (22-24). The light radius derived is the same as that for the equatorial case. 5. 5.

We discussed the non-equatorial separatrix orbits, which asymptote to the unstable spherical radius sharing the same $E$ and $L$ values with the eccentric bound orbit. The radius of this unstable spherical radius for the separatrix orbit exists between $MBSO$ and $ISSO$. We write the exact forms for the eccentricity and inverse-latus rectum ($e_{s}$, $\mu_{s}$) for the non-equatorial separatrix orbits as functions of $r_{s}$, $a$, and $Q$, given by Eqs. (29d). 6. 6.

We use our general trajectory solutions to derive the equations of motions for non-equatorial separatrix orbits, given by Eqs. (34c), and find that the radial part of the solutions can be completely reduced to the form containing only trigonometric and logarithmic functions. We also show the reduction of these trajectories to the equatorial case which is also a new and useful form and match the solutions with the previously known result derived in [20]. Separatrix trajectories are essential in the study of gravitational waves from EMRIs, where our analytic solutions are directly applicable. These results are summarized in Table 3. 7. 7.

We discuss families of allowed bound orbit trajectories like non-equatorial eccentric, non-equatorial separatrix, zoom whirl, and spherical orbits around a rotating black hole using our analytic solution for the trajectories. Homoclinic trajectories have their applications in the gravitational wave astronomy as these trajectories are the boundaries between bound eccentric and plunge orbits. Separatrix/homoclinic orbits were studied for the equatorial case in [20, 28]. In this paper, we describe non-equatorial homoclinic and zoom-whirl trajectories, which is the more generalized case for the application to gravitational astronomy. 8. 8.

We derived the closed-form expressions for the fundamental frequencies in terms of elliptic integrals, Eqs. (37a, 37b, 37c), using the long time average method and without using Mino time, $\lambda$. We show that these expressions match with those derived in [11] using Hamilton-Jacobi formulation, which were left in the quadrature form, and we have obtained a closed form using elliptic integrals. These expressions are summarized in Table 5. We present the consistency of our trajectory solution by reducing it to the equatorial separatrix case and also show that the frequency ratio, $\nu_{\phi}/\nu_{\theta}$, matches with the standard expression derived [9] for the spherical orbits.

The results include novel aspects given in (i) and (iii)-(vii), listed above, and alternate new forms of the known formulae, given in the points (ii) and (viii) above. The equations and tables providing these results are indicated in the points above.

9 Discussion and Conclusions

There are several notable results in the vast literature discussing various aspects of dynamics in Kerr geometry such as the quadrature formulae for the trajectories ([6, 7]), circular orbit formulae [8], conditions for spherical orbits [9], expressions in terms of quadratures for the oscillation frequencies [11], formulae for trajectories in terms of quadratures for spherical polar motion [13], trajectories for non-spherical polar motion [12], and expressions for the trajectories and oscillation frequencies [16] in terms of Mino time [14]. Besides these key results there are other useful expressions reported for example on separatrix orbits [20], and on eccentric equatorial bound orbits ([6, 18]).

We discuss below the utility of the results in our paper:

The recipe for calculating frequencies and trajectories by [16] is as follows: The operative equations are $\{\phi(\lambda),t(\lambda),r(\lambda)\}$, Eqs. (6), (23)-(33), (35)-(45), which require linear combinations of many other equations. The analogy to $r(\chi)$ or $\chi(r)$ is $r(\lambda)$ or $\lambda(r)$ (Eqs. (26, 27)); the latter is non-trivial, whereas the former is simple. Given $\lambda$, ($\phi,t$) are calculated subsequently inverting linear combinations of many other elliptic integrals. We have numerically matched our frequency formulae with that given in [16] and we find that there is a minor typo in their expression of $\Gamma$ below Eq. (20) in section 3.3, where there is a factor of $E/2$ missing in the term $(r_{1}-r_{3})(r_{2}-r_{4})E(k_{r})$. However, the correct factor has been applied to calculate the numbers in their Tables (1, 2, 3) given in [16]. By using the set of equations in [16] and comparing with our expressions, it is found that our calculation is easier to implement and numerically faster by $\sim 20$, in the equatorial case, for example.

The novel results listed in (i) and (iii)-(vii) of the summary: translation conditions of $\{E,L\}\rightarrow\{e,\mu\}$, bound orbit conditions, $\{E(r_{s}),L(r_{s})\}$, ISSO, MBSO, and light radius formulae, besides new form for equatorial trajectories, are useful for various applications and simulations related to astrophysical scenarios involving relativistic precession like QPOs and accretion disks. We have also derived the locus of the $Q\neq 0$ separatrix curve in the $e-\mu$ plane besides providing the form of the trajectories. Using this, further studies can be carried for chaotic motion and study of gravitational waves from zoom-whirl orbits which can be set-up by locating them near the separatrix locus, in the same spirit, as was done for the equatorial case [20, 21].

The analytic results presented in this paper have direct applications in astrophysics for example, the study of non-equatorial separatrix orbits which has not been discussed before. They also help in understanding the highly eccentric behaviour of trajectories seen in numerical simulations [30] just before plunging onto the massive black hole in the case of EMRIs which is possibly related to the eccentric and inclined homoclinic orbits, besides relativistic precession in other astrophysical systems like binary pulsars and black holes, spin precession of gyroscopes around rotating black holes for the test of general relativity, and the study of chaotic orbits in the phase space.

Appendix A Deriving translation formulae between ($E$, $L$) and ($e$, $\mu$) parameters

The turning points of radial motion around rotating black hole are derived from Eq. (25a), which can be factorized into two quadratics

[TABLE]

and the comparison of coefficients on both sides yield the following relations

[TABLE]

Assuming that the first quadratic in Eq. (38) have the turning points of the orbit, $u_{1}=\mu\left(1-e\right)$ and $u_{2}=\mu\left(1+e\right)$, implies that ${a_{1}}^{{}^{\prime}}=-\left(u_{1}+u_{2}\right)=-2\mu$ and ${b_{1}}^{{}^{\prime}}=\mu^{2}\left(1-e^{2}\right)$. Using Eqs. (39d) and (39a) to replace ${a_{2}}^{{}^{\prime}}$ and ${b_{2}}^{{}^{\prime}}$ in Eqs. (39b) and (39c) and substituting for ${a_{1}}^{{}^{\prime}}=-\left(u_{1}+u_{2}\right)=-2\mu$ and ${b_{1}}^{{}^{\prime}}=\mu^{2}\left(1-e^{2}\right)$ yields

[TABLE]

The substitution of $a^{{}^{\prime}}$, $b^{{}^{\prime}}$, $c^{{}^{\prime}}$ and $d^{{}^{\prime}}$ and rearrangement of the terms in above equations gives

[TABLE]

where

[TABLE]

Next, the substitution of Eq. (41b) in Eq. (41a) gives

[TABLE]

which is further solved for $x^{2}$ to obtain Eqs. (7b-7e). These relations also reduce to the equatorial case, $Q=0$, which was first derived in [18].

Appendix B Solution of integrals $I_{1}$-$I_{8}$

In this appendix, we show the solutions of the integrals given in §2.2. First, we derive the radial integrals $I_{1}-I_{8}$, given by the Eqs. (8a-8h). We make the substitution $1/r^{{}^{\prime}}=\mu\left(1+e\cos\chi\right)$, which reduces the integrals to

[TABLE]

where $y=\left(1+e\cos\chi\right)$. Further, we implement the partial fraction method to reduce the integrals to

[TABLE]

[TABLE]

Next, we make the substitution, $\displaystyle{\cos\chi=2\cos^{2}\frac{\chi}{2}-1}$ and $\displaystyle{\psi=\frac{\chi}{2}-\frac{\pi}{2}}$, which reduces the integrals to

[TABLE]

where the constants, $C_{3}$-$C_{8}$, $n^{2}$, $m^{2}$, ${p_{1}}^{2}$, ${p_{2}}^{2}$, and ${p_{3}}^{2}$ are defined by Eq. (9) in §2.2. First, we solve the integrals $I_{6}$, $I_{3}$ or $I_{4}$ which are of the form given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Now, the substitution given by

[TABLE]

reduces the integrals in Eq. (48) to

[TABLE]

Hence, integrals given by Eqs. (8, 8, 8f) reduce to the forms given above. Next, we solve for $I_{5}$, which is of the form

[TABLE]

where the substitution given by Eq. (49) reduces the second integral in the above equation to

[TABLE]

where

[TABLE]

[TABLE]

the substitution given by Eq. (49) reduces the above expression to

[TABLE]

and multiplying and dividing the first integral in above equation by $\sqrt{1-k^{2}\sin^{2}\alpha}$ gives

[TABLE]

[TABLE]

which can be further reduced to

[TABLE]

where

[TABLE]

By defining integrals having a general form given by

[TABLE]

we use the following identity [26]:

[TABLE]

[TABLE]

The integral $I_{d}$ has a form similar to $T_{2}$ with $\left(h=1,g=\frac{p^{2}+m^{2}}{1-m^{2}},w=k\right)$, which yields

[TABLE]

where

[TABLE]

and

[TABLE]

were substituted. See Table 2 for the summary of the final expressions of $I_{1}-I_{7}$ derived using the method given in this section.

Appendix C Reduction to the equatorial plane ($Q=0$)

In this appendix, we reduce the integrals of motion to the case of equatorial plane. We start with the final expressions of $I_{3}$, $I_{4}$, $I_{5}$, and $I_{6}$ (in §2.2) and take the limit $Q\rightarrow 0$. As shown in §2.4, that for $Q\rightarrow 0$, we have $n^{2}\rightarrow 0$ which gives $\displaystyle{k^{2}=\frac{-m^{2}}{1-m^{2}}}$. Now, we first reduce the expressions of $I_{3}$ and $I_{4}$, Eqs. (8, 8) under the limit $Q\rightarrow 0$. We use the following identity (cf. [46], Eq. 160.02) to write

[TABLE]

where

[TABLE]

$\alpha_{2}^{2}=\alpha_{1}^{2}{k_{2}^{{}^{\prime}}}^{2}+k_{2}^{2}$, $k_{2}=\frac{k_{1}}{\sqrt{1+k_{1}^{2}}}$ and $k_{2}^{{}^{\prime}}=\frac{k_{2}}{k_{1}}$, which reduces $I_{3}$ and $I_{4}$ to the forms

[TABLE]

[TABLE]

The above expressions can be directly obtained if $n^{2}=0$ is substituted in the definition of $I_{3}$ and $I_{4}$ at the intermediate step, Eq. (47a). However, we aim to directly validate the final forms of ($I_{3}-I_{7}$). Hence, we apply the reduced expressions of $I_{3}$ and $I_{4}$ to $\left(\phi-\phi_{0}\right)$ for the equatorial plane to obtain Eq. (16a).

Next, we reduce the coordinate time integral. The expression for $I_{6}$ also reduces in similar way to $I_{3}$ and $I_{4}$ by applying the identity, Eq. (62), to

[TABLE]

which again can be directly obtained by substituting $n^{2}=0$ in the definition of $I_{6}$ in Eq. (47b). Next, we reduce the expression of $I_{5}$, Eq. (8e). We substitute for $\displaystyle{\sin\alpha=\frac{\sqrt{1-m^{2}}\sin\psi}{\sqrt{1-m^{2}\sin^{2}\psi}}}$, $\displaystyle{\cos\alpha=\frac{\cos\psi}{\sqrt{1-m^{2}\sin^{2}\psi}}}$, $\displaystyle{k^{2}=\frac{-m^{2}}{1-m^{2}}}$ and $\displaystyle{s^{2}=-\frac{p_{1}^{2}+m^{2}}{1-m^{2}}}$ into the expression of $I_{7}$, Eq. (8g), and use the following identities (cf. [46], Eq. 160.02)

[TABLE]

and Eq. (62) to transform $F\left(\alpha,k^{2}\right)$, $K\left(\alpha,k^{2}\right)$, and $\Pi\left(s^{2},\alpha,k^{2}\right)$ to

[TABLE]

The substitution of the above equations reduces $I_{7}$ to

[TABLE]

[TABLE]

Now, the substitution of $I_{7}$ from the above equation and $\sin\alpha$, $\cos\alpha$, $k^{2}$, $s^{2}$ into Eq. (8e) reduces $I_{5}$ to

[TABLE]

The above expression of $I_{5}$ can also be directly obtained by substituting $n^{2}=0$ in its definition given in Eq. (47b) and applying the identity given by Eq. (58). Hence, the substitution of $I_{3}$, $I_{4}$, $I_{5}$ and $I_{6}$ from Eqs. (64), (65), (73) and (66) in Eq. (12b) gives the expression for coordinate time, $t$, which simplifies to Eq. (16b).

See Table 2 for the final expressions of $(\phi,t)$ for the equatorial plane derived using the method given in this section.

Appendix D Innermost stable and marginally bound spherical radii

We discussed in §3.2 that ISSO and MBSO radii represent the end points of the separatrix curve defined by Eq. (20). Hence, we use this to write equations for ISSO and MBSO.

For the case of ISSO, we substitute $e=0$ and $\mu=1/r_{s}$ into Eq. (20), which gives

[TABLE]

which further expands, by the substitution of $x^{2}$ for spherical orbits, to

[TABLE]

The above equation has a complicated explicit expression of order twelve in $r_{s}$, which factorizes to

[TABLE]

where the second factor corresponds to the light radius [8] and the equation for ISSO is given by

[TABLE]

For the equatorial plane, $Q=0$ gives the following equation for the ISCO radius ($Z$)

[TABLE]

this equation can be factorized as

[TABLE]

where the first bracket gives the solution for light radius and the second bracket,

[TABLE]

gives the solution for ISCO [8] given by

[TABLE] 2. 2.

The condition for MBSO is derived by substituting $e=1$ and $\mu=1/(2Y)$ in Eq. (20), which yields

[TABLE]

where $r_{a}$ of the orbit reaches infinity and $r_{s}$ in the above equation corresponds to the unstable radius. By substituting $x^{2}$ for spherical orbits, the above equation reduces to

[TABLE]

The above expression expands to an equation of eleventh order in $r_{s}$, which factorizes to

[TABLE]

where the second factor corresponds to the light radius ($X$) [8] and the equation for MBSO is given by

[TABLE]

In the equatorial plane, $Q=0$, this reduces to

[TABLE]

which gets further factorized to

[TABLE]

where the last bracket in above equation gives the solution for light radius and the first two brackets give retrograde and prograde solutions for marginally bound orbit in the equatorial plane [8]. According to the sign convention used in this paper ($-1<a<1$), both retrograde and prograde cases are covered by the formula given by

[TABLE]

Appendix E Solution of ($e_{s}$, $\mu_{s}$) in the equatorial separatrix case

The expressions for eccentricity and inverse-latus rectum for the non-equatorial separatrix orbits, given by Eqs. (29a) and (29b) can be written in the form

[TABLE]

where the expressions for $a^{{}^{\prime}}$ and $d^{{}^{\prime}}$ are given by Eq. (25b). We solve for the factor in the denominator of $e_{s}$ by substituting $a^{{}^{\prime}}$ and $d^{{}^{\prime}}$ and taking the limit $Q\rightarrow 0$, we find

[TABLE]

[TABLE]

Substitution of Eq. (98) in (95, 96) gives

[TABLE]

Now, by substituting $x^{2}=\left(L-aE\right)^{2}$ and $E^{2}$ from Eq. (19) into the above equations, we find

[TABLE]

and

[TABLE]

which are the expressions in Eq. (30) given in §3.2.

Appendix F Reducing the radial integrals for the case of non-equatorial separatrix trajectories

Here, we reduce the radial integrals in our general trajectory solutions presented in §2.2, for the case of separatrix orbits with $Q\neq 0$. As shown in §3.3 that $k^{2}=1$ for the separatrix orbits, we use the identities given by Eq. (33) [26] to reduce the corresponding Elliptic integrals.

The substitution of these reduced form of the Elliptic integrals further reduces the integrals of the form $S_{3}$ and $S_{4}$, Eqs. (8, 8), to the general form

[TABLE]

where $p^{2}$ equals ${p_{2}}^{2}$ and ${p_{3}}^{2}$ for $S_{3}$ and $S_{4}$ respectively. Hence, from Eq. (8a) the expression of $I_{1}$ reduces in terms of $S_{3}$ and $S_{4}$ for the separatrix trajectories to

[TABLE]

Further, Eqs. (33) reduce the integrals $S_{5}-S_{8}$, Eqs. (8e-8h), to

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, the expression of $I_{2}$, Eq. (8b), reduces for the separatrix orbits in terms of $S_{3}-S_{6}$ to

[TABLE]

[TABLE]

See Table 3 for the summary of radial integrals for the separatrix trajectories.

Appendix G Derivations for the consistency check with the previous results

Equatorial separatrix orbits: Here, we show the reduction of the separatrix trajectory formulae to the equatorial case. We substitute $Q=0$ in the expression of the azimuthal angle given by Eq. (34a), which yields

[TABLE]

We saw that $n^{2}=0$ for the equatorial orbits and $m^{2}$ reduces to 1 for separatrix orbits, which gives $\alpha=0$ by definition, given by Eq. (49). Hence, we solve the integrals in the form of $\psi$ variable. First, we solve integrals $I_{3}$ and $I_{4}$ for separatrix orbits using

[TABLE]

Using the substitution $p_{2}\sin\psi=z$ for $I_{3}$ and $p_{3}\sin\psi=z$ for $I_{4}$, and applying the method of partial fractions, we find

[TABLE]

Hence, when the above equations for $I_{3}$ and $I_{4}$ are substituted into the expressions for azimuthal angle, Eq. (109), yields

[TABLE]

[TABLE]

Now, plugging in $Q=0$, reduces the expression for coordinate time of the separatrix orbits, Eq. (34b), to

[TABLE]

[TABLE]

And similarly, the integrals $I_{5}$ and $I_{6}$ reduce to

[TABLE]

Hence, the substitution of $I_{3}$, $I_{4}$, $I_{5}$, $I_{6}$ into Eq. (113) yields the expression for coordinate time for equatorial separatrix orbits to be

[TABLE]

[TABLE]

[TABLE]

We have numerically matched these expressions of the azimuthal angle and coordinate time for equatorial separatrix orbits, Eqs. (112) and (115), with those derived in [20]. 2. 2.

Spherical orbits:

Now, we verify our frequency formulae with the spherical orbit case. From Eqs. (37c) and (37b), we have

[TABLE]

In [9], the azimuthal to polar motion frequency ratio, $\nu_{\phi}/\nu_{\theta}$, for maximally rotating black hole, $a=1$, was found to be

[TABLE]

where $P=\left[E\left(r^{2}+1\right)-L\right]\left(r-1\right)^{2}-E$ for $a=1$.

For the case of spherical orbits, the limit $e\rightarrow 0$ reduces the ratio $\left(I_{1}/I_{8}\right)$ to the ratio of their integrands, which yields

[TABLE]

this reduces $\nu_{\phi}/\nu_{\theta}$ from, Eq. (116) to Eq. (117), and establishes the consistency of our frequency formulae with the spherical orbits case.