Revisiting critical orbits of test particles traveling in a black hole background

Critical orbits are defined as those trajectories for which the radial kinetic energy function admits either a double or triple real root, corresponding to an unstable circular orbit or motion that neither falls into the black hole nor escapes to infinity. Critical orbits are powerful tools for analytically studying black hole shadows and accretion processes. However, as a specialized subclass of geodesics, they are rarely examined as stand-alone entities in the literature. In some literature, critical orbits may be overlooked in the broader discussion. The primary objective of this paper is to conduct a comprehensive review of the key properties of critical orbits, aiming to renew researchers' attention toward this significant topic in the field.

The study of geodesic motion in black hole spacetimes has a long and rich history, closely intertwined with the development of the solutions themselves. In spherically symmetric spacetimes, the earliest analytical solution for geodesic motion in Schwarzschild spacetime can be traced back to Droste [1] in 1917, who expressed his solution in terms of the Weierstrass elliptic function. Thirteen years later, Hagihara conducted a comprehensive classification of all possible test particle motions in Schwarzschild spacetime, and his work [2] has remained a classic reference. Simultaneously with the development of solutions based on Weierstrass functions, a significant body of research [3−8] by Forsyth, Greenhill, Darwin, and Scharf et. al. succeeded in expressing Schwarzschild geodesics using Jacobi elliptic functions and Legendre integrals, with publication years spanning several decades. In Reissner–Nordström (RN) spacetime, the motion of test particles exhibits richer behavior due to the presence of charge; for instance, a particle crossing the event horizon may not necessarily end at the singularity but can traverse the Cauchy horizon and emerge into another universe [9]. Grunau and Kagramanova [9] provided an analytical solution for the motion of electrically and magnetically charged test particles in this background.

In rotating, axisymmetric spacetimes, the theory of Kerr geodesics developed rapidly after Carter's fundamental discovery in 1968 that the Hamilton–Jacobi equation is completely separable [10, 11]. In this work, he identified a conserved quantity now universally known as the Carter constant, the geometrical nature of which was later understood by Walker and Penrose in terms of Killing tensors [12]. Early analyses focused on equatorial orbits [13−15], while non-equatorial geodesics proved to be much more involved. Wilkins [16] provided an early account of bound orbits, and vortical orbits were first discussed by De Felice and Calvani [17, 18]. A major technical breakthrough came in 2003, when Mino [19] introduced the so-called Mino time, a parameter that completely decouples the radial and latitudinal equations of motion, enabling a straightforward application of elliptic functions. This facilitated the derivation of analytic solutions for bound timelike orbits by Fujita and Hikida [20]. The phase-space structure of Kerr geodesics, including homoclinic orbits and the separatrix, has been investigated by Levin and collaborators [21, 22]. More recently, null geodesics in the Kerr exterior have been revisited by Gralla and Lupsasca [23, 24], providing convenient classifications and analytic solutions relevant to observational problems like black hole lensing and photon rings. And Cieślik et. al. [57] derive analytical solutions describing timelike and null geodesics by using Weierstrass elliptic functions in the Kerr spacetime. For the Kerr–Newman spacetime, which incorporates charge, the geodesic motion of both neutral and charged particles has been extensively studied. Early work by Johnston and Ruffini [26] examined timelike equatorial and spherical orbits of uncharged particles, while Young [27] analyzed the last stable orbit for charged particles. Bičák et al. [28, 29] provided a systematic study of charged particle motion, including radial motion and motion along the symmetry axis. Later, Kovář et al. [30] discovered unstable off-equatorial circular orbits for charged particles, and Pugliese et al. [31] used equatorial circular orbits to distinguish between black holes and naked singularities. Hackmann and Xu [32] completely classifies the colatitudinal and radial motion of charged test particles in the Kerr-Newman spacetime and presents analytical solutions in terms of elliptic functions that are valid for all types of orbits. A comprehensive analysis of photon orbits in Kerr–Newman spacetime was presented by Calvani and Turolla [33]. For more recent literature, see [34−38].

In the case of spacetimes with a cosmological constant, Hackmann and Lämmerzahl achieved a significant breakthrough by performing analytical integrations for Schwarzschild–de Sitter [39, 40] and Reissner–Nordström–de Sitter geometries [41] using hyperelliptic θ and σ functions, based on the Jacobi inversion problem restricted to the θ-divisor. This method was subsequently extended to higher-dimensional Schwarzschild and Reissner–Nordström spacetimes with a cosmological constant [41]. Similar hyperelliptic function methods have also been successfully applied to the axially symmetric Taub-NUT [44] and Kerr–de Sitter [43] spacetimes, where the types of orbits are classified and extensively studied. In addition to the literature mentioned above, a number of books also discuss geodesic motion in black hole spacetimes in considerable detail, e.g. [45−48].

