© ESO, 2014

1. Introduction

In the vicinity of the black hole and any other compact objects, the gravitational field is extremely strong and the spacetime is significantly warped and twisted. Thus, the general relativity effects cannot be ignored. The motion of photons and test particles in this curved spacetime is along geodesics if we do not consider the external forces or perturbations exerting on them. The assumption that photons and particles propagate along geodesic trajectories is valid in most astrophysical contexts. The fast calculation of the null and time-like geodesics in curved spacetime is significantly important and has been widely used in astrophysical research (e.g., Cunningham & Bardeen 1973; Luminet 1979; Rauch & Blandford 1994; Hackmann 2010).

The calculation and applications of null geodesics in a curved spacetime, especially in a Kerr spacetime, have been discussed by many authors in different attempts to date (Dexter & Agol 2009; Hackmann 2010; Hackmann & Xu 2013; Chan et al. 2013; Yang & Wang 2013,and the references therein). To compute the geodesics it is possible to integrate a set of second-order differential equations in any relativistic spacetime directly, or to evaluate a set of elliptic integrals of motion in a Kerr-Newman (K-N) spacetime. In the present paper, we focus on the latter approach. There are four constants for any geodesic motions in a K-N spacetime (Carter 1968), which makes the reduction of the order of motion equations possible.

To obtain the optical appearance of a star orbiting around an extreme Kerr black hole, Cunningham & Bardeen (1973) calculated the null geodesics in a Kerr spacetime based on the elliptic integral method and proposed the impact parameters for the first time. After that, a method called ray-tracing was developed (e.g., Luminet 1979). Rauch & Blandford (1994) researched the optical caustics in a Kerr spacetime with an attempt to explain rapid X-ray variability in AGN. As a bonus they presented, in tabular form, cases that need to be considered for the calculation of both the null and time-like geodesics in a Kerr spacetime. Similar discussions and results are also given by Li et al. (2005) in their appendix.

The cases discussed by the above authors are very detailed and very complicated. As discussed in Yang & Wang (2013), this sophisticated situation can be significantly simplified by the introductions of Mino time p (Mino 2003) and the Weierstrass elliptic integrals and functions. Hackmann (2010) and Hackmann & Xu (2013) systematically discuss how to solve equations of geodesic motion, instead of restricting their research to the Kerr or K-N spacetime by Mino time and many types of elliptic functions. Carlson’s elliptic approach is quite suitable and efficient for evaluating elliptic integrals and functions, which has been demonstrated by Dexter & Agol (2009) and Yang & Wang (2013).

Motivated by the above discussions and the fact that there is currently no public code available to calculate time-like geodesics in a K-N spacetimes for all coordinates (including the proper times) at the same time, in this paper, we extend the scheme of Yang & Wang (2013) from null geodesics in a Kerr spacetime to time-like geodesics in a K-N spacetime. As a result, we developed a new public code, named ynogkm.

Analogous to ynogk, in ynogkm we also express the Boyer-Lidquist (B-L) coordinates r, θ, φ, t and the proper time σ as functions of the Mino time p semianalytically by using Weierstrass’ and Jacobi’s elliptic functions and integrals. Such treatment allows us to handle the practical applications conveniently. The Mino time p is an integral value along a particular geodesic. All of the elliptical integrals are computed using Carlson’s approach. Similarly to ynogk, we discuss how to compute the constants of motion from the initial conditions, i.e., the initial four-momentum of the particles measured under the local nonrotating frame (LNRF, Bardeen et al. 1972). For a massive particle with electric charge in a K-N spacetime, the number of constants of motion becomes 4. For a photon whose rest mass μ_m and electric charge ϵ are both zero, the number of constants of motion is 2. When taking μ_m and ϵ to be zero, the discussions here are reduced to those for null geodesics.

The paper is organized as follows. In Sect. 2, we provide the equations of motion for an electrically charged massive particle in a K-N spacetime. In Sect. 3, we discuss the expressions of the B-L coordinates and proper time as functions of parameter p. Then we reduce all of the elliptic integrals to standard forms that are evaluated by Carlson’s approach. Next, we discuss the calculation of constants of motion from initial conditions in Sect. 4. A brief introduction and discussion about the code are given in Sect. 5. In Sect. 6, we demonstrate the applications of our code to toy problems in the literature. Finally, a brief summary is presented in Sect. 7. Throughout this paper, the natural unit is used in which the constants G = c = 1. The mass of central black hole M is also taken to be 1, unless otherwise stated.

2. Equations of motion for time-like geodesics

We assume that the spin and electric charge of the black hole are a and e, respectively. Using the notation of Bardeen et al. (1972), we can write the K-N metric under the B-L coordinate as $\begin{array}{ccc}\mathrm{d}{\mathit{s}}^{\mathrm{2}}\mathrm{=}\mathrm{-}{\mathit{e}}^{\mathrm{2}\mathit{\nu }}\mathrm{d}{\mathit{t}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}\mathit{\psi }}\mathrm{(}\mathrm{d}\mathit{\phi }\mathrm{-}\mathit{\omega }\mathrm{d}\mathit{t}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}}\mathrm{d}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}}\mathrm{d}{\mathit{\theta }}^{\mathrm{2}}\mathit{,}& & \end{array}$(1)where $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{e}}^{\mathrm{2}\mathit{\nu }}\mathrm{=}\frac{\mathrm{\Sigma \Delta }}{\mathit{A}}\mathit{,} {\mathit{e}}^{\mathrm{2}\mathit{\psi }}\mathrm{=}\frac{{\sin }^{\mathrm{2}}\mathit{\theta A}}{\mathrm{\Sigma }}\mathit{,} {\mathit{e}}^{\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}}\mathrm{=}\frac{\mathrm{\Sigma }}{\mathrm{\Delta }}\mathit{,} \begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{e}}^{\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}}\mathrm{=}\mathrm{\Sigma }\mathit{,} \mathit{\omega }\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{r}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathit{a}}{\mathit{A}}\mathit{,}\end{array}$(2)and $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathrm{\Delta }\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{r}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}}\mathit{,} \mathrm{\Sigma }\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta ,} \begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathit{A}\mathrm{=}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{\Delta }{\mathit{a}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }\mathit{.}\end{array}$(3)Carter (1968) gave the first-order differential equations of motion for electrically charged massive particles as follows: $\begin{array}{ccc}\mathrm{\Sigma }\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{\lambda }}& \mathrm{=}& \mathrm{\pm }\sqrt{{\mathit{R}}_{\mathit{r}}}\mathit{,}\\ \mathrm{\Sigma }\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}\mathit{\lambda }}& \mathrm{=}& \mathrm{\pm }\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}\mathit{,}\\ \mathrm{\Sigma }\frac{\mathrm{d}\mathit{\phi }}{\mathrm{d}\mathit{\lambda }}& \mathrm{=}& \mathrm{-}\left(\mathit{aE}\mathrm{-}\frac{\mathit{L}}{{\sin }^{\mathrm{2}}\mathit{\theta }}\right)\mathrm{+}\frac{\mathit{aT}}{\mathrm{\Delta }}\mathit{,}\\ \mathrm{\Sigma }\frac{\mathrm{d}\mathit{t}}{\mathrm{d}\mathit{\lambda }}& \mathrm{=}& \mathrm{-}\mathit{a}\mathrm{(}\mathit{aE}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\mathit{L}\mathrm{)}\mathrm{+}\frac{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}\mathit{T}}{\mathrm{\Delta }}\mathit{,}\end{array}$where

$\begin{array}{ccc}& & \mathit{T}\mathrm{=}\mathit{E}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}\mathrm{-}\mathit{La}\mathrm{+}\mathit{eϵr,}\\ & & {\mathit{R}}_{\mathit{r}}\mathrm{=}{\mathit{T}}^{\mathrm{2}}\mathrm{-}\mathrm{\Delta }\mathrm{[}{\mathit{\mu }}_{\mathit{m}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{L}\mathrm{-}\mathit{aE}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{Q}}^{\mathrm{]}}\mathit{,}\\ & & {\mathrm{\Theta }}_{\mathit{\theta }}\mathrm{=}\mathit{Q}\mathrm{-}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{[}{\mathit{a}}^{\mathrm{2}}\mathrm{(}{\mathit{\mu }}_{\mathit{m}}^{\mathrm{2}}\mathrm{-}{\mathit{E}}^{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{L}}^{\mathrm{2}}\mathit{/}{\sin }^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{]}}\mathit{,}\end{array}$and λ = τ/μ_m, τ is the proper time, μ_m is the rest mass of the particle, Q is the Carter constant, E is the energy tested by an observer at infinity, L is the angular momentum of the particle about the black hole spin axis, and ϵ is the electric charge of the particle. From Eqs. (4)–(7) we can obtain the expression of the four-momentum for a particle $\begin{array}{ccc}{\mathit{p}}_{\mathit{\mu }}\mathrm{=}{\mathit{g}}_{\mathit{\mu \nu }}\frac{\mathrm{d}{\mathit{x}}^{\mathit{\nu }}}{\mathrm{d}\mathit{\lambda }}\mathrm{=}\left(\mathrm{-}\mathit{E,}\mathrm{\pm }\frac{\sqrt{{\mathit{R}}_{\mathit{r}}}}{\mathrm{\Delta }}\mathit{,}\mathrm{\pm }\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}\mathit{,L}\right)\mathit{.}& & \end{array}$(11)Equivalently, the equations of motion with integral forms can be written as:

$\begin{array}{ccc}& & \mathrm{\pm }{\mathrm{\int }}^{\mathit{\theta }}\frac{\mathrm{d}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathrm{=}\mathrm{\pm }{\mathrm{\int }}^{\mathit{r}}\frac{\mathrm{d}\mathit{r}}{\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathit{,}\\ & & \mathit{\sigma }\mathrm{=}\frac{\mathit{\tau E}}{{\mathit{\mu }}_{\mathit{m}}}\mathrm{=}{\mathrm{\int }}^{\mathit{\theta }}\frac{\mathit{E}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathrm{d}\mathit{\theta }\mathrm{+}{\mathrm{\int }}^{\mathit{r}}\frac{\mathit{E}{\mathit{r}}^{\mathrm{2}}}{\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathrm{d}\mathit{r,}\\ & & \mathit{t}\mathrm{=}\mathit{\sigma }\mathrm{+}\mathrm{2}{\mathrm{\int }}^{\mathit{r}}\frac{{\mathit{N}}_{\mathit{r}}}{\mathrm{\Delta }\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathrm{d}\mathit{r,}\\ & & \mathit{\phi }\mathrm{=}{\mathrm{\int }}^{\mathit{\theta }}\frac{\mathit{L}{\csc }^{\mathrm{2}}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathrm{d}\mathit{\theta }\mathrm{+}\mathit{a}{\mathrm{\int }}^{\mathit{r}}\frac{\mathit{r}\mathrm{(}\mathrm{2}\mathit{E}\mathrm{-}\mathit{eϵ}\mathrm{)}\mathrm{-}\mathrm{(}\mathit{E}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\mathit{L}{\mathit{a}}^{\mathrm{)}}}{\mathrm{\Delta }\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathrm{d}\mathit{r,}\end{array}$where $\begin{array}{ccc}{\mathit{N}}_{\mathit{r}}& \mathrm{=}& \mathrm{(}\mathrm{2}\mathit{E}\mathrm{+}\mathit{eϵ}\mathrm{)}{\mathit{r}}^{\mathrm{3}}\mathrm{-}\mathit{E}{\mathit{e}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{+}[\mathrm{2}\mathit{a}\mathrm{(}\mathit{Ea}\mathrm{-}\mathit{L}\mathrm{)}\\ & & \mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathit{eϵ}]\mathit{r}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathit{a}\mathrm{(}\mathit{Ea}\mathrm{-}\mathit{L}\mathrm{)}\mathit{.}\end{array}$(16)Here σ is a new variable, which is related to the proper time τ of the particle. For a photon, it becomes an affine parameter.

In many cases, we only need the equations of motion with integral forms. But in two special cases, i.e., the equatorial plane motion and the spherical motion, we need the differential equations of motion. In the case of the former, the particle is confined to the equatorial plane, one has Q = 0, θ ≡ π/2, and thus Θ_θ ≡ 0. Then the equations of motion with integral forms become invalid, since Θ_θ ≡ 0 appears in the denominator. But, from the differential equations, we can obtain the right equations to describe the plane motion. From Eq. (4), we have $\begin{array}{ccc}\mathit{\sigma }\mathrm{=}\mathit{\lambda E}\mathrm{=}{\mathrm{\int }}^{\mathit{r}}\frac{\mathit{E}{\mathit{r}}^{\mathrm{2}}}{\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathrm{d}\mathit{r}\mathit{.}& & \end{array}$(17)Dividing both sides of Eq. (7) by Eq. (4), we obtain $\begin{array}{ccc}\mathit{t}\mathrm{=}\mathit{\sigma }\mathrm{+}{\mathrm{\int }}^{\mathit{r}}\frac{{\mathit{N}}_{\mathit{r}}}{\mathrm{\Delta }\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathrm{d}\mathit{r}\mathit{.}& & \end{array}$(18)Similarly, from Eqs. (4) and (6), we obtain $\begin{array}{ccc}\mathit{\phi }\mathrm{=}\mathit{L}{\mathrm{\int }}^{\mathit{r}}\frac{\mathrm{d}\mathit{r}}{\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathrm{+}\mathit{a}{\mathrm{\int }}^{\mathit{r}}\frac{\mathit{r}\mathrm{(}\mathrm{2}\mathit{E}\mathrm{-}\mathit{eϵ}\mathrm{)}\mathrm{-}\mathrm{(}\mathit{E}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\mathit{L}{\mathit{a}}^{\mathrm{)}}}{\mathrm{\Delta }\sqrt{{\mathit{R}}_{\mathit{r}}}}\mathrm{d}\mathit{r}\mathit{.}& & \end{array}$(19)For spherical motion, we have R_r ≡ 0, thus the equations of motion with integral forms become invalid, because R_r appears in the denominator. Similarly, from Eq. (5) we obtain $\begin{array}{ccc}\mathit{\sigma }\mathrm{=}\mathit{E}\mathrm{\int }\mathrm{d}\mathit{\lambda }\mathrm{=}\mathit{E}{\mathit{a}}^{\mathrm{2}}{\mathrm{\int }}^{\mathit{\theta }}\frac{{\cos }^{\mathrm{2}}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathrm{d}\mathit{\theta }\mathrm{+}\mathit{E}{\mathit{r}}^{\mathrm{2}}{\mathrm{\int }}^{\mathit{\theta }}\frac{\mathrm{d}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathrm{·}& & \end{array}$(20)From Eqs. (5) and (7), we obtain $\begin{array}{ccc}\mathit{t}\mathrm{=}\mathit{\sigma }\mathrm{+}\frac{{\mathit{N}}_{\mathit{r}}}{\mathrm{\Delta }}{\mathrm{\int }}^{\mathit{\theta }}\frac{\mathrm{d}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathit{.}& & \end{array}$(21)And from Eqs. (5) and (6), we find $\begin{array}{ccc}\mathit{\phi }\mathrm{=}\mathit{L}{\mathrm{\int }}^{\mathit{\theta }}\frac{{\csc }^{\mathrm{2}}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathrm{d}\mathit{\theta }\mathrm{+}\mathit{a}\frac{\mathit{r}\mathrm{(}\mathrm{2}\mathit{E}\mathrm{-}\mathit{eϵ}\mathrm{)}\mathrm{-}\mathrm{(}\mathit{E}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\mathit{L}{\mathit{a}}^{\mathrm{)}}}{\mathrm{\Delta }}{\mathrm{\int }}^{\mathit{\theta }}\frac{\mathrm{d}\mathit{\theta }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}}\mathrm{·}& & \end{array}$(22)With Eqs. (12)–(22), we can calculate the geodesics by evaluating the elliptical integrals, instead of solving the differential equations of motion, and we can also express the B-L coordinates and proper time as functions of a parameter p. In the next section, we discuss how to get these functions semianalytically with elliptical functions and integrals.

3. Expressions of B-L coordinates and proper time as functions of p

3.1. The turning points

As discussed in Yang & Wang (2013), when we introduce a new parameter p with following definition from Eq. (12) $\begin{array}{ccc}\mathit{p}\mathrm{=}\mathrm{\pm }{\mathrm{\int }}^{\mathit{r}}\frac{\mathrm{d}\mathit{r}}{\sqrt{\mathit{R}}}\mathrm{=}\mathrm{\pm }{\mathrm{\int }}^{\mathit{\mu }}\frac{\mathrm{d}\mathit{\mu }}{\sqrt{{\mathrm{\Theta }}_{\mathit{\mu }}}}\mathit{,}& & \end{array}$(23)we can obtain functions r(p),μ(p),φ(p),t(p), and σ(p) from the equations of motion with integral forms, where μ = cosθ and

$\begin{array}{ccc}& & \mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{R}}_{\mathit{r}}}{{\mathit{E}}^{\mathrm{2}}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{m}}^{\mathrm{2}}\mathrm{)}{\mathit{r}}^{\mathrm{4}}\mathrm{+}\mathrm{2}\mathrm{(}{\mathit{m}}^{\mathrm{2}}\mathrm{+}\mathit{e}{\mathit{\epsilon }}^{\mathrm{)}}{\mathit{r}}^{\mathrm{3}}\mathrm{-}[\mathit{q}\mathrm{+}{\mathit{\lambda }}^{\mathrm{2}}\\ & & \mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{(}{\mathit{m}}^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}\mathrm{+}{\mathit{e}}^{\mathrm{2}}\mathrm{(}{\mathit{m}}^{\mathrm{2}}\mathrm{-}{\mathit{\epsilon }}^{\mathrm{2}}\mathrm{)}]{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{2}[\mathit{q}\mathrm{+}{\mathrm{(}\mathit{a}\mathrm{-}{\mathit{\lambda }}^{\mathrm{)}}}^{\mathrm{2}}\\ & & \mathrm{+}\mathit{ea\epsilon }\mathrm{(}\mathit{a}\mathrm{-}{\mathit{\lambda }}^{\mathrm{)}}]\mathit{r}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{(}\mathit{a}\mathrm{-}{\mathit{\lambda }}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathit{q,}\\ & & {\mathrm{\Theta }}_{\mathit{\mu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Theta }}_{\mathit{\theta }}{\sin }^{\mathrm{2}}\mathit{\theta }}{{\mathit{E}}^{\mathrm{2}}}\\ & & \hspace{0.17em}\mathrm{=}{\mathit{a}}^{\mathrm{2}}\mathrm{(}{\mathit{m}}^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}{\mathit{\mu }}^{\mathrm{4}}\mathrm{-}\mathrm{[}\mathit{q}\mathrm{+}{\mathit{\lambda }}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}{\mathrm{(}{\mathit{m}}^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{]}}{\mathit{\mu }}^{\mathrm{2}}\mathrm{+}\mathit{q}\mathit{.}\end{array}$Since the signs before the integrals are the same with dr and dθ, the parameter p monotonously increases along a particular geodesic. Here λ = L/E, q = Q/E², m = μ_m/E, and ε = ϵ/E, which are defined as constants of motion throughout this paper. Because R and Θ_μ are quartic (when m ≠ 1) or cubic (when m = 1) polynomials, the integrals about r and μ are elliptical integrals, which are reduced to the Weierstrass standard elliptic integrals (or those of Legendre’s when equation R(r) = 0 has no real roots).

Since both R and Θ_θ appear under the radical sign in the equations of motion, they must be nonnegative. The critical points satisfying R(r) = 0 or Θ_θ(θ) = 0 are so-called turning points in which the corresponding coordinate velocity is zero. When the motion of a particle is bounded in r or in θ coordinate, two turning points exist for the coordinate. We use r_tp1,r_tp2 and θ_tp1, θ_tp2 to denote the coordinates of these points, and assume that r_tp1 ≤ r_tp2, θ_tp1 ≤ θ_tp2. Since ${\mathit{p}}_{\mathit{r}}\mathrm{=}\mathrm{\pm }\sqrt{{\mathit{R}}_{\mathit{r}}}\mathit{/}\mathrm{\Delta }\mathit{,}{\mathit{p}}_{\mathit{\theta }}\mathrm{=}\mathrm{\pm }\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}$, when p_r = 0 or p_θ = 0 at the initial point, we have R_r(r_ini) = 0,Θ_θ(θ_ini) = 0, implying that the initial point is a turning point and r_ini (or θ_ini) is equal to one of r_tp1,r_tp2 (or θ_tp1, θ_tp2). Then we find r_ini ∈  [r_tp1,r_tp2] and θ_ini ∈  [θ_tp1,θ_tp2].

When r_tp2 does not exist at all (or equivalently, r_tp2 = ∞) and r_tp1 > r_h (r_h is the radius of the event horizon), the particle will eventually go to infinity far away. When r_tp1 does not exist (or r_tp1 < r_h) and r_tp2 > r_h, then the particle will eventually fall into the event horizon of the black hole. If (1) both r_tp1 and r_tp2 do not exist, this case equivalently corresponds to the fact that the equation R(r) = 0 has no real roots; or (2) r_tp2 does not exist and r_tp1 exists, but r_tp1 < r_h, for both cases the particle can move from infinity to the event horizon freely.

To obtain the θ coordinate of a turning point, we usually solve the equation Θ_μ(μ) = 0 to get μ_tp (=cosθ_tp) instead of solving the equation Θ_θ(θ) = 0. The roots of two equations are exactly the same except these special cases with constant λ = 0. The equation Θ_μ(μ) = 0 with λ = 0 has real roots ± 1, or 0,π, which are not the roots of equation Θ_θ(θ) = 0, indicating that a particle with λ = 0 can move from 0 to π freely and can go through the spin axis due to non-zero poloidal velocity p_θ =  ±$\sqrt{{\mathrm{\Theta }}_{\mathit{\theta }}}$ at the spin axis. Meanwhile, the particle changes the sign of its angular velocity dθ/dλ, and its azimuthal coordinate jumps from φ to φ ± π (Shakura 1987) instantaneously.

3.2. Coordinates μ and r

In this section, we express μ and r as functions of parameter p, i.e., μ = μ(p),r = r(p). The procedure to obtain these explicit expressions for electrically charged massive particles is quite tedious, but similar to the procedure for photons (refer to the discussions in Yang & Wang 2013). Thus, there is no need to present the details of the procedure. For reference, we present expressions of μ(p),r(p) in tabular form. See Tables 1 and 2.

For r, there are five cases:

• 1.

  m ≠ 1 and equation R(r) = 0 has at least one real root.

• 2.

  1 − m² > 0 (or | E |  > μ_m) and R(r) = 0 has no real roots.

• 3.

  1 − m² < 0 (or | E |  < μ_m) and R(r) = 0 has no real roots.

• 4.

  | m |  = 1 (or | E |  = μ_m) and the geodesic is unbounded.

• 5.

  | m |  = 1 (or | E |  = μ_m) and the geodesic is bounded.

In case 1, R(r) = 0 with m ≠ 1 has at least one real root. We are not concerned with the number of real roots in the equation or with their practical distribution. When the equation R(r) = 0 has real roots, turning point r_tp1 (or r_tp2) does exist, and can be easily determined with given r_ini¹.

In cases 2 and 3, R(r) = 0 has no real roots, but two pairs of complex conjugate roots are written as: $\begin{array}{ccc}{\mathit{r}}_{\mathrm{1}}& \mathrm{=}& \mathit{u}\mathrm{+}\mathrm{i}\mathit{v,} {\mathit{r}}_{\mathrm{2}}\mathrm{=}\mathit{u}\mathrm{-}\mathrm{i}\mathit{v,}\\ {\mathit{r}}_{\mathrm{3}}& \mathrm{=}& \mathit{w}\mathrm{+}\mathrm{i}\mathit{s,} {\mathit{r}}_{\mathrm{4}}\mathrm{=}\mathit{w}\mathrm{-}\mathrm{i}\mathit{s}\mathit{.}\end{array}$To avoid dealing with a complex integral, we use Jacobi’s elliptic function to express r instead of that of Weierstrass.

For cases 4 and 5, 1 − m² = 0, thus R(r) reduces to $\begin{array}{ccc}\mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \mathrm{2}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e\epsilon }\mathrm{)}{\mathit{r}}^{\mathrm{3}}\mathrm{-}\mathrm{[}\mathit{q}\mathrm{+}{\mathit{\lambda }}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\epsilon }}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{[}\mathit{q}\mathrm{+}\mathrm{(}\mathit{a}\mathrm{-}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}\\ & & \mathrm{+}\mathit{ea\epsilon }\mathrm{(}\mathit{a}\mathrm{-}\mathit{\lambda }\mathrm{)}\mathrm{]}\mathit{r}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{(}\mathit{a}\mathrm{-}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathit{q}\mathit{.}\end{array}$(41)From the expression of r(p) in Eq. (30), we know that when p + Π_r = ω′, ℘(p + Π_r;g₂,g₃) = ∞, i.e., r(p) = ∞, meaning that no matter what initial value p_r takes, the particle shall reach infinity eventually. As a result, the particular geodesic is unbounded.

3.3. Coordinates φ, t and the proper time σ

As discussed in Yang & Wang (2013), the expressions φ,t, and σ as functions of parameter p can be converted to evaluate the elliptic integrals appeared in the equations of motion with a given p. We divide the process into two steps. In the first step, the path and limits of the integrals are determined. In the second step, the integrals are reduced to standard forms, which are evaluated by Carlson’s approach.

3.3.1. The path and limits of integrals

For convenience, we use F_r and F_μ to represent the complicated integrands in integrals of r and μ, respectively.

The path is not monotonic when one or more than one turning points exist for r and μ. The path is divided into several parts in which each one has the maximum monotonic length, and the integrals are the sum of all individual parts. In Fig. 1, the integral path of r coordinate for a particular bounded geodesic is illustrated schematically. The motion of the particle is confined between two turning points, r_tp1 and r_tp2.

There are four important points in a particular path for r (or for μ), which are related to the integral limits. They are: 1. the initial position r_ini (or μ_ini); 2. the two turning points, r_tp1 and r_tp2 (or μ_tp1, μ_tp2); 3. the position corresponding to a given p, r_p, and μ_p. The z values of these points are z_ini, z_tp1, z_tp2, and z_p. The former three can be evaluated from z(r) functions listed in the right column of Table 3. The value of z_p can be evaluated from function z(p) = ℘(p + Π_r;g₂,g₃) for cases 1, 4, 5 and $\mathit{z}\mathrm{=}\mathrm{sn}\mathrm{(}\mathit{s}\sqrt{{\mathit{\lambda }}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{m}}^{\mathrm{2}}\mathrm{)}}\mathit{p}\mathrm{+}{\mathrm{\Pi }}_{\mathit{r}}\mathrm{|}{\mathit{k}}^{\mathrm{2}}\mathrm{)}$ for case 2, and $\mathit{z}\mathrm{=}\mathrm{sn}\mathrm{(}\mathit{s}\sqrt{{\mathit{\lambda }}_{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{m}}^{\mathrm{2}}\mathrm{)}}\mathit{p}\mathrm{+}{\mathrm{\Pi }}_{\mathit{r}}\mathrm{|}{\mathit{k}}^{\mathrm{2}}\mathrm{)}$ for case 3. For μ, the value of z_p can be computed from z = ℘(p + Π_μ;g₂,g₃). It is noted that the functions z(r) are monotonously decreasing, we have z_tp1 ≥ z_ini ≥ z_tp2 and z_tp1 ≥ z_p ≥ z_tp2. For μ, since the function z(μ) given in Table 1 is monotonously increasing, the two relationships are still valid.

In addition to z_p, for a given p, we also obtain the number of times that the particle meets the two turning points. We assume that the particle meets r_tp1 (or μ_tp1) for Nt₁ times, and r_tp2 (or μ_tp2) for Nt₂ times. If r_tp1 or r_tp2 does not exist, Nt₁ or Nt₂ is zero. To get Nt₁ and Nt₂ for a given p, we define the following five integrals with the help of Table 4: $\begin{array}{ccc}\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(42)where W(z) represents the integrands in Table 4. Apparently, we have p₁ ≥ 0 and p₂ ≥ 0, and $\begin{array}{ccc}{\mathit{p}}_{\mathrm{1}}\mathrm{=}{\mathit{I}}_{\mathrm{1}}\mathrm{-}{\mathit{p}}_{\mathrm{0}}\mathit{,} {\mathit{p}}_{\mathrm{2}}\mathrm{=}{\mathit{I}}_{\mathrm{2}}\mathrm{+}{\mathit{p}}_{\mathrm{0}}\mathit{.}& & \end{array}$(43)With the above definitions, we get the following identity $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{p}\mathrm{=}\mathrm{-}\mathrm{sign}\mathrm{\left(}{\mathit{p}}_{\mathit{\beta }}\mathrm{\right)}{\mathit{p}}_{\mathrm{0}}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{1}}{\mathit{p}}_{\mathrm{1}}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{2}}{\mathit{p}}_{\mathrm{2}}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\mathrm{-}\mathrm{\left[}\mathrm{sign}\mathrm{\right(}{\mathit{p}}_{\mathit{\beta }}\mathrm{)}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{1}}\mathrm{-}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{2}}\mathrm{]}{\mathit{p}}_{\mathrm{0}}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{1}}{\mathit{I}}_{\mathrm{1}}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{2}}{\mathit{I}}_{\mathrm{2}}\mathit{,}\end{array}& & \end{array}$(44)where β = r or θ, and p_β is the initial value of r or θ component of the four-momentum. We obtain Nt₁ and Nt₂ from the above equations by trial and error because Nt₁ and Nt₂ increase regularly as the particle moves, i.e., when p_β > 0 (or, p_β = 0 and z_ini = z_tp1), Nt₁ and Nt₂ increase as: $\begin{array}{ccc}\mathit{N}{\mathit{t}}_{\mathrm{1}}\mathrm{=}\mathrm{0}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3...}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& & \\ \mathit{N}{\mathit{t}}_{\mathrm{2}}\mathrm{=}\mathrm{0}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{4...}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& & \end{array}$When p_β < 0 (or, p_β = 0 and z_ini = z_tp2), Nt₁ and Nt₂ increase as: $\begin{array}{ccc}\mathit{N}{\mathit{t}}_{\mathrm{1}}\mathrm{=}\mathrm{0}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{4...}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& & \\ \mathit{N}{\mathit{t}}_{\mathrm{2}}\mathrm{=}\mathrm{0}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{3...}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& & \end{array}$Note the path and limits of integrals in t, φ, and σ are exactly the same with those of p. We introduce the following definitions: $\begin{array}{ccc}\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(45)where α = t, φ, σ. Similarly, we have $\begin{array}{ccc}{\mathit{I}}_{\mathrm{1}\mathit{\alpha \beta }}\mathrm{=}\mathit{I}{\mathit{T}}_{\mathrm{1}\mathit{\alpha \beta }}\mathrm{-}{\mathit{I}}_{\mathrm{0}\mathit{\alpha \beta }}\mathit{,} {\mathit{I}}_{\mathrm{2}\mathit{\alpha \beta }}\mathrm{=}\mathit{I}{\mathit{T}}_{\mathrm{2}\mathit{\alpha \beta }}\mathrm{+}{\mathit{I}}_{\mathrm{0}\mathit{\alpha \beta }}\mathit{.}& & \end{array}$(46)Then the integrals in t, φ, and σ can be written as $\begin{array}{ccc}{\mathit{I}}_{\mathit{\alpha \beta }}& \mathrm{=}& \mathrm{-}\mathrm{sign}\mathrm{\left(}{\mathit{p}}_{\mathit{\beta }}\mathrm{\right)}{\mathit{I}}_{\mathrm{0}\mathit{\alpha \beta }}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{1}}{\mathit{I}}_{\mathrm{1}\mathit{\alpha \beta }}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{2}}{\mathit{I}}_{\mathrm{2}\mathit{\alpha \beta }}\mathit{,}\\ & \mathrm{=}& \mathrm{-}\mathrm{\left[}\mathrm{sign}\mathrm{\right(}{\mathit{p}}_{\mathit{\beta }}\mathrm{)}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{1}}\mathrm{-}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{2}}\mathrm{]}{\mathit{I}}_{\mathrm{0}\mathit{\alpha \beta }}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{1}}\mathit{I}{\mathit{T}}_{\mathrm{1}\mathit{\alpha \beta }}\mathrm{+}\mathrm{2}\mathit{N}{\mathit{t}}_{\mathrm{2}}\mathit{I}{\mathit{T}}_{\mathrm{2}\mathit{\alpha \beta }}\mathit{.}\end{array}$(47)Finally, we have $\begin{array}{ccc}\mathit{t}\mathrm{=}{\mathit{I}}_{\mathit{tr}}\mathrm{+}{\mathit{I}}_{\mathit{t\theta }}\mathit{,} \mathit{\phi }\mathrm{=}{\mathit{I}}_{\mathit{\phi r}}\mathrm{+}{\mathit{I}}_{\mathit{\phi \theta }}\mathit{,} \mathit{\sigma }\mathrm{=}{\mathit{I}}_{\mathit{\sigma t}}\mathrm{+}{\mathit{I}}_{\mathit{\sigma \theta }}\mathit{.}& & \end{array}$(48)

