# Spherical particle orbits around a rotating black hole in massive gravity

Mohsen Fathi<EMAIL_ADDRESS><EMAIL_ADDRESS>Departamento de
Física, Universidad de Santiago de Chile, Avenida Víctor Jara 3493, Estación
Central, 9170124, Santiago, Chile  J.R. Villanueva<EMAIL_ADDRESS>Instituto de Física y Astronomía, Universidad de Valparaíso, Avenida Gran
Bretaña 1111, Valparaíso, Chile  Norman Cruz<EMAIL_ADDRESS>Departamento de Física, Universidad de Santiago de Chile, Avenida Víctor Jara
3493, Estación Central, 9170124, Santiago, Chile

###### Abstract

In this paper, we present a rotating de Rham-Gabadadze-Tolley black hole with
a positive cosmological constant in massive gravity, achieved by applying a
modified Newman-Janis algorithm. The black hole exhibits stable orbits of
constant radii, prompting a numerical study on the behavior of the solutions
to a nonic equation governing the radii of planar orbits around the black
hole. Additionally, we investigate the stability of orbits near the event
horizon and provide a comprehensive analytical examination of the solutions to
the angular equations of motion. This is followed by simulating some spherical
particle orbits around the black hole.

keywords: Particle geodesics, massive gravity, cosmological constant, ISCO

PACS numbers: 04.20.Fy, 04.20.Jb, 04.25.-g

###### Contents

  1. I Introduction and Motivation
  2. II Massive theory of gravity and its static black hole solution in the cosmological background
  3. III Spherical particle orbits
    1. III.1 Radii of planar orbits
      1. III.1.1 Non-monotonic behavior of the solutions
  4. IV Analytical solutions for the spherical particle orbits
    1. IV.1 The latitudinal motion
    2. IV.2 The azimuth motion
    3. IV.3 Explicit examples of non-planar orbits
  5. V Summary and conclusions
  6. A The roots of the polynomial equation $\Delta(r)=0$

##  I Introduction and Motivation

The significance of investigating stationary black hole spacetimes derived
from general relativity is undeniable. Among these, the Kerr-like black holes
have garnered significant attention, particularly following the groundbreaking
detections of gravitational waves from a merger by the Laser Interferometer
Gravitational-Wave Observatory (LIGO) PhysRevLett.116.061102  and the
remarkable imaging of black holes M87* and Sgr A* by the Event Horizon
Telescope (EHT) EventHorizonTelescope:2019ths ; EventHorizonTelescope:2022xnr
. Subsequently, there has been a surge in the study of geodesics involving
massless and massive particles in the vicinity of Kerr-like black holes, as
evidenced by recent publications (see for example, Refs.
rana_astrophysically_2019 ; kapec_particle_2020 ; gralla_null_2020 ;
stein_location_2020 ; compere_near-horizon_2020 ; rana_geometric_2020 ).
However, the study of particle motion around Kerr black holes dates back to
1968 with the introduction of Carter’s method for separating the Hamilton-
Jacobi equation Carter:1968 . This method allows orbits to be classified as
either planar or non-planar based on variations in their polar component.
While the resulting first-order differential equations of motion can be easily
solved numerically, finding exact analytical solutions for the coordinate
evolution is challenging due to their nonlinear nature. Extensive research has
been conducted to address this issue. For example, Mino demonstrated that by
introducing a new time parameter, two of the four equations of motion can be
decoupled Mino:2003 , leading to solutions expressed in terms of elliptic
integrals fujita_analytical_2009  (also discussed in the review article
nicolini_analytical_2016  and references therein). Moreover, when the radial
coordinate is held constant, the resulting orbits are either circular on the
equatorial plane or completely non-planar, forming spherical orbits. In fact,
the study of such orbits is valuable in astrophysical research as they help
determine the specific domains where light and particles are captured by the
black hole from different polar angles. Numerous studies have focused on
determining spherical particle orbits characterized by time-like constant-
radius geodesics, which follow the pioneering work by Wilkins PhysRevD.5.814 .
Since then, extensive investigations have explored various aspects of these
orbits that provide insights into their properties and behavior (see for
example, Refs. rana_astrophysically_2019 ; stein_location_2020 ; compere_near-
horizon_2020 ; PhysRevD.10.2324 ; stoghianidis_polar_1987 ;
hughes_evolution_2000 ; hughes_evolution_2001 ; kraniotis_precise_2004 ;
fayos_geometrical_2008 ; hackmann_analytical_2010 ; grossman_harmonic_2012 ;
hod_marginally_2013 ; teo_spherical_2021 ; Tavlayan:2021 ;
battista_geodesic_2022 ). In these studies, one often encounters algebraically
complex equations, such as the polynomial equations that govern the radii of
spherical orbits. Consequently, both rigorous analytical and numerical
approaches are necessary to determine reliable ranges for the motion
parameters that satisfy the orbital conditions in general relativistic
spacetimes, which are crucial for unveiling accurate insights into the
behavior and constraints of spherical orbits in the context of general
relativity.

However, while general relativity and its black hole solutions have been
successful in numerous observational tests, there are still unanswered
questions, particularly regarding the mysterious dark side of the universe
RevModPhys.75.1433 , which remains one of the major enigmas in modern
cosmology. In response, some scientists suggest that extending general
relativity in the right way could provide an alternative approach that may
eliminate the need for dark matter and dark energy, offering new perspectives
and potential solutions to these intriguing phenomena. Such extensions include
the $f(R)$ RevModPhys.82.451 ; de_felice_fr_2010  and scalar-tensor theories
of gravity quiros_selected_2019 ; fujii_maeda_2003 , in the four-dimensional
spacetimes. In particular, massive gravity has received significant attention,
tracing back to the Lorentz invariant Fierz-Pauli massive spin-2 theory in
1939 doi:10.1098/rspa.1939.0140 . Since then, this theory has been further
developed and generalized by ed Rham, Gabadadze, and Tolley (dRGT)
PhysRevD.82.044020 ; PhysRevLett.106.231101  (see also the review
de_rham_massive_2014 ). In Ref. ghosh_class_2016 , a spherically symmetric
black hole solution within the dRGT theory was presented. By generalizing the
Schwarzschild spacetime, this latter solution incorporates the new parameter
$\gamma$ that can potentially account for the flat galactic rotation curves,
in the sense that the massive gravitons constitute a dark matter halo
Panpanich:2018 . Furthermore, the aforementioned solution is endowed with a
positive cosmological constant, which serves to compensate for the accelerated
expansion of the universe. In fact, the presence of a positive cosmological
constant has been argued to have significant gravitational effects, not only
on black hole dynamics but also on the cosmos as a whole Ashtekar_2017 . In a
related study, it was demonstrated in Ref. Visser:2020  that the inclusion of
a positive cosmological constant enables black hole spacetimes to exhibit
innermost and outermost stable circular orbits (ISCOs and OSCOs), which have
important implications for the behavior and stability of particles around
black holes.

In line with the same research interest, our paper investigates a black hole
within the framework of dRGT massive gravity, specifically considering the
form proposed in Ref. Panpanich:2018 . This black hole solution incorporates a
positive cosmological constant, allowing for the formation of an ISCO and an
OSCO, as recently studied in Ref. universe7080278 . However, as highlighted in
the opening of this section, it is worth noting that stable circular orbits
represent only a subset of the broader category of spherical particle orbits,
which forms the primary focus of our investigation. In this regard, we
construct a rotating counterpart of the aforementioned dRGT spacetime and
explore various facets of spherical time-like geodesics within the black
hole’s exterior region. Consequently, this paper centers on two key
objectives; examining the radii of spherical orbits and their corresponding
orbit stability, and presenting analytical solutions to the angular equations
of motion for a comprehensive analysis. To achieve these objectives, the paper
is organized as follows: In Sect. II, we provide a concise introduction to the
dRGT massive gravity theory and its static black hole solution. Subsequently,
we present the rotating counterpart of this black hole spacetime and analyze
its causal structure. Moving on to Sect. III, we initiate the investigation of
spherical orbits using Carter’s method of separation of the Hamilton-Jacobi
equation. By applying the general criteria for the formation of spherical
orbits, we obtain a polynomial equation of order fourteen. However, for the
purpose of this investigation, we focus specifically on planar orbits and
numerically analyze the solutions to the corresponding characteristic
equation. In Section IV, we tackle the problem by solving the first-order
nonlinear differential equations governing the polar and azimuth angles. This
allows us to derive exact analytical solutions for the constant-radius
equations of motion that describe spherical orbits. These solutions are
expressed in terms of Weierstrassian elliptic functions, which are then
employed to showcase illustrative examples of orbits around the black hole.
Finally, in Sect. V, we summarize our findings. Throughout this study, we
adopt geometrized units where $G=c=1$, employ the sign convention $(-++\,+)$
for the spacetime line element, and denote differentiations with respect to
$r$ or $x$-coordinates using primes.

##  II Massive theory of gravity and its static black hole solution in the
cosmological background

Considering the Riemannian manifold $(\mathcal{M},g_{\mu\nu})$, the
gravitational action of the dRGT massive theory of gravity PhysRevD.82.044020
; PhysRevLett.106.231101 , can be re-expressed as Panpanich:2018

$\mathcal{S}=\frac{M_{\mathrm{Pl}}^{2}}{2}\int\mathrm{d}x^{4}\,\sqrt{-g}\Big{[}R+m_{g}^{2}~{}\mathcal{U}(\bm{g},\bm{f})\Big{]}+\mathcal{S}_{\mathrm{m}},$
(1)

in which $\mathcal{S}_{\mathrm{m}}$ is the matter action, $M_{\mathrm{Pl}}$
represents the reduced Planck mass, $m_{g}$ is the graviton’s mass, and
$\mathcal{U}$ is the gravitons’ potential, which to avoid the Boulware-Deser
ghost, must obey the following self-interaction

$\mathcal{U}=\mathcal{U}_{2}+\alpha_{3}\mathcal{U}_{3}+\alpha_{4}\mathcal{U}_{4},$
(2)

where

$\displaystyle\mathcal{U}_{2}=[\mathcal{K}]^{2}-[\mathcal{K}^{2}],$ (3a)
$\displaystyle\mathcal{U}_{3}=[\mathcal{K}]^{3}-3[\mathcal{K}][\mathcal{K}^{2}]+2[\mathcal{K}^{3}],$
(3b)
$\displaystyle\mathcal{U}_{4}=[\mathcal{K}]^{4}-6[\mathcal{K}]^{2}[\mathcal{K}^{2}]+3[\mathcal{K}^{2}]^{2}+8[\mathcal{K}][\mathcal{K}^{3}]-6[\mathcal{K}^{4}],$
(3c)

in which $[\mathcal{K}]=\mathcal{K}^{\mu}_{\mu}$ and
$(\mathcal{K}^{i})^{\mu}_{\nu}=\mathcal{K}^{\mu}_{\rho_{1}}\mathcal{K}^{\rho_{1}}_{\rho_{2}}\cdots\mathcal{K}^{\rho_{i}}_{\nu}$,
considering the definition

$\mathcal{K}^{\mu}_{\nu}=\delta^{\mu}_{\nu}-\sqrt{g^{\mu\lambda}\partial_{\lambda}\varphi^{a}\partial_{\nu}\varphi^{b}f_{ab}}\,.$
(4)

In the above expression, it is important to distinguish between the physical
metric $\bm{g}$ and the reference metric $\bm{f}$, that both act on the
Stückelberg field $\bm{\varphi}$. If the unitary gauge
$\bm{\varphi}=\bm{x}\cdot\bm{\delta}$ is taken into account, then one can
recast
$\sqrt{\bm{g}\cdot\bm{\partial\varphi}\cdot{\bm{\partial\varphi}}\cdot\bm{f}}=\sqrt{\bm{g}\cdot\bm{f}}$
in Eq. (4). Consequently, the field equations are obtained as

