On Innermost Stable Spherical Orbits near a Rotating Black Hole: A Numerical Study of the Particle Motion near the Plunging Region

According to general relativity, astrophysical black holes are described by a small number of parameters. Apart from the mass of the black hole (M), among the most interesting characteristics is the spin (a), which determines the degree of rotation, i.e., the angular momentum of the black hole. The latter is observationally constrained by the spectral and timing properties of the radiation signal emerging from an accretion disk of matter orbiting near the event horizon. In the case of the planar (standard, equatorial) accretion disk, this is the location of the innermost stable circular orbit that determines the observable radiation characteristics and allows us to measure the spin. In this paper, we discuss a more general case of the innermost stable spherical orbits (ISSOs) extending above and below the equatorial plane. To this end, we study the nonequatorial geodesic motion of particles following inclined, spherical, relativistically precessing trajectories with the aim of exploring the boundary between the regions of stable (energetically bound) and escaping (energetically unbound) motion. The concept of the radius of the ISSO should play a role in determining the inner rim of a tilted or geometrically thick accretion flow. We demonstrate that the region of inclined bound orbits has a complicated structure due to enhanced precession near the inner rim. We also explore the fate of particles launched below the radius of the marginally bound spherical orbit: these may either plunge into the event horizon or escape to radial infinity.


INTRODUCTION
A growing number of studies have been published during the recent decade indicating that black holes are ubiquitous in different types of cosmic objects, ranging from stellar-mass black holes in compact binary systems to supermassive black holes that reside in cores of galaxies (Meier 2012).The processes of accretion and ejection of matter onto black holes turn out to be essential for our ability to identify cosmic black holes and to explore their interaction with the surrounding gaseous environment and the radiation field (Frank et al. 2002;Shakura 2018, e.g.).Strong-gravity effects operate in the regions close to the black hole event horizon, and so the General Relativity framework has to be employed in order to describe the motion of particles and fluids (Chandrasekhar 1998;Kato et al. 2008;Misner et al. 2017).
The gravitational field near black holes can be described in terms of Kerr metric (Kerr 1963) with mass and spin being the only two astrophysically relevant parameters.The motion of test particles around classical black holes is regular, i.e., without signatures of chaos (Carter 1968).The special importance has been attributed to the innermost stable circular orbit of the equatorial motion (ISCO Bardeen et al. 1972;Mummery & Balbus 2022), when test particles transit from the phase-space region of energetically stable circulation into plunging trajectories.Even in the planar case, complete taxonomy is very rich, as can be seen in a systematic classification by Levin & Perez-Giz (2008).
Pressure and turbulence must affect the corresponding structure of gaseous flows, nevertheless, ISCO plays an important role in the standard accretion disc scenario, and it may also influence gaseous flows to a certain extent.This raises a long-standing question about a boundary between stable (energetically bound), plunging, and escaping (unbound) trajectories for non-equatorial motion; despite its apparent relevance for geometrically thick and/or inclined accretion flows, the topology of that boundary has been so far explored only indirectly, and the results about a potential analogy of ISCO outside the equatorial plane have been presented only in an implicit form (Will 2012;Rana & Mangalam 2019;Teo 2021;Tavlayan & Tekin 2021).
In this context, we revisit the issue and give a detailed account of the problem of stability and energetic binding of spherical orbits including a numerical analysis of an unstable case occurring in an inner region of an accreting black hole system.The paper is organized as follows.In Section 2 we review the equations of motion and we introduce a two-dimensional effective potential for the geodesic motion of free test particles in Kerr spacetime.In Section 3 we review the properties of spherical orbits using two alternative parameterizations (in terms of Carter constant or inclination angle) and we study the locations of the innermost stable spherical orbits and marginally bound spherical orbits with respect to parameters of the system.The results are numerically verified in Section 4 with the technique of escape-boundary plots.Furthermore, we focus here on the unstable orbits below ISSO and employ this method to demonstrate, how the parameters (namely, spin and Carter constant) affect the instability and what effect they may have on the probability of the plunge or escape.Results are summarized and briefly discussed in Section 5. Details on numerical methods employed for the analysis are provided in Appendix.

EQUATIONS OF MOTION AND EFFECTIVE POTENTIAL
Kerr metric describes the geometry of the spacetime around the rotating black hole (Kerr 1963).It can be expressed in Boyer-Lindquist coordinates x µ = (t, r, θ, φ) as follows (Misner et al. 2017): where ∆ ≡ r 2 − 2M r + a 2 , Σ ≡ r 2 + a 2 cos 2 θ. (2) A coordinate singularity at ∆(r) = 0 corresponds to the position of the outer and the inner event horizons of the black hole, respectively: r ± = M ± √ M 2 − a 2 .The rotation of the black hole is measured by the spin parameter a ∈ ⟨−M, M ⟩; we can assume a ≥ 0 without any loss of generality.
We note that geometrized units will be used throughout the paper, making the gravitational constant and the speed of light equal to one; G = c = 1.For the rest of the paper, we switch to dimensionless units scaled by the mass of the black hole, i.e., M = 1 is set in all formulas.
The Hamiltonian H of a particle with rest mass m in the spacetime with contravariant metric components g µν may be defined as (Misner et al. 2017): where p µ is the generalized (canonical) momentum which in this case corresponds to the kinematical momentum.The equations of motion are expressed as: where λ ≡ τ /m is a dimensionless affine parameter (τ denotes the proper time).The conserved value of the Hamiltonian is given as H = −m 2 /2.The system is stationary and the time component of canonical momentum π t is an integral of motion, which equals the (negatively taken) energy of the test particle, p t ≡ −E.Moreover, the system is axially symmetric which assures the conservation of the axial component of the angular momentum L ≡ p φ .Analysis of the geodesics in Kerr spacetime using the Hamilton-Jacobi formalism reveals (Carter 1968) that besides the above obvious integrals E, L corresponding to relevant Killing vectors, there exists another independent integral of motion, namely, the Carter constant Q, related to the hidden symmetry of the system corresponding to the existence of the Killing tensor.Carter constant may be expressed as: (5) In some cases, the constant is being expressed in an alternative way as K ≡ Q + (L − aE) 2 , however, for our analysis, the form given by Eq. ( 5) remains more appropriate.Endowed with the four independent integrals of motion which are in involution (i.e., they commute in terms of Poisson brackets) the system becomes fully integrable in the Liouville sense and its Hamilton-Jacobi equation is completely separable (e.g., Goldstein et al. 2002).Besides other properties that follow from integrability, we stress that all trajectories in such a system are regular, i.e., no chaotic orbits are present in the whole phase space.
In the rest of the paper, we switch to specific quantities E/m → E, L/m → L, and Q/m 2 → Q which corresponds to setting the rest mass of the particle m = 1 in the formulas.
Effective potential expressing the minimal allowed energy of a particle might be expressed as a function of coordinates r, θ and parameters a, L in the following form (Kopáček & Karas 2018a): where Equipotential curves V eff = E drawn in the (r, θ) plane connect simultaneous turning points of the particle with energy E in coordinates r and θ.Analysis of the two-dimensional effective potential may thus be used to localize circular orbits not only in the equatorial plane but also off this plane.Although in the pure Kerr spacetime, the circular geodesics are allowed only in the equatorial plane, if we, for instance, consider electrically charged particles near weakly magnetized Kerr black hole, we also find circular orbits above and below equatorial plane (Kovář et al. 2008(Kovář et al. , 2010)), and the method of two-dimensional effective potential also allows to discuss their stability (Tursunov et al. 2016).

