Charged particle motions near non-Schwarzschild black holes with external magnetic fields in modified theories of gravity

A small deformation controlled by four free parameters to the Schwarzschild metric could be referred to a nonspinning black hole solution in alternative theories of gravity. Because such a non-Schwarzschild metric can be changed into a Kerr-like black hole metric via a complex coordinate transformation, the recently proposed time-transformed explicit symplectic integrators for the Kerr type spacetimes are suitable for a Hamiltonian system describing the motion of charged particles around the non-Schwarzschild black hole surrounded with an external magnetic field. The obtained explicit symplectic methods are based on a time-transformed Hamiltonian split into seven parts, whose analytical solutions are explicit functions of new coordinate time. Numerical tests show that such explicit symplectic integrators for intermediate time-steps perform good long-term performance in stabilizing Hamiltonian errors regardless of regular or chaotic orbits. One of the explicit symplectic integrators with the techniques of Poincar\'{e} sections and fast Lyapunov indicators is applied to investigate the effects of the parameters including the four free deformation parameters on the orbital dynamical behavior. From the global phase-space structure, chaotic properties are typically strengthened under some circumstances as any one of the energy and the magnitudes of magnetic parameter and negative deformation parameters increases. However, they are weakened when each of the angular momentum and positive deformation parameters increases.


