Applying explicit symplectic integrator to study chaos of charged particles around magnetized Kerr black hole

In a recent work of Wu, Wang, Sun and Liu, a second-order explicit symplectic integrator was proposed for the integrable Kerr spacetime geometry. It is still suited for simulating the nonintegrable dynamics of charged particles moving around the Kerr black hole embedded in an external magnetic field. Its successful construction is due to the contribution of a time transformation. The algorithm exhibits a good long-term numerical performance in stable Hamiltonian errors and computational efficiency. As its application, the dynamics of order and chaos of charged particles is surveyed. In some circumstances, an increase of the dragging effects of the spacetime seems to weaken the extent of chaos from the global phase-space structure on Poincare sections. However, an increase of the magnetic parameter strengthens the chaotic properties. On the other hand, fast Lyapunov indicators show that there is no universal rule for the dependence of the transition between different dynamical regimes on the black hole spin. The dragging effects of the spacetime do not always weaken the extent of chaos from a local point of view.


Introduction
The Kerr spacetime [1] as well as the Schwarzschild spacetime is integrable because the Carter constant [2] appears as a fourth constant of motion. However, a test charge motion in a magnetic field around the Kerr black hole or the Schwarzschild black hole is nonintegrable due to the absence of the fourth constant of motion. In some cases, the electromagnetic perturbation induces the onset of chaos and islands of regularity [3][4][5][6][7][8][9][10]. In a e-mail: wuxin_1134@sina.com particular, Takahashi & Koyama [11] found the connection between chaoticness of the motion and the black hole spin. That is, the chaotic properties are weakened with an increase of the black hole spin.
It is worth pointing out that the detection of the chaotical behavior needs reliable results on the trajectories. Usually, it needs long-time calculations, especially in the case of the computation of Lyapunov exponents. Thus, the adopted computational scheme is required to have good stability, high precision and small cost of computational time. Because these general relativistic curved spacetimes correspond to Hamiltonian systems with the symplectic nature, a low order symplectic integration scheme which respects the symplectic nature of Hamiltonian dynamics [12][13][14] is naturally regarded to as the most appropriate solver. The Hamiltonian systems for these curved spacetimes are not separable to the phase-space variables, therefore, no explicit symplectic integrators but implicit (or implicit and explicit combined) symplectic integrators [15][16][17][18][19][20][21][22] can be available in general. However, the implicit methods are more computationally demanding than the explicit ones at same order. To solve this problem,  split the Hamiltonians of Schwarzschild, Reissner-Nordström and Reissner-Nordström -(anti)-de Sitter spacetime geometries into four, five and six integrable separable parts, whose analytical solutions are explicit functions of proper time. In this way, explicit symplectic integrators were designed for the spacetime geometries. Unfortunately, the dragging effects of the spacetime by a rotating black hole make the Hamiltonian of Kerr spacetime have no desired splitting form similar to the Hamiltonian splitting form of Schwarzschild spacetime, and then lead to the difficulty in the construction of explicit symplectic integrators. This obstacle was successfully overcame by the introduction of time transformations of Mikkola [26] to the Hamiltonian of Kerr spacetime geometry. The obtained timetransformed Hamiltonian exists a desired splitting, and explicit symplectic integrators were established by Wu et al. [27]. At this point, a note is worthwhile on such time-transformed explicit symplectic integrators. They are constant time-steps in the new coordinate time and the algorithmic symplecticity can be maintained. Although they are variant time-steps in the proper coordinate time, the varying proper time-steps have only small differences and are approximately constant when the new integration time is not very long. In other words, the time transformation plays a main role in obtaining separable time-transformed Hamiltonian for the use of explicit symplectic integrators rather than adaptive time step control.
One of the main purposes in this paper is to explore whether the proposed explicit symplectic integrators for the Kerr spacetime geometry work well in the simulation of the chaotic motion of charged particles near the Kerr black hole immersed in a weak, asymptotically uniform external magnetic field. As another important purpose, we apply the established secondorder explicit symplectic method to check the result of Takahashi & Koyama [11] on the dragging effects weakening the chaotic properties, and to survey how a small change of the magnetic parameter affects a dynamical transition from order to chaos. Besides the Poincaré map method, the methods of Lyapunov exponents [28,29] and fast Lyapunov indicators [30][31][32] are employed to distinguish between regular and chaotic orbits and to quantify the dependence of transition from regular to chaotic dynamics on a certain parameter.
For the sake of our purposes, we introduce a Hamiltonian system for the description of charged particles moving around the Kerr black hole surrounded with an external magnetic field in Sect. 2. A second-order time-transformed explicit symplectic algorithm is constructed and used to quantify the transition from regular to chaotic dynamics by means of different methods finding chaos in Sect. 3. Finally, the main results are concluded in Sect. 4.