Another important goal of this paper is to lay the groundwork for a 3+1 formalism describing the accretion of a Vlasov gas onto a Kerr–Newman black hole. In 2017, Rioseco and Sarbach [49] developed a modern analytical theory for relativistic collisionless Vlasov gas accretion onto a Schwarzschild black hole. In their model, the tangential pressure exceeds the radial pressure at the horizon [50], offering a partial explanation for the issue of low accretion rates, which has attracted considerable attention. Subsequently, this model was quickly extended to describe accretion onto a moving Schwarzschild [51, 52], Reissner–Nordström [53] and Kerr [54−56] black holes. In a related development, Mach et al. constructed a numerical framework using Monte Carlo methods to simulate the motion of collisionless gas in a Schwarzschild background [57, 58]. Since our focus is on accretion models, we are particularly interested in the scenario where test particles fall toward the black hole from infinity. Similar studies of unbound orbits in Schwarzschild spacetime have already been discussed, see [59, 60]. It should be noted that the accretion of neutral particles onto a black hole does not induce charge evolution. To investigate the effect of accretion on black hole charge, we have recently extended this model to describe the accretion of a collisionless Fermi gas onto a Reissner–Nordström black hole [61].

In the accretion model proposed by Rioseco and Sarbach [49], test particles are assumed to be in thermal equilibrium at infinity. Under gravitational influence, these particles follow geodesics into the finite region. Consequently, they naturally separate into three categories: those that fall into the black hole and are absorbed, those that are scattered by the black hole and return to infinity, and—lying between these two scenarios—a critical case in which particles neither enter the black hole nor escape to infinity, but instead asymptotically approach an unstable bound orbit. Therefore, examining this critical behavior is essential in accretion theory, as it defines the boundary in parameter space between absorbed and scattered particles. A similar line of reasoning applies to the analysis of the shadow boundary of a black hole. Therefore, this paper provides a systematic review of methods for handling critical geodesics in both spherically symmetric and axisymmetric rotating spacetimes. It aims to enable efficient identification of the relevant parameters associated with critical geodesics and to facilitate their analytical or numerical analysis.

In this paper, we systematically review the classical results related to critical geodesics for massless and neutral massive particles in the double root case, and present the relevant results for the triple root case and for charged particles in the Reissner–Nordström and Kerr–Newman backgrounds. And the paper is organized as follows. In section II, we review the derivation of the equations of motion for charged test particles in Kerr–Newman spacetime using the separability of the Hamilton–Jacobi equation. Sections III through VI present a detailed analytical treatment of critical orbits — both null and timelike — in Schwarzschild, Reissner–Nordström, Kerr, and Kerr–Newman black hole spacetimes, respectively. The final section provides conclusions and discussion.

In this section, we review the equations of motion for charged particles moving in the Kerr-Newman spacetime. We work in the Boyer-Lindquist coordinates $ (t,r,\theta,\varphi) $. The Kerr-Newman metric is well known and can be written as

where

and M is the mass, Q is the charge, and $ a = J/M $ is the angular momentum per unit mass of the black hole. The outer and inner horizons are the roots of $ \Delta(r) = 0 $ and are given by

In general, the condition for a black hole to possess an event horizon is given by $ 0\leq a^2+Q^2\leq M^2 $. The electromagnetic four-pontential in this background is $ A_{\mu} = (-\frac{Q r}{\rho^2},\frac{Q r}{\rho^2}a\sin^2\theta,0,0) $. The inverse metric can be expressed as

which satisfies $ g^{\mu \alpha}g_{\alpha \nu} = \delta^{\mu}_{\nu} $.

Consider a test particle with mass m and charge q moving in the background of a rotating, charged black hole, the Lagrangian for its motion can be written as

where $ e = \frac{q}{m} $ is the charge-to-mass ratio, λ is the affine parameter, and the dot denotes differentiation with respect to λ. The canonical momenta $ p_{\mu} $ are defined by

The Hamiltonian is obtained by performing a Legendre transformation

There are four constants of motion: the rest mass m, the energy E, the Carter constant D, and the angular momentum $ L_z $ in the z-direction, which are defined by

With these constants, the equations of motion $ \dot{x}^{\mu} = \frac{\partial {\cal{H}}}{\partial p_{\mu}},\dot{p}_{\mu} = -\frac{\partial {\cal{H}}}{\partial p_{\mu}} $ can be rewritten as a set of decoupled, first-order equations

where we have introduced the variables $ R = \rho^4\dot{r}^2 $ and $ \Theta = \rho^4\dot{\theta}^2 $ for brevity.

We now introduce Hamilton's principal function S via

and the Hamilton–Jacobi equation for S then takes the form

By demonstrating that the Hamilton–Jacobi equation (21) for test particles is separable, Carter [10] first established the existence of a fourth constant of geodesic motion in Kerr spacetime – now known as the Carter constant D. As a consequence of this separability, the abbreviated action S is constrained to take a particular form

Beacause D is constant, generic geodesic motion is not planar; thus, the orbits lie in the three-dimensional coordinate space $ (r,\theta,\varphi) $. The equations of motion are obtained by requiring that the partial derivatives of the abbreviated action S with respect to the four constants $ D,E,m,L_z $ vanish. Setting $ \frac{\partial S}{\partial D} = 0 $, one yields the constraint satisfied by the $ (r,\theta) $ coordinates

