Construction of Explicit Symplectic Integrators in General Relativity. I. Schwarzschild Black Holes

Black holes and gravitational waves were predicted in Einstein's theory of general relativity (Einstein 1915; Einstein & Sitzungsber 1916). The Schwarzschild solution was obtained from the field equations of a nonrotating black hole (Schwarzschild 1916). The Kerr solution was provided for a rotating black hole (Kerr 1963). The recent detection of gravitational waves (GW150914) from a binary black hole merger (Abbott et al. 2016) and the images of a supermassive black hole candidate at the center of the giant elliptical galaxy M87 (EHT Collaboration et al. 2019) provide powerful evidence for confirming the two predictions.

Although the relativistic equations of motion for test particles in the Schwarzschild and Kerr metrics are highly nonlinear, they are separable in variables and solved analytically in the Hamiltonian–Jacobi equation. Thus, they are integrable and the motions of particles near the two black holes are strictly regular. This integrability is attributed to the existence of four independent constants of motion, namely, energy, angular momentum, the four-velocity relation of particles, and the Carter constant (Carter 1968). However, no additional information regarding the solutions except the integrability of spacetimes is known because the solutions are expressed in terms of quadratures rather than elementary functions. Good numerical methods for computing these geodesics are highly desirable. In particular, when magnetic fields are included in curved spacetimes, the separation of variables in the Hamiltonian–Jacobi equation, associated with the equations of charged particle motion, is generally highly improbable. This condition may lead to the non-integrability of systems and the chaotic behavior of motion (Takahashi & Koyama 2009; Kopáček et al. 2010; Kopáček & Karas 2014; Kološ et al. 2015; Azreg-Aïnou 2016; Stuchlík & Kološ 2016; Tursunov et al. 2016; Li & Wu 2019). Numerical methods play an important role in analyzing the properties of these non-integrable problems.

Supposedly, good numerical methods are integrators that provide reliable results, particularly in the case of long-term integrations. In addition, the preservation of structural properties, such as symplectic structures, integrals of motion, phase-space volume, and symmetries, is desired. Such structure-preserving algorithms belong to a class of geometric integrators (Hairer et al. 1999). Among the properties, the most important ones are the preservation of energy and symplecticity.

In many cases, checking energy accuracy is a basic reference for testing the performance of numerical integration algorithms, although energy conservation does not necessarily yield high-precision numerical solutions. To demonstrate this scenario, we present a two-body problem as an example. Energy errors from the truncation or discretization errors of Runge–Kutta-type algorithms in the two-body problem typically increase linearly with integration time (Rein & Spiegel 2015). The growth speeds of in-track errors (Huang & Innanen 1983), which correspond to errors along the tangent to a trajectory in phase space, directly depend on the relative error in Keplerian energy (Avdyushev 2003). Accordingly, the Keplerian orbit is Lyapunov's instability, which leads to an increase in various errors. However, the stabilization or conservation of energy along the orbit is more efficient in eliminating Lyapunov's instability and the fast drifting of in-track errors than that of other integrals. The energy stabilization method of Baumgarte (1972, 1973) includes known integrals (such as an energy integral) in the equations of motion. The stabilization in the perturbed two-body and restricted three-body problems of satellites, asteroids, stars, and planets has been demonstrated to improve the accuracy of numerical integrations by several orders of magnitude (Avdyushev 2003). In contrast with Baumgarte's method, the manifold correction or projection method of Nacozy (1971) applies a least-squares procedure to add a linear correction vector to a numerical solution. This vector is computed from the gradient vectors of the integrals. The application of Nacozy's method is generalized to quasi-Keplerian motions of perturbed two-body or N-body problems with the aid of the integral invariant relation of slowly varying individual Kepler energies (Wu et al. 2007; Ma et al. 2008a). Some projection methods (Fukushima 2003a, 2003b, 2003c, 2004; Ma et al. 2008b; Wang et al. 2016, 2018; Deng et al. 2020) for rigorously satisfying integrals, including Kepler energy in a two-body problem, have been proposed and extended to perturbed two-body problems, N-body systems, non-conservative elliptic restricted three-body problems, and dissipative circular restricted three-body problems. In addition to explicit projection methods that exactly preserve the energy integral, exact energy-preserving implicit integration methods that discretize Hamiltonian gradients in terms of the average Hamiltonian difference terms have been specifically designed for conservative Hamiltonian systems (Feng & Qin 2009; Bacchini et al. 2018a, 2018b; Hu et al. 2019).

Although energy-preserving integrators and some projection methods exactly conserve energy, they are non-symplectic. Symplectic algorithms (Wisdom 1982; Ruth 1983; Feng 1986; Suzuki 1991; McLachlan & Atela 1992; Chin 1997; Omelyan et al. 2002a, 2002b, 2003) do not exactly conserve the energy of a Hamiltonian system, but they cause energy errors to oscillate and become bounded as evolution time increases. In this manner, these algorithms are also considered to conserve energy efficiently over long-term integrations. Moreover, they preserve the symplectic structure of Hamiltonian flows. Given these two advantages, symplectic integrators are widely used in long-term studies on solar system dynamics. The most popular algorithms in solar system dynamics are the second-order symplectic integrator of Wisdom & Holman (1991) and its extensions (Wisdom et al. 1996; Chambers & Murison 2000; Laskar & Robutel 2001; Hernandez & Dehnen 2017). Notably, the explicit symplectic algorithms in a series of references (Suzuki 1991; Chin 1997; Omelyan et al. 2002a, 2002b, 2003) require the integrated Hamiltonian to be split into two parts with analytical solutions as explicit functions of time. However, the two splitting parts from the Hamiltonian in Wisdom & Holman (1991), Wisdom et al. (1996), Chambers & Murison (2000), and Laskar & Robutel (2001) should be the primary and secondary parts. For the latter, the analytical solutions can be given in explicit functions of time. The former also has explicit analytical solutions, but eccentric anomaly is calculated using an iteration method, such as the Newton–Raphson method.

