Spherical particle orbits around a rotating black hole in massive gravity

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.


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) [1] and the remarkable imaging of black holes M87* and Sgr A* by the Event Horizon Telescope (EHT) [2,3].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.[4][5][6][7][8][9]).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 [10].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 [11], leading to solutions expressed in terms of elliptic integrals [12] (also discussed in the review article [13] and references therein).Moreover, when the radial coordinate is held constant, the resulting orbits are either circular on the equatorial plane or completely nonplanar, 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 [14].Since then, extensive investigations have explored various aspects of these orbits that provide insights into their properties and behavior (see for example, Refs.[4,7,8,[15][16][17][18][19][20][21][22][23][24][25][26]).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 [27], 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) [28,29] and scalar-tensor theories of gravity [30,31], 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 [32].Since then, this theory has been further developed and generalized by ed Rham, Gabadadze, and Tolley (dRGT) [33,34] (see also the review [35]).In Ref. [36], a spherically symmetric black hole solution within the dRGT theory was presented.By generalizing the Schwarzschild spacetime, this latter solution incorporates the new parameter γ that can potentially account for the flat galactic rotation curves, in the sense that the massive gravitons constitute a dark matter halo [37].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 [38].In a related study, it was demonstrated in Ref. [39] 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. [37].This black hole solution incorporates a positive cosmological constant, allowing for the formation of an ISCO and an OSCO, as recently studied in Ref. [40].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 (M, g µν ), the gravitational action of the dRGT massive theory of gravity [33,34], can be re-expressed as [37] in which S m is the matter action, M Pl represents the reduced Planck mass, m g is the graviton's mass, and U is the gravitons' potential, which to avoid the Boulware-Deser ghost, must obey the following self-interaction where in which In the above expression, it is important to distinguish between the physical metric g and the reference metric f , that both act on the Stückelberg field ϕ.If the unitary gauge ϕ = x • δ is taken into account, then one can recast 4).Consequently, the field equations are obtained as in which G is the Einstein tensor, T (m) is the matter energy-momentum tensor, and [36,41,42] is the massive graviton tensor.Now by assuming the ansatz f µν = diag 0, 0, C 2 , C 2 sin 2 θ with C > 0, the terms of order O(K 4 ) are eliminated, and the static spherically symmetric vacuum solution to the dRGT massive gravity will be characterized by the line element in the usual Schwarzschild coordinates, where the lapse function is given by [37] in which Λ performs the role of the cosmological constant, as in the common Schwarzschild-de Sitter spacetime, and the parameters γ and ζ stem in the massive theory of gravity.The three terms {Λ, γ, ζ} obey the relationships [37] The flat space with ζ = 0 is obtained in terms of the conditions α = −3β and β = 1/2 + with 0 < 1 [37], which is respected in this study to guarantee the positiveness of Λ and γ.Now, to obtain the rotating counterpart of this black hole spacetime, a modified version of the Newman-Janis algorithm [43], proposed by Azreg-Aïnou [44] 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 (U, r, θ, φ).This can be achieved by introducing the transformation dU = dt − dr/B(r), which produces By introducing the null tetrad set Z µ A = (l µ , n µ , m µ , mµ ), in which m is the complex conjugate of m, and considering l µ = δ µ r , n µ = δ µ U − 1 2 B(r)δ µ r , and m µ = δ µ θ + iδ µ φ / sin θ /( √ 2r), one can express the inverse of the metric as g µν = −l µ n ν − l ν n µ + m µ mν + m ν mµ .Note that, the tetrad vector fields obey the conditions l Now to incorporate the spin parameter a of the stationary counterpart, we apply the transformation δ µ θ → δ µ r + ia sin θ(δ µ U − δ µ r ), while all the others remain unchanged [44].This procedure leads to the transformations B(r) → B(r, a, θ) and r 2 → H(r, a, θ).Consequently, the null tetrad assumes the form and the new inverse metric becomes g µν = −l µ n ν − l ν n µ + m µ m ν + m ν m µ in the Eddington-Finkelstein coordinates.To transition to the desired Boyer-Lindquist coordinates, we adopt the transformations dU = dt + ω(r)dr and dφ = dφ + χ(r)dr, and we choose [44] ω The remaining steps involve setting B(r, a, θ) = H(r, a, θ) −2 [r 2 B(r) + a 2 cos 2 θ] in order to eliminate the cross term dtdr in the line element.Additionally, we select H(r, a, θ) = r 2 + a 2 cos 2 θ to satisfy G rθ = 0. Finally, this algorithm generates the stationary spacetime of the rotating dRGT (termed as RdRGT) black hole in the form in which the spin parameter is given by a = J/M , with J being the black hole's angular momentum.We have also defined 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 ∆(r) = 0, which results in the quartic equation admitting the discriminant Since Λ 1, the above discriminant can be positive within a certain range for the values of γ, and hence, the above quartic can possess four real roots, one negative and three positives.In Fig. 1, the region inside which, δ ∆ > 0, has been presented when the black hole's spin parameter varies from zero to one.Accordingly, for γ min < γ < γ max , the positivity of the discriminant is guaranteed for each value of a, and the cases of γ = γ min and γ = γ max correspond to extremal black holes.The figure shows that the width of this domain decreases with the increase of the spin parameter.Therefore, the condition δ ∆ = 0 is satisfied for both fast and slow-rotating black holes if we choose |γ| 1.Note that, the extremal limits γ min max correspond to a vanishing discriminant, for which the equation ∆(r) = 0 has one negative and three positive roots, which two of them are degenerate (coinciding horizons).On the other hand, the equation δ ∆ = 0 has one positive, one negative, and two complex conjugate roots, and hence, one can recast the discriminant in Eq. ( 16) as This quartic equation can be solved explicitly for γ min max , but we find it unnecessary to present their analytical expressions.On the 0.0 0.2 0.4 0.6 0.8 1.0 The region plot of δ∆ > 0 for Λ = 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 γmin < γ < γmax.In this case, the shaded region represents the domain corresponding to a = 0.9M .For all values of a, it holds that γmin < 0, and the width of the region δ∆ > 0 significantly increases as the spin parameter decreases.It is worth noting that for the exact values of γmin and γmax, δ∆ = 0, and possessing these values results in extremal black holes.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 γ, the solutions (A1)-(A4) 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 γ = γ max , the extremal black hole is obtained as the Cauchy and event horizons coincide, whereas for the case of γ = γ min , extremality corresponds to the unification of the event and cosmological horizons.In Fig. 2, the behavior of ∆(r) has been plotted for different values of the γ-parameter, where the solutions to ∆ = 0 and the extremal limit have been indicated.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 st++ = r 1 , r st+ = r 2 and r st− = r 3 , as those given in Eqs.(A1)-(A4) with the same form for the included parameters, considering the only replacement a 2 → a 2 cos 2 θ in Eq. (A6c).Hence, it is straightforward to see that for θ = 0, the static limits and the black hole horizons coincide.Also as expected, these radii satisfy the conditions 0 < r st− < r − and r + < r st+ < r 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. 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 = const., when the θ-coordinate oscillates between two turning points (for the case of Kerr black holes, see for example Refs.[45][46][47]).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 [10], we can determine the time-like geodesic equations in the exterior geometry of the RdRGT black hole, which are given as where we have used the dimensionless Mino time, λ, defined as Σdλ = M dτ [11], in which τ 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, defining the Carter's constant Q, as the third constant of motion.Since we are interested in the motion of massive particles, we set = 1.Note that, from Eqs. ( 17) and (18), one can infer that the particles' motion is performed under the mutual conditions R(r) ≥ 0 and Θ(θ) ≥ 0. Also for the sake of convenience, the positive segments of these equations are adopted.It is important to accentuate the role of Q in the classification of the orbits; orbits on the equatorial plane (planar orbits that correspond to θ = π/2) are obtained for Q = 0.These orbits can be regarded as the boundaries of the spherical orbits in their general case, for which, Q ≥ 0. In fact, circular orbits satisfy the conditions R(r) = 0 = R (r), which provides the critical quantities Furthermore, the effective inclination angle can be defined by means of the relation [25] Using this relation, and by eliminating L between the equations R(r) = 0 and R (r) = 0, we get to the single equation as the general characteristic equation for the radii of spherical particle orbits in the RdRGT spacetime, where in which we have used the dimensionless quantities and the definition ν = sin 2 i, for which, the Carter's constant takes the form [25] In general, the polynomial equation ( 25) has fourteen solutions of the form x(u, ν, 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. [25].One can also make a primary classification of orbits, based on the value of the k-parameter.In this regard, for 0 ≤ 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 [24].
To proceed with our discussion, we restrict ourselves to the case of planar orbits at the vicinity of the black hole's event horizon.

