Construction of Explicit Symplectic Integrators in General Relativity. IV. Kerr Black Holes

Gravitational waves and black holes are two fundamental predictions of the theory of general relativity of Einstein. The predictions have been frequently confirmed by a number of observations of gravitational waves from binary black hole or neutron star mergers (e.g., GW 150914 (Abbott et al. 2016) and GW 190521 (Abbott et al. 2020)) and images of the supermassive black hole candidate in the center of the giant elliptical galaxy M87 (EHT Collaboration et al. 2019). The observed image provides powerful evidence for the presence of a rotating black hole, which was derived by Kerr (1963) from the field equations of general relativity. There are also other solutions for the field equations, such as a nonrotating Schwarzschild black hole.

The geodesics of the Schwarzschild, Reissner–Nordström, Reissner–Nordström–(anti)–de Sitter, and Kerr spacetime geometries are highly nonlinear but are integrable. This integrability does not mean that their solutions can be expressed in terms of elementary functions but quadratures. Thus, numerical integration schemes are still necessarily used to solve these geodesic equations. The popular fourth-order Runge–Kutta explicit integration method (RK4) is applicable for arbitrary metrics (i.e., arbitrary spacetimes in arbitrary coordinates). Recently, the RAPTOR code with RK4 was applied to study accretion models of supermassive Kerr black holes so as to produce physically accurate images of black hole accretion disks (Bronzwaer et al. 2018, 2020). The RK4 integrator is an accurate and efficient integrator for a short integration time. However, it would yield a secular drift in energy errors and would provide unreliable numerical results for a long-term integration. The manifold correction of Nacozy (1971) and its extensions (Fukushima 2003; Wu et al. 2007, Ma et al. 2008; Wang et al. 2016, 2018; Deng et al. 2020) are helpful to compensate for the defect of RK4 by pulling the integrated orbit back to the original integral hypersurface.

The RK4 integrator combined with a manifold correction scheme is regarded as one of the geometric integration algorithms that preserve structures, integrals, symmetries, reversing symmetries, and phase-space volumes (Hairer et al. 1999). The manifold correction scheme can strictly satisfy the energy integral but does not conserve the symplecticity. A class of discrete Hamiltonian gradient schemes that conserve energy to machine precision can be used for any globally hyperbolic spacetimes with six dimensions (Bacchini et al. 2018, 2019) and eight- and ten-dimensional Hamiltonian problems (Hu et al. 2019, 2021). These energy-conserving discrete gradient integrators are implicit, nonsymplectic, and do not preserve other integrals in general. ⁴ Their constructions become complex as the dimensionality increases. No higher-order accuracy but only first- and second-order accuracies can be given to numerical solutions in the existing energy-conserving discrete gradient schemes.

Symplectic integrators (Swope et al. 1982; Wisdom 1982; Ruth 1983; Forest & Ruth 1990; Wisdom & Holman 1991; Chambers & Murison 2000; Laskar & Robutel 2001; Omelyan et al. 2003) are the best geometric integrators for studying the long-term evolution of the geodesics and other Hamiltonian problems. Although they are different from the manifold correction schemes and energy-conserving integrators that can exactly preserve the energy integral, they remain bounded and show no secular growth in energy errors. Other constants of motion and symplectic geometrical structure are also conserved in the course of numerical integrations. Because of the difficulty of the separation of variables in curved spacetimes, the standard explicit symplectic integrators become useless. Instead, implicit symplectic methods (Feng 1986; Brown 2006; Tsang et al. 2015) or implicit and explicit mixed symplectic methods (Preto & Saha 2009; Kopáček et al. 2010; Lubich et al. 2010; Zhong et al. 2010; Mei et al. 2013a, 2013b) are used.

It is well known that an explicit algorithm is generally superior to the same-order implicit method in computational efficiency. Noting this fact, Pihajoki (2015) presented explicit leapfrog integration schemes for inseparable Hamiltonian systems including arbitrary metrics by doubling the phase space and introducing a new Hamiltonian with two variable-separating parts equal to the original Hamiltonian. Of course, the leapfrog algorithms are symplectic in the extended phase space. However, the two parts are coupled and dependent through the derivatives, and therefore, the two numerical flows for the two parts diverge with time. Mixing maps that act as feedback between the two solutions would address the problem of the divergence of both solutions. If the mixing maps are not symplectic, then the extended phase-space leapfrogs are not, either. Even if the mixing maps are symplectic, the extended phase-space leapfrogs are not symplectic when projection maps are used to project a vector in extended phase space back to the original phase space in any case. In spite of this, the leapfrogs still preserve the original Hamiltonian without secular growth in the error because the mixing maps and projection maps are operated only in obtaining outputs and do not alter the state in the extended phase space. In this sense, these algorithms are viewed as extended phase-space explicit symplectic-like integrators. Permutations of momenta were shown to be the best mixing maps. Liu et al. (2016) proposed fourth-order extended phase-space explicit symplectic-like methods and found that the sequent permutations of coordinates and momenta are superior to the permutations of momenta. Luo et al. (2017) showed that the best choice for permuted maps should be midpoint permutations. These extended phase-space explicit symplectic-like integrators are applicable for various inseparable problems (e.g., a vast family of spacetimes and post-Newtonian problems; Li & Wu 2017; Wu & Wu 2018; Li & Wu 2019). On the other hand, Tao (2016) replaced the mixing maps with a third part as an artificial restraint on the binding of the two copies of the original system with mixed-up positions and momenta. In this way, the extended phase-space leapfrogs are symplectic. More recently, FANTASY based on this idea was applied to allow for the integration of geodesics in arbitrary spacetimes with automatic differentiation (Christian & Chan 2021). However, there is an open problem on how to determine the most appropriate constant for controlling the binding of the two copies. This choice is not given in a universal method but relies on many numerical tests. Wu & Wu (2018) reported that Tao's method with an appropriate choice of control constant is not better than the method with midpoint permutations in accuracy. Another problem is that the extended phase-space leapfrogs are not symplectic when restricted to the original phase space in any case, as Pihajoki claimed.

Are the standard explicit symplectic methods not applicable for general relativistic metrics? No is the key to this question. Recently, Wang et al. (2021a, Paper I) successfully separated the Hamiltonian of Schwarzschild black hole into four integrable parts with analytical solutions as explicit functions of proper time and used these explicit analytical solutions to compose the standard second- and fourth-order explicit symplectic integrators. When the Hamiltonian of a Reissner–Nordström black hole has five similar splitting parts, the standard explicit symplectic integrators were easily set up in Paper II (Wang et al. 2021b). The standard explicit symplectic integrators were also designed for the Hamiltonian of a Reissner–Nordström–(anti)–de Sitter black hole with six integrable separable parts in Paper III (Wang et al. 2021c).

Unfortunately, the standard explicit symplectic integrators become useless if the Hamiltonian of a Kerr black hole is split according to the splitting techniques of the Hamiltonians of nonrotating black holes in Papers I, II, and III. This is because Σ as a function of r and θ exists in the denominators of the Hamiltonian and leads to inseparable variables or splitting Hamiltonian parts without the desired analytical solutions. To overcome this difficulty, we use the time transformation method introduced by Mikkola (1997) to obtain a time-transformed Hamiltonian in which the denominators do not contain the function Σ. In this way, the standard explicit symplectic algorithms can be available for the time-transformed Hamiltonian. This is the main aim of this paper.

The rest of this paper is organized as follows. In Section 2 we introduce the symplectic integrators with adaptive time steps of Mikkola (1997). Then, the Kerr geometry is described in Section 3. Explicit symplectic algorithms are designed for the Kerr geometry in Section 4. We check the performance of the proposed algorithms in Section 5. Finally, the main results are concluded in Section 6. The Carter constant and the parameters and initial conditions of unstable spherical photon orbits in the Kerr spacetime are described in Appendix A. Codes of the new second-order method are given in Appendix B. Other choices of the time transformation function are presented in Appendix C.

Consider a perturbed two-body problem with the Hamiltonian

$$\begin{eqnarray}&&H({\boldsymbol{p}},{\boldsymbol{r}},\tau )={H}_{0}+R({\boldsymbol{r}},\tau ),{H}_{0}=\displaystyle \frac{{{\boldsymbol{p}}}^{2}}{2}-\displaystyle \frac{\mu }{r},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{p}}$ denotes a momentum vector, τ is a physical time, and μ is a constant associated with the constant of gravity and masses of the two bodies. The Kepler part H₀ is integrable, and so is the perturbing function R.

Take $\tau ={q}_{0}$ as a new coordinate, which corresponds to a conjugate momentum ${p}_{0}=-H$. The phase space is extended by ${\boldsymbol{Q}}=({q}_{0},{\boldsymbol{r}})$ and ${\boldsymbol{P}}=({p}_{0},{\boldsymbol{p}})$. By introducing a fictitious time variable w through the relation

$$\begin{eqnarray}&&d\tau =g({\boldsymbol{r}},{q}_{0}){dw},\end{eqnarray} \tag{ 2 }$$

Mikkola (1997) obtained an extended phase-space Hamiltonian