Similarly, by setting$ \frac{\partial S}{\partial m} = 0 $, $ \frac{\partial S}{\partial E} = 0 $, and $ \frac{\partial S}{\partial L_z} = 0 $, we find

The above expressions are obtained by simplification using Eq.(23).

Consider test particles incident from infinity. Some of these particles will fall into the black hole, whereas others will be scattered back to infinity. Between these two cases lies a critical regime in which particles are neither captured by the black hole nor do they escape to infinity; instead, they asymptotically approach a critical radius at a finite distance from the black hole. We next provide a detailed analysis of the critical orbital solutions for test-particle motion, including both null and timelike orbits, in Schwarzschild, Reissner-Nordström, Kerr, and Kerr-Newman spacetimes.

In this section, we consider the motion of neutral particles around a Schwarzschild black hole. In this case, the parameters satisfy $ a = Q = e = 0 $. The equations of motion (16) - (19) reduce to

where $ L^2 = D+L_z^2 $ denotes the squared angular momentum. In spherically symmetric spacetimes, the Carter constant D is related to the conserved angular momentum $ L^2 $, which can be further expressed as

where $ d\Omega^2 = d\theta^2+\sin^2\theta d\varphi^2 $. When the square of the angular momentum $ L^2 $ is constant, it follows that the geodesics in this spacetime are confined to a single plane. This conclusion also holds for any spherically symmetric spacetime. To simplify the calculations, we restrict our discussion to orbital motion in the equatorial plane, $ \theta = \frac{\pi}{2} $. Orbits in other planes can be obtained by appropriately rotating the equatorial plane. In the equatorial case, there are $ \dot{\theta} = 0 $ and $ L = L_z $. We are not concerned with the solution for the time coordinate t, and the orbital equations (28) - (30) further reduce to

Since the condition $ \left( \frac{dr}{d\varphi} \right)^2 \geq 0 $ must hold, particles arriving from infinity must satisfy $ E \geq m $. The case $ E = m $ corresonds to particles falling freely from rest and is not a critical orbit. Therefore, we restrict our analysis to $ E \gt m $. Defining $ u = 1/r $, Eq. (32) can be rewritten as

In the case of a null geodesic $ (m = 0) $, the function $ f(u) $ simplifies to:

The orbital types are classified based on the roots of $ f(u) = 0 $. Since $ f(0) = \frac{E^2}{L^2}>0 $, there is always a negative root. The three possible orbital types — absorption, scattering, and critical — correspond to the nature of the remaining roots; see Fig. 1. For the critical orbit, the remaining two positive roots coincide at $ u_c $, that is:

Figure 1.  We consider all possible root configurations of $ f(u) = 0 $. In Case 1, where there are no real roots for $ u>0 $, a photon incident from infinity is absorbed by the black hole. In Case 2, characterized by two distinct real roots for $ u>0 $, the photon is scattered. Case 3 corresponds to the critical orbit, where a double root occurs at $ u_c>0 $. In this critical case, as we will discuss later, the photon is neither absorbed nor scattered but instead approaches a circular orbit at a finite radius $ r_c = \frac{1}{u_c} $.

where prime denotes differentiation with respect to the variable. Solving these equations yields

The critical radius $ r_c = \frac{1}{u_c} = 3M $ defines the photon sphere, whereas the critical impact parameter $ R_c = \frac{L_c}{E_c} = 3\sqrt{3}M $ corresponds to the radius of the black hole's shadow as seen by an observer at infinity.

Substituting the critical impact parameter into (34), the orbital equation reduces to

where the negative root is $ u_{*} = -\frac{u_c}{2} $. The solution is given by:

where $ \varphi_0 $ is an integration constant. The critical null geodesic is shown in Fig. 2, depicting the trajectory of a massless particle originating at infinity and spiraling inward, asymptotically approuching the circular at $ r_c = 3M $.

Figure 2.  Critical orbits of null geodesics in the Schwarzschild metric are shown. The mass is set to $ M = 1 $, with the constant $ \varphi_0 = 0 $. The dashed line indicates the location of the photon sphere.

Unlike the case of null geodesics, the critical radius for timelike geodesics depends on the parameters E and L. ¹ When $ E^2>m^2 $, we have $ f(0) = \frac{E^2-m^2}{L^2}>0 $, which imples the existence of a single negative root. If the other two real roots coincide, then the conditions (35) are also satisfied, namely

The solutions are given by

The conditions $ L_c^2>0 $ and $ E_c^2-m^2>0 $ imply $ u_c\in(\frac{1}{4M},\frac{1}{3M}) $. Substituting the critical impact parameter into equation (33) reduces orbital equation to

where the negative root is $ u_{*} = \frac{1}{2M}-2u_c $ and $ u_c $ lies in the interval $ (\frac{1}{4M},\frac{1}{3M}) $. Mathematical analysis shows that