$G^{\mu}_{\nu}+m_{g}^{2}X^{\mu}_{\nu}=8\pi GT^{\mu(\mathrm{m})}_{\nu},$ (5)

in which $\bm{G}$ is the Einstein tensor, $\bm{T}^{(\mathrm{m})}$ is the
matter energy-momentum tensor, and ghosh_class_2016 ; PhysRevD.85.044024 ;
PhysRevD.87.064001

$X^{\mu}_{\nu}=\mathcal{K}^{\mu}_{\nu}-[\mathcal{K}]\delta^{\mu}_{\nu}-\alpha\left[(\mathcal{K}^{2})^{\mu}_{\nu}-[\mathcal{K}]\mathcal{K}^{\mu}_{\nu}+\frac{1}{2}\delta^{\mu}_{\nu}\left([\mathcal{K}]^{2}-[\mathcal{K}^{2}]\right)\right]\\\
+3\beta\left[(\mathcal{K}^{3})^{\mu}_{\nu}-[\mathcal{K}](\mathcal{K}^{2})^{\mu}_{\nu}+\frac{1}{2}\mathcal{K}^{\mu}_{\nu}\left([\mathcal{K}]^{2}-[\mathcal{K}^{2}]\right)-\frac{1}{6}\delta^{\mu}_{\nu}\left([\mathcal{K}]^{3}-3[\mathcal{K}][\mathcal{K}^{2}]+2[\mathcal{K}^{3}]\right)\right],$
(6)

is the massive graviton tensor. Now by assuming the ansatz
$f_{\mu\nu}=\mathrm{diag}\left(0,0,C^{2},C^{2}\sin^{2}\theta\right)$ with
$C>0$, the terms of order $\mathcal{O}(\mathcal{K}^{4})$ are eliminated, and
the static spherically symmetric vacuum solution to the dRGT massive gravity
will be characterized by the line element

$\mathrm{d}s^{2}=-B(r)\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{B(r)}+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2}\right),$
(7)

in the usual Schwarzschild coordinates, where the lapse function is given by
Panpanich:2018

$B(r)=1-\frac{2M}{r}-\frac{1}{3}\Lambda r^{2}+\gamma r+\zeta,$ (8)

in which $\Lambda$ performs the role of the cosmological constant, as in the
common Schwarzschild-de Sitter spacetime, and the parameters $\gamma$ and
$\zeta$ stem in the massive theory of gravity. The three terms
$\\{\Lambda,\gamma,\zeta\\}$ obey the relationships Panpanich:2018

$\displaystyle\Lambda=-3m_{g}^{2}(1+\alpha+\beta),$ (9a)
$\displaystyle\gamma=-m_{g}^{2}C(1+2\alpha+3\beta),$ (9b)
$\displaystyle\zeta=m_{g}^{2}C^{2}(\alpha+3\beta).$ (9c)

The flat space with $\zeta=0$ is obtained in terms of the conditions
$\alpha=-3\beta$ and $\beta=1/2+\epsilon$ with $0<\epsilon\ll 1$
Panpanich:2018 , which is respected in this study to guarantee the
positiveness of $\Lambda$ and $\gamma$.

Now, to obtain the rotating counterpart of this black hole spacetime, a
modified version of the Newman-Janis algorithm Newman:1965 , proposed by
Azreg-Aïnou Azreg:2014  is applied. This method employs a non-complexification
procedure to generate stationary spacetimes from their static counterparts. In
order to accomplish this, it is necessary to express the line element (7)
using the Eddington-Finkelstein coordinates $(\mathfrak{U},r,\theta,\phi)$.
This can be achieved by introducing the transformation
$\mathrm{d}\mathfrak{U}=\mathrm{d}t-\mathrm{d}r/B(r)$, which produces

$\mathrm{d}s^{2}=-B(r)\mathrm{d}\mathfrak{U}^{2}-2\mathrm{d}\mathfrak{U}\,\mathrm{d}r+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2}\right).$
(10)

By introducing the null tetrad set
$\mathfrak{Z}_{A}^{\mu}=\left(l^{\mu},n^{\mu},m^{\mu},\bar{m}^{\mu}\right)$,
in which $\bar{\bm{m}}$ is the complex conjugate of $\bm{m}$, and considering
$l^{\mu}=\delta_{r}^{\mu}$,
$n^{\mu}=\delta_{\mathfrak{U}}^{\mu}-\frac{1}{2}B(r)\delta_{r}^{\mu}$, and
$m^{\mu}=\big{(}\delta_{\theta}^{\mu}+\mathrm{i}\delta_{\phi}^{\mu}/\sin\theta\big{)}/(\sqrt{2}r)$,
one can express the inverse of the metric as
$g^{\mu\nu}=-l^{\mu}n^{\nu}-l^{\nu}n^{\mu}+m^{\mu}\bar{m}^{\nu}+m^{\nu}\bar{m}^{\mu}$.
Note that, the tetrad vector fields obey the conditions
$\bm{l}\cdot\bm{l}=\bm{n}\cdot\bm{n}=\bm{m}\cdot\bm{m}=\bm{l}\cdot\bm{m}=\bm{n}\cdot\bm{m}=0$,
and $\bm{l}\cdot\bm{n}=-\bm{m}\cdot\bar{\bm{m}}=-1$. Now to incorporate the
spin parameter $a$ of the stationary counterpart, we apply the transformation
$\delta_{\theta}^{\mu}\rightarrow\delta_{r}^{\mu}+\mathrm{i}a\sin\theta(\delta_{\mathfrak{U}}^{\mu}-\delta_{r}^{\mu})$,
while all the others remain unchanged Azreg:2014 . This procedure leads to the
transformations $B(r)\rightarrow\mathfrak{B}(r,a,\theta)$ and
$r^{2}\rightarrow\mathfrak{H}(r,a,\theta)$. Consequently, the null tetrad
assumes the form

$\displaystyle l^{\prime\mu}=\delta_{r}^{\mu},$ (11a) $\displaystyle
n^{\prime\mu}=\delta_{\mathfrak{U}}^{\mu}-\frac{\mathfrak{B}(r,a,\theta)}{2}\delta_{r}^{\mu},$
(11b) $\displaystyle
m^{\prime\mu}=\frac{1}{\sqrt{2\mathfrak{H}(r,a,\theta)}}\left[\mathrm{i}a\sin\theta\left(\delta_{\mathfrak{U}}^{\mu}-\delta_{r}^{\mu}\right)+\delta_{\theta}^{\mu}+\frac{\mathrm{i}}{\sin\theta}\delta_{\phi}^{\mu}\right],$
(11c)

and the new inverse metric becomes
$g^{\mu\nu}=-l^{\prime\mu}n^{\prime\nu}-l^{\prime\nu}n^{\prime\mu}+m^{\prime\mu}\bar{m}^{\prime\nu}+m^{\prime\nu}\bar{m}^{\prime\mu}$
in the Eddington-Finkelstein coordinates. To transition to the desired Boyer-
Lindquist coordinates, we adopt the transformations
$\mathrm{d}\mathfrak{U}=\mathrm{d}t^{\prime}+\omega(r)\mathrm{d}r$ and
$\mathrm{d}\phi=\mathrm{d}\phi^{\prime}+\chi(r)\mathrm{d}r$, and we choose
Azreg:2014

$\displaystyle\omega(r)=-\frac{r^{2}+a^{2}}{r^{2}B(r)+a^{2}},$ (12a)
$\displaystyle\chi(r)=-\frac{a}{r^{2}B(r)+a^{2}}.$ (12b)

The remaining steps involve setting
$\mathfrak{B}(r,a,\theta)=\mathfrak{H}(r,a,\theta)^{-2}[r^{2}B(r)+a^{2}\cos^{2}\theta]$
in order to eliminate the cross term $\mathrm{d}t\mathrm{d}r$ in the line
element. Additionally, we select
$\mathfrak{H}(r,a,\theta)=r^{2}+a^{2}\cos^{2}\theta$ to satisfy
$G_{r\theta}=0$. Finally, this algorithm generates the stationary spacetime of
the rotating dRGT (termed as RdRGT) black hole in the form

$\mathrm{d}s^{2}=-\frac{\Delta-a^{2}\sin^{2}\theta}{\Sigma}\mathrm{d}t^{2}+\frac{\Sigma}{\Delta}\mathrm{d}r^{2}-2a\sin^{2}\theta\left(1-\frac{\Delta-a^{2}\sin^{2}\theta}{\Sigma}\right)\mathrm{d}t\mathrm{d}\phi+\Sigma\mathrm{d}\theta^{2}\\\
+\sin^{2}\theta\left[\Sigma+a^{2}\sin^{2}\theta\left(2-\frac{\Delta-a^{2}\sin^{2}\theta}{\Sigma}\right)\right]\mathrm{d}\phi^{2},$
(13)

in which the spin parameter is given by $a=J/M$, with $J$ being the black
hole’s angular momentum. We have also defined

$\displaystyle\Delta(r)=a^{2}+r^{2}B(r)=\gamma
r^{3}+r^{2}+a^{2}-2Mr-\frac{\Lambda}{3}r^{4},$ (14a)
$\displaystyle\Sigma(r,\theta)=r^{2}+a^{2}\cos^{2}\theta.$ (14b)

The stationary black holes are surrounded by several hyper-surfaces that
configure their exterior spacetime. First of all, the horizons of the RdRGT
black hole are located by the roots of $\Delta(r)=0$, which results in the
quartic equation

$\Lambda r^{4}-3r^{2}\left(1+\gamma r\right)+6Mr-3a^{2}=0,$ (15)

admitting the discriminant

$\delta_{\Delta}=\frac{1}{27}\left(-256a^{6}\Lambda^{3}-729a^{4}\gamma^{4}-1296a^{4}\gamma^{2}\Lambda-384a^{4}\Lambda^{2}+1152a^{4}\gamma\Lambda^{2}M-108a^{2}\gamma^{2}\right.\\\
\left.-144a^{2}\Lambda+216a^{2}\gamma^{2}\Lambda
M^{2}+1728a^{2}\Lambda^{2}M^{2}-972a^{2}\gamma^{3}M-1440a^{2}\gamma\Lambda
M-1296\Lambda^{2}M^{4}\right.\\\ \left.+864\gamma^{3}M^{3}+1296\gamma\Lambda
M^{3}+108\gamma^{2}M^{2}+144\Lambda M^{2}\right).$ (16)

Since $\Lambda\ll 1$, the above discriminant can be positive within a certain
range for the values of $\gamma$, and hence, the above quartic can possess
four real roots, one negative and three positives. In Fig. 1, the region
inside which, $\delta_{\Delta}>0$, has been presented when the black hole’s
spin parameter varies from zero to one. Accordingly, for
$\gamma_{\min}<\gamma<\gamma_{\max}$, the positivity of the discriminant is
guaranteed for each value of $a$, and the cases of $\gamma=\gamma_{\min}$ and
$\gamma=\gamma_{\max}$ correspond to extremal black holes.

Figure 1: The region plot of $\delta_{\Delta}>0$ for $\Lambda=10^{-6}\,M^{-2}$
and $0<a<1$ is shown. According to the diagrams, for each value of $a$, the
positivity of the discriminant is guaranteed inside a particular domain
$\gamma_{\min}<\gamma<\gamma_{\max}$. In this case, the shaded region
represents the domain corresponding to $a=0.9M$. For all values of $a$, it
holds that $\gamma_{\min}<0$, and the width of the region $\delta_{\Delta}>0$
significantly increases as the spin parameter decreases. It is worth noting
that for the exact values of $\gamma_{\min}$ and $\gamma_{\max}$,
$\delta_{\Delta}=0$, and possessing these values results in extremal black
holes.

