{'fraction_non_alphanumeric': 0.046614291058735505, 'fraction_numerical': 0.025614166354907094, 'mean_word_length': 3.8030666027791087, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Black Holes With Many Horizons in the Theories of Nonlinear Electrodynamics

In this article, we construct exact black hole solutions with many horizons (more than number two) in the Einstein-nonlinear electrodynamic theories. In particular, we acquire the explicit expression of nonlinear electrodynamic Lagrangian for the 3-horizon black holes. Then we make the investigations of 3-horizon black holes on the horizons, the null and timelike geodesics, the Love numbers and the thermodynamics.


I. INTRODUCTION
In order to remove the inconsistencies of infinite selfenergy of the point charges in the Maxwell's electrodynamics, Max Born firstly proposed the nonlinear electrodynamic (NLE) theory [1]. The theory was subsequently extended in collaboration with Leopard Infield [2] to the so-called Born-Infield theory. It is amazing that the Born-Infield theory was rediscovered in the low energy limit of the string theory [3] half a century later. It is found that the Born-Infield parameter b is related to the string tension α ′ via 2πα ′ = 1/b [4]. On the other hand, not long after the proposal of Born and Infield, Warner K. Heisenberg and Khans H. Euler [5] achieved the one-loop correction of quantum electrodynamics to Maxwell's Lagrangian.
Take into account gravity, one suppose the NLE theories may erase the black hole singularities. However, it is not the case. It is found neither electrically charged Einstein-Born-Infield black holes [6][7][8][9] nor electrically charged Einstein-Euler-Heisenberg black holes [10,11] are regular. The first and remarkable regular black hole solution was written down by Bardeen [12] by hand. It was interpreted very latter by Ayon-Beat and Garcia [13,14] as the solution of a particular Einstein-NLE theory. Later on, a lot of regular black hole solutions were found within Einstein-NLE theories based on various Lagrangian functions, such as logarithmic [15], hyperbolic tangent [16], power [17,18], exponential [19], de Sitter or anti-de Sitter asymptotic [20] and so on. Finally, the black holes with charged scalar hairs in Einstein-NLE is considered in [21]. In the respect of regular black holes in NLE, an important point is the existence of a no-go theorem [22,23]. It is shown that there is no such Lagrangian function which has the Maxwell weak field limit that the resulting black hole has the regular center. In order to circumvent the no-go theorem, an interesting proposal is to consider a kind of phase transition on a certain sphere, * gaocj@nao.cas.cn outside of which there is a pure electric field but inside of which the field is pure magnetic [24].
On the other hand, it is well-known that the Schwarzschild black hole has one event horizon. The Reissner-Nordstrom (RN) black hole, the Kerr black hole and the Kerr-Newman black hole have two horizons among which one is the event horizon and the other is the inner Cauchy horizon. So one interesting question that one would ask is " Can a black hole have horizons more than two? " The answer is yes. In fact, regular and non-regular multi-horizon (more than two horizons) black holes in Einstein-NLE theories have been presented in Ref. [25,26].
It is found that the multi-horizon black holes have very rich physics. For instance, different from the standard one-horizon black holes, the multi-horizon black holes show up not only the Hawking evaporation but also antievaporation (or related instability phenomenon). In detail, one can consider the limit where the radius of one horizon coincides with that of another horizon for multihorizon black holes. This is called the Nariai limit. Usually, the radius of the horizon decreases by the Hawking radiation. However, in case of the Nariai limit, the radius can increase due to the quantum effects. This is the well-known anti-evaporation effect [27][28][29].
However, to the best of our knowledge, none of th solutions in Ref. [25,26] brings us with an analytic expression for NLE Lagrangian. In this article, starting from the most general, analytic form of Lagrangian (the sum of infinite series of Maxwell invariant) for NLE, we shall seek for the multi-horizon black hole solutions. In particular, we find and investigate the solution for a 3-horizon black hole with explicit form of NLE Lagrangian.
The paper is organized as follows. In Sec. II, we report that we find the solutions for black holes with many horizons, in particular with three horizons, in NLE. In Sec. III, we give a detailed analysis on the horizons for the 3-horizon black holes. In Sec. IV, we derive the geodesic equation in the equatorial plane for both null and timelike geodesics in the background of 3-horizon black hole. Sec. V and Sec. VI are devoted to the investigations on radial geodesics, general geodesics, stable circular orbits and innermost stable circular orbit (ISCO). In Sec. VII, we consider the motion of test charged particles in the background of 3-horizon black hole. In Sec. VIII, we show the Love number of 3-horizon black hole is vanishing. This reveals the 3-horizon black hole is totally rigid. In Sec. IX, we make a study on the thermodynamics for 3-horizon black holes. Finally, we give the conclusion and discussion in Sec. X. Throughout this paper, we adopt the system of units in which G = c = = 1 and the metric signature (−, +, +, +).

MANY HORIZONS
We consider the Einstein theory coupled with the nonlinear electromagnetic field which has the action with Here R is the Ricci scalar and L EM is the extended Maxwell Lagrangian. A µ is the Maxwell field. α i are dimensional constants and have the dimension of lengh 2(i−1) . i is a positive integer. The physical motivation for taking the extended Maxwell Lagrangian Eq. (2) for the gauge field is that it can cover nearly all the known proposals for NLE, for example, the Born-Infield-Lagrangian [1,2], the Euler-Heisenberg Lagrangian [5], the power-law Maxwell Lagrangian [17,18] and so on. When α 1 = 1 and α i = 0 (for i > 1), it reduces to the Einstein-Maxwell theory. The variation of the action with respect to the metric gives the Einstein equations On the other hand, the variation of the action with respect to the field A µ gives the generalized Maxwell equations We shall look for the static, spherically symmetric black hole solutions in the theory. To this end, we take the ansatz of the metric and the Maxwell field Then the Einstein and Maxwell equations give Here Q, as an integration constant, is nothing but the electric charge of the RN black hole. The prime denotes the derivative with respect to r. Solving Eq. (7), we obtain It seems rather difficult to solve remaining equations. But in fact, the solutions with electric charge have been solved in a general form by Pellicer and Torrence in [30]. The magnetic counterpart was solved by Bronnikov in [31]. Here we will look for the solutions in the form of series (the most general form). So we expand U and Φ as follows such that Φ is asymptotically vanishing and the spacetime is asymptotically Minkowski. Here b i and c i are constants. Taking Eq. (11) into account and substituting Eqs. (12) into Eqs. (8)(9)(10), we obtain the non-vanishing constants when α 1 = 1 and These form a solution of power series. For this solution, there are only two integration constants, Q and c 1 . Q is the electric charge and c 1 = −2M (M is the mass of the black hole). Basing on this solution, we have the following conclusions.

A. 2-horizon black hole
When we obtain Let c 1 = −2M , then it is the RN black hole. M and Q represent the mass and electric charge of the black hole, respectively. In general, RN spacetime has two horizons. The corresponding Lagrangian L EM is given by

B. 3-horizon black hole
When we obtain the solution Given the expressions of Φ, U above and using the Lagrangian-generation-method in [32], we can construct the corresponding Lagrangian L EM with We note that the Lagrangian is dependent on not the mass M and charge Q, but only the coupling constant α 2 . Expand Eq. (27) in series of F 2 , we obtain L EM = F 2 + α 2 F 2 2 + 4α 2 2 F 2 3 + 24α 3 2 F 2 4 +176α 4 2 F 2 5 + 1456α 5 2 F 2 6 + · · · , (29) which are consistent with the coefficients in Eq. (25). In order that the root in Eq. (28) always makes a sense, we should require because of F 2 = −2Φ ′ 2 < 0. We shall find shortly later this spacetime can have three horizons when α 2 is negative. As an example, we plot the locations of three horizons for M = 1, Q = 0.7, α 2 = −0.001 in Fig. 1. We point that there is no physical reasons for the choice of above parameters. We simply choose the parameters by hand.

C. 4-horizon black hole
When we obtain the solution Using the Lagrangian-generation-method in [32], we find the corresponding Lagrangian L EM is where η is determined by Similarly, the Lagrangian is dependent on not the mass M and charge Q, but only the coupling constants α 2 and α 3 . If the parameters, M, Q, α 2 and α 3 are properly chosen, the spacetime can have four horizons. As an example, we plot the locations of four horizons for M = 1, Q = 0.74, α 2 = −0.0009, α 3 = 3.14 · 0 −6 in Fig. 2.
Same as the scenario of 3-horizon black holes, there is no physical reasons for the choice of above parameters. They are chosen simply by hand. Similarly, we can obtain black holes with 5-horizon, 6-horizon and so on. But in the next sections, we shall focus on 3-horizon black holes. Our lagrangians have intrinsics advantage over previous studies [25,26]. They are not dependent on the mass M , charge Q, or the ratio between charge to mass Q/M .

III. HORIZONS OF 3-HORIZON BLACK HOLE
The positions of horizons for 3-horizon black holes are determined by where V is defined by In order to determine the number of horizons, we make an analysis on V in the first place. In Fig. 3, we plot five critical curves for V with Q 2 > 25 24 M 2 , Q 2 = 25 24 M 2 , M 2 < Q 2 < 25 24 M 2 , Q 2 = M 2 and Q 2 < M 2 , from up to down, respectively. We note that r = 0 is the curvature singularity. When Q 2 ≥ 25 24 M 2 , V is an increasing function. When Q 2 < 25 24 M 2 , there are in general two local extremals, the local maximum V max = V | r=r− and local minimum V min = V | r=r+ at and respectively. In particular, when M 2 < Q 2 < 25 24 M 2 , the local minimum is positive. On the other hand, when Q 2 < M 2 , the local minimum is negative. Finally, when Q 2 = M 2 , the local minimum is zero. Therefore, we have the following conclusions.
In this case, V is an increasing function of r. Therefore, if α 2 < 0, there is only one black hole event horizon. But if α 2 > 0, there is no horizon and the singularity is naked.
In this case, we have V max > 0 and V min > 0. Thus we have the following conclusions.

If
there is only one black hole event horizon. We conclude there are two horizons. We conclude α 2 < 0.

If
there are two horizons. We conclude α 2 < 0.
In this case, we have V max > 0 and V min = 0.
there is only one black hole event horizon. We conclude there are two horizons. We conclude α 2 < 0.

If
there is one horizon. Since V min = 0, we conclude α 2 = 0. It is nothing but the extreme Reissner-Nordstrom black hole.

If
there is no horizon. We conclude α 2 > 0.
In this case, we have V max > 0 and V min < 0. Thus we have the following conclusions.
1. If there is only one black hole event horizon. We conclude there are two horizons. We conclude α 2 < 0. 3. If there are three horizons. We conclude α 2 < 0.

If
there is one horizon. We conclude α 2 > 0.

If
there are two horizons. We conclude α 2 > 0. 6. If there is no horizon and the singularity is naked. We conclude α 2 > 0. In summary, when α 2 < 0, the spacetime can have three horizons. We shall mainly focus on α 2 < 0 in the next sections.

SPACETIME OF 3-HORIZON BLACK HOLE
In this section, we give the equation of motion for null and timelike geodesics. In this section and later sections, we shall closely follow the way of Chandrasekhar [33]. The equations governing the geodesics in a spacetime with the line element can be derived from the Lagrangian where τ is the some affine parameter along the geodesic. For timelike geodesics, τ can be identified with the proper time, s, of the co-moving observer along the geodesic. For the 3-horizon black hole spacetime, the Lagrangian is where the dot denotes the derivative with respect to the τ . Using the Euler-Lagrange equation, we obtain the equation of motion for the geodesic in the equatorial plane, θ = π/2, where E is the energy (inclusive of the rest energy) and L is the angular momentum about the axis to the equatorial plane for the particle. By rescaling the affine parameter τ , we can arrange that 2L has the value of +1 for timelike geodesics. For null geodesics, L has the value of zero. By considering r as the function of φ and letting u = 1/r, we obtain the basic equation of the problem as in the analysis of Keplerian orbit in the Newtonian theory This equation determines the geometry of the geodesics in the equatorial plane.

A. Radial null geodesics
The equations governing the radial null-geodesics can be obtained by setting L = 0 and L = 0 in Eq. (60) and Eq. (61), thuṡ Therefore, In order to solve this equation, we define the tortoise coordinate and rewrite U as follows Here r 3 , r 2 , r 1 denote the radii of three horizons. We assume r 3 > r 2 > r 1 > 0 and r 4 < 0. r 5 and r 6 are two conjugated complex numbers. By using this form of U , we obtain As defined, The solution to Eq. (64) is Therefore, for the in-going null-rays, the coordinate time t increases from −∞ to +∞ as r decreases from +∞ to r 3 , decreases from +∞ to −∞ as r further decreases from r 3 to r 2 , increases from −∞ to +∞ as r continually decreases from r 2 to r 1 , and decreases finally from a finite value to −∞ as r decreases from r 1 to the black hole center. It should be noted that dt/dτ tends to +∞ for r → r 3 +0, to −∞ for r → r 2 +0 and to +∞ for r → r 1 + 0. Therefore, the radiation observed from infinity would appear infinitely red-shifted at the crossing of horizons r 3 and r 1 , and infinitely blue-shifted at the crossing of horizon r 2 .

B. General null geodesics
Turning to the consideration of Eq. (62) with L = 0 for the general null geodesics. We have where denotes the impact parameter. The property of orbits is determined by the number of roots of eight-degree equation, f (u) = 0. In general, the equation has four distinguished and positive roots, u 1 , u 2 , u 3 , u 4 . But for some cases, double roots can arise. For example, u 1 and u 2 coincide. In this case, we denote u 1 = u 2 with u 12 .
Taking into account of α 2 < 0, we find we must consider 14 cases as distinguished below: no real roots .
In Fig. (4), we present the plots of u vs f (u) for these cases. These different cases lead to the following conclusions. Case (1): The eight-degree equation, f (u) = 0, allows four positive real roots with 0 < u 1 < u 2 < u 3 < u 4 . We must distinguish between orbits of the two kinds: the orbits of the first kind restricted to the interval, 0 < u ≤ u 1 and u 2 ≤ u ≤ u 3 , and the second kind with u ≥ u 4 . The orbits in the interval, 0 < u ≤ u 1 , are the analogues of the hyperbolic orbits for massive particles in the Newtonian theory. The orbit in the interval, u 2 ≤ u ≤ u 3 , oscillates between two extreme values of u = u 2 and u = u 3 and they behave as the relativistic analogues of the Keplerian orbits for massive particles. The orbit of the second kind, starting at a certain aphelion distance, u = u 4 , plunges into the singularity.
Case (2): In this case, we have 0 < u 1 = u 2 = u 12 < u 3 < u 4 . The hyperbolic orbits and the oscillating orbits coalesce as they approach, asymptotically, a common circle from opposite sides by spiralling around it an infinite number of times. The orbit in the interval, u ≥ u 4 , starting at a certain aphelion distance, u = u 4 , plunges into the singularity.
Case (3): In this case, we have 0 < u 1 < u 2 = u 3 = u 23 < u 4 . The orbit at u = u 23 is a stable circular orbit with zero eccentricity. The orbits at intervals, 0 < u ≤ u 1 and u ≥ u 4 , are the same as the case (1).
Case (4): In this case, we have 0 < u 1 < u 2 < u 3 = u 4 = u 34 . The orbit of oscillation starts at a certain aphelion distance, u 2 , and approaches the circle of radius u −1 4 , asymptotically, by spiraling around it an infinite number of times. The orbit of the second kind spirals away from the same circle and eventually plunges into the central singularity. The orbit at the interval, 0 < u ≤ u 1 is the same that in the case (1).
Case (5): In this case, we have 0 < u 1 < u 2 with u 3 , u 4 a pair of complex-conjugate roots. The orbits belong to the hyperbolic and the plunging kind, respectively.
Case (6): In this case, we have 0 < u 3 < u 4 with u 1 , u 2 a pair of complex-conjugate roots. The orbits belong to the hyperbolic and the plunging kind, respectively.
Case (7): In this case, we have 0 < u 1 = u 2 = u 12 < u 3 = u 4 = u 34 . Then the hyperbolic orbits and the oscillating orbits coalesce as they approached, asymptotically, a common circle at radius, u −1 1 (= u −1 2 ), from opposite sides by spiralling round it an infinite number of times. On the other hand, the oscillating orbits and the plunging orbits also coalesce as they approached, asymptotically, a common circle at radius, u −1 3 (= u −1 4 ), from opposite sides by spiralling round it an infinite number of times.
Case (8): In this case, we have 0 < u 1 < u 4 with u 2 , u 3 a pair of complex-conjugate roots. The orbits belong to the hyperbolic and the plunging kind, respectively.
Case (9): In this case, we have 0 < u 1 < u 2 = u 3 = u 4 = u 234 . The orbits include the hyperbolic orbits and the plunging orbits. A remarkable fact of the plunging orbits is that they have an unstable circular orbit with radius, 1/u 234 .
Case (10): In this case, we have 0 < u 1 = u 2 = u 3 = u 123 < u 4 . The orbits include the hyperbolic orbits and the plunging orbits. A remarkable fact of the hyperbolic orbits is that they have an unstable circular orbit with radius, 1/u 123 .
Case (11): In this case, we have 0 < u 1 = u 2 = u 3 = u 4 = u 1234 . Then the hyperbolic orbits and the plunging orbits coalesce as they approached, asymptotically, a common circle at radius, u −1 1234 , from opposite sides by spiralling round it an infinite number of times.
Case (12): In this case, we have 0 < u 1 = u 2 = u 12 with u 3 , u 4 a pair of complex-conjugate roots. The hyperbolic orbits and the plunging orbits coalesce as they approached, asymptotically, a common circle at radius, u −1 12 , from opposite sides by spiralling round it an infinite number of times.
Case (13): In this case, we have 0 < u 3 = u 4 = u 34 with u 1 , u 2 a pair of complex-conjugate roots. The hyperbolic orbits and the plunging orbits coalesce as they approached, asymptotically, a common circle at radius, u −1 34 , from opposite sides by spiralling round it an infinite number of times.
Case (14): In this case, f (u) = 0 has no real roots or u 1 , u 2 a pair of complex-conjugate roots and u 3 , u 4 a pair of complex-conjugate roots. The resulting orbits can be considered as belonging to imaginary eccentricities with the remarkable fact that they are unbound orbits. Like the bound orbits, these unbound orbits similarly fall into the central singularity, but are allowed to start from infinity rather than from finite aphelion distance.

C. Stable circular orbits
Case (3) presents us a stable null circular orbit. We are interested in this orbit. The conditions for the occurrence of stable circular orbits are determined by where the comma denotes the derivative with respect to u. We note that these conditions are necessary but not sufficient. Observing the plot of case (3) in Fig. (4), we assign the roots as u 1 = 1, u 23 = 2, u 4 = 3. Then we have Solving these equations, we obtain With above parameters, the equation of horizon, U = 0, tells us there is only one event horizon with radius, r EH = 0.462. The radius of stable circular orbit is 1/u 23 = 0.5. It is outside of the event horizon. So we have a stable photon sphere outside the black hole.

D. Null innermost stable circular orbit
Different from Newtonian mechanics, the presence of an innermost stable circular orbit (ISCO) in General Relativity is a purely relativistic effect. The difference comes from the fact that the velocity of massive test particles cannot be equal or exceeding the speed of light. The ISCO represents the boundary between test particles orbiting the black hole and test particles falling into the black hole. As a result, it marks the inner edge of the accretion disk in the accretion disk model of Shakura and Sunjaev [34,35]. In accretion disk physics which can be compared to Event Horizon Telescopes observations, the model is frequently used as the starting point [36]. Therefore, it is valuable to compute the ISCO for 3-horizon black holes.
The radius of null ISCO, r SI = 1/u SI is determined by where prime denotes the derivative with respect to u. Now we have four equations, but five unknown functions. So we are left with one freedom, for example, u I . Thus we obtain Define the ratio for radii of null ISCO to event horizon, ǫ ≡ r I /r EH , we find from above equation which is larger than the Schwarzschild black hole (with ǫ = 3/2). The radius of null ISCO determines the boundary of the black hole shadow. Thus we conclude that the shadow of a 3-horizon black hole is larger than that of a Schwarzschild black hole although they have the same size of event horizon.

A. The radial geodesics
In this subsection, we shall consider the radial timelike geodesics with zero angular momentum. The equations governing these geodesics are Since U > 0 in the interval r 1 < r < r 2 , it is obvious that E 2 − U will vanish for some finite value of r 1 < r < r 2 provided that E 2 is smaller than the maximum of U for the interval r 1 < r < r 2 . We conclude that the trajectory will have a turning point in the interval r 1 < r < r 2 in this case. On the other hand, if E 2 is larger than the maximum of U for the interval r 1 < r < r 2 , the radial timelike trajectory would reach the singularity.
Case (2): In this case, u 1 and u 2 coincide and there exists stable circular orbit with radius 1/u 12 except for the oscillating and plunging orbits.
Case (3): In this case, u 2 and u 3 coincide and there exists unstable circular orbit with radius 1/u 23 . The two oscillating orbit coalesce as they approach, asymptotically, the common circle with radius 1/u 23 , from opposite sides by spiralling around it an infinite number of times.
Case (4): In this case, u 3 and u 4 coincide and there exists stable circular orbit with radius 1/u 34 except for the oscillating and plunging orbits.
Case (5): In this case, u 4 and u 5 coincide and there exists unstable circular orbit with radius 1/u 45 . The oscillating orbit and the plunging orbits coalesce as they approach, asymptotically, the common circle with radius 1/u 45 , from opposite sides by spiralling around it an infinite number of times.
Case (6): In this case, u 1 , u 2 and u 3 coincide and there exists unstable circular orbit with radius 1/u 123 except for the oscillating and plunging orbits.
Case (7): In this case, u 1 and u 2 coincide while u 3 and u 4 coincide. There exist two stable circular orbits with radius 1/u 12 and 1/u 34 , respectively, except for the plunging orbits.
Case (8): In this case, u 1 and u 2 coincide while u 4 and u 5 coincide. There exists a stable circular orbit with radius 1/u 12 and am unstable circular orbit with radius 1/u 45 , respectively. The oscillating orbit in the interval, u 3 ≤ u ≤ u 45 and the plunging orbits in the interval, u ≥ u 45 coalesce as they approach, asymptotically, the common circle with radius 1/u 45 , from opposite sides by spiralling around it an infinite number of times.
Case (9): In this case, u 2 , u 3 and u 4 coincide and there exists unstable circular orbit with radius 1/u 234 except for the oscillating and plunging orbits.
Case (10): In this case, u 2 and u 3 coincide while u 4 and u 5 coincide. There exist two unstable circular orbits with radius 1/u 23 and 1/u 45 , respectively, except for the oscillating and the plunging orbits.
Case (11): In this case, u 3 , u 4 and u 5 coincide and there exists an unstable circular orbit with radius 1/u 345 except for the oscillating orbit in the interval, u 1 ≤ u ≤ u 2 and the plunging orbits in the interval, u ≥ u 345 .
Case (12): In this case, there are three real roots, u 1 , u 2 and u 5 . u 3 and u 4 are a pair of complex-conjugate roots. So there exist the oscillating orbit in the interval, u 1 ≤ u ≤ u 2 and the plunging orbits in the interval, u ≥ u 5 .
Case (13): In this case, there are three real roots, u 3 , u 4 and u 5 . u 1 and u 2 are a pair of complex-conjugate roots. Similar to case (12), there exist the oscillating orbit in the interval, u 3 ≤ u ≤ u 4 and the plunging orbits in the interval, u ≥ u 5 .
Case (14): In this case, there are three real roots, u 1 , u 2 and u 3 . u 4 and u 5 are a pair of complex-conjugate roots. Similar to case (12), there exist the oscillating orbit in the interval, u 1 ≤ u ≤ u 2 and the plunging orbits in the interval, u ≥ u 3 .
Case (15): In this case, there are three real roots, u 1 , u 4 and u 5 . u 2 and u 3 are a pair of complex-conjugate roots. Similar to case (12), there exist the oscillating orbit in the interval, u 1 ≤ u ≤ u 4 and the plunging orbits in the interval, u ≥ u 5 .
Case (16): In this case, u 1 , u 2 , u 3 , u 4 coincide and there exists a stable circular orbit with radius 1/u 1234 except for the plunging orbits in the interval, u ≥ u 5 .
Case (17): In this case, u 1 , u 2 and u 3 coincide while u 4 and u 5 coincide. There exist two unstable circular orbits with radius 1/u 123 and 1/u 45 , respectively, except for the oscillating and the plunging orbits.
Case (18): In this case, u 1 and u 2 coincide while u 3 , u 4 and u 5 coincide. There exists a stable circular orbit with radius 1/u 12 and am unstable circular orbit with radius 1/u 345 , respectively, except for the plunging orbits.
Case (19): In this case, u 2 , u 3 , u 4 , u 5 coincide. There exists an unstable circular orbit with radius 1/u 2345 except for the oscillating and the plunging orbits.
Case (20): In this case, there are two real roots with double toots u 12 and u 5 . u 3 and u 4 are a pair of complexconjugate roots. There is a stable circular orbit with radius 1/u 12 and the plunging orbits in the interval, u ≥ u 5 .
Case (21): In this case, there are two real roots with double toots u 34 and u 5 . u 1 and u 2 are a pair of complexconjugate roots. Similar to case (20), there is a stable circular orbit with radius 1/u 34 and the plunging orbits in the interval, u ≥ u 5 .
Case (22): In this case, there are two real roots with double toots u 23 and u 1 . u 4 and u 5 are a pair of complexconjugate roots. The oscillating orbits stars at a certain aphelion distance, 1/u 1 , and approaches the circle of radius 1/u 23 ,asymptotically, by spiralling around it an infinite number of times. The orbit of the plunging is, in some sense, a continuation of the oscillating orbits in that it spirals away from the circle and then plunge eventually into the central singularity.
Case (23): In this case, there are two real roots with double toots u 45 and u 1 . u 2 and u 3 are a pair of complexconjugate roots. The properties of orbits are similar to the case (22).
Case (24): In this case, there is only one real root, u 123 or u 345 or u 12345 . There exists an unstable circular orbit with radius 1/u 123 , 1/u 345 or 1/u 12345 except for the plunging orbits.
Case (25): In this case, there is only one real root, u 1 or u 5 (the other four roots are complex-conjugate roots). We have only one class of orbits. Namely, all the orbits plunge into the singularity after starting from certain aphelion distances.
C. Stable timelike circular orbits Case (7) presents us two stable timelike circular orbits with radius 1/u 12 and 1/u 34 , respectively. We consider the orbits in this subsection. The conditions for the occurrence of stable circular orbits are determined by f (u) = 0 , and f ,u = 0 , where the comma denotes the derivative with respect to u. We assign the roots as u 12 = 1, u 34 = 2, u 5 = 3.
Then we have Solving these equations, we obtain With above parameters, the equation of horizon, U = 0, tells us there is only one event horizon with radius, r EH = 0.264. The radius of two stable circular orbits are 1/u 12 = 1 and 1/u 34 = 0.5, respectively. They are all larger than the event horizon. So we have two stable matter spheres outside the black hole.

D. Timelike innermost stable circular orbit
The radius of timelike ISCO, r SI = 1/u SI is determined by where prime denotes the derivative with respect to u. We have five equations, but six unknown functions. We choose u I as the remaining freedom. Thus we obtain Notice that the energy of the particles is independent on the inverse radius u I of ISCO. In order to compare the ISCO with black hole event horizon, we insert Eqs. (90),(91) and (92) into the horizon equation and eventually find which is smaller than the Schwarzschild black hole (with ǫ = 3). This means the inner edge of accretion disk for the 3-horizon black holes is pushed more closer to the event horizon compared with the Schwarzschild black hole.

VII. THE MOTION OF CHARGED PARTICLES
The motion of a test particle with net charge is determined by the Lagrangian where q denotes the charge per unit mass of the test particle. The equations of motion in the equatorial plane following from this Lagrangian are and in place of Eq. (62) we obtain du dφ We can do orbital analysis as in previous sections. But we will not do it here. One novel feature of Eq. (98) is that when the test particle has a turning point at it arrives at the event horizon, its energy will be (100) In the absence of α 2 , the energy is negative if and only if the charges q and Q have different signs. However, in the presence of α 2 , the energy can be negative even if the charges q and Q have the same signs provided that α 2 is sufficiently negative. The negative of energy E leads to the Penrose process that one can extract energy from the black hole by using charged particles.

NUMBERS
In this section, we shall study the static response of the 3-horizon black holes to external scalar field. The quantities of response are intrinsic and contribute to the form of gravitational waves. Therefore, they can be in principle discovered [37] in the observation of gravitational waves. The response of an object to a long-wavelength tidal field is encoded in the so-called Love numbers [38], which describe the deformability or rigidity of the object. The Love numbers of Schwarzschild black hole are exactly vanishing in four dimensional spacetime [39][40][41][42][43][44][45]. This reveals the Schwarzschild black hole is totally rigid. However, it is not case when one consider higher dimensions [42,47], anti-de Sitter asymptotical [48], the presence of higher-curvature terms [37], or many alternative theories of gravity [49]. It is found that the Love numbers of those cases are non-vanishing. Here we compute the Love numbers of 3-horizon black hole caused by scalar tidal field. The 3-horizon black hole can be thought as the extension of RN black hole and to our knowledge, one did not compute the Love numbers for RN black holes. So in the next, we shall begin from the studying of RN black holes.

Reissner-Nordstrom black holes
The extremal RN metric is where ∆ = r 2 − 2M r + M 2 with M the mass of the black hole. This spacetime has only one event horizon with the radius r EH = M . A static, massless scalar field in this background satisfies the equation of motion [45,46] where φ l is the radial component of the scalar field Φ which has been decomposed as Φ = φ l (r)Y lm (θ, φ). Y lm (θ, φ) is the spherical harmonic function. l is the angular quantum number. Make the change of variables, r → z as follows namely, The region of r ∈ [r EH , +∞) is mapped into z ∈ [+∞, 0). Then the equation of motion becomes Solving this equation, we obtain or where c 1 and c 2 are integration constants. Now we require two boundary conditions to specify completely the solution. The first boundary condition is that the scalar field is finite at the event horizon. The second boundary condition is to fix the normalization of the growing mode solution at spatial infinity. Then one can read off the induced sub-leading fall-off at spatial infinity, which plays the role of the linear response to the external field. The Love number, describing the response of the black hole to external perturbations, is defined as the ratio between the coefficients of the decaying and growing modes at infinity. The first boundary condition reveals c 1 = 0. The second boundary condition tells us we have φ l = c 2 r l in spatial infinity. We can understand this point by imagining that we are applying a scalar field, in which the black hole is immersed, that scales like r l as r → +∞ with angular structure given by l = 1, 2, 3, · · · harmonics. Therefore, the long-wavelength external scalar field one can apply is This solution is a pure growing mode and the sub-leading induced fall-off mode does not exist. This means the Love numbers of extreme RN black holes are vanishing.
B. The Love numbers of non-extremal

Reissner-Nordstrom black holes
In this case, we have ∆ = (r − r + )(r − r 1 ) with r ± = M ± M 2 − Q 2 . Make the change of variables, r → z as follows namely, The region of r ∈ [r + , +∞) is mapped into z ∈ [+∞, 0). The equation of motion becomes Solving this equation, we obtain where L is defined by c 1 and c 2 are integration constants. We have z = +∞ on the event horizon r = r + . When z → +∞, we have The first boundary tells us c 2 term must be dropped. Therefore, the long-wavelength external scalar field can be applied is It is a growing mode with respect to r and the decaying mode does not exist. Thus the Love number of nonextremal RN black hole is vanishing.
Make the change of variables, r → z as follows Then the region of r ∈ [r 3 , +∞) is mapped into z ∈ [+∞, 0). The equation of motion becomes Here ∆ is understood as the function of z. Now we look for the fall-off solution at r → +∞, or z → 0. When r → +∞, we have ∆ = 1/z 2 and r = 1/z. The equation of motion, Eq. (117), reduces to exactly that for the extremal RN black hole. The corresponding Love number is zero. On the other hand, when r → r 3 , or z → +∞, we have ∆ = c 3 (r − r 5 ) with c 3 a positive constant and r = r 3 + e −c3z . The resulting equation of motion is The solution is The c 2 term is divergent on the horizon. Therefore it should be dropped. We eventually find that It is a growing mode with φ L = c 1 when r = r 3 . There is no the decaying fall-off mode. In other words, the Love number of 3-horizon black hole is also vanishing.

IX. THERMODYNAMICS
Finally, we make an investigation of the thermodynamics for the 3-horizon black hole. Concretely, we shall derive the Smarr formula and the first law of thermodynamics. It is worth noting that many studies in NLE in this respect have been carried out, for example in [50][51][52][53][54][55]. Hawking showed that for the outermost, event horizon in an asymptotically flat spacetime, the temperature of black hole is where the surface gravity κ EH is defined by evaluating on the event horizon. Here l µ is the future-directed null generator of the event horizon, which coincides with a Killing vector K µ on the horizon. The metric for the 3-horizon black holes is static. So if the Killing vector is adopted as K µ = ∂/∂t, then we have As a result, the black hole temperature takes the form The radius r EH of event horizon is determined by Thus the temperature is eventually determined by the mass M , charge Q and coupling constant α 2 . The entropy of black holes generally satisfies the area law which states that the entropy is a quarter of the area of black hole event horizon [56][57][58]. Therefore we have the entropy of the black hole The electrostatic potential on the event horizon is We find if we define the pressure P and the thermodynamic volume V as follows, respectively, V ≡ ∂M ∂P S,Q = − α 2 2 Q 4 5r 5

EH
. (129) Then the Smarr formula is satisfied. We could make an examination on whether the thermal quantities fulfill the requirement of the first law of thermodynamics. To this end, we treat the mass M , the entropy S and the pressure P as the function of r EH , Q, P . Then we have dM = M ,rEH dr EH + M ,Q dQ + M ,P dP , dS = S ,rEH dr EH .
After computation, we find the first law of thermodynamics is indeed satisfied.

X. CONCLUSION AND DISCUSSION
In summary, starting from the NLE Lagrangian with infinite series of Maxwell invariant and using the method of infinite series, we find the black hole solutions with many horizons. To be specific, we present the solutions for 3-horizon and 4-horizon black holes. In particular, the explicit and analytic expression for NLE Lagrangian of the 3-horizon black hole is obtained. To our knowledge, one did not yet get 3-horizon black holes with an analytic NLE Lagrangian. On the other hand, our lagrangians have intrinsics advantage over previous studies [25,26]. They are not dependent on the mass M , charge Q, or the ratio between charge to mass Q/M . There are only related to the coupling constants, α i . We find that there are three physical parameters, the mass M , charge Q and the coupling constant α 2 for the 3-horizon black hole. For negative coupling constant α 2 , the spacetime can have three horizons. When the coupling constant α 2 vanishes, it reduces to the RN spacetime. What is more, the charge Q can be much larger than the mass M while the central singularity remains dressed by an event horizon. For black holes with N + 1 horizons, there are N + 1 physical parameters, the mass M , charge Q and the coupling constants α i with i running over from 2 to N . In this sense, the black holes have N + 1 "hairs".
We find that, for the multi-horizon black holes, both the null geodesics and the timelike geodesics are considerable rich. For null geodesics, there are: (1) the hyperbolic orbits; (2) the oscillating orbits; (3) the plunging orbits; (4) the stable circular orbit; (5) the unstable circular orbit; (6) the hyperbolic and plunging coalescing orbits; (7) the hyperbolic and the oscillating coalescing orbits; (8) the plunging and the oscillating coalescing orbits; (9) the hyperbolic, the oscillating and the plunging coalescing orbits. It is found there is one stable circular orbit for the null geodesics. The ratio of radii for null ISCO to event horizon is ǫ ≃ 1.5055 which is larger than the Schwarzschild black hole (with ǫ ≃ 1.5). Above conclusion are also applicable to timelike geodesics when E 2 > 1. When E 2 < 1 (the bound orbits), the structure of timelike geodesics is richer than that of null geodesics. For example, one can have two oscillating orbits in case (1) and two stable circular orbits in case (7). The ratio of radii for timelike ISCO to event horizon is ǫ ≃ 2.0678 which is smaller than the Schwarzschild black hole (with ǫ ≃ 3).
For the RN black hole, the Penrose process occurs if and only if the charge of test particle q has different sign from the black hole charge Q. However, it is not the case for 3-horizon black hole. It is found the Penrose process can also occur when the two charges, q and Q have the same signs provided that α 2 is sufficiently negative. This is an interesting property for multi-horizon black holes.
We also calculate the Love numbers of RN black holes and 3-horizon black holes. We find they are all vanishing. This reveals both the RN black holes and the 3-horizon black holes are rigid. It is the same as the Schwarzschild black holes. Finally, the thermodynamics for 3-horizon black holes is developed. It is found that the inverse of coupling constant α 2 plays the role of thermal pressure.