A. Radii of planar orbits
Planar orbits correspond to trajectories that occur in the equatorial plane, for which Q = 0 (i.e., i = 0, π, or ν = 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.
in which 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 = 1, 9), for which, we have presented some examples in the following.In Fig. 4, the k-profiles of the solutions within the domains where x i ∈ 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→j we notate the values of k, at which, the solutions x i and x j (i ≤ j) coincide.
As it is seen in the left panel of the figure, x 3 , x 4 < x + for 0 < k < k 5→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→3 = 1.838 and k 5→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→5 = 0.894 < k < k 5→3 and k > k 5→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→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→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 ∈ R have been shown for this particular case.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→3 = 0.0003 < k < k 5→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→3 and k 5→4 = 2.546.The solution x 5 has three branches which are ramified as follows.For k 4→3 < k < k 6→5 = 0.250, it is x 5 < x + , whereas for k 6→5 = 0.731 < k < k 5→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 ∈ R occurs at x = 6.840.Beyond this point, x 5 > x + is real-valued only in the region k > k 5→4 .The solution x 6 is divided into two branches outside the event horizon, which correspond to the domains k 4→3 < k < k 6→5 = 0.250 and k 6→5 = 0.731 < k < k 6 = 1.137, within x u = 0.3, h = 0.08 x 5 x 6 x 7 x 8 (a)  x u = 0.9, h = 0.08 x 5 = x 6 x 6 ∈  limit x 4 = x 5 x 8 ∈  limit FIG. 5.The x − k diagram for the real parts of the solutions xi, plotted for h = 0.08, l = 10 −6 , and u = 0.9.
which, the upper limit of x 6 ∈ 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) ≡ x(u, k 0 , h 0 , l 0 ), so that the uprofiles 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 ext = 0.923, for which, x + = x − = x ext = 0.902, giving the extremal black hole horizon where the event and Cauchy horizons coincide.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→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→4 = 1.06 > u 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→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→6 = 0.013 < u < u 8→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 ∈ 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 ext , up to u 3→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→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 (x) ≥ 0. From Eq. (21a) and the changes of variables in Eqs.(27), one can infer that 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 (x) for these two solutions are restricted to the range u 7→6 < u < u 8→7 , and in both cases R (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→5 = 0.894 ≤ 1.  I. 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 .

Non-monotonic behavior of the solutions
The above examples provide numerical solutions for determining the radii of prograde and retrograde circular orbits.Table I 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.
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 up to first order of h, that provides the critical value for which, δ ∆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 Furthermore, for k = 1 and k = 2 we obtain Now to find the desired radii, we first consider examples of h 0 and equal them to h i crit (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, {k 0 = 0, 1, 2} , 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 , {k 0 = 0, 1, 2} , 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 prograde = 1.542 and x retrograde = 2.367, corresponding to u 0 = 0.319, which are in conformity with the data given in Table I.
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 ν = 0 in the main equation (25).For the investigation of polar orbits, characterized by ν = 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 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 nonplanar 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 ν = 0 and ν = 0 respectively.

A. The latitudinal motion
Considering the definitions in Eqs. ( 27) and (28), one can re-write the evolution equation (18) as in which Z = cos θ, and for bound orbits of massive particles (i.e.k ≤ 1), we specify Naturally, for planar orbits with ν = 0, the above expression reduces to Θ Z = (1 − k)uZ 4 − (L 2 + u − uk)Z 2 , whereas for polar orbits with L = 0, it becomes Θ 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 Θ Z = 0.This equation provides the values with Z 0 = (1−k)(1−ν)u, that confine the Z-parameter.In this way, the polar angle θ oscillates in the domain θ ∈ [θ min , θ max ], in which θ min = arccos(Z max ) and θ max = arccos(Z min ).Note that, there is a sign change for the azimuth angle φ, where dφ/dλ = 0, and can be determined by means of Eq. (19).After doing the substitutions, this angular turning point is obtained as 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 ), is the Weierstrassian ℘ elliptic function with the invariants where The above solution can also be used to determine the Mino time at which the θ-profile periodically traverses the equatorial plane and generates nodes.Thus, for each full oscillation of the θ-coordinate, the profile encounters two nodes.According to Eq. ( 42) we get for θ = π/2, where 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) → (k − 1) in the quantities in Eqs.(44).