The figure shows that the width of this domain decreases with the increase of
the spin parameter. Therefore, the condition $\delta_{\Delta}=0$ is satisfied
for both fast and slow-rotating black holes if we choose $|\gamma|\ll 1$. Note
that, the extremal limits $\gamma_{\begin{subarray}{c}{\min}\\\
{\max}\end{subarray}}$ correspond to a vanishing discriminant, for which the
equation $\Delta(r)=0$ has one negative and three positive roots, which two of
them are degenerate (coinciding horizons). On the other hand, the equation
$\delta_{\Delta}=0$ has one positive, one negative, and two complex conjugate
roots, and hence, one can recast the discriminant in Eq. (16) as
$\delta_{\Delta}=27a^{4}(\gamma_{\max}-\gamma)(\gamma-\gamma_{\min})(\gamma^{2}-|\gamma_{1}|^{2})$,
in which $\gamma_{1}\in\mathbb{C}$. This quartic equation can be solved
explicitly for $\gamma_{\begin{subarray}{c}{\min}\\\ {\max}\end{subarray}}$,
but we find it unnecessary to present their analytical expressions. On the
other hand, the suppressed quartic (15) has the four analytical roots
$r_{1},r_{2},r_{3}$, and $r_{4}$, as they have been presented in appendix A.
It can be verified that for the aforementioned domain of $\gamma$, the
solutions (55)–(58) admit $r_{4}<0$ and $r_{1}>r_{2}>r_{3}>0$. This way, the
black hole will possess three horizons, which from smaller to larger are
namely, the Cauchy horizon $r_{-}=r_{3}$, the event horizon $r_{+}=r_{2}$, and
the cosmological horizon $r_{++}=r_{1}$. Once $\gamma=\gamma_{\max}$, the
extremal black hole is obtained as the Cauchy and event horizons coincide,
whereas for the case of $\gamma=\gamma_{\min}$, extremality corresponds to the
unification of the event and cosmological horizons. In Fig. 2, the behavior of
$\Delta(r)$ has been plotted for different values of the $\gamma$-parameter,
where the solutions to $\Delta=0$ and the extremal limit have been indicated.

(a)  (b)

Figure 2: The behavior of $\Delta(r)$ for a slow and a fast-rotating black
hole, plotted for different values of the $\gamma$-parameter. The thick and
thin solid curves represent, respectively, the negative and positive values,
and the dashed curves correspond to the extremal cases.

The other hypersurfaces that characterize the stationary spacetimes are those
that correspond to the static limits, and together, they form the black holes’
ergoregions. These regions, inside which no static observer can exist, are
identified by the equation $g_{tt}=0$, in accordance with the line element
(13). For the RdRGT black hole this results in another quartic, whose
solutions are exactly the same positive solutions
$r_{\mathrm{st}_{++}}=r_{1}$, $r_{\mathrm{st}_{+}}=r_{2}$ and
$r_{\mathrm{st}_{-}}=r_{3}$, as those given in Eqs. (55)–(58) with the same
form for the included parameters, considering the only replacement
$a^{2}\rightarrow a^{2}\cos^{2}\theta$ in Eq. (60c). Hence, it is
straightforward to see that for $\theta=0$, the static limits and the black
hole horizons coincide. Also as expected, these radii satisfy the conditions
$0<r_{\mathrm{st}_{-}}<r_{-}$ and
$r_{+}<r_{\mathrm{st}_{+}}<r_{\mathrm{st}_{++}}<r_{++}$, and comprise the
three boundaries of the interior and exterior ergoregions. In Fig. 3, the
radial profile of the function $g_{tt}$ has been plotted for the same
parameters as in Fig. 2.

(a)  (b)

Figure 3: The radial profile of $g_{tt}$ plotted for $\theta=\pi/4$ for the
two cases of $a=0.3M$ and $a=0.9M$, and for the same values of the
$\gamma$-parameter as in Fig. 2. The first two roots of $g_{tt}=0$ are shown
for all of the cases, and for the extremal black holes, only two static limits
are available.

To focus on the motion of particles in the exterior geometry of the RdRGT
black hole, we limit our analysis to the region defined by $r_{+}<r<r_{++}$.
Within this domain, we utilize Carter’s geodesic equations of motion to study
the behavior of particles in orbits with constant radii.

##  III Spherical particle orbits

Black holes with non-zero spin parameter can form several spherical particle
trajectories around them, each of which, stays on a sphere with
$r=\mathrm{const.}$, when the $\theta$-coordinate oscillates between two
turning points (for the case of Kerr black holes, see for example Refs.
Bardeen:1972a ; Bardeen:1973b ; Chandrasekhar:1998 ). These spherical orbits
exist in a certain interval and only a subclass of them, such as the stable
circular orbits that are confined by the ISCO and OSCO are on the equatorial
plane, while all the other orbits are non-planar. In this section, we employ
the standard geodesic equations for massive particles to investigate the
planar and non-planar $r$-constant particle orbits around the RdRGT black
hole.

Applying the Carter’s method of separation of the Hamilton-Jacobi equation
Carter:1968 , we can determine the time-like geodesic equations in the
exterior geometry of the RdRGT black hole, which are given as

$\displaystyle M\frac{\mathrm{d}r}{\mathrm{d}\lambda}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\mathcal{R}(r)},$ (17) $\displaystyle
M\frac{\mathrm{d}\theta}{\mathrm{d}\lambda}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\Theta(\theta)},$ (18) $\displaystyle
M\frac{\mathrm{d}\phi}{\mathrm{d}\lambda}$ $\displaystyle=$
$\displaystyle\frac{a}{\Delta}\left[E\left(r^{2}+a^{2}\right)-aL\right]-\left(aE-\frac{L}{\sin^{2}\theta}\right),$
(19) $\displaystyle M\frac{\mathrm{d}t}{\mathrm{d}\lambda}$ $\displaystyle=$
$\displaystyle\frac{r^{2}+a^{2}}{\Delta}\left[E\left(r^{2}+a^{2}\right)-aL\right]-a\left(aE\sin^{2}\theta-L\right),$
(20)

where we have used the dimensionless Mino time, $\lambda$, defined as
$\Sigma\mathrm{d}\lambda=M\mathrm{d}\tau$ Mino:2003 , in which $\tau$ is the
geodesics’ affine parameter. In Eqs.(17)–(20), $E$ and $L$ are the constants
of motion related to the Killing symmetries of the spacetime, which are
termed, respectively, as the energy and angular momentum of the particles.
Furthermore,

$\displaystyle\mathcal{R}(r)=\Big{[}E\left(r^{2}+a^{2}\right)-aL\Big{]}^{2}-\Delta\Big{[}\epsilon^{2}r^{2}+\left(L-aE\right)^{2}+\mathscr{Q}\Big{]},$
(21a)
$\displaystyle\Theta(\theta)=\mathscr{Q}-\cos^{2}\theta\left[a^{2}\left(\epsilon^{2}-E^{2}\right)+\frac{L^{2}}{\sin^{2}\theta}\right],$
(21b)

defining the Carter’s constant $\mathscr{Q}$, as the third constant of motion.
Since we are interested in the motion of massive particles, we set
$\epsilon=1$. Note that, from Eqs. (17) and (18), one can infer that the
particles’ motion is performed under the mutual conditions $\mathcal{R}(r)\geq
0$ and $\Theta(\theta)\geq 0$. Also for the sake of convenience, the positive
segments of these equations are adopted. It is important to accentuate the
role of $\mathscr{Q}$ in the classification of the orbits; orbits on the
equatorial plane (planar orbits that correspond to $\theta=\pi/2$) are
obtained for $\mathscr{Q}=0$. These orbits can be regarded as the boundaries
of the spherical orbits in their general case, for which, $\mathscr{Q}\geq 0$.
In fact, circular orbits satisfy the conditions
$\mathcal{R}(r)=0=\mathcal{R}^{\prime}(r)$, which provides the critical
quantities

$L_{c}(r)=\frac{1}{3a\Big{[}6M+r\left(4\Lambda r^{2}-9\gamma
r-6\right)\Big{]}}\Bigg{\\{}3E\bigg{[}a^{2}\bigg{(}6M+r^{2}(4\Lambda
r-9\gamma)\bigg{)}-6Mr^{2}+r^{4}(2\Lambda r-3\gamma)\bigg{]}\\\
-\bigg{[}\Lambda r^{4}-3\bigg{(}a^{2}+r\left(r-2M+\gamma
r^{2}\right)\bigg{)}\bigg{]}\sqrt{6r\bigg{[}6M+r\left(4\Lambda r^{2}-9\gamma
r-6\right)+6E^{2}r\bigg{]}}\,\Bigg{\\}},$ (22)
$\mathscr{Q}_{c}(r)=\frac{r^{2}}{3a^{3}\Big{[}6M+r\left(4\Lambda r^{2}-9\gamma
r-6\right)\Big{]}^{2}}\Bigg{\\{}ar\Bigg{[}27E^{2}r\bigg{[}20M^{2}-4Mr(3\gamma
r+4)+r^{2}\Big{(}\gamma r(5\gamma r+8)+4\Big{)}\bigg{]}\\\
+6\Lambda^{2}r^{6}\Big{[}18M-r(11\gamma r+10)+4E^{2}r\Big{]}+36\Lambda
r^{3}\bigg{[}E^{2}r\Big{(}2M-r(3\gamma r+2)\Big{)}+\Big{(}2M-r(\gamma
r+1)\Big{)}\\\ \times\Big{(}6M-r(5\gamma
r+4)\Big{)}\bigg{]}+54\left(r-2M+\gamma r^{2}\right)^{2}\Big{(}2M-r(3\gamma
r+2)\Big{)}+8\Lambda^{3}r^{9}\Bigg{]}\\\
+3a^{3}\Big{(}6M+r^{2}(3\gamma-2\Lambda r)\Big{)}\Big{(}6M+r\left(4\Lambda
r^{2}-9\gamma
r-6\right)+12E^{2}r\Big{)}+2\sqrt{6}E\Big{(}6M+r^{2}(3\gamma-2\Lambda
r)\Big{)}\\\ \times\bigg{[}3a^{3}+ar\Big{(}r\left(3-\Lambda r^{2}+3\gamma
r\right)-6M\Big{)}\bigg{]}\sqrt{r\Big{(}6M+r\left(4\Lambda r^{2}-9\gamma
r-6\right)+6E^{2}r\Big{)}}\Bigg{\\}}.$ (23)

Furthermore, the effective inclination angle can be defined by means of the
relation Tavlayan:2021

$\cos i=\frac{L}{\sqrt{L^{2}+\mathscr{Q}}}.$ (24)

Using this relation, and by eliminating $L$ between the equations
$\mathcal{R}(r)=0$ and $\mathcal{R}^{\prime}(r)=0$, we get to the single
equation

$p_{14}(x)=\sum_{j=0}^{14}m_{j}x^{j}=0,$ (25)

as the general characteristic equation for the radii of spherical particle
orbits in the RdRGT spacetime, where