$$\begin{eqnarray}&&{\rm{\Gamma }}={{\rm{\Gamma }}}_{0}+{{\rm{\Gamma }}}_{1},\,\,{{\rm{\Gamma }}}_{0}=g({H}_{0}+{p}_{0}),\,\,{{\rm{\Gamma }}}_{1}={gR}.\end{eqnarray} \tag{ 3 }$$

g and Γ are called the time transformation function and Hamiltonian, respectively. g is a constant along the orbit over the fictitious time step from the beginning of each step to the end of each step, but it changes at the beginning and at the end of each step. Even if Γ is explicitly dependent on the physical time coordinate q₀, it is identical to zero (i.e., ${\rm{\Gamma }}\equiv 0$) for any time w. ${{\rm{\Gamma }}}_{1}$ does not depend on any momenta and is thus easily solvable. However, ${{\rm{\Gamma }}}_{0}$ is difficult to analytically solve due to the inseparable variables. To overcome this problem, Mikkola used the time transformation $1/g$ (e.g., g = r) to transform ${{\rm{\Gamma }}}_{0}$ back to the physical time:

$$\begin{eqnarray}&&{\bar{H}}_{0}=\displaystyle \frac{1}{g}({{\rm{\Gamma }}}_{0}-\varepsilon )=\displaystyle \frac{{{\boldsymbol{p}}}^{2}}{2}-\displaystyle \frac{\mu +\varepsilon }{r}+{p}_{0},\end{eqnarray} \tag{ 4 }$$

where $\varepsilon ={{\rm{\Gamma }}}_{0}$ is the value of ${{\rm{\Gamma }}}_{0}$ at the beginning of the next step. ε is a constant along the orbit over the physical time step during the beginning and end of each step. Equation (4) is still a Kepler problem with a modified mass $\kappa =\mu +\varepsilon$. Of course, the mass at the beginning of one step is unlike that at the end of this step. Because ${dw}=d\tau /g=d\tau /r$, $h=w={\int }^{\tau }(1/r)d\tau$ is a time step. Based on Stumpff's form of Kepler's equation, the physical time step τ can be expressed in terms of h. The positions and velocities of ${\bar{H}}_{0}$ are also functions of h. Take ${\mathscr{A}}$ as a differential operator with respect to ${\bar{H}}_{0}$. The analytical solutions (i.e., the momentum jumps) of ${{\rm{\Gamma }}}_{1}$ are easily obtained by ${\rm{\Delta }}{\boldsymbol{p}}=-h\partial ({gR})/\partial {\boldsymbol{r}}$ and ${\rm{\Delta }}{p}_{0}=-h\partial ({gR})/\partial {q}_{0}$. ${\mathscr{B}}$ is viewed as another differential operator with respect to ${{\rm{\Gamma }}}_{1}$. In this way, Mikkola established the second-order symplectic leapfrog of Wisdom & Holman (1991):

$$\begin{eqnarray}&&{{\mathscr{S}}}_{2}^{{\rm{\Gamma }}}={\mathscr{B}}\left(\displaystyle \frac{h}{2}\right)\circ {\mathscr{A}}(h)\circ {\mathscr{B}}\left(\displaystyle \frac{h}{2}\right).\end{eqnarray} \tag{ 5 }$$

The operators ${\mathscr{A}}$ and ${\mathscr{B}}$ can also compose the fourth-order symplectic method of Forest & Ruth (1990)

$$\begin{eqnarray}\begin{array}{rcl}{{\mathscr{S}}}_{4}^{{\rm{\Gamma }}} & = & {\mathscr{B}}\left(\displaystyle \frac{h}{2{b}_{1}}\right)\circ {\mathscr{A}}\left(\displaystyle \frac{h}{{b}_{1}}\right)\circ {\mathscr{B}}\left(\displaystyle \frac{1-{b}_{2}}{2{b}_{1}}h\right)\circ {\mathscr{A}}\left(-\displaystyle \frac{{b}_{2}}{{b}_{1}}h\right)\\ & & \times \ \circ {\mathscr{B}}\left(\displaystyle \frac{1-{b}_{2}}{2{b}_{1}}h\right)\circ {\mathscr{A}}\left(\displaystyle \frac{h}{{b}_{1}}\right)\circ {\mathscr{B}}\left(\displaystyle \frac{h}{2{b}_{1}}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

where ${b}_{2}={2}^{1/3}$ and ${b}_{1}=2-{b}_{2}$. These adaptive integrators, Equations (5) and (6), nearly preserve the Hamiltonian Γ, i.e., the Hamiltonian (1). They demonstrate good qualitative properties of symplectic integrators with constant time steps. In addition, the use of variable steps causes the accuracy and efficiency to be significantly improved.

The above symplectic algorithms are implemented in the new time w. They rely on the analytical solutions of the Keplerian motion ${\bar{H}}_{0}$ and require that the physical time τ be given in an expressional form of h. It is a hard task. However, when the Hamiltonian (1) in the extended phase space is split into kinetic energy $T={{\boldsymbol{p}}}^{2}/2+{p}_{0}$ and potential energy $U=R({\boldsymbol{r}},{q}_{0})$, such similar explicit symplectic integrators are easily available for the logarithmic Hamiltonian method proposed by Mikkola & Tanikawa (1999) and the time transformation function suggested by Preto & Tremaine (1999). Extensions and applications of such adaptive time-step symplectic integrators were considered by Mikkola & Aarseth (2002), Emel'yanenko (2007), Preto & Saha (2009), Mikkola & Tanikawa (2013), Ni & Wu (2014), Li & Wu (2017), and Wang & Nitadori (2020).

Because the fixed time step h is used in the new fictitious time w, the good long-time behavior properties of symplectic methods are not lost. In addition, the physical time steps vary in different positions of an orbit. In particular, they become smaller as the perturbing body is closer to the central body. This leads to increasing precision. The two points are what the symplectic methods with time transformations satisfy.

The Kerr black hole is a rotating black hole. Its gravitational field is described by the spacetime metric

$$\begin{eqnarray}&&-d{\tau }^{2}={{ds}}^{2}={g}_{\alpha \beta }{{dx}}^{\alpha }{{dx}}^{\beta }.\end{eqnarray} \tag{ 7 }$$

In the standard Boyer-Lindquist coordinates $(t,r,\theta ,\phi )$, this metric is an axially symmetric stationary and has covariant nonzero components (Kerr 1963; Takahashi & Koyama 2009):

$$\begin{eqnarray*}\begin{array}{rcl}{g}_{{tt}} & = & -\left(1-\displaystyle \frac{2r}{{\rm{\Sigma }}}\right),\,\,{g}_{t\phi }=-\displaystyle \frac{2{ar}{\sin }^{2}\theta }{{\rm{\Sigma }}}={g}_{\phi t},\\ {g}_{{rr}} & = & \displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}},{g}_{\theta \theta }={\rm{\Sigma }},\\ {g}_{\phi \phi } & = & \left({\rho }^{2}+\displaystyle \frac{2{{ra}}^{2}}{{\rm{\Sigma }}}{\sin }^{2}\theta \right){\sin }^{2}\theta .\end{array}\end{eqnarray*}$$

Its contravariant nonzero components are

$$\begin{eqnarray*}\begin{array}{rcl}{g}^{{tt}} & = & -\displaystyle \frac{A}{{\rm{\Delta }}{\rm{\Sigma }}},\,\,{g}^{t\phi }=-\displaystyle \frac{2{ar}}{{\rm{\Delta }}{\rm{\Sigma }}}={g}^{\phi t},\\ {g}^{{rr}} & = & \displaystyle \frac{{\rm{\Delta }}}{{\rm{\Sigma }}},\,\,{g}^{\theta \theta }=\displaystyle \frac{1}{{\rm{\Sigma }}},\\ {g}^{\phi \phi } & = & \displaystyle \frac{{\rm{\Sigma }}-2r}{{\rm{\Delta }}{\rm{\Sigma }}{\sin }^{2}\theta }.\end{array}\end{eqnarray*}$$

Note that ${\rm{\Sigma }}={r}^{2}+{a}^{2}{\cos }^{2}\theta$, ${\rm{\Delta }}={\rho }^{2}-2r$, ${\rho }^{2}={r}^{2}+{a}^{2}$, and $A={\rho }^{4}-{\rm{\Delta }}{a}^{2}{\sin }^{2}\theta$. τ is a proper time, and t is a coordinate time. The speed of light c and the gravitational constant G use geometrized units, $c=G=1$. The mass of the black hole M also takes one unit, M = 1. In such a unit system, a denotes the angular momentum of the rotating body. The body should be slowly rotating, namely, $| a| \leqslant 1$.

Based on the spacetime metric, a Lagrangian system is given by

$$\begin{eqnarray}&&{\ell }=\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{ds}}{d\tau }\right)}^{2}=\displaystyle \frac{1}{2}{g}_{\mu \nu }{\dot{x}}^{\mu }{\dot{x}}^{\nu }.\end{eqnarray} \tag{ 8 }$$

Four-velocity ${\dot{x}}^{\mu }=(\dot{t},\dot{r},\dot{\theta },\dot{\phi })={U}^{\mu }={\boldsymbol{U}}$ or ${\boldsymbol{U}}={U}_{\mu }={g}_{\mu \nu }{U}^{\nu }={g}_{\mu \nu }{\dot{x}}^{\nu }$ satisfies the identical relation