B. The azimuth motion
The solution to the evolution equation ( 19) can be written as for spherical particle orbits on the radii x i , where  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 ν = 1.Additionally, we observe that for the same inclination, an increase in the spin parameter reduces the difference between consecutive values of λ 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 ∆(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. [48].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* [49].The methods used in this paper can also be applied to study the evolution of the ISCOs in dynamical spacetimes [50].Additionally, the reflection symmetry of the RdRGT black hole can be explored, which may provide insights into the existence of curved accretion disks [51].Moreover, numerical analyses of the radii of planar orbits can be conducted for various axially symmetric spacetimes, serving as tests for the analytical solutions [52].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 [53,54].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 [55].All of these investigations can be further extended to explore regular rotating black holes and solitons [56].Recently, a new concept of particle surfaces has been defined in Ref. [57].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.

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

FIG. 3 .
FIG.3.The radial profile of gtt plotted for θ = π/4 for the two cases of a = 0.3M and a = 0.9M , and for the same values of the γ-parameter as in Fig.2.The first two roots of gtt = 0 are shown for all of the cases, and for the extremal black holes, only two static limits are available.

FIG. 4 .
FIG.4.The x − k diagrams for the real parts of the solutions xi 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 x7 and x8 are real, has been magnified in panel (b).The color-coding of the solutions xi that is used here will be also applied in all of the forthcoming diagrams within the paper.

FIG. 8 .
FIG.8.Some examples of spherical particle orbits in accordance with the data presented in TableII.The sphere indicates the closure of points swept by the radii xi, which is cut to halves by a circle on the θ = π/2 surface.
k u x prograde x retrograde k u x prograde x retrograde k u x prograde x retrograde