3.3.2. The computation of elliptic integrals by Carlson’s approach

In this section, we discuss how to compute the elliptic integrals appeared in t, φ, and σ using Carlson’s approach. First, we reduce these integrals to the standard forms. Before the reductions we introduce two notations J_k(h) and I_k(h) with the following definition: $\begin{array}{ccc}{\mathit{J}}_{\mathit{k}}\mathrm{\left(}\mathit{h}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\\ \mathit{y}\end{array}\frac{\mathrm{d}\mathit{t}}{\mathrm{(}\mathit{t}\mathrm{-}\mathit{h}{\mathrm{)}}^{\mathit{k}}\sqrt{\mathrm{4}{\mathit{t}}^{\mathrm{3}}\mathrm{-}{\mathit{g}}_{\mathrm{2}}\mathit{t}\mathrm{-}{\mathit{g}}_{\mathrm{3}}}}\mathit{,}\\ {\mathit{I}}_{\mathit{k}}\mathrm{\left(}\mathit{h}\mathrm{\right)}& \mathrm{=}& {\mathrm{\int }}_{\mathit{y}}^{\mathit{x}}\frac{\mathrm{d}\mathit{r}}{\mathrm{(}\mathit{r}\mathrm{-}\mathit{h}{\mathrm{)}}^{\mathit{k}}\sqrt{\mathrm{\left[}\mathrm{\right(}\mathit{r}\mathrm{-}\mathit{u}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}^{\mathrm{2}}\mathrm{\left]}\mathrm{\right[}\mathrm{(}\mathit{r}\mathrm{-}\mathit{w}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{s}}^{\mathrm{2}}\mathrm{]}}}\mathit{,}\end{array}$where k =  −2, − 1,0,1,2. From Eqs. (13)–(15), we have

$\begin{array}{ccc}& & {\mathit{\sigma }}_{\mathit{\mu }}\mathrm{=}{\mathit{t}}_{\mathit{\mu }}\mathrm{=}{\mathit{a}}^{\mathrm{2}}\left[\frac{{\mathit{b}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{16}}{\mathit{J}}_{\mathrm{2}}\left(\frac{{\mathit{b}}_{\mathrm{1}}}{\mathrm{4}}\right)\mathrm{+}\frac{{\mathit{b}}_{\mathrm{0}}{\mathit{\mu }}_{{\mathrm{tp}}_{\mathrm{1}}}}{\mathrm{2}}{\mathit{J}}_{\mathrm{1}}\left(\frac{{\mathit{b}}_{\mathrm{1}}}{\mathrm{4}}\right)\mathrm{+}{\mathit{\mu }}_{{\mathrm{tp}}_{\mathrm{1}}}^{\mathrm{2}}\mathit{p}\right]\mathit{,}\\ & & {\mathit{\phi }}_{\mathit{\mu }}\mathrm{=}\mathit{\lambda }\left[\frac{\mathit{p}}{\mathrm{1}\mathrm{-}{\mathit{\mu }}_{{\mathrm{tp}}_{\mathrm{1}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{C}}_{\mathrm{-}}}{\mathrm{2}}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}{\mathit{z}}_{\mathrm{-}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{C}}_{\mathrm{+}}}{\mathrm{2}}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}{\mathit{z}}_{\mathrm{+}}\mathrm{\right)}\right]\mathit{,}\end{array}$where $\begin{array}{ccc}& & {\mathit{C}}_{\mathrm{\pm }}\mathrm{=}\frac{{\mathit{b}}_{\mathrm{0}}}{\mathrm{4}\mathrm{(}\mathrm{\pm }\mathrm{1}\mathrm{-}{\mathit{\mu }}_{{\mathrm{tp}}_{\mathrm{1}}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}\\ & & {\mathit{z}}_{\mathrm{\pm }}\mathrm{=}\frac{{\mathit{b}}_{\mathrm{0}}}{\mathrm{4}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{\pm }\mathrm{1}\mathrm{-}{\mathit{\mu }}_{{\mathrm{tp}}_{\mathrm{1}}}\mathrm{)}}\mathrm{+}\frac{{\mathit{b}}_{\mathrm{1}}}{\mathrm{4}}\mathrm{·}\end{array}$Noting the definition of parameter p, we have replaced J₀ with p in the above equations.

The reduced standard forms of integrals for r are given in Tables 5 and 6. The standard forms for the special cases, such as the equatorial plane motion and the spherical motion, can be obtained directly and are not given here anymore.

Carlson (1988, 1989, 1991, 1992) developed a new approach to compute elliptic integrals (Press et al. 2007). He gave new definitions of the standard elliptic integrals of the first and third kind

$\begin{array}{ccc}& & {\mathit{R}}_{\mathit{F}}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{t}}{\sqrt{\mathrm{(}\mathit{t}\mathrm{+}\mathit{x}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{+}\mathit{y}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{+}\mathit{z}\mathrm{)}}}\mathit{,}\\ & & {\mathit{R}}_{\mathit{J}}\mathrm{\left(}\mathit{x,y,z,p}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{t}}{\mathrm{(}\mathit{t}\mathrm{+}\mathit{p}\mathrm{)}\sqrt{\mathrm{(}\mathit{t}\mathrm{+}\mathit{x}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{+}\mathit{y}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{+}\mathit{z}\mathrm{)}}}\mathit{,}\end{array}$and the degenerate cases of R_C(x,y) = R_F(x,y,y) and R_D(x,y,z) = R_J(x,y,z,z). Case R_D can be regarded as the standard elliptic integral of the second kind. Carlson denotes the elliptic integrals with a symbol with the following definition: $\begin{array}{ccc}\mathrm{\left[}{\mathit{p}}_{\mathrm{1}}\mathit{,}\mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}\mathit{,}{\mathit{p}}_{\mathit{k}}\mathrm{\right]}\mathrm{=}{\mathrm{\int }}_{\mathit{y}}^{\mathit{x}}\underset{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}{\overset{\mathit{k}}{\mathrm{􏽙}}}\mathrm{(}{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{b}}_{\mathit{i}}\mathit{t}{\mathrm{)}}^{{\mathit{p}}_{\mathit{i}}\mathit{/}\mathrm{2}}\mathrm{d}\mathit{t}\mathit{.}& & \end{array}$(66)If a_i + b_it is complex, then its complex conjugate $\overline{{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{b}}_{\mathit{i}}\mathit{t}}$ must exist and guarantee the integral is real. And, $\mathrm{(}{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{b}}_{\mathit{i}}\mathit{t}\mathrm{)}\mathrm{\left(}\overline{{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{b}}_{\mathit{i}}\mathit{t}}\mathrm{\right)}\mathrm{=}\mathit{f}\mathrm{+}\mathit{gt}\mathrm{+}\mathit{h}{\mathit{t}}^{\mathrm{2}}$. Thus, $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathrm{\left[}{\mathit{p}}_{\mathrm{1}}\mathit{,}{\mathit{p}}_{\mathrm{1}}\mathit{,}{\mathit{p}}_{\mathrm{3}}\mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}\mathit{,}{\mathit{p}}_{\mathit{k}}\mathrm{\right]}\mathrm{=}\hspace{0.17em}{\mathrm{\int }}_{\mathit{y}}^{\mathit{x}}\mathrm{(}\mathit{f}\mathrm{+}\mathit{gt}\mathrm{+}\mathit{h}{\mathit{t}}^{\mathrm{2}}{\mathrm{)}}^{{\mathit{p}}_{\mathrm{1}}\mathit{/}\mathrm{2}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{\times }\underset{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{3}}{\overset{\mathit{k}}{\mathrm{􏽙}}}\mathrm{(}{\mathit{a}}_{\mathit{j}}\mathrm{+}{\mathit{b}}_{\mathit{j}}\mathit{t}{\mathrm{)}}^{{\mathit{p}}_{\mathit{j}}\mathit{/}\mathrm{2}}\mathrm{d}\mathit{t}\mathit{.}\end{array}& & \end{array}$(67)For a particular elliptic integral, there is an unique formula to evaluate it. We give a simple example here: $\begin{array}{ccc}\mathrm{[}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathrm{]}\mathrm{=}\mathrm{2}{\mathit{R}}_{\mathit{F}}\mathrm{\left(}{\mathit{U}}_{\mathrm{12}}^{\mathrm{2}}\mathit{,}{\mathit{U}}_{\mathrm{13}}^{\mathrm{2}}\mathit{,}{\mathit{U}}_{\mathrm{14}}^{\mathrm{2}}\mathrm{\right)}\mathit{,}& & \end{array}$(68)where

$\begin{array}{ccc}& & {\mathit{U}}_{\mathit{ij}}\mathrm{=}\mathrm{(}{\mathit{X}}_{\mathit{i}}{\mathit{X}}_{\mathit{j}}{\mathit{Y}}_{\mathit{k}}{\mathit{Y}}_{\mathit{m}}\mathrm{+}{\mathit{Y}}_{\mathit{i}}{\mathit{Y}}_{\mathit{j}}{\mathit{X}}_{\mathit{k}}{\mathit{X}}_{\mathit{m}}\mathrm{)}\mathit{/}\mathrm{(}\mathit{x}\mathrm{-}\mathit{y}\mathrm{)}\mathit{,}\\ & & {\mathit{X}}_{\mathit{i}}\mathrm{=}\sqrt{{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{b}}_{\mathit{i}}\mathit{x}}\mathit{,} {\mathit{Y}}_{\mathit{i}}\mathrm{=}\sqrt{{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{b}}_{\mathit{i}}\mathit{y}}\mathit{.}\end{array}$The elliptic integrals we need to evaluate in this paper are: J₁,J₂,I₂,I_-1,I_-2, which can be recast with Carlson’s notations. When equation 4t³ − g₂t − g₃ = 0 has three real roots denoted by e₁,e₂,e₃, we obtain (Carlson 1988) $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{J}}_{\mathit{k}}\mathrm{\left(}\mathit{h}\mathrm{\right)}\mathrm{=}{\mathit{s}}_{\mathit{h}}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathit{y}}^{\mathit{x}}\frac{\mathrm{d}\mathit{t}}{\sqrt{\mathrm{(}\mathit{t}\mathrm{-}{\mathit{e}}_{\mathrm{1}}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{e}}_{\mathrm{2}}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{e}}_{\mathrm{3}}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{-}\mathit{h}{\mathrm{)}}^{\mathrm{2}\mathit{k}}}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}{\mathit{s}}_{\mathit{h}}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{[}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{2}\mathit{k}\mathrm{]}\mathit{,}\end{array}& & \end{array}$(69)where s_h = sign [(y − h)^k]. When equation 4t³ − g₂t − g₃ = 0 has one pair of complex roots and one real root e₁, we obtain (Carlson 1991) $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{J}}_{\mathit{k}}\mathrm{\left(}\mathit{h}\mathrm{\right)}\mathrm{=}{\mathit{s}}_{\mathit{h}}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathit{y}}^{\mathit{x}}\frac{\mathrm{d}\mathit{t}}{\sqrt{\mathrm{(}\mathit{t}\mathrm{-}{\mathit{e}}_{\mathrm{1}}\mathrm{)}\mathrm{(}{\mathit{t}}^{\mathrm{2}}\mathrm{+}\mathit{gt}\mathrm{+}\mathit{f}\mathrm{)}\mathrm{(}\mathit{t}\mathrm{-}\mathit{h}{\mathrm{)}}^{\mathrm{2}\mathit{k}}}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}{\mathit{s}}_{\mathit{h}}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{[}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{2}\mathit{k}\mathrm{]}\mathit{.}\end{array}& & \end{array}$(70)The integral I_k(h) corresponds to the case that equation R_r(r) = 0 has no real roots and can be expressed as (Carlson 1992): $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{I}}_{\mathit{k}}\mathrm{\left(}\mathit{h}\mathrm{\right)}\mathrm{=}{\mathit{s}}_{\mathit{h}}{\mathrm{\int }}_{\mathit{y}}^{\mathit{x}}\frac{\mathrm{d}\mathit{r}}{\sqrt{\mathrm{\left[}\mathrm{\right(}\mathit{r}\mathrm{-}\mathit{u}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}^{\mathrm{2}}\mathrm{\left]}\mathrm{\right[}\mathrm{(}\mathit{r}\mathrm{-}\mathit{w}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{s}}^{\mathrm{2}}\mathrm{]}\mathrm{(}\mathit{r}\mathrm{-}\mathit{h}{\mathrm{)}}^{\mathrm{2}\mathit{k}}}}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}{\mathit{s}}_{\mathit{h}}\mathrm{[}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{-}\mathrm{2}\mathit{k}\mathrm{]}\mathit{.}\end{array}& & \end{array}$(71)Up to now, we have expressed all coordinates and proper time as functions of parameter p semi-analytically. As discussed in Yang & Wang (2013), such treatment is very convenient for massive particles whose geodesics can be bounded, and the number of times that the particle meets the turning points can be arbitrary both for r and μ coordinates. In addition to p, one needs to prescribe the constants of motion. In the next section, we discuss how to get them from the initial four-momentum of a particle.

