Charged Particle Motions near Non-Schwarzschild Black Holes with External Magnetic Fields in Modified Theories of Gravity

Abstract

A small deformation to the Schwarzschild metric controlled by four free parameters could be referred to as a nonspinning black hole solution in alternative theories of gravity. Since 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 well long-term when stabilizing Hamiltonian errors, regardless of regular or chaotic orbits. One of the explicit symplectic integrators with the techniques of Poincaré 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 the magnitude of the magnetic parameter or any one of the negative deformation parameters increases. However, they are weakened when the angular momentum or any one of the positive deformation parameters increases.

1. Introduction

A Schwarzschild solution describing a nonrotating black hole and a Kerr solution describing a rotating black hole are two exact solutions of Einstein’s field equations of general relativity in a 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 an 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]. By 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 and 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 $\gamma$ metric (or Zipoy–Voorhees metric) describing a static and axially symmetric field [11,12] is also a parameterizing deviation from the Schwarzschild solution for $\gamma \ne 1$, it is an exact solution of Einstein’s equations in a vacuum.

In addition to the above-mentioned simply modified theories of gravity, such as scalar-tensor gravity, many other forms of modified theories of gravity can be found in the literature. Some examples are scalar-tensor theories, such as the Brans–Dicke theory [13,14], general scalar-tensor theories [15,16,17,18,19], Einstein-ther theories [20], Bimetric theories [21,22], tensor-vector-scalar theories [23,24], Einstein–Cartan–Sciama–Kibble theory [25,26], scalartensor-vector theory [27], $f\left(R\right)$ theories [28,29,30], $f\left(G\right)$ 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 should also see 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 and 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\left(R\right)$ 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 the innermost stable circular orbits are useful for understanding 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, such as magnetic fields, are 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 strong-field gravity features of general relativity, the motions of charged particles in the modified solutions with or without perturbations of weak external sources are naturally taken into account by some authors. The authors of [10] focused on the question of 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 $\gamma$ spacetime in the presence of an external magnetic field were studied in [66]. In fact, the $\gamma$ spacetime is nonintegrable and can allow for the onset of chaos if the external magnetic field is not included [67]. The authors of [68] gave some insights into the effect of one deformation parameter on the chaos of charged particles in the vicinity of a non-Schwarzschild black hole with an external magnetic field.