$\displaystyle m_{0}=36k^{2}\nu^{2}u^{4},$ (26a) $\displaystyle
m_{1}=72\left(k-1\right)k\nu^{2}u^{4},$ (26b) $\displaystyle m_{2}=36\nu
u^{3}\left(-4k^{2}\nu+k^{2}\nu u+10k\nu-2k\nu u-2k+\nu u\right),$ (26c)
$\displaystyle m_{3}=72\nu u^{2}\left(2k^{2}\nu u-2k^{2}u-8k\nu+k\nu
u+3ku+4k-3\nu u-u\right),$ (26d) $\displaystyle
m_{4}=36u^{2}\left(-5hk\nu^{2}u-12k^{2}\nu^{2}+10k^{2}\nu+4k^{2}\nu^{2}u+12k\nu^{2}-28k\nu-8k\nu^{2}u+9\nu^{2}+14\nu+4\nu^{2}u+1\right),$
(26e) $\displaystyle m_{5}=-36u\left(5hk\nu^{2}u^{2}-16hk\nu^{2}u+8hk\nu
u-5h\nu^{2}u^{2}+12k^{2}\nu u-2kl\nu^{2}u^{2}-16k\nu-12k\nu^{2}u-24k\nu
u+24\nu\right.$ $\displaystyle\left.\qquad+12\nu^{2}u+12\nu u+8\right),$ (26f)
$\displaystyle m_{6}=12\left(45hk\nu u^{2}-45h\nu^{2}u^{2}-39h\nu
u^{2}+12k^{2}\nu^{2}u^{2}+6k^{2}\nu u^{2}+36k^{2}\nu
u+6kl\nu^{2}u^{3}-20kl\nu^{2}u^{2}+10kl\nu u^{2}\right.$
$\displaystyle\left.\qquad-24k\nu^{2}u^{2}-12k\nu u^{2}-138k\nu
u-30ku-6l\nu^{2}u^{3}+12\nu^{2}u^{2}+6\nu u^{2}+120\nu u+24u+48\right),$ (26g)
$\displaystyle m_{7}=-12\left(30hk\nu^{2}u^{2}+9hk\nu u^{2}+48hk\nu
u-30h\nu^{2}u^{2}-9h\nu u^{2}-132h\nu u-24hu+36k^{2}\nu u+12k^{2}u\right.$
$\displaystyle\left.\qquad+18kl\nu u^{2}-102k\nu
u-18ku-72k-18l\nu^{2}u^{2}-16l\nu u^{2}+66\nu u+6u+96\right),$ (26h)
$\displaystyle m_{8}=3\left(75h^{2}\nu^{2}u^{2}+372hk\nu u+72hku-540h\nu
u-36hu-384h+48k^{2}\nu u+108k^{2}+48kl\nu^{2}u^{2}+16kl\nu u^{2}\right.$
$\displaystyle\left.\qquad+80kl\nu u-96k\nu u-384k-48l\nu^{2}u^{2}-16l\nu
u^{2}-216l\nu u-40lu+48\nu u+288\right),$ (26i) $\displaystyle
m_{9}=-12\left(60h^{2}\nu u+33hk\nu u+72hk+15hl\nu^{2}u^{2}-33h\nu
u-132h+18k^{2}+38kl\nu u+8klu-42k-56l\nu u\right.$
$\displaystyle\left.\qquad-4lu-40l+24\right),$ (26j) $\displaystyle
m_{10}=6\left(45h^{2}\nu u+96h^{2}+102hk+98hl\nu u-120h+6k^{2}+28kl\nu
u+60kl-12k+6l^{2}\nu^{2}u^{2}-28l\nu u-112l\right.$
$\displaystyle\left.\qquad+6\right),$ (26k) $\displaystyle
m_{11}=-12\left(36h^{2}+9hk+19hl\nu u+40hl-9h+22kl+10l^{2}\nu u-26l\right),$
(26l) $\displaystyle m_{12}=81h^{2}+372hl+48kl+48l^{2}\nu u+100l^{2}-48l,$
(26m) $\displaystyle m_{13}=-8l(9h+10l),$ (26n) $\displaystyle
m_{14}=16l^{2},$ (26o)

in which we have used the dimensionless quantities

$\displaystyle x=\frac{r}{M},$ (27a) $\displaystyle u=\frac{a^{2}}{M^{2}},$
(27b) $\displaystyle k=E^{2},$ (27c) $\displaystyle l=\frac{1}{3}\Lambda
M^{2},$ (27d) $\displaystyle h=\gamma M,$ (27e)

and the definition $\nu=\sin^{2}i$, for which, the Carter’s constant takes the
form Tavlayan:2021

$\mathscr{Q}=\frac{\nu L^{2}}{1-\nu}.$ (28)

In general, the polynomial equation (25) has fourteen solutions of the form
$x(u,\nu,k,h,l)$, among which, those real roots that reside in the domain
$x_{+}<x<x_{++}$ are of our interest. On the other hand, and in accordance
with the Abel-Ruffini theorem, this equation cannot be solved analytically in
terms of finite radicals. We, therefore, consider some specific cases in what
follows in the rest of this section. Note that, for the case of Kerr black
hole with $h=l=0$, the above equation reduces to a polynomial equation of
order ten, whose analytical solutions to the radii of spherical time-like
orbits has been investigated in Ref. Tavlayan:2021 .

One can also make a primary classification of orbits, based on the value of
the $k$-parameter. In this regard, for $0\leq k<1$, the test particles travel
on bound orbits, and accordingly, the case of $k=1$ corresponds to marginally
bound orbits. On the other hand, $k>1$ provides unbound orbits. For unbound
spherical orbits the particles travel on unstable trajectories that are
deflected from the black hole upon small radial perturbations. This is while
for the bound orbits, the particle trajectories can still be stable or
unstable. In the latter case, the mentioned perturbations will generate
eccentric orbits that oscillate between two finite radial distances
teo_spherical_2021 .

To proceed with our discussion, we restrict ourselves to the case of planar
orbits at the vicinity of the black hole’s event horizon.

###  III.1 Radii of planar orbits

Planar orbits correspond to trajectories that occur in the equatorial plane,
for which $\mathscr{Q}=0$ (i.e., $i=0$, $\pi$, or $\nu=0$). These orbits are
of great importance in black hole astrophysics since they can also
characterize the orbit of particles within the accretion disk which are
confined within the boundaries composed by the ISCO and OSCO. The above
condition reduces the order of the polynomial equation (25) down to ten. A
further reduction relies on the smallness of $h$ and $l$, and their negligible
impacts at the vicinity of the black hole event horizon. In this sense, we can
take into account only up to the first order of these parameters for our
purpose. Hence, the characteristic polynomial equation reduces to the nonic
equation.

$p_{9}(x)=x^{9}+\sum_{j=0}^{8}\bar{m}_{j}x^{j}=0,$ (29)

in which

$\displaystyle\bar{m}_{0}=-\frac{u^{2}}{2hl},$ (30a)
$\displaystyle\bar{m}_{1}=\frac{4u}{hl},$ (30b)
$\displaystyle\bar{m}_{2}=\frac{(5k-4)u-8}{hl},$ (30c)
$\displaystyle\bar{m}_{3}=\frac{16-u\left(4h-2k^{2}+3k-1\right)-12k}{hl},$
(30d)
$\displaystyle\bar{m}_{4}=\frac{3\left(32h-9k^{2}+32k-24\right)+u(10l-18hk+9h)}{6hl},$
(30e)
$\displaystyle\bar{m}_{5}=\frac{\left(36hk-66h+9k^{2}-21k-20l+12\right)+2(2k-1)lu}{3hl},$
(30f) $\displaystyle\bar{m}_{6}=\frac{-51hk+60h-3k^{2}-30kl+6k+56l-3}{6hl},$
(30g) $\displaystyle\bar{m}_{7}=\frac{9hk+40hl-9h+22kl-26l}{6hl},$ (30h)
$\displaystyle\bar{m}_{8}=\frac{-31h-4k+4}{6h}.$ (30i)

To study the roots of this equation, $x(u,k,h,l)$, let us first fix $u=u_{0}$,
$h=h_{0}$, and $l=l_{0}$. This way, the solutions will be of the form
$x_{i}(k)$ $(i=\overline{1,9})$, for which, we have presented some examples in
the following.

(a)  (b)

Figure 4: The $x-k$ diagrams for the real parts of the solutions $x_{i}$ of
Eq. (29), plotted for $h=0.08$, $l=10^{-6}$, and $u=0.3$. In panel (a), the
whole range of $k$ has been shown, whereas the range at which $x_{7}$ and
$x_{8}$ are real, has been magnified in panel (b). The color-coding of the
solutions $x_{i}$ that is used here will be also applied in all of the
forthcoming diagrams within the paper.

In Fig. 4, the $k$-profiles of the solutions within the domains where
$x_{i}\in\mathbb{R}$, have been plotted for a slowly rotating black hole with
$u_{0}=0.3$, $h_{0}=0.08$ and $l_{0}=10^{-6}$, for which $x_{+}=1.607$. In
these diagrams, by $k_{i\rightarrow j}$ we notate the values of $k$, at which,
the solutions $x_{i}$ and $x_{j}$ $(i\leq j)$ coincide. As it is seen in the
left panel of the figure, $x_{3},x_{4}<x_{+}$ for $0<k<k_{5\rightarrow 3}$.
After this point, they are of real values, only for $x_{3},x_{4}>x_{+}$. Both
of these radii connect to $x_{5}$, respectively, after $k_{5\rightarrow
3}=1.838$ and $k_{5\rightarrow 4}=2.388$, and at $x=2.512$ and $x=4.934$.
Furthermore, $x_{1}$ and $x_{2}$ are real-valued only inside the event
horizon. $x_{5}$ has two significant branches in the domains $k_{6\rightarrow
5}=0.894<k<k_{5\rightarrow 3}$ and $k>k_{5\rightarrow 4}$, that respect the
condition $x_{5}>x_{+}$. At its upper limit, $x_{5}$ connects $x_{6}$ at
$x=6.318$. Note that, these radii correspond to unbound orbits. In fact, bound
orbits can happen inside the event horizon for the case of
$x_{5},x_{6}<x_{+}$, within the domain $0<k<k_{6\rightarrow 5}=0.144$. These
two branches coincide at $x=1.502$. The other possibility corresponds to the
solutions $x_{7}$ and $x_{8}$, which have been shown, separately, in the right
panel of Fig. 4. These radii are located outside the event horizon and are
real-valued in the domain $0<k<k_{8\rightarrow 7}=0.008$, for which, they
coincide at $x=2.45$. Since this latter case corresponds to bound orbits
outside the event horizon, it is worth discussing the stability of the orbits
in these radii, which will be referred to further in this subsection. But
before that, let us consider the $k$-profile of the solutions for a fast-
rotating black hole with $u=0.9$ and the same values for the other parameters,
for which, $x_{+}=1.603$. In Fig. 5, the regions where $x_{i}\in\mathbb{R}$
have been shown for this particular case.

Figure 5: The $x-k$ diagram for the real parts of the solutions $x_{i}$,
plotted for $h=0.08$, $l=10^{-6}$, and $u=0.9$.

The same as before, the real parts of $x_{8}$ and $x_{9}$ coincide. According
to the figure, $x_{3},x_{4}<x_{+}$ for $k_{4\rightarrow
3}=0.0003<k<k_{5\rightarrow 3}=1.615$. Beyond this domain, $x_{3}$ and $x_{4}$
switch, respectively, to $x_{1}$ and $x_{2}$ inside the event horizon, and
outside it, they have real values after $k_{5\rightarrow 3}$ and
$k_{5\rightarrow 4}=2.546$. The solution $x_{5}$ has three branches which are
ramified as follows. For $k_{4\rightarrow 3}<k<k_{6\rightarrow 5}=0.250$, it
is $x_{5}<x_{+}$, whereas for $k_{6\rightarrow 5}=0.731<k<k_{5\rightarrow 4}$,
it holds $x_{5}>x_{+}$. A particular domain for this case corresponds to
$0.731<k<k_{8}=0.901$, within which, bound orbits can form outside the event
horizon. In this domain, the lower limit of $x_{8}\in\mathbb{R}$ occurs at
$x=6.840$. Beyond this point, $x_{5}>x_{+}$ is real-valued only in the region
$k>k_{5\rightarrow 4}$. The solution $x_{6}$ is divided into two branches
outside the event horizon, which correspond to the domains $k_{4\rightarrow
3}<k<k_{6\rightarrow 5}=0.250$ and $k_{6\rightarrow 5}=0.731<k<k_{6}=1.137$,
within which, the upper limit of $x_{6}\in\mathbb{R}$ occurs at $x=3.40$.
Again, bound orbits can form inside the domain $0.731<k<k_{8}$. Note that, the
real part of the solution $x_{7}$ coincides with that of $x_{6}$.

To proceed further, it is also important to distinguish between the directions
of planar orbits outside the event horizon, within the energy domains
discussed above. For this reason, we construct the set of solutions
$x_{i}(u)\equiv x(u,k_{0},h_{0},l_{0})$, so that the $u$-profiles of the roots
of the nonic equation can be checked. Considering the previous values for
$h_{0}$ and $l_{0}$, the critical value of the spin parameter is
$u_{\mathrm{ext}}=0.923$, for which, $x_{+}=x_{-}=x_{\mathrm{ext}}=0.902$,
giving the extremal black hole horizon where the event and Cauchy horizons
coincide.

