Chaos of Particle Motion near the Black Hole with Quasi-topological Electromagnetism

We explore the chaotic behavior of particle motion in a black hole with quasi-topological electromagnetism. The chaos bound is found to be violated in the higher order expansion of the metric function and the electric potential near the horizon. We draw the Poincare sections of particle motion corresponding to the chaos bound violated and non-violated cases, respectively. Then we study the relationship between the maximal Lyapunov exponent \lambda_s defined by the static equilibrium and the Lyapunov exponent of the particle geodesic motion near the Reissner-Nordstrom(RN) black hole and the black hole with quasi-topological electromagnetism. We find an interesting relationship between the Lyapunov exponent \lambda_{ph} of photon's radial falling into the black hole and the maximal Lyapunov exponent \lambda_s. For the black holes whose metric function increases monotonically with radius outside horizon, this leads to \lambda_{ph} \geq 2\lambda_s.


I. INTRODUCTION
Chaos is an important nonlinear phenomenon, which describes the violent response of the dynamic system to perturbations. Chaos of particle motion near the black hole has been studied for a long time, such as using the Melnikov method to find the chaotic orbit [1], distinguishing the chaotic orbits of particles or strings from periodic orbits [2][3][4][5][6][7][8][9][10][11][12][13][14] and probing the instability of particles motion [15][16][17][18][19]. Recently, much attention has been paid to the quantum chaos originated from the AdS/CFT correspondence [20]. In 2015, Maldacena, Shenker and Stanford pointed out in quantum field theory that for a quantum system like a black hole, its Lyapunov exponent λ should have an upper bound [20] λ ≤ 2πT, (1) where T is the Hawking temperature of the black hole in the natural unit = 1. Note that the Hawking temperature T is related to the surface gravity κ via κ = 2πT . In turn, the chaos bound can be recast as λ ≤ κ. The chaos bound is closely related to the chaotic behavior in the strongly correlated system and the physics of black holes. The chaos bound can also be studied through examining the geodesic motion of probing particles near black holes. Susskind proposed that for a neutral particle falling radially toward a black hole, its Rindler momentum will increase exponentially with the exponent equaling to the surface gravity κ [21], indicating the chaos bound is saturated. On the other hand, when a particle in its static equilibrium near the horizon with an external potential, the "maximal" Lyapunov exponent λ s can be obtained [22] with the relation λ s = κ at the horizon. However, the chaos bound could be violated for charged particles outside charged black holes [23] by studying the near-horizon expansion.
Our motivations of this paper are to study the Poincaré section and the chaos of particle motion for a very special black hole solution with quasi-topological electromagnetism obtained in [24]. We hope to study the Lyapunov exponent of particle motion near this black hole to understand the chaos bound and black hole chaos more. Interestingly, we find that the chaos bound can be violated for this special black hole. And in the system where chaos bound can be violated, the chaos is indeed strengthened.
Subsequently, we calculate the Lyapunov exponent λ c of the circular geodesic motion of natural particles near black hole, and find that λ c is not constrained by the "maximal" Lyapunov exponent λ s . Meanwhile, we are inspired by [25][26][27] in which the authors pointed out that as the speed of particle cannot reach the speed of light, the growth rate of the particle Rindler momentum will not saturate the chaos bound, or even increase exponentially. There may be a connection between particles' radial falling, the speed of light and chaos bound. So we study the radial falling of particles and find an interesting relation between the Lyapunov exponent λ ph of photon's radial falling and the "maximal" Lyapunov exponent λ s . For the black hole whose metric function increases monotonically outside horizon, we can see there is λ ph ≥ 2λ s . This paper is organized as follows. In Section II, we briefly review the black hole chaos and motivate why the study of black hole chaos is interesting. In Section III, we provide some information about the black hole with quasi-topolopical electromagnetism [24]. In Section IV, we discuss the static equilibrium of charged particles near the special dyonic black holes to check whether the chaos bound λ ≤ κ will be violated, and perform a numerical analysis of a toy model. In Section V, we study the geodesic motion of neutral particles, including circular geodesic motion and radially falling. In Section VI, we briefly summarize the main results. In Appendix A, the static equilibrium of a charged particle outside the RN black hole is studied by the fast Lyapunov indicator. In Appendix B, we discuss the static equilibrium with only gravity. In Appendix C, we calculate the "maximal" Lyapunov exponent λ s by the Jacobian matrix method. In Appendix D, we consider some details about the "maximal" Lyapunov exponent λ s outside these black holes. In Section E, we discuss in detail the parameters for the Poincaré section in a toy model.