Numerical integration methods are vital to detecting the chaotic behavior of charged particles in the vicinity of the standard general-relativistic or modified black hole solutions with or without perturbations from weak external sources. They should have good stability and high precision so as to provide reliable results when detecting the chaotic behavior. The most appropriate long-term integration solvers for Hamiltonian systems constitute a class of symplectic integrators which respect the symplectic structures of Hamiltonian dynamics [69,70]. The motions of charged particles near the black holes with or without 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 group with explicit analytical solutions and another group with implicit solutions, explicit and implicit combined symplectic methods can be constructed [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 with or without 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 and 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 the question of 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 address that question in this paper. In addition, we mainly pay attention to the effects of the four free deformation parameters on the chaotic behavior. The present work is unlike the study in [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 presented in Section 5.

2. Deformed Schwarzschild Metric

In Schwarzschild coordinates $(t,r,\theta ,\phi )$, a Schwarzschild-like metric $d{s}^2=g_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }$ is written in [7,10] as

$$\begin{array}{ccc}\hfill d{s}^2& =& -f(1+h)d{t}^2+f^{-1}(1+h)d{r}^2+r^{2}d{\theta }^2+r^2{sin}^2\theta d{\phi }^2,\hfill \\ \hfill f& =& 1-\frac{2M}{r},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill h& =& k_0+\frac{k_{1}M}{r}+\frac{k_2{M}^2}{r^2}+\frac{k_3{M}^3}{r^3}.\hfill \end{array}$$

(2)

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 in the Schwarzschild metric are obtained from the modified field equations and the scalar field’s equation in dynamical theory. Clearly, Equation (1) with $h=0$ corresponds to the Schwarzschild metric. When $h\ne 0$, Equation (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:

$$\begin{array}{c}\hfill A_{\mu }=\frac{1}{2}{\delta }_{\mu }^{\phi }B{r}^2{sin}^2\theta ,\end{array}$$

(3)

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.

$$\begin{array}{c}\hfill H=\frac{1}{2m}{g}^{\mu \nu }(p_{\mu }-q{A}_{\mu })(p_v-q{A}_{\nu }).\end{array}$$

(4)

where $p_{\mu }$ is a generalized momentum, which is determined by

$$\begin{array}{c}\hfill {\dot{x}}^{\mu }=\frac{\partial H}{\partial p_{\mu }}=\frac{1}{m}{g}^{\mu \nu }(p_{\nu }-q{A}_{\nu }),\end{array}$$

(5)

equivalently,

$$\begin{array}{c}\hfill p_{\mu }=m{\dot{x}}^{\nu }{g}_{\mu \nu }+q{A}_{\mu }.\end{array}$$

(6)

The 4-velocity ${\dot{x}}^{\mu }$ is a derivative of the coordinate $x^{\mu }$ with respect to proper time $\tau$. As the Hamiltonian equations satisfy Equation (5) and

$$\begin{array}{c}\hfill {\dot{p}}_{\mu }=-\frac{\partial H}{\partial x^{\mu }};\end{array}$$

(7)

$p_t$ and $p_{\phi }$ are two constants of motion:

$$\begin{array}{ccc}\hfill p_t& =& m\dot{t}{g}_{tt}=-m\dot{t}f(1+h)=-E,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill p_{\phi }& =& m\dot{\phi }{g}_{\phi \phi }+q{A}_{\phi }=m{r}^2\dot{\phi }{sin}^2\theta +\frac{1}{2}qB{r}^2{sin}^2\theta =L.\hfill \end{array}$$

(9)

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\to tM$, $\tau \to \tau M$, $r\to rM$, $B\to B/M$, $E\to mE$, $p_r\to m{p}_r$, $L\to mML$, $p_{\theta }\to mM{p}_{\theta }$, $q\to mq$ and $H\to mH$. In this way, M and m in Equations (1)–(9) are taken as geometric units; $m=M=1$. The Hamiltonian (4) has two degrees of freedom $(r,\theta )$ in a four-dimensional phase space $(r,\theta ,p_r,p_{\theta })$, and can be rewritten as a dimensionless form:

$$\begin{array}{c}\hfill H=-\frac{E^2}{2f(1+h)}+\frac{1}{2r^2{sin}^2\theta }{(L-\frac{1}{2}Q{r}^2{sin}^2\theta )}^2+\frac{f{p}_r^2}{2(1+h)}+\frac{p_{\theta }^2}{2r^2},\end{array}$$

(10)

where $Q=Bq$.

Besides the two constants (8) and (9), the conserved Hamiltonian quantity

$$\begin{array}{c}\hfill H=-\frac{1}{2}\end{array}$$

(11)

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

3. Explicit Symplectic Integrations

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

3.1. 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 $\tau$. Since the Schwarzschild-like metric (1) can correspond to a Kerr-like metric via some coordinate transformation [12], the time-transformed explicit symplectic methods for the Kerr-type spacetimes proposed in [73] are guessed to be applicable to the system (10). The implementations of the algorithms are detailed below.

By extending the phase-space variables $(p_r,p_{\theta };r,\theta )$ of the Hamiltonian (10) to $(p_r,p_{\theta },p_0;$$r,\theta ,q_0)$, where $\tau$ is viewed as a new coordinate $q_0=\tau$ and its corresponding momentum is $p_0$ with $p_0=-H=1/2\ne p_t$, we have an extended phase-space Hamiltonian:

$$\begin{array}{c}\hfill J=H+p_0.\end{array}$$

(12)

It is clear that J is always identical to zero, $J=0$. By taking a time transformation

$$\begin{array}{ccc}\hfill d\tau & =& g\left(r\right)dw,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill g\left(r\right)& =& 1+h,\hfill \end{array}$$

(14)

we get a new time transformation Hamiltonian:

$$\begin{array}{ccc}\hfill \mathcal{H}=g\left(r\right)J& =& -\frac{E^2}{2f}+\frac{(1+h){(L-\frac{1}{2}Q{r}^2{sin}^2\theta )}^2}{2r^2{sin}^2\theta }+\frac{p_r^2}{2}-\frac{p_r^2}{r}\hfill \\ & & +\frac{(1+k_0)p_{\theta }^2}{2r^2}+\frac{k_1{p}_{\theta }^2}{2r^3}+\frac{k_2{p}_{\theta }^2}{2r^4}+\frac{k_3{p}_{\theta }^2}{2r^5}+p_{0}g\left(r\right).\hfill \end{array}$$

(15)

The Hamiltonian $\mathcal{H}$ has new coordinate time variable w and the phase-space variables $(p_r,p_{\theta },p_0;r,\theta ,q_0)$. As $J=0$, $\mathcal{H}=0$.

Similarly to the Hamiltonians of the Schwarzschild-type spacetimes in references [80,81,82], the time-transformed Hamiltonian $\mathcal{H}$ can be split in the following way:

$$\begin{array}{c}\hfill \mathcal{H}={\mathcal{H}}_1+{\mathcal{H}}_2+{\mathcal{H}}_3+{\mathcal{H}}_4+{\mathcal{H}}_5+{\mathcal{H}}_6+{\mathcal{H}}_7,\end{array}$$

(16)

where sub-Hamiltonians read

$$\begin{array}{ccc}\hfill {\mathcal{H}}_1& =& -\frac{E^2}{2f}+\frac{(1+h){(L-\frac{1}{2}Q{r}^2{sin}^2\theta )}^2}{2r^2{sin}^2\theta }+p_0(1+h),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\mathcal{H}}_2& =& \frac{p_r^2}{2},\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\mathcal{H}}_3& =& -\frac{p_r^2}{r},\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\mathcal{H}}_4& =& \frac{(1+k_0)p_{\theta }^2}{2r^2},\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\mathcal{H}}_5& =& \frac{k_1{p}_{\theta }^2}{2r^3},\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\mathcal{H}}_6& =& \frac{k_2{p}_{\theta }^2}{2r^4},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\mathcal{H}}_7& =& \frac{k_3{p}_{\theta }^2}{2r^5}.\hfill \end{array}$$