4. Constants of motion

As mentioned above, the constants of motion throughout this paper are defined as $\begin{array}{ccc}\mathit{\lambda }\mathrm{=}\frac{\mathit{L}}{\mathit{E}}\mathit{,q}\mathrm{=}\frac{\mathit{Q}}{{\mathit{E}}^{\mathrm{2}}}\mathit{,m}\mathrm{=}\frac{{\mathit{\mu }}_{\mathit{m}}}{\mathit{E}}\mathit{,\epsilon }\mathrm{=}\frac{\mathit{ϵ}}{\mathit{E}}\mathit{,}& & \end{array}$(72)which can be obtained from the initial four-momentum of a particle given in an LNRF reference. But to handle more complicated applications, we want to specify the initial four-momentum in the reference of an assumed emitter, instead of an LNRF reference directly. However, the initial four-momentum is finally transformed into an LNRF reference by a Lorentz transformation.

Now we introduce the LNRF reference, which is also called zero angular momentum observers (ZAMO, Bardeen et al. 1972). The orthonormal tetrad is given by ${{e}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\left(}\mathrm{LNRF}\mathrm{\right)}\mathrm{=}{\mathit{e}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathit{\nu }}{\mathit{\partial }}_{\mathit{\nu }}\mathit{,}$(73)where $\begin{array}{ccc}{\mathit{e}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathit{\nu }}\mathrm{=}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{,}& & \end{array}$(74)and the dual form of which is ${{e}}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\left(}\mathrm{LNRF}\mathrm{\right)}\mathrm{=}{\mathit{e}}_{\mathit{\nu }}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{d}{\mathit{x}}^{\mathit{\nu }}\mathit{,}$(75)where $\begin{array}{ccc}{\mathit{e}}_{\mathit{\nu }}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{=}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{.}& & \end{array}$(76)We assume that the particle is emitted by an emitter at the initial position, where the emitter has coordinate velocities ṙ = dr/dt, , and , then its physical velocities υ_r,υ_θ,υ_φ, with respect to the LNRF fixed at the same point, can be written as (Bardeen et al. 1972): $\begin{array}{ccc}{\mathit{\upsilon }}_{\mathit{r}}\mathrm{=}{\mathit{e}}^{{\mathit{\mu }}_{\mathrm{1}}\mathrm{-}\mathit{\nu }}\mathit{ṙ}\mathit{,} {\mathit{\upsilon }}_{\mathit{\theta }}\mathrm{=}{\mathit{e}}^{{\mathit{\mu }}_{\mathrm{2}}\mathrm{-}\mathit{\nu }}\mathit{\theta ̇}\mathit{,} {\mathit{\upsilon }}_{\mathit{\phi }}\mathrm{=}{\mathit{e}}^{\mathit{\psi }\mathrm{-}\mathit{\nu }}\mathrm{(}\mathrm{\Omega }\mathrm{-}\mathit{\omega }\mathrm{)}\mathit{.}& & \end{array}$(77)The orthonormal tetrad of the emitter can be obtained by rotating the tetrad of the LNRF reference in the local four dimension spacetime. The rotation is simply a Lorentz transformation. We denote the matrix of the rotation by Λ, and find $\begin{array}{ccc}{{e}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\left(}\mathrm{em}\mathrm{\right)}\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathrm{\left(}\mathit{b}\mathrm{\right)}}{{e}}_{\mathrm{\left(}\mathit{b}\mathrm{\right)}}\mathrm{\left(}\mathrm{LNRF}\mathrm{\right)}\mathit{,}& & \end{array}$(78)where (Misner et al. 1973) $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathrm{\left(}\mathit{b}\mathrm{\right)}}\mathrm{=}(\begin{array}{c}\\ & & \\ & & \\ & & \\ & & \end{array}& & \\ \begin{array}{c}\\ & \\ & \\ & \\ & \end{array})\mathit{,}& & \end{array}$(79)where $\mathit{\gamma }\mathrm{=}\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{)}{\mathrm{]}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$, and the covariant tetrad of the emitter is ${{e}}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\left(}\mathrm{em}\mathrm{\right)}\mathrm{=}\mathrm{\left(}{\mathrm{\Lambda }}^{-1}{\mathrm{)}}_{\mathrm{\left(}\mathit{b}\mathrm{\right)}}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{e}}^{\mathrm{\left(}\mathit{b}\mathrm{\right)}}\mathrm{\right(}\mathrm{LNRF}\mathrm{)}\mathit{,}$(80)where $\begin{array}{ccc}\mathrm{(}{\mathrm{\Lambda }}^{-1}{\mathrm{)}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathrm{\left(}\mathit{b}\mathrm{\right)}}\mathrm{=}(\begin{array}{c}\\ & & \\ & & \\ & & \\ & & \end{array}& & \\ \begin{array}{c}\\ & \\ & \\ & \\ & \end{array})\mathit{.}& & \end{array}$(81)Equivalently, we have