However, a relativistic gravitational Hamiltonian system, such as the Schwarzschild spacetime, is inseparable or has no two separable parts with analytical solutions being explicit functions of proper time. This condition leads to the difficulty in applying explicit symplectic integrators. By extending the phase space of such an inseparable Hamiltonian system, Pihajoki (2015) obtained a new Hamiltonian consisting of two sub-Hamiltonians equal to the original Hamiltonian, where one sub-Hamiltonian is a function of the original coordinates and new momenta, and the other is a function of the original momenta and new coordinates. The two sub-Hamiltonians are separable in variables; therefore, standard explicit symplectic leapfrog splitting methods are applicable to the new Hamiltonian. Mixing maps of feedback between the two sub-Hamiltonian solutions and a map for projecting a vector in the extended phase space back to the original number of dimensions are necessary and require a suitable choice. Liu et al. (2016) confirmed that sequent permutations of coordinates and momenta achieve good results in preserving the original Hamiltonian without an increase in secular errors compared with the permutations of momenta suggested by Pihajoki (2015). Luo et al. (2017) found that midpoint permutations exhibit the best results. However, mixing maps generally destroy symplecticity in extended phase space. In addition, extended phase space leapfrogs are not symplectic for the use of any projection map. Despite the absence of symplecticity, mixing and projection maps are used only as output and exert no influence on the state in extended phase space. Consequently, leapfrogs, such as partitioned multistep methods, can exhibit good long-term behavior in stabilizing the original Hamiltonian (Liu et al. 2017; Luo & Wu 2017; Wu & Wu 2018). Thus, extended phase-space leapfrog methods, including extended phase-space logarithmic Hamiltonian methods (Li & Wu 2017), are called explicit symplectic-like integrators. In addition to the two copies of the original system with mixed-up positions and momenta, a third sub-Hamiltonian, as an artificial restraint to the divergence between the original and extended variables, was introduced by Tao (2016). Neither mixing nor projection maps are used in Tao's method, and thus explicit leapfrog methods are still symplectic in the extended phase space. Two problems exist. (i) A binding constant for controlling divergence has an optimal choice. This choice cannot be given theoretically but requires considerable values to test which one minimizes the original Hamiltonian error. (ii) Whether the original variables in the newly extended Hamiltonian coincide with those in the original Hamiltonian is unclear.

To date, no standard explicit symplectic leapfrogs, only implicit symplectic methods, have been established in a relativistic Hamiltonian problem because of the difficulty in separating variables. The second-order implicit midpoint method (Feng 1986) is the most common choice among implicit symplectic methods. It can function as a variational symplectic integrator for constrained Hamiltonian systems (Brown 2006). To save computational cost, explicit and implicit combined symplectic algorithms have been provided by some authors (Liao 1997; Preto & Saha 2009; Lubich et al. 2010; Zhong et al. 2010; Mei et al. 2013a, 2013b). Notably, the symplectic integration scheme for the post-Newtonian motion of a spinning black hole binary (Lubich et al. 2010) is non-canonical because of the use of non-canonical spin variables. However, this scheme can become canonical when canonically conjugated cylindrical-like spin coordinates (Wu & Xie 2010) are used. The symplectic implicit Gauss–Legendre Runge–Kutta method has been applied to determine the regular and chaotic behavior of charged particles around a Kerr black hole immersed in a weak, asymptotically uniform magnetic field (Kopáček et al. 2010). Implicit symmetric schemes with adaptive step size control that effectively conserve the integrals of motion are appropriate for studying geodesic orbits in curved spacetime backgrounds (Seyrich & Lukes-Gerakopoulos 2012). Symplectic integrators for general non-conservative systems (Tsang et al. 2015) can share many benefits of traditional symplectic integrators.

In general, implicit symplectic methods are numerically more expensive to solve than same-order explicit symplectic integrators. The latter algorithms should be used if possible. Accordingly, we intend to address the difficulty in constructing explicit symplectic integrators for Schwarzschild-type spacetimes similar to the standard explicit symplectic leapfrogs for Hamiltonian problems in solar system dynamics. If the Hamiltonians of Schwarzschild-type spacetimes are separated into two parts that resemble the splitting form of Hamiltonian systems in the construction of standard symplectic leapfrogs, then no explicit symplectic algorithms are available. The conditions for constructing explicit symplectic schemes may require Hamiltonians to be split into more parts with analytical solutions as explicit functions of proper time.

The remainder of this paper is organized as follows. In Section 2, we briefly introduce the standard explicit symplectic leapfrog and its extensions for a separable Hamiltonian system. The Hamiltonian of charged particles moving around a Schwarzschild black hole with an external magnetic field is described in Section 3. Explicit symplectic schemes are designed for curved Schwarzschild spacetimes in Section 4. The performance of explicit symplectic integrators is tested numerically in Section 5. Section 6 gives the major results. A discrete difference scheme of the new second-order explicit symplectic integrator is presented in Appendix A. Explicit and implicit combined symplectic methods and extended phase-space explicit symplectic-like methods are provided in Appendix B.

Set q as an N-dimensional coordinate vector. Its corresponding generalized momentum is p. Let Z = (p, q) be a 2N-dimensional phase-space variable. Consider the following Hamiltonian:

$$\begin{eqnarray}&&H({\boldsymbol{p}},{\boldsymbol{q}})={H}_{1}({\boldsymbol{p}},{\boldsymbol{q}})+{H}_{2}({\boldsymbol{p}},{\boldsymbol{q}}),\end{eqnarray} \tag{ 1 }$$

where the two separable parts H₁ and H₂ are supposed to be independently integrable. A typical splitting form of H takes H₁ as kinetic energy T(p) and H₂ as potential V(q).

Two differential operators are defined as follows:

$$\begin{eqnarray*}\begin{array}{rcl}{ \mathcal A } & = & \displaystyle \sum _{i=1}^{N}\left(\displaystyle \frac{\partial {H}_{1}}{\partial {{\boldsymbol{p}}}_{i}}\displaystyle \frac{\partial }{\partial {{\boldsymbol{q}}}_{i}}-\displaystyle \frac{\partial {H}_{1}}{\partial {{\boldsymbol{q}}}_{i}}\displaystyle \frac{\partial }{\partial {{\boldsymbol{p}}}_{i}}\right),\\ { \mathcal B } & = & \displaystyle \sum _{i=1}^{N}\left(\displaystyle \frac{\partial {H}_{2}}{\partial {{\boldsymbol{p}}}_{i}}\displaystyle \frac{\partial }{\partial {{\boldsymbol{q}}}_{i}}-\displaystyle \frac{\partial {H}_{2}}{\partial {{\boldsymbol{q}}}_{i}}\displaystyle \frac{\partial }{\partial {{\boldsymbol{p}}}_{i}}\right).\end{array}\end{eqnarray*}$$

System (1) has the following formal solution:

$$\begin{eqnarray}&&{\boldsymbol{Z}}(h)={ \mathcal C }(h){\boldsymbol{Z}}(0),\end{eqnarray} \tag{ 2 }$$

where Z(0) denotes the value of Z in the beginning of time step h. The differential operator ${ \mathcal C }={ \mathcal A }+{ \mathcal B }$ is approximately expressed as a series of products of ${ \mathcal A }$ and ${ \mathcal B }$:

$$\begin{eqnarray}&&{ \mathcal C }(h)\approx {{\rm{\Pi }}}_{j=1}^{e}{ \mathcal A }(h{\alpha }_{j}){ \mathcal B }(h{\beta }_{j})+O({h}^{d+1}),\end{eqnarray} \tag{ 3 }$$

where coefficients α_j and β_j are determined by the conditions of order d. In this manner, symplectic numerical integrators of arbitrary orders are built.

If d = 2, then Equation (3) is the Verlet algorithm (Swope et al. 1982)

$$\begin{eqnarray}&&{{ \mathcal S }}_{2}(h)={ \mathcal A }\left(\displaystyle \frac{h}{2}\right){ \mathcal B }(h){ \mathcal A }\left(\displaystyle \frac{h}{2}\right).\end{eqnarray} \tag{ 4 }$$

This algorithm is an explicit standard symplectic leapfrog method. When d = 4, Equation (3) corresponds to the explicit symplectic algorithm of Forest & Ruth (1990):

$$\begin{eqnarray}\begin{array}{rcl}{FR}4(h) & = & { \mathcal A }\left(\displaystyle \frac{\gamma }{2}h\right){ \mathcal B }(\gamma h){ \mathcal A }\left(\displaystyle \frac{1-\gamma }{2}h\right){ \mathcal B }((1-2\gamma )h)\\ & & \circ \,{ \mathcal A }\left(\displaystyle \frac{1-\gamma }{2}h\right){ \mathcal B }(\gamma h){ \mathcal A }\left(\displaystyle \frac{\gamma }{2}h\right),\end{array}\end{eqnarray} \tag{ 5 }$$

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

Evidently, the construction of these explicit symplectic integrators is based on the Hamiltonian with an analytically integrable decomposition. Can such an operator-splitting technique be available in strictly general relativistic systems, such as a Schwarzschild spacetime? The succeeding discussions answer this question.

A Schwarzschild black hole with mass M is a non-rotating black hole. In spherical-like coordinates (t, r, θ, ϕ), the Schwarzschild metric is described by

$$\begin{eqnarray}\begin{array}{rcl}-{c}^{2}d{\tau }^{2} & = & {{ds}}^{2}={g}_{\alpha \beta }{{dx}}^{\alpha }{{dx}}^{\beta }\\ & = & -\left(1-\displaystyle \frac{2{GM}}{{{rc}}^{2}}\right){c}^{2}{{dt}}^{2}+{\left(1-\displaystyle \frac{2{GM}}{{{rc}}^{2}}\right)}^{-1}\\ & & \cdot \,{{dr}}^{2}+{r}^{2}d{\theta }^{2}+{r}^{2}{\sin }^{2}\theta d{\phi }^{2},\end{array}\end{eqnarray} \tag{ 6 }$$

where τ, c, and G denote proper time, the speed of light, and constant of gravity, respectively. In general, c and G use geometrized units, c = G = 1. M also has one unit, M = 1. This unit mass can be obtained via scale transformations to certain quantities: t → tM, r → rM, and τ → τM. In this manner, this metric is transformed into a dimensionless form as follows:

$$\begin{eqnarray}\begin{array}{rcl}-d{\tau }^{2}={{ds}}^{2} & = & -\left(1-\displaystyle \frac{2}{r}\right){{dt}}^{2}+{\left(1-\displaystyle \frac{2}{r}\right)}^{-1}{{dr}}^{2}\\ & & +{r}^{2}d{\theta }^{2}+{r}^{2}{\sin }^{2}\theta d{\phi }^{2}.\end{array}\end{eqnarray} \tag{ 7 }$$

This metric corresponds to a Lagrangian system

$$\begin{eqnarray}&&{ \mathcal L }=\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 }$$

where ${\dot{x}}^{\mu }={\boldsymbol{U}}$ is a four-velocity. A covariant generalized momentum p is defined in the following form:

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

This Lagrangian does not explicitly depend on t and ϕ, and thus two constant momentum components exist. They are

$$\begin{eqnarray}&&{p}_{t}=-\left(1-\displaystyle \frac{2}{r}\right)\dot{t}=-E,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{p}_{\phi }={r}^{2}{\sin }^{2}\theta \dot{\phi }={\ell },\end{eqnarray} \tag{ 11 }$$

where E and ℓ are the energy and angular momentum of a test particle moving around a black hole, respectively.

In accordance with classical mechanics, a Hamiltonian derived from the Lagrangian is expressed as

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal H } & = & {\boldsymbol{U}}\cdot {\boldsymbol{p}}-{ \mathcal L }=\displaystyle \frac{1}{2}{g}^{\mu \nu }{p}_{\mu }{p}_{\nu }=-\displaystyle \frac{1}{2}{\left(1-\displaystyle \frac{2}{r}\right)}^{-1}{E}^{2}\\ & & +\,\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{2}{r}\right){p}_{r}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{{p}_{\theta }^{2}}{{r}^{2}}+\displaystyle \frac{1}{2}\displaystyle \frac{{{\ell }}^{2}}{{r}^{2}{\sin }^{2}\theta }.\end{array}\end{eqnarray} \tag{ 12 }$$

This Hamiltonian governs the motion of a test particle around the Schwarzschild black hole.

A point is worth noting. A magnetic field arises due to the relativistic motion of charged particles in an accretion disk around the central black hole (Borm & Spaans 2013). It also leads to generating gigantic jets along the magnetic axes. The magnetic field is too weak to change the gravitational background and alter the metric tensor of the Schwarzschild black hole spacetime. However, it can exert a considerable influence on the motion of charged test particles. Considering this point, we suppose that the particle has a charge q and the black hole is immersed into an external asymptotically uniform magnetic field. The magnetic field is parallel to the z-axis, and its strength is B. The electromagnetic four-vector potential A^α in the Lorentz gauge is a linear combination of the time-like and space-like axial Killing vectors ${\xi }_{(t)}^{\alpha }$ and ${\xi }_{(\phi )}^{\alpha }$ (Abdujabbarov et al. 2013; Shaymatov et al. 2015; Tursunov et al. 2016; Benavides-Gallego et al. 2019):

$$\begin{eqnarray}&&{A}^{\alpha }={C}_{1}{\xi }_{(t)}^{\alpha }+{C}_{2}{\xi }_{(\phi )}^{\alpha }.\end{eqnarray} \tag{ 13 }$$

In Felice & Sorge (2003), the constants are set as C₁ = 0 and C₂ = B/2. In this manner, the four-vector potential has only one non-zero covariant component