(a)  (b)

Figure 6: The $x-u$ diagrams for the real parts of the solutions $x_{i}$,
plotted for $h=0.08$, $l=10^{-6}$, (a) $k=0.003$ and (b) $k=3.2$.

We first consider $k_{0}=0.003$, which lies in the energy domain of Fig. 4(b),
and corresponds to the evolution of $x_{7}$ and $x_{8}$. The $u$-profiles of
the solutions have been shown in the left panel of Fig. 6. Like before, we now
let $u_{i\rightarrow j}$ to notate the values of the spin parameter, at which,
$x_{i}=x_{j}$. As expected from the $x-k$ diagrams, $x_{1}$ and $x_{2}$ are
totally complex-valued for this particular case, so they are absent from the
profiles. Furthermore, the real parts of $x_{8}$ and $x_{9}$ coincide. In this
diagram, we have also shown how the event and Cauchy horizons evolve until
they reach the extremal radius. The solutions $x_{3},x_{4},x_{5}$ and $x_{6}$
are inside the event horizon for this energy, but they are real-valued up to
$u=u_{5\rightarrow 4}=1.06>u_{\mathrm{ext}}$, so that will be available, as
well, around a naked singularity. At this point, the solutions $x_{4}$ and
$x_{5}$ meet at $x=0.574$. The two other solutions $x_{7}$ and $x_{8}$ exist
outside the event horizon, and by increasing the spin parameter up to
$u_{8\rightarrow 7}=0.463$, they finally coincide at $x=2.505$. Note that,
$x_{7}$ and $x_{8}$ correspond, respectively, to the retrograde and prograde
planar orbits for this choice of energy, and are confined in the range
$1.972<x<2.505$ within the domain $u_{7\rightarrow 6}=0.013<u<u_{8\rightarrow
7}=0.463$.

Let us now consider another example by choosing $k_{0}=3.2$, which corresponds
to the unbound orbits in Fig. 4(a). In this case,
$x_{6},x_{7},x_{8},x_{9}\in\mathbb{C}$, and the real-valued parts of $x_{4}$
and $x_{5}$ coincide. The solutions $x_{1},x_{2},x_{3}<x_{+}$ are extended
beyond $u_{\mathrm{ext}}$, up to $u_{3\rightarrow 2}=0.980$, at which
$x_{2}=x_{3}=0.984$. Hence they are available around the naked singularity.
The other two solutions $x_{3},x_{4}>x_{+}$ coincide at $x=3.367$ for
$u_{2\rightarrow 1}=0.004$, the spin parameter value at which, it is also
$x_{1}=x_{2}=0.001$. Note that, $x_{3}$ and $x_{4}$ correspond, respectively,
to prograde and retrograde planar orbits.

Now that we have established that bound orbits outside the event horizon can
be present on $x_{7}$ and $x_{8}$ for the above examples, we can now inspect
their stability. In fact, the stability condition for spherical particle
orbits respects the additional condition $R^{\prime\prime}(x)\geq 0$. From Eq.
(21a) and the changes of variables in Eqs. (27), one can infer that

$R^{\prime\prime}(x)=30lx^{4}-20hx^{3}+12\Big{[}klu-(1-k)\Big{]}x^{2}+6(2-hku)x-2u(1-k).$
(31)

Figure 7: The $u$-profile of $R^{\prime\prime}(x)$ plotted for $x_{7}$ and
$x_{8}$, in accordance with Fig. 6(a).

In Fig.7, the $u$-profile of the above function is plotted for the two radii
$x_{7}(u)$ and $x_{8}(u)$ in agreement with the parameter values considered in
Fig.6(a). As expected, the values of $R^{\prime\prime}(x)$ for these two
solutions are restricted to the range $u_{7\rightarrow 6}<u<u_{8\rightarrow
7}$, and in both cases $R^{\prime\prime}(x)>0$. Hence, the values of $x_{7}$
and $x_{8}$ correspond to stable circular orbits located outside the event
horizon. These orbits can, therefore, serve as generators of the ISCO for the
RdRGT black hole when $k=0.003$ and $u=0.3$. On the other hand, to obtain the
OSCO for this spin parameter, we can choose either $x_{5}$ or $x_{6}$, while
ensuring that the energy range satisfies $k_{6\rightarrow 5}=0.894\leq 1$.

####  III.1.1 Non-monotonic behavior of the solutions

The above examples provide numerical solutions for determining the radii of
prograde and retrograde circular orbits. Table 1 presents additional examples
obtained by solving the nonic equation (29) for the cases of $k=0,1$, and $2$,
in accordance with the parameter values shown in Figs. 4, 5, and 6. This table
offers a clearer understanding of the variations in these radii outside the
event horizon for the aforementioned cases. It should be noted that, under
certain circumstances, prograde or retrograde trajectories with nonzero
energies are only possible inside the event horizon, which is why they are not
included in the table.

$k$ | $u$ | $x_{\mathrm{prograde}}$ | $x_{\mathrm{retrograde}}$ | $k$ | $u$ | $x_{\mathrm{prograde}}$ | $x_{\mathrm{retrograde}}$ | $k$ | $u$ | $x_{\mathrm{prograde}}$ | $x_{\mathrm{retrograde}}$
---|---|---|---|---|---|---|---|---|---|---|---
0 | 0.1 | 1.762 | 2.213, 2.517 | 1 | 0.1 | – | 4.047 | 2 | 0.1 | 2.918 | –
0 | 0.2 | 1.652 | 2.295, 2.533 | 1 | 0.2 | – | 4.096 | 2 | 0.2 | 2.677 | –
0 | 0.3 | 1.562 | 2.357, 2.549 | 1 | 0.3 | – | 4.144 | 2 | 0.3 | 2.473 | –
0 | 0.4 | 1.480 | 2.408, 2.564 | 1 | 0.4 | – | 4.193 | 2 | 0.4 | 2.289 | –
0 | 0.5 | 1.403 | 2.452, 2.578 | 1 | 0.5 | – | 4.242 | 2 | 0.5 | 2.115 | –
0 | 0.6 | 1.328 | 2.491, 2.593 | 1 | 0.6 | – | 4.290 | 2 | 0.6 | 1.944 | –
0 | 0.7 | 1.252 | 2.526, 2.606 | 1 | 0.7 | – | 4.338 | 2 | 0.7 | 1.769 | –
0 | 0.8 | 1.175 | 2.559, 2.619 | 1 | 0.8 | 1.730 | 2.269, 4.386 | 2 | 0.8 | 1.578 | –
0 | 0.9 | 1.092 | 2.590, 2.632 | 1 | 0.9 | 1.335 | 2.437, 4.432 | 2 | 0.9 | 1.345 | –
0 | 1.0 | 1.0 | 2\. 618, 2.644 | 1 | 1.0 | – | 2.529, 4.478 | 2 | 1.0 | – | –
Table 1: The radii of prograde and retrograde orbits outside the event horizon
for different energies and spin parameters, obtained by assuming $h=0.08$ and
$l=10^{-6}$.

Here, we argue that these radii can also be calibrated by finding the critical
values of the spacetime parameters. To fulfill this task, let us calculate the
discriminant of the nonic equation (29) with respect to the $u$-parameter,
which is obtained as

$\delta_{\Delta_{9}}=h\bigg{[}-432k^{3}x^{4}+4k^{2}x\left(-72lx^{5}+216x^{3}-414x^{2}\right)-4k\left(-72lx^{6}+186lx^{5}+108x^{4}-414x^{3}+360x^{2}\right)\bigg{]}\\\
+144k^{4}x^{3}+48k^{3}x^{2}\left(4lx^{3}-9x+15\right)+4k^{2}x\left(16l^{2}x^{6}-96lx^{4}+180lx^{3}+108x^{2}-360x+288\right)\\\
-4k\bigg{[}16l^{2}x^{7}-40l^{2}x^{6}-48lx^{5}+180lx^{4}-6(26l-6)x^{3}-180x^{2}+288x-144\bigg{]}+4l^{2}(5-2x)^{2}x^{5},$
(32)

up to first order of $h$, that provides the critical value

$h_{\mathrm{crit}}=\frac{1}{6kx^{2}\Big{[}12(k-1)lx^{4}+18(k-1)^{2}x^{2}+69(k-1)x+31lx^{3}+60\Big{]}}\bigg{[}36k^{4}x^{3}+12k^{3}x^{2}\left(4lx^{3}-9x+15\right)\\\
+4k^{2}x\left(4l^{2}x^{6}-24lx^{4}+45lx^{3}+27x^{2}-90x+72\right)-4k\Big{(}4l^{2}x^{7}-10l^{2}x^{6}-12lx^{5}+45lx^{4}\\\
+(9-39l)x^{3}-45x^{2}+72x-36\Big{)}+l^{2}(5-2x)^{2}x^{5}\bigg{]}$ (33)

for which, $\delta_{\Delta_{9}}=0$. By means of this critical value, one can
determine the turning points of the polynomial (29). For the zero-energy
particles with $k=0$, we get

$h_{\mathrm{crit}}^{0}=\frac{2\Big{[}4l^{2}x^{7}-10l^{2}x^{6}-12lx^{5}+45lx^{4}+(9-39l)x^{3}-45x^{2}+72x-36\Big{]}}{3x^{2}\left(12lx^{4}-31lx^{3}-18x^{2}+69x-60\right)}.$
(34)

Furthermore, for $k=1$ and $k=2$ we obtain

$\displaystyle
h_{\mathrm{crit}}^{1}=\frac{l^{2}(2x+5)^{2}x^{5}+156lx^{3}+144}{6x^{2}\left(31lx^{3}+60\right)},$
(35) $\displaystyle
h_{\mathrm{crit}}^{2}=\frac{36l^{2}x^{7}+60l^{2}x^{6}+l(25l+96)x^{5}+360lx^{4}+24(13l+3)x^{3}+360x^{2}+576x+288}{12x^{2}\left(12lx^{4}+31lx^{3}+18x^{2}+69x+60\right)}.$
(36)

Now to find the desired radii, we first consider examples of $h_{0}$ and equal
them to $h_{\mathrm{crit}}^{i}\,\,(i=0,1,2)$. Using the Eqs. (34)–(36), the
values of $x_{0}$ can be obtained for each of the energy cases $k=0,1$ and
$2$. The nonic equation
$p_{9}(x_{0},u,\left\\{k_{0}=0,1,2\right\\},h_{0},l_{0})=0$, is then solved
for $u$, to find the corresponding $u_{0}$ for each of the cases. This spin
parameter is again substituted in the nonic equation
$p_{9}(x,u_{0},\left\\{k_{0}=0,1,2\right\\},h_{0},l_{0})=0$, to find the
$x$-solutions. For example, for $k_{0}=0$, $h_{0}=0.08$ and $l_{0}=10^{-6}$,
the prograde and retrograde radii are obtained as
$x_{\mathrm{prograde}}=1.542$ and $x_{\mathrm{retrograde}}=2.367$,
corresponding to $u_{0}=0.319$, which are in conformity with the data given in
Table 1.

In this section, we have examined the characteristic equations governing the
radii of planar orbits around the RdRGT black hole. These orbits correspond to
the case of $\nu=0$ in the main equation (25). For the investigation of polar
orbits, characterized by $\nu=1$ in Eq. (25), once again a polynomial equation
of fourteenth order is obtained. However, instead of directly analyzing this
equation, we can exploit the fact that for polar orbits, $L_{c}=0$. By
utilizing Eq. (22) and performing some algebraic manipulations, we arrive at
the octic equation