$\begin{array}{ccc}& & \\ & & {{e}}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\left(}\mathrm{LNRF}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{\Lambda }{\mathrm{)}}_{\mathrm{\left(}\mathit{b}\mathrm{\right)}}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{e}}^{\mathrm{\left(}\mathit{b}\mathrm{\right)}}\mathrm{\right(}\mathrm{em}\mathrm{)}\mathit{.}\end{array}$In the rest frame of the emitter, the components of the four-momentum of the particle are denoted by ${\mathit{p}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathrm{\prime }}$, which can be regarded as the projections of the momentum p on the corresponding basis vectors, i.e., ${\mathit{p}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathrm{\prime }}\mathrm{=}{p}\mathrm{·}{{e}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\left(}\mathrm{em}\mathrm{\right)}\mathit{.}$(84)Multiplying both sides of Eq. (82) by p, we find $\begin{array}{ccc}{\mathit{e}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathit{\nu }}{\mathit{p}}_{\mathit{\nu }}\mathrm{=}\mathrm{(}{\mathrm{\Lambda }}^{-1}{\mathrm{)}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathrm{\left(}\mathit{b}\mathrm{\right)}}{\mathit{p}}_{\mathrm{\left(}\mathit{b}\mathrm{\right)}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(85)By assuming that the physical velocities of the particle with respect to the emitter are: ${\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}$, we find $\begin{array}{ccc}{\mathit{p}}_{\mathrm{\left(}\mathit{a}\mathrm{\right)}}^{\mathrm{\prime }}\mathrm{=}{\mathit{\gamma }}^{\mathrm{\prime }}{\mathit{\mu }}_{\mathit{m}}\mathrm{(}\mathrm{-}\mathrm{1}\mathit{,}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{)}\mathit{,}& & \end{array}$(86)where ${\mathit{\gamma }}^{\mathrm{\prime }}\mathrm{=}\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }\mathrm{2}}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }\mathrm{2}}\mathrm{)}{\mathrm{]}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$ is the Lorentz factor. Equation (85) can be expanded explicitly with Eqs. (11) and (74):