The critical timelike geodesic in the Schwarzschild spacetime is shown in Fig. 3.

Figure 3.  Critical orbits for timelike geodesics in the Schwarzschild metric. The parameters are chosen as $ M = 1,\varphi_0 = 0 $ and $ u_c = \frac{7}{24} $. The dashed line marks the critical radius $ r_c = \frac{1}{u_c} $.

In the accretion model of Rioseco and Sarbach [49], the key observable quantities are the particle current density $ J_{\mu} $ and the energy-momentum tensor $ T_{\mu\nu} $, which are defined by

where $ f(x,p) $ is the distribution function and $ \mathrm{dvol}_x(p) = \sqrt{-\det[g^{\mu\nu}]}dp_tdp_rdp_{\theta}dp_{\varphi} $ is the volume element. Introducing the action-angle momenta $ (m,E,l_z,l) $, the volume element can be re-expressed as $ \mathrm{dvol}_x(p)\propto dmdEdLdL_z $. For the integral over $ dm $, it is generally assumed that the Vlasov gas consists of a large number of identical point particles with mass $ m_0 $, and therefore the distribution function includes a factor $ \delta(m-m_0) $. For the integral over $ dE $, unbound orbits statisfy $ E>m_0 $. For the integral over $ dL_z $, the integration limits are $ L_z\in[-L\sin\theta, L\sin\theta] $. The remaining integrals over $ dL $ separate the observable quantities $ J_{\mu} $ and $ T_{\mu\nu} $ into three parts: the absorbed, the scattered, and the critical parts, which are:

Since the critical part occurs at $ L = L_c $, we have $ \int_{L_c}^{L_c}f(l)dl = 0 $, so the critical part does not contribute substantially. We can re-express (41) and (42) as

For the absorption component, the interval is $ L\in [0,L_c] $. For the scattering component, the interval is $ L\in [L_c, L_{m}] $, where $ L_{m} $ is chosen so that $ f>0 $ throughout. Solving $ f = 0 $ yields

Thus, we have obtained all integration intervals. In fact, we do not need to express them explicitly in the form $ L_c(E) $, as the relation is already fully specified by equations (41) and (42). Therefore, we will not discuss the integration interval of the accretion model separately in the subsequent sections.

In this section, we study the motion of charged particles in the Reissner-Nordström black hole spacetime with $ a = 0 $. The equations (16) - (19) governing this motion reduce to

As in Schwarzschild spacetime, the square of the angular momentum $ L^2 $ is conserved. Therefore, the worldlines are confined to a single plane. Note that a charged test particle in Reissner-Nordström spacetime does not, in general, follow geodesics. The orbital equation in the plane is given by $ \left(\frac{du}{d\varphi}\right)^2 = f(u) $, where

Now, $ f(u) $ is a quartic function in u. Given that $ f(0) \gt 0 $, the equation $ f(u) = 0 $ must have one negative real root $ u_{*}^{-} $ and one positive real root $ u_{*}^{+} $. Therefore, the zeros of the function $ f(u) $ may fall into one of the following cases: (i) four real roots, three of which coincide; (ii) four real roots, two of which coincide; (iii) two real roots and two complex roots. This paper primarily concerns with cases (i) and (ii). Both cases satisfy $ f(u_c) = 0,f'(u_c) = 0 $.

The null geodesic $ (m = e = 0) $ in the Reissner-Nordström geometry is governed by the function

In the critical case $ f(u_c) = 0 $ and $ f'(u_c) = 0 $, it is straightforward to obtain

Photons incident from infinity will first encounter the larger critical radius $ r_c = \frac{1}{u_c} $; this is the case we focus on in the subsequent analysis.

As in the Schwarzschild case, the shadow radius of a Reissner–Nordström black hole, as seen by a distant observer, is given by $ R_c = \frac{L_c}{E_c} $. For small values of the charge Q, a Taylor series expansion yields

Thus, relative to an uncharged black hole, a charged black hole has a smaller photon sphere and therefore casts a smaller shadow. In the extremal case $ Q = M $, the photon sphere and shadow attain their minimal radii

Given the condition for the existence of a horizon $ Q \lt M $, the case of three coincident real roots cannot occur. Thus, in the critical case, the function $ f(u) $ can only degenerate into

where

and $ u_{*}^{-}<0<u_c<u_{*}^{+} $. Solving the equation $ \frac{du}{d\varphi} = \sqrt{f(u)} $, we obtain

The plot of this solution closely resembles Fig. 2 and is therefore omitted here.

The conditions for timelike critical orbits are as follow:

If the critical orbits have radius $ r_c = \frac{1}{u_c} $, then the energy E and the angular momentum are given by

where

These equations require that $ 4m^2Y_c+e^2Q^2u_c^2>0 $, which implies $ u_{c}^{-}<u_c<u_{c}^{+} $, where

A fundamentally new scenario arises—the case of triply coincident roots. For this to occur, the additional constraint $ f''(u) = 0 $ must be satisfied, namely,