Kerr black hole with external magnetic field
Let us consider a particle with charge q moving around the Kerr black hole, surrounded by an external asymptotically uniform magnetic field with strength B. The particle motion is described by the following Hamiltonian formalism [33] Contravariant Kerr metric g µν has nonzero components A µ is an electromagnetic four-vector potential [33, 34] where ξ µ (t) = (1, 0, 0, 0) and ξ µ (φ) = (0, 0, 0, 1) are timelike and spacelike axial Killing vectors. Obviously, the four-vector potential has two nonzero covariant components The above covariant metric's components are 2ar sin 2 θ Σ = g φt , g φφ = (r 2 + a 2 + 2ra 2 Σ sin 2 θ) sin 2 θ. p µ represents a covariant generalized four-momentum, determined by the relatioṅ Note thatẋ µ is a derivative of coordinate x µ = (t, r, θ, φ) with respect to proper time τ , called as a 4-velocity. Eq.
(1) is rewritten as where F is a function of r and θ as follows: In Eq. (14), we take the speed of light c and the gravitational constant G as geometrized units, c = G = 1. The black hole's mass M is also one unit, M = 1. In practice, dimensionless operations are given to the related quantities via a series of scale transformations. That is, r → rM , t → tM , τ → τ M , a → aM , E → Em, p r → mp r , L → mM L, p θ → mM p θ , q → mq, B → B/M and H → m 2 H, where m denotes the particle's mass. Hereafter, we take β = qB. For any rotating black hole in general relativity, its spin angular momentum always satisfies |a| ≤ 1.
Because the timelike 4-velocity (5) always satisfies the relationẋ µẋ µ = −1, the Hamiltonian (14) itself is identical to a given constant Besides this value and E and L, other constants do not exist in the system. Thus, the Hamiltonian is nonintegrable.

Numerical investigations
Following the work of Wu et al. [27], we establish a second-order explicit symplectic integrator for a timetransformed Hamiltonian of the system (14). For comparison, a second-order implicit and explicit mixed symplectic method is applied to solve the Hamiltonian (14) in Sect. 3.1. Then, the explicit symplectic integrator is used to investigate the dependence of the regular and chaotic dynamics of the time-transformed Hamiltonian on a certain parameter by means of several methods finding chaos in Sect. 3.2.