$\begin{array}{ccc}& & \\ & & \\ & & \\ & & \end{array}$where $\begin{array}{ccc}{\mathit{k}}_{\mathrm{\left(}\mathit{t}\mathrm{\right)}}& \mathrm{=}& \mathit{\gamma }\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\upsilon }}_{\mathit{r}}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\theta }}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\phi }}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{)}\mathit{,}\\ {\mathit{k}}_{\mathrm{\left(}\mathit{r}\mathrm{\right)}}& \mathrm{=}& \mathit{\gamma }{\mathit{\upsilon }}_{\mathit{r}}\mathrm{+}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\upsilon }}_{\mathit{r}}}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{(}{\mathit{\upsilon }}_{\mathit{r}}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\theta }}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\phi }}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{)}\mathit{,}\\ {\mathit{k}}_{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}& \mathrm{=}& \mathit{\gamma }{\mathit{\upsilon }}_{\mathit{\theta }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\upsilon }}_{\mathit{\theta }}}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{(}{\mathit{\upsilon }}_{\mathit{r}}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\theta }}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\phi }}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{)}\mathit{,}\\ {\mathit{k}}_{\mathrm{\left(}\mathit{\phi }\mathrm{\right)}}& \mathrm{=}& \mathit{\gamma }{\mathit{\upsilon }}_{\mathit{\phi }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\upsilon }}_{\mathit{\phi }}}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{(}{\mathit{\upsilon }}_{\mathit{r}}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\theta }}{\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{+}{\mathit{\upsilon }}_{\mathit{\phi }}{\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{)}\mathit{.}\end{array}$If we solve Eqs. (87) and (90) for λ, we obtain $\begin{array}{ccc}\mathit{\lambda }\mathrm{=}\frac{{\mathit{k}}_{\mathrm{\left(}\mathit{\phi }\mathrm{\right)}}\sin \mathit{\theta }}{{\mathit{k}}_{\mathrm{\left(}\mathit{t}\mathrm{\right)}}\sqrt{\mathrm{\Delta }}\mathrm{\Sigma }\mathit{/}\mathit{A}\mathrm{+}{\mathit{k}}_{\mathrm{\left(}\mathit{\phi }\mathrm{\right)}}\mathit{\omega }\sin \mathit{\theta }}\mathrm{·}& & \end{array}$(95)With λ, from Eq. (87), we find $\begin{array}{ccc}\mathit{m}\mathrm{=}\frac{{\mathit{e}}^{\mathrm{-}\mathit{\nu }}}{{\mathit{\gamma }}^{\mathrm{\prime }}{\mathit{k}}_{\mathrm{\left(}\mathit{t}\mathrm{\right)}}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\lambda \omega }\mathrm{)}\mathit{.}& & \end{array}$(96)With λ and m, from Eq. (89), we find $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{q}\mathrm{=}\hspace{0.17em}\left[{\left(\frac{{\mathit{k}}_{\mathrm{\left(}\mathit{\phi }\mathrm{\right)}}}{{\mathit{k}}_{\mathrm{\left(}\mathit{t}\mathrm{\right)}}\sqrt{\mathrm{\Delta }}\mathrm{\Sigma }\mathit{/}\mathit{A}\mathrm{+}{\mathit{k}}_{\mathrm{\left(}\mathit{\phi }\mathrm{\right)}}\mathit{\omega }\sin \mathit{\theta }}\right)}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{(}{\mathit{m}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}\right]{\cos }^{\mathrm{2}}\mathit{\theta }\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{+}{\mathit{m}}^{\mathrm{2}}\mathrm{\Sigma }{\mathit{\gamma }}^{\mathrm{\prime }\mathrm{2}}{\mathit{k}}_{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}^{\mathrm{2}}\mathit{.}\end{array}& & \end{array}$(97)With λ, m, and q, from Eq. (88), we obtain a quadratic equation for ε, $\begin{array}{ccc}{\mathit{a}}_{\mathrm{1}}{\mathit{\epsilon }}^{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{1}}\mathit{\epsilon }\mathrm{+}{\mathit{c}}_{\mathrm{1}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(98)where $\begin{array}{ccc}{\mathit{a}}_{\mathrm{1}}& \mathrm{=}& {\mathit{e}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathit{,}\\ {\mathit{b}}_{\mathrm{1}}& \mathrm{=}& \mathrm{2}\mathit{er}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{-}\mathit{a\lambda }\mathrm{)}\mathit{,}\\ {\mathit{c}}_{\mathrm{1}}& \mathrm{=}& \mathrm{(}\mathrm{1}\mathrm{-}{\mathit{m}}^{\mathrm{2}}\mathrm{)}{\mathit{r}}^{\mathrm{4}}\mathrm{+}\mathrm{2}{\mathit{m}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{3}}\mathrm{-}\mathrm{[}\mathit{q}\mathrm{+}{\mathit{\lambda }}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{(}{\mathit{m}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{e}}^{\mathrm{2}}{\mathit{m}}^{\mathrm{2}}\mathrm{]}{\mathit{r}}^{\mathrm{2}}\\ & & \mathrm{+}\mathrm{2}\mathrm{[}\mathit{q}\mathrm{+}\mathrm{(}\mathit{a}\mathrm{-}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}\mathrm{]}\mathit{r}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{(}\mathit{a}\mathrm{-}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathit{q}\\ & & \mathrm{-}\mathrm{\Delta \Sigma }\mathrm{(}\mathit{m}{\mathit{\gamma }}^{\mathrm{\prime }}{\mathit{k}}_{\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{)}}^{\mathrm{2}}\mathit{.}\end{array}$We then find two values, $\begin{array}{ccc}{\mathit{\epsilon }}_{\mathrm{\pm }}\mathrm{=}\frac{\mathrm{-}{\mathit{b}}_{\mathrm{1}}\mathrm{\pm }\sqrt{{\mathit{b}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}\mathrm{4}{\mathit{a}}_{\mathrm{1}}{\mathit{c}}_{\mathrm{1}}}}{\mathrm{2}{\mathit{a}}_{\mathrm{1}}}\mathit{,}& & \end{array}$(102)where ε₊ represents positive electric charge, and ε₋ represents negative electric charge. In Fig. 2, we plot a set of geodesic orbits of particles emitted isotropically in an LNRF reference. These particles are confined in the equatorial plane. The particles have ε₊ in the top panel, and ε₋ in the bottom panel.

5. A brief introduction to the code

According to the discussions above, we have developed a new public code for computing null and time-like geodesics in a K-N spacetime². We name the code ynogkm, which is written in fortran 95, and the object-oriented method is used. The code consists of several independent modules, in which each one completes a special goal. The most important two modules are ellfunction and blcoordinates. The former contains the supporting functions and subroutines computing elliptic integrals by Carlson’s approach. The latter contains the functions and routines computing the B-L coordinate functions: r(p), μ(p), φ(p), and t(p), as well as proper time function σ(p). In blcoordinates, we provide a subroutine named ynogkm to compute all coordinates and proper times simultaneously for a given p. We also provide two functions named radius and mucos to compute r(p) and μ(p), respectively. In an axis-symmetry case, we only need to compute r and μ.

Before calling these functions and subroutines to compute the B-L coordinates, we need to provide the constants of motion, namely, λ, q, m, and ε. As discussed in the above section, we have provided a set of formulae to compute these constants from ${\mathit{\upsilon }}_{\mathit{i}}^{\mathrm{\prime }}$, which are the physical velocities of the particle with respect to an assumed emitter, which also has physical velocities υ_i with respect to an LNRF reference. According to these formulae, we provide a subroutine named lambdaqm to calculate λ, q, m, ε, and k_(a) defined by Eqs. (91)–(94). Except for a factor γ′μ_m, k_(a) actually are exactly equal to the initial four-momentum of the test particle given in an LNRF. Thus, k_(a) can be used to determine the signs in front of Π_r or Π_μ. The other initial parameters that need to be specified include: (1) the initial coordinates of the particle, r_ini, θ_ini, φ_ini, and t_ini. The latter two are usually set at zero; (2) the physical velocities of the assumed emitter with respect to an LNRF, υ_r, υ_θ, and υ_φ; (3) the physical velocities of the particle with respect to the assumed emitter, ${\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}$, ${\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}$, and ${\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}$; (4) the spin parameter a and the electric charge e of the black hole. With a given p and those initial parameters, we can do the calculations directly without giving the number of times that the particle meets the two turning points, namely Nt₁ and Nt₂.

In our code, the parameter p is an independent variable, which is always positive and monotonously increasing along a particular geodesic. When the geodesic is unbounded, it has a termination, either at infinity or the event horizon. The value of p corresponding to the termination is a finite number, denoted by p_max. We provide a subroutine named ptotal to calculate this number. Apparently, when a p given by the user is bigger than p_max, it has no meaning and the code resets it to be p_max mandatorily. When a geodesic is bounded, its termination does not exist at all and p can take any positive value.

For a more detailed introduction, refer to the README³ file. In the next section, we provide the results of our code for toy problems.

6. Applications for toy problems

To show the utility of our code, we apply it to toy problems. The results for six such examples are illustrated in this section.

6.1. Geodesic orbits of massive particles

The most important application of the code is to compute the geodesics of massive particles in a K-N spacetime. As in the first application, we use the code to compute the orbits of a set of test particles that are emitted isotropically in the local rest frame of an assumed emitter. The particles have a constant speed ${\mathit{\upsilon }}_{\mathrm{pt}}^{\mathrm{\prime }}$, but different directions in the local reference, and ${\mathit{\upsilon }}_{\mathrm{pt}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0.6}$. The orientation of the velocity is described by ϑ and ϕ, thus the components of the velocity under the reference of the emitter are $\begin{array}{ccc}{\mathit{\upsilon }}_{\mathit{r}}^{\mathrm{\prime }}& \mathrm{=}& {\mathit{\upsilon }}_{\mathrm{pt}}^{\mathrm{\prime }}\sin \mathit{\vartheta }\cos \mathit{\varphi ,}\\ {\mathit{\upsilon }}_{\mathit{\theta }}^{\mathrm{\prime }}& \mathrm{=}& {\mathit{\upsilon }}_{\mathrm{pt}}^{\mathrm{\prime }}\sin \mathit{\vartheta }\sin \mathit{\varphi ,}\\ {\mathit{\upsilon }}_{\mathit{\phi }}^{\mathrm{\prime }}& \mathrm{=}& {\mathit{\upsilon }}_{\mathrm{pt}}^{\mathrm{\prime }}\cos \mathit{\vartheta }\mathit{.}\end{array}$The physical velocities of the emitter with respect to the LNRF reference are υ_r = 0, υ_θ = 0, in which only the φ component is not zero and can take on different values. We demonstrate the results in Fig. 3. It is shown that as the speed of emitter increases, more particles become unbounded, and the beaming effect becomes more significant.

6.2. The orbits of spherical motion

The circular orbits in the Kerr spacetime have significant applications in the standard geometrically thin accretion disk systems. A particle in the accretion flow loses its angular momentum by viscosity and moves inward slowly. Its angular velocity is far greater than its radial velocity. Thus, to be a good approximation, the particle moves in a circular orbit. The inner radius of the disk is usually located at the innermost stable circular orbit (ISCO). Based on this assumption, one can measure the black hole spin by fitting the line profiles or the continuous spectra. Actually, the circular orbit, in which the particle is confined in the equatorial plane of the black hole, can be regarded as a special case of spherical orbit. The radial velocity and acceleration of the particle in a spherical orbit are vanished, leading to two conditions: dr/dτ = 0 and d²r/dτ² = 0. Using Eq. (4), these conditions reduce to (Bardeen et al. 1972; Wilkins 1972): $\begin{array}{ccc}{\mathit{R}}_{\mathit{r}}\mathrm{=}\mathrm{0}\mathit{,} \frac{\mathrm{d}{\mathit{R}}_{\mathit{r}}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{0.}& & \end{array}$(106)We use θ_∗ to denote the coordinate of one of the θ turning points, therefore we have Θ_θ(θ_∗) = 0. Using the same strategy discussed in Shakura (1987), we can obtain the angular velocity of the particle at θ_∗$\begin{array}{ccc}{\mathrm{\Omega }}^{\mathrm{\ast }}\mathrm{=}\frac{\sqrt{\mathit{P}}}{\sin {\mathit{\theta }}_{\mathrm{\ast }}\mathrm{(}\mathrm{\pm }\mathrm{\Sigma }\sqrt{\mathit{r}}\mathrm{+}\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\sqrt{\mathit{P}}\mathrm{)}}\mathit{,}& & \end{array}$(107)where P = M(r² − a²cos²θ_∗) − e²r, and the constants of motion: $\begin{array}{ccc}\frac{\mathit{E}}{{\mathit{\mu }}_{\mathit{m}}}& \mathrm{=}& \frac{\mathrm{\pm }\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\sqrt{\mathit{P}}\mathrm{+}\mathrm{(}\mathrm{\Delta }\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}\sqrt{\mathit{r}}}{\sqrt{\mathrm{\Sigma }}\sqrt{\mathrm{-}\mathit{P}\mathrm{+}\mathrm{(}\mathrm{\Delta }\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}\mathit{r}\mathrm{\pm }\mathrm{2}\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\sqrt{\mathit{rP}}}}\mathit{,}\\ \frac{\mathit{L}}{{\mathit{\mu }}_{\mathit{m}}}& \mathrm{=}& \frac{\sin {\mathit{\theta }}_{\mathrm{\ast }}\mathrm{[}\mathrm{\pm }\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}\sqrt{\mathit{P}}\mathrm{-}\sqrt{\mathit{r}}\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\mathrm{(}\mathrm{2}\mathit{Mr}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}}{\sqrt{\mathrm{\Sigma }}\sqrt{\mathrm{-}\mathit{P}\mathrm{+}\mathrm{(}\mathrm{\Delta }\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}\mathit{r}\mathrm{\pm }\mathrm{2}\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\sqrt{\mathit{rP}}}}\mathit{,}\end{array}$$\begin{array}{ccc}\frac{\mathit{Q}}{{\mathit{\mu }}_{\mathit{m}}^{\mathrm{2}}}& \mathrm{=}& {\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\left[{\mathit{a}}^{\mathrm{2}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{E}}^{\mathrm{2}}}{{\mathit{\mu }}_{\mathit{m}}^{\mathrm{2}}}\right)\mathrm{+}\frac{\mathrm{1}}{{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}}\frac{{\mathit{L}}^{\mathrm{2}}}{{\mathit{\mu }}_{\mathit{m}}^{\mathrm{2}}}\right]\mathit{,}\\ & \mathrm{=}& \frac{\mathit{r}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{[}{\mathit{A}}_{\mathit{Q}}\mathrm{\mp }\mathrm{2}\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\mathrm{(}\mathrm{2}\mathit{Mr}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathit{r}{\sqrt{\mathit{rP}}}^{\mathrm{]}}}{\mathrm{\Sigma }\mathrm{[}\mathrm{-}\mathit{P}\mathrm{+}\mathrm{(}\mathrm{\Delta }\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}\mathit{r}\mathrm{\pm }\mathrm{2}\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\sqrt{\mathit{rP}}\mathrm{]}}\mathit{,}\end{array}$(110)where $\begin{array}{ccc}{\mathit{A}}_{\mathit{Q}}& \mathrm{=}& \mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathit{Mr}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\{\mathrm{[}\mathit{Mr}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{}\\ & & \mathrm{+}\hspace{0.17em}{\mathit{e}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{]}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{-}\mathrm{(}\mathrm{2}\mathit{Mr}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\}\mathit{.}\end{array}$(111)In these formulae, the upper sign refers to the prograde orbits (i.e., corotating with L > 0), while the lower sign refers to retrograde orbits (counter rotating with L < 0).

In Fig. 4, we plot the orbit of a particle in a spherical motion. Given the parameters: a,e,r, and θ_∗, from Eqs. (108)–(110), we can obtain the constants of motion: λ,q,m, with which from equation R_r = 0 (or dR_r/dr = 0), we can obtain the final constant ε. For simplicity, we let both e and ε be zero in this figure. In comparison with circular motion, the most significant effect of spherical motion is the precession of the orbit.

Correspondingly, the spherical motion also has three kinds of marginal orbits, which are:

(1): The photon orbit r_ph, which is the innermost boundaryof the spherical orbits for particles, occurs when the de-nominator of Eqs. (108)–(110) vanishes, i.e.,$\begin{array}{ccc}\mathrm{-}\mathit{P}\mathrm{+}\mathrm{(}\mathrm{\Delta }\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}\mathit{r}\mathrm{\pm }\mathrm{2}\mathit{a}\sin {\mathit{\theta }}_{\mathrm{\ast }}\sqrt{\mathit{rP}}\mathrm{=}\mathrm{0.}& & \end{array}$(112)

(2): The marginally bound spherical orbit r_mb, which occurs when E/μ_m = 1.

(3): The innermost stable spherical orbit (ISSO) r_ms. The stable condition requires that d${}^{\mathrm{2}}\mathit{R}_{\mathit{r}}\mathit{/}\mathrm{d}{\mathit{r}}^{\mathrm{2}}\mathrm{\le }\mathrm{0}$, which yields the equivalent condition, $\begin{array}{ccc}\mathrm{1}\mathrm{-}{{\left(\frac{\mathit{E}}{{\mathit{\mu }}_{\mathit{m}}}\right)}^{\mathrm{2}}|}_{\mathit{e}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}\mathrm{\ge }\frac{\mathrm{2}\mathit{Mr}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{3}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}}{\mathrm{3}{\mathit{r}}^{\mathrm{4}}\mathrm{-}{\mathit{a}}^{\mathrm{4}}{\cos }^{\mathrm{4}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{-}\mathrm{6}{\mathit{r}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}}\mathit{,}& & \end{array}$(113)

or r ≥ r_ms. For simplicity, we have let e be zero in the above equation. If θ_∗ = π/2, namely the circular orbits, this condition reduces to the same form of Eq. (2.20) of Bardeen et al. (1972).

In our code, we provide three functions named r_ms,r_mb, and r_ph to compute the radii of these orbits with given a,θ_∗,e. In Fig. 5, we plot the radii of inner most stable spherical orbits as functions of a for various θ_∗. For simplicity, we also let e = 0 in this figure. One can see that as θ_∗ increases, the radii become smaller for a > 0 and larger for a < 0. For a = 0, the radii remain unchanged. When θ_∗ = 0°, the curve becomes symmetry for a > 0 and a < 0. Similar properties can be obtained for the radii of photon orbits and marginally bound orbits.

6.3. Orbits inside r_ms

The region inside ISSO is usually called the plug region in which a particle moves along geodesics with constants of motion of the marginally stable spherical geodesic (Cunningham 1975) when its initial radial perturbation velocity υ_r = 0. With the results presented in the above section, we obtain the constants of motion for the marginally stable spherical orbits $\begin{array}{ccc}{\left(\frac{\mathit{E}}{{\mathit{\mu }}_{\mathit{m}}}\right)}^{\mathrm{2}}& \mathrm{=}& \mathrm{1}\mathrm{-}\frac{\mathrm{2}\mathit{M}{\mathit{r}}_{\mathrm{ms}}\mathrm{(}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{2}}\mathrm{-}\mathrm{3}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}}{{\mathit{D}}_{\mathrm{ms}}}\mathit{,}\\ {\left(\frac{\mathit{L}}{{\mathit{\mu }}_{\mathit{m}}}\right)}^{\mathrm{2}}& \mathrm{=}& \mathrm{2}\mathit{M}{\mathit{r}}_{\mathrm{ms}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{[}\mathrm{3}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{4}}\mathrm{-}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}\\ & & \mathrm{+}\mathrm{3}{\mathit{a}}^{\mathrm{4}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{]}\mathit{/}{\mathit{D}}_{\mathrm{ms}}\mathit{,}\end{array}$$\begin{array}{ccc}\frac{\mathit{Q}}{{\mathit{\mu }}_{\mathit{m}}^{\mathrm{2}}}\hspace{0.17em}\hspace{0.17em}& \mathrm{=}& \frac{\mathrm{2}\mathit{M}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{3}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{(}\mathrm{3}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{)}}{{\mathit{D}}_{\mathrm{ms}}}\mathit{,}\\ {\mathit{D}}_{\mathrm{ms}}& \mathrm{=}& \mathrm{3}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{4}}\mathrm{-}{\mathit{a}}^{\mathrm{4}}{\cos }^{\mathrm{4}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{-}\mathrm{6}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathit{,}\end{array}$and

$\begin{array}{ccc}& & {\mathit{R}}_{\mathit{r}}\mathrm{=}\mathrm{2}\mathit{M}\mathrm{(}{\mathit{r}}_{\mathrm{ms}}\mathrm{-}\mathit{r}{\mathrm{)}}^{\mathrm{3}}\mathrm{[}\mathit{r}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{3}}\mathrm{-}\mathrm{3}{\mathit{a}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{ms}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{(}{\mathit{r}}_{\mathrm{ms}}\mathrm{+}\mathit{r}\mathrm{)}\\ & & \hspace{0.17em}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{a}}^{\mathrm{4}}{\cos }^{\mathrm{4}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{]}\mathit{/}{\mathit{D}}_{\mathrm{ms}}\mathit{,}\\ & & {\mathrm{\Theta }}_{\mathit{\mu }}\mathrm{=}\mathrm{2}\mathit{M}{\mathit{r}}_{\mathrm{ms}}\mathrm{(}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\mathrm{\left[}\mathrm{\right(}\mathrm{3}{\mathit{a}}^{\mathrm{4}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{-}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{)}{\mathit{\mu }}^{\mathrm{2}}\\ & & \hspace{0.17em}\hspace{0.17em}\mathrm{-}\hspace{0.17em}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{\ast }}\mathrm{+}\mathrm{3}{\mathit{r}}_{\mathrm{ms}}^{\mathrm{4}}\mathrm{]}\mathit{/}{\mathit{D}}_{\mathrm{ms}}\mathit{.}\end{array}$From the above two equations, we know that both r_ms and θ_∗ are the turning points. With these expressions, we can find the constants of motion immediately to compute the geodesic orbits inside r_ms. In Fig. 6, we plot such an orbit. We take a = 0.998 and θ_∗ = 0, namely the particle goes through the spin axis of the black hole. Using the function r_ms(θ_∗,a,e) in our code, we find r_ms = 5.2781 r_g. From this figure, one can see that the orbit is almost the same with a spherical motion because the radial velocity is much smaller than the poloidal and azimuthal velocities.

6.4. Accretion flow of disk

Now we use our code to construct a toy model for mimicking the accretion flows of a skewed geometrically thin disk. The flows are composed of non-interacting particles, which fall freely into the black hole along the geodesic trajectories. This implies that we make a ballistic treatment to the fluid flow, and the dynamics and the structure of the disk are uniquely determined by the gravitational field of the black hole. It is also convenient to regard the accretion flow as a collection of test particles with the same mass. The boundary conditions of the disk are assumed to be a ring at r = r₀, from which the test particles are continuously injected. The plane of the ring has an inclination angle β with respect to the spin axis.

To describe the initial conditions of the test particles, namely the velocities, an orthonormal tetrad is established on the ring. We choose three spacial basis vectors of the tetrad as : e_x,e_y,e_z, where e_x is the tangent vector of the ring, e_y is aligned along the radial direction and pointed inward, e_z is the normal vector of the plane of the ring, and these vectors satisfy right hand rule. In the local rest frame of the tetrad, the physical velocities of the particle are ${\mathit{\upsilon }}_{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\upsilon }}_{\mathit{y}}^{\mathrm{\prime }}\mathit{,}{\mathit{\upsilon }}_{\mathit{z}}^{\mathrm{\prime }}$. According to the discussions in Sect. 4, to compute the constants of motion from these velocities, we need to transform them into the LNRF reference for getting υ_r,υ_θ,υ_φ.

We denote the transformation matrix by T(θ,φ,β). The explcit expression of T is presented in Appendix A. Hence, we have ${\mathit{\upsilon }}_{\mathit{i}}\mathrm{=}{\mathit{T}}_{\mathit{i}}^{\mathit{j}}{\mathit{\upsilon }}_{\mathit{j}}^{\mathrm{\prime }}$. Since the tetrad is attached on the ring, θ and φ are not independent variables, actually they satisfy the following equation, $\begin{array}{ccc}\sin \mathit{\theta }\mathrm{=}\frac{\cos \mathit{\beta }}{\sqrt{\mathrm{1}\mathrm{-}\sin {\mathit{\phi }}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\beta }}}\mathrm{·}& & \end{array}$(120)Therefore, we have ${\mathit{\upsilon }}_{\mathit{i}}\mathrm{=}{\mathit{T}}_{\mathit{i}}^{\mathit{j}}\mathrm{\left(}\mathit{\phi ,\beta }\mathrm{\right)}{\mathit{\upsilon }}_{\mathit{j}}^{\mathrm{\prime }}$. Taking φ as an independent variable that varies from 0 to 2π, we can set the initial conditions for all test particles. In Fig. 7, we plot such a set of geodesic orbits of test particles with the same mass. The physical velocities are ${\mathit{\upsilon }}_{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0.01}$, ${\mathit{\upsilon }}_{\mathit{y}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0.5}$, and ${\mathit{\upsilon }}_{\mathit{z}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathrm{0.01}$, respectively. From the figure, one can see that all trajectories form a smooth but curved surface.

Using the ray-tracing approach (Luminet 1979), we can image the curved surface. In Yang & Wang (2013), we presented a new public code named ynogk to compute the null geodesics in a Kerr spacetime and a more general method to image a target object. The method requires one to provide the function describing the surface, i.e., F(r,θ,φ) = 0, or F(x,y,z) = 0. For this curved surface formed by geodesic orbits of test particles, we cannot write out its explicit form, and only use the interpolation approach. To approximate the surface, we take N particles with N geodesic orbits. We take M points in each orbit, and get a total of N × M points. We can get the coordinates of each point easily and write them as: z_ij = z_ij(x_i,y_j). Using the interpolation approaches provided in Press et al. (2007), we can obtain an approximation function z = z(x,y) that describes the surface.

In Fig. 8, we plot the images of a skewed accretion disk that is composed of test particles falling freely into a black hole along geodesics viewed from different inclinations. The initial tiled angle of the disk is β = 30°. Due to the frame drag effect, the particles drift into the black hole along spiral orbits, instead of a straight lines. The orbits make a gradual transition into the equatorial plane. The disk is significantly warped as it moves inward. In the figure, the false color represents the redshift of emission coming from the disk surface. The approaching and receding sides of the disk are no longer the left and right sides, but the regions are farthest and nearest with respect to the observer, respectively.

6.5. Stationary axisymmetric accretion flow

Tejeda et al. (2013) presented an analytic toy model to mimic the stationary axisymmetric accretion flow of a rotating cloud of non interacting particles falling onto a Kerr black hole in which the streamlines are described analytically in terms of time-like geodesics. Thus, they solve the equations of motion with integral forms with elliptic functions. However, their results are completely different from ours. In addition, they just obtain the solutions for r and μ.

As a check of the validation of our code, we model a similar accretion flow. Since the flow is axisymmetric, we just need to consider the spacial projections onto the r − θ plane. The boundary of the flow is a spherical shell at r = r₀ from which test particles are continuously injected. On the shell, the four-velocity of particles are taken to be constants, i.e., $\begin{array}{ccc}{\mathit{u}}^{\mathit{r}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,\theta }\mathrm{\right)}\mathrm{=}{\mathit{u}}_{\mathrm{0}}^{\mathit{r}}\mathit{,} {\mathit{u}}^{\mathit{\theta }}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,\theta }\mathrm{\right)}\mathrm{=}{\mathit{u}}_{\mathrm{0}}^{\mathit{\theta }}\mathit{,} {\mathit{u}}^{\mathit{\phi }}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,\theta }\mathrm{\right)}\mathrm{=}{\mathit{u}}_{\mathrm{0}}^{\mathit{\phi }}\mathit{,}& & \end{array}$(121)with which and the identity g_μνu^μu^ν =  −1, we obtain the time component ${\mathit{u}}_{\mathrm{0}}^{\mathit{t}}$. The four-velocity ${\mathit{u}}_{\mathrm{0}}^{\mathit{\mu }}$ at r₀ are taken as the initial conditions for the calculation of geodesics. Using Eq. (77) the physical velocities υ_r,υ_θ,υ_φ can be computed from ${\mathit{u}}_{\mathrm{0}}^{\mathit{\mu }}$. Then the constants of motion are uniquely determined. The results are illustrated in Fig. 9, which agree well with those of Tejeda et al. (2013).

6.6. The tidal disruption of a ball

As the final application of our code, we model a tidal disruption event of a ball, which falls freely to the central black hole. To use the code, we have to assume that the ball consists of a set of equal mass test particles without any interactions. At the initial point the ball has a kick velocity, the physical components of which measured under the LNRF reference are υ_r,υ_θ,υ_φ. All of the particles share the same initial velocities, but different positions. Then with the given initial conditions, each particle freely falls inward along a geodesic trajectory.

In Fig. 10, we show the deformed images of the ball for five different coordinate times. Initially, we assume that the shape of the ball is a regular sphere and the center of the ball is located in the equatorial plane of the black hole. The velocities are υ_r =  −0.1,υ_θ = 0,υ_φ = 0.1. The shape of the ball is significantly deformed and stretched as it approaches the central black hole due to the strong tidal disruption force. The former part is stretched more seriously than the latter part. The debris of the ball orbits around the black hole along a spiral trajectory and goes inward slowly instead of falling into the black hole directly because of the frame drag effect. However, the picture illustrated by this example is definitely a toy model.

7. Discussion and conclusion

We have developed a new fast public code named ynogkm for calculating time-like geodesics under a K-N spacetime, which is a direct extension of Yang & Wang (2013). In ynogkm, we adopt the same strategies used in ynogk, i.e., expressing all coordinates and proper times as functions of a parameter p and calculating all elliptic integrals with Carlson’s approach. The former guarantees the convenience of the code in practical application and the latter guarantees the fast speed of the code. The extension is involved in many more complicated cases.

In the expressions, we also use the Weierstrass elliptic function ℘(z;g₂,g₃) and integral as investigated by many authors in the literature. They not only investigate the geodesic motion itself but also the properties of the spacetime, while what we discussed in this paper focuses on the potential real applications of the calculation of geodesic orbits in astrophysics. To avoid complex integrals, we also adopt the Jacobi’s elliptic functions sn(z | k²),cn(z | k²) when equation R(r) = 0 has no real roots.

Since ynogkm uses the same strategies as ynogk and their speeds are almost the same, we do not present the speed test results. As discussed in Chan et al. (2013), a powerful approach to improve the speed of tracing the trajectories of billions of photons in a curved spacetime is based on the massively parallel algorithm and GPU graphic cards. Their results show that this approach is two orders of magnitude faster than the CPU-based tracing codes. Therefore, a future work will involve the extension of ynogkm from a serial program to a parallel program.

To demonstrate the utility of ynogkm, we apply it to six toy problems and present the results simply. Its application to more complicated and practical cases will be provided in a future work.

---

Acknowledgments

We acknowledge the anonymous referee for his/her valuable comments and advice, which significantly improved the manuscript. We acknowledge financial support from the National Natural Science Foundation of China 11133006, 11163006, 11173054, the National Basic Research Program of China (973 Program 2009CB824800), and the Policy Research Program of Chinese Academy of Sciences (KJCX2-YW-T24).

References

Appendix A: Transformation matrix

Here we discuss how to get the explicit expression for the matrix T, which transforms the physical velocities ${\mathit{\upsilon }}_{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\upsilon }}_{\mathit{y}}^{\mathrm{\prime }}\mathit{,}{\mathit{\upsilon }}_{\mathit{z}}^{\mathrm{\prime }}$ of a particle specified in the reference of the tetrad ${{e}}_{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{{e}}_{\mathit{y}}^{\mathrm{\prime }}\mathit{,}{{e}}_{\mathit{z}}^{\mathrm{\prime }}$ into the LNRF reference whose origin is fixed at the same point. As shown in Fig. A.1, we have four references, i.e., R′: {p, ${{e}}_{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{{e}}_{\mathit{y}}^{\mathrm{\prime }}\mathit{,}{{e}}_{\mathit{z}}^{\mathrm{\prime }}$}, R′′: {O,x′′,y′′,z′′}, R: {O,x,y,z}, and R_rθφ: {p, e_r,e_θ,e_φ}.

The matrix of transformation T₁ from R′ into R′′ can be obtained directly, $\begin{array}{ccc}\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathit{.}& & \end{array}$(A.1)The transformation T₂ from R′′ into R is given by, $\begin{array}{ccc}\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathit{.}& & \end{array}$(A.2)The transformation T₃ from R into R_rθφ is given by, $\begin{array}{ccc}\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathit{.}& & \end{array}$(A.3)Thus the transformation T from R′ into R_rθφ is given by $\begin{array}{ccc}\mathit{T}\mathrm{=}{\mathit{T}}_{\mathrm{3}}{\mathit{T}}_{\mathrm{2}}{\mathit{T}}_{\mathrm{1}}\mathit{,}& & \end{array}$(A.4)and noting that ψ = γ + φ, one has $\begin{array}{ccc}{\mathit{T}}_{\mathrm{3}}{\mathit{T}}_{\mathrm{2}}\mathrm{=}(\begin{array}{c}\\ & \\ & \\ & \end{array}& & \\ \begin{array}{c}\\ \\ \\ \end{array})\mathit{.}& & \end{array}$Finally, we obtain $\begin{array}{ccc}\mathit{T}\mathrm{=}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{.}& & \end{array}$(A.5)In the reduction, the following identities are used

$\begin{array}{ccc}& & \sin \mathit{\theta }\mathrm{=}\frac{\cos \mathit{\beta }}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{\phi }{\sin }^{\mathrm{2}}\mathit{\beta }}}\mathit{,} \cos \mathit{\theta }\mathrm{=}\mathrm{-}\frac{\cos \mathit{\phi }\sin \mathit{\beta }}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{\phi }{\sin }^{\mathrm{2}}\mathit{\beta }}}\mathit{,}\\ & & \sin \mathit{\varphi }\mathrm{=}\frac{\sin \mathit{\phi }\cos \mathit{\beta }}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{\phi }{\sin }^{\mathrm{2}}\mathit{\beta }}}\mathit{,} \cos \mathit{\varphi }\mathrm{=}\frac{\cos \mathit{\phi }}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{\phi }{\sin }^{\mathrm{2}}\mathit{\beta }}}\mathrm{·}\end{array}$