$$\begin{eqnarray}&&{A}_{\phi }=\displaystyle \frac{B}{2}{g}_{\phi \phi }=\displaystyle \frac{B}{2}{r}^{2}{\sin }^{2}\theta .\end{eqnarray} \tag{ 14 }$$

The charged particle motion is described by the Hamiltonian system

$$\begin{eqnarray}\begin{array}{rcl}K & = & \displaystyle \frac{1}{2}{g}^{\mu \nu }({p}_{\mu }-{{qA}}_{\mu })({p}_{\nu }-{{qA}}_{\nu })\\ & = & -\displaystyle \frac{1}{2}{\left(1-\displaystyle \frac{2}{r}\right)}^{-1}{E}^{2}+\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{2}{r}\right){p}_{r}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{{p}_{\theta }^{2}}{{r}^{2}}\\ & & +\,\displaystyle \frac{1}{2{r}^{2}{\sin }^{2}\theta }{\left(L-\displaystyle \frac{\beta }{2}{r}^{2}{\sin }^{2}\theta \right)}^{2},\end{array}\end{eqnarray} \tag{ 15 }$$

where β = qB. The energy E is still determined using Equation (10). However, the expression of angular momentum is dissimilar to that of Equation (11) and is presented as

$$\begin{eqnarray}&&L={r}^{2}{\sin }^{2}\theta \dot{\phi }+\displaystyle \frac{\beta }{2}{r}^{2}{\sin }^{2}\theta .\end{eqnarray} \tag{ 16 }$$

A point is illustrated here. The dimensionless Hamiltonian (15) is obtained after scale transformations of B → B/M, E → mE, p_r → mp_r, q → mq, L → mML, p_θ → mMp_θ, and K → m²K, where m is the particle's mass. In addition, the Schwarzschild solution with an external magnetic field is the Hamiltonian (15), and it no longer has a background solution to general relativity.

The Hamiltonians ${ \mathcal H }$ and K always remain at a given constant as follows:

$$\begin{eqnarray}&&{ \mathcal H }=-\displaystyle \frac{1}{2},\end{eqnarray} \tag{ 17 }$$

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

They are attributed to the four-velocity relation U · U = −1. In addition, a second integral (i.e., the Carter constant) can be easily found in the Hamiltonian ${ \mathcal H }$ by performing the separation of variables in the Hamilton–Jacobi equation. Thus, this Hamiltonian is integrable and has formal analytical solutions. However, the perturbation from the external magnetic field leads to the absence of a second integral. In this case, no formal analytical solutions exist in the Hamiltonian K.

Suppose the Hamiltonian (12) is similar to the Hamiltonian (1) and has two splitting parts:

$$\begin{eqnarray}&&{ \mathcal H }={ \mathcal T }+{ \mathcal V },\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{ \mathcal T }=\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{2}{r}\right){p}_{r}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{{p}_{\theta }^{2}}{{r}^{2}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{ \mathcal V }=-\displaystyle \frac{1}{2}{\left(1-\displaystyle \frac{2}{r}\right)}^{-1}{E}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{{{\ell }}^{2}}{{r}^{2}{\sin }^{2}\theta }.\end{eqnarray} \tag{ 21 }$$

The ${ \mathcal V }$ part is analytically integrable, and its analytical solutions p_r and p_θ are explicit functions of proper time τ. Although the ${ \mathcal T }$ part exhibits no separation of variables, it is still analytically integrable. However, its analytical solutions r and p_r are not explicit functions of proper time τ but are implicit functions. In this case, the explicit symplectic integrators in Equations (4) and (5) are unsuitable for the Hamiltonian-splitting form (19). Consequently, implicit symplectic integrators rather than explicit ones can be constructed in relativistic Hamiltonian systems, such as Equation (12), in the general case. The ${ \mathcal V }$ part is more complicated and is not a separation of variables in most cases in general relativity. Thus, the construction of explicit symplectic methods becomes more difficult.

From the preceding demonstrations, the key for constructing explicit symplectic integrators requires the integrated Hamiltonian to exist as an analytically integrable decomposition. In particular, the obtained analytical solutions for each splitting part should be explicit functions of proper time τ. In summary, the two points must be satisfied for constructing explicit symplectic integrators. The Hamiltonian (12) with the two analytically integrable splitting parts fails to construct any explicit symplectic scheme. Subsequently, we focus on the Hamiltonian with more analytically integrable splitting parts.

We split the Hamiltonian ${ \mathcal H }$ into four pieces:

$$\begin{eqnarray}&&{ \mathcal H }={{ \mathcal H }}_{1}+{{ \mathcal H }}_{2}+{{ \mathcal H }}_{3}+{{ \mathcal H }}_{4},\end{eqnarray} \tag{ 22 }$$

where these sub-Hamiltonians are

$$\begin{eqnarray}&&{{ \mathcal H }}_{1}=\displaystyle \frac{1}{2}\displaystyle \frac{{{\ell }}^{2}}{{r}^{2}{\sin }^{2}\theta }-\displaystyle \frac{1}{2}{\left(1-\displaystyle \frac{2}{r}\right)}^{-1}{E}^{2},\end{eqnarray} \tag{ 23 }$$

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

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

$$\begin{eqnarray}&&{{ \mathcal H }}_{4}=\displaystyle \frac{{p}_{\theta }^{2}}{2{r}^{2}}.\end{eqnarray} \tag{ 26 }$$

For the sub-Hamiltonian ${{ \mathcal H }}_{1}$, its canonical equations are $\dot{r}=\dot{\theta }=0$ and

$$\begin{eqnarray}&&\displaystyle \frac{{{dp}}_{r}}{d\tau }=-\displaystyle \frac{\partial {{ \mathcal H }}_{1}}{\partial r}=\displaystyle \frac{{{\ell }}^{2}}{{r}^{3}{\sin }^{2}\theta }-\displaystyle \frac{{E}^{2}}{{\left(r-2\right)}^{2}},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dp}}_{\theta }}{d\tau }=-\displaystyle \frac{\partial {{ \mathcal H }}_{1}}{\partial \theta }=\displaystyle \frac{{{\ell }}^{2}\cos \theta }{{r}^{2}{\sin }^{3}\theta }.\end{eqnarray} \tag{ 28 }$$

Evidently, r and θ are constants when proper time goes from τ₀ to τ₁ = τ₀ + τ. Thus, p_r and p_θ can be solved analytically from Equations (27) and (28) They are explicit functions of τ in the following form