$4lx^{8}-\left(3h+14l\right)x^{7}+\left(11h+4k+6lu+10l-2\right)x^{6}\\\
-\left(5hu+8h+8k+6lu-10\right)x^{5}+\left(5hu+4ku+6k-8u-16\right)x^{4}\\\
+4\left(2-2ku+3u\right)x^{3}+2\left(ku^{2}+2ku-u^{2}-4u\right)x^{2}+2u^{2}x-2ku^{2}=0.$
(37)

Indeed, the methods employed to analyze the behavior of the solutions of the
aforementioned equation are quite similar to those discussed earlier for the
nonic equation (29). Consequently, we can close this section at this point and
proceed with our study by calculating the analytical solutions for the orbits.
With these solutions available, it becomes feasible to perform simulations of
spherical orbits around the RdRGT black hole.

However, before proceeding to the next section, it is crucial to highlight
some important features of stable circular orbits. In particular, solutions
such as $x_{7}$ and $x_{8}$ correspond to extremums in the radial effective
potential experienced by particles as they approach the black hole. These
radii give rise to circular orbits within the equatorial plane. When these
extremums represent minima in the effective potential (as in the cases of
$x_{7}$ and $x_{8}$), the stability of these orbits implies that deviations
from these minima result in unstable orbits. In such cases, particles exhibit
Keplerian motion, oscillating between two points known as the periapsis and
apoapsis of their trajectories. It is important to note that, unlike circular
orbits, these Keplerian trajectories are non-planar and correspond to
planetary bound orbits. However, when deviations from the potential minima
become significant and the turning points approach the maximum of the
effective potential, marginally bound orbits occur. At this point, particles
have the opportunity to escape the black hole. The maximum of the effective
potential exhibits a certain degree of stability, allowing for the existence
of spherical orbits. Nevertheless, these orbits are highly sensitive to radial
perturbations, which can cause the transition from bound to unbound motion.
Hence, particles can deflect and escape from the black hole when approaching
the maximum of the effective potential. In summary, the bound and marginally
bound circular orbits encompass two additional motion possibilities; particles
can either remain confined to the black hole or escape from its gravitational
field.

Now that we have discussed the classifications of particle motion in black
hole spacetimes, we can now continue with the derivation of the analytical
solutions for the equations of motion governing spherical orbits.

##  IV Analytical solutions for the spherical particle orbits

In this section, we derive the analytical solutions for the evolution
equations governing the polar and azimuth angles of spherical orbits, as
described in Section III. As discussed previously, the classification of
geodesics is based on the energy values, which determine whether the particle
motion is bound or unbound. The solutions encompass both planar and non-planar
orbits, corresponding to the cases of $\nu=0$ and $\nu\neq 0$ respectively.

###  IV.1 The latitudinal motion

Considering the definitions in Eqs. (27) and (28), one can re-write the
evolution equation (18) as

$-\frac{\mathrm{d}Z}{\mathrm{d}\lambda}=\sqrt{\Theta_{Z}},$ (38)

in which $Z=\cos\theta$, and for bound orbits of massive particles (i.e.
$k\leq 1$), we specify

$\Theta_{Z}=\left(1-k\right)uZ^{4}-\left[\left(1-k\right)u-\frac{L^{2}}{1-\nu}\right]Z^{2}+\frac{L^{2}\nu}{1-\nu}.$
(39)

Naturally, for planar orbits with $\nu=0$, the above expression reduces to
$\Theta_{Z}=(1-k)uZ^{4}-(L^{2}+u-uk)Z^{2}$, whereas for polar orbits with
$L=0$, it becomes $\Theta_{Z}=(1-k)uZ^{2}(Z^{2}-1)$. Note that, Eq. (39) can
help us determine the maximum and minimum latitudes reachable by particles,
which are the angular values where $\Theta_{Z}=0$. This equation provides the
values

$Z_{\begin{subarray}{c}{\max}\\\
{\min}\end{subarray}}^{2}=\frac{1}{2Z_{0}}\bigg{[}L^{2}+Z_{0}\pm\sqrt{\left(L^{2}+Z_{0}\right)^{2}-4L^{2}\nu
Z_{0}}\bigg{]},$ (40)

with $Z_{0}=(1-k)(1-\nu)u$, that confine the $Z$-parameter. In this way, the
polar angle $\theta$ oscillates in the domain
$\theta\in[\theta_{\min},\theta_{\max}]$, in which
$\theta_{\min}=\arccos(Z_{\max})$ and $\theta_{\max}=\arccos(Z_{\min})$. Note
that, there is a sign change for the azimuth angle $\phi$, where
${\mathrm{d}\phi}/{\mathrm{d}\lambda}=0$, and can be determined by means of
Eq. (19). After doing the substitutions, this angular turning point is
obtained as

$Z_{t}^{2}=\frac{\Big{(}\sqrt{ku}-L\Big{)}\Big{(}lx^{4}-hx^{3}+2x\Big{)}+Lx^{2}}{\sqrt{ku}\Big{(}lx^{4}-hx^{3}+2x\Big{)}-lu}.$
(41)

Hence, one can infer that the physically acceptable polar segments are where
$|Z_{t}|<|Z_{\max}|$. Now, to obtain the exact analytical solution to the
evolution of the polar angle, we directly integrate Eq. (38), which after
proper substitutions results in

$\theta(\lambda)=\arccos\left(Z_{\max}-\frac{3}{\psi_{2}-12\wp(\sqrt{\psi}\,\lambda)}\right),$
(42)

in which $\wp(\cdots)\equiv\wp(\cdots;\bar{g}_{2},\bar{g}_{3})$, is the
Weierstrassian $\wp$ elliptic function with the invariants

$\displaystyle\bar{g}_{2}=-\frac{1}{4}\left(\psi_{1}-\frac{\psi_{2}^{2}}{3}\right),$
(43a)
$\displaystyle\bar{g}_{3}=-\frac{1}{16}\left(\frac{2\psi_{2}^{3}}{27}-\frac{\psi_{1}\psi_{2}}{3}-\psi_{0}\right),$
(43b)

where

$\displaystyle\psi_{0}=\frac{(1-k)u}{\psi},$ (44a)
$\displaystyle\psi_{1}=\frac{4Z_{\max}(1-k)u}{\psi},$ (44b)
$\displaystyle\psi_{2}=-\frac{1}{\psi}\left[(1-k)u-\frac{L^{2}}{1-\nu}-6Z_{\max}^{2}(1-k)u\right],$
(44c)
$\displaystyle\psi=Z_{\max}\left[(1-k)uZ_{\max}^{2}+(1-k)u-\frac{L^{2}}{1-\nu}\right]-Z_{\max}\left[\frac{L^{2}}{1-\nu}-(1-k)u\right]+3(1-k)uZ_{\max}^{3}.$
(44d)

The above solution can also be used to determine the Mino time at which the
$\theta$-profile periodically traverses the equatorial plane and generates
nodes. Thus, for each full oscillation of the $\theta$-coordinate, the profile
encounters two nodes. According to Eq. (42) we get

$\lambda_{\mathrm{nod}}=\frac{1}{\sqrt{\psi}}\wp^{-1}\left(\frac{\psi_{2}}{12}-\frac{1}{4Z_{\max}}\right),$
(45)

for $\theta=\pi/2$, where
$\wp^{-1}(\cdots)\equiv\wp^{-1}(\cdots;\bar{g}_{2},\bar{g}_{3})$. Note that,
for unbound orbits with $k>1$, the form of the analytical solution (42)
remains the same, however, the Weierstrass invariants (43) are affected
through the exchange $(1-k)\rightarrow(k-1)$ in the quantities in Eqs. (44).

###  IV.2 The azimuth motion

The solution to the evolution equation (19) can be written as

$\phi(\lambda)=\mathcal{C}(x_{i})\Phi_{\theta_{1}}(\lambda)+\Phi_{\theta_{2}}(\lambda),$
(46)

for spherical particle orbits on the radii $x_{i}$, where

$\displaystyle{\mathcal{C}}(x_{i})=\frac{\sqrt{uk}}{\Delta(x_{i})}\Big{[}x_{i}^{2}+u-u\,\xi_{c}(x_{i})-\Delta(x_{i})\Big{]},$
(47a)
$\displaystyle\Phi_{\theta_{1}}(\lambda)=\int^{\theta(\lambda)}\frac{\mathrm{d}\theta}{\sqrt{\Theta(\theta)}},$
(47b)
$\displaystyle\Phi_{\theta_{2}}(\lambda)=\xi_{c}\int^{\theta(\lambda)}\frac{\mathrm{d}\theta}{\sin^{2}\theta\sqrt{\Theta(\theta)}},$
(47c)

in which $\xi_{c}(x)=L_{c}(x)/\sqrt{k}$, is the impact parameter of particles.
The integral (47b) can be directly calculated from the solution (42), which
provides

$\Phi_{\theta_{1}}(\lambda)=\frac{1}{\sqrt{\psi}}\wp^{-1}\left(\frac{\psi_{2}}{4}-\mathcal{U}_{\theta}\right),$
(48)

where

$\mathcal{U}_{\theta}=\frac{\psi_{2}}{6}+\frac{1}{4\left[Z_{\max}-\cos\theta\right]}.$
(49)

In addition, the integral (47c) can be solved by the same methods, yielding
the following:

$\Phi_{\theta_{2}}(\lambda)=\mathcal{K}_{0}\left[\mathcal{K}_{1}\mathscr{F}_{1}(\mathcal{U}_{\theta})-\mathcal{K}_{2}\mathscr{F}_{2}(\mathcal{U}_{\theta})-\sqrt{\psi}\lambda\right],$
(50)

where

$\displaystyle\mathcal{K}_{0}=\frac{\xi_{c}^{2}}{u(1-Z_{\max})(1+Z_{\max})\sqrt{2Z_{\max}\left(Z_{\max}^{2}+Z_{0}^{2}\right)}},$
(51a) $\displaystyle\mathcal{K}_{1}=\frac{1+Z_{\max}}{8(1-Z_{\max})},$ (51b)
$\displaystyle\mathcal{K}_{2}=\frac{1-Z_{\max}}{8(1+Z_{\max})},$ (51c)

and

$\mathscr{F}_{j}(\mathcal{U}_{\theta})=\frac{1}{\wp^{\prime}(\upsilon_{j})}\left[\ln\left(\frac{\sigma\left(\ss(\mathcal{U}_{\theta})-\upsilon_{j}\right)}{\sigma\left(\ss(\mathcal{U}_{\theta})+\upsilon_{j}\right)}\right)+2\ss(\mathcal{U}_{\theta})\zeta(\upsilon_{j})\right],\qquad
j=1,2,$ (52)

in which $\wp^{\prime}(\upsilon)=\mathrm{d}\wp(\upsilon)/\mathrm{d}\upsilon$.
Furthermore, $\sigma(\cdots)$ and $\zeta(\cdots)$ are the Sigma and Zeta
Weierstrassian elliptic functions, and the Weierstrass invariants are the same
as those given in Eq. (43). In Eq. (52), we have defined

$\displaystyle\upsilon_{1}=\wp^{-1}\left(\frac{\psi_{2}}{12}-\frac{1}{4|1-Z_{\max}|}\right),$
(53a)
$\displaystyle\upsilon_{2}=\wp^{-1}\left(\frac{\psi_{2}}{12}-\frac{1}{4|1+Z_{\max}|}\right).$
(53b)

In order to simulate the orbits, it is convenient choosing $\phi(Z_{\max})=0$.
Then by interpolating $\theta\rightarrow\theta(\lambda)$ in the expression
(48), the evolution of the $\phi$-coordinate is obtained explicitly.

###  IV.3 Explicit examples of non-planar orbits

We can now apply the analytical solutions obtained above, to make some
illustrative examples of spherical particle orbits around the RdRGT black
hole. To proceed with this, we first have to consider some specific initial
conditions for the black hole parameters. We then span the analytical
solutions (42) and (46) onto a Cartesian system defined as

$\displaystyle\mathcal{X}={x}\sin\theta\cos\phi,$ (54a)
$\displaystyle\mathcal{Y}={x}\sin\theta\sin\phi,$ (54b)
$\displaystyle\mathcal{Z}=x\cos\theta.$ (54c)

