# Construction of explicit symplectic integrators in general relativity. I.
Schwarzschild black holes

Ying Wang1,2, Wei Sun1, Fuyao Liu1, Xin Wu1,2,3,† 1\. School of Mathematics,
Physics and Statistics, Shanghai University of Engineering Science, Shanghai
201620, China
2\. Center of Application and Research of Computational Physics, Shanghai
University of Engineering Science, Shanghai 201620, China
3\. Guangxi Key Laboratory for Relativistic Astrophysics, Guangxi University,
Nanning 530004, China  Emails<EMAIL_ADDRESS>(Y. W.),
<EMAIL_ADDRESS>(W. S<EMAIL_ADDRESS>(F. L.); ${\dagger}$
Corresponding Author<EMAIL_ADDRESS>(X. W.)

###### Abstract

Symplectic integrators that preserve the geometric structure of Hamiltonian
flows and do not exhibit secular growth in energy errors are suitable for the
long-term integration of N-body Hamiltonian systems in the solar system.
However, the construction of explicit symplectic integrators is frequently
difficult in general relativity because all variables are inseparable.
Moreover, even if two analytically integrable splitting parts exist in a
relativistic Hamiltonian, all analytical solutions are not explicit functions
of proper time. Naturally, implicit symplectic integrators, such as the
midpoint rule, are applicable to this case. In general, these integrators are
numerically more expensive to solve than same-order explicit symplectic
algorithms. To address this issue, we split the Hamiltonian of Schwarzschild
space-time geometry into four integrable parts with analytical solutions as
explicit functions of proper time. In this manner, second- and fourth-order
explicit symplectic integrators can be easily available. The new algorithms
are also useful for modeling the chaotic motion of charged particles around a
black hole with an external magnetic field. They demonstrate excellent long-
term performance in maintaining bounded Hamiltonian errors and saving
computational cost when appropriate proper time steps are adopted.

_Unified Astronomy Thesaurus concepts_ : Black hole physics (159);
Computational methods (1965); Computational astronomy (293); Chaos (222)

##  1 Introduction

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 given to 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, four-
velocity relation of particles, and Carter constant (Carter 1968). However, no
additional information regarding the solutions but only the integrability of
space-times 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 space-times, the separation of variables in the
Hamiltonian-Jacobi equation, associated to 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;
Stuchlík $\&$ Kološ 2016; Tursunov et al. 2016; Azreg-Aïnou 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, should be 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 that 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, 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 involving the total energy. 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, nonconservative 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
the 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 secondary
part, the analytical solutions can be given in explicit functions of time. The
primary part 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 space-time, 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 have 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 but 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 in
some references (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 noncanonical because of the use of noncanonical 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 space-time backgrounds (Seyrich $\&$ Lukes-
Gerakopoulos 2012). Slimplectic integrators for general nonconservative
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 space-
times similar to the standard explicit symplectic leapfrogs for Hamiltonian
problems in solar system dynamics. If the Hamiltonians of Schwarzschild type
space-times 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 space-times in Section 4. The performance of explicit symplectic
integrators is tested numerically in Section 5. Section 6 concludes 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.

##  2 Standard explicit symplectic integrators for a separable Hamiltonian

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

$H(\mathbf{p},\mathbf{q})=H_{1}(\mathbf{p},\mathbf{q})+H_{2}(\mathbf{p},\mathbf{q}),$
(1)

where the two separable parts $H_{1}$ and $H_{2}$ are supposed to be
independently integrable. A typical splitting form of $H$ takes $H_{1}$ as
kinetic energy $T(\mathbf{p})$ and $H_{2}$ as potential $V(\mathbf{q})$.

Two differential operators are defined as follows:

$\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\sum^{N}_{i=1}(\frac{\partial
H_{1}}{\partial\mathbf{p}_{i}}\frac{\partial}{\partial\mathbf{q}_{i}}-\frac{\partial
H_{1}}{\partial\mathbf{q}_{i}}\frac{\partial}{\partial\mathbf{p}_{i}}),$
$\displaystyle\mathcal{B}$ $\displaystyle=$
$\displaystyle\sum^{N}_{i=1}(\frac{\partial
H_{2}}{\partial\mathbf{p}_{i}}\frac{\partial}{\partial\mathbf{q}_{i}}-\frac{\partial
H_{2}}{\partial\mathbf{q}_{i}}\frac{\partial}{\partial\mathbf{p}_{i}}).$

System (1) has the following formal solution

$\mathbf{Z}(h)=\mathcal{C}(h)\mathbf{Z}(0),$ (2)

where $\mathbf{Z}(0)$ denotes the value of $\mathbf{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}$:

$\mathcal{C}(h)\approx\Pi^{e}_{j=1}\mathcal{A}(h\alpha_{j})\mathcal{B}(h\beta_{j})+O(h^{d+1}),$
(3)

where coefficients $\alpha_{j}$ and $\beta_{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)

$\mathcal{S}_{2}(h)=\mathcal{A}(\frac{h}{2})\mathcal{B}(h)\mathcal{A}(\frac{h}{2}).$
(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)

$\displaystyle FR4(h)$ $\displaystyle=$
$\displaystyle\mathcal{A}(\frac{\gamma}{2}h)\mathcal{B}(\gamma
h)\mathcal{A}(\frac{1-\gamma}{2}h)\mathcal{B}((1-2\gamma)h)$ (5)
$\displaystyle\circ\mathcal{A}(\frac{1-\gamma}{2}h)\mathcal{B}(\gamma
h)\mathcal{A}(\frac{\gamma}{2}h),$

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 space-time? The succeeding discussions answer
this question.

##  3 Schwarzschild black holes

A Schwarzschild black hole with mass $M$ is a nonrotating black hole. In
spherical-like coordinates $(t,r,\theta,\phi)$, the Schwarzschild metric is
described by

$\displaystyle-c^{2}d\tau^{2}$ $\displaystyle=$ $\displaystyle
ds^{2}=g_{\alpha\beta}dx^{\alpha}dx^{\beta}$ (6) $\displaystyle=$
$\displaystyle-(1-\frac{2GM}{rc^{2}})c^{2}dt^{2}+(1-\frac{2GM}{rc^{2}})^{-1}$
$\displaystyle\cdot dr^{2}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\phi^{2},$

where $\tau$, $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\rightarrow tM$, $r\rightarrow
rM$ and $\tau\rightarrow\tau M$. In this manner, this metric is transformed
into a dimensionless form as follows:

$\displaystyle-d\tau^{2}=ds^{2}$ $\displaystyle=$
$\displaystyle-(1-\frac{2}{r})dt^{2}+(1-\frac{2}{r})^{-1}dr^{2}$ (7)
$\displaystyle+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\phi^{2}.$

This metric corresponds to a Lagrangian system

$\mathcal{L}=\frac{1}{2}(\frac{ds}{d\tau})^{2}=\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu},$
(8)

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

$p_{\mu}=\frac{\partial\mathcal{L}}{\partial\dot{x}^{\mu}}=g_{\mu\nu}\dot{x}^{\nu}.$
(9)

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

$\displaystyle p_{t}$ $\displaystyle=$
$\displaystyle-(1-\frac{2}{r})\dot{t}=-E,$ (10) $\displaystyle p_{\phi}$
$\displaystyle=$ $\displaystyle r^{2}\sin^{2}\theta\dot{\phi}=\ell,$ (11)

where $E$ and $\ell$ 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

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\mathbf{U}\cdot\mathbf{p}-\mathcal{L}=\frac{1}{2}g^{\mu\nu}p_{\mu}p_{\nu}=-\frac{1}{2}(1-\frac{2}{r})^{-1}E^{2}$
(12)
$\displaystyle+\frac{1}{2}(1-\frac{2}{r})p^{2}_{r}+\frac{1}{2}\frac{p^{2}_{\theta}}{r^{2}}+\frac{1}{2}\frac{\ell^{2}}{r^{2}\sin^{2}\theta}.$

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 disc 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 space-
time. 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^{\alpha}$ in
the Lorentz gauge is a linear combination of the time-like and space-like
axial Killing vectors $\xi^{\alpha}_{(t)}$ and $\xi^{\alpha}_{(\phi)}$
(Abdujabbarov et al. 2013; Shaymatov et al. 2015; Tursunov et al. 2016;
Benavides-Gallego et al. 2019):

$A^{\alpha}=C_{1}\xi^{\alpha}_{(t)}+C_{2}\xi^{\alpha}_{(\phi)}.$ (13)

In Felice $\&$ Sorge (2003), the constants are set as $C_{1}=0$ and
$C_{2}=B/2$. In this manner, the four-vector potential has only one nonzero
covariant component

$A_{\phi}=\frac{B}{2}g_{\phi\phi}=\frac{B}{2}r^{2}\sin^{2}\theta.$ (14)

The charged particle motion is described by the Hamiltonian system

$\displaystyle K$ $\displaystyle=$
$\displaystyle\frac{1}{2}g^{\mu\nu}(p_{\mu}-qA_{\mu})(p_{\nu}-qA_{\nu})$ (15)
$\displaystyle=$
$\displaystyle-\frac{1}{2}(1-\frac{2}{r})^{-1}E^{2}+\frac{1}{2}(1-\frac{2}{r})p^{2}_{r}+\frac{1}{2}\frac{p^{2}_{\theta}}{r^{2}}$
$\displaystyle+\frac{1}{2r^{2}\sin^{2}\theta}(L-\frac{\beta}{2}r^{2}\sin^{2}\theta)^{2},$

where $\beta=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

$L=r^{2}\sin^{2}\theta\dot{\phi}+\frac{\beta}{2}r^{2}\sin^{2}\theta.$ (16)

A point is illustrated here. The dimensionless Hamiltonian (15) is obtained
after scale transformations of $B\rightarrow B/M$, $E\rightarrow mE$,
$p_{r}\rightarrow mp_{r}$, $q\rightarrow mq$, $L\rightarrow mML$,
$p_{\theta}\rightarrow mMp_{\theta}$ and $K\rightarrow m^{2}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:

$\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle-\frac{1}{2},$ (17)
$\displaystyle K$ $\displaystyle=$ $\displaystyle-\frac{1}{2}.$ (18)

They are attributed to the four-velocity relation
$\mathbf{U}\cdot\mathbf{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 such case, no formal analytical solutions exist in the
Hamiltonian $K$.

##  4 Construction of explicit symplectic integrators for Schwarzschild space-
times

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

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\mathcal{T}+\mathcal{V},$ (19) $\displaystyle\mathcal{T}$
$\displaystyle=$
$\displaystyle\frac{1}{2}(1-\frac{2}{r})p^{2}_{r}+\frac{1}{2}\frac{p^{2}_{\theta}}{r^{2}},$
(20) $\displaystyle\mathcal{V}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(1-\frac{2}{r})^{-1}E^{2}+\frac{1}{2}\frac{\ell^{2}}{r^{2}\sin^{2}\theta}.$
(21)

The $\mathcal{V}$ part is analytically integrable, and its analytical
solutions $p_{r}$ and $p_{\theta}$ are explicit functions of proper time
$\tau$. 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 $\tau$ but are implicit
functions. In such 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
$\tau$. 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:

$\mathcal{H}=\mathcal{H}_{1}+\mathcal{H}_{2}+\mathcal{H}_{3}+\mathcal{H}_{4},$
(22)

where these sub-Hamiltonians are

$\displaystyle\mathcal{H}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{\ell^{2}}{r^{2}\sin^{2}\theta}-\frac{1}{2}(1-\frac{2}{r})^{-1}E^{2},$
(23) $\displaystyle\mathcal{H}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}p^{2}_{r},$ (24) $\displaystyle\mathcal{H}_{3}$
$\displaystyle=$ $\displaystyle-\frac{1}{r}p^{2}_{r},$ (25)
$\displaystyle\mathcal{H}_{4}$ $\displaystyle=$
$\displaystyle\frac{p^{2}_{\theta}}{2r^{2}}.$ (26)

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

$\displaystyle\frac{dp_{r}}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{\partial\mathcal{H}_{1}}{\partial
r}=\frac{\ell^{2}}{r^{3}\sin^{2}\theta}-\frac{E^{2}}{(r-2)^{2}},$ (27)
$\displaystyle\frac{dp_{\theta}}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{\partial\mathcal{H}_{1}}{\partial\theta}=\frac{\ell^{2}\cos\theta}{r^{2}\sin^{3}\theta}.$
(28)

Evidently, $r$ and $\theta$ are constants when proper time goes from
$\tau_{0}$ to $\tau_{1}=\tau_{0}+\tau$. Thus, $p_{r}$ and $p_{\theta}$ can be
solved analytically from Equations (27) and (28). They are explicit functions
of $\tau$ in the following forms

$\displaystyle p_{r}(\tau)$ $\displaystyle=$ $\displaystyle
p_{r0}+\tau[\frac{\ell^{2}}{r^{3}_{0}\sin^{2}\theta_{0}}-\frac{E^{2}}{(r_{0}-2)^{2}}],$
(29) $\displaystyle p_{\theta}(\tau)$ $\displaystyle=$ $\displaystyle
p_{\theta 0}+\tau\frac{\ell^{2}\cos\theta_{0}}{r^{2}_{0}\sin^{3}\theta_{0}},$
(30)

where $r_{0}$, $\theta_{0}$, $p_{r0}$ and $p_{\theta 0}$ represent values of
$r$, $\theta$, $p_{r}$ and $p_{\theta}$ at the proper time $\tau_{0}$; and
$p_{r}(\tau)$ and $p_{\theta}(\tau)$ denote the values of $p_{r}$ and
$p_{\theta}$ at proper time $\tau_{1}$. A differential operator for solving
$\mathcal{H}_{1}$ is labeled as $\psi^{\mathcal{H}_{1}}_{\tau}$.

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

$\displaystyle\mathcal{H}_{2}:~{}\frac{dr}{d\tau}$ $\displaystyle=$
$\displaystyle p_{r},~{}~{}\dot{p}_{r}=0;$ (31)
$\displaystyle\mathcal{H}_{3}:~{}\frac{dr}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{2}{r}p_{r},~{}\frac{dp_{r}}{d\tau}=-\frac{p^{2}_{r}}{r^{2}};$
(32) $\displaystyle\mathcal{H}_{4}:~{}\frac{d\theta}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{p_{\theta}}{r^{2}},~{}\frac{dp_{r}}{d\tau}=\frac{p^{2}_{\theta}}{r^{3}},~{}\dot{r}=\dot{p}_{\theta}=0.$
(33)

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

$\displaystyle\psi^{\mathcal{H}_{2}}_{\tau}:~{}r(\tau)$ $\displaystyle=$
$\displaystyle r_{0}+\tau p_{r0};$ (34)
$\displaystyle\psi^{\mathcal{H}_{3}}_{\tau}:~{}r(\tau)$ $\displaystyle=$
$\displaystyle[(r^{2}_{0}-3\tau p_{r0})^{2}/r_{0}]^{1/3},$ $\displaystyle
p_{r}(\tau)$ $\displaystyle=$ $\displaystyle p_{r0}[(r^{2}_{0}-3\tau
p_{r0})/r^{2}_{0}]^{1/3};$ (35)
$\displaystyle\psi^{\mathcal{H}_{4}}_{\tau}:~{}\theta(\tau)$ $\displaystyle=$
$\displaystyle\theta_{0}+\tau p_{\theta 0}/r^{2}_{0},$ $\displaystyle
p_{r}(\tau)$ $\displaystyle=$ $\displaystyle p_{r0}+\tau p^{2}_{\theta
0}/r^{3}_{0}.$ (36)

It is clear that these solutions are explicit functions of proper time $\tau$.
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$, $\theta$ and $p_{r}$ for the sum cannot be
expressed as explicit functions of proper time $\tau$. 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^{\mathcal{H}}_{h}$ of the Hamiltonian (12) over time step $h$
is approximately given by the symmetric composition of these operators

$\displaystyle\psi^{\mathcal{H}}_{h}\approx S^{\mathcal{H}}_{2}(h)$
$\displaystyle=$
$\displaystyle\psi^{\mathcal{H}_{4}}_{h/2}\circ\psi^{\mathcal{H}_{3}}_{h/2}\circ\psi^{\mathcal{H}_{2}}_{h/2}\circ\psi^{\mathcal{H}_{1}}_{h}$
(37)
$\displaystyle\circ\psi^{\mathcal{H}_{2}}_{h/2}\circ\psi^{\mathcal{H}_{3}}_{h/2}\circ\psi^{\mathcal{H}_{4}}_{h/2}.$

The above construction is a second order explicit symplectic integrator marked
as $S^{\mathcal{H}}_{2}$. 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

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

where $\delta=1-2\gamma$.

The Hamiltonian (15) exhibits the following splitting form

$K=K_{1}+K_{2}+K_{3}+K_{4},$ (39)

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

$\displaystyle K_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2r^{2}\sin^{2}\theta}(L-\frac{\beta}{2}r^{2}\sin^{2}\theta)^{2}$
(40) $\displaystyle-\frac{1}{2}(1-\frac{2}{r})^{-1}E^{2}.$

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

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