$$\begin{eqnarray}&&{p}_{r}(\tau )={p}_{r0}+\tau \left[\displaystyle \frac{{{\ell }}^{2}}{{r}_{0}^{3}{\sin }^{2}{\theta }_{0}}-\displaystyle \frac{{E}^{2}}{{\left({r}_{0}-2\right)}^{2}}\right],\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{p}_{\theta }(\tau )={p}_{\theta 0}+\tau \displaystyle \frac{{{\ell }}^{2}\cos {\theta }_{0}}{{r}_{0}^{2}{\sin }^{3}{\theta }_{0}},\end{eqnarray} \tag{ 30 }$$

where r₀, θ₀, p_r0, and p_θ0 represent values of r, θ, p_r and p_θ at the proper time τ₀, and p_r(τ) and p_θ(τ) denote the values of p_r and p_θ at proper time τ₁. A differential operator for solving ${{ \mathcal H }}_{1}$ is labeled as ${\psi }_{\tau }^{{{ \mathcal H }}_{1}}$.

The canonical equations of the sub-Hamiltonians ${{ \mathcal H }}_{2}$, ${{ \mathcal H }}_{3}$, and ${{ \mathcal H }}_{4}$ are

$$\begin{eqnarray}&&{{ \mathcal H }}_{2}:\,\displaystyle \frac{{dr}}{d\tau }={p}_{r},\,\,{\dot{p}}_{r}=0;\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{3}:\,\displaystyle \frac{{dr}}{d\tau }=-\displaystyle \frac{2}{r}{p}_{r},\,\displaystyle \frac{{{dp}}_{r}}{d\tau }=-\displaystyle \frac{{p}_{r}^{2}}{{r}^{2}};\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{4}:\,\displaystyle \frac{d\theta }{d\tau }=\displaystyle \frac{{p}_{\theta }}{{r}^{2}},\,\displaystyle \frac{{{dp}}_{r}}{d\tau }=\displaystyle \frac{{p}_{\theta }^{2}}{{r}^{3}},\,\dot{r}={\dot{p}}_{\theta }=0.\end{eqnarray} \tag{ 33 }$$

Let ${\psi }_{\tau }^{{{ \mathcal H }}_{2}}$, ${\psi }_{\tau }^{{{ \mathcal H }}_{3}}$, and ${\psi }_{\tau }^{{{ \mathcal H }}_{4}}$ be three operators. We obtain the solutions for Equations (31)–(33)as follows:

$$\begin{eqnarray}&&{\psi }_{\tau }^{{{ \mathcal H }}_{2}}:\,r(\tau )={r}_{0}+\tau {p}_{r0};\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\psi }_{\tau }^{{{ \mathcal H }}_{3}}:\,r(\tau ) & = & {\left[{\left({r}_{0}^{2}-3\tau {p}_{r0}\right)}^{2}/{r}_{0}\right]}^{1/3},\\ {p}_{r}(\tau ) & = & {p}_{r0}{\left[({r}_{0}^{2}-3\tau {p}_{r0})/{r}_{0}^{2}\right]}^{1/3};\end{array}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\psi }_{\tau }^{{{ \mathcal H }}_{4}}:\,\theta (\tau ) & = & {\theta }_{0}+\tau {p}_{\theta 0}/{r}_{0}^{2},\\ {p}_{r}(\tau ) & = & {p}_{r0}+\tau {p}_{\theta 0}^{2}/{r}_{0}^{3}.\end{array}\end{eqnarray} \tag{ 36 }$$

It is clear that these solutions are explicit functions of proper time τ. If the sum of ${{ \mathcal H }}_{2}$ and ${{ \mathcal H }}_{3}$ is regarded as an independent sub-Hamiltonian, then it is analytically solved. However, the analytical solutions of r, θ, and p_r for the sum cannot be expressed as explicit functions of proper time τ. Thus, such a composed sub-Hamiltonian is not considered. Equation (22) is a possible Hamiltonian splitting for satisfying this requirement. Other appropriate splitting forms may be provided to the Hamiltonian (12).

The flow ${\psi }_{h}^{{ \mathcal H }}$ of the Hamiltonian (12) over time step h is approximately given by the symmetric composition of these operators:

$$\begin{eqnarray}\begin{array}{rcl}{\psi }_{h}^{{ \mathcal H }}\approx {S}_{2}^{{ \mathcal H }}(h) & = & {\psi }_{h/2}^{{{ \mathcal H }}_{4}}\circ {\psi }_{h/2}^{{{ \mathcal H }}_{3}}\circ {\psi }_{h/2}^{{{ \mathcal H }}_{2}}\circ {\psi }_{h}^{{{ \mathcal H }}_{1}}\\ & & \circ \,{\psi }_{h/2}^{{{ \mathcal H }}_{2}}\circ {\psi }_{h/2}^{{{ \mathcal H }}_{3}}\circ {\psi }_{h/2}^{{{ \mathcal H }}_{4}}.\end{array}\end{eqnarray} \tag{ 37 }$$

The above construction is a second-order explicit symplectic integrator marked as ${S}_{2}^{{ \mathcal H }}$. Its difference scheme is provided in Appendix A.

The order of algorithm (37) can be lifted to four by using the composition scheme of Yoshida (1990). That is, a fourth-order symplectic composition construction is

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

where δ = 1–2γ.

The Hamiltonian (15) exhibits the following splitting form:

$$\begin{eqnarray}&&K={K}_{1}+{K}_{2}+{K}_{3}+{K}_{4},\end{eqnarray} \tag{ 39 }$$

where ${K}_{2}={{ \mathcal H }}_{2}$, ${K}_{3}={{ \mathcal H }}_{3}$, ${K}_{4}={{ \mathcal H }}_{4}$, and the inclusion of A_ϕ only changes ${{ \mathcal H }}_{1}$ as

$$\begin{eqnarray}\begin{array}{rcl}{K}_{1} & = & \displaystyle \frac{1}{2{r}^{2}{\sin }^{2}\theta }{\left(L-\displaystyle \frac{\beta }{2}{r}^{2}{\sin }^{2}\theta \right)}^{2}\\ & & -\displaystyle \frac{1}{2}{\left(1-\displaystyle \frac{2}{r}\right)}^{-1}{E}^{2}.\end{array}\end{eqnarray} \tag{ 40 }$$

When ${{ \mathcal H }}_{1}$ gives place to K₁, the explicit symplectic integrators S₂ and S₄ are still suitable for the non-integrable Hamiltonian K of the Schwarzschild solution with an external magnetic field, labeled as S^K₂ and S^K₄.

In summary, when the Hamiltonians (12) and (15) are split into four analytically integrable parts, their explicit symplectic integrators are easily constructed.