SPHERICAL ORBITS IN KERR SPACETIME: THEIR STABILITY AND BINDING
Spherical orbits are a special class of bound geodesics in Kerr spacetime with a constant radial coordinate r = const and latitudinal angle θ varying around the equatorial plane between the turning points that are positioned symmetrically at θ ⋆ and π − θ ⋆ .The value of θ ⋆ is determined by the value of Carter constant Q and, in particular, for Q = 0 the spherical orbit reduces to circular (Keplerian) trajectory in the equatorial plane with θ ⋆ = π/2.An illustration of spherical orbits with different values of Carter constant is provided in Fig. 1.They can represent trajectories of test particles near the event horizon.
Spherical orbits of massive particles in the Kerr metric were analyzed by Wilkins (1972) for the special case of extreme Kerr black hole with a = 1, while Goldstein (1974) performed the early numerical calculations of such geodesics (see also Dymnikova 1986).Energy and angular momentum of the spherical orbits may be expressed explicitly as a function of radius, spin, and latitudinal turning point θ ⋆ (Shakura 1987): where Σ ⋆ = r 2 + a 2 cos 2 θ ⋆ and q ⋆ = (r 2 − a 2 cos 2 θ ⋆ ) 1/2 .Upper signs correspond to co-rotating (direct) orbits while the lower signs are valid for counter-rotating orbits.With the inclination θ ⋆ = π/2 the above expressions reduce to familiar formulas for circular (Keplerian) orbits in the equatorial plane (Bardeen et al. 1972): Example of a set of spherical orbits launched at the radius r = 5 from the equatorial plane of a Kerr black hole (a = 0.8) differing solely by the value of Carter constant Q which specifies the initial value of p θ = Q 1/2 .In particular, we compare the trajectories with Setting the other limiting value, θ ⋆ = 0, in Eqs. ( 8) and ( 9) we obtain polar orbits with L = 0 and energy given as: At a given radius, polar orbits periodically cross both poles while being continuously dragged due to the rotation of the black hole as illustrated in Fig. 2.
a=0.1 a=0.5 a=0.8 a=1 Alternatively to the parameterization by the value of the turning point θ ⋆ , it is also possible to parameterize the integrals of motion E and L by the Carter constant Q (Teo 2021): where Υ = r 5 − Q(r − 3)r 3 + a 2 Q 2 .Submitting Q = 0, we may again verify that the above formulas reduce to expressions for circular orbits given by Eqs. ( 10) and ( 11).
In Fig. 3 we plot the square root of the Carter constant Q 1/2 (which is equal to the equatorial value of p θ ) as a function of the angular momentum given by Eq. ( 14) for spherical orbits at r = 10 with several values of spin.In particular, we observe that for a given radius and spin, the maximal value of the Carter constant is attained by the counter-rotating spherical orbit (L < 0) while the value corresponding to polar orbit with L = 0 is slightly lower.Each value of Q 1/2 below its maximum corresponds to the pair of orbits with different L, which might be both negative (two counter-rotating orbits) for the Carter constant close to the maximal value, or with mixed signs of L (one co-rotating and one counter-rotating) for lower values.The described behavior of the function Q 1/2 (L) is not universal and it changes closer to the horizon.As we shall see, neither polar nor spherical orbits may generally exist all the way down to the horizon.
We note that alternative expressions for the energy and angular momentum of the spherical orbit have also been provided by Rana & Mangalam (2019) in Eqs.(18a), (18b) and (18c) of the referred paper.These equations, however, are not entirely correct.Although they work well for low values of Q for which they reproduce proper parameters of spherical orbits, with higher values of Q, they fail to do so and provide erroneous values of the angular momentum (see the discussion and Fig. 22 in Appendix).
Figure 4. Locations of the innermost stable spherical orbits (left panel) and marginally bound spherical orbits (right panel) with different values of the Carter constant and spin.On the vertical axis, parameter a[M] ≡ σa is shown, where a is the Kerr spin parameter, for which a ≥ 0 is assumed throughout the paper, and σ = ±1 is a switch that distinguishes prograde vs. retrograde motion.Positive values (σ = 1) correspond to orbits co-rotating with the black hole, while negative values (σ = −1) are for the counter-rotating orbits.These graphs generalize the well-known dependence of the equatorial ISCO radius to the case of orbits inclined with respect to the equatorial plane.
In the equatorial plane, stable circular geodesics are not possible below the radius of Innermost Stable Circular Orbit (ISCO), the position of which is given as (Bardeen et al. 1972): where 1 and the upper sign corresponds to co-rotating orbit and the lower sign to counter-rotating one.In fact, bound circular orbits exist only above the radius of marginally bound circular orbit (MBCO) given as (Bardeen et al. 1972): The concept of ISCO has profound importance in our understanding of how standard accretion disks operate in the black hole vicinity.Indeed, spectral characteristics of accreting black holes, both the supermassive ones in cores of Active Galactic Nuclei (AGN) as well as stellar-mass black holes in close binary systems, indicate that the accretion disk often proceeds down to ISCO, although a growing number of counter-examples have been indicated in various sources and related to numerical simulations of truncated accretion disks with strong magnetic fields (Narayan et al. 2003;Liska et al. 2019).Moreover, recently introduced "puffy" disks seem to proceed well below ISCO (Lančová et al. 2019;Wielgus et al. 2022).Nonetheless, much of the emerging flux is typically released in X-rays that are produced just above the critical ISCO radius.As the black hole spin grows, the ISCO location approaches the event horizon, and the radiation output from the accretion mechanism is also thought to grow in parallel.On the other hand, accretion does not have to be bound to the equatorial plane and, indeed, numerous examples have been reported of objects where the accretion flow extends above the equatorial plane as it proceeds via a non-planar (twisted) accretion disk (e.g.Kumar & Pringle 1985;Martin et al. 2007).This brings us to an important question about the analogy of the equatorial ISCO radius outside the equatorial plane.We can introduce the notion of Innermost Stable Spherical Orbit (ISSO) as a natural extension of ISCO.
In a close analogy with the case of circular orbits residing in the equatorial plane, for spherical (generally inclined) orbits the radius r = r s of ISSO, and the related radius r = r b of the Marginally Bound Spherical Orbit (MBSO) are given implicitly by the following algebraic relations (Rana & Mangalam 2019): The meaning of θ⋆ is the angular distance of the latitudinal turning point from the zenith or nadir, respectively (see the main text for further details).Dashed lines correspond to the counter-rotating spherical orbits.We notice that the vertical axes of the plots are inverted, i.e., their lower edges (θ⋆ = π/2) correspond to circular orbits in the equatorial plane, while polar orbits with θ⋆ = 0 are found on the upper edges.
and These expressions for ISSO and MBSO radii were obtained using the equation of the separatrix curve in the (e, ρ) parameter plane of eccentricity e and inverse-latus rectum ρ of the trajectory.In particular, the ISSO radius is derived by inserting e = 0 and ρ = 1/r s , while for the MBSO, the values e = 1 and ρ = 1/2r b need to be used.The polynomial that defines the separatrix is thus reduced to the degree of nine and eight, respectively.See Rana & Mangalam (2019) and especially Appendix D therein for further details on the derivation of these relations.
In the relevant range of spin a ∈ ⟨0, 1⟩, both equations ( 17) and ( 18) have two real roots such that r > r + .The higher value of each pair corresponds to the counter-rotating spherical orbit, while the smaller value locates the corotating orbit.In the Schwarzschild limit (a = 0), the values coincide at r s = 6 and r b = 4 corresponding to circular orbit in the static Schwarzschild spacetime.In Fig. 4, we plot the radii of ISSO and MBSO as functions of spin for different values of the Carter constant for both directions of orbits.In particular, we note that ISSO and MBSO radii of counter-rotating spherical orbits are located below those of circular orbits while the co-rotating spherical orbits have these radii always higher than circular orbits.
a=0.1 a=0.3 a=0.5 a=0.7 a=0.9 a=1 In Figs. 5 and 6, we provide an alternative discussion of ISSO and MBSO radii and corresponding ratio L/E with respect to the inclination angle θ ⋆ of spherical orbits with different values of spin.In particular, we notice that stable nor bound counter-rotating spherical orbits (dashed curves in Figs. 5 and 6) can never have arbitrarily low inclination θ ⋆ .Indeed, even if we start with the stable counter-rotating orbit (i.e. with a radius above the corresponding ISSO), we inevitably lose stability if we gradually decrease the inclination (i.e., increase the Carter constant of the particle) before the orbit becomes polar.On the other hand, in the case of co-rotating orbits, we can, in principle, have a stable orbit of an arbitrary inclination, provided that the radius is sufficiently large to remain above the corresponding ISSO curve.In particular, with the appropriate choice of the Carter constant, Q = Q p , we get the polar orbit with θ ⋆ = 0.
The value of Q p corresponding to a polar orbit may be expressed as: where E p is given by Eq. ( 12), while g tt and g θθ are contravariant components of the Kerr metric (1).Values of Q p for various spins are plotted in Fig. 7. Marginal stability and marginal binding of polar orbits may also be investigated using Eqs.( 17) and ( 18), if the value Q = Q p is inserted.The resulting radii of the marginally bound polar orbit (MBPO) and innermost stable polar orbit (ISPO) are plotted as a function of spin in Fig. 8.In particular, we observe that MBPO changes between r ≈ 3.38 (with a = 1) and r = 4 for a = 0 (circular orbit around Schwarzschild black hole).Boundary values of ISPO radii are r ≈ 5.28 (corresponding to a = 1) and r = 6 in the Schwarzschild limit.
In the above discussion, we have studied spherical orbits of massive particles using two alternative parameters (with one-to-one correspondence), namely, the Carter constant Q and the inclination angle θ ⋆ .In particular, radii of the innermost stable spherical orbit (ISSO) and marginally bound spherical orbit (MBSO) were discussed for different spin values using both complementary parameterizations.Special attention was paid to the limiting case of polar orbits θ ⋆ = 0 for which the corresponding radii of the innermost stable polar orbit (ISPO) and marginally bound polar orbit (MBPO) were also determined.In the following Section 4, we shall discuss the behavior of particle orbits below the ISSO (or ISPO), and we will investigate the nature of their instability numerically.
We note that spherical orbits in Kerr spacetime may also be followed by photons (Teo 2003).It appears, however, that in the region of our interest (i.e., above the outer horizon at r = r + ) these are always unstable against perturbations in the radial direction.The discussion of spherical photon orbits is considerably simplified compared to the case of massive particles, nevertheless, both cases still share some interesting common features.Moreover, photon orbits near rotating black holes are also astrophysically relevant as they define the apparent shapes (shadows) of the black holes for the distant observers (Bardeen 1973;Young 1976) Unlike massive particles, spherical photon orbits are obtained as a one-parameter family of orbits, i.e., for a given black hole spin a the orbit is fully defined by radius r, while in the case of timelike geodesics, we always have one more parameter to select (see Eqs. ( 8)-( 9) and Eqs. ( 13)-( 14)).Moreover, it appears that for null geodesics described by integrals E, L and Q, only the two ratios ϕ ≡ L/E and η ≡ Q/E 2 are really independent (Bardeen 1973).These ratios are directly related to the impact parameters of photons received by a distant observer and are thus relevant for the analysis of black hole shadows and gravitational lensing.In the case of spherical photon orbits their values are given as (Teo 2003): Spherical orbits are only allowed in the range of radii r 1 ≤ r ≤ r 2 , where r 1,2 = 2 1 + cos 2 3 arccos (∓a) are the radii of unstable circular photon orbits in the equatorial plane with η = 0. Radius r 1 corresponds to a co-rotating (prograde) orbit while at r 2 we get a counter-rotating (retrograde) circular orbit.The radial range between r 1 and r 2 is divided by an intermediate radius corresponding to the polar orbit with ϕ = 0.For r 1 ≤ r < r 3 the orbits are co-rotating (ϕ is positive) while for r 3 < r ≤ r 2 we get counter-rotating orbits with negative ϕ.The value of η is maximized within the range of radii of counter-rotating orbits which is analogical to the case of massive particles where also the maximal value of the Carter constant corresponds to the counter-rotating orbit (see Fig. 3).In the limit of maximally rotating black hole, a = 1, the boundary radii of spherical photon orbits become r 1 = 1, r 2 = 4 and r 3 = 1 + √ 2. For a = 0, they reduce to r 1 = r 2 = r 3 = 3, which defines the photon sphere of the Schwarzschild black hole.
For more details on spherical photon orbits outside the Kerr black holes, we refer to Teo (2003).General analysis of photon orbits in Kerr and Kerr-Newman spacetimes is given by Gal 'tsov & Kobialko (2019).A detailed review on black hole shadows is provided by Perlick & Tsupko (2022).In the rest of the article, we further analyze the spherical orbits of massive particles only.