(23)

Each of the seven sub-Hamiltonians is analytically solvable, and its solutions are explicit functions of the new coordinate time w. $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$, $\mathcal{D}$, $\mathcal{E}$, $\mathcal{F}$ and $\mathcal{G}$ are differential operators, which correspond to ${\mathcal{H}}_1$, ${\mathcal{H}}_2$, ${\mathcal{H}}_3$, ${\mathcal{H}}_4$, ${\mathcal{H}}_5$, ${\mathcal{H}}_6$ and ${\mathcal{H}}_7$, respectively. These operators are written as

$$\begin{array}{ccc}\hfill \mathcal{A}& =& -\frac{\partial {\mathcal{H}}_1}{\partial r}\frac{\partial }{\partial p_r}-\frac{\partial {\mathcal{H}}_1}{\partial \theta }\frac{\partial }{\partial p_{\theta }}+\frac{\partial {\mathcal{H}}_1}{\partial p_0}\frac{\partial }{\partial q_0}\hfill \\ & =& f_1\frac{\partial }{\partial p_r}+f_2\frac{\partial }{\partial p_{\theta }}+(1+h)\frac{\partial }{\partial q_0},\hfill \\ \hfill f_1& =& \frac{k_1}{2r^2}+\frac{k_2}{r^3}+\frac{3k_3}{2r^4}-\frac{E^2}{r^2{(\frac{2}{r}-1)}^2}\hfill \\ \hfill & & +\frac{{(L-\frac{Q{r}^2{sin}^2\theta }{2})}^2(\frac{k_1}{r^2}+\frac{2k_2}{r^3}+\frac{3k_3}{r^4})}{2r^2{sin}^2\theta }\hfill \\ \hfill & & +(L-\frac{Q{r}^2{sin}^2\theta }{2})[Q+\frac{(L-\frac{Q{r}^2{sin}^2\theta }{2})}{r^2{sin}^2\theta }]\hfill \\ \hfill & & ·(\frac{k_0+1}{r}+\frac{k_1}{r^2}+\frac{k_2}{r^3}+\frac{k_3}{r^4})-p_0\frac{\partial h}{\partial r},\hfill \\ \hfill f_2& =& (L-\frac{Q{r}^2{sin}^2\theta }{2})[Q+\frac{(L-\frac{Q{r}^2{sin}^2\theta }{2})}{r^2{sin}^2\theta }]\hfill \\ \hfill & & ·(k_0+1+\frac{k_1}{r}+\frac{k_2}{r^2}+\frac{k_3}{r^3})cot\theta ,\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill \mathcal{B}& =& p_r\frac{\partial }{\partial p_r},\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill \mathcal{C}& =& -\frac{2}{r}{p}_r\frac{\partial }{\partial r}-\frac{p_r^2}{r^2}\frac{\partial }{\partial p_r},\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill \mathcal{D}& =& \frac{(1+k_0)p_{\theta }}{r^2}\frac{\partial }{\partial \theta }-\frac{(1+k_0)p_{\theta }^2}{r^3}\frac{\partial }{\partial p_r},\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill \mathcal{E}& =& \frac{k_1{p}_{\theta }}{r^3}\frac{\partial }{\partial \theta }-\frac{3}{2}\frac{k_1{p}_{\theta }^2}{r^4}\frac{\partial }{\partial p_r},\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill \mathcal{F}& =& \frac{k_2{p}_{\theta }}{r^4}\frac{\partial }{\partial \theta }-2\frac{k_1{p}_{\theta }^2}{r^5}\frac{\partial }{\partial p_r},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \mathcal{G}& =& \frac{k_3{p}_{\theta }}{r^5}\frac{\partial }{\partial \theta }-\frac{5}{2}\frac{k_3{p}_{\theta }^2}{r^6}\frac{\partial }{\partial p_r}.\hfill \end{array}$$