Numerical integration scheme
Clearly, the Hamiltonian (14) is inseparable to the variables. Implicit symplectic methods can be applied to it without doubt. The first term F in Eq. (14) to obtain a time-transformed Hamiltonian to the Kerr geometry. The time-transformed Hamiltonian can be solved by explicit symplectic integrators. In the present problem, the time-transformed Hamiltonian is p 0 = −H = 1/2 is a momentum with respect to coordinate q 0 = τ . Thus, K is always identical to zero, K = 0. The time-transformed Hamiltonian is split into five parts where the sub-Hamiltonians read K 1 seems to be the same as that of Wu et al. [27] from the expressional form, but is unlike that of Wu et al. Their solutions correspond to operators K 1 , K 2 , K 3 , K 4 and K 5 . According to the idea of Wu et al., these operators can compose second-order explicit symplectic integrator for K: where h represents a new coordinate time-step. The use of fixed new coordinate time-step maintains the symplecticity of the algorithms for the time-transformed Hamiltonian K. However, the proper time τ takes varying step-sizes. Such symplectic integrators are adaptive time-steps [26,35,36]. However, the adaptive proper time step control is almost absent in the present problem because the time transformation function 1 < g ≤ 1 + 1/r 2 ≈ 1 for r 1. This indicates that the time transformation function mainly provides the desirable separable time-transformed Hamiltonian.
Based on the numerical results of Wu et al., the second-order explicit symplectic method S2 performs good computational efficiency and stabilizing error behavior, compared with the fourth-order explicit symplectic method S4. Thus, we focus on the application of S2 to the system (19). For comparison, the second-order implicit and explicit mixed symplectic integrator IE2 is applied to H. Let us take the time-step h = 1. Note that h = 1 denotes a proper time-step for IE2 and a new coordinate time-step for S2. The parameters are given by E = 0.995, L = 4.6, a = 0.5 and β = 0.001. The initial conditions are θ = π/2 and p r = 0. The initial separations r 0 = 11 and 75 are respectively given to Orbits 1 and 2. The initial value of p θ > 0 is determined by Eqs. (14) and (16). It is shown in Fig. 1(a) that Hamiltonian errors ∆ = K for the explicit symplectic algorithm S2 and ∆ = −H −1/2 for the second-order implicit and explicit mixed symplectic method IE2 remain stable and bounded. Both algorithms give no explicit error differences to the same orbit. The errors have the same order for the two orbits. If IE2 is applied to K, the results are almost the same as those for the application of IE2 to H. Although S2 is not explicitly better than IE2 in numerical accuracy, it is in computational efficiency, as claimed in the work of Wu et al.
It is worth pointing out that IE2 and S2 work in different time systems. The integration time for IE2 is the proper time τ = 10 8 . However, the integration time for S2 is the new coordinate time w = 10 8 . The new coordinate time w = 10 8 and its corresponding real proper time have no large differences. This result is shown clearly in Table 1. In fact, differences between the solutions of H and those of K are negligible before the integration time arrives at 10 5 . So are the differences between the proper times of K and the coordinate times of K. These facts confirm that the time transformation does not provide adaptive time-steps but a separable time-transformed Hamiltonian. Of course, these differences between the solutions of the two systems get slowly large as the integration time spans 10 5 and increases. In this case, the solutions of K at a given proper time are obtained with the aid of some interpolation method.
In fact, the tested orbits have different dynamical behaviors, as shown on Poincaré sections at the plane θ = π/2 and p θ > 0 in Fig. 1(b). The phase-space structure of Orbit 1 is a torus, regarded as the characteristic of an ordered orbit. However, the phase-space structure of Orbit 2 has many discrete points that are randomly filled in an area, regarded as the characteristic of a chaotic orbit. It can be seen from Fig. 1(a) that the performance of S2 (with IE2) is independent of the regularity and chaoticity of orbits. The phasespace structures in Fig. 1(b) are described by S2. They are also given by IE2.

Dynamical transition with parameters varying
In what follows, we use the algorithm S2 to trace the orbital dynamics of the system K. Meantime, methods of Poincaré sections, Lyapunov exponents and fast Lyapunov indicators are employed to find chaos.