Introduction
A Schwarzschild solution describing a nonrotating black hole and a Kerr solution describing a rotating black hole are two exact solutions of the Einstein's field equations of general relativity in vacuum. According to the no-hair theorem, astrophysical (Kerr) black holes have their masses and spins as their unique characteristics. The theoretical prediction of the existence of black holes has been confirmed frequently by a wealth of observational evidence, such as X-ray binaries [1,2], detections of gravitational waves [3,4] and event-horizon-scale images of M87 [5,6].
Observational tests of strong-field gravity features cannot be based on a priori hypothesis about the correctness of general relativity. Instead, such tests must allow ansatz metric solutions to deviate from the general relativistic black hole scenarios predicted by the no-hair theorem. These metric solutions often come from perturbations of the usual Schwarzschild (or Kerr) black hole or exact solutions in alternative (or modified) theories of gravity. A small deformation to the Schwarzschild metric describing a nonspinning black hole (i.e., a modified Schwarzschild metric) [7] could be required to satisfy the modified field equations in dynamical Chern-Simons modified gravity [8,9]. Applying the Newman-Janis algorithm and a complex coordinate transformation, Johannsen and Psaltis [10] transformed such a Schwarzschild-like metric with several free deformation parameters into a Kerr-like metric, including a set of free deformation parameters as well as mass and spin. This Kerr-like metric, which is a parametric deformation of the Kerr solution and is not a vacuum solution, is regular everywhere outside of the event horizon. These metric deformations away from the Schwarzschild or Kerr metric by one or more parameters contain modified multipole structures. Although the γ metric (or Zipoy-Voorhees metric) describing a static and axially symmetric field [11,12] is also a parameterizing deviation from the Schwarzschild solution for γ = 1, it is an exact solution of Einstein's equations in vacuum.
In addition to the above-mentioned simply modified theories of gravity like scalar-tensor gravity, many other forms of modified theories of gravity can be found in the literature. Some examples are scalar-tensor theories including Brans-Dicke theory [13,14] and general scalar-tensor theories [15][16][17][18][19], Einstein-aether theories [20], Bimetric theories [21,22], tensor-vector-scalar theories [23,24], Einstein-Cartan-Sciama-Kibble theory [25,26], scalar-tensor-vector theory [27], f (R) theories [28][29][30], f (G) theory [31,32], Hořava-Lifschitz gravity [33][34][35] and higher dimensional theories of gravity [36][37][38]. Researchers and students in cosmology and gravitational physics see also review articles [29,30,[39][40][41] for more information on these modified gravity theories. Black-hole solutions in modified theories of gravity are generally unlike those in general relativity, and include many additional free parameters as well as the parameters predicted by the no-hair theorem in general relativity. Although a solution in a modified gravity model can be mathematically equivalent to a scalar field model, this mathematical correspondence does not always mean physical equivalence. The two corresponding solutions may have different physical behaviors. Corrections to the classical Einsteinian black hole entropy are necessary so as to constrain the viability of modified gravity theories in the study of Schwarzschild-de Sitter black holes by the use of the Noether charge method [42]. However, Not all black-hole solutions in modified theories of gravity must necessarily dissatisfy the Einstein field equations. For example, a stationary black-hole solution of the Brans-Dicke field equations must be that of the Einstein field equations [43]; this result is still present if no symmetries apart from stationarity are assumed [44]. The Kerr metric also remains a solution of certain f (R) theories [45].
A deep understanding of the relevant properties of the standard general relativistic black hole solutions and particle motions in the vicinity of the black holes is important to study accretion disk structure, gravitational lensing, cosmology and gravitational wave theory. Observational data from the vicinity of the circular photon orbits or the innermost stable circular orbits could be used as tests of the no-hair theorem. The properties of innermost stable circular orbits are useful to understand the energetic processes of a black hole. For this reason, radial effective potentials and (innermost) stable circular orbits of charged particles in electromagnetic fields surrounding a black hole have been extensively investigated in a large variety of papers (see, e.g., [46][47][48][49][50][51]). The motions of charged particles in the equatorial plane sound simple, but off-equatorial motions of charged particles in the magnetic fields become very complicated. In a stationary and axisymmetric black hole solution, there are three conserved quantities including the energy, angular momentum and rest mass of a charged particle. The fourth invariable quantity related to the azimuthal motion of the particle is destroyed in general when an electromagnetic field is included around the black hole. Thus, the particle motion in the spacetime background is not an integrable system. Chaos describing a dynamical system with sensitive dependence on initial conditions can occur in some circumstances. Various aspects of chaotic motions of charged particles around the standard general relativistic black holes perturbed by weak external sources like magnetic fields were discussed in many references (see, e.g., [52][53][54][55][56][57][58][59]).
Thanks to the importance of the deformed (or modified) black hole solutions in tests of strongfield gravity features of general relativity, the motions of charged particles in the modified solutions without or with perturbations of weak external sources are naturally taken into account by some authors. The authors of [10] focused on a question how the radii of the innermost stable circular orbits and circular photon orbits vary with increasing values of the spin and deviation parameters in a Kerr-like metric of a rapidly rotating black hole. They demonstrated that their Kerr-like metric is suitable for strong-field tests of the no-hair theorem in the electromagnetic spectrum. Charged particle motions around non-Schwarzschild or non-Kerr black hole immersed in an external uniform magnetic field were considered in [60][61][62]. The influence of a magnetic field on the radial motion of a charged test particle around a black hole surrounded with an external magnetic field in Hořava-Lifshitz gravity was investigated in [63][64][65]. The radial motions of charged particles in the γ spacetime in the presence of an external magnetic field was studied in [66]. In fact, the γ spacetime is nonintegrable and can allow for the onset of chaos if the external magnetic field is not included in [67]. The authors of [68] gave some insight into the effect of one deformation parameter on chaos of charged particles in the vicinity of non-Schwarzschild black hole with external magnetic field.
Numerical integration methods are vital to detect the chaotical behavior of charged particles in the vicinity of the standard general-relativistic or modified black hole solutions without or with perturbations from weak external sources. They should have good stability and high precision so as to provide reliable results to detecting the chaotical behavior. The most appropriate long-term integration solvers for Hamiltonian systems are a class of symplectic integrators which respect the symplectic structures of Hamiltonian dynamics [69,70]. The motions of charged particles near the black holes without or with weak external sources can be described by Hamiltonian systems, and thus allow for the applicability of symplectic methods. If the Hamiltonian systems are split into two parts, explicit symplectic integration algorithms are not available in general. However, implicit symplectic integrators, such as the implicit midpoint rule [71,72] and implicit Gauss-Legendre Runge-Kutta symplectic schemes [54,73,74], are always suitable for their applications to these Hamiltonian systems that do not need any separable forms. When the Hamiltonians are separated into one part with explicit analytical solutions and another part with implicit solutions, explicit and implicit combined symplectic methods could be constructed in [75][76][77][78][79]. The implicit algorithms are more computationally demanding than the explicit ones in general, therefore, the explicit symplectic integrations should be developed as much as possible. Recently, the authors of [80][81][82] successfully constructed the explicit symplectic integrators for the Schwarzschild type black holes without or with external magnetic fields by splitting the corresponding Hamiltonians into several parts having analytical solutions as explicit functions of proper time. More recently, the time-transformed explicit symplectic integrators were designed for the Kerr family spacetimes [83][84][85].
The idea for constructing the time-transformed explicit symplectic integrators as well as the explicit symplectic integrators introduced in [80][81][82][83] allows for the applicability of many standard general-relativistic or modified black hole solutions with or without perturbations of weak external sources. In spite of this, there is no universal rule on how to construct explicit symplectic integrators for Hamiltonians corresponding to the spacetimes. Specific Hamiltonian problems have different separations, or different choices of time-transformed Hamiltonians and their splitting forms. As is claimed above, the non-Schwarzschild metric with four free deformation parameters could produce a Kerr-like metric through a complex coordinate transformation [10]. Now, there is a question whether the time-transformed explicit symplectic integrators for the Kerr type spacetimes [83] are applicable to such a deformed non-Schwarzschild black hole immersed in an external magnetic field. We plan to address this question in this paper. In addition, we mainly pay attention to the effects of the four free deformation parameters on the chaotical behavior. The present work is unlike Ref. [68] in which one deformation parameter is added to the non-Schwarzschild metric and no explicit symplectic integrators are considered.
The remainder of this paper is organized as follows. A metric deformation to the Schwarzschild spacetime is introduced in Section 2. Time-transformed explicit symplectic integrators are described in Section 3. Orbital dynamical properties are discussed in Section 4. Finally, the main results are concluded in Section 5.

Deformed Schwarzschild metric
In Schwarzschild coordinates (t, r, θ , ϕ), a Schwarzschild-like metric ds 2 = g αβ dx α dx β is written in [7,10] as Here, M denotes a mass of the black hole. The speed of light c and the gravitational constant G are taken as geometric units, c = G = 1. Deformation function h is a perturbation to the Schwarzschild metric, where k 0 , k 1 , k 2 and k 3 are deformation parameters. It comes from modified multipole structures related to spherical deformations of the star. When the action through algebraic, quadratic curvature invariants coupled to scalar fields is modified, such small deformations to the Schwarzschild metric are obtained from the modified field equations and the scalar field's equation in the dynamical theory. Clearly, Eq. (1) with h = 0 corresponds to the Schwarzschild metric. When h = 0, Eq. (1) looks like the Schwarzschild metric but can be transformed into a Kerr-like black-hole metric by the Newman-Janis algorithm [86] and a complex coordinate transformation [10].
Suppose the black hole is immersed in an external electromagnetic field with a four-vector potential where B is a constant strength of the uniform magnetic field. The motion of a test particle with mass m and charge q is described in the following Hamiltonian Here, p µ is a generalized momentum, which is determined bẏ equivalently, The 4-velocityẋ µ is a derivative of the coordinate x µ with respect to proper time τ. Because the Hamiltonian equations satisfy Eq. (5) anḋ p t and p ϕ are two constants of motion: E is an energy of the particle, and L is an angular momentum of the particle.
For simplicity, dimensionless operations are given to the related quantities as follows: t → tM, In this way, M and m in Eqs. (1)-(9) are taken as geometric units, m = M = 1. The Hamiltonian (4) has two degrees of freedom (r, θ ) in a four-dimensional phase space (r, θ , p r , p θ ), and can be rewritten as a dimensionless form where Q = Bq.
Besides the two constants (8) and (9), the conserved Hamiltonian quantity is a third constant of the system (10). The third constant of motion exists due to the invariance of the 4-velocity or the rest mass of the particle in the time-like spacetime (1). Given Q = 0, the system (10) holds a fourth constant of motion and therefore is integrable and nonchaotic. When Q = 0, the system (10) has no fourth constant and then becomes nonintegrable. In this case, analytical solutions cannot be given to the system (10), but numerical solutions can.

Explicit symplectic integrations
First, time-transformed explicit symplectic methods for the system (10) is introduced. Then, their performance is numerically evaluated.

Design of algorithms
As is claimed above, the metric (1) seems to be the Schwarzschild metric but the system (10) is not suitable for the application of the explicit symplectic methods suggested in [80][81][82] because the Hamiltonian (10) is unlike the Hamiltonians of the Schwarzschild type spacetimes (including the Reissner-Nordström metric, the Reissner-Nordström-(anti)-de Sitter solution and these spacetimes perturbed by external magnetic fields), which can be separated into several parts having analytical solutions as explicit functions of proper time τ. Since the Schwarzschild-like metric (1) can correspond to a Kerr-like metric via some coordinate transformation [12], the timetransformed explicit symplectic methods for the Kerr type spacetimes proposed in [73] are guessed to be applicable to the system (10). The implementation of the algorithms is detailed below.
Extending the phase-space variables (p r , p θ ; r, θ ) of the Hamiltonian (10) to (p r , p θ , p 0 ; r, θ , q 0 ), where τ is viewed as a new coordinate q 0 = τ and its corresponding momentum is p 0 with p 0 = −H = 1/2 = p t , we have an extended phase-space Hamiltonian It is clear that J is always identical to zero, J = 0. Taking a time transformation we get a new time transformation Hamiltonian The Hamiltonian H has new coordinate time variable w and the phase-space variables (p r , p θ , p 0 ; r, θ , q 0 ).
Similar to the Hamiltonians of the Schwarzschild type spacetimes in Refs. [80][81][82], the timetransformed Hamiltonian H can be split in the following form where sub-Hamiltonians read Each of the seven sub-Hamiltonians is analytically solvable and its solutions are explicit functions of the new coordinate time w. A , B, C , D, E , F and G are differential operators, which correspond to H 1 , H 2 , H 3 , H 4 , H 5 , H 6 and H 7 , respectively. These operators are written as The solutions z = (r, θ , q 0 , p r , p θ ) T for the time-transformed Hamiltonian H advancing a new coordinate time step ∆w = σ from the initial solutions z(0) = (r 0 , θ 0 , q 00 , p r0 , p θ 0 ) T can be given by where S H 2 is symmetric products of exponents of the seven operators and has the expressional form Such symmetric products are a component of symplectic operators at second order. The symplectic method S 2 is an extension to the works of [83][84][85] regarding the time-transformed explicit symplectic methods for the Kerr spacetimes. Of course, such symmetric products of order 2 easily yield a fourth-order construction of Yoshida [86] where γ = 1/(1 − 3 √ 2) and δ = 1 − 2γ.
Given the time step σ = 1, the errors of the Hamiltonian J for the second-order method S 2 and the fourth-order method S 4 solving Orbit 1 have no secular drifts. The errors are three orders of magnitude smaller for S 4 than for S 2 before the integration time w = 10 7 , as shown in Figure 1(a). With the integration spanning this time and tending to w = 10 8 , the errors still remain bounded for S 2 , but exhibit long-term growths for S 4 . The secular drifts of the Hamiltonian errors for S 4 are due to roundoff errors. When the number of integration steps is small, the truncation errors are more important than the roundoff errors. As the integration is long enough, the roundoff errors are dominant errors and cause the Hamiltonian errors to grow with time. However, such error drifts for S 4 lose when a larger time step σ = 4 is adopted. If Orbit 1 is replaced with Orbit 2, the Hamiltonian errors for each of the two methods are not explicitly altered.
In what follows, S 4 with the time step σ = 4 is used. Figure 1(b) describes the relationship between the proper time τ and the new coordinate time w when Orbit 1 is tested. Clearly, w is almost equal to τ. This result coincides with the theoretical result g ≈ 1 + k 0 ≈ 1 when r ≫ 2 and k 0 ≈ 0. Therefore, the time transformation function g in Eq. (14) mainly plays an important role in implementing the desired separable form of the time-transformed Hamiltonian H rather than adaptive control to time steps.

Regular and chaotic dynamics of orbits
The regularity of Orbit 1 and the chaoticity of Orbit 2 are clearly shown through the Poincaré section at the plane θ = π/2 with p θ > 0 in Figure 1(c). The phase-space of Orbit 1 is a Kolmogorov-Arnold-Moser (KAM) torus, which belongs to the characteristic of a regular quasi-periodic orbit. For Orbit 2, many discrete points are densely, randomly filled with an area and are regarded as the characteristic of a chaotic orbit. The Hamiltonian errors for S 4 acting on Orbit 1 are approximately same as those for S 4 acting on Orbit 2. This fact indicates that the algorithmic performance in the Hamiltonian error behavior is regardless of the regularity or chaoticity of orbits. Now, we continue to use the technique of Poincaré section to trace the orbital dynamical evolution. The parameters are the same as those in Figure 1 but Q = 8 × 10 −4 , k 0 = 10 −4 , and different values E are given. When E = 0.991 in Figure 2(a), the plotted seven orbits are ordered. As the energy increases, e.g., E = 0.9925, three of the orbits are chaotic in Figure 2(b). For E = 0.9975 in Figure 2(c), chaos is present almost elsewhere in the whole phase space. These results indicate an increase of the energy leading to enhancing the strength of chaos from the global phase-space structure. However, the chaotic properties are weakened as the particle angular momenta L increase, as shown in Figure 3.
Besides the technique of Poincaré section, Lyapunov exponents for measuring an exponential rate of the separation between two nearby orbits with time are often used to detect chaos from order. The largest Lyapunov exponent is defined in [88] by where d0 is the starting separation between the two nearby orbits and d(w) is the distance between the two nearby orbits at time w. However, it takes long enough time to obtain stabilizing values of the Lyapunov exponents. Instead, a fast Lyapunov indicator (FLI), as a quicker method to distinguish between the ordered and chaotic two cases, is often used. It comes from a slightly modified version of the largest Lyapunov exponent, and is calculated in [88] by An exponential growth of FLI with time log 10 w means the chaoticity of an bounded orbit, whereas a power law growth of FLI shows the regularity of an bounded orbit. When the integration time arrives at 10 6 , the FLIs in Figure 4(a) can clearly identify the regular and chaotic properties of three energies corresponding to the orbits with the initial separation r = 15 in Figure 2. The regular and chaotic properties of three angular momenta corresponding to the orbits with the initial separation r = 70 in Figure 3 are also described the FLIs in Figure 4(b). Clearly, the angular momentum L = 4.4 corresponds to the regularity, whereas the angular momenta L = 3.85 and L = 4 correspond to chaos. Chaos is stronger for L = 3.85 than for L = 4. As far as the Poincaré sections and FLIs are concerned, they are two popular methods to detect chaos from order. The technique of Poincaré sections can clearly, intuitively describe the global phase-space structure, but is mainly applicable to conservative systems with two degrees of freedom or four dimensional phase spaces. The method of FLIs is suitable for any dimensions.
Taking the parameters L = 4, k 0 = 10 −4 , k 1 = 10 −2 , k 2 = 10 −1 and k 3 = 1, we employ the technique of Poincaré sections to plot the global phase-space structures with E = 0.9915 for three positive values of the magnetic parameter Q in Figures 5 (a)-(c). When Q = 5 × 10 −4 , all orbits are regular KAM tori in Figure 5(a). Given Q = 8 × 10 −4 in Figure 5(b), many tori are twisted and a few orbits can be chaotic. When Q = 10 −3 in Figure 5 Figure 5(e). More orbits can be chaotic when Q = −10 −3 in Figure 5(f). That is to say, the chaotic properties from the global phase-space structures are typically strengthened with an increase of the absolute value of the negative magnetic parameter. In short, chaos becomes stronger as the magnitude of the positive or negative magnetic parameter (|Q|) varies from small to large. This result is also supported by the FLIs in Figure 6. Here, the FLI for a given value of Q is obtained after the integration time w = 2 × 10 6 . All FLIs that are not less than 6 correspond to the onset of chaos, while those that are less than this value turn out to indicate the regularity of orbits. When Q > 8.5 × 10 −4 in Figure  6(a) or Q < −7.5 × 10 −4 in Figure 6 [69] for describing the effect of deformation parameter k 3 on the chaotic properties.
The above demonstrations clearly show how small changes of these parameters affect the dynamical transitions from order to chaos. The main result is that chaos in the global phase space is strengthened as any one of the energy E, magnetic parameter |Q|, absolute values of the negative deformation parameters (|k 0 |, |k 1 |, |k 2 | and |k 3 |) increases, but weakened when any one of the angular momentum L and positive deformation parameters k 0 , k 1 , k 2 and k 3 increases. Here, an interpretation is given to the result. Expanding 1/ f in the Taylor series, we rewrite Eq. (17) at the equatorial plane θ = π/2 as The second term corresponds to the black hole gravity to the particles. The third term yields an attractive force from a contribution of the magnetic field regardless of whether Q > 0 or Q < 0. The fourth term provides an inertial centrifugal force due to the particle angular momentum. The fifth, sixth and seventh terms come from coupled interactions among the metric deformation perturbations, angular momentum and magnetic field. For 1 − LQ ≈ 1, they have repulsive force effects to the charged particles when k 1 > 0, k 2 > 0 and k 3 > 0, but attractive force effects when k 1 < 0, k 2 < 0 and k 3 < 0. A small increase of the energy E or the magnetic field |Q| means enhancing the attractive force effects and therefore the motions of particles can become more chaotic in some circumstances. As the angular momentum L increases, the repulsive force effects are strengthened and chaos is weakened. With a minor increase of relatively small positive deformation parameter k 0 , the magnetic field attractive force and the centrifugal force will increase, but the centrifugal force has a larger increase than the magnetic field force for the parameters chosen in Figure 7. This leads to weakening the strength of chaos. However, as the absolute value |k 0 | with k 0 < 0 increases, the centrifugal force has a larger decrease than the magnetic field force and chaos becomes stronger. Increases of the other positive deformation parameters k 1 , k 2 and k 3 cause the repulsive forces to increase, and chaos to get weaker. However, the attractive force effects are enhanced and chaos gets stronger as the magnitudes of negative deformation parameters k 1 , k 2 and k 3 increase.

Conclusions
When a nonrotating compact object has spherical deformations, it is suffered from metric deformation perturbations. Such small deformation perturbations to the Schwarzschild metric could be regarded as a nonrotating black hole solution departure from the standard Schwarzschild spacetime in modified theories of gravity. The non-Schwarzschild spacetime with four free deformation parameters is integrable. However, the dynamics of charged particles moving around the Schwarzschild-like black hole is nonintegrable when the inclusion of an external asymptotically uniform magnetic field destroys the fourth invariable quantity related to the azimuthal motion of the particles.
Although the deformation perturbation metric looks like the Schwarzschild metric, it can be changed into a Kerr-like black hole metric via some appropriate coordinate transformation. Therefore, the time-transformed explicit symplectic integrators for the Kerr type spacetimes introduced in [84] should be similarly applicable to the deformation perturbation Schwarzschild black hole surrounded with external magnetic field. In fact, we can design explicit symplectic methods for a time-transformed Hamiltonian, which is split into seven parts with analytical solutions as explicit functions of new coordinate time. A main role for the time transformation function is the implementation of such desired separable form of the time-transformed Hamiltonian rather than that of adaptive time-steps control. It is shown numerically that the obtained time-transformed explicit symplectic integrators perform good long-term stable error behavior regardless of regular or chaotic orbits when intermediate time-steps are chosen.
One of the obtained time-transformed explicit symplectic integrators combined with the techniques of Poincaré sections and FLIs is used to well show how small changes of the parameters affect the dynamical transitions from order to chaos. Chaos in the global phase space can be strengthened under some circumstances as any one of the energy and the absolute values of the (positive or negative) magnetic parameter and negative deformation parameters increases. However, it is weakened as any one of the angular momentum and positive deformation parameters increases.