Table 1: Dependence of stable (S) or unstable (U) behavior of Hamiltonian errors for the seven algorithms on step size $h$. Chaotic Orbit 3 in Figure 2 is integrated until proper time $\tau=10^{8}$. 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
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
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

##  5 Numerical evaluations

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 $\beta=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 $\ell$
(or $L$) =4.6, and proper time step size $h=1$. Initial conditions are $r=11$,
$\theta=\pi/2$ and $p_{r}=0$. The initial value of $p_{\theta}$ ($>0$) is
determined by 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^{-6}$ to Hamiltonian errors $\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^{-8}$.

The four fourth-order algorithms, namely, S4, EI4, EE4 and RK4, yield the
Hamiltonian errors in Figures 1(b) and 1(c). The algorithms S4, EI4 and EE4
are accurate to an order of $10^{-8}$. 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^{-6}$ 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 $\theta=\pi/2$ and $p_{\theta}>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 the
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, respectively. 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 $\beta\neq 0$

When an external magnetic field with parameter $\beta=8.9\times 10^{-4}$ 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.

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)-2(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 such case, accuracy is maintained with
an order of $10^{-5}$. 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\epsilon$, where
$\epsilon\sim 10^{-16}$ 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 in an order of $10^{-12}$. The
total errors in the energy are stabilized at the order of magnitude when
$N<10^{4}$, but grow linearly as $N\gg 10^{4}$. If a symplectic method has a
truncation energy error higher than the order of $10^{-8}$, then the total
errors in the energy remain bounded and approach the order of truncation
errors when $N<10^{8}$, whereas increase linearly as $N\gg 10^{8}$. 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 in the order of
$10^{-9}$ for $h=1$ but the roundoff errors are $10^{-8}$ after $10^{8}$
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 in the order of $10^{-5}$ for $h=10$. They are
larger than the roundoff errors after $10^{8}$ 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.

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 (B.2) should be iterated. Such
partially implicit constructions are faster to compute than the completely
implicit integrators.

##  6 Conclusions

The major contribution of this study is the successful construction of
explicit symplectic integration algorithms in general relativistic
Schwarzschild type space-time geometries. The construction is based on an
appropriate splitting form of the Hamiltonian corresponding to this space-
time. 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 space-time 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 nonlinear electromagnetic source (Jawad et al.
2016). Such class of explicit symplectic integration algorithms will be
developed to address other black hole gravitational problems, such as the
Reissner-Nordström space-time.

## APPENDIX

##  Appendix A Discrete difference scheme of algorithm $S^{\mathcal{H}}_{2}$

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

$\displaystyle\theta^{\mathcal{H}4}$ $\displaystyle=$
$\displaystyle\theta_{n-1}+\frac{h}{2}p_{\theta,n-1}/r^{2}_{n-1},$
$\displaystyle p^{\mathcal{H}4}_{r}$ $\displaystyle=$ $\displaystyle
p_{r,n-1}+\frac{h}{2}p^{2}_{\theta,n-1}/r^{3}_{n-1};$ $\displaystyle
r^{\mathcal{H}3}$ $\displaystyle=$
$\displaystyle[(r^{2}_{n-1}-\frac{3}{2}hp^{\mathcal{H}4}_{r})^{2}/r_{n-1}]^{1/3},$
$\displaystyle p^{\mathcal{H}3}_{r}$ $\displaystyle=$ $\displaystyle
p^{\mathcal{H}4}_{r}[(r^{2}_{n-1}-\frac{3}{2}hp^{\mathcal{H}4}_{r})/r^{2}_{n-1}]^{1/3};$
$\displaystyle r^{\mathcal{H}2}$ $\displaystyle=$ $\displaystyle
r^{\mathcal{H}3}+\frac{h}{2}p^{\mathcal{H}3}_{r};$ $\displaystyle
p^{\mathcal{H}1}_{r}$ $\displaystyle=$ $\displaystyle
p^{\mathcal{H}3}_{r}+h[\frac{\ell^{2}}{(r^{\mathcal{H}2})^{3}\sin^{2}\theta^{\mathcal{H}4}}-\frac{E^{2}}{(r^{\mathcal{H}2}-2)^{2}}],$
$\displaystyle p_{\theta n}$ $\displaystyle=$ $\displaystyle
p_{\theta,n-1}+h\frac{\ell^{2}\cos\theta^{\mathcal{H}4}}{(r^{\mathcal{H}2})^{2}\sin^{3}\theta^{\mathcal{H}4}};$
$\displaystyle r^{*\mathcal{H}2}$ $\displaystyle=$ $\displaystyle
r^{\mathcal{H}2}+\frac{h}{2}p^{\mathcal{H}1}_{r};$ $\displaystyle r_{n}$
$\displaystyle=$
$\displaystyle[((r^{*\mathcal{H}2})^{2}-\frac{3}{2}hp^{\mathcal{H}1}_{r})^{2}/r^{*\mathcal{H}2}]^{1/3},$
$\displaystyle p^{*\mathcal{H}3}_{r}$ $\displaystyle=$ $\displaystyle
p^{\mathcal{H}1}_{r}[((r^{*\mathcal{H}2})^{2}-\frac{3}{2}hp^{\mathcal{H}1}_{r})/(r^{*\mathcal{H}2})^{2}]^{1/3};$
$\displaystyle\theta_{n}$ $\displaystyle=$
$\displaystyle\theta^{\mathcal{H}4}+\frac{h}{2}p_{\theta n}/(r_{n})^{2},$
$\displaystyle p_{rn}$ $\displaystyle=$ $\displaystyle
p^{*\mathcal{H}3}_{r}+\frac{h}{2}(p_{\theta n})^{2}/(r_{n})^{3}.$

In this manner, the solutions $(r_{n},\theta_{n},p_{rn},p_{\theta n})$ at the
$n$th step are presented. Let the integration continue from the $n$th step to
the $(n+1)$th step.

##  Appendix B Descriptions of algorithms EI4 and EE4

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

$K=K_{1}+\Lambda,$ (B1)

where $\Lambda=K_{2}+K_{3}+K_{4}$. The sub-Hamiltonian $K_{1}$ does not depend
on momenta $p_{r}$ and $p_{\theta}$, and thus, it is easily, explicitly and
analytically solved, and then labeled as operator $\psi^{K_{1}}_{h}$. Another
sub-Hamiltonian $\Lambda$ 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_{2}$ in Equation (4), a second-order explicit and implicit mixed symplectic
integrator is symmetrically composed of two explicit and implicit operators by

$EI2(h)=\psi^{K_{1}}_{h/2}\circ IM2(h)\circ\psi^{K_{1}}_{h/2}.$ (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 four-order explicit and implicit mixed
symplectic integrator EI4 can be obtained by substituting EI2 into
$S^{\mathcal{H}}_{2}$ 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,\theta,p_{r},p_{\theta})$ 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.,

$\Gamma=\kappa_{1}(r,\theta,\tilde{p}_{r},\tilde{p}_{\theta})+\kappa_{2}(\tilde{r},\tilde{\theta},p_{r},p_{\theta}),$
(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 $\kappa_{1}$ and $\kappa_{2}$ are
independently, explicitly and analytically solved, and then labeled as
operators $\psi^{\kappa_{1}}_{h}$ and $\psi^{\kappa_{2}}_{h}$. 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 $\kappa_{1}$ and $\kappa_{2}$ 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), Luo et al. (2017), Liu 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):