NUMERICAL ANALYSIS OF STABLE AND UNSTABLE SPHERICAL ORBITS
From the analysis performed in Section 3 we learn two basic facts regarding the stability of spherical orbits with given values of Q and a: i) the radii r s of innermost stable spherical orbits (ISSO) below which the spherical orbits become unstable, and, ii) the radii r b of marginally bound spherical orbits (MBSO) below which they become unbound (their energy exceeds the rest energy, i.e., E > 1 in dimensionless units) and may escape to infinity.Generally, r b ≤ r s , while the both radii only coincide in the case of co-rotating circular orbits (Q = 0) around maximally spinning black hole for which r b = r s = r ± = a = 1.6) of the particle with angular momentum given by Eq. ( 14) while the equipotential curve of corresponding energy given by Eq. ( 13) is marked by the dashed black line.The solid black line marks the horizon of the black hole.The blue curve shows the projected trajectory.In the right panel, the same trajectory is shown in three dimensions with the horizon marked by the grey sphere.Following parameters were employed: r0 = 6.5, Q 1/2 = 0.75 and a = 0.5.
An unstable spherical orbit maintains its initial radius and keeps evolving as a spherical orbit only if it remains completely unperturbed.However, in a real system (physical or numerically simulated), the fluctuations that perturb the dynamics are always present, and unstable orbits are therefore eradicated.As a result, a perturbed spherical orbit with initial radius r b < r < r s may, in principle, plunge into the horizon or quasiperiodically oscillate in both radial and latitudinal directions.On the other hand, particles launched with r + < r ≤ r b may also escape to infinity as they have sufficient energy to surpass the attraction of the center.
In Figs.9-12 we present typical examples of the above-mentioned trajectories.Each trajectory is visualized in two ways.In the left panels, we show the projections of the trajectories to the poloidal plane (x, z) being bounded by the relevant zero velocity curves (given by the equipotential curves of the effective potential Eq. ( 6)) as well as the shape of the potential surface in its neighborhood, while the right panels of the respective figures show the full 3D trajectory above the black hole horizon.Examples in Figs.9-12 share the same value of the Carter constant and spin parameter (Q 1/2 = 0.75 and a = 0.5), while they differ in the initial radii (and thus the values of E and L given by Eq. ( 13) and Eq. ( 14), respectively).
A stable spherical orbit above ISSO is shown in Fig. 9.For an unstable orbit launched between radii of ISSO and MBSO, there are two options, i.e., plunge into the horizon (Fig. 10) or quasiperiodic oscillations in radial and latitudinal directions (Fig. 11) while it largely depends on the initial perturbation of an unstable equilibrium which one will be realized.The particle launched below MBSO also has two options depending on the perturbation of the unstable orbit; it may plunge into the horizon or escape to infinity (both cases are shown in Fig. 12).The behavior of particles set on unstable geodesics below ISSO and their possible stabilization on quasiperiodic orbit due to small perturbations represent astrophysically interesting problems that may be tackled by various approaches.In general, the relevant phase space consists of regions corresponding to different modes of motion which are separated by boundary denoted as separatrix.In particular, here we identify three different modes of behavior, i.e., radially bounded motion, plunge onto the horizon, and escape to infinity.Especially the separatrix between plunging and stable regions of Kerr geodesics is particularly relevant for the discussion of Extreme Mass Ratio Inspirals (EMRIs) which represent one of the key observational targets for the future LISA mission.In this context, it has been recently systematically studied using the parameterization of orbits by eccentricity and semilatus rectum (Stein & Warburton 2020).
The class of escaping orbits of particles forming an outflow of matter from the vicinity of the spinning black hole is also of particular astrophysical interest, especially in the context of observed high-energy cosmic rays.While the acceleration mechanism beyond cosmic rays is still a matter of debate and particular mechanisms are being discussed, active galactic nuclei powered by spinning supermassive black holes represent the main suspects for accelerating observed extra-galactic cosmic rays (e.g., Rodrigues et al. 2021;Tursunov et al. 2020;Kološ et al. 2021).While the cosmic rays are composed dominantly of charged particles (protons and electrons, in particular), the electrically neutral component is also present, although probably not reaching ultra-high energies (Navarro & Pierre Auger Collaboration 2013).Particles escaping from an unstable equilibrium of spherical orbit below MBSO (illustrated in Fig. 12) may produce an outflow of electrically neutral particles if they are perturbed outward in the radial direction (δ(p r ) > 0).
In our previous studies (Kopáček & Karas 2020, 2018b) we investigated a different dynamical system with escapes (Karas & Kopáček 2021;Stuchlík & Kološ 2016) by a straightforward, yet effective, approach based on escape-boundary plots.Such graphic representation of the dynamics is inspired by standard basin-boundary plots which visualize basins of attraction of attractors in systems with dissipation.Although conservative systems do not have real attractors we may consider escaping orbits (if present) as being attracted by an attractor at infinity (Contopoulos & Harsoula 2010;Contopoulos 2006).If we pursue this analogy, we may also consider the horizon of the black hole as an attractor for plunging orbits and explore which parts of the phase space belong to which attractor and determine the initial conditions that belong to their basins of attraction.If we visualize a particular set of trajectories (differing in initial conditions or other parameters) in a way that distinguishes them (e.g., by assigning different colors) by their final states (i.e., plunge, escape, or bounded motion, in our case) we obtain an escape-boundary plot as an analogy of basin-boundary plot for the conservative system with escapes.
The method of escape-boundary plots has been successfully applied in a previously studied non-integrable system of charged particles in magnetized Kerr background (Kopáček & Karas 2018b;Karas et al. 2013).Due to the nonintegrability of this system (e.g., Mukherjee et al. 2023), we have encountered fractal structures of the boundaries with self-similar patterns which were associated with chaotic dynamics (Kopáček & Karas 2020) in which case specific methods, e.g., recurrence analysis (Marwan et al. 2007), become especially useful (Kopáček & Karas 2020; Lukes-Gerakopoulos & Kopáček 2018).Nevertheless, in a current study we examine a completely integrable system of Kerr geodesics where the deterministic chaos cannot develop, and the boundaries between domains of stable motion, plunge, and escape should therefore be regular (i.e., with non-fractal dimension) and, in particular, they should respect ISSO and MBSO radii determined in Section 3.
In order to numerically verify the above-given assumptions, we construct the escape-boundary plots of two types: i) with coordinates r 0 and Q 1/2 and fixed values of spin shown in Fig. 13, and ii) with coordinates r 0 and a and fixed Q presented in Fig. 14.The special case of polar orbits is treated separately in Fig. 15.To obtain these plots we use the grid with resolution 400 × 400, where each point represents a spherical trajectory launched at r 0 in the equatorial plane with constants E and L given by Eqs. ( 13) and ( 14).We evolve each trajectory until λ fin (see Appendix 5 for details) and assign the color to each point of the grid according to the following convention: yellow for escape (r fin ≥ r e ), red for radially bounded trajectory with r + < r fin < r e (i.e., not necessarily with constant r), and blue for the plunge.The choice of the escape threshold radius r e is discussed in Appendix 5.In addition, we use black color for the black hole interior (r ≤ r + ) and grey for the parts of the parameter space where the spherical orbits are not defined (values given by Eqs. ( 13) and ( 14) become complex).We also present MBSO (cyan line) and ISSO (green line) radii in these plots.Moreover, for the case of escape boundary plots in (r 0 , Q 1/2 ) plane shown in Fig. 14, we highlight the parameters' values corresponding to polar orbits by a white line.The case of polar orbits (L = 0; θ ⋆ = 0) is presented in Fig. 15, where the energy of the test particles is set according to Eq. ( 12) and the value of the Carter constant is given by Eq. ( 19).
Inspecting the escape-boundary plots presented in Figs.13-15 we find that they confirm the previous analysis of ISSO and MBSO locations presented in Section 3. Indeed, we observe that: i) the regions above ISSO lines are entirely red since spherical orbits launched here are stable against small numerical perturbations caused by integration errors, and ii) yellow zones of escaping particles only appear below MBSO lines where the energy of particles becomes sufficient for the escape.However, in regions between MBSO and ISSO lines and below the MBSO line, we observe an irregular mixture of different colors; blue and red dots in the former case, and blue and yellow dots in the latter case.
For the case of circular orbits (upper left panel of Fig. 13) the distribution of different orbits in these regions appears random while it becomes more organized for spherical orbits (Q > 0) and large escape zones also develop here.Especially in Figs. 14 and  .Left panel: a sketch of the dynamics in the neighborhood of a saddle-type unstable hyperbolic fixed point in a two-dimensional dynamical system with phase-space coordinates (q, p) = (0, 0).Separatrix (red) and flow lines (blue) in different phase-space regions are shown.Right panel: example of actual dynamics near an unstable spherical orbit between MBSO and ISSO radii with parameters r0 = 4.1, Q 1/2 = 3 and a = 0.9 (light-blue asterisk).The blue region corresponds to plunging orbits, red region consists of radially bounded quasiperiodic orbits.Poincaré surface of section (taken at θ = π/2 for θ > 0) with canonical coordinates (r, pr) of several quasiperiodic trajectories within the bound region differing in initial radii is shown with black color.
Karas 2020).Nevertheless, since here we investigate a fully integrable system, these unexpected features associated with deterministic chaos are necessarily of a numerical origin, i.e., caused by numerical errors of the integration.Although the integration errors may be largely decreased by an appropriate choice of the integrator (discussed in Appendix 5), they can never be completely avoided.The reason why these small numerical perturbations are so crucial in the region below ISSO is obviously the fact that in this zone we are dealing with the dynamics of an unstable equilibrium.Dynamical systems close to such an instability must be studied with special caution and appropriate tools must be used.Artificial (numerically induced) fractal structures observed in Figs.[13][14][15] show that for the regions below ISSO the escape-boundary plots of this type are not an optimal method and may even become misleading.
Based on the analysis of escape-boundary plots in Figs.13-15 we may only conclude, that parameters a and Q strongly affect the behavior of numerically perturbed unstable spherical orbits.In particular, regarding the red-blue zone between MBSO and ISSO, these plots suggest that for co-rotating orbits the Carter constant increases the probability that a perturbed orbit becomes radially bounded, while for counter-rotating orbits it slightly increases the probability of the plunge.The spin parameter appears to have the opposite effect; slightly increasing the probability of plunge for co-rotating orbits while supporting the stabilization of perturbed counter-rotating orbits into quasiperiodic radially bounded orbits.Nevertheless, since the numerical perturbation is not random and non-trivially depends on the integration scheme, such probabilistic interpretation is necessarily vague and needs to be confirmed by other means.Moreover, regarding the yellow-blue zone below MBSO we may draw hardly any conclusion here and irregular fractal-like structures present in the plots make the applicability of this approach in this zone even more questionable.Q 1/2 =2 6 6.05 6.1 6.15 Figure 17.Effect of the Carter constant on the dynamics near unstable spherical orbits (marked by light-blue asterisks in the plots) located between MBSO and ISSO radii.The upper row shows the co-rotating orbits near the unstable orbit at r0 = 3.24 with spin a = 0.7, while the bottom row reveals the counter-rotating orbits around unstable orbit r0 = 6 with the spin a = 0.8.The range of pr values (vertical axis) is different in each row in order to make the effect of Q 1/2 more apparent.
In order to assess the stability of a periodic orbit (or a corresponding fixed point of the Poincaré map) and the dynamics in its neighborhood, the usual procedure is to linearize the system (considering only the first-order terms of the corresponding Taylor expansion) and solve the eigenvalue problem.(for details see, e.g., Strogatz 2019;Tabor 1989).Based on the eigenvalues of the relevant Jacobian matrix (also denoted as stability matrix in this context), which in our case consists of the second derivatives of the Hamiltonian (3) with respect to the phase-space variables x µ and p µ , we may introduce a following classification and terminology of the fixed points (Wiggins 2006).Namely, if all the eigenvalues have nonzero real parts, the fixed point is denoted as hyperbolic.If some, but not all, of the eigenvalues have real parts greater than zero (resp., moduli greater than one) and the rest of the eigenvalues have real parts less than zero (resp., moduli less than one), then the hyperbolic point is called a saddle.If all the eigenvalues have negative real parts (resp., moduli less than one), then the hyperbolic fixed point is called a stable node or sink.If all of the eigenvalues have positive real parts (resp., moduli greater than one), then the hyperbolic fixed point is called an unstable node or source.For each type of hyperbolic fixed point, the structure of the (linearized) flow in its vicinity significantly differs.Moreover, for the case of a purely imaginary set of eigenvalues of the stability matrix, we obtain a fixed point of different class, namely, an elliptic fixed point.For further details and precise mathematical treatment of the problem of linearized stability, we refer to Wiggins (2006).For our purpose, the basic classification of the fixed points associated with the periodic spherical orbits remains sufficient.
In particular, regarding the zone between MBSO and ISSO the above analysis of linearized dynamics indicates that initial conditions of a spherical orbit given by Eqs. ( 13) and ( 14) correspond to an unstable saddle-type hyperbolic fixed point of the relevant Poincaré map.A sketch of the dynamics near such a point in a phase space of a two- a=-0.9 5.9 5.95 6 6.05 a=-1 Figure 18.Effect of the spin parameter on the dynamics near unstable spherical orbits (marked by blue asterisks in the plots) located between MBSO and ISSO radii.The upper row shows the co-rotating orbits near the unstable orbit at r0 = 4.1, with Q 1/2 = 3, while the bottom row reveals the counter-rotating orbits around the unstable orbit at r0 = 5.9 with Q 1/2 = 1.8.
dimensional dynamical system is presented in the left panel of Fig. 16.In that case, the trajectories near the fixed point are exactly hyperbolic (blue curves in the sketch) while in the fixed point itself (located at the origin of the plot), the stable (attracting) and unstable (repelling) manifolds (indicated by red lines) intersect.Stable manifolds of the fixed point consist of orbits that reach it asymptotically (as time t → ∞) while those of unstable manifolds do so as t → −∞.Together they form a boundary that separates different types of trajectories (different modes of behavior of a dynamical system) and is thus usually denoted as separatrix.For precise definitions of the above terms and more detailed treatment of the topic, we refer to standard textbooks on nonlinear dynamical systems (e.g., Ott 1993;Lichtenberg & Lieberman 1992).
In the case of saddle-type unstable fixed points, we may observe a special class of orbits, namely homoclinic orbits, which are found at the intersections of the stable and unstable manifolds and connect the fixed point to itself (while the heteroclinic orbit connects different fixed points).For the investigated system, however, only the former is relevant.In particular, for the case of circular geodesics in Kerr spacetime, it has been previously shown, that there is a oneto-one correspondence between energetically-bound unstable circular orbits (i.e., circular orbits with radii between MBCO and ISCO) and homoclinic orbits (Levin & Perez-Giz 2009;Perez-Giz & Levin 2009).Moreover, it has been demonstrated that homoclinic orbits represent a limiting case of zoom-whirl behavior of orbits with extreme perihelion precession occurring in a strong-field regime.
The zoom-whirl behavior of geodesics in Kerr background is obviously not restricted to unstable circular orbits and may also be observed in connection with unstable spherical orbits (see, e.g., Fig. 11).To our knowledge, the properties of homoclinic orbits associated with energetically bound unstable spherical orbits (i.e., with radii between MBSO and ISSO) have not been studied in detail.Nevertheless, there are further qualitative indications based on numerical analysis which allow us to establish such correspondence.In particular, in the right panel of Fig. 16, we combine the escape-boundary technique to distinguish plunging initial conditions (blue) from those leading to radially bounded motion (red), and thus locate the separatrix of the fixed point (light-blue asterisk), with the method of Poincaré surface of section recording the one-way intersections of trajectories with the equatorial plane θ = π/2 and θ > 0 (black).These trajectories, differing in r 0 , are launched from the equatorial plane with u r (0) = p r (0) = 0 and they share the same values of energy and angular momentum given by Eqs. ( 13) and ( 14) at the unstable point at r 0 = 4.1.For sufficiently long integration, the intersection points draw closed curves typical for quasiperiodic orbits which are characterized by incommensurable frequencies.Only at r 0 ≈ 5.1 we get the stable periodic orbit of constant r, i.e., stable spherical orbit above ISSO which corresponds to a fixed point of a Poincaré map.
In order to determine the role of the Carter constant and spin in the evolution of unstable spherical orbit under perturbation, we employ the technique of escape-boundary plots once again.However, unlike the application in Q 1/2 =3.5 Figure 19.The effect of the Carter constant on the dynamics near unstable spherical orbits (marked by blue asterisks in the plots) that are located below MBSO radius.The upper row shows co-rotating orbits near the unstable orbit at r0 = 1.7 for the spin a = 0.9.The bottom row reveals the counter-rotating orbits around the unstable orbit r0 = 4.5 with the spin a = 0.6.Note that the range of pr values (vertical axis) is different in each row.
Figs. 13-15, now we shall study the vicinity of the particular unstable orbit in the plane (r 0 , p r ).In Fig. 17 we observe how the slope of separatrix, which divides the blue plunging zone from the red wedge of bounded motion, changes with different values of the Carter constant (while all other parameters of the orbits are kept constant).In the first row of Fig. 17 we show an unstable co-rotating orbit (r 0 = 3.24, a = 0.7) with different values of the Carter constant.The slope of the separatrix rises and the opening angle of the red stable zone gradually and significantly increases as the Carter constant increases.If we consider the direction of the perturbation of the unstable equilibrium to be random in the (r 0 , p r ) plane, then the probability of the stabilization of the perturbed co-rotating orbit clearly increases with the Carter constant.On the other hand, in the case of counter-rotating orbits (with parameters r 0 = 6 and a = 0.8), the probability of stabilization is decreased by the Carter constant, although the effect is not as prominent as in the co-rotating case (see the second row of Fig. 17).
In Fig. 18 we perform an analogous discussion of the role of the spin parameter.Therefore for a particular unstable orbit, we fix the Carter constant (here we choose a co-rotating orbit with r 0 = 4.1, Q 1/2 = 3, and counter-rotating orbit with r 0 = 5.9, Q 1/2 = 1.8) and vary the spin.As previously, for each case, we evolve a grid of 400 × 400 initial conditions in (r 0 , p r ) plane around the unstable point.Unlike the effect of the Carter constant, the spin is observed to decrease the probability of the stabilization of co-rotating orbits (first row of Fig. 18) and increase it for counterrotating orbits (second row).As for the Carter constant, the effect is stronger for the co-rotating orbits, although the difference is not that significant.
We may conclude that both parameters (the Carter constant and the black hole spin) play an important role in the evolution of unstable spherical orbits with radii between MBSO and ISSO.In the case of co-rotating orbits, the Carter constant increases the probability of stabilization, while the spin suppresses it.For the counter-rotating orbits, the role of the parameters is reversed -the Carter constant supports the tendency to plunge while the spin tends to stabilize the perturbed orbits.Of course, the probabilistic interpretation of Figs. 17 and Fig. 18 is based on the assumption that the direction of the perturbation in the (r 0 , p r ) plane is random, and thus the probability of stabilization in the form of a quasiperiodic eccentric orbit is proportional to the opening angle of the red wedge of stable motion.These observations directly and more clearly confirm the assumptions previously inferred from Figs. 13-15.
The above results apply to the region between MBSO and ISSO.Below MBSO the orbits become energetically unbound allowing the particles to escape to infinity and the nature of instability changes from the saddle type to an unstable node.Nevertheless, we may use the analogous approach as previously to study this part of the phase space.In Fig. 19 we study the effect of the Carter constant on a particular spherical orbit below MBSO.Namely, in the upper row of Fig. 19 we pick co-rotating orbits with r 0 = 1.7 and a = 0.9, while the bottom row reveals the counter-rotating orbits with r 0 = 4.5 and a = 0.6.It appears that in this case, the Carter constant has a rather negligible effect on the destiny of particles.Indeed, we observe that it only slightly rotates the separatrix between the yellow escape zone and the blue plunging region.With increasing Carter constant the rotation of the separatrix is clockwise for corotating orbits and counterclockwise for counter-rotating orbits.Nevertheless, the separatrix remains straight and the probability of the escape or plunge for the random perturbation therefore does not change.The fractal-like structures observed below MBSO in Figs.13-15 were therefore a product of numerical noise lacking physical significance.