$$\begin{eqnarray}&&{\boldsymbol{U}}\cdot {\boldsymbol{U}}={U}_{\mu }{U}^{\mu }={g}_{\mu \nu }{\dot{x}}^{\mu }{\dot{x}}^{\nu }=-1.\end{eqnarray} \tag{ 9 }$$

According to classical mechanics, a covariant generalized four-momentum is defined as

$$\begin{eqnarray}&&{p}_{\mu }=\displaystyle \frac{\partial {\ell }}{\partial {\dot{x}}^{\mu }}={g}_{\mu \nu }{\dot{x}}^{\nu }.\end{eqnarray} \tag{ 10 }$$

Clearly, each component of the metric tensor ${g}_{\mu \nu }$ does not explicitly depend on the coordinates t and ϕ. Based on the Euler–Lagrangian equations, the four-momentum exists with two constant components:

$$\begin{eqnarray}\begin{array}{rcl}{p}_{t}={g}_{{tt}}\dot{t}+{g}_{t\phi }\dot{\phi } & = & -\left(1-\displaystyle \frac{2r}{{\rm{\Sigma }}}\right)\dot{t}-\displaystyle \frac{2{ra}{\sin }^{2}\theta }{{\rm{\Sigma }}}\dot{\phi }\\ & = & -E,\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\begin{array}{rcl}{p}_{\phi }={g}_{\phi \phi }\dot{\phi }+{g}_{t\phi }\dot{t} & = & \left({\rho }^{2}+\displaystyle \frac{2{{ra}}^{2}}{{\rm{\Sigma }}}{\sin }^{2}\theta \right){\sin }^{2}\theta \dot{\phi }\\ & & -\ \displaystyle \frac{2{ra}{\sin }^{2}\theta }{{\rm{\Sigma }}}\dot{t}=L.\end{array}\end{eqnarray} \tag{ 12 }$$

E represents the energy of a test particle moving around the rotating body, and L is the angular momentum of the particle. Equations (11) and (12) can be rewritten as

$$\begin{eqnarray}&&\dot{t}=\displaystyle \frac{{dt}}{d\tau }=-{f}_{1}E-{f}_{2}L,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\dot{\phi }=\displaystyle \frac{d\phi }{d\tau }={f}_{2}E+{f}_{3}L,\end{eqnarray} \tag{ 14 }$$

where f₁, f₂, and f₃ are functions of r and θ as follows:

$$\begin{eqnarray}&&{f}_{1}=\displaystyle \frac{{g}_{\phi \phi }}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{f}_{2}=\displaystyle \frac{{g}_{t\phi }}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{f}_{3}=\displaystyle \frac{{g}_{{tt}}}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}}.\end{eqnarray} \tag{ 17 }$$

As in classical mechanics, this Lagrangian is strictly equivalent to the Hamiltonian

$$\begin{eqnarray}\begin{array}{rcl}H & = & {\boldsymbol{U}}\cdot {\boldsymbol{p}}-{\ell }=\displaystyle \frac{1}{2}{g}^{\mu \nu }{p}_{\mu }{p}_{\nu }\\ & = & F+\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}}{{\rm{\Sigma }}}{p}_{r}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{{p}_{\theta }^{2}}{{\rm{\Sigma }}},\end{array}\end{eqnarray} \tag{ 18 }$$

where F is a function of r and θ and reads

$$\begin{eqnarray}\begin{array}{rcl}F & = & \displaystyle \frac{1}{2}({g}^{{tt}}{E}^{2}+{g}^{\phi \phi }{L}^{2})-{g}^{t\phi }{EL}\\ & = & \displaystyle \frac{1}{2}({f}_{1}{E}^{2}+{f}_{3}{L}^{2})+{f}_{2}{EL}\\ & = & -\displaystyle \frac{{{AE}}^{2}}{2{\rm{\Delta }}{\rm{\Sigma }}}+\displaystyle \frac{{L}^{2}({\rm{\Sigma }}-2r)}{2{\rm{\Delta }}{\rm{\Sigma }}{\sin }^{2}\theta }+\displaystyle \frac{2{ar}}{{\rm{\Delta }}{\rm{\Sigma }}}{EL}.\end{array}\end{eqnarray} \tag{ 19 }$$

The motion of the particle around the Kerr black hole is determined by the Hamiltonian (18) with two degrees of freedom in a four-dimensional phase space. Aside from the two constants E and L, the Hamiltonian itself is always equal to a known constant:

$$\begin{eqnarray}&&H=-\displaystyle \frac{1}{2}.\end{eqnarray} \tag{ 20 }$$

This constant is due to the four-velocity relation, Equation (9).

It is clear that the Hamiltonian system (18) has two degrees of freedom with four phase-space variables $({p}_{r},{p}_{\theta };r,\theta )$. Aside from the Hamiltonian itself as an integral, a second integral ⁵ can be found by use of the separation of variables in the Hamilton–Jacobi equations. This is the so-called Carter constant (Carter 1968), which is given in Appendix A. Thus, the system H is integrable and formally analytically solvable. In spite of this, no elementary functions but quadratures can be given to the formal solutions. In this case, it is still necessary to employ good numerical methods to study the long-term evolution of geodesics of the Kerr geometry.

In the previous papers (Wang et al. 2021a, 2021b, 2021c), we designed explicit symplectic integrators for the Hamiltonians of nonrotating black holes like the Reissner–Nordström black hole. This successful construction is completely based on the splitting parts of each Hamiltonian that exists as analytical solutions of explicit functions of proper time τ. If a splitting technique similar to that of the Hamiltonian of a Reissner–Nordström black hole in Paper II (Wang et al. 2021b) is used on the Hamiltonian (18), we have

$$\begin{eqnarray}\begin{array}{rcl}H & = & {H}_{1}+{H}_{2}+{H}_{3}+{H}_{4}+{H}_{5},\\ {H}_{1} & = & F,\\ {H}_{2} & = & \displaystyle \frac{1}{2}\displaystyle \frac{{r}^{2}}{{\rm{\Sigma }}}{p}_{r}^{2},\\ {H}_{3} & = & \displaystyle \frac{1}{2}\displaystyle \frac{{a}^{2}}{{\rm{\Sigma }}}{p}_{r}^{2},\\ {H}_{4} & = & -\displaystyle \frac{r}{{\rm{\Sigma }}}{p}_{r}^{2},\\ {H}_{5} & = & \displaystyle \frac{1}{2}\displaystyle \frac{{p}_{\theta }^{2}}{{\rm{\Sigma }}}.\end{array}\end{eqnarray} \tag{ 21 }$$

H₁ is easily solved. However, it is difficult to give analytical solutions to any one of H₂, H₃, H₄, and H₅ because none of the Hamiltonians are separable to the variables. Even if these sub-Hamiltonians can be solved analytically, their solutions are not explicit functions of proper time τ. For example, the evolution of r with τ in the sub-Hamiltonian H₃ is described by

$$\begin{eqnarray}&&\displaystyle \frac{r}{2}\sqrt{{\rm{\Sigma }}}+\displaystyle \frac{{a}^{2}{\cos }^{2}\theta }{2}\mathrm{ln}(r+\sqrt{{\rm{\Sigma }}})={c}_{1}{a}^{2}\tau +{c}_{2},\end{eqnarray} \tag{ 22 }$$

where c₁ and c₂ are integral constants. Because r is only one implicit function of τ, the splitting Hamiltonian method fails to construct an explicit symplectic integrator. In fact, this failure is directly due to Σ in the denominators of the p_r and ${p}_{\theta }$ terms acting as a function of r and θ. To successfully construct explicit symplectic integrators for the above-mentioned Hamiltonian, we must eliminate the function Σ in the denominators by constructing a time-transformed Hamiltonian like Equation (3).

In the present case, we take the time transformation (2) in the form

$$\begin{eqnarray}&&d\tau =g(r,\theta ){dw},\end{eqnarray} \tag{ 23 }$$

where τ is still the proper time and w is a new coordinate time unlike the original coordinate time t. Now the proper time τ is referred to as a coordinate ${q}_{0}=\tau$, and its corresponding momentum is p₀. Note that ${p}_{0}\ne {p}_{t}$ but ${p}_{0}=-H=1/2$. Then, the original phase-space variables $({p}_{r},{p}_{\theta };r,\theta )$ are extended to a set of new phase-space variables $({p}_{0},{p}_{r},{p}_{\theta };{q}_{0},r,\theta )$. Similar to Equation (3), a time-transformed Hamiltonian is given to the Hamiltonian (18) in the form

$$\begin{eqnarray}&&{\mathscr{H}}=g(H+{p}_{0}).\end{eqnarray} \tag{ 24 }$$

When the time transformation function in Equation (23) takes

$$\begin{eqnarray}&&g(r,\theta )=\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}},\end{eqnarray} \tag{ 25 }$$

the time-transformed Hamiltonian in Equation (24) is

$$\begin{eqnarray}&&{\mathscr{H}}=\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}}(F+{p}_{0})+\displaystyle \frac{{\rm{\Delta }}}{2{r}^{2}}{p}_{r}^{2}+\displaystyle \frac{1}{2{r}^{2}}{p}_{\theta }^{2}.\end{eqnarray} \tag{ 26 }$$