(30)

The solutions $\mathbf{z}={(r,\theta ,q_0,p_r,p_{\theta })}^T$ for the time-transformed Hamiltonian $\mathcal{H}$ advancing a new coordinate time step $\mathsf{\Delta }w=\sigma$ from the initial solutions $\mathbf{z}\left(0\right)={(r_0,{\theta }_0,q_{00},p_{r0},p_{\theta 0})}^T$ can be given by

$$\begin{array}{c}\hfill \mathbf{z}=S_2^{\mathcal{H}}\left(\sigma \right)\mathbf{z}\left(0\right),\end{array}$$

(31)

where $S_2^{\mathcal{H}}$ represents symmetric products of exponents of the seven operators and has the expressional form

$$\begin{array}{ccc}\hfill S_2^{\mathcal{H}}\left(\sigma \right)& =& e^{\frac{\sigma }{2}\mathcal{G}}\times e^{\frac{\sigma }{2}\mathcal{F}}\times e^{\frac{\sigma }{2}\mathcal{E}}\times e^{\frac{\sigma }{2}\mathcal{D}}\times e^{\frac{\sigma }{2}\mathcal{C}}\times e^{\frac{\sigma }{2}\mathcal{B}}\times e^{\sigma \mathcal{A}}\times e^{\frac{\sigma }{2}\mathcal{B}}\hfill \\ & & \times e^{\frac{\sigma }{2}\mathcal{C}}\times e^{\frac{\sigma }{2}\mathcal{D}}\times e^{\frac{\sigma }{2}\mathcal{E}}\times e^{\frac{\sigma }{2}\mathcal{F}}\times e^{\frac{\sigma }{2}\mathcal{G}}.\hfill \end{array}$$

(32)

Such symmetric products are a component of symplectic operators of 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 [87]:

$$\begin{array}{c}\hfill S_4^{\mathcal{H}}=S_2^{\mathcal{H}}\left(\gamma \sigma \right)\times S_2^{\mathcal{H}}\left(\delta \sigma \right)\times S_2^{\mathcal{H}}\left(\gamma \sigma \right),\end{array}$$

(33)

where $\gamma =1/(1-\sqrt[3]{2})$ and $\delta =1-2\gamma$.

3.2. Numerical Evaluations

Let us choose parameters $E=0.9965$, $L=4$, $Q=6\times {10}^{-4}$, $k_0={10}^{-3}$, $k_1={10}^{-2}$, $k_2={10}^{-1}$ and $k_3=1$. The initial conditions are $p_r=0$ and $\theta =\pi /2$. The initial value $r=15$ for Orbit 1, and $r=50$ for Orbit 2. The initial values $p_{\theta }>0$ for the two orbits are determined by Equation (11).

Given the time step $\sigma =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 1a. 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 $\sigma =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 $\sigma =4$ is used. Figure 1b describes the relationship between the proper time $\tau$ and the new coordinate time w when Orbit 1 is tested. Clearly, w is almost equal to $\tau$. This result coincides with the theoretical result $g\approx 1+k_0\approx 1$ when $r\gg 2$ and $k_0\approx 0$. Therefore, the time transformation function g in Equation (14) mainly plays an important role in implementing the desired separable form of the time-transformed Hamiltonian $\mathcal{H}$ rather than adaptive control to time steps.

4. 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 $\theta =\pi /2$ with $p_{\theta }>0$ in Figure 1c. 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 for the Hamiltonian error behavior is not related to 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\times {10}^{-4}$, $k_0={10}^{-4}$ and different values E are given. When $E=0.991$ in Figure 2a, the plotted seven orbits are ordered. As the energy increases, e.g., $E=0.9925$, three of the orbits are chaotic in Figure 2b. For $E=0.9975$ in Figure 2c, chaos is present almost elsewhere in the whole phase space. These results indicate an increase in the energy enhances the strength of chaos from the global phase-space structure. However, the chaotic properties are weakened as the particle’s angular momentum L increasing, 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 distinguish chaos from order. The largest Lyapunov exponent is defined in [88] by

$$\begin{array}{c}\hfill \lambda =\underset{w\to \infty }{lim}\frac{1}{w}ln\frac{d\left(w\right)}{d0},\end{array}$$