In Fig. 8, we have used the data given in Table 2 to simulate some exemplary
bound and unbound spherical orbits on slow and fast-rotating RdRGT black
holes. The radii of orbits have been obtained by numerically solving the
general equation (25), for different inclinations.

name | $u$ | $\nu$ | $h$ | $k$ | $x_{i}$ | $\xi_{c}(x_{i})$
---|---|---|---|---|---|---
(a) | 0.3 | 0.1 | 0.08 | 0.003 | 2.342 | 7.093
(b) | 0.9 | 0.1 | 0.08 | 0.890 | 1.921 | 2.208
(c) | 0.85 | 0.3 | 0.04 | 0.91 | 5.038 | 6.590
(d) | 0.5 | 0.4 | 0.004 | 1.1 | 2.533 | 3.938
(e) | 0.6 | 0.6 | 0.0004 | 0.93 | 5.803 | 7.916
(f) | 0.8 | 0.7 | 0.01 | 1.3 | 4.109 | 2.711
(g) | $u_{\mathrm{ext}}=0.990$ | 0.9 | 0.001 | 0.9 | 2.922 | 1.067
(h) | 0.85 | 1 (polar) | 0.0001 | 1 | 3.497 | 0.0045
Table 2: The information for the exemplary cases outside the event horizon,
considered for $l=10^{-6}$.

(a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)

Figure 8: Some examples of spherical particle orbits in accordance with the
data presented in Table 2. The sphere indicates the closure of points swept by
the radii $x_{i}$, which is cut to halves by a circle on the $\theta=\pi/2$
surface.

Referring to the table, case (a) is related to the solution $x_{7}$ in Fig.
4(b), but here corresponds to a retrograde bound orbit above the equatorial
plane. In the same sense, case (b) is related to the solution $x_{6}$, which
for the selected energy value corresponds to a bound orbit above the
equatorial plane. The other non-planar orbits are assumed as follows: Case (c)
considers a retrograde bound orbit for a fast-rotating black hole. In case
(d), a prograde unbound orbit is assumed. Case (e) once again assumes a
retrograde bound orbit, which is here derived for low energies. Case (f)
corresponds to a retrograde unbound orbit for a fast-rotating black hole. In
case (g), we present an example of prograde bound orbits around an extremal
black hole, and finally, in case (h), a polar bound orbit around a fast-
rotating black hole is considered. As expected, increasing the inclination
leads to orbits that span broader regions of the spherical surface,
accompanied by an increase in the frequency of polar angle oscillations. This
behavior persists until the orbits become entirely polar when $\nu=1$.
Additionally, we observe that for the same inclination, an increase in the
spin parameter reduces the difference between consecutive values of
$\lambda_{\mathrm{nod}}$, as the orbits more frequently intersect the
equatorial plane. Furthermore, we discover that a lower massive gravity
parameter results in more uniform orbits.

Given that the solutions are fully analytic, these aforementioned findings can
be replicated in different examples, allowing for simulations of various types
of orbits. Consequently, we leave the discussion at this point and proceed to
summarize our results in the following section.

##  V Summary and conclusions

In this paper, we presented a comprehensive analysis of the spherical orbits
of massive particles around an RdRGT black hole. Our investigation encompassed
both analytical and numerical approaches, with a focus on providing exact
analytical solutions whenever feasible. Commencing our investigation, we
provided a concise overview of the dRGT massive theory of gravity and its
static black hole solution. Employing a modified version of the Newman-Janis
algorithm, we proceeded to construct a rotating analog of the spacetime,
specifically termed as the RdRGT solution. We then examined the causal
structure of this spacetime by studying the fourth-order polynomial equation
$\Delta(r)=0$. We found that, under certain conditions, the equation’s
discriminant is positive, implying the presence of three horizons for the
black hole. By considering particular parameter cases, we were able to
determine specific ranges for the massive gravity parameter for the RdRGT
black hole with three horizons. We obtained the analytical expressions for the
radial positions of the horizons and complemented our findings with a
numerical investigation. This numerical study aimed to illustrate the radial
evolution of the black hole horizons and the static limit hypersurfaces for
various values of the massive gravity parameter. We continued our
investigation by presenting the nonlinear differential equations of motion for
massive particles and derived the conditions for the motion invariants, that
lead to spherical orbits. Subsequently, we derived the complete characteristic
polynomial equation governing the radii of spherical particle orbits in the
RdRGT spacetime. This equation is of the fourteenth order and encompasses all
possible inclinations. To focus our analysis, we specifically examined planar
orbits confined to the equatorial plane, resulting in a nonic equation.
Numerical analysis of the solutions to this nonic equation was conducted,
considering various aspects such as their evolution in terms of energy and
spin parameter for specific cases. We observed certain discontinuities in the
energy profiles; however, we found that the bound and unbound solutions remain
connected under specific criteria. By calculating the connecting energies, we
argued that there exist potentially bound and stable solutions beyond the
event horizon. Additionally, we conducted a detailed numerical examination of
the solutions’ evolution in relation to variations in the spin parameter,
distinguishing between prograde and retrograde particle orbits. Interestingly,
we discovered that the solutions remain valid even in the presence of naked
singularities. In all cases, we numerically calculated the intersection points
for the energy and spin parameter, accompanied by their corresponding radial
distances. Furthermore, through a comprehensive stability analysis, we
established the existence of at least two stable solutions outside the event
horizon, which could potentially serve as candidates for the ISCO around the
RdRGT black hole. We terminated our study on the radii of spherical orbits by
examining the critical values of the massive gravity parameter, at which, the
nonic equation becomes degenerate. We demonstrated that by utilizing these
critical values, we can obtain the solutions, effectively by regenerating the
desired spin parameters associated with each of the roots. The obtained
results were consistent with the previously determined radii of prograde and
retrograde orbits for each specific case. In the final section of the paper,
we focused on obtaining analytical solutions for the angular equations of
motion. This involved solving the evolution equations for both the polar and
azimuth angles in their most general forms. By directly integrating the
original differential equations, we expressed the solutions in terms of the
Weierstrassian elliptic functions, which covered both the bound and unbound
orbits at constant radii. To further explore their practical implications, we
considered specific cases involving definite values for the black hole
parameters and the initial energy of the test particles. Utilizing the derived
analytical solutions, we simulated the motion of the test particles on
spherical orbits. This approach includes a gradual increase in the inclination
while the overall massive gravity parameter decreases. We considered different
spin parameters and energies to encompass both bound and unbound orbits around
regular and extremal states of the RdRGT black hole. Our investigations
revealed a variety of behaviors, ranging from nearly planar orbits to polar
orbits. In this paper, we conducted numerical investigations on the solutions
of a nonic equation and performed analytical studies on the equations of
motion. Since the black hole spacetime of the RdRGT black hole extends that of
the Kerr black hole, the results and methods presented in this study can be
applied to different types of stationary spacetimes derived from other
extended theories of gravity that involve different matter properties.

For example, one intriguing aspect that can be studied is the redshift and
blueshift of particles trapped by the gravitational pull of the RdRGT black
hole, using a similar approach as in Ref. kraniotis_gravitational_2021 .
Furthermore, it is possible to investigate the shape of RdRGT black holes by
examining the distribution of accreting material along orbits of constant
radius. This can lead to predictions regarding the observed shapes of
astrophysical black holes, such as M87* and Sgr A* universe6090154 . The
methods used in this paper can also be applied to study the evolution of the
ISCOs in dynamical spacetimes Song:2021 . Additionally, the reflection
symmetry of the RdRGT black hole can be explored, which may provide insights
into the existence of curved accretion disks chen_curved_2022 . Moreover,
numerical analyses of the radii of planar orbits can be conducted for various
axially symmetric spacetimes, serving as tests for the analytical solutions
ospino_all_2022 . These methods can also be extended to investigate spherical
particle orbits and accretion in stationary black holes with quintessential
parameters or minimally coupled scalar fields PhysRevD.105.124039 ;
PhysRevD.106.024034 . In addition, the study of gravitomagnetic phenomena,
such as the Lense-Thirring effect, can be carried out for orbits of constant
radius in stationary spacetimes baines_constant-r_2022 . All of these
investigations can be further extended to explore regular rotating black holes
and solitons dymnikova_orbits_2023 . Recently, a new concept of particle
surfaces has been defined in Ref. 2023arXiv230612888B . These surfaces
represent the collection of points swept by particles on spherical orbits.
This notion can be generalized to other types of stationary spacetimes using
the methods provided in our paper and can be used to delimit the regions of
particles around black holes. These are just a few examples of potentially
interesting subjects that can be investigated in future works, utilizing the
numerical and analytical methods discussed in this paper.

## Acknowledgements

MF and NC acknowledge Universidad de Santiago de Chile for financial support
through the Proyecto POSTDOC-DICYT, Código 042331CM$\\_$Postdoc. JRV is
partially supported by the Centro de Astrofísica de Valparaíso (CAV).

##  Appendix A The roots of the polynomial equation $\Delta(r)=0$

By employing the method outlined in Appendix B of Ref. fathi_spherical_2023 ,
the solutions to the suppressed quartic equation (15) are derived as follows:

$\displaystyle
r_{1}=\frac{1}{\Lambda}\left[\tilde{\mathcal{A}}+\sqrt{\tilde{\mathcal{A}}^{2}-\tilde{\mathcal{B}}}\,\right]+\frac{3\gamma}{4\Lambda},$
(55) $\displaystyle
r_{2}=\frac{1}{\Lambda}\left[\tilde{\mathcal{A}}-\sqrt{\tilde{\mathcal{A}}^{2}-\tilde{\mathcal{B}}}\,\right]+\frac{3\gamma}{4\Lambda},$
(56) $\displaystyle
r_{3}=\frac{1}{\Lambda}\left[-\tilde{\mathcal{A}}+\sqrt{\tilde{\mathcal{A}}^{2}-\tilde{\mathcal{C}}}\,\right]+\frac{3\gamma}{4\Lambda},$
(57) $\displaystyle
r_{4}=\frac{1}{\Lambda}\left[-\tilde{\mathcal{A}}-\sqrt{\tilde{\mathcal{A}}^{2}-\tilde{\mathcal{C}}}\,\right]+\frac{3\gamma}{4\Lambda},$
(58)

in which

$\displaystyle\tilde{\mathcal{A}}=\sqrt{\tilde{\mathcal{U}}-\frac{\mathcal{A}}{6}},$
(59a)
$\displaystyle\tilde{\mathcal{B}}=2\tilde{\mathcal{A}}^{2}+\frac{\mathcal{A}}{2}+\frac{\mathcal{B}}{4\tilde{\mathcal{A}}},$
(59b)
$\displaystyle\tilde{\mathcal{C}}=2\tilde{\mathcal{A}}^{2}+\frac{\mathcal{A}}{2}-\frac{\mathcal{B}}{4\tilde{\mathcal{A}}},$
(59c)

where

$\displaystyle\mathcal{A}=-\left(\frac{27\gamma^{2}}{8\Lambda^{3}}+\frac{3}{\Lambda^{2}}\right)\Lambda^{3},$
(60a)
$\displaystyle\mathcal{B}=\left(\frac{6M}{\Lambda}-\frac{27\gamma^{3}}{8\Lambda^{3}}-\frac{9\gamma}{2\Lambda^{2}}\right)\Lambda^{3},$
(60b)
$\displaystyle\mathcal{C}=\left(\frac{2M\gamma}{2\Lambda}-\frac{27\gamma^{2}}{16\Lambda^{2}}-\frac{243\gamma^{4}}{256\Lambda^{3}}-3a^{2}\right)\Lambda^{3},$
(60c)

and

$\tilde{\mathcal{U}}=\sqrt{\frac{g_{2}}{3}}\cosh\left(\frac{1}{3}\operatorname{arccosh}\left(3g_{3}\sqrt{\frac{3}{g_{2}^{3}}}\right)\right),$
(61)