CONCLUSIONS
Spherical orbits in Kerr spacetime are a special class of geodesics with constant radii.They are defined with respect to a suitably selected radial coordinate and generalize the equatorial circular orbits by allowing non-zero values of the Carter constant.These orbits can be interpreted as inclined (by the angle θ ⋆ with respect to the rotation axis) circular geodesics; the orbital plane is continuously rotated by the frame-dragging effect.The latter effect due to the spinning black hole has no classical analog and it represents a purely relativistic phenomenon.The spherical orbits were first studied in the early 1970s (Wilkins 1972) and have recently attracted a renewed interest (Teo 2021;Stein & Warburton 2020;Rana & Mangalam 2019) due to their relevance for EMRIs (e.g., Amaro-Seoane 2018), which represent a key target for the future gravitational observatory LISA (Amaro-Seoane et al. 2023).In this paper, we have revisited the topic, focusing on the dynamics of both stable and unstable spherical orbits.
In Section 3 we first discussed the radii of the Innermost Stable Spherical Orbit (ISSO) and the Marginally Bound Spherical Orbit (MBSO) with respect to the values of the spin parameter and the Carter constant.We introduced the parameterization of orbits by the (or inclination) angle θ ⋆ as an alternative to the Carter constant.The parameterization by θ ⋆ provides an intuitive picture of the geometric shape of the orbit, and the two limiting special cases, namely, polar orbits (θ ⋆ = 0) and circular equatorial orbits (θ ⋆ = π/2) are easily identified.For a given orientation of the orbit (co-rotating or counter-rotating), there is a one-to-one correspondence between θ ⋆ and the Carter constant.The analysis shows that for co-rotating orbits, both ISSO and MBSO radii are gradually shifted to higher radii compared to ISCO and MBCO of the corresponding circular orbit as the Carter constant increases.On the other hand, for counter-rotating spherical orbits, ISSO and MBSO are always smaller than those of circular orbits, and they gradually shrink as the Carter constant is increased.
Spherical orbits below ISSO become unstable, and the fate of a particle, i.e., whether it plunges to the black hole, stabilizes on a quasiperiodic orbit, or escapes (from below MBSO), depends primarily on the phase-space direction of the perturbation of the unstable equilibrium.In Section 4 we numerically study the nature of the instability using the technique of escape-boundary plots (Kopáček & Karas 2023, 2020) to reveal how the parameters, namely the spin and the Carter constant, affect the probability of each possible outcome.It appears that a straightforward application of this method may become questionable since the system near instability is prone to numerical noise.Therefore, we adapt the method and plot the neighborhood of the unstable spherical orbit in the plane of the canonical variables (r, p r ), which allows us to combine the plot with the Poincaré sections of particular trajectories and, more importantly, we can thus directly trace the separatrix between the regions of the plunge and escape (occurring below MBSO) or between the plunge and radially bounded motion (above MBSO).Comparing the results for different values of the parameters (and assuming that the perturbation has a random direction in the (r, p r ) plane), allows us to conclude that the spin parameter decreases the probability of the stabilization for co-rotating unstable spherical orbits while it increases this probability for counter-rotating orbits.The Carter constant has an opposite role: it contributes to the stabilization of co-rotating orbits while it makes counter-rotating orbits more prone to plunge.Nevertheless, below the MBSO radii, the nature of the instability changes.Namely, it becomes an unstable node for which the probability of plunge or escape does not change with these parameters since the separatrix remains a straight line that only slightly rotates as the spin or Carter constant increases.
Our results have been obtained with a general numerical approach of escape-boundary plots, whose applicability is not limited to a fully integrable problem of Kerr geodesics and can also be used to study a chaotic system under non-integrable perturbations (e.g.Kopáček & Karas 2020, 2018b).On the other hand, the integrable case allows to proceed analytically.In particular, with an appropriate parameterization, it is possible to find an implicit formula for the separatrix between the plunging and quasiperiodic regions of the phase space (Stein & Warburton 2020).More recently, an analytical solution for plunging geodesics in Kerr has also been discussed (Dyson & van de Meent 2023).In this context, the numerical method used in our analysis provides a complementary and intuitive approach that allows be especially sensitive to errors.In Fig. 20, we observe that secular numerical error may introduce here an artificial plunging zone above ISSO as an obvious numerical artifact which grows with the integration time λ fin .Nevertheless, the integration time must be sufficient for the escaping orbits to reach the escape threshold r e which is another crucial parameter of escape-boundary plots which needs to be set with caution.In Fig. 21 we observe, how the improper choice of r e may lead to incorrect identification of orbits.In the first two panels the underestimated value of r e leads to the false identification of highly eccentric quasiperiodic orbits near MBPO boundary as escaping.On the other hand, an overestimated escape threshold may cause misinterpretation of escaping orbits below MBPO as quasiperiodic (shown in the third panel of Fig. 21).
After testing different values of the parameters λ fin and r e , we found the optimal setting for the escape-boundary plots in Figs.13-15 to be λ fin = 2 × 10 4 with r e = 1000.For the survey of the vicinity of the unstable fixed point in Figs.17-19, a shorter integration remains sufficient, for which also the threshold needs to be adjusted.Namely, we set λ fin = 2000 and r e = 100 for these calculations.
In addition to the numerical problems mentioned above, we have also come across an erroneous result in the literature.Namely, the expression of the angular momentum of a spherical orbit as a function of the Carter constant published by Rana & Mangalam (2019) in formulas (18a)-(18c) of the cited paper gives incorrect results for the co-rotating nearpolar spherical orbits, although it works well for counter-rotating orbits and for co-rotating orbits with lower values of the Carter constant (see Fig. 22).Therefore, for our calculations, we employ an alternative expression Eq. ( 14) by Teo (2021), which works correctly regardless of the value of the Carter constant (see Fig. 3).