(34)

where $d0$ is the starting separation between the two nearby orbits and $d\left(w\right)$ 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

$$\begin{array}{c}\hfill \mathrm{FLI}={log}_{10}\frac{d\left(w\right)}{d0}.\end{array}$$

(35)

An exponential growth of FLI with time ${log}_{10}w$ means that the bounded orbit is chaotic, whereas a power law growth of FLI shows the bounded orbit is regular. When the integration time arrives at ${10}^6$, the FLIs in Figure 4a 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 4b. 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 distinguish 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 Figure 5a–c. When $Q=5\times {10}^{-4}$, all orbits are regular KAM tori in Figure 5a. Given $Q=8\times {10}^{-4}$ in Figure 5b, many tori are twisted and a few orbits can be chaotic. When $Q={10}^{-3}$ in Figure 5c, the number of chaotic orbits increases and the strength of chaos is enhanced. In other words, an increase in the positive magnetic parameter is helpful to induce the occurrence of chaos. How does a negative magnetic parameter affect the chaotic behavior as the magnitude of the negative magnetic parameter increases? The key to this question can be found in Figure 5d–f with $E=0.9975$. No chaos exists for $Q=-{10}^{-4}$ in Figure 5d. Three chaotic orbits are plotted for $Q=-8\times {10}^{-4}$ in Figure 5e. More orbits can be chaotic when $Q=-{10}^{-3}$ in Figure 5f. That is to say, the chaotic properties from the global phase-space structures are typically strengthened as the absolute value of the negative magnetic parameter increases. In short, chaos becomes stronger as the magnitude of the positive or negative magnetic parameter ($\left|Q\right|$) 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\times {10}^6$. All FLIs that are not less than 6 correspond to the onset of chaos, and those that are less than this value turn out to indicate the regularity of orbits. When $Q>8.5\times {10}^{-4}$ in Figure 6a or $Q<-7.5\times {10}^{-4}$ in Figure 6b, a dynamical transition from order to chaos occurs.

Now, let us focus on the dependence of chaos on the deformation parameters. Chaos becomes weaker when the deformation parameter $k_0$ is positive and increases in Figure 7a–c. However, it gets stronger when the deformation parameter $k_0$ is negative and its magnitude increases in Figure 7d–f. The effects of the deformation parameter $k_0$ on chaos described by the technique of Poincaré sections are consistent with those described by the method of FLIs in Figure 8. The effects of the other deformation parameters on chaos are shown through the methods of Poincaré sections and FLIs in Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14. They are similar to the effect of the deformation parameter $k_0$ on chaos. Precisely speaking, an increase in any one of the positive deformation parameters $k_1$, $k_2$ and $k_3$ weakens the chaotic properties, and an increase in any of the magnitudes of the negative deformation parameters $k_1$, $k_2$ and $k_3$ strengthens the chaotic properties. The result regarding the effects of the four deformation parameters on the chaotic properties is similar to the result of [68] 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 the energy E, magnetic parameter $\left|Q\right|$ or an absolute value of one of the negative deformation parameters ($|k_0|$, $|k_1|$, $|k_2|$ and $|k_3|$) increases, but weakened when the angular momentum L or any one of the 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 Equation (17) at the equatorial plane $\theta =\pi /2$ as

$$\begin{array}{ccc}\hfill {\mathcal{H}}_1& \approx & \frac{1}{2}[(1+k_0)(1-LQ)-E^2+\frac{k_2}{4}{Q}^2]-\frac{E^2}{r}+\frac{Q^2}{8}(1+k_0)r^2\hfill \\ & & +\frac{L^2}{2r^2}(1+k_0)+\frac{L^2{k}_1}{2r^3}+\frac{1-LQ}{2}(\frac{k_1}{r}+\frac{k_2}{r^2}+\frac{k_3}{r^3})+\frac{k_3}{2r}{Q}^2.\hfill \end{array}$$

(36)

The second term corresponds to the black hole gravity acting on 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’s 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\approx 1$, they have repulsive force effects on 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 in the energy E or the magnetic field $\left|Q\right|$ 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 in 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.

5. Conclusions

When a nonrotating compact object has spherical deformations, it has suffered from metric deformation perturbations. Such small deformation perturbations in 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 [83] should be similarly applicable to the deformed Schwarzschild black hole surrounded with an 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-step control. It was shown numerically that the obtained time-transformed explicit symplectic integrators perform well long-term in terms of 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 was used to 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, if the energy or the absolute value of the (positive or negative) magnetic parameter or any of the negative deformation parameters increases. However, it is weakened as the angular momentum or any one of the positive deformation parameters increases.