Poincaré sections
It is worth pointing out that the chaoticity of Orbit 2 is due to the charged particle suffered from the electromagnetic field interaction. As claimed in the Introduction, the Kerr spacetime without the electromagnetic field interaction holds the Carter constant as the fourth integral of motion and therefore is integrable and nonchaotic. If the electromagnetic field is included in the Kerr spacetime, then the Carter constant is no longer present and the dynamics of charged particles becomes nonintegrable. This nonintegrability is an insufficient but necessary condition for the occurrence of chaos. When the electromagnetic field interaction acts as a relatively small perturbation, the gravitational force from the black hole is a dominant force and the dynamical system is nearly integrable. In this case, the Kerr spacetime with the electromagnetic field interaction can exhibit the similar dynamical features of the Kerr spacetime without the electromagnetic field interaction. That is to say, no chaos occurs and all orbits are regular Kolmogorov-Arnold-Moser (KAM) tori. This result is suitable for the magnetic parameter β = 4 × 10 −4 in Fig. 2(a). The tori in the perturbed case unlike those in the unperturbed case are twisted in the shapes; e.g., Orbit 1 corresponds to a peculiar torus, and Orbit 4 with initial separation r 0 = 110 yields a triangle torus. When β = 6 × 10 −4 in Fig. 2(b), the tori are still present, but are to a large degree different from those for β = 4 × 10 −4 in the shapes. There are typical differences of Orbits 1, 2 and 3 between Figs. 2 (a) and (b). In particular, Orbit 4 in Fig. 2(a) is evolved to a twisted figure-eight orbit with a hyperbolic point in Fig. 2(b) due to the electromagnetic field interaction. Such a hyperbolic point has a stable direction and another unstable direction, and easily induces chaos. As the electromagnetic field interaction increases and can generally match with the black hole gravitational force, chaos occurs. The hyperbolic point causes Orbit 4 to be chaotic for β = 7×10 −4 in Fig. 2(c). When the magnetic field strength is further enhanced, e.g., β = 8 × 10 −4 in Fig. 2(d), β = 1 × 10 −3 in Fig. 1(b) and β = 1.2 × 10 −3 in Fig. 2(e), chaos becomes stronger and stronger. Particular for β = 3 × 10 −3 in Fig. 2(f), chaos is existent almost everywhere in the whole phase space. In a word, it can be concluded clearly from Figs. 1(b) and 2 that an increase of the magnetic parameter is easier to induce chaos and enhances the strength of chaos from the global phase-space structure.
On the other hand, we are also interested in knowing what the dynamical transition with an increase of the black hole angular momentum a is. To answer this question, we consider some values of a. For a = 0 corresponding to the Schwarzschild case in Fig. 3(a), chaos densely exists almost everywhere in the whole phase space. This does not mean that any regular orbits are ruled out. In fact, regular Orbit 2 and black orbit with initial separation r 0 = 25 in Fig. 1(b) are still ordered in Fig. 3(a). Given a = 0.2, chaotic Orbits 1 and 3 in Fig. 3(a) become regular in Fig. 3(b), whereas ordered Orbit 2 in Fig. 3(a) becomes chaotic in Fig. 3(b). As the spin frequently increases [e.g., a = 0.5 in Fig. 1(b), a = 0.8 in Fig. 3(c) and a = 1 in Fig. 3(d)], the regular region seems to gradually increase and the chaotic region seems to decrease. In other words, an increase of a seems to weaken the extent of chaos from the global phase-space structure. This seems to support the result of [11] about the black hole spin weakening the chaotic properties. However, this does not mean that some individual orbits must be more chaotic as the dragging effects of the spacetime by a rotating black hole decrease. For example, Orbit 2 is regular for a=0, 1 but chaotic for a=0.2, 0.5, 0.8.