Eliminating $ L^2 $ using Eq. (68), we obtain:

Note that, while the above equation may admit multiple solutions, only those with $ u_c $ satisfying $ u_{c}^{-}<u_c<u_{c}^{+} $ and corresponding to a triple root are valid. In particular, when $ e = 0 $, the above equation reduces to

This equation determines the minimum radius of a stable circular orbit in Reissner-Nordström geometry. In the critical case (i), the planar orbital equation takes the form:

where the negative root $ u_{*}^{-} = \frac{2M-3Q^2u_c}{Q^2} $ shows that the triple root satisfies $ u_c>\frac{2M}{3Q^2} $. The solution is

An example of this orbit is plotted in Fig. 4.

Figure 4.  Shown are the critical orbits for case (i) of a timelike worldline in the Reissner–Nordström metric. The parameters are chosen as $ u_c = \frac{2}{3M},e = 1.5m,Q = \frac{3 m M}{\sqrt{e^2+8 m^2}},M = 1,\varphi_0 = 0 $. The dashed line indicates the location of the critical radius $ r_c = \frac{1}{u_c} $.

In the critical case (ii), the planar orbital equation takes the following form:

where

where $ L_c^2 $ is given by Eq. (68). This equation has the same form as the null-geodesic orbital equation (62); hence, its solution can be written as Eq. (64), but with different parameters $ u_c $ and $ u_{*}^{\pm} $.

In this section, we study the critical orbits of neutral particles in the Kerr geometry. For $ Q = e = 0 $, equations (16) - (19) reduce to

Here, Eq. (19) reduces to Eq. (83).

In the following discussion, we introduce new variables:

Then, the functions R and θ can be rewritten as follow:

We begin by considering the solution in the θ-direction, before discussing the r-motion. The θ-motion is governed by the integral $ \int\frac{d\theta}{\sqrt{\Theta}} $. For this integral to be real-valued, the condition $ \eta+(a-\xi)^2 \geq 0 $ must be satisfied. Let $ \cos\theta = \mu $ be the variable of integration; then $ \int\frac{d\theta}{\sqrt{\Theta}} = -\int\frac{d\mu}{\sqrt{\tilde{\Theta}}} $ where

The behavior of the θ-motion is categorized into three cases: $ \eta>0, \eta<0 $, and $ \eta = 0 $.

1. $ \eta>0 $: We can re-express $ \tilde{\Theta} $ as follow:

where

The requirement $ \tilde{\Theta}\geq 0 $ implies that $ 0\leq \mu^2\leq \mu_2^2 $. Here, $ \mu^2 = 0 $ corresponds to the equatorial plane $ \theta = \frac{\pi}{2} $. When $ \mu^2 = \mu_2^2 $, the orbits attain the turning angles at $ \theta = \frac{\pi}{2}\pm \theta_m $, where $ \theta_m = \arccos \mu_2 $. In other words, these orbits cross the equatorial plane and oscillate symmetrically about it. The integration yields

where $ \text{EllipticF}[\phi,m] $ denotes the incomplete elliptic integral of the first kind, and C is an integration constant.

2. $ \eta = 0 $: In this special case, we have

where $ \mu_m^2 = 1-\frac{\xi^2}{a^2} $, and $ 0\leq \mu^2\leq\mu_m^2\leq 1 $. The integration yields:

3. $ \eta<0 $: In this case, the function $ \tilde{\Theta} $ can be re-expressed as

where

Thus, the parameter $ \mu^2 $ is restricted to $ 0\leq\mu_1^2\leq\mu^2\leq \mu_2^2\leq 1 $. In fact, if the $ \eta<0 $ case actually happens, it requires very stringent conditions: $ |\eta|<a^2 $ and $ 0<|\xi|\leq|a| $. In such cases, the test particle will inevitably enter the region of negative r, and its worldline will terminate this negative r region. However, to date our analysis has not identified a physically feasible set of parameters that would realize this scenario. We include it here solely as a mathematical possibility for completeness. The integration yields

We then consider the radial motion. For massless particles traveling in Kerr geometry ($ m = 0 $ and $ \chi = 1 $), R reduces to

We consider the critical orbits $ R = 0 $ and $ R' = 0 $

The equations above yield two distinct pairs of solutions for ξ and η; one is given below

and the other is

Substituting solution (100) into equation (86), we find that $ \Theta = -\frac{\rho^2 E^2}{a^2 \sin^2\theta} \leq 0 $. Since physical motion requires $ \Theta \geq 0 $, this condition can only be satisfied if $ \theta = \theta_0 $ remains constant, corresponding to shear-free null congruences. Our primary interest lies in the second set of solutions (101).

Celestial coordinates were introduced in [45]

Here, $ r_0 $ is the distance from the black hole to the observer, and $ \theta_0 $ is the observer's viewing angle. If the observer lies in the equatorial plane, $ \theta_0 = \frac{\pi}{2} $, the celestial coordinates can be written as $ \alpha = -\xi $ and $ \beta = \pm\sqrt{\eta} $. The apparent shape of the black hole is determined by the critical coordinates $ (\alpha_c,\beta_c) $. Figure 5 depicts the shadow of the Kerr black hole as viewed from the equatorial plane.