II. BLACK HOLE CHAOS: A BRIEF REVIEW
Gravity theory is a nonlinear theory, and chaos should naturally exist in it. As the black hole is an important research object in gravity theory, the chaos near the black hole has high research value. Especially after the AdS/CFT duality is proposed, studying the connection between black hole theory and quantum chaos can help us further understand gravity theory and the essence of chaos. In this section, we briefly review and disscuss black hole and the related issues.

A. Some indicators for chaos around black holes
In the early studies on black hole chaos, searching for the behavior of chaos is the main task. The non-integrability of dynamical equation in black holes was studied by some analysis methods [1]. Some indicators such as the Lyapunov exponent, the fast Lyapunov indicator and the Poincaré section were also proposed to check the existence of chaos near the black hole.

Lyapunov exponent
The Lyapunov exponent describes the average rates of expansion and contraction of two adjacent orbits in the classical phase space, which can be defined by where x is the distance between the two orbits. The Lyapunov exponent can be used to describe the perturbation's exponential increasion in chaotic motion. The positive Lyapunov exponent means the existence of chaos. For the equilibrium of particles outside the black hole, the Lyapunov exponent can be obtained by the Jacobian matrix [15][16][17][18][19]. The equation of motion of particles can be schematically written as X i is the coordinates and F i (X j ) is a function to be determined. Considering a certain orbit, we can linearize the equation of perturbation where is the Jacobian matrix. When the equilibrium of particles outside black holes is considered, the Lyapunov exponent can be given by λ = det(K ij ).

The fast Lyapunov indicator (FLI)
The fast Lyapunov indicator (FLI) is a more effective tool to search for chaos [28]. FLI is usually defined by considering the difference between two nearby trajectories of particles, which is the so-called two-particle method or two-nearby-trajectory method. For chaotic motion, even for weak chaotic case, FLI will grow at an exponential ratio. In general relativity, F LI(τ ) can be described as [29,30] where d(τ ) = |g µν x µ x ν |, x µ is the deviation vector between two nearby trajectories at proper time τ . To avoid two orbits expand too fast, the sequential number k of renormalization is considered (the value of k can take 0,1,2...). In this paper, we use FLI to analyse the perturbation's growth in the static equilibrium of charged particles in Appendix A.

Poincaré section
The Poincaré section is the most commonly used tool to analyze dynamical systems, which automatically follows the non-integrability analysis in examining the existence of chaos. It can be defined as the intersection of a given hypersurface and motion trajectory in high-dimensional phase space (d≥3). Using the Poincaré section, the periodic orbits, quasi-periodic orbits and chaotic orbits can be clearly distinguished in dynamic systems.
Using these indicators, many chaotic phenomena in black holes have been studied, such as the motion of spining particle [8], scalar particle [11], and even strings [12][13][14] and the chaos in the background of a massive magnetic dipole [9], rotating black holes [5,10], the black holes with halo [5,6]. However, there were only the chaotic phenomenon near the black hole shown. The essential characteristics of black hole chaos are still to be discovered.

B. More on black hole chaos
Black hole chaos is an important subject in the study of black hole theory, which can be related to many physical problems. Actually, there are many physical problems related to black hole chaos. There are many phenomena related to black hole observation that have attracted much attention, such as photon sphere [31], shadow and gravitational lensing [32]. It is natural to discuss black hole chaos in these problems. The quasinormal modes (QNMs) of black holes describe the damping and oscillation of the gravitational waves emitted by the perturbation of black hole [33][34][35][36][37][38][39]. In [15][16][17][18][19], the Lyapunov exponent of circular null geodesic motion was related to QNMs, which gives an easier method to calculate the QNMs.
From the momentum-size duality proposed in [21], we can see that the black hole chaos is a bridge between gravity and quantum theory. Follow this duality, the connections between momentum and quantum complexity were explored in [40][41][42][43], where the quantum complexity [44] is a very meaningful concept that can be related to many physical topics such as the action of black holes [45,46], the shock wave geometry [47], the accelerated expansion of the universe [48] and the partition function [49]. There are also many more problems related to black hole chaos, such as black hole thermodynamics [50,51], acoustic black holes [52] and gauge/gravity correspondence [53]. Therefore, it is worthwhile to study the chaos of particles near the black hole with quasi-topological electromagnetism.

