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 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.

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.

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$.

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.

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.

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.