In this section, we focus on checking the numerical performance of the proposed integrators. For comparison, a conventional fourth-order Runge–Kutta integrator (RK4), second- and fourth-order symplectic algorithms consisting of explicit and implicit mixed methods (EI2 and EI4), and second- and fourth-order extended phase-space explicit symplectic-like methods (EE2 and EE4) are used. The details of EI2, EI4, EE2, and EE4 are provided in Appendix B.

5.1. Case of β = 0

When no charges are assigned to test particles, the system (15) is transformed to the Schwarzschild problem (12). We consider parameters E = 0.995 and ℓ (or L) = 4.6, and proper time step size h = 1. Initial conditions are r = 11, θ = π/2 and p_r = 0. The initial value of p_θ (>0) is determined using Equation (17). We conduct our numerical experiments by applying each of the aforementioned algorithms to solve the Hamiltonian (12). As shown in Figure 1(a), the three second-order methods, namely, S2, EI2, and EE2, provide an order of 10⁻⁶ to Hamiltonian errors ${\rm{\Delta }}H=1+2{ \mathcal H }$ from Equation (17) at the end of integration time. Differences exist among the algorithmic errors. The new symplectic algorithm S2 and the explicit and implicit mixed symplectic method EI2 have nearly the same errors, which remain bounded and stable. This result indicates the superiority of S2 in the conservation of the long-term stable behavior of energy (or Hamiltonian) errors. However, the extended phase-space method EE2 exhibits an increase in secular errors. This increase can be prevented if a small time size h = 0.1 is used. In such case, the errors (not plotted) can be stabilized within an order of 10⁻⁸.

Zoom In Zoom Out Reset image size

The four fourth-order algorithms, namely, S4, EI4, EE4, and RK4, yield the Hamiltonian errors in Figures 1(b) and (c). The algorithms S4, EI4, and EE4 are accurate to an order of 10⁻⁸. The new method S4 and the extended phase-space method EE4 have stable and bounded errors. The explicit and implicit mixed symplectic method EI4 causes the errors to become bounded. Meanwhile, RK4 provides the lowest accuracy with an order of 10⁻⁶ and its errors increase linearly with time. This result is expected because RK4 is not a geometric integrator.

The considered orbit, called Orbit 1, can be observed from the Poincaré section map on the plane θ = π/2 and p_θ > 0. The map relates to a two-dimensional plane, which exhibits intersections of the particles' trajectories with the surface of section in phase space (Lichtenberg & Lieberman 1983). If the plotted points form a closed curve, then the motion is regular. This result is based on a regular trajectory moving on a torus in phase space and the curve being a cross-section of the torus. By contrast, if the plotted points are distributed randomly, then the motion is chaotic. With the aid of the distribution of the points in the Poincaré map, we can determine the phase-space structure, indicating whether the motion is chaotic. The Kolmogorov–Arnold–Moser (KAM) torus in the section in Figure 1(d) is provided by the new method S2 and indicates the regularity of Orbit 1. In addition, the structure of Orbit 1, and those of Orbits 2 and 3 with initial separations r = 70 and 110 are described. The numerical performance of the aforementioned algorithms acting on Orbit 1 is approximately consistent with those acting on Orbits 2 and 3.

5.2. Case of β ≠ 0

When an external magnetic field with parameter β = 8.9 × 10⁻⁴ is included within the vicinity of a black hole, the system is non-integrable. The magnetic field causes the three orbits in Figure 1(d) to have different phase-space structures in Figure 2(a). Although Orbit 1 remains a simply closed torus, it is shrunk drastically and becomes a small torus. By contrast, Orbit 2 becomes a more complicated KAM torus, consisting of seven small loops wherein the successive points jump from one loop to the next. These small loops belong to the same trajectory and form a chain of islands (Hénon & Heiles 1964). Such a torus is regular but easily induces the occurrence of resonance and chaos. In particular, Orbit 3, which is a small loop in Figure 1(d), is considerably enlarged and densely filled in the phase space. This result indicates the onset of strong chaoticity.

Zoom In Zoom Out Reset image size

Although the loop of Orbit 1 is considerably smaller under the interaction of the electromagnetic forces in Figure 2(a) than in the case without electromagnetic forces in Figure 1(d), each algorithm exhibits nearly the same performance in the two cases because the tori of Orbit 1 in the two cases belong to the same category of trajectories, namely, simple single regular loops. Orbits 2 and 3 exhibit completely different dynamical behavior, but correspond to approximately the same Hamiltonian errors for each integration method. Figures 2(b)–(d) plot the errors for chaotic Orbit 3. The errors of the second-order methods for chaotic Orbit 3 shown in Figure 2(b) are approximately consistent with those for regular Orbit 1 shown in Figure 1(a). The fourth-order algorithms S4 and EE4 exhibit no dramatic differences in errors in Figure 2(c), similar to that in Figure 1(b). This result indicates that orbital chaoticity does not explicitly affect algorithmic accuracy. However, the explicit and implicit mixed method EI4 presents a secular drift in errors due to roundoff errors. The increase in errors can be prevented when a large time size h = 10 is adopted. In this case, accuracy is maintained with an order of 10⁻⁵. EI4 exhibits secular drift in the Hamiltonian errors for the smaller time step h = 1 but does not for the larger time size h = 10. The following is a simple analysis. The errors of a symplectic integrator mostly consist of truncation and roundoff errors. When truncation errors are more than roundoff errors, the symplectic integrator causes the Hamiltonian errors to remain bounded and to exhibit no secular drift in appropriate situations. Roundoff errors increase with an increase in the number N of calculations. They are approximately estimated using N, where  ∼ 10⁻¹⁶ demonstrates machine precision in double floating-point precision. When roundoff errors completely dominate total errors, the Hamiltonian or energy errors increase linearly with time. Assume that a symplectic method has a truncation energy error on an order of 10⁻¹². The total errors in the energy are stabilized at the order of magnitude when N < 10⁴, but grow linearly as N ≫ 10⁴. If a symplectic method has a truncation energy error higher than the order of 10⁻⁸, then the total errors in the energy remain bounded and approach the order of truncation errors when N < 10⁸, whereas they increase linearly as N ≫ 10⁸. These results have been confirmed by numerical experiments on N-body problems in the solar system (Wu et al. 2003; Deng et al. 2020). In the present numerical simulations, the truncation Hamiltonian errors of EI4 are on the order of 10⁻⁹ for h = 1 but the roundoff errors are 10⁻⁸ after 10⁸ integration steps. Given that the former errors are smaller than the latter ones, secular drift exists in the Hamiltonian errors. However, the truncation Hamiltonian errors of EI4 are on the order of 10⁻⁵ for h = 10. They are larger than the roundoff errors after 10⁸ integration steps. Therefore, no secular drift occurs in the Hamiltonian errors.