Figure 3 .
Figure 3. Square root of the Carter constant of spherical orbits at r = 10 as a function of angular momentum L for several values of spin.In the right panel, we zoom the section of the left plot around maximal values.

Figure 5 .
Figure5.Locations of the innermost stable spherical orbits (left panel) and marginally bound spherical orbits (right panel) with respect to the inclination θ⋆ for different values of the spin parameter.The meaning of θ⋆ is the angular distance of the latitudinal turning point from the zenith or nadir, respectively (see the main text for further details).Dashed lines correspond to the counter-rotating spherical orbits.We notice that the vertical axes of the plots are inverted, i.e., their lower edges (θ⋆ = π/2) correspond to circular orbits in the equatorial plane, while polar orbits with θ⋆ = 0 are found on the upper edges.

Figure 6 .
Figure 6.Ratio L/E (a.k.a.impact parameter) as a function of the latitudinal turning point θ⋆ for spherical orbits at ISSO radius (left panel) and at MBSO radius (right panel) for different values of spin.Solid lines correspond to co-rotating orbits, whereas the dashed lines are for counter-rotating orbits.

Figure 7 .
Figure7.Square root of Carter constant Qp of polar orbits is plotted with respect to their radii for various values of the spin parameter.Solid lines are used for stable orbits above ISPO, while dashed lines correspond to unstable polar orbits below ISPO.
. The concept of black hole shadow has recently been observationally confirmed with spectacular images from Event Horizon Telescope showing the supermassive black hole in M87 galaxy (Event Horizon Telescope Collaboration et al. 2024) and Sagittarius A* in the center of the Milky Way (Event Horizon Telescope Collaboration et al. 2022).