Fast Lyapunov indicators
Apart from the Poincaré map method, the technique of Lyapunov exponents is very efficient to detect the onset of chaos and provides a quantitative measure to chaos. An average exponential deviation of two nearby orbits in the system K is measured by the largest Lyapunov exponent where d(0) and d(w) represent the separations between the two nearby orbits at times 0 and w, respectively. Precisely speaking, d(w) is written as where (r 1 , θ 1 , p r1 , p θ1 ) is a set of phase-space variables of one orbit at time w and (r 2 , θ 2 , p r2 , p θ2 ) is a set of phase-space variables of its nearby orbit. Three points should be noted. (i) The initial distance d (0)  A positive Lyapunov exponent indicates the chaoticity of a bounded orbit, and zero Lyapunov exponent describes the regularity of a bounded orbit. This allows us to detect chaos from order. It is clearly seen from the Lyapunov exponents in Fig. 4(a) that Orbits 1 and 2 in Fig. 1 are ordered and chaotic, respectively. Orbit 3 in Fig. 2(d) seems to be regular till w = 10 7 because its Lyapunov exponent like Orbit 1's Lyapunov exponent decreases and does not tend to a stabilizing value. However, the Lyapunov exponent of Orbit 3 unlike that of Orbit 1 tends to a stabilizing positive value when the integration time spans from 10 7 to 10 8 ; therefore, Orbit 3 is chaotic. Because the Lyapunov exponent of Orbit 3 is smaller than that of Orbit 2, the chaoticity of Orbit 3 is weaker than that of Orbit 2. The result is consistent with that described by the method of Poincaré sections in Figs. 1(b) and 2(d).
In general, a long enough integration time is necessary to obtain a stabilizing value of Lyapunov exponent. In particular, the chaoticity of Orbit 3 cannot be distinguished from the regularity of Orbit 1 in terms of Lyapunov exponents until w = 10 8 . Compared with the method of Lyapunov exponents, the fast Lyapunov indicator (FLI) of Froeschlé & Lega [31] is based on computations of tangent vectors and is regarded as a quicker technique to separate chaotic orbits from regular ones. Wu et al. [32] modified the FLI with tangent vectors as the FLI of two nearby orbits: Here, d(0) = 10 −9 is admissible in the double-precision case, and a few but not many renormalizations to the distance d(w) are still necessary to avoid the saturation of orbits. Regular and chaotic orbits can be identified according to completely different time rates of growth of FLIs. An exponential increase of FLI of a bounded with time log 10 w indicates the chaoticity of the orbit; an algebraical increase of FLI of a bounded with time log 10 w indicates the regularity of the orbit. In other words, if the FLI of a bounded orbit is much larger than that of another bounded orbit for a given time, the former orbit is chaotic, but the latter orbit is regular.
On the basis of this point, the properties of orbits in Fig. 4(b) can be known easily. Without doubt, Orbit 2 is strongly chaotic because its FLI arrives at 12 when the integration time w = 6 × 10 5 , and is larger than 200 when w = 10 7 . Unlike the Lyapunov exponents, the FLIs are easy to distinguish between the regularity of Orbit 1 and the chaoticity of Orbit 3 when w = 10 7 because the FLI of Orbit 1 is smaller than 6, but the FLI of Orbit 3 is 13.5. As claimed in Ref.
[32], the FLI is a good tool to allow us to identify the transition between different dynamical regimes with a variation of a certain dynamical parameter. Letting a = 0.5 be fixed and β run from 0 to 1.5 × 10 −3 , we trace the dependence of the dynamical transition on the magnetic parameter. For a given value β, the integration lasts to the time w = 10 7 , but it ends if the FLI of some orbit is 20. There are three areas: regular area with FLIs smaller than 6, weak chaos area with FLIs no less than 6 but no more than 20 and strong chaos area with FLIs larger than 20. When β < 0.7 × 10 −3 , no chaos occurs in Figs. 5 (a) and (b). This is because the magnetic field force to charged particles is smaller than the black hole gravity to charged particles, as is mentioned above. β < 0.7 × 10 −3 is a regular area of the magnetic parameter. Chaos will get easier to a large degree as β increases. When β > 0.8 × 10 −3 , many values of β correspond to the onset of chaos and many other values indicate the presence of order in Fig. 5(a). In particular, chaos and order are present in the neighbourhood of β = 1.5 × 10 −3 . Unlike in Fig. 5(a), all values of β > 1.1 × 10 −3 can induce strong chaos in Fig. 5(b).
Of course, few values of the magnetic parameter correspond to weak chaos in the two panels. In short, the technique of FLIs and the method of Poincaré sections consistently show that an increase of the magnetic parameter leads to a strong probability for inducing the occurrence of chaos and strengthening the extent of chaos.
Figs. 6 (a)-(d) plot the dependence of FLIs on the black hole spin a when β = 0.001 is fixed. Strong chaos occurs as a → 0 but no chaos appears as a → 1 for the initial separation r 0 = 40 in Fig. 6(a). For the most values of a ∈ [0.25, 0.65], chaos is absent, too. The transition has some differences for the initial separation r 0 = 60 in Fig. 6(b). Chaos exists when a → 1. However, there is no chaos at the left end a = 0 and the right end a = 1 in Figs. 6 (c) and (d). These facts show that different combinations of initial conditions and other dynamical parameters affect the dependence of the dynamical transition on the black hole spin. Thus, no universal rule can be given to the relation between the dynamical transition and the black hole spin.
A notable point is that the aforementioned results are completely based on mathematical calculations. To justify the mathematical correctness of the obtained results, we employ an eighth-and ninth-order Runge-Kutta-Fehlberg integrator (RKF89) with adaptive step sizes to integrate the system (14) or (19). This algorithm can provide a machine precision to the Hamiltonians if roundoff errors are ignored. Of course, the higherorder method is more expensive in computations than the second-order explicit scheme S2. Consequently, the results obtained from RKF89 are consistent with those given by S2. It is very perfect and idealized to make a comparison of the obtained theoretical results with real observational results in the relativistic astrophysics. However, the known observational results have not sufficiently supported the study of chaotic behavior of black holes. There are two main reasons. On one hand, it costs a long enough time span to distinguish between the regular and chaotic two cases. The time for the determination of chaos is 10 7 M or 10 8 M from the numerical computations. Nevertheless, it is impossibly available from the observations. On the other hand, it costs massive data to detect chaos from order. This fact can be seen clearly from the numerical simulations. However, only a small amount of observational data are given. Because of these reasons, the study of chaos and islands of regularity of black holes in all the known publications such as [3-11, 23-25, 33 and 34] is restricted to the theoretically numerical results rather than the real observational results.