A conclusion can be drawn from Figures 1 and 2 that the stable behavior and magnitude of the Hamiltonian errors for each algorithm mostly depend on the choice of step sizes. To demonstrate this fact clearly, we list them in Tables 1 and 2, where chaotic Orbit 3 is used as a test orbit. The two second-order symplectic integrators S2 and EI2 can make the errors bounded for the three time steps, h = 0.1, 1, 10. A larger time step is also suitable for the two fourth-order symplectic integrators S4 and EI4. However, a smaller time step is suitable for the extended phase-space methods. The reason why EE2 does not produce stable errors for h = 1 but does for h = 0.1 (or EE4 does not produce stable errors for h = 10 but does for h = 1) differs from why S4 does not provide stable errors for h = 0.1 but does for h = 1. The error stability or instability for the former case is mostly dependent on permutations, which are frequently required in appropriately small times. However, it is primarily related to the roundoff errors for the latter case. Such a smaller time step is also necessary for RK4 to obtain higher accuracy, although RK4 does not remain at a stable or bounded value of energy errors.

Table 1.  Dependence of Stable (S) or Unstable (U) Behavior of Hamiltonian Errors for the Seven Algorithms on Step Size h

Method	S2	EI2	EE2	S4	EI4	EE4	RK4
h = 0.1	S	S	S	U	U	S	U
h = 1.0	S	S	U	S	U	S	U
h = 10	S	S	U	S	S	U	U

Note. Chaotic Orbit 3 in Figure 2 is integrated until proper time τ = 10⁸.

Download table as:  ASCIITypeset image

Table 2.  Same as Table 1, but Dependence of the Largest Absolute Values of Hamiltonian Errors on h

Method	S2	EI2	EE2	S4	EI4	EE4	RK4
h = 0.1	4e-8	4e-8	3e-8	7e-9	3e-12	1e-12	4e-12
h = 1.0	6e-6	5e-6	2e-6	3e-8	7e-9	2e-8	4e-7
h = 10	8e-4	6e-3	6e-3	4e-4	7e-5	4e-3	3e-2

Download table as:  ASCIITypeset image

Computational costs are listed in Table 3. Given the smaller step sizes, several differences among CPU times exist for the same-order methods. The proposed explicit symplectic integrators achieve the best computational efficiency compared with the other algorithms at the same order and time step. The explicit and implicit mixed symplectic methods require smaller additional computational labor than the same-order new integrators because only the solutions of r and p_r in IM2 of Equation (B2) should be iterated. Such partially implicit constructions are faster to compute than the completely implicit integrators.

Table 3.  Same as Table 1, but Dependence of Computational Cost, i.e., CPU Times (Minute: Second), on h

Method	S2	EI2	EE2	S4	EI4	EE4	RK4
h = 0.1	9:13	10:13	14:22	27:42	30:33	33:35	17:48
h = 1.0	0:56	1:03	1:26	2:46	3:09	3:21	1:46
h = 10	0:05	0:07	0:07	0:16	0:20	0:19	0:10

Download table as:  ASCIITypeset image

The major contribution of this study is the successful construction of explicit symplectic integration algorithms in general relativistic Schwarzschild-type spacetime geometries. The construction is based on an appropriate splitting form of the Hamiltonian corresponding to this spacetime. The Hamiltonian exists four integrable separable parts with analytical solutions as explicit functions of proper time. The solutions from the four parts are symmetrically composed of second- and fourth-order explicit symplectic integrators, similar to the standard explicit symplectic leapfrog methods that split the considered Hamiltonian into two integrable parts with analytical solutions as explicit functions of time. The proposed algorithms are still valid for an external magnetic field included within the vicinity of the black hole.

Numerical tests show that the newly proposed integration schemes effectively control Hamiltonian errors without secular changes when appropriate step sizes are adopted. They are well-behaved in the simulation of the long-term evolution of regular orbits with single or many loops and weakly or strongly chaotic orbits. Appropriately larger step sizes are acceptable for such explicit symplectic integrators to maintain stable or bounded energy (or Hamiltonian) errors. Explicit constructions are generally superior to same-order implicit methods in computational efficiency.

In summary, the new methods achieve long-time performance. Therefore, they are highly appropriate for the long-term numerical simulations of regular and chaotic motions of charged particles in the present non-integrable magnetized Schwarzschild spacetime background (Felice & Sorge 2003; Kološ et al. 2015; Yi & Wu 2020). The methods are also useful for studying the chaotic motion of a charged particle in a tokamak magnetic field (Cambon et al. 2014). They are suitable for investigating the capture cross-section of magnetized particles and the magnetized particles' acceleration mechanism near a black hole with an external magnetic field (Abdujabbarov et al. 2014). These methods are applicable to the simulation of the dynamics of charged particles around a regular black hole with a non-linear electromagnetic source (Jawad et al. 2016). This class of explicit symplectic integration algorithms will be developed to address other black hole gravitational problems, such as the Reissner–Nordström spacetime.

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

From an (n − 1)th step to an nth step, algorithm ${S}_{2}^{{ \mathcal H }}$ has the following discrete difference scheme:

$$\begin{eqnarray*}\begin{array}{rcl}{\theta }^{{ \mathcal H }4} & = & {\theta }_{n-1}+\displaystyle \frac{h}{2}{p}_{\theta ,n-1}{/r}_{n-1}^{2},\\ {p}_{r}^{{ \mathcal H }4} & = & {p}_{r,n-1}+\displaystyle \frac{h}{2}{p}_{\theta ,n-1}^{2}{/r}_{n-1}^{3};\\ {r}^{{ \mathcal H }3} & = & {\left[{\left({r}_{n-1}^{2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{{ \mathcal H }4}\right)}^{2}{/r}_{n-1}\right]}^{1/3},\\ {p}_{r}^{{ \mathcal H }3} & = & {p}_{r}^{{ \mathcal H }4}{\left[\left({r}_{n-1}^{2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{{ \mathcal H }4}\right){/r}_{n-1}^{2}\right]}^{1/3};\\ {r}^{{ \mathcal H }2} & = & {r}^{{ \mathcal H }3}+\displaystyle \frac{h}{2}{p}_{r}^{{ \mathcal H }3};\\ {p}_{r}^{{ \mathcal H }1} & = & {p}_{r}^{{ \mathcal H }3}+h\left[\displaystyle \frac{{{\ell }}^{2}}{{\left({r}^{{ \mathcal H }2}\right)}^{3}{\sin }^{2}{\theta }^{{ \mathcal H }4}}-\displaystyle \frac{{E}^{2}}{{\left({r}^{{ \mathcal H }2}-2\right)}^{2}}\right],\\ {p}_{\theta n} & = & {p}_{\theta ,n-1}+h\displaystyle \frac{{{\ell }}^{2}\cos {\theta }^{{ \mathcal H }4}}{{\left({r}^{{ \mathcal H }2}\right)}^{2}{\sin }^{3}{\theta }^{{ \mathcal H }4}};\\ {r}^{* { \mathcal H }2} & = & {r}^{{ \mathcal H }2}+\displaystyle \frac{h}{2}{p}_{r}^{{ \mathcal H }1};\\ {r}_{n} & = & {\left[{\left({\left({r}^{* { \mathcal H }2}\right)}^{2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{{ \mathcal H }1}\right)}^{2}{/r}^{* { \mathcal H }2}\right]}^{1/3},\\ {p}_{r}^{* { \mathcal H }3} & = & {p}_{r}^{{ \mathcal H }1}{\left[\left({\left({r}^{* { \mathcal H }2}\right)}^{2}-\displaystyle \frac{3}{2}{{hp}}_{r}^{{ \mathcal H }1}\right)/{\left({r}^{* { \mathcal H }2}\right)}^{2}\right]}^{1/3};\\ {\theta }_{n} & = & {\theta }^{{ \mathcal H }4}+\displaystyle \frac{h}{2}{p}_{\theta n}/{\left({r}_{n}\right)}^{2},\\ {p}_{{rn}} & = & {p}_{r}^{* { \mathcal H }3}+\displaystyle \frac{h}{2}{\left({p}_{\theta n}\right)}^{2}/{\left({r}_{n}\right)}^{3}.\end{array}\end{eqnarray*}$$

In this manner, the solutions (r_n, θ_n, p_rn, p_θn) at the nth step are presented. Let the integration continue from the nth step to the (n + 1)th step.

Algorithm EI4 was discussed by Lubich et al. (2010), Zhong et al. (2010), and Mei et al. (2013a, 2013b). Here, it is used to solve the Hamiltonian (15). Its construction requires splitting this Hamiltonian into two parts:

$$\begin{eqnarray}&&K={K}_{1}+{\rm{\Lambda }},\end{eqnarray} \tag{ B1 }$$

where Λ = K₂ + K₃ + K₄. The sub-Hamiltonian K₁ does not depend on momenta p_r and p_θ, and thus it is easily, explicitly, and analytically solved, and then labeled as operator ${\psi }_{h}^{{K}_{1}}$. Another sub-Hamiltonian Λ exhibits difficulty in providing explicit analytical solutions, but can be integrated using the second-order implicit midpoint rule (Feng 1986), labeled as operator IM2(h). Similar to the explicit algorithm S₂ in Equation (4), a second-order explicit and implicit mixed symplectic integrator is symmetrically composed of two explicit and implicit operators by

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

Such a mixed symplectic method demonstrates an explicit advantage over the implicit midpoint method acting on the complete Hamiltonian K in terms of computational efficiency. The fourth-order explicit and implicit mixed symplectic integrator EI4 can be obtained by substituting EI2 into ${S}_{2}^{{ \mathcal H }}$ in Equation (38).

Algorithm EE4 is based on the idea of Pihajoki (2015). Its construction relies on extending the four-dimensional phase-space variables (r, θ, p_r, p_θ) of the Hamiltonian K to eight-dimensional phase-space variables $(r,\theta ,\tilde{r},\tilde{\theta },{p}_{r},$ ${p}_{\theta },{\tilde{p}}_{r},{\tilde{p}}_{\theta })$ of a new Hamiltonian, i.e.,

$$\begin{eqnarray}&&{\rm{\Gamma }}={\kappa }_{1}(r,\theta ,{\tilde{p}}_{r},{\tilde{p}}_{\theta })+{\kappa }_{2}(\tilde{r},\tilde{\theta },{p}_{r},{p}_{\theta }),\end{eqnarray} \tag{ B3 }$$

where ${\kappa }_{1}(r,\theta ,{\tilde{p}}_{r},{\tilde{p}}_{\theta })={\kappa }_{2}(\tilde{r},\tilde{\theta },{p}_{r},{p}_{\theta })=K(r,\theta ,{p}_{r},{p}_{\theta })$. Evidently, the two sub-Hamiltonians κ₁ and κ₂ are independently, explicitly, and analytically solved, and then labeled as operators ${\psi }_{h}^{{\kappa }_{1}}$ and ${\psi }_{h}^{{\kappa }_{2}}$. The two operators are used to yield the second-order symplectic method ${{ \mathcal S }}_{2}$ and the Forest–Ruth fourth-order algorithm FR4, which are respectively given by Equations (4) and (5) but ${ \mathcal A }$ and ${ \mathcal B }$ are respectively replaced with ${\psi }^{{\kappa }_{1}}$ and ${\psi }^{{\kappa }_{2}}$.

If the two independent Hamiltonians κ₁ and κ₂ have the same initial conditions, then they should have the same solutions, i.e., $r=\tilde{r}$, $\theta =\tilde{\theta }$, ${\tilde{p}}_{r}={p}_{r}$ and ${\tilde{p}}_{\theta }={p}_{\theta }$. However, these solutions are not equal because of their couplings in the methods ${{ \mathcal S }}_{2}$ and FR4. To make them equal, Pihajoki (2015), Liu et al. (2016, 2017), Luo et al. (2017), Luo & Wu (2017), Li & Wu (2017), and Wu & Wu (2018) introduced permutations between the original variables and their corresponding extended variables after the implementation of ${{ \mathcal S }}_{2}$ or FR4. A good choice is the midpoint permutation method (Luo et al. 2017):

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal 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{ B4 }$$

By adding the midpoint permutation map ${ \mathcal M }$ after ${{ \mathcal S }}_{2}$ or FR4, Luo et al. (2017) obtained algorithms EE2 and EE4 as follows:

$$\begin{eqnarray}&&{EE}2={ \mathcal M }\otimes {{ \mathcal S }}_{2},\,\,{EE}4={ \mathcal M }\otimes {FR}4.\end{eqnarray} \tag{ B5 }$$

The inclusion of ${ \mathcal M }$ destroys the symplecticity of ${{ \mathcal S }}_{2}$ and FR4, but EE2 and EE4, similar to the symplectic schemes ${{ \mathcal S }}_{2}$ and FR4, still exhibit good long-term stable behavior in energy errors because of their symmetry. Thus, they are called explicit symplectic-like algorithms for the newly extended phase-space Hamiltonian Γ.