Figure 5.  The shadow of a Kerr black hole, as observed by a distant observer in the equatorial plane, is shown. The mass of the black hole is set to $ M = 1 $.

There are two types of critical orbits, defined by the root structures of $ R = 0 $: (i) triple coincident roots and (ii) double coincident roots. In addition to Eqs. (98) and (99), the case of triple coincident roots requires that the condition $ R'' = 0 $ also be satisfied, i.e.,

Define $ \mu_0^2 = \frac{M r_c}{3M^2-3M r_c+r_c^2} $; one can then obtain the solution

Rewriting relation (105) yields

This expression reveals a monotonic relationship: $ a^{*}_c $ increases with $ r_c $. Given the constraints for the existence of the horizon, $ 0\leq a\leq M $, and that the photon sphere lies outside the horizon, $ r_{+}\leq r_c $, it follows that the case of triply coincident roots occurs only when $ a^{*}_c = M $ and $ r_c = r_{+} = M $, which implies $ \mu_0^2 = 1 $, with $ \xi_c = 2M, \eta_c = 3M^2 $. Substituting these solutions into R and $ \tilde{\Theta} $ yields

Mathematically, we obtain

where $ C_1, C_2 $ are integration constants. Since $ \eta = 3M^2>0 $, the θ-motion oscillates between $ (\frac{\pi}{2}-\theta_{m},\frac{\pi}{2}+\theta_{m}) $, where $ \theta_{m} = \arccos \mu_{m} $. The constraint on $ (r,\theta) $ given in (23) can now be integrated to yield

where $ r_i, \mu_i $ are the initial conditions. Here, the limits $ \mu_m^{n} $ alternate between $ \mu_m $ and $ -\mu_m $, namely $ \mu_m^{2k+1} = \pm\mu_m $ and $ \mu_m^{2k} = \mp\mu_m $. The number of integrals n on the right-hand side equals the number of oscillations, which in turn depends on the final radial coordinate r. The solution for φ can be written as

where $ \text{ElliptiPi}[n,\varphi,m] $ denotes the incomplete elliptic integral of the third kind. The numerical solution is shown in Fig. 6.

Figure 6.  Critical null orbits with triple coincident roots in Kerr spacetime. The left panel illustrates the trajectory of a null geodesic in the $ (r,\theta) $ plane; the right panel shows the trajectory in three-dimensional space. The mass of the black hole is set to $ M = 1 $, and the initial conditions are $ r_{i} = 10,\mu_{i} = 0 $, and $ \varphi_i = 0 $. The red dashed line (ball) denotes the phonton sphere, and the black dashed lines correspond to $ \mu = \pm \mu_{m} $.

The critical null geodesics with double coincident roots can be analyzed in a similar manner. Substituting $ \xi = \xi_c, \eta = \eta_c $ into Eqs. (97), we can rewrite R as

Define $ x = \frac{1}{r-r_c} $, we obtain

where $ c = 3r_c^2-\frac{a^2\eta_c}{r_c^2} $. The dependence between the coordinates r and θ is given by Eq. (23), and a detailed analysis of the θ motion is provided in the previous subsection. The solution for φ is obtained by numerical integration. Critical null geodesics with double root are shown in Fig. 7.

Figure 7.  Critical null orbits corresponding to double-coincident root in Kerr spacetime. The left diagram illustrates the trajectory of a null geodesic in the $ (r,\theta) $ plane, while the right diagram shows the trajectory of a null geodesic in three-dimensional space. The mass of the black hole is set to $ M = 1 $, and the initial conditions are $ r_{i} = 10,\mu_{i} = 0,\varphi_i = 0 $. The parameters for each row of images (from top to bottom) are $ a = 0.1,r_c = 2.9; a = 0.5,r_c = 2.5; $ and $ a = 0.9,r_c = 2.0 $. The red dashed line (ball) is located at $ r = r_c $, and the black dashed lines are at $ \mu = \pm \mu_{m} $.

For timelike geodesics with $ m\neq 0 $, we show in Appendix A that no triple-root unbound orbits exist when $ r_c \gt r_+ $. ² The conditions for a double root, $ R = 0 $ and $ R' = 0 $, yield

The function R can then be expressed as

Defining $ x = \frac{1}{r-r_c} $, one can evaluate the integrals to obtain

where

The numerical results are plotted in Fig. 8.

Figure 8.  Critical timelike geodesics in Kerr spacetime. The left diagram illustrates the trajectory of a timelike geodesic in the $ (r,\theta) $ plane, while the right diagram shows the trajectory in three-dimensional space. The mass of the black hole is set to $ M = 1 $, and the initial conditions are $ r_{i} = 10,\mu_{i} = 0,\varphi_i = 0 $. The parameters for each row of images (from top to bottom) are $ a = 0.1,r_c = 2.9,\chi = 1.02; a = 0.5,r_c = 2.5,\chi = 1.4; a = 0.9,r_c = 2.0,\chi = 4 $. The red dashed line (ball) is at $ r = r_c $, and the black dashed line denotes $ \mu = \pm \mu_{m} $.