with

$\displaystyle g_{2}=\frac{\mathcal{A}^{2}}{48}+\frac{\mathcal{C}}{4},$ (62a)
$\displaystyle
g_{3}=\frac{\mathcal{A}^{3}}{864}-\frac{\mathcal{A}\mathcal{C}}{24}+\frac{\mathcal{B}^{2}}{64}.$
(62b)

## References

  * (1) B. P. Abbott et al., “Observation of gravitational waves from a binary black hole merger,” Phys. Rev. Lett., vol. 116, p. 061102, Feb 2016.
  * (2) K. Akiyama et al., “First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole,” Astrophys. J. Lett., vol. 875, no. 1, p. L4, 2019.
  * (3) K. Akiyama et al., “First Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way,” Astrophys. J. Lett., vol. 930, no. 2, p. L12, 2022.
  * (4) P. Rana and A. Mangalam, “Astrophysically relevant bound trajectories around a Kerr black hole,” Classical and Quantum Gravity, vol. 36, p. 045009, Feb. 2019.
  * (5) D. Kapec and A. Lupsasca, “Particle motion near high-spin black holes,” Classical and Quantum Gravity, vol. 37, p. 015006, Jan. 2020.
  * (6) S. E. Gralla and A. Lupsasca, “Null geodesics of the Kerr exterior,” Physical Review D, vol. 101, p. 044032, Feb. 2020.
  * (7) L. C. Stein and N. Warburton, “Location of the last stable orbit in Kerr spacetime,” Physical Review D, vol. 101, p. 064007, Mar. 2020.
  * (8) G. Compère and A. Druart, “Near-horizon geodesics of high spin black holes,” Physical Review D, vol. 101, p. 084042, Apr. 2020.
  * (9) P. Rana and A. Mangalam, “A Geometric Origin for Quasi-periodic Oscillations in Black Hole X-Ray Binaries,” The Astrophysical Journal, vol. 903, p. 121, Nov. 2020.
  * (10) B. Carter, “Global structure of the kerr family of gravitational fields,” Physical Review, vol. 174, pp. 1559–1571, Oct 1968.
  * (11) Y. Mino, “Perturbative approach to an orbital evolution around a supermassive black hole,” Phys. Rev. D, vol. 67, p. 084027, Apr 2003.
  * (12) R. Fujita and W. Hikida, “Analytical solutions of bound timelike geodesic orbits in Kerr spacetime,” Classical and Quantum Gravity, vol. 26, p. 135002, July 2009.
  * (13) C. Lämmerzahl and E. Hackmann, “Analytical Solutions for Geodesic Equation in Black Hole Spacetimes,” in 1st Karl Schwarzschild Meeting on Gravitational Physics (P. Nicolini, M. Kaminski, J. Mureika, and M. Bleicher, eds.), vol. 170, pp. 43–51, Cham: Springer International Publishing, 2016.
  * (14) D. C. Wilkins, “Bound geodesics in the kerr metric,” Phys. Rev. D, vol. 5, pp. 814–822, Feb 1972.
  * (15) M. Johnston and R. Ruffini, “Generalized wilkins effect and selected orbits in a kerr-newman geometry,” Phys. Rev. D, vol. 10, pp. 2324–2329, Oct 1974\.
  * (16) E. Stoghianidis and D. Tsoubelis, “Polar orbits in the Kerr space-time,” General Relativity and Gravitation, vol. 19, pp. 1235–1249, Dec. 1987.
  * (17) S. A. Hughes, “Evolution of circular, nonequatorial orbits of Kerr black holes due to gravitational-wave emission,” Physical Review D, vol. 61, p. 084004, Mar. 2000.
  * (18) S. A. Hughes, “Evolution of circular, nonequatorial orbits of Kerr black holes due to gravitational-wave emission. II. Inspiral trajectories and gravitational waveforms,” Physical Review D, vol. 64, p. 064004, Aug. 2001\.
  * (19) G. V. Kraniotis, “Precise relativistic orbits in Kerr and Kerr–(anti) de Sitter spacetimes,” Classical and Quantum Gravity, vol. 21, pp. 4743–4769, Oct. 2004.
  * (20) F. Fayos and C. Teijón, “Geometrical locus of massive test particle orbits in the space of physical parameters in Kerr space–time,” General Relativity and Gravitation, vol. 40, pp. 2433–2460, Nov. 2008.
  * (21) E. Hackmann, C. Lämmerzahl, V. Kagramanova, and J. Kunz, “Analytical solution of the geodesic equation in Kerr-(anti-) de Sitter space-times,” Physical Review D, vol. 81, p. 044020, Feb. 2010.
  * (22) R. Grossman, J. Levin, and G. Perez-Giz, “Harmonic structure of generic Kerr orbits,” Physical Review D, vol. 85, p. 023012, Jan. 2012.
  * (23) S. Hod, “Marginally bound (critical) geodesics of rapidly rotating black holes,” Physical Review D, vol. 88, p. 087502, Oct. 2013.
  * (24) E. Teo, “Spherical orbits around a Kerr black hole,” General Relativity and Gravitation, vol. 53, p. 10, Jan. 2021.
  * (25) A. Tavlayan and B. Tekin, “Radii of spherical timelike orbits around kerr black holes,” Phys. Rev. D, vol. 104, p. 124059, Dec 2021.
  * (26) E. Battista and G. Esposito, “Geodesic motion in Euclidean Schwarzschild geometry,” The European Physical Journal C, vol. 82, p. 1088, Dec. 2022\.
  * (27) W. L. Freedman and M. S. Turner, “Colloquium: Measuring and understanding the universe,” Rev. Mod. Phys., vol. 75, pp. 1433–1447, Nov 2003.
  * (28) T. P. Sotiriou and V. Faraoni, “$f(r)$ theories of gravity,” Rev. Mod. Phys., vol. 82, pp. 451–497, Mar 2010.
  * (29) A. De Felice and S. Tsujikawa, “f(R) Theories,” Living Reviews in Relativity, vol. 13, p. 3, Dec. 2010.
  * (30) I. Quiros, “Selected topics in scalar–tensor theories and beyond,” International Journal of Modern Physics D, vol. 28, p. 1930012, May 2019.
  * (31) Y. Fujii and K.-i. Maeda, The Scalar-Tensor Theory of Gravitation.  Cambridge Monographs on Mathematical Physics, Cambridge University Press, 2003.
  * (32) M. Fierz, W. E. Pauli, and P. A. M. Dirac, “On relativistic wave equations for particles of arbitrary spin in an electromagnetic field,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, vol. 173, no. 953, pp. 211–232, 1939.
  * (33) C. de Rham and G. Gabadadze, “Generalization of the fierz-pauli action,” Phys. Rev. D, vol. 82, p. 044020, Aug 2010.
  * (34) C. de Rham, G. Gabadadze, and A. J. Tolley, “Resummation of massive gravity,” Phys. Rev. Lett., vol. 106, p. 231101, Jun 2011.
  * (35) C. De Rham, “Massive Gravity,” Living Reviews in Relativity, vol. 17, p. 7, Dec. 2014.
  * (36) S. G. Ghosh, L. Tannukij, and P. Wongjun, “A class of black holes in dRGT massive gravity and their thermodynamical properties,” The European Physical Journal C, vol. 76, p. 119, Mar. 2016.
  * (37) S. Panpanich and P. Burikham, “Fitting rotation curves of galaxies by de rham-gabadadze-tolley massive gravity,” Phys. Rev. D, vol. 98, p. 064008, Sep 2018.
  * (38) A. Ashtekar, “Implications of a positive cosmological constant for general relativity,” Reports on Progress in Physics, vol. 80, p. 102901, aug 2017\.
  * (39) P. Boonserm, T. Ngampitipan, A. Simpson, and M. Visser, “Innermost and outermost stable circular orbits in the presence of a positive cosmological constant,” Phys. Rev. D, vol. 101, p. 024050, 2020.
  * (40) A. Rincón, G. Panotopoulos, I. Lopes, and N. Cruz, “Iscos and oscos in the presence of a positive cosmological constant in massive gravity,” Universe, vol. 7, no. 8, 2021.
  * (41) L. Berezhiani, G. Chkareuli, C. de Rham, G. Gabadadze, and A. J. Tolley, “On black holes in massive gravity,” Phys. Rev. D, vol. 85, p. 044024, Feb 2012\.
  * (42) Y.-F. Cai, D. A. Easson, C. Gao, and E. N. Saridakis, “Charged black holes in nonlinear massive gravity,” Phys. Rev. D, vol. 87, p. 064001, Mar 2013\.
  * (43) E. T. Newman and A. I. Janis, “Note on the kerr spinning‐particle metric,” Journal of Mathematical Physics, vol. 6, no. 6, pp. 915–917, 1965.
  * (44) M. Azreg-Aïnou, “Generating rotating regular black hole solutions without complexification,” Phys. Rev. D, vol. 90, p. 064041, Sep 2014.
  * (45) J. M. Bardeen, W. H. Press, and S. A. Teukolsky, “Rotating Black Holes: Locally Nonrotating Frames, Energy Extraction, and Scalar Synchrotron Radiation,” Astrophysical Journal, vol. 178, pp. 347–370, Dec. 1972.
  * (46) J. Bardeen, “Timelike and null geodesics in the Kerr metric,” in Les Houches Summer School of Theoretical Physics: Black Holes, pp. 215–240, 1973\.
  * (47) S. Chandrasekhar, The mathematical theory of black holes.  Oxford classic texts in the physical sciences, Oxford University Press, 1998.
  * (48) G. V. Kraniotis, “Gravitational redshift/blueshift of light emitted by geodesic test particles, frame-dragging and pericentre-shift effects, in the Kerr–Newman–de Sitter and Kerr–Newman black hole geometries,” The European Physical Journal C, vol. 81, p. 147, Feb. 2021\.
  * (49) V. I. Dokuchaev and N. O. Nazarova, “Visible shapes of black holes m87* and sgra*,” Universe, vol. 6, no. 9, 2020.
  * (50) Y. Song, “The evolutions of the innermost stable circular orbits in dynamical spacetimes,” European Physical Journal C, vol. 81, p. 875, Oct. 2021.
  * (51) C.-Y. Chen and H.-Y. K. Yang, “Curved accretion disks around rotating black holes without reflection symmetry,” The European Physical Journal C, vol. 82, p. 307, Apr. 2022.
  * (52) J. Ospino, J. L. Hernández-Pastora, and L. A. Núñez, “All analytic solutions for geodesic motion in axially symmetric space-times,” The European Physical Journal C, vol. 82, p. 591, July 2022.
  * (53) W. Cao, W. Liu, and X. Wu, “Integrability of kerr-newman spacetime with cloud strings, quintessence and electromagnetic field,” Phys. Rev. D, vol. 105, p. 124039, Jun 2022.
  * (54) I. Bogush, D. Gal’tsov, G. Gyulchev, K. Kobialko, P. Nedkova, and T. Vetsov, “Photon surfaces, shadows, and accretion disks in gravity with minimally coupled scalar field,” Phys. Rev. D, vol. 106, p. 024034, Jul 2022.
  * (55) J. Baines, T. Berry, A. Simpson, and M. Visser, “Constant-r geodesics in the Painlevé–Gullstrand form of Lense–Thirring spacetime,” General Relativity and Gravitation, vol. 54, p. 79, Aug. 2022.
  * (56) I. Dymnikova and A. Dobosz, “Orbits of Particles and Photons around Regular Rotating Black Holes and Solitons,” Symmetry, vol. 15, p. 273, Jan. 2023.
  * (57) I. Bogush, K. Kobialko, and D. Gal’tsov, “Glued massive particles surfaces,” arXiv e-prints, p. arXiv:2306.12888, June 2023.
  * (58) M. Fathi, M. Olivares, and J. R. Villanueva, “Spherical photon orbits around a rotating black hole with quintessence and cloud of strings,” The European Physical Journal Plus, vol. 138, p. 7, Jan. 2023.