The denominators in the second and third terms of the new Hamiltonian ${\mathscr{H}}$ just eliminate the function Σ, compared with those in the original Hamiltonian H of Equation (18).

An operator-splitting technique is easily given to the Hamiltonian ${\mathscr{H}}$. However, it is unlike that of the Hamiltonian Γ with two separable integrable parts in Equation (3). The separable form of ${\mathscr{H}}$ is almost the same as that of the Hamiltonian of Reissner–Nordström black hole in Paper II. The Hamiltonian (26) takes five splitting parts

$$\begin{eqnarray}&&{\mathscr{H}}={{\mathscr{H}}}_{1}+{{\mathscr{H}}}_{2}+{{\mathscr{H}}}_{3}+{{\mathscr{H}}}_{4}+{{\mathscr{H}}}_{5},\end{eqnarray} \tag{ 27 }$$

where the sub-Hamiltonians ${{\mathscr{H}}}_{i}\,(i=1,\ldots ,5)$ are

$$\begin{eqnarray}&&{{\mathscr{H}}}_{1}=\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}}(F+{p}_{0}),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{2}=\displaystyle \frac{1}{2}{p}_{r}^{2},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{3}=-\displaystyle \frac{1}{r}{p}_{r}^{2},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{4}=\displaystyle \frac{{a}^{2}}{2{r}^{2}}{p}_{r}^{2},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{5}=\displaystyle \frac{1}{2{r}^{2}}{p}_{\theta }^{2}.\end{eqnarray} \tag{ 32 }$$

The equations of motion for the sub-Hamiltonian ${{\mathscr{H}}}_{1}$ in the new coordinate time read ${{dp}}_{0}/{dw}={dr}/{dw}=d\theta /{dw}=0$, and

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\tau }{{dw}} & = & \displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}}=g(r,\theta ),\\ \displaystyle \frac{{{dp}}_{r}}{{dw}} & = & \displaystyle \frac{2}{r}\left(\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}}-1\right)(F+{p}_{0})-\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}}\displaystyle \frac{\partial F}{\partial r}={P}_{r}(r,\theta ),\\ \displaystyle \frac{{{dp}}_{\theta }}{{dw}} & = & \displaystyle \frac{{a}^{2}}{{r}^{2}}(F+{p}_{0})\sin (2\theta )-\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}}\displaystyle \frac{\partial F}{\partial \theta }={P}_{\theta }(r,\theta ).\end{array}\end{eqnarray} \tag{ 33 }$$

The equations of motion for the other sub-Hamiltonians are written as

$$\begin{eqnarray}&&{{\mathscr{H}}}_{2}:\displaystyle \frac{{dr}}{{dw}}={p}_{r},\,\,{\dot{p}}_{r}=0;\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{3}:\displaystyle \frac{{dr}}{{dw}}=-\displaystyle \frac{2}{r}{p}_{r},\,\displaystyle \frac{{{dp}}_{r}}{{dw}}=-\displaystyle \frac{{p}_{r}^{2}}{{r}^{2}};\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{4}:\displaystyle \frac{{dr}}{{dw}}=\displaystyle \frac{{a}^{2}}{{r}^{2}}{p}_{r},\,\displaystyle \frac{{{dp}}_{r}}{{dw}}=\displaystyle \frac{{a}^{2}}{{r}^{3}}{p}_{r}^{2};\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{5}:\displaystyle \frac{d\theta }{{dw}}=\displaystyle \frac{{p}_{\theta }}{{r}^{2}},\,\displaystyle \frac{{{dp}}_{r}}{{dw}}=\displaystyle \frac{{p}_{\theta }^{2}}{{r}^{3}},\,\dot{r}={\dot{p}}_{\theta }=0.\end{eqnarray} \tag{ 37 }$$

Each equation is independently solved in an analytical method. Let $({r}_{0},{\theta }_{0},{p}_{r0},{p}_{\theta 0})$ be the values at the beginning of one step and $(r,\theta ,{p}_{r},{p}_{\theta })$ denote the analytical solutions at the end of the step over time w. Analytical solutions are easily given for the five pieces, unlike those in Equation (4). They are labeled as operators ${\tilde{{\mathscr{H}}}}_{1}(w)$, ${\tilde{{\mathscr{H}}}}_{2}(w)$, ${\tilde{{\mathscr{H}}}}_{3}(w)$, ${\tilde{{\mathscr{H}}}}_{4}(w)$, and ${\tilde{{\mathscr{H}}}}_{5}(w)$. In fact, they are explicit functions of the new coordinate time w:

$$\begin{eqnarray}\begin{array}{rcl}\tau (w) & = & {\tau }_{0}+{wg}({r}_{0},{\theta }_{0}),\\ {\tilde{{\mathscr{H}}}}_{1}:\,{p}_{r}(w) & = & {p}_{r0}+{{wP}}_{r}({r}_{0},{\theta }_{0}),\end{array}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}\begin{array}{rcl}{p}_{\theta }(w) & = & {p}_{\theta 0}+{{wP}}_{\theta }({r}_{0},{\theta }_{0});\\ {\tilde{{\mathscr{H}}}}_{2}:\,r(w) & = & {r}_{0}+{{wp}}_{r0};\end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{{\mathscr{H}}}}_{3}:\,r(w) & = & {\left[{\left({r}_{0}^{2}-3{{wp}}_{r0}\right)}^{2}/{r}_{0}\right]}^{\displaystyle \frac{1}{3}},\\ {p}_{r}(w) & = & {p}_{r0}{\left[({r}_{0}^{2}-3{{wp}}_{r0})/{r}_{0}^{2}\right]}^{\displaystyle \frac{1}{3}};\end{array}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{{\mathscr{H}}}}_{4}:\,r(w) & = & {\left(2{a}^{2}{p}_{r0}w/{r}_{0}+{r}_{0}^{2}\right)}^{\displaystyle \frac{1}{2}},\\ {p}_{r}(w) & = & \displaystyle \frac{{p}_{r0}}{{r}_{0}}{\left(2{a}^{2}{p}_{r0}w/{r}_{0}+{r}_{0}^{2}\right)}^{\displaystyle \frac{1}{2}};\end{array}\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{{\mathscr{H}}}}_{5}:\,\theta (w) & = & {\theta }_{0}+{{wp}}_{\theta 0}/{r}_{0}^{2},\\ {p}_{r}(w) & = & {p}_{r0}+{{wp}}_{\theta 0}^{2}/{r}_{0}^{3}.\end{array}\end{eqnarray} \tag{ 42 }$$

Any one of the three compositions, involving Equations (29) and (30), Equations (29)–(31), and Equations (29)–(32), is analytically solvable. However, the analytical solutions are implicit functions of w.

Using the splitting operators with a new coordinate time step h, we design an explicit second-order symplectic integrator for the system ${\mathscr{H}}$

$$\begin{eqnarray}\begin{array}{rcl}{S}_{2}^{{\mathscr{H}}}(h) & = & {\tilde{{\mathscr{H}}}}_{5}\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{4}\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{3}\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{2}\left(\displaystyle \frac{h}{2}\right)\\ & & \times \ \circ {\tilde{{\mathscr{H}}}}_{1}(h)\circ {\tilde{{\mathscr{H}}}}_{2}\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{3}\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{4}\left(\displaystyle \frac{h}{2}\right)\\ & & \times \ \circ {\tilde{{\mathscr{H}}}}_{5}\left(\displaystyle \frac{h}{2}\right).\end{array}\end{eqnarray} \tag{ 43 }$$

Its detailed expression is given in Appendix B.

A symmetric composition of three second-order methods yields the fourth-order explicit symplectic scheme of Yoshida (1990),

$$\begin{eqnarray}&&{S}_{4}^{{\mathscr{H}}}(h)={S}_{2}^{{\mathscr{H}}}(\gamma h)\circ {S}_{2}^{{\mathscr{H}}}(\delta h)\circ {S}_{2}^{{\mathscr{H}}}(\gamma h),\end{eqnarray} \tag{ 44 }$$

where $\delta =1-2\gamma$ and $\gamma =1/(2-\sqrt[3]{2})$.

Seen from the construction of the two explicit symplectic integrators for the Hamiltonian ${\mathscr{H}}$, the time transformation function g plays an important role in successfully eliminating Σ in the denominators of the p_r and ${p}_{\theta }$ terms of the inseparable Hamiltonian H. This is successful to overcome an obstacle to the application of such explicit symplectic integrators to the Hamiltonian ${\mathscr{H}}$. Because $g=1+{a}^{2}{\cos }^{2}\theta /{r}^{2}\leqslant 1\,+{a}^{2}/{r}^{2}\leqslant 1+1/{r}^{2}\approx 1$, we have ${\rm{\Delta }}\tau \sim g{\rm{\Delta }}w\approx {\rm{\Delta }}w=h$ in Equation (23). ⁶ When the fixed time step h is given to the coordinate time w, the proper time step ${\rm{\Delta }}\tau$ slightly depends on the radial distance r and even is almost the same fixed coordinate time step h. The symplectic structure of the time-transformed Hamiltonian flow ${\mathscr{H}}$ is preserved for the fixed coordinate time step ${\rm{\Delta }}w$ but is not for the slightly variable proper time step ${\rm{\Delta }}\tau$. A large difference between the two sets of time steps may exist when the time transformation function g is altered. Other choices of time transformation function are given in Appendix C.