In this section, we investigate the critical orbits of charged particles moving in a Kerr-Newman background. We define

and rewritten R as

The evaluation of $ \int\frac{d\theta}{\sqrt{\Theta}} $ has been addressed in detail in the preceding section; here we focus exclusively on solving for r.

In the case of null geodesics, with $ \tilde{e} = m = 0 $ and $ \chi = 1 $, the function R is further reduced to

Solving $ R = 0 $ and $ R' = 0 $, we obtains

As in the previous section, we compute the black hole shadow as seen by an observer in the equatorial plane. Figure 9 illustrates how the new parameter Q affects the black hole shadow.

Figure 9.  The shadow of a Kerr-Newman black hole as seen by a distant observer in the equatorial plane. The mass of the black hole is set to $ M = 1 $, and the spin parameter to $ a = 0.5 $.

For critical orbits with triply coincident roots $ R = 0,R' = 0,R'' = 0 $, we obtain

One can re-express $ r_c = M-(M(M^2-a^2-Q^2 ) )^{1/3} $. Since the photon sphere lies outside the event horizon, with $ r_+ \leq r_c $, we have $ r_c = r_+ = M $. It follows that $ \xi^{*}_c = \frac{M^2(3M^2-4Q^2)}{M^2-Q^2},\eta^{*}_c = \frac{2M^2-Q^2}{\sqrt{M^2-Q^2}},a^{*}_c = \sqrt{M^2-Q^2} $. Then, $ \frac{R}{E^2} $ reduces to the expression (108), and

where $ \mu_{1,2}^2 = \frac{2M\sqrt{3M^2-Q^2}\pm 3M^2}{M^2-Q^2} $. Thus, the charge Q affects the geodesic by altering $ \theta_m $. The integral $ \int \frac{dr}{\sqrt{R}} $ is given by Eq. (109), and $ \int\frac{d\mu}{\sqrt{\tilde{\Theta}}} $ is given by Eq. (91) for the case of $ m = 0 $. The numerical results for this case are shown in Fig. 10.

Figure 10.  Critical null geodesics with triply coincident roots in Kerr-Newman spacetime. The left diagram illustrates the trajectory of a null geodesic in the $ (r,\theta) $ plane, while the right diagram shows the trajectory in three-dimensional space. The mass of the black hole is set to $ M = 1 $, and the initial conditions are $ r_{i} = 10,\mu_{i} = 0,\varphi_i = 0 $. The parameters for each row of images (from top to bottom) are $ Q = 0.5, Q = 0.86, $ and $ Q = 0.866 $. The red dashed line (ball) is located at $ r = r_c $, and the black dashed lines are at $ \mu = \pm \mu_{m} $.

In the case of coincident double roots, one substitutes $ \xi = \xi_c,\eta = \eta_c $ into Eq. (124) to obtain

Defining $ x = \frac{1}{r-r_c} $, one can also obtain the expression (114) with a different value of c, namely $ c = 3r_c^2-\frac{a^2\eta_c}{r_c^2}-\frac{4Q^2\Delta(r_c)}{(r_c-M)^2} $. Figure 11 illustrates the solutions for critical orbits with a double root for different values of Q in the Kerr-Newman background.

Figure 11.  Critical null orbits corresponding to double coincident roots in Kerr-Newman spacetime are shown. The left diagram illustrates the trajectory of a null geodesic in the $ (r,\theta) $ plane, while the right diagram shows the trajectory of a null geodesic in three-dimensional space. The mass of the black hole is set to $ M = 1 $, and the initial conditions are $ r_{i} = 10,\mu_{i} = 0,\varphi_i = 0 $. The parameters for each row of images (from top to bottom) are $ Q = 0.3, Q = 0.5, $ and $ Q = 0.866 $ with $ a = 0.5, $ and $ r_c = 2.5 $. The red dashed line (ball) marks $ r = r_c $ and the black dashed line denotes $ \mu = \pm \mu_{m} $.

For timelike critical orbits, analytical calculations become increasingly challenging. Owing to the additional parameters, the expressions involved in such computations are highly complex. We begin by examining the case of a triple root. Solving the system $ R = 0 $, $ R' = 0 $, and $ R'' = 0 $ directly is difficult, so we can only solve each equation individually. We can rewrite $ R'' = 0 $ as

and obtains

Eliminating η in $ R' = 0 $, one obtains

Thus, we obtain

Eliminating ξ and η from the equation $ R = 0 $, yields a quartic algebraic equation in χ

