Applicability of the 0–1 test for chaos in magnetized Kerr–Newman spacetimes

The dynamics of electrically neutral or charged particles around a magnetized Kerr–Newman black hole immersed in an external electromagnetic field can be described by a dimensionless Hamiltonian system. This Hamiltonian is given an appropriate time transformation, which allows for construction of explicit symplectic integrators. Selecting one of the integrators with good accuracy, long-term stabilized Hamiltonian error behavior and less computational cost, we employ the 0–1 binary test correlation method to distinguish between regular and chaotic dynamics of electrically neutral or charged particles. The correlation method is almost the same as the techniques of Poincaré map and fast Lyapunov indicators in identifying the regular and chaotic two cases. It can well describe the dependence of the transition from regularity to chaos on varying one or two dynamical parameters. From a statistical viewpoint, chaos occurs easily under some circumstances with an increase of the external magnetic field strength and the particle electric charge and energy or a decrease of the black hole spin and the particle angular momentum. A small change of the black hole electric charge does not very sensitively affect the dynamics of neutral particles. With the black hole electric charge increasing, positively charged particles do not easily yield chaotic motions, but negatively charged particles do. On the other hand, the effect of a small change of the black hole magnetic charge on the dynamical transition from order to chaos has no universal rule.


Introduction
It is known that there are strong magnetic fields near supermassive black holes at the center of galaxies [1].When the external magnetic fields are close to the upper limit a e-mail: wuxin_1134@sina.com(corresponding author) 10 19 M /M Gauss, where M and M are the masses of the Sun and black holes, they make the black holes strongly magnetized.They can not only influence the motion of charged particles drastically but also change the background black hole geometry.They can affect the dynamics of accretion, and are helpful to explain the observed emission from the black holes.Some examples of magnetized black hole spacetimes are Schwarzschild, Reissner-Nordström (RN), and Kerr-Newman black holes embedded in the Melvin universe [2,3].The vector potential for the magnetic field was given by Wald [4].The authors of [5] provided a complicated expression of the magnetic vector potential at the magnetized Kerr-Newman solution.These magnetized black hole solutions can be obtained by applying a Harrison magnetizing transformation to the non-magnetized black hole solutions.The magnetized Kerr-Newman solution is an axisymmetric solution of Einstein-Maxwell equations because the rotation axis is aligned to the magnetic field.The magnetic fields in these types of spacetime geometries destroy the existence of a fourth constant of motion (analogous to Carter's constant in the Kerr metric), and thus lead to the non-integrability of the magnetized black hole spacetimes.This fact was confirmed by finding chaotic motion of electrically neutral particles around the magnetized Schwarzschild black hole [6,7].There are imprints of a chaotic pattern in the shadow images and the lensing of the magnetized Schwarzschild black hole [8].Self-similar fractal structures in the shadow of magnetized Kerr black hole show the chaotic scattering motion of photons near the magnetized Kerr black hole [9].The authors of [10] studied chaos of electrically charged and neutral particles in the magnetized RN spacetime.
When the strength of an external magnetic field surrounding a black hole is much smaller than the upper limit, it does not influence the background black hole geometry.However, it can strongly affect charged particle motion in the black hole vicinity for the particle with large charge to mass ratio.Under some circumstances, the charged particle motion can be chaotic, as several authors claimed in [11][12][13][14][15][16][17][18][19].The chaotic scatter of charged particles in the combined gravitational and magnetic fields leads to charged particle acceleration along the magnetic field lines [20].Combining the chaotic motion with regular motion in the final part of the motion can also lead to the possibility of extreme accelerations [21].
The detection of the orbital chaotical behavior requires the use of fast and accurate numerical tools that provide reliable results.Because the motion of electrically charged and neutral particles in many curved spacetimes can be described in terms of Hamiltonian systems, the most appropriate solvers in the case of long-term integration are symplectic schemes which preserve the symplectic nature of Hamiltonian dynamics.Such Hamiltonian systems have no separation of variables or two explicitly integrable splitting terms in general.In this case, the construction and application of explicit symplectic methods are somewhat difficult.Implicit symplectic methods are always applicable to these problems.They are, e.g., the Gauss-Legendre Runge-Kutta (GLRK) method [12,22], and slimplectic integrator [23].These completely implicit symplectic methods have poor computational efficiency compared with combined implicit and explicit symplectic methods [24][25][26][27][28][29].In order to further reduce the cost in the computational time, our group have recently developed a class of explicit symplectic integrators for the Hamiltonian systems of many curved spacetimes.Such algorithmic constructions are based on the Hamiltonians or time-transformed versions of the Hamiltonians split into more than two explicitly integrable pieces [30][31][32][33][34].
The detection of the orbital chaotical behavior also relies on efficient methods.The chaos detection methods are the Poincaré map, frequency map analysis of Laskar [35], Lyapunov exponents [36], fast Lyapunov indicators (FLIs) [37], smaller alignment index (SALI) [38], generalized alignment index (GALI) [39], and 0-1 test [40].The Poincaré map is mainly applied to clearly distinguish between ordered and chaotic motion in conservative Hamiltonian systems with four-dimensional phase space.The other methods are appropriate for Hamiltonian systems with any dimensions.The SALI, GALI and FLI can much rapidly determine the chaotic or regular nature of orbits than the maximum Lyapunov exponents.Unlike the maximal Lyapunov exponent whose computation requires phase space reconstruction, the 0-1 test method can be directly applied to the time series data.The latter method is superior to the former method in the case of underlying quasiperiodic dynamics with amounts of measurement noisy data.The 0-1 test correlation method of Gottwald and Melbourne [41] significantly increases the sensitivity for characterizing the chaotic or regular nature of orbits compared with the 0-1 test regression method of Gottwald and Melbourne [40].The 0-1 test has been applied successfully for the distinction between regular and chaotic dynamics in some Hamiltonian systems [42][43][44][45][46][47].
In the present work, we focus on the application of the 0-1 test correlation method combined an explicit symplectic integrator to a dynamical system in curved spacetime.This system relates to the motion of electrically neutral and charged particles around the magnetized Kerr-Newman black hole immersed in an external electromagnetic field [5].For the sake of our purpose, we simply introduce the dynamical model in Sect. 2. Explicit symplectic integrator are designed in Sect.3. In Sect.4, we check the validity of the correlation method for testing order and chaos.The 0-1 test is used to investigate effects of the parameters on chaos in Sect. 5. Finally, the main results are concluded in Sect.6.

Magnetized Kerr-Newman black holes
A magnetized Kerr-Newman black hole spacetime is introduced.Then, its corresponding Hamiltonian formulation is provided.

Black hole metric
Suppose that a rotating black hole has mass m, rotating angular momentum a, electric charge q and magnetic charge p, and is immersed in an external magnetic field B. In the Boyer-Lindquist coordinates x α = (t, r, θ, φ), such a rotating black hole can be described in terms of a Kerr-Newmanlike metric [5]: where the nonzero covariant metric components are The related notations are given as follows: Functions H (1) , H (2) , H (3) and H (4) are expressed as where q is given by The geometric unit systems C = G = 1 are given to the speed of light C and the gravitational constant G.
There is an external electromagnetic field, which corresponds to the vector potential Functions Φ 0 , Φ 3 and ω have complicated expressions.Φ 0 reads as where functions ω 1 , ω 2 , ω 3 and ω 4 are ( Φ 0 is written as where The expression of Φ 3 is where In a word, the four-vector potential (12) has two nonzero covariant components If B = 0 and p = 0, the metric (1) represents the Kerr-Newman black hole.For the case of B = 0, the magnetic field plays an important role in the metric (1).In other words, the Kerr-Newman spacetime geometry is typically changed by the magnetic field.Only when q = −am B, is the metric (1) believed to be asymptotic to the static Melvin metric [48], as was shown in [5].

Hamiltonian system
The motion of a test particle with electric charge e and mass m 0 around the magnetized Kerr-Newman black hole (1) is governed by the following Hamiltonian where function H 1 reads as The nonzero contravariant components of the metric (1) are The covariant 4-momentum p α = ( p t , p r , p θ , p φ ) is defined by where the 4-velocity ẋα = dx α /dτ is a derivative of the coordinate x α with respect to the proper time τ , and the Lagrangian system is This 4-momentum is also determined by a set of canonical equations for the Hamiltonian ( 29) Another set of canonical equations for the Hamiltonian (29) are This shows that p t and p φ are constants of motion, which correspond to the energy of the particle E = −p t and the angular momentum of the particle L = p φ .They are In this case, the Hamiltonian system ( 29) is axially symmetric, and also stationary.In other words, the existence of the two constants of motion is attributed to these background symmetries.A third constant of motion is the conserved Hamiltonian for a time-like geodesic orbit because the 4-velocity satisfies the relation g αβ ẋα ẋβ = −1.In fact, the third constant is due to the rest mass of the particle.There is no fourth constant of motion in the Hamiltonian system (29) due to the strong external magnetic field changing the Kerr-Newman spacetime structure.This fact was similarly shown by finding chaos of neutral particles in the magnetized Schwarzschild spacetime [6,7] and magnetized Reissner-Nordström spacetime [10].
For simplicity, the Hamiltonian system (29) uses dimensionless operations.Namely, scale transformations are given to the related variables and parameters: e, and B → B/m.In this case, m 0 and m appear as mass factors in the above equations and can be eliminated.

Explicit symplectic integrators
The Hamiltonian system (29) cannot be split into more than two integrable terms with analytical solutions as explicit functions of the proper time.However, its appropriate time transformation Hamiltonian can, as the time transformation Hamiltonians for the Kerr spacetime and some other curved spacetimes mentioned in the papers [10,33,34] can.In this way, explicit symplectic integrators are easily available for such time transformation Hamiltonians.In what follows, we introduce how to implement explicit symplectic methods for the Hamiltonian (29).

Construction of numerical schemes
The phase space (r, θ, p r , p θ ) of Hamiltonian ( 29) is extended to a new phase space (r, θ, p r , p θ , p 0 , q 0 ), where q 0 = τ and p 0 = −H = 0.5 are new coordinate and momentum, respectively.The Hamiltonian in the expanded phase space is It is clear that H * = 0 for any proper time.Taking the time transformation dτ = g(r, θ)dw (39) with the time transformation function we have a time transformation Hamiltonian Obviously, J = 0 for any new time w.
The time transformation Hamiltonian has five splitting pieces where these sub-Hamiltonians are It is easy to check that each of the sub-Hamiltonians is analytically solvable, and its analytical solutions are all explicit functions of the new time w.The five sub-Hamiltonians correspond to solvers T 1 , T 2 , T 3 , T 4 , and T 5 in sequence.Setting h as a time step, we obtain two first-order operators Symmetrically combining the two operators, we have a second-order explicit symplectic method for solving the Hamiltonian (42) Symmetrically combining three second-order explicit symplectic methods, we obtain a fourth-order method [49] where γ = 1/(2 − 3 √ 2) and δ = 1 − 2γ .The Hamiltonian (42) is suitable for the application of an optimized fourthorder partitioned Runge-Kutta explicit symplectic integrator [50,51] where the time coefficients are listed in [50,51] as follows: In short, S 2 , S 4 and P RK 6 4 are the explicit symplectic integrators we have designed for the time transformation Hamiltonian (42).They adopt fixed time steps for the new time w, but seem to employ variable ones for the proper time τ .

Numerical tests
The time step h = 1 is used in the new time.The parameters are given by the magnetic field strength B = 4 × 10 −4 , black hole electric charge q = 0.6, black hole rotating angular momentum a = 0.1, particle electric charge e = 0.5, particle angular momentum L = 4.6, and particle energy E = 0.998.The initial conditions are the angle θ = π/2, p r = 0, and p θ > 0 determined by Eq. (37).The initial separations are r = 25 for Orbit 1, and r = 85 for Orbit 2.
In Fig. 1a, b, the second-order method S 2 and fourth-order method S 4 solving the two tested orbits exhibit good longterm performance.There is no secular drift with the new time w in the errors of the time transformation Hamiltonian (41).In this sense, the two methods have one of the advantages which symplectic methods own.However, the method P RK 6 4 gives secular drift to the Hamiltonian errors due to the rapid accumulation of roundoff errors.
When r 2, g(r, θ) ≈ 1 in Eq. ( 40).Thus, τ ≈ w in Eq. ( 39).This result is also confirmed numerically.In fact, the new time w uses a fixed step size, and the proper time τ almost adopts the same fixed step size.In other words, the proper time step is not adaptive if the new time step is fixed.
The phase space structures of Orbits 1 and 2 are plotted in Fig. 1c regarding the Poincaré map in the two-dimensional r − p r plane.In fact, the structures consist of many intersections of the particles' trajectories with the surface of section θ = π/2 and p θ > 0. The distribution of the points from an orbit in the Poincaré map can show the regularity or chaoticity of the motion.Because Orbit 1 is a closed curve, it is regular nonchaotic.Orbit 2 has many randomly distributed points and therefore is irregular or chaotic.Seen from Fig. 1c, the accuracies in the energy conservation of the algorithms S 2 and S 4 for the regular orbit 1 are almost the same as those for the chaotic orbit 2. This fact shows that the algorithms S 2 and S 4 well conserve the energy irrespective of whether the tested orbit is regular or chaotic.
Although the algorithm S 4 is three orders of magnitude smaller than the method S 2 in the energy accuracy, the former is more computationally demanding than the latter.In addition, the method S 2 reaches an order of 10 −6 as an acceptable accuracy.Thus, the method S 2 is used in the later computations

Correlation method on the 0-1 test for chaos
Apart from the technique of Poincaré map, there are several methods finding chaos, such as the 0-1 test.Considering that the mean square displacements have different growth rates for nonchaotic and chaotic dynamics, Gottwald and Melbourne [40] used a regression method to establish a binary test yielding 0 or 1 for distinguishing between the two cases.The method for testing chaos is superior to tangent space methods for computing maximal Lyapunov exponents [43].Gottwald and Melbourne proposed computing the correlation coefficient as the correlation of the mean square displacement with time [41].The correlation method significantly increases its sensitivity for detecting weakly chaotic from regular dynamics, compared with the regression method.

Description of the correlation method
Let ψ( j) for j = 1, 2, 3, . . ., n 1 be a set of observation data from a deterministic dynamical system.For a given parameter c ∈ (0, π), two translation variables are defined as where n 1 = 1, 2, . . ., N 1 .Then, the mean square displacement of the translation variables P c and Q c is given by In principle, n 2 ≤ ñ N 2 is required, but ñ ≤ N 2 /10 is acceptable in practical computations.The time step is h = 1.The two orbits have the dynamical parameters, which involve the magnetic field strength B = 4×10 −4 , the black hole's electric charge q = 0.6, magnetic charge p = 0.05 and spin a = 0.1, and the particle's electric charge e = 0.5, energy E = 0.998 and angu-lar momentum L = 4.6.They have the initial conditions p r = 0 and θ = π/2.The initial separations are r = 25 for Orbit 1 and r = 85 for Orbit 2. The initial values p θ > 0 are given by Eq. ( 11).The secondorder method S 2 and fourth-order method S 4 show no secular drift in the errors.c Poincaré sections in the r − p r plane.Orbit 1 is regular, while Orbit 2 is chaotic The mean square displacement M c (n 2 ) is modified as where V osc (c, n 2 ) is an oscillatory term with an explicit expression If M c (n 2 ) has a linear term, then D c (n 2 ) remains the linear term of M c (n 2 ), but eliminates the oscillatory term of M c (n 2 ).In fact, the asymptotic growth rate of D c (n 2 ) with time number n 2 is the same as that of M c (n 2 ).Suppose that the vectors ξ = (1, 2, ..., ñ) and Δ c = (D c (1), D c (2), ..., D c ( ñ)).Their mean values are Variances on the two vectors are A covariance of ξ and Δ c is Now, the correlation coefficient of the two vectors is computed by From a statistical viewpoint, the mean square displacement D c (n 2 ) is bounded for any number of observation data n 2 associated with time in the case of nonchaotic dynamics, whereas grows linearly in the related time number n 2 for chaotic dynamics.Thus, K c = 0 for regular dynamics and K c = 1 for chaotic dynamics.In other words, K c = 0 indicates the presence of regular dynamics, while K c = 1 indicates the presence of chaotic dynamics.K c is the 0-1 chaos test from the correlation method, which was introduced by Gottwald and Melbourne [41].There is a minor difference between the description of K c in Eq. (62) and that in [41].That is, n 1 and n 2 were taken as n, and N 1 and N 2 were written as N in [41].In fact, a larger value N 2 can cause K c to much accurately determine regular or chaotic dynamics.However, we do not need such many points (P c , Q c ) when the visual of P c and Q c is used to underly regular or chaotic dynamics.This is why N 1 < N 2 is required.A note is worthwhile on the choice of c.In fact, the value of K c depends on the choice of c.Most values of c ∈ (0, π) can make the values of K c yield consistent results in signifying nonchaotic or chaotic dynamics.However, isolated values of c lead to the occurrence of resonances and then the values of K c do not signify the actual dynamics.In this case, the authors of [41] suggested that the median of the computed values of K c rather than the mean of the computed values of K c should be used against outliers associated with resonances.A theoretical justification of the 0-1 test for chaos was given by Gottwald and Melbourne [42].

Validity of the correlation method for testing order and chaos
In our present computations, an orbit is integrated until the integration time w = 10 7 .For computations of data in a plot on P c and Q c , we take c = 1.7,N 1 = 4500, and ψ( j) = r j , where j corresponds to the time w j = j × w/N 1 .When each value of K c is computed, we take N 2 = 10,000, and give 100 values to c, i.e. c i = i × π/100 (i = 1, 2, . . ., 100).For each orbit, the 100 values of c correspond to those of K c .The median of the 100 values of K c is taken as the value of K c for this orbit.In this way, the values of K c for Orbits 1 and 2 in Fig. 1 can be obtained as follows: K c = 0.00804 ≈ 0 for Orbits 1 in Fig. 2a, and K c = 0.98801 ≈ 1 for Orbits 2 in Fig. 2d.The two values 0 and 1 of K c clearly characterize the regularity of Orbits 1 and the chaoticity of Orbits 2, respectively.The results can also be shown from visual plots on the evolution behavior of P c and Q c .The points (P c , Q c ) form a bounded region such as a ring surface for underlying regular dynamics in Fig. 2a, but behave as an asymptotic unbounded Brownian motion for underlying chaotic dynamics in Fig. 2d.The regular and chaotic dynamics of Orbits 1 and 2 described by the 0-1 binary test and the visuals of P c and Q c is consistent with that given by the method of Poincaré map in Fig. 1c.It is also confirmed by the largest Lyapunov exponents in Fig. 2c, fast Lyapunov indicators (FLIs) in Fig. 2e, power spectra in Fig. 2c, f.There are tangent space methods and direct methods for computing the maximal Lyapunov exponents.The direct methods for computing the maximal Lyapunov exponents directly measure the divergence of nearby trajectories with time [52].They were modified in [53,54] as where d(0) and d(w) represent the proper distances between the nearby trajectories at time w and the initial time, respectively.A bounded orbit is chaotic if λ > 0, while it is regular if λ ≤ 0. The FLIs, which use completely different time rates of the lengths of deviation vectors to distinguish between the ordered and chaotic two cases, are from modified versions of the Lyapunov exponents.The tangent space methods for the FLIs were introduced by Froeschlé [37] and Froeschlé and Lega [55].The direct method for computing the FLI was given in [54] by An exponential growth of the FLI with time log 10 w describes chaotic dynamics of a bounded orbit, but an algebraical growth shows regular dynamics.The FLIs are a faster and more sensitive tool to detect chaos from order than the Lyapunov exponents.The power spectrum has discrete peaks at the harmonics and subharmonics for a periodic or quasiperiodic system, while has a broadband nature for a chaotic signal.These characteristics can be used to distinguish periodic, quasiperiodic and chaotic motions.Besides regular Orbit 1 as a KAM (Kolmogorov-Arnold-Moser) torus in Fig. 1, some regular trajectories consisting of several distinct islands are used to evaluate the performance of P c , Q c and K c in Fig. 3.The appearance of such islands in each trajectory means the existence of resonance.The values of K c are less than 0.12, and show the regular dynamics of the three resonant orbits.For the three-island orbit and seven-island orbit produced through the Poincaré surface of section in Fig. 3a, c, their corresponding visuals of P c and Q c in Fig. 3d, f are similar to the visual for the single-torus orbit 1 in Fig. 2a.The visual of P c and Q c for the five-island orbit in Fig. 3b is a torus surface in Fig. 3e.These results show that the visual shapes of P c and Q c for such resonant orbits have no typically different characteristics, compared with those for these single-torus orbits.

Effects of the parameters on chaos
In this section, the impacts of the parameters on a transition from regularity to chaos are considered.At first, we detect chaos from order by tracing the visual shapes of (P c , Q c ) and the values of K c .Then, we find chaos through the dependence of K c on varying one or two parameters.

Finding chaos by the evolution behavior of P c , Q c and K c
Let us take the parameters E = 0.998, L = 4.5, a = p = 0.1, q = 0.5 and e = 0.The initial separation is r = 60.
Seen the visuals of P c and Q c from Fig. 4, the magnetic field strengths B = 1.5×10 −4 , 3×10 −4 and 7×10 −4 correspond to order, chaos and chaos, respectively.Their corresponding values of K c = 0.13989, 0.99389 and 0.99598 also show these dynamical behaviors.In the case of the particle without electric charge, chaos can be present because the magnetic field acts as a gravitational effect destroying the integrability of the Kerr-Newman spacetime.In fact, the onset of chaos of neutral particles or photons moving around the magnetized Schwarzschild, Reissner-Nordström, and Kerr black holes was reported in [6][7][8][9][10].The FLIs with the values of show that the extent of chaos increases with the magnetic field strength increasing.Now, we consider three distinct energies E = 0.996, 0.9967 and 0.998 in Fig. 5, where the other parameters are B = 5 × 10 −4 , L = 4.7, a = p = 0.1, q = 0.5 and e = 0, and the initial separation is r = 65.The Poincaré surfaces of section describe that the three energies correspond to a regular KAM torus, weak chaos and strong chaos.These regular and chaotic dynamical behaviors can also be seen from the visuals of (P c , Q c ).The values of K c = 0.1922, 0.98038, 0.99668 tell us a path for the transition from regularity to weak and strong chaos.
Figure 6 plots the section surfaces, visuals of (P c , Q c ) and values of K c for three different particle angular momenta L = 4.5, 5.3 and 6.2.The other parameters are B = 4.5 × 10 −4 , E = 0.998, a = p = 0.1, q = 0.5 and e = 0, and the initial separation is r = 50.When the three values of the angular momentum are changed from small to large, strongly chaotic dynamics, weakly chaotic dynamics and regular dynamics occur, as shown via the section surfaces and values of K c .The extent of chaos is weakened with the particle angular momentum increasing under some circumstances.This result is also suitable for an increase of the black hole spin a in  To investigate the effects of the black hole electric charges q on the dynamical behavior, we consider three cases of the particle charge e =, >, or < 0. For the case of e = 0 in Fig. 9, the visuals of (P c , Q c ) and values of K c show that the black hole electric charges q = 0.2, 0.5 and 0.8 correspond to the regularity in Fig. 9a-c, where B = 2 × 10 −4 is given.However, the black hole charges q = 0.2, 0.5 and 0.8 correspond to the onset of chaos when B = 7 × 10 −4 in Fig. 9d-f.If e = 0.5 in Fig. 10, the values of K c are 0.99661 for q = 0.2, 0.99192 for q = 0.6 and 0.01599 for q = 0.9, and seem to show that an increase of the black hole charges q weakens the extent of chaos.When e = −0.8 in Fig. 11a-c, the black hole charges q = 0.2, 0.45 and 0.8 correspond to the values of K c = 0.00467, 0.97603 and 0.99178.This seems to show that chaos is induced easily as the black hole charge q increases.If q = 0.8 is fixed in Fig. 11c-e, chaos seems to be weakened as the magnitude of the particle negative charge q decreases.
Figure 12 plots the visuals of (P c , Q c ) and values of K c , which show how a variation of the black hole magnetic charge p affects the dynamics for the three cases of the particle charge e =, >, or < 0. In the three cases, the transition from regularity to chaos is not very sensitive dependence on the   Figure 13 plots the dependence of K c on varying one parameter.Here, we take N 2 = 10,000.For a given parameter, the value of K c is not obtained until the integration time reaches 10 7 .The values of 0 ≤ K c < 0.5 indicate the regularity, whereas those of 0.5 ≤ K c ≤ 1 correspond to the chaoticity.When the magnetic field B spans 3.6 × 10 −4 in Fig. 13a, the transition from regularity to chaos occurs.As the black hole spin a is larger than 0.69, chaos disappears in Fig. 13b.The values of e ≥ 0.21 induce chaos almost anywhere in Fig. 13c.These facts show that chaos is easily induced as the magnetic field strength and particle positive charge increase or the black hole spin decreases.On the other hand, chaos is always present for any one of the black hole magnetic charges p when a smaller angular momentum L = 4.0 is taken in Fig. 13d.However, only several values of p admit chaos when a larger angular momentum L = 6.5 is used.
We can also use the values of K c to exhibit the distribution of order and chaos in the space of two parameters.In Fig. 14, the particles are electrically neutral, and the two parameters are q and one of B, E and L. Clearly, the regular and chaotic dynamics of neutral particles is not very sensitive dependence on the variation of q.Under some circumstances, the increase of B and E or the decrease of L can easily induce chaos from the statistical computations.When the particles have negative charges in Fig. 15a, chaos occurs easily as q increases.However, the chaotic motion of positive charge particles occurs easily as q decreases in Fig. 15b.The dependence of chaos on the black hole charge in the magnetized Kerr-Newman spacetime is similar to that in the magnetized RN spacetime [10].There is no universal rule on the contribution of a variation of p to the transition from order to chaos in Fig. 15c, d.
As a point to be illustrated, the parameters corresponding to regular and chaotic dynamics described by the values of K c in Figs. 13, 14 and 15 are completely consistent with those shown by the FLIs.

Summary
An external magnetic field which surrounds the Kerr-Newman black hole is so strong that it changes the spacetime geometry.The strong magnetic field in the spacetime acts as a gravitational effect to electrically neutral particles.It leads to the absence of the Carter constant.As a result, such a magnetized Kerr-Newman spacetime is nonintegrable.When test particles moving around the black hole are electrically charged, they are subject to the interaction of an external electromagnetic field as well as the magnetized Kerr-Newman spacetime.In this case, the dynamical model of electrically charged particles becomes more complicated.
The motion of electrically neutral or charged particles around the magnetized Kerr-Newman black hole immersed in the external electromagnetic field is described in terms of a dimensionless Hamiltonian system.The Hamiltonian has several physical parameters, including the external magnetic field strength, the black hole's rotating angular momentum, electric charge and magnetic charge, and the particle's charge, angular momentum and energy.It has no separation of variables, and cannot be split into several explicitly integrable terms, either.However, an appropriate time transformation to the Hamiltonian allows for the existence of five explicitly integrable splitting parts.In this way, explicit symplectic integrators are easily available.It is shown numerically that the second-and fourth-order explicit symplectic integrators perform high accuracy in the conserved Hamiltonian quantity and show no secular drift in the Hamiltonian error.
Considering computational efficiency, we select the secondorder explicit symplectic method as an integration tool to provide some insight into the dynamics of electrically neutral or charged particles.In addition, the 0-1 binary test correlation method is mainly used to distinguish between regular and chaotic dynamics.The correlation coefficient K c = 0 corresponds to the regular dynamics, while the correlation coefficient K c = 1 represents the chaotic dynamics.A bounded region like a ring surface, as the visual of two translation variables P c and Q c , can characterize the regularity.If the visual of two translation variables P c and Q c behave as an asymptotic unbounded Brownian motion, it can characterize the chaoticity.The correlation method is effective to detect chaos from order, as the other techniques of Poincaré map, power spectrum, Lyapunov exponents and FLIs are.It is more sensitive to find chaos than the Lyapunov exponents.
From a statistical viewpoint, the transition from regularity to chaos occurs easily under some circumstances as both the external magnetic field strength and the particle charge and energy increase.Chaos is also easily induced when the black hole spin and the particle angular momentum decrease.The effects of the black hole charge on the regular and chaotic dynamics are dependent on the sign of the particle charge.The variation of the black hole charge does not very sensitively affect the dynamics of neutral particles.With the black hole charge increasing, the motions of positively charged particles are not easily chaotic, whereas those of negatively charged particles are.On the other hand, no universal rule can be given to the effect of a small change of the black hole magnetic charge on the dynamical transition from order to chaos.

Fig. 1 a
Fig. 1 a Errors of the Hamiltonian J for the three symplectic integrators acting on Orbit 1. b Same as a but Orbit 1 replaced with Orbit 2.The time step is h = 1.The two orbits have the dynamical parameters, which involve the magnetic field strength B = 4×10 −4 , the black hole's electric charge q = 0.6, magnetic charge p = 0.05 and spin a = 0.1, and the particle's electric charge e = 0.5, energy E = 0.998 and angu-

Fig. 2 a
Fig. 2 a, d The visuals of two translation variables P c and Q c with the values of K c for Orbits 1 and 2. b The largest Lyapunov exponents λ for Orbits 1 and 2. e Fast Lyapunov indicators (FLIs) for Orbits 1 and 2. c, f Power spectra for Orbits 1 and 2. The 0-1 binary test correla-

Fig. 4 a
Fig. 4 a-c The visuals of P c and Q c with the values of K c for three magnetic field strengths B = 1.5 × 10 −4 , B = 3 × 10 −4 , and B = 7 × 10 −4 .The other parameters are E = 0.998, L = 4.5,

Fig. 8 aFig. 9 Fig. 7 .
Fig. 8 a-c The visuals of P c and Q c with the values of K c for three values of the particle charges e = 0.2, e = 0.5, and e = 0.9.The other parameters are B = 4 × 10 −4 , E = 0.998, L = 5.7 and

Fig. 12 a 6 Fig. 13
Fig. 12 a-d The visuals of P c and Q c with the values of K c and FLIs for three values of the black hole magnetic charges p.The other parameters are e = 0, B = 4.5 × 10 −4 , E = 0.998, L = 6.5, a = 0.1, and the initial separation is r = 45.e-h Same as a-d, but the other

Fig. 14 45 Fig. 15
Fig.14 The dependence of K c on two varying parameters for the case of electrically neutral particles.The magnetic charge p = 0.1 and spin a = 0.1 are adopted.a The two parameters relate to the black hole charge q and magnetic field strength B. There are the particle energy E = 0.998, angular momentum L = 4.5, and initial separation r = 30.b The two parameters are the black hole charge q and particle energy E.