$\displaystyle\mathcal{M}:\frac{r+\tilde{r}}{2}$ $\displaystyle\rightarrow$
$\displaystyle
r=\tilde{r},~{}~{}~{}~{}~{}~{}~{}~{}\frac{\theta+\tilde{\theta}}{2}\rightarrow\theta=\tilde{\theta};$
$\displaystyle\frac{p_{r}+\tilde{p}_{r}}{2}$ $\displaystyle\rightarrow$
$\displaystyle
p_{r}=\tilde{p}_{r},~{}~{}\frac{p_{\theta}+\tilde{p}_{\theta}}{2}\rightarrow
p_{\theta}=\tilde{p}_{\theta}.$ (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:

$EE2=\mathcal{M}\otimes\mathcal{S}_{2},~{}~{}EE4=\mathcal{M}\otimes FR4.$ (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 $\Gamma$.

## Acknowledgments

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 (C0035736), 11803020, 41807437, U2031145]
and the Natural Science Foundation of Guangxi (Grant Nos. 2018GXNSFGA281007
and 2019JJD110006).

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phy. Rev. Lett., 116, 061102
  * Abdujabbarov et al. (2013) Abdujabbarov, A. A., Ahmedov, B. J., & Jurayeva, N. B. 2013, Phys. Rev. D, 87, 064042
  * Abdujabbarov et al. (2014) Abdujabbarov, A., Ahmedov, B., Rahimov, O., & Salikhbaev., U. 2014, Phys. Scr., 89, 084008
  * Avdyushev (2003) Avdyushev, E. A. 2003, Celestial Mechanics and Dynamical Astronomy, 87, 383
  * Azreg-Aïnou (2016) Azreg-Aïnou, M. 2016, Eur. Phys. J. C, 76, 414
  * Bacchini et al. (2018a) Bacchini, F., Ripperda, B., Chen, A. Y., & Sironi, L. 2018a, Astropys. J. Suppl., 237, 6
  * Bacchini et al. (2018b) Bacchini, F., Ripperda, B., Chen, A. Y., & Sironi, L. 2018b, Astropys. J. Suppl., 240, 40
  * Baumgarte (1972) Baumgarte, J. 1972, Comp. Math. Appl. Mech. Eng., 1, 1
  * Baumgarte (1973) Baumgarte, J. 1973, Celest. Mech., 5, 490
  * Benavides-Gallego et al. (2019) Benavides-Gallego, C. A., Abdujabbarov, A., Malafarina, D., Ahmedov, B., & Bambi, C. 2019, Phys. Rev. D, 99, 044012
  * Borm & Spaans (2013) Borm, C. V., & Spaans, M. 2013, Astron. Astrophys., 553, L9
  * Brown (2006) Brown, J. D. 2006, Phys. Rev. D, 73, 024001
  * Cambon et al. (2014) Cambon, B., Leoncini, X., Vittot, M., Dumont, R., & Garbet, X. 2014, Chaos, 24, 033101
  * Carter (1968) Carter, B. 1968, Phy. Rev., 174, 1559
  * Chambers (2000) Chambers, J. E., & Murison, M. A. 2000, AJ, 119, 425
  * Chin (1997) Chin, S. A. 1997, Phys. Lett. A, 226, 344
  * Deng et al. (2020) Deng, C., Wu, X., & Liang, E. 2020, MNRAS, 496, 2946
  * Collaboration et al. (2019) EHT Collaboration et al. 2019, ApJL, 875, L1
  * Einstein (1915) Einstein, A. 1915, Sitzungsberichte der Köiglich Preu$\beta$schen Akademie der Wissenschaften (Berlin: Deutsche Akademie der Wissenschaften zu Berlin)
  * Einstein & Sitzungsber (1996) Einstein, A., & Sitzungsber, K. 1996, Preuss. Akad. Wiss., 1, 688
  * Felice & Sorge (2003) Felice, D. d, & Sorge, F. 2003, Class. Quantum Grav., 20, 469
  * Feng (1986) Feng, K. 1986, Journal of Computational Mathematics, 44, 279
  * Feng et al. (2009) Feng, K., & Qin, M. Z. 2009, Symplectic Geometric Algorithms for Hamiltonian Systems (Hangzhou, New York: Zhejiang Science and Technology Publishing House, Springer)
  * Forest & Ruth (1990) Forest, E., & Ruth, R. D. 1990, Physica D, 43, 105
  * Fukushima (2003a) Fukushima, T. 2003a, AJ, 126, 1097
  * Fukushima (2003b) Fukushima, T. 2003b, AJ, 126, 2567
  * Fukushima (2003c) Fukushima, T. 2003c, AJ, 126, 3138
  * Fukushima (2004) Fukushima, T. 2004, AJ, 127, 3638
  * Hairer et al. (1999) Hairer E., Lubich C. & Wanner G. 1999, Geometric Numerical Integration, Springer-Verlag, Berlin
  * Hénon & Heiles (1964) Hénon, M., & Heiles, C. 1964, AJ, 69, 73
  * Hernandez & Dehnen (2017) Hernandez, D. M., & Dehnen, W. 2017, MNRAS, 468, 2614
  * Hu et al. (2019) Hu S., Wu X., Huang G., & Liang E. 2019, ApJ, 887, 191
  * Huang & Innanen (1983) Huang, T. Y. & Innanen, K. 1983, AJ, 88, 870
  * Jawad et al. (2016) Jawad, A., Ali, F., Jamil, M., & Debnath, U. 2016, Commun. Theor. Phys., 66, 509
  * Kerr (1963) Kerr, R. P. 1963, Phy. Rev. Lett., 11, 237
  * Kološ et al. (2015) Kološ, M., Stuchlík, Z., & Tursunov, A. 2015, Class. Quantum Grav., 32, 165009
  * Kopáček et al. (2010) Kopáček, O., Karas, V., Kovář, J., & Stuchlík, Z. 2010, ApJ, 722, 1240
  * Kopáček & Karas (2014) Kopáček, O., & Karas, V. 2014, ApJ, 787, 117
  * Laskar & Robutel (2001) Laskar J., & Robutel, P. 2001, Celest. Mech. Dyn. Astron., 80, 39
  * Li & Wu (2017) Li D., & Wu, X. 2017, Mon. Not. R. Astron. Soc., 469, 3031
  * Li & Wu (2019) Li D., & Wu, X. 2019, Eur. Phys. J. Plus, 134, 96
  * Liao (1997) Liao, X. H. 1997, Celest. Mech. Dyn. Astron., 66, 243
  * Lichtenberg & Lieberman (1983) Lichtenberg, A. J., & Lieberman, M. A. 1983, Regular and Chaotic Dynamics (Springer-Verlag, New York)
  * Liu et al. (2016) Liu, L., Wu, X., Huang, G. Q., & Liu, F. 2016, Mon. Not. R. Astron. Soc., 459, 1968
  * Liu et al. (2017) Liu, L., Wu, X., & Huang, G. Q. 2017, Gen. Relativ. Gravit., 49, 28
  * Lubich et al. (2010) Lubich, C., Walther, B., & Brügmann, B. 2010, Phys. Rev. D, 81, 104025
  * Luo et al. (2017) Luo, J., Wu, X., Huang, G., & Liu, F. 2017, ApJ, 834, 64
  * Luo & Wu (2017) Luo, J., & Wu, X. 2017, Eur. Phys. J. Plus, 132, 485
  * Ma et al. (2008a) Ma, D. Z., Wu, X., & Zhong, S. Y. 2008a, ApJ, 687, 1294
  * Ma et al. (2008b) Ma, D. Z., Wu, X., & Zhu, J. F. 2008b, New Astrom., 13, 216
  * McLachlan & Atela (1992) McLachlan, R. I., & Atela, P. 1992, Nonlinearity, 5, 541
  * Mei et al. (2013a) Mei, L., Ju, M., Wu, X., & Liu, S. 2013a, Mon. Not. R. Astron. Soc., 435, 2246
  * Mei et al. (2013b) Mei, L., Wu, X., & Liu, F. 2013b, Eur. Phys. J. C, 73, 2413
  * Nacozy (1971) Nacozy, P. E. 1971, Astrophys. Space Sci, 14, 40
  * Omelyan et al. (2002a) Omelyan, I. P., Mryglod, I. M., & Folk, R. 2002a, Comput. Phys. Commun., 146, 188
  * Omelyan et al. (2002b) Omelyan, I. P., Mryglod, I. M., & Folk, R. 2002b, Phys. Rev. E, 66, 026701
  * Omelyan et al. (2003) Omelyan, I. P., Mryglod, I. M., & Folk, R. 2003, Comput. Phys. Commun., 151, 272
  * Pihajoki (2015) Pihajoki, P. 2015, Celest. Mech. Dyn. Astron., 121, 211
  * Preto & Saha (2009) Preto, M., & Saha, P. 2009, ApJ, 703, 1743
  * Rein & Spiegel (2015) Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424
  * Ruth (1983) Ruth, R. D, 1983, IEEE Trans. Nucl. Sci. NS 30, 2669
  * Schwarzschild (1916) Schwarzschild, K. 1916, Sitzungsberichte der Köiglich Preu$\beta$schen Akademie der Wissenschaften (Berlin: Deutsche Akademie der Wissenschaften zu Berlin)
  * Seyrich & Lukes-Gerakopoulos (2012) Seyrich, J., & Lukes-Gerakopoulos, G. 2012, Phys. Rev. D, 86, 124013
  * Shaymatov et al. (2015) Shaymatov, S., Patil, M., Ahmedov, B., & Joshi, P. S. 2015, Phys. Rev. D, 91, 064025
  * Stuchlík & Kološ (2016) Stuchlík, Z., & Kološ, M. 2016, Eur. Phys. J. C, 76, 32
  * Suzuki (1991) Suzuki, M. 1991, J. Math. Phys. (N.Y.), 32, 400
  * Swope et al. (1982) Swope, W. C., Andersen, H. C., Berens, P. H., & Wilson, K. R. 1982, J. Chem. Phys. 76, 637
  * Takahashi & Koyama (2012) Takahashi, M., & Koyama, H. 2009, ApJ, 693, 472
  * Tao (2016) Tao, M. 2016, J. Comput. Phys., 327, 245
  * Tsang et al. (2015) Tsang, D., Galley, C. R., Stein, L. C., & Turner, A. 2015, ApJL, 809, L9
  * Tursunov et al. (2016) Tursunov, A., Stuchlík, Z., & Kološ, M. 2016, Phys. Rev. D, 93, 084012
  * Wang et al. (2016) Wang, S. C., Wu, X., & Liu, F. Y. 2016, MNRAS, 463, 1352
  * Wang et al. (2018) Wang, S. C., Huang, G. Q., & Wu, X. 2018, AJ, 155, 67
  * Wisdom (1982) Wisdom, J. 1982, AJ, 87, 577
  * Wisdom & Holman (1991) Wisdom, J., & Holman, M. 1991, AJ, 102, 1528
  * Wisdom et al. (1996) Wisdom, J., Holman, M., & Touma, J. 1996, Fields Inst. Commun., 10, 217
  * Wu et al. (2003) Wu, X., Huang, T. Y., & Wan, X. S. 2003, Chinese Astronomy and Astrophysics, 27, 114
  * Wu et al. (2007) Wu, X., Huang, T. Y., Wan, X. S., & Zhang, H. 2007, AJ, 133, 2643
  * Wu & Xie (2010) Wu, X., & Xie, Y. 2010, Phys. Rev. D, 81, 084045
  * Wu & Wu (2018) Wu, Y. L., & Wu, X. 2018, International Journal of Modern Physics C, 29, 1850006
  * Wu & Wu (2020) Yi, M., & Wu, X. 2020, Phys. Scr., Phys. Scr., 95, 085008
  * Yoshida (1990) Yoshida, H. 1990, Phys. Lett. A, 150, 262
  * Zhong et al. (2010) Zhong, S. Y., Wu, X., Liu, S. Q., & Deng, X. F. 2010, Phys. Rev. D, 82, 124040

Figure 1: (a)-(c) Hamiltonian errors $\Delta H=1+2\mathcal{H}$ from Eq. (17)
for several algorithms solving the Schwarzschild problem (12). The adopted
algorithms are the new second-order explicit symplectic integrator S2 in
Equation (37), the second-order explicit and implicit mixed symplectic method
EI2 in Equation (B.2), the second-order explicit extended phase-space
symplectic-like algorithm EE2, the new fourth-order explicit symplectic
integrator $S_{4}$ in Equation (38), the fourth-order explicit and implicit
mixed symplectic method EI4, the fourth-order explicit extended phase-space
symplectic-like algorithm EE4 in Equation (B.5), and the fourth-order Runge-
Kutta scheme RK4. The energy and angular momentum of particles are $E=0.995$
and $\ell$ (or $L)$=4.6, respectively, and the proper time-step is $h=1$. The
integrated orbit (called Orbit 1) has initial conditions $r=11$,
$\theta=\pi/2$ and $p_{r}=0$. The initial value of $p_{\theta}$ $(>0)$ is
given by Equation (17). (d) Poincaré sections on the plane $\theta=\pi/2$ and
$p_{\theta}>0$. Apart from Orbit 1, Orbits 2 and 3 with initial separations
$r=70$ and 110, respectively, are plotted. The initial values of $\theta$ and
$p_{r}$ for Orbits 2 and 3 are the same as those for Orbit 1. The three orbits
are regular tori because of the integrability of the system (12).

Figure 2: Same as Figure 1, but an external magnetic field with parameter
$\beta=8.9\times 10^{-4}$ is included within the vicinity of the black hole.
(a) Poincaré sections. Orbit 1 is still a regular torus, Orbit 2 has many
islands, and Orbit 3 is strongly chaotic. (b)-(d) Hamiltonian errors $\Delta
K=1+2K$ from Equation (18) for the algorithms solving the three orbits in the
system (15).