In the above-mentioned three time transformation functions, the time transformation functions g₁ and g₂ play important roles in step-size selection procedures with adaptive proper time steps when an invariant coordinate time step is adopted. The choice of g₁ exerts a larger influence on the adjustment of proper time steps than that of g but has a smaller influence than that of g₂. Here, we are mainly interested in applying a time transformation to implement the construction of explicit symplectic integrators for the Kerr spacetime geometry. Considering this purpose, we take the time transformation function g in the following discussions.

Compared with the newly proposed integrators, several existing numerical methods are considered. They are RK4, fourth-order implicit and explicit mixed symplectic algorithm (IE4), and the fourth-order extended phase-space explicit symplectic-like method (EP4).

5.1. Integrating the Original Hamiltonian H with the Existing Algorithms

The sum of the second and third terms in Equation (18) is labeled as K. It is solved in terms of the second-order implicit midpoint rule (Feng 1986). Its corresponding operator is ${IM}2(h)$, where h is a proper time step. F in Equation (18) is easily solvable and corresponds to an operator ${\psi }_{h}^{F}$. The two operators symmetrically compose a second-order implicit and explicit mixed symplectic integrator for the original Hamiltonian H in Equation (18):

$$\begin{eqnarray}&&{IE}2(h)={\psi }_{h/2}^{K}\circ {IM}2(h)\circ {\psi }_{h/2}^{K}.\end{eqnarray} \tag{ 45 }$$

The computational efficiency of the method IE2 acting on H is superior to that of the algorithm IM2 acting on H. More details on the implicit and explicit mixed symplectic methods were given in the references (Lubich et al. 2010; Zhong et al. 2010; Mei et al. 2013a, 2013b). Similar to the algorithm ${S}_{4}^{{\mathscr{H}}}(h)$ in Equation (44), a fourth-order implicit and explicit mixed symplectic algorithm is expressed as

$$\begin{eqnarray}&&{IE}4={IE}2(\gamma h)\circ {IE}2(\delta h)\circ {IE}2(\gamma h).\end{eqnarray} \tag{ 46 }$$

On the other hand, the four-dimensional phase-space variables $(r,\theta ,{p}_{r},{p}_{\theta })$ of H are extended to eight-dimensional phase-space variables $(r,\theta ,\tilde{r},\tilde{\theta },{p}_{r},$ ${p}_{\theta },{\tilde{p}}_{r},{\tilde{p}}_{\theta })$ in a new Hamiltonian,

$$\begin{eqnarray}&&{\mathbb{H}}={{\mathbb{H}}}_{1}(r,\theta ,{\tilde{p}}_{r},{\tilde{p}}_{\theta })+{{\mathbb{H}}}_{2}(\tilde{r},\tilde{\theta },{p}_{r},{p}_{\theta }),\end{eqnarray} \tag{ 47 }$$

where ${{\mathbb{H}}}_{1}(r,\theta ,{\tilde{p}}_{r},{\tilde{p}}_{\theta })={{\mathbb{H}}}_{2}(\tilde{r},\tilde{\theta },{p}_{r},{p}_{\theta })=H(r,\theta ,{p}_{r},{p}_{\theta })$. ${{\mathbb{H}}}_{1}$ and ${{\mathbb{H}}}_{2}$ are independently solved in analytical methods. They correspond to operators ${\mathscr{A}}$ and ${\mathscr{B}}$ similar to those in Section 2. The second- and fourth-order methods ${{\mathscr{S}}}_{2}^{{\mathbb{H}}}$ and ${{\mathscr{S}}}_{4}^{{\mathbb{H}}}$ like those in Equations (5) and (6) can be obtained. For the same initial conditions, the two independent Hamiltonians should have the same solutions: $r=\tilde{r}$, $\theta =\tilde{\theta }$, ${\tilde{p}}_{r}={p}_{r}$, and ${\tilde{p}}_{\theta }={p}_{\theta }$. However, the two solutions are not the same due to their couplings in ${{\mathscr{S}}}_{2}^{{\mathbb{H}}}$ and ${{\mathscr{S}}}_{4}^{{\mathbb{H}}}$. Permutations between the original variables and their corresponding extended variables are necessarily used after the implementation of ${{\mathscr{S}}}_{2}^{{\mathbb{H}}}$ and ${{\mathscr{S}}}_{4}^{{\mathbb{H}}}$ so as to frequently make the two solutions equal (Pihajoki 2015; Liu et al. 2016; Luo et al. 2017; Li & Wu 2017; Wu & Wu 2018). The midpoint permutation method (Luo et al. 2017)

$$\begin{eqnarray}\begin{array}{rcl}{\mathscr{M}}:\displaystyle \frac{r+\tilde{r}}{2} & \to & r=\tilde{r},\\ \displaystyle \frac{\theta +\tilde{\theta }}{2} & \to & \theta =\tilde{\theta };\\ \displaystyle \frac{{p}_{r}+{\tilde{p}}_{r}}{2} & \to & {p}_{r}={\tilde{p}}_{r},\\ \displaystyle \frac{{p}_{\theta }+{\tilde{p}}_{\theta }}{2} & \to & {p}_{\theta }={\tilde{p}}_{\theta }\end{array}\end{eqnarray} \tag{ 48 }$$

is a good choice. The algorithms ${{\mathscr{S}}}_{2}^{{\mathbb{H}}}$ and ${{\mathscr{S}}}_{4}^{{\mathbb{H}}}$ combined with the midpoint permutation are

$$\begin{eqnarray}&&{EP}2={\mathscr{M}}\otimes {{\mathscr{S}}}_{2}^{{\mathbb{H}}},\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{EP}4={\mathscr{M}}\otimes {{\mathscr{S}}}_{4}^{{\mathbb{H}}}.\end{eqnarray} \tag{ 50 }$$

The symplecticity of ${{\mathscr{S}}}_{2}^{{\mathbb{H}}}$ and ${{\mathscr{S}}}_{4}^{{\mathbb{H}}}$ is destroyed due to the inclusion of ${\mathscr{M}}$. However, the symmetry makes the methods IE2 and IE4 exhibit good long-term stable behavior in Hamiltonian errors. In this sense, the methods IE2 and IE4 are explicit symplectic-like algorithms for the newly extended phase-space Hamiltonian ${\mathbb{H}}$ in Equation (47).

Taking the parameters $E=0.995$, $L=4.6$, and $a=0.5$, we choose a test orbit with the initial conditions r = 11, $\theta =\pi /2$, and ${p}_{r}=0$. The starting value of ${p}_{\theta }\gt 0$ is given by Equation (17). The proper time step is h = 1 and the integration time arrives at 10⁸. As shown in Figure 1, RK4 shows a secular drift in Hamiltonian error ΔH = −1 – 2H in Equation (20). So does the second-order extended phase-space method EP2. However, this drift is absent in EP2* when the proper time step h = 1 is replaced with $h=0.1$. It is because the extended phase-space method maintaining the boundness of errors requires appropriate permutations. When the step size becomes smaller, more permutations are used. This results in the improvement of the algorithm's accuracy. In fact, the accuracy of EP2* has an order of ${\mathscr{O}}({10}^{-8})$ and is two orders of magnitude higher than that of EP2. The fourth-order extended phase-space method EP4 with the proper time step h = 1 also gives the same order errors and therefore makes the errors bounded. If the Hamiltonian truncation errors are not larger than the order of ${\mathscr{O}}({10}^{-9})$, 10⁸ integration steps yield a great number of round-off errors that cause the Hamiltonian errors to have secular growth. ⁷ The fourth-order implicit and explicit mixed symplectic algorithm IE4 is suitable for this case. When a larger step size h = 4 is used, the Hamiltonian errors for IE4* are similar to those for the second-order implicit and explicit mixed symplectic algorithm IE2 with the proper time step h = 1 and remain stable and bounded.

Zoom In Zoom Out Reset image size

In short, the algorithms' performance tested in the above Kerr geometry is very similar to that checked in the Schwarzschild spacetime in Paper I (Wang et al. 2021a).

5.2. Solving the Time-transformed Hamiltonian ${\mathscr{H}}$

Without question, the algorithms RK4, IE4, and EP4 are also suited for solving the time-transformed Hamiltonian ${\mathscr{H}}$ in Equation (26). In addition to them, the new explicit symplectic methods ${S}_{2}^{{\mathscr{H}}}$ (labeled S2) and ${S}_{4}^{{\mathscr{H}}}$ (labeled S4) are respectively used to solve the Hamiltonian splitting form (27). Now, h is a coordinate time step. It is shown by comparing Figures 1 and 2(a) and (b) that the existing algorithms for the time-transformed Hamiltonian ${\mathscr{H}}$ and those for the original Hamiltonian H in Equation (18) have no typical differences in the behavior of the Hamiltonian error. Methods S2 and S4 exhibit good long-term stabilized error behavior in Figure 2(c). No secular error change exists for S4 because S4 gives the order of ${\mathscr{O}}({10}^{-8})$ to the Hamiltonian errors rather than the order of ${\mathscr{O}}({10}^{-9})$.