III. A SPECIAL DYONIC BLACK HOLE
Let us consider a black hole with quasi-topological electromagnetism proposed by [24]. The authors in [24] considered a 4-dimensional gravitational theory including the pure cosmological constant Λ 0 and the minimum coupling electromagnetic interaction. Its Lagrangian is given by [24] where α 1 and α 2 are coupling constants and The corresponding Maxwell's field equation is The Einstein's field equation can be written as This theory yields a static dyonic black hole solution where the metric dΩ 2 2, corresponds to a two-dimensional hyperbolic, torus and sphere, respectively with takes values -1, 0, 1. The metric function can be expressed as [24] where M is the mass of black hole, Λ 0 is the cosmological constant, q is the electric charge, p is the magnetic charge and 2 F 1 is a hypergeometric function.
The electric and magnetic charges can be obtained from electromagnetic tensorF µν which is defined by Maxwell's field equation Eq.(8) So the electric and magnetic potential functions are given by [24] Φ e (r) = When the constant parameters (M, p, q, Λ 0 , α 1 , α 2 ) take proper values, different types of black hole solutions can be obtained. Set α 1 = 1, Λ 0 = 0, = 1, two types of dyonic black holes can be obtained as the constants (α 2 , q, p) in the general solution Eq.(11) take appropriate values. Note that in these black hole solutions, as M takes different values, the black hole will also have different properties.
• Case 1: black holes with at most two horizons Consider a special condition of Eq.(11) with the parameters The corresponding function and electric potential functions become It is one of the first type black hole solutions with a constant M 0 (its value is related to the value of q, p, α 2 and the concrete definition of M 0 can be found in [24] there is a naked singularity at r = 0. For Eq. (14), the corresponding constant M 0 = 1.6372384. In the subsequent calculation process, for the first type of black hole solutions, we will consider an example with two horizons. In particular, M = 1.996, two horizons locate at r − = 0.28617, r + = 3.00002.
For this black hole, the metric function f (r) is monotonically increasing outside the outer horizon and the rate df (r) dr is monotonically decreasing.
• Case 2: black holes with at most four horizons There is another special case for Eq.(11) with parameters: .
The corresponding metric and electric potential functions can be given by The second type black hole solutions with three constants: M − , M + , M 0 (their values depend on the values of q, p, α 2 and the concrete definition of them can be found in [24]). Similarly, different situations depend on the value of M:  For this black hole, the metric function f (r) is not monotonically increasing outside the outer horizon, and the Newton potential has equilibriums at r = 2.8480 and r = 6.1421. So it is possible for particles to maintain a static equilibrium at these positions without any external forces, where the Lyapunov exponent is discussed in Appendix B.
Case 2-2: Black hole with cosmological horizons where r − , r + are the two horizons of this black hole, and r 1 , r 2 are the so-called cosmological horizons. The metric function f (r) increases firstly and then decreases between the outer event horizon and the inner cosmology horizon. This black hole also has an equilibrium of the Newton potential at r = 3.06723. For this black hole, the metric function f (r) is monotonically increasing outside the outer horizon, but the rate df (r) dr increases first and then decreases outside the outer horizon.
Actually, when the parameters of black holes satisfy different constraint equations, the black holes can be classified into these different types. For the same type of black holes, they have same properties. It is reasonable to follow the parameter values in [24].

IV. STATIC EQUILIBRIUM OF CHARGED PARTICLES OUTSIDE THE BLACK HOLE
In this section, we consider the general discussion of charged particles outside the black hole and derive the Lyapunov exponent. We consider a 4-dimensional static spherically symmetric black hole At the horizon, the surface gravity κ is where the prime " " denotes derivative with respect to r. The Lagrangian of a particle near this black hole can be written as 1 where m is the mass of the particle and V (r) is the external potential in the radial direction. When V (r) < 0, the potential V (r) can provide the particle with a repulsive force away from the black hole, so that the particle may not fall into the black hole and maintain static equilibrium near the horizon. In this section, we take the static gauge τ = t, so the dot "·" denotes derivative with respect to the proper time " t " in this section. When the particle stays in its static equilibrium near the horizon, the particle's Lagrangian Eq.(18) describing its radial motion can be reduced to Note that when the particle maintains static equilibrium, we haveṙ 1. After expanding the Lagrangian in terms ofṙ, we can obtain the effective Lagrangian where V ef f is the effective potential and at the particle's equilibrium position there is V ef f = 0. After expanding V ef f at the particle's static equilibrium position r = r 0 , we have the effective Lagrangian satisfying the relation where λ is the Lyapunov exponent obeying In the equilibrium position r = r 0 , V ef f =0. If V ef f has the maximum V ef f < 0, the static equilibrium is unstable and λ 2 > 0. If V ef f has the minimum as V ef f > 0, the static equilibrium is stable and λ 2 < 0. From the effective Lagrangian Eq. (21), we can see that the particle should follow the equation of motion at the equilibrium positionr This equation of motion tells us the particle's trajectory in the radial direction, that is to say If λ is real, the exponential increasion of r may imply the existence of chaos. Actually, when there is no perturbation in other directions, the position r = r 0 (where V ef f has the maximum) will be a separatrix of the phase space, which can provide a "maximal" Lyapunov exponent here [22,[54][55][56]. When the external potential is strong enough, the position where V ef f has the maximum can be near the horizon.

A. Chaos bound at the horizon
Returning to the black hole solution Eq.(11), we calculate the Lyapunov exponent when a charged particle maintains its static equilibrium near the horizon of these black holes and further verify the chaos bound λ ≤ κ. Here, the external potential we are considering is provided by the electric field, and its form is where e and m are the charge and mass of the particle, and Φ e (r) is the electric potential function. For the particle's static equilibrium at r = r 0 , the Lyapunov exponent λ satisfies the expression [23] When the static equilibrium position is infinitely close to the black hole outer horizon r = r + , the Lyapunov exponent λ should have an upper bound, that is, the chaos bound λ ≤ κ. The content of this section is mainly to calculate the Lyapunov exponent λ of the charged particle near the outer horizon of these black holes and to check whether they satisfy the chaos bound. Some interesting phenomena would occur when we pay attention to the near-horizon behavior shown in [23]. The metric function and potential function of the black hole can be expanded to the second-order on the black hole horizon r = r + : Then, substitute Eq. (27) into Eq. (17) and Eq.(26) Expanding λ 2 in Eq.(28) at the horizon r = r + , we can obtain which can be rewritten as [23] From Eq. (30), it can be seen that when γ > 0, the chaos bound λ ≤ κ could be violated. The main reason we consider this dyonic black hole here is that it can provide some special conditions. Since the calculation is very complicated, here we only show the main results of the calculations in the tabular form as shown in Table I, II. • Case 1 We consider the black hole solution given in Eq. (14) and substitute it into Eq.(17) and Eq.(26) to calculate the surface gravity of the particles at the black hole horizon and the Lyapunov exponent λ. Then, the parameter γ was calculated through Eq.(30). The calculation results are shown in Table I.
Similarly, we also calculate the black hole solution given in Eq. (15), and show the results in Table II. There is always λ = κ at the outer horizon r+ for ever case. But as we see, γ is positive for the black hole in Case 2-2, which means that there may be a violation of the chaos bound λ ≤ κ.
The above are the results of all our calculations about λ, κ and γ for these dyonic black holes in Eq. (14) and Eq. (15). These calculations show that λ = κ is universally established on the black hole outside horizon, which means that the chaos bound of λ ≤ κ is satisfied.
However, when considering the near-horizon behavior of the black hole and expanding the correlation function to the second order, we find the case of γ > 0, which indicates that there is a violation of the chaos bound. In order to examine this situation, we need to do more general study about the second-order expansion of these black holes.

B. More discussion on the expansion
To more clearly show the anomalous behavior of γ in the second-order expansion, we will discuss these two sets of metric functions and their corresponding potential functions above (Eq. (14) and Eq.(15)), and then plot γ as a function of M in Figure 1. There is always γ < 0 in Figure 1(a), which means the chaos bound λ ≤ κ is not violated for Case 1. However, as shown in Figure 1(b), for the black holes with M ∈ (6.73442, 6.9135) in Case 2, the chaos bound can be violated in the near horizon region.

C. Numerical analysis
To study the effect of the violation of chaos bound more, we consider a toy model that describes the near-horizon geometry with external potentials and study the Poincaré section of the particle motion in this model. The Lagrangian of particles can be written as where "·" denotes derivative with respect to the coordinate time t, the metric function f (x) = 2κ(x − x h ), κ is the surface gravity, x h is the radius of black hole horizon. (Note that the near-horizon expansion of the metric Eq. (11) yields the form f (x) = 2κ(x − x h ).) A(x) and B(y) are external potentials, ω is the coupling coefficient of particle with external potential 2 . From the Lagrangian Eq.(31), we can derive the generalized momentum The Hamiltonian of particle (that is, the energy E) is Correspondingly, the equation of motion can be written as Similarly as in [3,22], we consider the external potential similar to the harmonic oscillator potential where x c is the center position of A(x). We set x h = 0 and x c = 1, then two cases that the chaos bound is violated or not can be decided by adjusting the parameters a and b. To avoid the particle falling into horizon, the particle energy E should have an upper bound E max (See Appendix E for the detailed derivation of a, b and E max ). Here, we consider the values of a and b as following: 1. Chaos bound can be violated

Chaos bound not violated
The other parameters are chosen as κ = 1 2 and ω = E max = 30. The equation of motion Eq.(34) can be solved numerically. About the initial condition, we set y(0) = 0, P x (0) = 0, and the same x(0) for two cases, and the corresponding P y (0) can be obtained by Eq. (33). Then we can define the Poincaré section in (x, P x ) plane with y = 0 and P y < 0. In the following figures, we put the Poincaré section of the system (a=1.3, b=-0.3) where the chaos bound can be violated on the left, and the non-violation system (a=1.1, b=-0.1) on the right. In Figure 2, we set E = 18 and show the Poincaré section. The particle motion seems to be more closer to chaos in the system where the chaos bound cannot be violated since the section in Figure 2(b) is more irregular than in Figure 2(a). We should point out that the reason of this difference maybe the system(a = 1.3, b = −0.3) has a bigger potential A(x) than another. For the same E, the particle will have smaller kinetic energy which may weaken the chaos in particle motion. As for why the chaos is not strengthened in the system where the chaos bound can be violated, we speculate that the low energy results in it. The greater the energy E, the closer the particle can approach the horizon until it falls into the black hole. The violation of chaos bound we studied is a near-horizon behavior, so its reinforcement of chaos will work when the particles are close to horizon. To verify this, we examined the Poincaré section with higher energy. The situation of E = 25 is shown in Figure 3. In this figure, we can see the chaos is strengthened in the system with a = 1.3 and b = −0.3. The points of x(0) = 0.5 and x(0) = 1.25 form closed curves in Figure 3(a) which indicates quasi-periodic orbits, and in Figure 3(b) they are single points which represent periodic orbits. As we expected, the chaos in system (a = 1.3, b = −0.3) has been strengthened, because of the bigger E results in particle motion approaching horizon. Meanwhile, there is an opposite result of points with x(0) = 1.5, as the system (a = 1.3, b = −0.3) with stronger A(x). These points in Figure 3(a) form a closed curve, and they are dispersed which means it is a chaotic orbit in the system (a = 1.1, b = −0.1).
To further examine our intuition, we set E=29.9 in Figure 4, which will make the particles move close enough to the horizon to make the enhancement of chaos in system (a = 1.3, b = −0.3) more obvious. At the same time, we calculated more data. As shown in Figure 4, there are many chaotic orbits in Figure 4 values of x(0) in Figure 4(b) are regular which means these orbits are not chaotic. The chaos in the system (a = 1.3, b = −0.3) has been significantly strengthened. Although the stronger potential A(x) in the system (a = 1.3, b = −0.3) will weaken the chaos, we still see that the chaos is strengthened in Figure 4, which undoubtedly shows in a system where the chaos bound can be violated, chaos can indeed be strengthened.

V. GEODESIC MOTION OF NEUTRAL PARTICLES
In this section, we study the Lyapunov exponent of the geodesic motion of neutral particles near the black hole, including circular geodesic motion and radial falling. We extend the static equilibrium of charged particles in Section IV beyond the horizon, and take the "maximal" Lyapunov exponent as λ s which can be obtained from the static equilibrium. From Eq.(26), we obtain the formula of λ s for charged black holes Arbitrary location outside the horizon where the effective potential V ef f has a maximum, the Lyapunov exponent λ s obtained from static equilibrium is at its maximal value. When we consider λ s at the horizon, the value of λ s will return to the surface gravity κ. We propose a hypothesis here: the "maximal" Lyapunov exponent λ s is an inherent property determined by the nature of the black hole 3 . It is worthwhile to explore the connection between the "maximal" Lyapunov exponent λ s with the Lyapunov exponent of particle geodesic motion.

A. Circular geodesic motion of neutral particles
A 4-dimensional spherically symmetric black hole is considered here. Its metric is When we focus on a neutral test particle moving on the equatorial plane of this black hole (θ = π 2 ,θ = 0), its Lagrangian can be written as where the dot " . " denotes derivative with respect to the proper time " τ ". According to the generalized momentum expression p q = ∂L ∂q , we can obtain its generalized momentum as where E is the particle's energy and L is the angular momentum of the particle. Using the normalization with four-velocity g µν u µ u ν = η we can obtainṙ Note that η = +1, −1, 0 corresponds to space-like, time-like and null geodesics respectively. For the circular geodesic motion of neutral particles, we have the Euler-Lagrange equation Setting the phase space variables X i (t) = (p r , r) and considering the particles moving in a circular orbit with a radius of r = r 0 , we have two equations dp r dτ = ∂L ∂r and dr dτ = p r g rr .
Then, the Jacobian matrix of particle motion can be obtained from Eq.(4) The eigenvalues of this matrix can give the expression of the proper time Lyapunov exponent of circular geodesic motion λ p , and we observe that λ p satisfies [16] From the Lagrange's equation of geodesic motion and the formula d dτ the expression of ∂L ∂r can be written as For the circular geodesic motion, we have the circular geodesic condition [57] r 2 = ṙ 2 = 0. (49) Under this condition, we substitute Eq.(48) into Eq. (45), then the proper time Lyapunov exponent λ p in Eq. (45) can be reduced to If we consider an alternative form of Eq.(42) we can express the Lyapunov exponent λ c in term of coordinate time [15] In Eq. (39),ṫ = E f (r) , so we obtain the relation between Eq.(50) and Eq.(52) Both f (r) and E are real, so the properties of λ p and λ c are closely related. It is evident that the Lyapunov exponents λ p and λ c are the most important parameters directly verifying the stability of the motion. Only when λ p and λ c are both real, they result in unstable geodesic motion. In contrast, when one of the λ p and λ c is imaginary, the circular geodesic motion is stable; when λ p = 0 or λ c = 0, the circular geodesic motion is marginal, which means the circular geodesic motion can be easily broken. The Lyapunov exponent is a measure of the deviation in the time evolution of two adjacent trajectories in phase space, obviously it depends on the time coordinate. In the following, we will study the coordinate time Lyapunov exponent λ c of the neutral particle's geodesic motion.

Circular time-like geodesic
For time-like geodesic with η = −1, we havė The circular geodesic conditionṙ 2 = ṙ 2 = 0 yields [15] After substituting the metric functions of the four black holes given by Eq. (14) and Eq. (15) where M is the mass of black hole, and Q is the charge. For RN black hole, λ 2 c and λ 2 s can be calculated by Eq.(52) and Eq.(36) We want to explore the relationship between the circular time-like geodesic motion's Lyapunov exponent λ c and "maximal" Lyapunov exponent λ s defined from the static equilibrium, so we plot figures of λ 2 c and λ 2 s in the region where the circular geodesic motion exists. and Q = 1, λ 2 c decreases from the radius of the innermost circular orbit (the beginning of the coordinate axis). λ 2 s decreases too, but λ 2 s is always positive in contrast to λ 2 c . Near the innermost orbital radius, we have λ 2 c > λ 2 s . Away from the black hole, λ 2 c is negative, which means the circular geodesic motion is stable.
Firstly, we focus on the stability of circular motion which can be expressed by λ 2 c . For the RN black hole, as we show in Figure 5, λ 2 c decreases monotonically from the innermost circular orbit until it becomes zero at the innermost stable circular orbit, the circular time-like geodesic motion can only exist beyond the innermost stable circular orbit. In Figure 6, λ 2 c in Case 1 has a similar behavior as in RN black hole. For the black hole in Case 2-1, the circular geodesic motion exists in three discontinuous range: 2.19387 < r < 2.84801, 6.14209 < r < 6.44301 and r > 12.4354. Figure 7(a), the circular motion is always unstable. But in the Figure 7(b), there is λ 2 c < 0 corresponding to the existence of the stable circular time-like geodesic orbit in the range of 6.14209 < r < 6.44301 which is closer to the horizon than in Figure 7(c). Compared with RN black hole, the black hole given by Case 2-1 has a similar evolution trend of λ 2 c in Figure 7(c) with r > 12.4354. From Figure 7(b) and Figure 7(c), it seems to be a stable circular orbits closer to the horizon. We guess that this special property of Case 2-1 has a Newtonian potential equilibrium position that allows it to have a stable circular time-like geodesic closer to the horizon. For Case 2-2 in Figure 8, the circular orbits are unstable. In Figure 9, the behavior of λ 2 c is similar to that of RN black hole. , there is λ 2 c < 0 in the range r ∈ (6.14209, 6.44301) near the horizon, which means that in this range there can be a stable time-like circular geodesic motion. In (c), λ 2 c has a RN-like behavior.

As shown in
Then, we discuss the relationship between the circular motion's Lyapunov exponent λ c and the "maximal" Lyapunov exponent λ s . In [22], the authors regarded the innermost circular geodesic motion as an equilibrium with a repulsive potential given by particle circular motion. But as shown in Figure 5, Figure 6 and Figure 9, we can see that when the circular time-like geodesic motion approach the innermost circular orbit, λ s will not constraint λ c any more. The "maximal" Lyapunov exponent λ s seems to constraint λ c in Figure 7(a) 4 and Figure 8 with λ s = λ c at the position where Newtonian potential has an equilibrium. We think λ s restricts λ c in Figure 7(a) and Figure 8 because the particle is close to the equilibrium position of the Newtonian potential where particle motion is slow. In other words, we speculate the main reason for λ c > λ s around the innermost circular orbit in Figure 5, Figure 6 and Figure 9 may be when the particle approaches the innermost circular orbit its speed is close to the speed of light. There may be another reason which is the movement in other directions makes Lyapunov exponent growth. To verify our conjecture, we discuss the particle's radial falling in Section V B.

Circular null geodesic and stable photon sphere
To study the motion of the photon outside the black hole, we set η = 0 in Eq. (41). Then, Eq.(41) can be reduced toṙ From the circular geodesic conditionṙ 2 = ṙ 2 = 0, we have constraints on the circular motion of photons 2f (r cn ) = r cn f (r cn ) , where r cn is the radius of the circular null geodesic, which is defined by Eq.(59a). After taking the black hole solutions we consider Eq.(59a) to obtain r cn , we substitute the radius r cn , Eq.(59b) and Eq.(58) into Eq.(50) and Eq.(52) to calculate the Lyapunov exponent. From the calculation results of the Lyapunov exponent as shown in Table III, we can judge the stability of the circular null geodesics and verify the existence of stable photon spheres.
From Table III, we can see there is an unstable photon sphere at r = 4.72028 for Case 1. For Case 2-1, there are three photon spheres at r = 2.19387, r = 6.44301, r = 12.4354, and the one at r = 6.44301 is stable which agrees with the result in [24]. For Case 2-2 and Case 2-3, the photon spheres are unstable.

B. Radial falling
In this subsection, we study the geodesic motion of particles falling radially into the black hole. The metric of a 4-dimensional spherically black hole can be written as Eq. (37). When considering that the particle falls freely towards the black hole in the radial direction, Eq.(41) can be reduced tȯ then we can obtain When we expand around each point on the particle trajectory (denoted as r f ), there is Similarly, we can get an exponential growth form of the coordinate r. Here, we discuss in two cases: massive particle For the massive particle falling from infinity, there is time-like geodesic. Setting η = −1, E = 1, we can obtain the Lyapunov exponent λ mp from Eq.(62) photon For the photon, there is null geodesic. With η = 0, we have the Lyapunov exponent λ ph from Eq.(62) We can see that there is λ mp ≈ λ ph = 2κ near the horizon 5 . We study the relationship between these two exponents outside the horizon and the "maximal" Lyapunov exponent λ s (λ s can be calculated via Eq.(26)).
For simplicity, we first study the RN black hole in Eq.(56), then Eq.(63), Eq.(64) and Eq.(26) can be rewritten as 5 Near the horizon r + , we have f = 0 for Eq.(62), then Eq.(62) can be rewritten as dr dt ∼ f (r − r + ). Because of the surface gravity κ = 1 2 f r=r + , we can obtain r − r + = e 2κt at the horizon agreeing with [27,52]. In Figure 10, we show the relationship between λ mp , λ ph and λ ph in RN black hole and the metric function of RN black hole. There are some interesting phenomena in particle motion near RN black hole. For massive particles, in the region near the horizon, the closer to the horizon, the closer λ mp approaches λ ph . The reason may be the speed of the massive particle approaches the speed of light. For photon which is radial falling into the RN black hole, there is always It is similar to λ OT OC ≥ 2λ chaos , which is a formula about chaos [58]. Considering the relationship between null geodesic and shock waves, we speculate that there may be some relationship between these two formulas, or even equivalent.
To study the property of Eq.(66), we perform the same calculation for the black holes given by Eq. (14) and Eq.(15) and show the results in Figure 11-14. From Figure 11 and Figure 14, there is always λ ph ≥ 2λ s . But it's not in Figure  12 and Figure 13. We think the reason is that the metric function f (r) is not monotonically increasing outside the horizon, and we could understand the relationship between λ ph ≥ 2λ s and λ OT OC ≥ 2λ chaos by studying more examples.

VI. CONCLUSION AND DISCUSSIONS
In summary, we mainly perform relevant calculations on the black hole solutions given in [24]. Firstly, we introduce black hole chaos to show why this is an interesting topic. Secondly, we review the definition of λ s which is the "maximal" Lyapunov exponent obtained from static equilibrium [22,23]. We calculate the Lyapunov exponent when the charged particle remains static equilibrium near the horizon of these black holes given in Section III. The results show that the chaos bound does not seem to be universally satisfied. When considering the higher-order expansion term, for black holes with specific parameters, the chaos bound can be violated. Thirdly, we consider a toy model with adjustable parameters to examine the difference between the chaos violated and non-violated cases. From the analysis of Poincaré section, we can see the chaos could be strengthened in the system where the chaos bound can be violated, which is a novel study about the black hole chaos. After studying the static equilibrium of the particles, we turn to the geodesic motion of the particles. We study two types of geodesic motion: circular geodesic motion and radial falling. For circular geodesic motion, we obtain the Lyapunov exponent λ c by the Jacobian matrix and there are some stable circular orbits. At the same time, we find that these Lyapunov exponents λ c of circular geodesic motion do not seem to be constrained by the "maximal" Lyapunov exponent λ s defined from static equilibrium. The reason may be that the particle speed approaches the speed of light. This is a topic worth discussing. Then we study the particle's radial falling, and find that the Lyapunov exponent λ mp of the mass particle will approach the Lyapunov exponent λ ph of the photon in the region near the horizon. We speculate the reason may be that the falling speed of the mass particle is close to the speed of light. Note that the discussions are carried out in the geodesic dynamics. It would be interesting to include the backreaction of particle motion to the background spacetime and study the many-body effect of the chaotic behavior. We find that for RN-like black holes or black hole examples where the metric function f (r) increases monotonically outside the horizon, the relation λ ph ≥ 2λ s is established.
Chaos near black holes is an important subject in contemporary physics research. As far as we know, the nature of the chaos near the black hole will depend on the black hole. Therefore, studying as many different black holes as possible may help us understand chaos more clearly. It may be worth studying the higher-dimensional black hole with quasi-topological electromagnetism given in [59]. ACKNOWLEDGEMENT We would like to thank Hong Lü, Shu Lin and Qing-Bing Wang for helpful discussions. The work was partially supported by NSFC, China (grant No.11875184). To study the perturbation of particles at static equilibrium more clearly, we use the fast Lyapunov indicator to analyze the perturbation growth. Taking the parameter k = 0 in Eq.(6), we have where d(τ ) = |g µν x µ x ν |. x µ is the deviation vector between two nearby trajectories at proper time τ , and in the computation, we choose two particle trajectories with the initial state at the equilibrium position r 1 = r 0 and the non-equilibrium position r 2 = r 0 − ( is the perturbation). We consider a 4d RN black hole metric here where the metric function is f (r) = 1 − 2M r + Q 2 r 2 , and the potential function is A t (r) = 2Q r . The trajectory of particles motion can be numerically solved from the equation of motion Eq.(C18). The parameters of RN black hole here we simply take the mass M = 2 and the charge Q = 1. For particle 1, which is initial at the equilibrium position r 1 = r 0 , its initial value condition is (t, π t , r, π r ) = (0, π t1 , r 1 , 0), where π t1 = − f (r 1 ) + f (r 1 ) 2 π 2 r − qA t (r 1 ). For particle 2, which is initial at the non-equilibrium position r 2 = r 0 − , its initial value condition is (t, π t , r, π r ) = (0, π t2 , r 2 , 0), where π t2 = − f (r 2 ) + f (r 2 ) 2 π 2 r − qA t (r 2 ). The charge of two particles are the same, as defined by the equilibrium At(r1) . The results are shown in Figure 15. As shown in Figure 15, for the perturbation on the static equilibrium of charged particles, the chaotic behavior can be found by FLI. The results show that the closer to the horizon, the faster the FLI increases, which is also consistent with the behavior of the "maximal" Lyapunov exponent λ s outside the black hole. From the analysis of the effective potential curve, we can know that V ef f has maximum and minimum values. We will study λ at the maximum position of V ef f shown in Figure 16. Since this position is far from the black hole horizon, the "maximal" Lyapunov exponent λ s should not be given by surface gravity κ. We also consider the calculation from the geodesic motion. For the massive particle making a static equilibrium with only gravity, Eq.(41) can be reduced toṙ With the static conditionṙ 2 = ṙ 2 , we can obtain the coordinate time Lyapunov exponent λ c from (52) These two calculations have the same result.
There is a static equilibrium condition for the charge of particle at the equilibrium position r = r 0 , Taken (r, P r ) as the phase space variables, for the dynamic system dr dt = F 1 (r, P r ), dP r dt = F 2 (r, P r ), the components of the Jacobian matrix K ij can be given (f (r) (P 2 r f (r) + 1)) 3/2 , K 21 = ∂F 2 ∂r = f (r) 2 − 2f (r) P 2 r f (r) + 1 2qA t (r) (r) f (r) (P 2 r f (r) + 1) + 2P 2 r f (r) + 1 f (r) 4 (f (r) (P 2 r f (r) + 1)) At the equilibrium location r = r 0 , where P r = dPr dt = 0, the Jacobian matrix K ij can be reduced to Substituting q = − ( √ f (r)) At(r) r=r0 into the Jacobian matrix K ij , we can obtain the Lyapunov exponent satisfies which is same as Eq.(36).

The equivalent form of Lagrangian
Considered another form of the Lagrangian when only the radial motion is considered, Eq.(C14) can be reduced to The generalized momentum can be defined as π r =ṙ f .
Appendix D: V ef f and the "maximal" Lyapunov exponent λs In studying the static equilibrium of particles to obtain the "maximal" Lyapunov exponent λ s , the most important point is to ensure that the static equilibrium position is just at the maximum value of the effective potential V ef f .
For the static equilibrium of charged particles, there is V ef f = 0 can be used to determine the static equilibrium position, and on the static equilibrium position, there is When we consider Eq.(D2), we can obtain Only when V ef f is negative, the effective potential V ef f will have a maximum value at the static equilibrium position, which corresponds to the "maximal" Lyapunov exponent λ s . As shown in Figure 17, for the black hole represented by Case 2-1, within a certain range r > 4.37690, there is always a minimum value for effective potential at the equilibrium position, where the "maximal" Lyapunov exponent cannot be defined.