Figure 9 .
Figure 9. Example of a stable spherical orbit above ISSO.The left panel shows the effective potential Eq. (6) of the particle with angular momentum given by Eq. (14) while the equipotential curve of corresponding energy given by Eq. (13) is marked by the dashed black line.The solid black line marks the horizon of the black hole.The blue curve shows the projected trajectory.In the right panel, the same trajectory is shown in three dimensions with the horizon marked by the grey sphere.Following parameters were employed: r0 = 6.5, Q 1/2 = 0.75 and a = 0.5.

Figure 13 .
Figure 13.Escape-boundary plots with coordinates (r0, a) and fixed values of the Carter constant.Spherical trajectories presented in these plots are color-coded as follows: yellow for escaping, red for bound, and blue for plunging orbits.Black color indicates the black hole interior and grey shows the regions above the horizon where spherical orbits cannot exist.Locations of MBSO are shown by a cyan line and ISSO by a green line.Negative spin values correspond to counter-rotating orbits.

Figure 14 .
Figure14.Escape-boundary plots with coordinates (r0, Q 1/2 ) and fixed values of spin parameter.Spherical trajectories presented in these plots are color-coded as in Fig.13.Besides the MBSO radii (cyan line) and ISSO radii (green line), we also show the locations of polar orbits by a white line.

Figure 15 .
Figure 15.Escape-boundary plot of polar orbits is shown in (r0, a) plane.Color coding as in Figs. 13 and 14.MBPO radii are indicated by a cyan line and ISPO radii by a green line.

Figure 21 .Figure 22 .
Figure21.Seemingly escaping polar orbits may appear above MBPO if the escape threshold radius re is set too low (first two panels).Overestimating its value may, on the other hand, lead to a false interpretation of escaping orbits as radially bounded (third panel).The same magnified section of Fig.15is shown in all panels which only differ in re.