Zoom In Zoom Out Reset image size

The results can be concluded from Figure 2. For the coordinate time step h = 1, the algorithms RK4, EP2, IE2, and S2 are approximately the same in Hamiltonian errors; the fourth-order methods EP4, IE4, and S4 are, too. The schemes without secular error drifts are IE2, S2, EP4, and S4. In addition, EP2* with a smaller coordinate time step $h=0.1$ approaches S4 in accuracy. IE4* and S4* with a larger coordinate time step h = 4 are close to S2. EP2*, IE*, and S4* maintain the boundness of the Hamiltonian errors. These similar results are still present when the dynamical parameters and orbits are altered.

Besides the Hamiltonian ${\mathscr{H}}$, the Carter constant K in Appendix A can be conserved by the algorithms S4, EP4, and IE4, but cannot be conserved by the RK4 method, as shown in Figure 3(a). RK4 is the poorest in the accuracy of the Carter constant, and EP4 is the best. IE4 is slightly better than S4. However, RK4 is the best in the accuracy of the solution for an integration time of w = 2000 in Figure 3(b), and S4 is better than EP4 or IE4. This shows that RK4 is accurate and efficient for such a short time integration.

Zoom In Zoom Out Reset image size

Testing the performance of these algorithms in Figures 1–3 is based on the boundness of massive particle orbits. What about the performance of these algorithms if the integrated orbits are unstable and unbounded? To answer this question, we adopt the unstable spherical photon orbit B with normalized angular momentum Φ = −6, normalized Carter constant $Q=-13+16\sqrt{2}$, and constant radius ${r}_{0}=1+2\sqrt{2}$ considered by Bacchini et al. (2018). Some details of the unstable spherical photon orbits are described in Appendix A. When w = 116 in Figure 4, the orbit begins to run to infinity for RK4, and the accuracies of the Hamiltonian ${\mathscr{H}}$ and Carter constant K become the worst. Although this orbit does not remain spherical due to orbital instability, S4, IE4, and EP4 can still work well in the conservation of the two constants. In fact, S4 exhibits the best accuracy.

Zoom In Zoom Out Reset image size

Apart from the accuracies of the constants and solutions for these algorithms, computational efficiency should be compared. Table 1 lists the CPU times of the algorithms in Figures 1 and 2. S2 has the best efficiency among the three second-order methods S2, IE2, and EP2. The efficiency of S4 is also the best among the efficiencies of the three fourth-order methods S4, IE4, and EP4. This shows the superiority of the application of explicit algorithms in computational efficiency.

Table 1. Computational Cost (CPU Times in Minute' Second'') for the Algorithms in Figures 1 and 2

Method	RK4	EP2	EP2*	IE2	EP4	IE4	IE4*	S2	S4	S4*
Figure 1	2'4''	1'33''	15'45''	1'53''	3'37''	5'48''	1'33''	/	/	/
Figure 2	1'51''	1'21''	13'43''	1'43''	3'8''	5'26''	1'27''	57''	2'50''	43''

Download table as:  ASCIITypeset image

Figure 5(a) plots the evolution curve of r with proper time τ (colored black dot), obtained by IE4 solving the original Hamiltonian system (18) with the proper time step h = 1. It also plots the evolution of r with coordinate time w (colored red line), given by S4 integrating the time-transformed Hamiltonian (26) with the coordinate time step h = 1. The two evolution curves coincide. Here is an explanation. The proper time steps are varied by using a fixed coordinate time step. However, the coordinate time w and the proper time τ have no typical differences when S4 solves the Hamiltonian (26) with the coordinate time step, as shown in Figure 5(b). Clearly, the time transformation function in Equation (25) mainly plays an important role in eliminating the function Σ in the denominators of the Hamiltonian (26). It then leads to successfully implementing the construction of explicit symplectic integrators for the time-transformed Hamiltonian (26). The time transformation function does not have any explicit effect on the step-size selection procedures.

Zoom In Zoom Out Reset image size

5.3. Discussions

The new explicit symplectic integrators are applicable for not only integrating geodesics in the Kerr spacetime, but also a Kerr black hole immersed in an external magnetic field (Kopáček et al. 2010; Kopáček & Karas 2014, 2018), Kerr–Newman black hole, and Kerr–Newman black hole with an external magnetic field. The construction of explicit symplectic integrators is based on the splitting of the time-transformed Hamiltonian when the denominators of the contravariant metric components ${g}^{{rr}}$ and ${g}^{\theta \theta }$ have Σ as a function of r and θ, and the time transformation function $g={\rm{\Sigma }}/{r}^{2}$ is taken. Such similar constructions are still possibly available if ${g}^{{rr}}$ and ${g}^{\theta \theta }$ have other expressions. In this case, the time transformation function g will need an appropriate choice. In fact, the splitting is valid if the time-transformed Hamiltonian becomes

$$\begin{eqnarray}&&{\mathscr{H}}={K}_{1}(r,\theta )+{{gg}}^{{rr}}{p}_{r}^{2}+{{gg}}^{\theta \theta }{p}_{\theta \theta }^{2},\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{gg}}^{{rr}} & = & \left({a}_{0}+{a}_{1}r\,+\cdots +\,{a}_{n}{r}^{n}+\,\displaystyle \frac{{a}_{-1}}{r}+\displaystyle \frac{{a}_{-2}}{{r}^{2}}\right.\\ & & \left.+\ \cdots +\,\displaystyle \frac{{a}_{-m}}{{r}^{m}}\right){K}_{2}(\theta ),\end{array}\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{gg}}^{\theta \theta } & = & \left({b}_{0}+{b}_{1}\theta \,+\cdots +\,{b}_{i}{\theta }^{i}+\displaystyle \frac{{b}_{-1}}{\theta }+\displaystyle \frac{{b}_{-2}}{{\theta }^{2}}\right.\\ & & \left.+\ \cdots +\,\displaystyle \frac{{b}_{-j}}{{r}^{j}}\right){K}_{3}(r),\end{array}\end{eqnarray} \tag{ 53 }$$

where K₁, K₂, and K₃ are arbitrary continuous differentiable functions when r and θ are restricted to $r\gt 2$ and $\theta \ne 0$, and a₀, ⋯, a_n , ${a}_{-1}$, ⋯, ${a}_{-m}$, b₀, ⋯, b_i , ${b}_{-1}$, ⋯, ${b}_{-j}$ are constant parameters. It is clear that the time-transformed Hamiltonian splitting is desired when the time transformation function g satisfies Equations (52) and (53). For example, the time transformation function takes $g={e}^{Q-P}$ for ${g}^{{rr}}={e}^{P-Q}(1\,-2/r)$ and ${g}^{\theta \theta }={e}^{P-Q}{r}^{-2}$, where P and Q are functions of r and θ in an axially symmetric core–shell system (Vieira & Letelier 1999). Another example is a five-dimensional black ring metric (Igata et al. 2011), where ${g}^{{xx}}={(x-y)}^{2}(1\,-{x}^{2})(1+\nu x)/[{R}^{2}(1+\lambda x)]$ and ${g}^{{yy}}=-{(x-y)}^{2}(1-{y}^{2})(1\,+\nu y)/[{R}^{2}(1+\lambda x)]$ in ring coordinates $(t,x,y,\phi ,\psi )$. Here, $| x| \leqslant 1$ and y ≤ −1 are coordinate variables; $R\gt 0$ and $0\lt \nu \leqslant \lambda \lt 1$ are constant parameters; and t, ϕ, and ψ do not explicitly appear in the metric. Obviously, the time transformation function $g=1+\lambda x$ is a good choice. Therefore, the idea for the construction of the present explicit symplectic integrators is not restricted to Kerr-type spacetimes and is applied to many other relativistic problems.

Although the application of the new algorithms is not wider than that of RAPTOR based on the RK4 scheme (Bronzwaer et al. 2018), the accuracies of the new algorithms have an advantage over those of the RK4 scheme in long-term integrations. The new second-order explicit symplectic integrator S2 and RK4 have the same order of magnitude in the Hamiltonian errors in Figures 2(a) and (c), but they have typical differences in simulations of unstable spherical photon orbits or long-term evolution of massive particle bounded orbits. The errors grow with the increase of integration time for RK4 but remain bounded for S2. When the integration time reaches 116, RK4 gives extremely larger errors for the constants for the unstable spherical photon orbit, but S4 shows smaller errors. When the integration time lasts 10¹⁰, RK4 does not work well and provides incorrect results for the massive particle bounded orbits, whereas S2 still shows no secular growth in the errors. This is an advantage of a symplectic scheme in long-term integrations. In particular, Deng et al. (2020) reported that RK4 combined with manifold correction exhibits poorer performance than the second-order symplectic leapfrog in a 10⁸ yr integration of the outer solar system involving the Sun and four outer planets. In addition, Figures 2(a) and (c) clearly show that the new fourth-order explicit symplectic integrator S4 is two orders of magnitude higher than RK4 in accuracy. S4 is typically better than RK4 in numerical performance for the integration of unstable spherical photon orbits.