where the coefficients $ a_0 $ – $ a_4 $ are listed in Appendix B. In principle, by solving equation (136), a complete set of parameter values $ \eta = \eta_c^{*} $, $ \xi = \xi_c^{*} $, and $ \chi = \chi_{c}^{*} $ can be obtained so that the equation $ R = 0 $ has a triple root. However, due to the complexity of these coefficients, it is practically impossible to derive an explicit expression for $ \chi = \chi_c^{*} $. It is also difficult to determine which parameter values satisfy the physical constraints.

Fortunately, we can adopt an alternative approach to analyze the triple-root scenario. In fact, we know that the function R can always be reduced to

By comparing with the expression (123), we find that $ r_1 = 2 \tilde{e} Q \chi_c^{*} -2 M \chi_c^{*}{} ^2+2 M-3 r_c $, where $ \chi = \chi_c^{*} $ satisfies Eq. (136). The integration yields

It is worth noting that, for this scenario to occur, the parameters must still satisfy reasonable conditions, including the following: (1) the existence of a black hole horizon; (2) $ r_c $ lying outside the horizon; (3) $ \chi>1 $; (4) $ r_1<r_c $; (5) $ \eta\geq 0 $; and so on. We have not found a fully consistent set of parameter values, although we cannot entirely rule out this possibility.

The case of double root is similar. The expressions for the solutions $ \eta = \eta_c $, $ \xi = \xi_c $ that satisfy the conditions $ R = 0 $ and $ R' = 0 $ are extremely complex, so we do not list them here. The function R then degenerates into

The integral $ \int\frac{dr}{\sqrt{E^2-m^2}\sqrt{R}} $ retains the form (119), albeit with different values of the parameters b and c:

The numerical results are shown in Fig. 12.

Figure 12.  Critical timelike orbits with coincident double roots in Kerr-Newman spacetime. The left diagram illustrates the trajectory of a timelike geodesic in the $ (r,\theta) $ plane, while the right diagram shows the corresponding trajectory in three-dimensional space. The initial conditions are $ r_{i} = 10,\mu_{i} = 0,\varphi_i = 0 $. The parameter sets for each row of images (from top to bottom) are $ a = 0.8, r_c = 1.8, Q = 0.5, \tilde{e} = 1.1, \chi = 1.1; a = 0.8, r_c = 2.1, Q = 0.3, \tilde{e} = 1.3, \chi = 1.3; a = 0.4, r_c = 2.5, Q = 0.6, \tilde{e} = 1.5, \chi = 1.5 $.

When a test particle moves from infinity toward a black hole, orbits that neither fall into the black hole nor are scattered by it are referred to as critical orbits. This corresponds to the equation $ R(r) = 0 $ for the radial kinetic energy having either a double or a triple real root. This paper investigates all possible cases of critical orbits, including their motion in four different black hole backgrounds and for three types of particles (massless, massive uncharged, and massive charged).

In spherically symmetric spacetimes, conservation of angular momentum confines orbits to a single plane. For critical orbits in the equatorial plane, the relation between the radius r and the angle φ can be expressed analytically. In axisymmetric spacetimes, however, the Carter constant supplements the conserved axial angular moentum, and generic orbits are three-dimensional. In the θ-direction, the orbit oscillates about the equatorial plane. We provide analytical expressions for $ (r, \theta, \varphi) $ and present the corresponding numerical results.

For null geodesics, the critical orbit lies at the boundary of accreting light rays. Consequently, its shape determines the appearance of the black hole as seen from infinity. For a Schwarzschild black hole, this boundary is a circle of radius $ 3\sqrt{3}M $. If the black hole is charged, this radius decreases. Moreover, as the spin parameter a increases, the shape becomes increasingly D-shaped.

In addition, we present a detailed discussion of the relationships between the physical parameters and the critical radius, which are crucial for modeling the accretion of collisionless Vlasov gas onto black holes. In future work, we will develop a more general theory of black hole accretion based on this study.

A triple root requires $ R = 0, R' = 0 $ and $ R'' = 0 $, which shows

Solving Eq. (144), one obtains

Inserting the solution into Eq. (143), one obtains

Then, the solution of $ \xi_c $ is given by

Inserting the expressions (145) and (147) into R and defining:

we can reexpress $ R = 0 $ as

Thus, we have

where $ B^2-AC = 16 M r_c^3 \Delta^3(r_c) $. Over the interval $ r_c\in[r_+,\infty) $, the function $ f(r_c) = B+\sqrt{B^2-AC} $ is monotonically decreasing, whereas $ C(r_c) $ is monotonically increasing. For $ a\in[0,M] $, one can show that $ f(r_+)\leq0 $ and $ C(r_+)\geq 0 $. Thus, for $ r_c> r_+ $, we have $ \chi_c^2<0 $, and no physically admissible parameter range exists.

Acknowledgement This work was supported in part by the National Natural Science Foundation of China (Grant No. 12505071) and by the Research Foundation of the Education Bureau of Hunan Province, China (Grant No. 25B0635). It was also supported in part by the Key Laboratory of Information Detection and Intelligent Processing Technology of the Hunan Provincial Department of Education, and by the Applied Characteristic Subject of Hunan Province, ”Electronic Science and Technology”.