Application of Explicit Symplectic Integrators in a Magnetized Deformed Schwarzschild Black Spacetime

Following the latest work of Wu et al., we construct time-transformed explicit symplectic schemes for a Hamiltonian system on the description of charged particles moving around a deformed Schwarzschild black hole with an external magnetic field. Numerical tests show that such schemes have good performance in stabilizing energy errors without secular drift. Meantime, tangent vectors are solved from the variational equations of the system with the aid of an explicit symplectic integrator. The obtained tangent vectors are used to calculate several chaos indicators, including Lyapunov characteristic exponents, fast Lyapunov indicators, a smaller alignment index, and a generalized alignment index. It is found that the smaller alignment index and generalized alignment index are the fastest indicators for distinguishing between regular and chaotic cases. The smaller alignment index is applied to explore the effects of the parameters on the dynamical transition from order to chaos. When the positive deformation factor and angular momentum decrease, or when the energy, positive magnetic parameter, and the magnitude of the negative deformation parameter increase, chaos easily occurs for the appropriate choices of initial conditions and the other parameters.


Introduction
Abundant highly nonlinear general relativistic dynamical models greatly enrich the content of chaotic dynamics. A double fixed-center problem can be integrable and nonchaotic in Newtonian mechanics. However, in the theory of relativity, it has strong chaotic behavior (Vieira & Letelier 1997;Contopoulos et al. 1999). This example shows a significant difference between the dynamics of general relativity and the dynamics of Newton's theory. Many scholars have focused on studying the chaotic phenomena in the black hole disk (or halo or ring) system (Vieira & Letelier 1999;Semerak & Sukova 2010). Takahashi & Koyama (2009) studied the dynamics of the charged particles around the Kerr black hole with an electromagnetic field, and found that the drag effect of the Kerr spacetime could weaken chaos. In the work on detecting gravitational waves, the post-Newtonian dynamics of two-and three-body problems of compact objects have been particularly focused upon (see, e.g., Levin 2000;Wu & Xie 2007Wang & Wu 2011;Huang & Wu 2014;Huang et al. 2016;Li et al. , 2020Li et al. , 2021.
Because the equations of motion in the theory of relativity are highly nonlinear and nonintegrable in general, they are often solved using reliable numerical integration methods. As far as the long-term evolution of Hamiltonian systems is concerned, geometric integrators have the advantage of preserving structure or motion integrals, and their use has been widely adopted. The geometric integrators mainly include symplectic algorithms, manifold correction methods, extended phase-space symplectic-like methods, and energy conservation algorithms (Huang & Mei 2020;Hu et al. 2021aHu et al. , 2021b. The manifold correction method of Nacozy (1971) is frequently developed.  used the least-squares method to design a velocity scaling factor, which approximately satisfies the Kepler energy relationship of individual celestial bodies. Ma et al. (2008) employed a scaling factor to act on the velocity, so that the scaling factor strictly meets the Kepler energy relationship of individual celestial bodies. Zhong & Wu (2009) proposed a speed-scale factor method with leastsquares correction for several motion integrals. The manifold correction methods are convenient for treating conservative N-body problems in the solar system , nonconversitive systems (Wang et al. 2016(Wang et al. , 2018, and relativistic post-Newtonian systems of spinning compact binaries (Zhong et al. 2010).
Symplectic integrators have been widely popular in celestial mechanics because they achieve good long-term stability. They are divided into explicit and implicit symplectic schemes. If a Hamiltonian system can be separated into momentum and coordinate terms, then constructing an explicit symplectic scheme is easy (Zhang et al. 2021). However, such a separation, for most Hamiltonian systems, such as the compact binary Hamiltonian system (Wu & Huang 2015), is impossible. Of course, implicit symplectic schemes are always employed (Yoshida 1990). Compared with the completely implicit methods, explicit and implicit combined symplectic integrators (Zhong et al. 2010;Mei et al. 2013aMei et al. , 2013b have higher computational efficiency. Pihajoki (2015) proposed an expanded phase-space method to construct explicit algorithms for inseparable Hamiltonian systems. Such algorithms are similar to the symplectic algorithms. Liu et al. (2016) stabilized the error of a leapfrog symplectic-like scheme through the sequent permutations of coordinates and momenta. Luo et al. (2017) found that the energy conservation and efficiency of the midpoint permutations are better than those of the coordinate and momentum permutations. In addition, the extended phase-space method with the midpoint permutations was applied to solve logarithmic Hamiltonian problems (Li & Wu 2017) and post-Newtonian Lagrangian problems (Pan et al. 2021). Recently, Wang et al. (2021aWang et al. ( , 2021bWang et al. ( , 2021c proposed explicit symplectic schemes by dividing Hamiltonian systems of nonrotating black holes into several sub-Hamiltonian systems with analytical solutions as explicit functions of proper time. However, such constructions are unsuitable for all Hamiltonian problems of black holes. For example, there is an obstacle in constructing the explicit symplectic scheme for the Kerr spacetime. Wu et al. (2021) introduced a time-transformed Hamiltonian to the Kerr metric and successfully constructed explicit symplectic integrators for the time-transformed Hamiltonian.  and  used the time-transformed explicit symplectic methods to study the dynamics of charged particles moving around the Kerr black hole.
The chaotic behavior of Newtonian dynamic systems can be identified by using different methods, such as the Poincaré section method, the Lyapunov index, the local Lyapunov index and its spectral distribution (Contopoulos et al. 1999;Szezech et al. 2005), the fast Lyapunov index (Froeschle et al. 1997;Froeschle & Lega 2000), the spectrum analysis method (Laskar 1994), a relatively limited-time Lyapunov index (Sandor et al. 2004), a smaller alignment index (Skokos et al. 2004), a power spectrum (Binney & Spergel 1982), and a fractal graph (Levin 2000). Each of these methods has its advantages and disadvantages, and should be selected carefully. The Poincaré section method can intuitively describe the phase-space structure, but it is only suitable for dynamic systems in which the phase space minus the number of motion integrals does not exceed 3. The Lyapunov index is important for measuring the separation ratio of two adjacent orbits. It is suitable for any dimensional systems, but it usually requires a long calculation time. The fast Lyapunov index and smaller alignment index are methods for quickly identifying chaos. The latter is more sensitive than the former.
Chaos is a physical phenomenon inherent in a system itself, which is independent of the choice of the coordinate system. An observable or physical quantity in general relativity should be defined as tensor or scalar. Similarly, the identification of ordered and chaotic indicators of the relativistic system should be independent of the choice of spacetime coordinates. The invariant Lyapunov index in the theory of relativity (Wu & Huang 2003) and the invariant fast Lyapunov index  can satisfy this requirement.
In the present work, a time-transformed explicit symplectic integrator from the idea of Wu et al. (2021) is used to study the chaotic dynamics of charged test particles moving around a deformed Schwarzschild black hole with an external magnetic field (Rayimbaev 2016). Explicit symplectic methods are also applied to the variational equations of the problem. The present paper is unlike the work of Yi & Wu (2020), in which the Poincaré section method was used to distinguish chaos from order in the system. The rest of this paper is organized as follows. In Section 2, an explicit symplectic scheme (S4) of the Hamiltonian describing the deformed Schwarzschild black hole and an explicit symplectic scheme S v 4 ( ) of the Hamilton variational equations are constructed. In Section 3, we test the performance of the explicit symplectic scheme S4, and use S v 4 to calculate a smaller alignment index (SALI), generalized alignment index (GALI), and Lyapunov characteristic exponents (LCE). Conclusions are then presented in Section 4.

Deformed Schwarzschild Black Hole
In Schwarzschild coordinates (t, r, θ, f), a nonspinning deformed Schwarzschild black hole with mass M has a line element (Rayimbaev 2016) Here, the speed of light c and the constant of gravity G are geometric units, c = G = 1. In addition, f = 1 − 2M/r, h = ε/r 3 is a deformation part, and the deformation factor ε indicates whether a black hole is deformed. If the deformation factor ε = 0, then this metric is the Schwarzschild metric; if ε ≠ 0, then this metric is a Schwarzschild-like metric in alternative theories of gravity. Dimensionless operations of t → tM, r → rM, τ → τM, and ε → εM 3 are performed, where τ denotes the proper time.
Suppose that the black hole is immersed in an electromagnetic field, which is described by a four-dimensional potential (Abdujabbarov et al. 2013): where B is a uniform magnetic field. The Hamiltonian of a nonspinning particle with a charge q moving around the black hole plus the external electromagnetic field is where a covariant generalized momentum scheme is defined as follows: In the above equation, x U  = m is a four-velocity. Based on the Hamiltonian equations, two constant momentum components are where Q = Bq, and E and l are the energy and angular momentum of the test particles moving around the black hole, respectively. The Hamiltonian (3) is rewritten as (7) is the Hamiltonian describing the motions of charged test particles around the deformed Schwarzschild black hole. The dimensionless operations are given as follows: where m is the mass of a test particle. Due to the rest mass, the Hamiltonian (7) is a given constant:

Explicit Symplectic Methods for Hamiltonian Equations
Because h is a function of r, it is difficult to use the Hamiltonian (7) to construct explicit symplectic integrators similar to those in the Schwarzschild-type spacetimes. Wu et al. (2021) successfully solved this problem by means of time transformations. Now, we construct a time-transformation Hamiltonian associated to . Mikkola & Tanikawa (1999) introduced a time-transformation function where g is a function of r and θ, and ω is the new time coordinate. The proper time is referred to a new coordinate τ = q 0 , and the corresponding momentum is p The new Hamiltonian H and the original Hamiltonian  are equivalent. Note that H = 0.
According to the Hamiltonian splitting forms introduced by Wang et al. (2021aWang et al. ( , 2021bWang et al. ( , 2021c, the Hamiltonian (11) can be split into the following sub-Hamiltonian parts , . 14 After obtaining the time-dependent analytical solutions of all sub-Hamiltonian parts, we can construct a second-order explicit symplectic (S2) form as follows: where w is a new coordinate time step. Then, a fourth-order explicit symplectic scheme (S4) (Yoshida 1990 Methods S2 and S4 are the application of the time-transformed explicit symplectic integrators proposed by Wu et al. (2021) in the present problem.
In order to test the performance of this explicit symplectic scheme, we will compare the implicit symplectic scheme iS4 with the explicit symplectic scheme. The Hamiltonian (11) (16), H A y and H B y are analytical solutions of H A and H B . As such, we can get second-order and fourth-order implicit leapfrogging symplectic schemes:

Symplectic Schemes of Hamiltonian Variational Equations
The phase-space variables of the Hamiltonian (11) are represented by a vector where q 1 = q 0 , q 2 = q r , q 3 = q θ , p 1 = p 0 , p 2 = p r , and p 3 = p θ , The equation of motion for the Hamiltonian (11) is also written where 0 3 is a 3 × 3 matrix with all of its elements equal to zero, and I 3 is a 3 × 3 identity matrix.
Then, a second-order explicit symplectic scheme S v 2 and a fourthorder explicit symplectic scheme S v 4 are constructed as follows:

Numerical Calculations
Now, we use the above fourth-order explicit symplectic scheme to study the dynamics of charged particles. If ε = 0, Equation (1) is the Schwarzschild metric. This metric has been investigated in many studies.  used the extended phase-space method to study the motion of charged particles around a black hole with a magnetic field. If ε ≠ 0, then the black hole is deformed.
In Figure 1(a), we use the RKF8(9) method to plot the Poincaré sections with four orbits on the plane p , 0 2 q = > p q . This figure describes the intersection points of the trajectories of the particles on this plane. If these points of an orbit are regularly encircled in one or several circles, then this orbit is ordered. If these points are irregularly and randomly filled with an area, then the orbit is chaotic. Orbit-1 is a strong chaotic orbit, orbit-2 is a weak chaotic orbit, and orbit-3 and orbit-4 are ordered orbits. In Figure 1(c), errors of the Hamiltonian (11) are obtained when the explicit symplectic scheme S4 is used to calculate the four orbits in Figure 1(a). The four error graphs remain stably bounded, and increasing the initial value of r will reduce the errors due to a larger initial separation r corresponding to a longer orbital period. Figure 1(b) presents the energy comparisons among different integration methods. The error of the RK4 method increases with time. Thus, RK4 cannot be used for long-term integrators. The error of iS4 remains bounded at an order of O(10 −5 ). Compared with iS 4 and RK4, the newly constructed explicit symplectic scheme S2 can better preserve energy, because the error of S2 is O(10 −6 ) and is smaller than that of iS4. The error of S4 can reach an order of O(10 −9 ). The calculation time of S4 is less than that of iS4. We compare the time consumptions of 10 7 among different integrators in Table 1. With the same integration step ω, the explicit symplectic schemes are better than the implicit symplectic scheme in terms of calculation accuracy and time efficiency. The calculation speed of RK4 is the fastest of the fourth-order methods, but the error of RK4 gradually increases with the length of the calculation time.
Thanks to the excellent performance of the explicit scheme S4, the algorithm S4 is employed to study the dynamics of the charged test particles in the deformed Schwarzschild black hole. In what follows, we consider several methods for determining the orbital dynamical behavior.
LCE are expressed as  where v 0 ( )  and v ( )  w are the tangent vectors at time 0 and time ω, respectively. The LCE of the ordered orbit will tend to 0, whereas that of chaotic bounded orbit is stabilized to a positive value. The definition of the fast Lyapunov indicator (FLI) is provided by ) and v ( )  w are the lengths of the two tangent vectors at times 0 and ω, respectively. If the FLI of a bounded orbit grows linearly with time log10Ø, this orbit is regular. If the FLI of a bounded orbit grows exponentially, this orbit is chaotic.
where v 1 ( )  w and v 2 ( )  w are the two tangent vectors at time ω.
The SALI of a chaotic orbit will be close to 0, and that of an ordered orbit will remain at a positive constant. The definition of GALI is given as follows.
The GALI of s deviation vectors is determined through the evolution of s deviation vectors with 2 s 2N, which is initially linearly independent of the deviation vectors is the "norm" of the wedge product of the normalized deviation vectors v 1 and v 2 .
In the case of chaotic orbits, GALI s approaches zero as SALI does: w µ s s s s s s w --+ -+ +where σ 1 , σ 2 , ... σ s are the approximations of the s Lyapunov exponents.
Let us compare the efficiency of LCEs, FLIs, and SALIs in distinguishing between order and chaos. Figures 2(a)-(c) show the time evolution of the LCEs, FLIs, and SALIs of orbit-1 to orbit-4 in Figure 1(a). In Figure 2(a), the LCEs of the ordered orbits tend to 0 over time, whereas the LCE of the strong chaotic orbit-1 will not. As shown in Figure 2(b), the FLI of the strong chaotic orbit-1 grows rapidly, and the FLIs of the weak chaotic orbit-2 and ordered orbit-3 and orbit-4 are all less than 5. It is difficult to distinguish whether the orbits are weak chaotic or ordered. The LCE and FLI cannot distinguish weak chaotic orbits. In Figure 2(c), the SALIs of ordered orbit-3 and orbit-4 oscillate in the horizontal and approximately remain at two constants. The SALI of the weak chaotic orbit-2 appears unstable, and that of the strong chaotic orbit-1 approaches zero.   The comparison of these three chaotic indexes shows that SALI has the best effect in distinguishing chaotic orbits, which can effectively distinguish weak chaotic orbits, but the FLI and LCE cannot do so. Moreover, SALI will obtain the result at ω = 10 4 , whereas the FLI and LEC need ω > 10 5 .
In Figure 3(a), the SALI of ordered orbit-1 is a line that oscillates up and down, and is stable and greater than 10 −3 . The SALI of chaotic orbit-2 is a downward line, and, if there is enough calculation time, SALI will be equal to 0. The GALI and SALI of an ordered orbit are similar. GALI is a generalization of SALI, and its effectiveness in distinguishing chaos is similar to that of SALI, as shown in Figure 3(b). GALI2, GALI3, and GALI4 are the same as the initial values of orbit-2 in Figure 3(a). GALI2, GALI3, and GALI4 are all decreasing with time. The smaller the value of GALI, the stronger the effectiveness in terms of distinguishing chaos. Figure 3(b) shows the ability to distinguish chaos, where GALI4 > GALI3 > GALI2.
To know how the parameters Q, E, l, and ε affect the orbital dynamics, we use SALIs and Poincaré sections to study the effects of the four parameters on the dynamical transition from order to chaos. Figures 4(a1), (b1), (c1), and (d1) show the xcoordinate axis corresponding to the four parameters, and the ycoordinate axis corresponding to SALI. The value of SALI becomes smaller from top to bottom. Figures 4(a2), (b2), (c2), and (d2) plot phase-space structures associated with the four parameter values on the left-hand side. We find in Figure 4(a1) that if the deformation factor ε < 4, then strong chaos occurs. If ε is close to 5, then the SALI quickly rises to 10 −2 ; ε = 5 corresponds to a weak chaotic orbit. If ε > 5.2, then the values of SALI are always greater than 10 −3 , meaning that the orbits are no longer chaotic. In Figure 4(a2), the ε values of the three orbits are ε 1 = 3, ε 2 = 5, and ε 3 = 7. Obviously, the orbit with ε = 3 is composed of some disordered points, which constitute a typical strong chaotic orbit. For ε = 5, these points form a smooth curved solid line, which seems to be a standard weak chaotic orbit. For ε = 7, these points form a smooth solid line, which is an ordered orbit. ε = 5 is a boundary value between regularity and chaoticity.
We also find that the magnetic field parameter Q can explicitly affect the dynamics of the particle. In Figure 4(b1), if Q < 0 and the SALI is greater than 10 −3 , then a negative value of Q will not cause chaos. When Q > 4.5 * 10 −4 , the SALI becomes smaller and the orbit becomes unstable. The larger the value of Q, the smaller the SALI gets. In Figure 4(b2), we present the strong chaoticity for Q 1 = 9 * 10 −4 , the weak chaoticity orbit-2 for Q 2 = 5 * 10 −4 , and regularity for Q 3 = 3 * 10 −4 . The magnetic field around the black hole can be regarded as a perturbation to the charged test particles; if the value of Q is larger, then the perturbation to the charged particles becomes larger, and the orbit becomes unstable or chaos. However, there is no chaos when Q < 0. Now, we want to know how the angular momentum of the test particles affects the dynamic transition from order to chaos. The dynamic transition is clearly shown in Figure 4(c1). The xcoordinate axis is the angular momentum l, which increases from left to right, and the corresponding SALI shows an upward trend. When l = 4.2, the minimum SALI represents the maximum strength of chaos. If l > 4.45, the upward trend slows down, and at l > 4.6, the SALI becomes stable and the orbital dynamics will no longer be chaotic. The orbit is weakly chaotic at 4.45 < l < 4.6. The Poincaré sections of Figure 4(c2) are used to verify the SALIs for three values of the angular momentum, l 1 = 4.2, l 2 = 4.57, and l 3 = 4.7. The orbit for l 1 = 4.2 is composed of some disordered points. It is located in the outermost layer of the three orbits, and is a strong chaotic orbit. Although the orbit for l 3 = 4.7 is also composed of some disordered points, it is regular. These points on the orbit with l 2 = 4.57 seem to be a smooth closed ring, but in fact they are not. This is a separation orbit between chaos and order.
Finally, we study the dynamic transition from order to chaos as the energy of the particle varies from small to large. Figure 4(d1) shows that if 0.990 < E < 0.9945, the SALIs greater than 10 −3 are stable; if 0.9945 < E < 0.9975, the SALIs show a downward trend. At approximately E = 0.9975, the corresponding SALI value is the smallest, which represents the strongest chaos in the orbits. For 1 > E > 0.9975, the corresponding SALI values show an upward trend, but the orbits are still chaotic. Next, we analyze Figure 4 (d2), in which the orbit for E 1 = 0.9955 is composed of disordered points, and its chaoticity is stronger than that of Figures 4 (a2), (b2), and (c2). In fact, very strong chaos occurs for E 1 . The orbits for E 2 = 0.9954 and E 3 = 0.9930 are composed of regular points, which correspond to the regularity. Even if we increase the value of E 2 slowly, the Poincaré section directly changes from order to strong chaos. This shows the large effect of energy E on the dynamics of the test particles.
In sum, several main results can be concluded from Figure 4. Chaos is weakened as the positive deformation parameter and angular momentum increase. It is also weakened with an increase of the absolute value of the negative magnetic field parameter. However, chaos is easily induced when the energy, positive magnetic field parameter, and the magnitude of the negative deformation parameter increase from small to large. The results described by the SALI are consistent with those obtained from the method of the Poincaré section in the work of Yi & Wu (2020). The dependence of the chaotic properties on the parameters was explained by Yi & Wu (2020). That is to say, the chaotic properties are due to the inclusion of electromagnetic fields destroying the integrability of the deformed Schwarzschild spacetime. When the positive magnetic field parameter increases, the perturbation to the particles from the magnetic field increases, and chaos is naturally enhanced. An increase of the positive deformation parameter leads to increasing the gravity of the black hole to the particle and to weakening the strength of chaos. However, an increase of the magnitude of the negative deformation parameter makes the magnetic field force stronger, and chaos becomes stronger.

Conclusions
In this paper, we mainly focus on the construction and application of explicit symplectic schemes in a Hamiltonian system describing the motion of charged particles around a deformed Schwarzschild black hole with an external magnetic field.
First, we develop the time-transformed explicit symplectic schemes for the Kerr-type spacetimes introduced by Wu et al. (2021) for the deformed Schwarzschild spacetime. It is shown that the established algorithms exhibit good numerical performance in long-term integrations. Then, we also construct explicit symplectic methods for the variational equations of the Hamiltonian system. The obtained tangent vectors are applied to several chaos indicators, such as Lyapunov characteristic exponents, fast Lyapunov indicators, a smaller alignment index, and a generalized alignment index. It is found that the smaller alignment index and generalized alignment index are the fastest indicators for distinguishing between the regular and chaotic cases.
The smaller alignment index is employed to investigate how the parameters affect the dynamical transition from order to chaos. The main results are given here. For appropriate choices of initial conditions and other parameters, chaos is strengthened when the energy, positive magnetic field parameter, and the magnitude of the negative deformation parameter increase from small to large. However, chaos is weakened as the positive deformation parameter, angular momentum, and the absolute value of the negative magnetic field parameter increase.