As shown in the table, the newly proposed algorithms are greatly superior to the implicit and explicit mixed symplectic methods without a doubt in terms of computational efficiency. They also have an advantage over the extended phase-space explicit symplectic-like methods EP2 and EP4 because they integrate over four dimensions, while the extended phase-space schemes require doubling the phase space (i.e., eight dimensions). Figures 2(a) and (c) describe how the extended phase-space method EP2 possesses larger errors than the proposed algorithm S2. Figure 3(b) shows that EP4 is poorer than S4 in the accuracy of the solutions. The methods EP2 and EP4 are in fact similar to the method of Tao (2016) or FANTASY (Christian & Chan 2021). Wu & Wu (2018) found that fourth-order methods like EP4 with midpoint permutations show better accuracies than Tao's method or FANTASY with a good choice of the control constant.

The newly proposed algorithms are inferior to the energy-conserving discrete gradient scheme of Bacchini et al. (2018) in energy conservation but have some explicit advantages over the latter method. The latter algorithms are implicit, nonsymplectic, and do not preserve other integrals in general (perhaps some integrals may be conserved, as claimed in footnote 1), but the proposed algorithms are explicit, symplectic, and preserve other integrals like the Carter constant. Therefore, the proposed algorithms are extremely superior to the latter algorithm in computational efficiency. In addition, fourth-, sixth-, and higher-order schemes are easily available according to the present idea, whereas they are not in terms of the construction method of Bacchini et al.

The function Σ depending on r and θ exists in the denominators of the Hamiltonian of Kerr geometry. Thus, explicit symplectic integrators are useless if the Hamiltonian is split into several parts like the Hamiltonian of Schwarzschild spacetimes in previous papers (Wang et al. 2021a, 2021b, 2021c).

To overcome this difficulty, we use an appropriate time transformation function to eliminate the function Σ in the denominators and to obtain a time-transformed Hamiltonian. This Hamiltonian can be separated into five parts with analytical solutions as explicit functions of the new coordinate time, and therefore explicit symplectic integrators are easily available. Numerical tests show that the explicit symplectic methods exhibit good long-term stable behavior in Hamiltonian errors when appropriate coordinate time steps are chosen. The use of fixed coordinate time steps maintains the symplecticity of the time-transformed Hamiltonian. Although the proper time sizes are variant, there is no explicit difference between the proper time and newly transformed coordinate time. In this sense, the chosen time transformation function does not mainly play a role in step-size selection procedures. The proposed algorithms are superior to the same-order extended phase-space explicit symplectic-like methods and implicit and explicit mixed symplectic algorithms in computational efficiency.

The new idea for the construction of such explicit symplectic integrators is not limited to the application of the Kerr metric and allows for the application of many spacetimes, given in Equation (51). This problem will be considered further in our future works.

The authors are very grateful to the referee for valuable comments and useful suggestions. This research has been supported by the National Natural Science Foundation of China (grant Nos. 11973020 (C0035736), 11533004, 11803020, 41807437, and U2031145), and the Natural Science Foundation of Guangxi (grant Nos. 2018GXNSFGA281007 and 2019JJD110006).

Noticing Equations (18)–(20), we rewrite Equation (18) as follows:

$$\begin{eqnarray}\begin{array}{rcl}-\chi {\rm{\Sigma }} & = & -\displaystyle \frac{{E}^{2}}{{\rm{\Delta }}}({\rho }^{4}-{\rm{\Delta }}{a}^{2}{\sin }^{2}\theta )+\displaystyle \frac{{L}^{2}({\rm{\Delta }}-{a}^{2}{\sin }^{2}\theta )}{{\rm{\Delta }}{\sin }^{2}\theta }\\ & & +\ \displaystyle \frac{4{ar}}{{\rm{\Delta }}}{EL}+{\rm{\Delta }}{p}_{r}^{2}+{p}_{\theta }^{2},\end{array}\end{eqnarray} \tag{ A1 }$$

where $\chi =1$ for massive particles and $\chi =0$ for massless photons. The equation has the separation of variables:

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\theta }^{2}+\displaystyle \frac{{L}^{2}}{{\sin }^{2}\theta }+{E}^{2}{a}^{2}{\sin }^{2}\theta +\chi {a}^{2}{\cos }^{2}\theta =-\chi {r}^{2}\\ +\displaystyle \frac{{E}^{2}}{{\rm{\Delta }}}{\rho }^{4}+\displaystyle \frac{{L}^{2}{a}^{2}}{{\rm{\Delta }}}-\displaystyle \frac{4{ar}}{{\rm{\Delta }}}{EL}-{\rm{\Delta }}{p}_{r}^{2}.\end{array}\end{eqnarray} \tag{ A2 }$$

Clearly, the left-hand side is a function of the variables θ and ${p}_{\theta }$, whereas the right-hand side is a function of the variables r and p_r . This equality is possible only when the two sides are equal to the same constant

$$\begin{eqnarray}&&K={p}_{\theta }^{2}+\left[\displaystyle \frac{{L}^{2}}{{\sin }^{2}\theta }+{a}^{2}(\chi -{E}^{2})\right]{\cos }^{2}\theta ,\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}\begin{array}{rcl}K & = & -{L}^{2}-{a}^{2}{E}^{2}-\chi {r}^{2}+\displaystyle \frac{{E}^{2}}{{\rm{\Delta }}}{\rho }^{4}+\displaystyle \frac{{L}^{2}{a}^{2}}{{\rm{\Delta }}}\\ & & -\ \displaystyle \frac{4{ar}}{{\rm{\Delta }}}{EL}-{\rm{\Delta }}{p}_{r}^{2}.\end{array}\end{eqnarray} \tag{ A4 }$$

Equation (A3) or (A4) is an expressional form of the Carter constant.

Equation (A4) is rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}^{2}{\left(\displaystyle \frac{{dr}}{d\tau }\right)}^{2} & = & -(K+{L}^{2}+{a}^{2}{E}^{2}+\chi {r}^{2}){\rm{\Delta }}\\ & & +\ {E}^{2}{\rho }^{4}+{L}^{2}{a}^{2}-4{arEL}=R(r).\end{array}\end{eqnarray} \tag{ A5 }$$

When $R(r)={dR}(r)/{dr}=0$ in Equation (A5), r corresponds to a constant radius of the spherical photon orbit. Teo (2003) studied such spherical photon orbits with normalized angular momentum Φ and normalized Carter constant Q:

$$\begin{eqnarray}&&{\rm{\Phi }}=\displaystyle \frac{L}{E}=-\displaystyle \frac{{r}^{3}-3{r}^{2}+{a}^{2}(r+1)}{a(r-1)},\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&Q=\displaystyle \frac{K}{{E}^{2}}=-\displaystyle \frac{{r}^{3}({r}^{3}-6{r}^{2}+9r-4{a}^{2})}{{\left[a(r-1)\right]}^{2}}.\end{eqnarray} \tag{ A7 }$$

These two expressions were also provided by Chan et al. (2018). Because ${d}^{2}R(r)/{{dr}}^{2}\gt 0$, these orbits would be unstable when they suffer from perturbations in the radial direction. Müller & Grave (2010) gave the expression of E as follows:

$$\begin{eqnarray}&&E={\left(\sqrt{\displaystyle \frac{A}{{\rm{\Sigma }}{\rm{\Delta }}}}-\displaystyle \frac{2{ra}{\rm{\Phi }}}{\sqrt{A{\rm{\Sigma }}{\rm{\Delta }}}}\right)}^{-1}.\end{eqnarray} \tag{ A8 }$$

Thus, the other two parameters are easily determined by

$$\begin{eqnarray}&&L=E{\rm{\Phi }},\qquad K={{QE}}^{2}.\end{eqnarray} \tag{ A9 }$$

For given parameters r and a = 1 with the initial value $\theta =\pi /2$, we use Equation (A3) to obtain the initial value

$$\begin{eqnarray}&&{p}_{\theta }=E\sqrt{Q}.\end{eqnarray} \tag{ A10 }$$

Codes of the method in Equation (43) are written in the following forms:

$$\begin{eqnarray*}\begin{array}{rcl}{\theta }^{* } & = & {\theta }_{0}+\displaystyle \frac{h}{2}\displaystyle \frac{{p}_{\theta 0}}{{r}_{0}^{2}},\\ {p}_{r}^{* } & = & {p}_{r0}+\displaystyle \frac{h}{2}\displaystyle \frac{{p}_{\theta 0}^{2}}{{r}_{0}^{3}};\\ {r}^{\star } & = & {\left(\displaystyle \frac{h}{{r}_{0}}{a}^{2}{p}_{r}^{* }+{r}_{0}^{2}\right)}^{1/2},\\ {p}_{r}^{\star } & = & \displaystyle \frac{{p}_{r}^{* }}{{r}_{0}}{\left(\displaystyle \frac{h}{{r}_{0}}{a}^{2}{p}_{r}^{* }+{r}_{0}^{2}\right)}^{1/2};\\ {r}^{\dagger } & = & {\left[{\left({r}^{\star 2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{\star }\right)}^{2}/{r}^{\star }\right]}^{1/3},\\ {p}_{r}^{\dagger } & = & {p}_{r}^{\star }{\left[({r}^{\star 2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{\star })/{r}^{\star 2}\right]}^{1/3};\\ {r}^{\unicode{x02021}} & = & {r}^{\dagger }+\displaystyle \frac{h}{2}{p}_{r}^{\dagger };\\ {p}_{r}^{\sharp } & = & {p}_{r}^{\dagger }+{{hP}}_{r}({r}^{\unicode{x02021}},{\theta }^{* }),\\ {\tau }_{1} & = & {\tau }_{0}+{hg}({r}^{\unicode{x02021}},{\theta }^{* }),\\ {p}_{\theta 1} & = & {p}_{\theta 0}+{{hP}}_{\theta }({r}^{\unicode{x02021}},{\theta }^{* });\\ {r}^{\natural } & = & {r}^{\unicode{x02021}}+\displaystyle \frac{h}{2}{p}_{r}^{\sharp };\\ {r}^{\lozenge } & = & {\left[{\left({r}^{\natural 2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{\sharp }\right)}^{2}/{r}^{\natural }\right]}^{1/3},\\ {p}_{r}^{\lozenge } & = & {p}_{r}^{\sharp }{\left[{\left({r}^{\natural 2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{\sharp }\right)}^{2}/{r}^{\natural 2}\right]}^{1/3};\\ {p}_{r}^{\unicode{x02661}} & = & \displaystyle \frac{{p}_{r}^{\lozenge }}{{r}^{\lozenge }}{\left(\displaystyle \frac{h}{{r}^{\lozenge }}{a}^{2}{p}_{r}^{\lozenge }+{r}^{\lozenge 2}\right)}^{1/2},\\ {r}_{1} & = & {\left(\displaystyle \frac{h}{{r}^{\lozenge }}{a}^{2}{p}_{r}^{\lozenge }+{r}^{\lozenge 2}\right)}^{1/2};\\ {p}_{r1} & = & {p}_{r}^{\unicode{x02661}}+\displaystyle \frac{h}{2}\displaystyle \frac{{p}_{\theta 1}^{2}}{{r}_{1}^{3}},\\ {\theta }_{1} & = & {\theta }^{* }+\displaystyle \frac{h}{2}\displaystyle \frac{{p}_{\theta 1}}{{r}_{1}^{2}}.\end{array}\end{eqnarray*}$$

In this way, the solutions $({\tau }_{1},{r}_{1},{\theta }_{1},{p}_{r1},{p}_{\theta 1})$ are outputted after the values $({\tau }_{0},{r}_{0},{\theta }_{0},{p}_{r0},{p}_{\theta 0})$ advance a fixed coordinate time step h.

In principle, there are various choices for the time transformation function. For instance, it is given by

$$\begin{eqnarray}&&{g}_{1}=\displaystyle \frac{{\rm{\Sigma }}}{r}=r+\displaystyle \frac{1}{r}{a}^{2}{\cos }^{2}\theta .\end{eqnarray} \tag{ C1 }$$

The new coordinate time in Equation (23) becomes w₁. Advancing the fixed coordinate time ${\rm{\Delta }}{w}_{1}=h$ means advancing a variable proper time ${\rm{\Delta }}{\tau }_{1}\approx {rh}$. The advancement of proper time is absolutely dominated by the radial separation r. When a particle approaches the central black hole, ⁸ ${\rm{\Delta }}{\tau }_{1}$ is smaller, but it becomes larger when the particle moves always from the central object. This kind of step-size selection procedure uses adaptive time steps and satisfies the realistic need for the improvement of numerical precisions. An important contribution of using the time transformation equation is that this transformation brings the implementation of symplectic integrators with adaptive time-step controls. This is the main result claimed in the paper of Mikkola (1997). In terms of the time transformation function g₁, another new time-transformed Hamiltonian has four separable parts:

$$\begin{eqnarray}&&{{\mathscr{H}}}^{* }={{\mathscr{H}}}_{1}^{* }+{{\mathscr{H}}}_{2}^{* }+{{\mathscr{H}}}_{3}^{* }+{{\mathscr{H}}}_{4}^{* },\end{eqnarray} \tag{ C2 }$$

where these parts are

$$\begin{eqnarray}&&{{\mathscr{H}}}_{1}^{* }=\displaystyle \frac{{\rm{\Sigma }}}{r}(F+{p}_{0}),\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{2}^{* }=\left(\displaystyle \frac{r}{2}-1\right){p}_{r}^{2},\end{eqnarray} \tag{ C4 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{3}^{* }=\displaystyle \frac{{a}^{2}}{2r}{p}_{r}^{2},\end{eqnarray} \tag{ C5 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{4}^{* }=\displaystyle \frac{1}{2r}{p}_{\theta }^{2}.\end{eqnarray} \tag{ C6 }$$

Like the method ${S}_{2}^{{\mathscr{H}}}$, another new explicit second-order symplectic integrator for ${{\mathscr{H}}}^{* }$ reads

$$\begin{eqnarray}\begin{array}{rcl}{S}_{2}^{{{\mathscr{H}}}^{* }}(h) & = & {\tilde{{\mathscr{H}}}}_{4}^{* }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{3}^{* }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{2}^{* }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{1}^{* }(h)\\ & & \times \circ {\tilde{{\mathscr{H}}}}_{2}^{* }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{3}^{* }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{4}^{* }\left(\displaystyle \frac{h}{2}\right).\end{array}\end{eqnarray} \tag{ C7 }$$

A fourth-order method ${S}_{4}^{{{\mathscr{H}}}^{* }}$ can also be obtained. The new explicit symplectic algorithms in the new coordinate time w₁ are similar to those in the proper time of the Schwarzschild spacetime geometry in Paper I (Wang et al. 2021a).

If the time transformation function takes the form

$$\begin{eqnarray}&&{g}_{2}(r,\theta )={\rm{\Sigma }},\end{eqnarray} \tag{ C8 }$$

the new coordinate time in Equation (23) is w₂. In this case, Equation (24) corresponds to a simpler time-transformed Hamiltonian,

$$\begin{eqnarray}&&{{\mathscr{H}}}^{\star }={\rm{\Sigma }}(F+{p}_{0})+\displaystyle \frac{{\rm{\Delta }}}{2}{p}_{r}^{2}+\displaystyle \frac{1}{2}{p}_{\theta }^{2}.\end{eqnarray} \tag{ C9 }$$

It has three separable parts,

$$\begin{eqnarray}&&{{\mathscr{H}}}^{\star }={{\mathscr{H}}}_{1}^{\star }+{{\mathscr{H}}}_{2}^{\star }+{{\mathscr{H}}}_{3}^{\star },\end{eqnarray} \tag{ C10 }$$

where the sub-Hamiltonians are

$$\begin{eqnarray}&&{{\mathscr{H}}}_{1}^{\star }={\rm{\Sigma }}(F+{p}_{0}),\end{eqnarray} \tag{ C11 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{2}^{\star }=\displaystyle \frac{1}{2}{r}^{2}{p}_{r}^{2},\end{eqnarray} \tag{ C12 }$$

$$\begin{eqnarray}&&{{\mathscr{H}}}_{3}^{\star }=\displaystyle \frac{1}{2}({a}^{2}-2r){p}_{r}^{2}+\displaystyle \frac{1}{2}{p}_{\theta }^{2}.\end{eqnarray} \tag{ C13 }$$

Obviously, the three parts are independently solved in analytical methods, and their analytical solutions are explicit functions of the new coordinate time w₂. Using the Lie operators ${\tilde{{\mathscr{H}}}}_{1}^{\star }$, ${\tilde{{\mathscr{H}}}}_{2}^{\star }$, and ${\tilde{{\mathscr{H}}}}_{3}^{\star }$ to represent the analytical solutions of ${{\mathscr{H}}}_{1}^{\star }$, ${{\mathscr{H}}}_{2}^{\star }$, and ${{\mathscr{H}}}_{3}^{\star }$, we easily obtain an explicit second-order symplectic integrator for the system ${{\mathscr{H}}}^{\star }$ as follows:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{2}^{{{\mathscr{H}}}^{\star }}(h) & = & {\tilde{{\mathscr{H}}}}_{3}^{\star }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{2}^{\star }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{1}^{\star }(h)\\ & & \times \ \circ {\tilde{{\mathscr{H}}}}_{2}^{\star }\left(\displaystyle \frac{h}{2}\right)\circ {\tilde{{\mathscr{H}}}}_{3}^{\star }\left(\displaystyle \frac{h}{2}\right).\end{array}\end{eqnarray} \tag{ C14 }$$

An explicit fourth-order symplectic integrator ${S}_{4}^{{{\mathscr{H}}}^{\star }}$ can be easily available, too. As a point to illustrate, the newly fixed coordinate time step ${\rm{\Delta }}{w}_{2}=h$ corresponds to a variant proper time step ${\rm{\Delta }}{\tau }_{2}\approx {r}^{2}h$.

For $r\gg 1$, the fixed coordinate time steps ${\rm{\Delta }}{w}_{1}\ll 1$ and ${\rm{\Delta }}{w}_{2}\ll 1$ should be required when the variant proper time steps ${\rm{\Delta }}{\tau }_{1}$ and ${\rm{\Delta }}{\tau }_{2}$ are 1. Such extremely small coordinate time steps would lead to fatal round-off errors. Thus, we take the time transformation function Equation (25) rather than Equations (C1) and (C8).