Appendix B: Taking t and σ to be the independent variable

In some practical applications, we prefer using t or σ as the independent variable to parameter p. Since we have expressed all B-L coordinates and proper times as functions of the parameter p, when a value of t or σ is given, there is an unique p corresponds to it. Namely, both the equations t(p) = t₀ and σ(p) = σ₀ have one and only one real root, which is denoted by p₀. Apparently, if we can solve these equations efficiently and precisely to get p₀, we can take t or σ as the independent variable, for p₀ is obtained, the other three coordinates r,μ, and φ are also uniquely determined.

Actually, we can solve both equations t(p) = t₀ and σ(p) = σ₀ by the bisection method or iterative method. From the expressions for t_r,t_μ and σ_r,σ_μ given in Sect. 3.3.2 we can rewrite functions t(p) and σ(p) as $\begin{array}{ccc}\mathit{\sigma }\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\overline{\mathit{\sigma }}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathit{\sigma }}\mathit{p,} \mathit{t}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\overline{\mathit{t}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathit{t}}\mathit{p,}& & \end{array}$(B.1)where the definitions of C_σ and C_t are provided in Table B.1. Thus for a given σ₀ and t₀, we have $\begin{array}{ccc}\mathit{p}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{0}}\mathrm{-}\overline{\mathit{\sigma }}\mathrm{\left(}\mathit{p}\mathrm{\right)}}{{\mathit{C}}_{\mathit{p}}}\mathrm{=}{\mathit{f}}_{\mathit{\sigma }}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathit{,} \mathit{p}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{0}}\mathrm{-}\overline{\mathit{t}}\mathrm{\left(}\mathit{p}\mathrm{\right)}}{{\mathit{C}}_{\mathit{t}}}\mathrm{=}{\mathit{f}}_{\mathit{t}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathit{.}& & \end{array}$(B.2)

We illustrate schematically how to solve these equations using an iterative method in Fig. B.1. To use the bisection method, we define two new functions $\begin{array}{ccc}{\mathit{F}}_{\mathit{\sigma }}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{\sigma }}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{-}\mathit{p,} {\mathit{F}}_{\mathit{t}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{t}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{-}\mathit{p}\mathit{.}& & \end{array}$(B.3)Figure B.1 shows that when p < p₀, F_σ(p) or F_t(p) > 0; when p > p₀, F_σ(p) or F_t(p) < 0. Thus through the use of bisection method, we can solve the equations F_σ(p) = 0 and F_t(p) = 0 immediately.

All Tables

All Figures