Conclusion
Because of the dragging effects of the spacetime by a rotating black hole, the Kerr geometry is complicated. In spite of this, it is integrable and nonchaotic. When an external magnetic field is included, the dynamics of charged particles becomes more complicated. In fact, it is nonintegrable and possibly chaotic. We find that the second-order explicit symplectic integrator for the integrable Kerr spacetime proposed in the work of Wu et al. [27] is still suited for the present nonintegrable problem. The successful application of the explicit symplectic integrator is owing to the contribution of the time transformation given in Eq. (17). The theoretical analysis and numerical results show that such a time transformation does not mainly bring adaptive proper time steps but a separable time-transformed Hamiltonian, which satisfies a need for the use of explicit symplectic integrator. The explicit symplectic algorithm performs good performance in numerical accuracy, stable error behaviour and computational cost regardless of whether tested orbits are regular or chaotic.
The explicit symplectic integrator is very suitable for studying the long-term qualitative evolution of orbits of charged particles moving around the magnetized Kerr black hole. The electromagnetic force interactions are mainly responsible for the nonintegrability and chaoticity of charged particle dynamics. In some circumstances, charged particle motions can be chaotic. An increase of the dragging effects of the spacetime seems to weaken the extent of chaos from the global phase-space structure on Poincaré sections. However, an increase of the magnetic parameter results in strengthening the chaotic properties. On the other hand, the fast Lyapunov indicators by scanning the black hole spin parameter show that no universal rule can be given to the dependence of the transition between different dynamical regimes on the black hole spin. This is because the transition is dependent on not only the spin parameter but also different combinations of initial conditions and other parameters. Unlike an increase of the black hole spin, that of the magnetic parameter may much easily induce chaos. Fig. 1 (a) Hamiltonian errors ∆ = K for the explicit symplectic algorithm S2 and ∆ = −H − 1/2 for the second-order implicit and explicit mixed symplectic method IE2. The parameters are E = 0.995, L = 4.6, a = 0.5 and β = 0.001. Orbits 1 and 2 have their initial separations r 0 = 11 and 75, respectively. The other initial conditions are θ = π/2, pr = 0 and p θ > 0 given by Eqs. (14) and (16). h = 1 is a proper time step for IE2 and a new coordinate time step for S2. Both algorithms make the errors have no secular drift and remain of the same order for Orbits 1 and 2. (b) Poincaré sections on the plane θ = π/2 and p θ > 0, given by S2. It is clear that Orbit 1 is regular, whereas Orbit 2 is chaotic.   Fig. 1(b), but for different values of the black hole's angular momentum a. Here, the magnetic parameter β = 1×10 −3 .
An increase of a seems to weaken the extent of chaos from the global phase space structure. The FLI of Orbit 2 reaches 12 when the integration time w = 6 × 10 5 , but is larger than 200 for the integration time w = 10 7 . Thus, this orbit is strongly chaotic. Weak chaos of Orbit 3 can be distinguished from regular Orbit 1 in terms of the growth of FLIs of the two orbits during the integration time w = 10 7 . The initial conditions θ = π/2 and pr = 0 are the same in the two panels, but the initial values r 0 and p θ > 0 are different. For the strong chaotic case, the integration ends when the FLI reaches 20, whereas the integration time w = 10 7 for the two other cases.