Quasi-Normal Modes and Stability of Einstein-Born-Infeld Black Holes in de Sitter Space

We study gravitational perturbations of electrically charged black holes in (3+1)-dimensional Einstein-Born-Infeld gravity with a positive cosmological constant. For the axial perturbations, we obtain a set of decoupled Schrodinger-type equations, whose formal expressions, in terms of metric functions, are the same as those without cosmological constant, corresponding to the Regge-Wheeler equation in the proper limit. We compute the quasi-normal modes (QNMs) of the decoupled perturbations using the Schutz-Iyer-Will's WKB method. We discuss the stability of the charged black holes by investigating the dependence of quasi-normal frequencies on the parameters of the theory, correcting some errors in the literature. It is found that all the axial perturbations are stable for the cases where the WKB method applies. There are cases where the conventional WKB method does not apply, like the three-turning-points problem, so that a more generalized formalism is necessary for studying their QNMs and stabilities. We find that, for the degenerate horizons with the"point-like"horizons at the origin, the QNMs are quite long-lived, close to the quasi-resonance modes, in addition to the"frozen"QNMs for the Nariai-type horizons and the usual (short-lived) QNMs for the extremal black hole horizons. This is a genuine effect of the branch which does not have the general relativity limit. We also study the exact solution near the (charged) Nariai limit and find good agreements even far beyond the limit for the imaginary frequency parts.


I. INTRODUCTION
The organization of this paper is as follow. In Sec. 2, we review the electrically charged black hole solution in EBI gravity with a positive cosmological constant. In Sec. 3, we consider the axial perturbations from the spherically symmetric EBI black hole background and obtain a set of decoupled Regge-Wheeler equations. In Sec. 4, we use the Shutz-Iyer-Will's third-order WKB method to compute the QNMs numerically and discuss the stabilities of the EBI black holes. In Sec. 5, we study the exact solution near the (charged) Nariai limit and compare it with the numerical results in Sec. 4. In Sec. 6, we conclude with several discussions. Throughout this paper, we use the conventional units for the speed of light c, the electric and magnetic constants ǫ 0 , µ 0 , and the Boltzman constant k B , c = 4µ 0 = 4/ǫ 0 = k B = 1, but keep the Newton constant G and the Planck constanth unless stated otherwise.

II. BACKGROUND SOLUTIONS
In this section, we briefly review the black hole solution in EBI gravity with a positive cosmological constant Λ in (3 + 1) dimensions, from which the gravitational perturbations will be considered. The EBI action is given by where the BI Lagrangian density L(F ) is given by Here, β is the BI's coupling constant with dimensions [length] −2 [5,6]. In the weak field or equivalently strong coupling limit, |F µν F µν | ≪ 2β 2 , L(F ) may be expanded as where F 2 ≡ F µν F µν and the usual Maxwell's electrodynamics is recovered at the lowest order. The correction terms, when F is comparable to β, represent the effect of the nonlinear BI fields at the short distance ∼ β −1 and the possible electromagnetic field strength is bounded by −F 2 ≤ 2β 2 . Taking 16πG = 1 for simplicity, the equations of motion are obtained as where the energy-momentum tensor for BI fields is given by Let us now consider a static and spherically symmetric solution with the metric ansatz ds 2 = −N 2 (r)dt 2 + 1 f (r) dr 2 + r 2 (dθ 2 + sin 2 θdφ 2 ).
For the static, electrically-charged case where the only non-vanishing component of the field strength tensor is F rt ≡ E r , the general solutions for the metric and the BI electric field are obtained as in terms of the hypergeometric function [12,14]. Here Q is the electric charge 1 and M is the ADM mass which is composed of the intrinsic mass C and (finite) self energy of a point charge M 0 , defined by In the large-distance limit r ≫ Q/β, the solutions become which show the usual Reissner-Nordstrom-(Anti) de Sitter (RN-(A)dS) black hole at the leading order and the short-distance corrections at the sub-leading orders. On the other hand, in the short-distance limit r ≪ Q/β, we have which show the regularized solutions near the origin, with the milder (curvature) singularity of the metric and the finite electric field at the origin [14,15] due to the non-linear BI fields for a finite coupling β.
The solution (8) can have three horizons generally, i.e., two (inner, outer) black hole horizons, r − , r + , and one cosmological horizon r ++ . The Hawking temperature for the outer The plots of the extremal horizons r * H and r * C (> r * H ) vs. βQ for varying β. We consider β = 2, 1, 1/2, 1/3 (right to left) with Λ = 0.2. The solid/dashed lines denote the '+/−' roots in (17) and there is no '+' root for the last case. black horizon r + is given by and plotted in Fig. 1. There exist two extreme limits of vanishing temperature, when the outer black hole horizon r + meets the inner horizon r − (extremal black holes) or the cosmological horizon r ++ (Nariai solution), at for a positive cosmological constant Λ (r * H < r * C ) (Figs. 1 and 2). At the extreme points, the ADM mass M, becomes [15,16] The first law of black hole thermodynamics is found as with the usual Bekenstein-Hawking formula for the black hole entropy and the scalar potential [12][13][14] A 0 (r) = Q r 2 F 1 The integrated first-law of black hole thermodynamics, called the generalized Smarr relation, is given by [13,14,17,18] 2 The generalized Smarr relation without the last term has been obtained for the planar black holes in [14,18]. If we consider topological EBI black holes with the topological parameter k = +1, 0, −1 for spherical, plane, and hyperbolic geometry generally, with the solid angle Ω k , one finds [13] the last term as (k/3) 4hS BH /Ω k G. The plots of the ADM mass M vs. the black hole horizon radius r + for varying βQ with a fixed Λ > 0. The marginal mass M 0 is given by the mass value at r + = 0. The top two curves in the left represent M 0 > M * , M 0 = M * for βQ > 1/2, βQ = 1/2, respectively, with the extremal mass M * , whereas the bottom two curves represent the cases where M * is absent for βQ < 1/2. We consider βQ = 2/3, 1/2, 2/4.5, 1/3 (top to bottom) with β = 2 and Λ = 0.2. The effect of cosmological constant is not significant for small r + (left) but important for large r + (right). In the latter, we compare the dS case of Λ = 0.2 (thick curve) with the flat Λ = 0 (medium curve) and the AdS case of Λ = −0.2 (thin curve), respectively. Now, from the relation (17), one can classify the black holes by the values of βQ [15].
(i) βQ > 1/2: In this case, there are generally three horizons (two smaller ones for black holes and the largest one for the cosmological horizon) for M * H ≤ M ≤ M * C (< M 0 ), depending on M and Λ (Fig. 3, left), where M * H and M * C denote the values of the extremal mass M * in (19) at r * H and r * C for the extremal black holes and the (charged) Nariai solutions, respectively, with the vanishing Hawking temperature (Fig.1). When the mass is outside this range, i.e., M < M * H or M > M * C , the singularity at r = 0 becomes naked as in the RN black hole.
(ii) βQ = 1/2: In this case, only the Schwarzschild -de Sitter (Sch-dS) -like black holes with a non-degenerate black hole horizon r + are possible for M 0 < M ≤ M * C (Fig. 3, center). It is peculiar that the black hole horizon r + shrinks to zero size, i.e., the point-like horizon, for the marginal case M = M 0 (Fig. 4). This corresponds to the extremal black hole, with the vanishing Hawking temperature (Fig.1), where the horizon degenerates at the origin. When the mass M is smaller than the marginal mass M 0 or larger than M * C , the singularity at r = 0 becomes naked. Note that, in this case, the GR limit β → ∞ does not exist.
(iii) βQ < 1/2: This case is similar to the case (ii), except that the marginal case M = M 0 has no (even a point) horizon so that its curvature singularity at r = 0 is naked always (Fig.  3, right), even though the (non-degenerate) black hole horizon can arbitrarily approach close to the point-like horizon r + = 0 in the limit M → M 0 , with the divergent Hawking temperature (Fig. 1) 3 .

III. PERTURBATION EQUATIONS
In this section, we consider the gravitational perturbations of an electrically charged black hole in EBI gravity with a positive cosmological constant Λ > 0, whose solutions for the metric and fields are given by (8) and (9), respectively. We follow the procedure of Chandrasekhar [19] for the perturbations of RN solution but in the different conventions which agree with those of Wald [20] 4 . We generalize the earlier computations of [21] to a positive cosmological constant but with some important corrections of errors. To accommodate the procedure of [19], it is useful to write the metric (7) as , e 2µ 3 ≡ r 2 , e 2ψ ≡ r 2 sin 2 θ.
We will consider (24) as a background metric and consider the first-order metric perturbations. In our linear perturbations of the spherically symmetric system, we may restrict ourselves to axisymmetric modes of perturbations, without loss of generality [19], whose metric may be generally written as where the metric functions ν, ψ, µ 2 , µ 3 , q 0 , q 2 , and q 3 are functions of only t, r, and θ, due to axisymmetry. There are two kinds of perturbations: The perturbations of q 0 , q 2 , and q 3 , which are called "axial" perturbations, induce a dragging effect and impart a rotation to the black hole, whereas the increments of δν, δψ, δµ 2 , and δµ 3 , which are called "polar" perturbations, do not impart such rotation. These two perturbations are decoupled and can be treated independently [19]. In this paper, we will focus only on the axial perturbations for simplicity. To obtain the solutions for the perturbed metric (26), we first need to consider the perturbations of BI and Einstein equations, (4) and (5). For this purpose, we adopt the tetrad formalism as in [19] and use the Roman indices a, b = 0, 1, 2, 3 for the tetrad indices and the Greek indices µ, ν = t, φ, r, θ for the curved coordinate indices. The tetrad basis e a µ satisfies e a µ e b ν η ab = g µν with η ab = diag(−1, +1, +1, +1). For the metric given by (26), e a µ can be obtained as its mass M has still a non-vanishing remnant M 0 . If this tiny black hole with a unit charge Q = √ 4πǫ 0 αhc (α is the fine-structure constant) is created at the Planck energy scale in the early universe, i.e., M 0 ≈ hc/G, we obtain β ≈ (9/2)Γ(3/4) 4 π −7/2 α −3/2 (ǫ 0h c) −1/2 G −1 . If we consider α P l ∼ 1.5 at the early universe, contrast to α ∼ 10 −2 at the current epoch, we have β P l ∼ 0.1 (ǫ 0h c) −1/2 G −1 and β P l Q ≈ 9Γ(3/4) 4 π −3 α −1 P l G −1 ∼ 0.4 G −1 . 4 We use a metric with signature (− + ++) and the Wald's conventions for the Riemann and Ricci tensors which differ from [19,21,22].
Tensors in the two bases are related by the tetrads, for example, for the field strength and the Ricci tensor. We will now consider equations of motion (4) and (5) in the tetrad basis.
Here, (37), (38), and (39) represent the Gauss's law and the Ampere's law (with the Maxwell's correction term) in the curved space. On the other hand, (40), which corresponds the φ-component of the Ampere's law, shows the induced source terms which correspond the electric currents J (e) ∼ F ab Q ab , due to non-linear interactions between the electromagnetic and gravitational fields, similar to (33). These induced source terms in (40) and (33) represent the back-reactions from the dynamical metric q 0 , q 2 and q 3 , which are also driven by T µν in the Einstein equations (5).
Up to now, we have not assumed any specific solution of the background BI fields. Now, using that the only non-zero components of F and D in the background solution (9) are and substituting ψ, ν, µ 2 , µ 3 with ψ + δψ, ν + δν, µ 2 + δµ 2 , µ 3 + δµ 3 in the above equations, we find the linearized versions of perturbed BI equations as

The perturbed Ricci tensor equations
To find the metric perturbations of the Einstein equation (5), we may conveniently write the Einstein equations (5) as 5 With a simple computation, one can easily find that the cosmological constant term does not appear explicitly in the perturbation of Einstein equations (48). The lineralized form of the perturbed Ricci tensor equations are then given by Since the only non-zero componsnt of the background F ab is F 02 , we replace δF ab by F ab except for δF 02 for simplification, following [19]. Explicit forms of the right-hand side of the perturbed Ricci tensors (51) are given by C. The wave equations for the axial perturbations The linearized Einstein equations can be obtained by equating δR ab of the left-hand side in the above equations with those which can be computed from the Ricci tensors for the metric (26) in the tetrad basis, whose explicit forms are given in the Appendix. As mentioned before, the perturbation equations can be categorized into two kinds. One is the axial perturbations which are characterized by the non-vanishing of q 0 , q 2 , q 3 , F 01 , F 12 , and F 13 . The other is the polar perturbations which are characterized by the non-vanishing of δF 02 , F 03 , F 23 , δν, δψ, δµ 2 , and δµ 3 . In this paper, we only consider the axial perturbations, which are relatively easier and can be treated independently of the polar perturbations. For this purpose, we first consider (57) and (58), which are now given by Eliminating D 12 and D 13 from (44) using (42) and (43), we obtain where With the substitution (60) and (61) take the forms where ∆ ≡ r 2 e 2ν . Assuming that the perturbed fields q 0 , q 2 , q 3 , and Q have a time dependence e −iωt and eliminating q 0 from (66) and (67), we obtain Similarly, eliminating (q 0,2 − q 2,0 ) ,t from (62) using (66), we obtain We can further separate the variables r and θ in (68) and (69) with and where C m n is the Gegenbauer function and we have used the recurrence relation, Substituting (70) and (71) in (68) and (69), we obtain two radial wave equations, where and l = 0, 1, 2, · · · denotes the orbital number. Introducing the tortoise coordinate r * , and further defining Q and D as we find a pair of coupled equations of H 1 and H 2 [21] as To find the one-dimensional Schrödinger-type wave equations, we rewrite the coupled equations (79) and (80) in the matrix form, where The matrix V ij can be diagonalized by the similarity transformation as in the RN case [19]. Then, the two decoupled, one-dimensional Schrödinger-type equations are given by, where the effective potentials, U 1 and U 2 , which are real-valued functions, are given by Here, the decoupled solutions Z 1 and Z 2 are related with the coupled solutions H 1 and H 2 , which are basically the perturbations of the metric variable Q 23 and the fields strength D 01 [19], where The explicit form of U 1 , U 2 , q 1 , and q 2 will be available upon substituting e 2ν and e ϕ with the background solutions (25) and (78), respectively, when we compute the QNMs in the next section.

IV. QUASI-NORMAL MODES FOR THE AXIAL PERTURBATIONS: THE WKB APPROACH
A. The WKB approximation Quasi-Normal Modes (QNMs) are defined as the solutions which are purely ingoing as e −i(ωt+kr * ) near the (outer) black hole horizon r + , where r * = −∞. In this paper, in order to compute QNMs, we will consider the WKB approximation method [9][10][11], which is simple to apply in our case though there are several inherent limitations, as will be explained below. The applicability depends also on several conditions, like the (asymptotic) boundary conditions, and our dS space-time satisfies this condition which is a technical reason of using the WKB approach in this paper. In this section, we briefly summarize the basic results of the WKB approach.
The master equation for the WKB approach may be written as the one-dimensional Schödinger-type equation 6 , where the potential function −Q(x) is assumed to be constant at the asymptotic boundaries x = ±∞, although not necessary the same at both ends, and has one maximum at a certain finite x 0 (Fig. 5). In the black hole case, Φ represents the radial part of the perturbation with the usual time dependence of a positive-frequency mode e −iωt , as well as the appropriate angular dependence. The coordinate x is related to the tortoise coordinate r * which ranges from −∞, at the (outer) black hole horizon r + , to +∞, at the spatial infinity r = ∞ for the FIG. 5: A typical plot of the real part of the potential function −Q(x) ≡ U − ω 2 with the effective potential U , which is independent of the frequency ω ≡ ω R −iω I , and the (complex-valued) "energy" and two turning points at x 1 and x 2 are assumed.
asymptotically flat [10] or at the cosmological horizon r ++ for the asymptotically dS case [22], as adopted in this paper.
In order to solve the one-dimensional potential barrier problem of (92), we adopt a modified WKB approach by matching the two exterior WKB solutions in regions I and III across the two turning points x 1 and x 2 simultaneously [23]. In the interior region II, we expand the potential function −Q(x) near its maximum at x = x 0 and use it to connect the two exterior WKB solutions. For the black hole case with the quasi-normal mode boundary conditions, i.e., purely "outgoing" (from the potential barrier), the result may be written as the simple semi-analytic formula, where the primes denote the derivatives in terms of x and Q ′′ 0 > 0 due to the extremum of the potential −Q at x 0 . Here, the subscripts 2k, 2k + 1 denote the order of the WKB approximations, N = N or N ±1, and Λ 2k , Ω 2k+1 are polynomials of the derivatives Q (m) It is important to note that even and odd-order terms become pure-imaginary and real-valued, respectively, so that, when considering Q 0 = ω 2 − U 0 in the black hole case, they play the distinct roles in obtaining the quasi-normal frequencies, ω = ω R − iω I 7 . Then, the result up to the third WKB orders [10,11] (86), is given by 7 This fact does not seem to be well emphasized in the literature. This causes some misleading results, for example, in [21] due to the related typos in the original work of the third-order terms in [10] where "i" factor in the even part of (93) is missing. Similar typos also appear in the 4th-order terms in [24] (Appendix A), in contrast to the correct i factors in the corresponding 4th-order terms in [10] (Appendix) and [25] (private communications: We thank the authors for kindly sharing their results). 8 The convergence in each order may not be guaranteed generally since the expansions are just asymptotic ones. However, the popular third-order result of [10] can also be obtained as the lower-order results of the phase integral approach which allows the convergent expansions with controlled accuracy [26] and works with the real-valued WKB correction terms, where α ≡ n+ 1 2 (n = 0, 1, 2, · · ·), Ω 3 ≡ Ω 3 /(n+1/2), and the primes and the superscript (m) of U (m) 0 denote the derivatives of the effective potential U(r * ) with respect to the tortoise coordinate r * , evaluated at its maximum point r 0 . The maximum point r 0 is determined by dU /dr * = e 2ν dU/dr = 0, which does not have analytic solutions in our case, where the effective potentials U i in (87) and (88) are quite complicated so that r 0 can be found only numerically. Then, QNMs can be found numerically by solving (94) from the numerical computations of (95) and (96) as a function of the parameters of the theory.  QNMs with n = 0 as functions of the BI parameter β. As can be seen in Fig. 6 (left), there are one degenerate event horizon for β = β min , and one degenerate black hole horizon with one cosmological horizon for β = β max . In between them, β min < β < β max , there are generally three horizons (two black hole horizons and one cosmological horizon), depending on β. Only for this dS black hole regime, the WKB formula (94) for QNMs may apply. Fig. well in many cases. For extensions to higher orders, see [10] (Appendix) (4th order), [24] (6th order), and [25] (13th order). See also [27] for a recent review.  7 shows that the black holes are stable for the range of β min < β ≤ β max since their metric perturbations are decaying (which is called also as the"ring-down" phase) with the decay time τ ∼ ω −1

I
for ω I > 0. Moreover, the results show ω ≈ 0 9 , i.e., the "frozen" mode, at β = β min ≈ 0.103 and no QNMs beyond β max ≈ 1.166, where the WKB formula (94) does not apply, as expected from the behaviors of the metric functions and the effective potentials in Fig. 6. As β increases, the negative imaginary parts of QNMs, ω I increase first monotonically and reach maximum points around βQ = 1/2 (here, β = 0.5 with Q = 1), which divides the Sch-type and the RN-type [15], and then decrease monotonically. Beyond the above range, β < β min or β > β max , one can not say anything about its stability since the WKB formula for QNMs can not apply due to the lack of the required boundary conditions and we need other analysis to test its stability. For example, the WKB formula does not say anything about the stability of the naked singularity for β > β max . On the other hand, the real parts of QNMs ω R increase monotonically and saturate to maximum values, corresponding to those of the RN-dS case [22,34], showing two decoupled oscillation modes with a rough relation, ω R(1) ∼ 2ω R (2) . This is in contrast to the imaginary parts with a rough relation, ω I(1) ∼ ω I(2) , so that the real-to-imaginary frequency ratios are ω R(1) /ω I(1) ∼ 2ω R(2) /ω I(2) ∼ 10.   show the effective potentials and corresponding lowest QNMs as functions of the electric charge Q. As can be seen in Fig. 8, there are one degenerate event horizon (charged Nariai solution) with one Cauchy horizon for Q = Q min , and one degenerate black hole horizon with one cosmological horizon for Q = Q max . In between them, Q min < Q < Q max , there are generally three (two black hole and one cosmological) horizons, depending on Q. The results show ω ≈ 0 at Q = Q min ≈ 0.9320 and no QNMs beyond Q max ≈ 1.0041, as expected from the behaviors of the metric functions and the effective potentials in Fig.  8. The behaviors of ω I and ω R are almost the same as those of Fig. 7 for β < 0.5 with the similar relations, ω R(1) ∼ 2ω R(2) , ω I(1) ∼ ω I (2) and ω R(1) /ω I(1) ∼ 2ω R(2) /ω I(2) ∼ 10, and other discussions are similar. Here, we divide into three cases depending on the values of βQ, which divides Sch-type and RN-type by βQ = 1/2 [15]. First of all, Figs. 10 and 11 show the effective potentials and corresponding lowest QNMs as functions of the mass M for a RN-type black hole with fixed βQ > 1/2. As can be seen in Fig. 10, there are one degenerate event horizon (chraged Nariai solution) with one Cauchy horizon for M = M max , and one degenerate black hole horizon       However, the behavior of ω R is similar to the case of βQ > 1/2 in Fig. 10, except the rapid but finite oscillation near M min .
Figs. 14 and 15 show the cases for a Sch-type black hole with βQ < 1/2. The basic difference from the case of βQ = 1/2 in Figs. 12 and 13 is that there is no (even a point) black hole horizon in addition to the cosmological horizon for M = M min , even though the (non-degenerate) black hole horizon can arbitrarily approach close to r + = 0 as M → M 0 , with the divergent Hawking temperature (Fig. 1). One more remarkable difference is that, as M → M min , ω is divergent as  Fig. 15 shows ω ≈ 0 at M max ≈ 0.9774, similarly to the case of βQ > 1/2 in Fig. 11.   10 In the numerical computations of WKB formula (93), we find that the "smallness" of U ′ 0 , which should be "0" by definition, can be a good barometer of the numerical errors. As M → M min , the (non-vanishing) magnitude of U ′ 0 is increasing as shown in Fig. 17. The plots of f = e 2ν(r) and the effective potentials U 1 (r) and U 2 (r) for varying the orbital number l with a fixed βQ = 1/2. Here, we consider from l = 0 to l = 5 (bottom to top, barrier region) when M = 0.95, β = 1/2, Q = 1, and Λ = 0.2 with two (one for black hole and one for cosmological) horizons.

E. QNMs vs. l
We also divide into three cases depending on the values of βQ. Figs. 18 and 19 show the effective potentials and corresponding lowest QNMs as functions of the orbital number l for a fixed βQ > 1/2. As shown in Fig. 18, there are three (two for black hole and one for cosmological) horizons independent of l ( Fig. 18 (left)). On the other hand, the effective potentials depend on l, representing the angular momentum barriers for higher l with a single maximum. However, for lower l, there is no local maximum (but a minimum) in the  effective potential (l = 0 case for U 2 ) or there are "three turning points" with one additional turning point within the usual two turning points, i.e., the outer black hole horizon at r + ≈ 1.2 and the cosmological horizon at r ++ ≈ 2.4 (l = 0 case for U 1 and l = 1 case for U 2 ) so that the usual WKB formula does not apply 11 . One remarkable thing in the result of Fig. 19 is that, in contrast to other varying parameters, ω I(1) and ω I(2) respond differently to the varying l, while ω R(1) and ω R(2) respond almost identically. In particular, in the large l limit, ω I approaches asymptotically to a limiting value ω I ≈ 0.03055, while ω R shows a linear dependence ω R = σl with σ ≈ 0.1092, similar to earlier results [11]. Figs. 20 and 21 show the case for βQ = 1/2. The main difference from βQ > 1/2 case in Figs. 18 and 19 is the "magnitude flip" between ω I(1) and ω I(2) in Fig. 21. In the large l limit, ω I approaches asymptotically to a limiting value ω I ≈ 0.03349, while ω R shows a linear dependence ω R = σl with σ ≈ 0.1006. On the other hand, for lower l, i.e., l = 0 for U 1 and l = 0, 1 for U 2 , the usual WKB formula does not apply, due to either "three turning points" (l = 0 for U 1 and l = 1 for U 2 ) or "no local maximum" (l = 0 for U 2 ) within the usual two turning points r + ≈ 1.4 and r ++ ≈ 2.4.
Figs. 22 and 23 show the cases for βQ < 1/2 but the results do not show any qualitative difference from βQ = 1/2 in Figs. 20 and 21. In the large l limit, ω I approaches asymptotically to a limiting value, ω I ≈ 0.03345, while ω R shows a linear dependence ω R = σl with σ ≈ 0.0910. For lower l, i.e., l = 0 for U 1 and l = 0, 1 for U 2 , the usual WKB formula does not apply, due to either "three turning points" (l = 0 for U 1 and l = 1 for U 2 ) or "no local maximum" (l = 0 for U 2 ) within the usual two turning points r + ≈ 1.5 and r ++ ≈ 2.4.    there are one degenerate event horizon (charged Nariai solution) with one Cauchy horizon for Λ = Λ max , and one degenerate black hole horizon with one cosmological horizon for Λ min . In between them, Λ min < Λ < Λ max , there are generally three horizons depending on Λ, similar to Fig. 10. Fig. 25 shows ω ≈ 0 at Λ max ≈ 0.2352 and no QNMs below Λ min ≈ 0.1796, similar to the case in Fig. 11. Figs. 26 and 27 show the case for βQ = 1/2. As can be seen in Fig. 26, there are one degenerate event horizon for Λ = Λ max , and one non-degenerate black hole horizon for Λ = "0", i.e., the flat case. In between them, 0 < Λ < Λ max , there are generally two (one for black hole and one for cosmological) horizons. The main difference in Fig. 27 from βQ > 1/2 case in Fig. 25 is the "magnitude flip" between ω I(1) and ω I (2) , as in Figs. 13 and 21.   Figs. 28 and 29 show the case for βQ < 1/2 but there is no qualitative difference in the result from βQ = 1/2 case in Fig. 27.

V. EXACT SOLUTION
Generally, finding analytic expressions for QNMs is difficult and one needs to consider their numerical computations. However, there exist some special cases where exact solutions can be found. In this section, we consider the exact solution near the (charged) Nariai solution where the black hole horizon r + and the cosmological horizon r ++ merge. To this end, we first note that the metric function f (r) = e 2ν(r) in (7) and (24) can be written [28,29], near the Nariai limit r + → r ++ , as where ǫ ≡ (r ++ − r + )/r + ≪ 1 and κ + ≡ (1/2)(df /dr)| r=r + is the surface gravity at the horizon r + , which is related to the Hawking temperature T H =hκ/(2π) in (16). In this limit, the tortoise coordinate r * , for the physical region r + ≤ r ≤ r ++ , can be obtained as where we have chosen the integration constant such that r * = −∞ at r = r + and r * = +∞ at r = r ++ . Then, one can invert the relation (99) to get and Now, substituting (100) and (101) in the effective potentials U i of (87) and (88), one can find where and Here, we have used f (r * ) ∼ O(ǫ 2 ) and df /dr = −2κ + tanh(κ + r * ) + O(ǫ 3 ) ∼ O(ǫ 2 ) from κ + ∼ O(ǫ) near the Nariai limit. The potential in (102) is known as Pöshl-Teller potential [30] and its QNMs can be solved analytically [31] as where n = 0, 1, 2, · · · is the overtone mode number. In Figs. 30 -33, considering the dependence on the Hawking temperature T H or black-hole horizon radius r + , the analytic result of QNMs in (108) are compared with the numerical results based on the WKB approximations in Sec. IV. The results show quite good agreements for the imaginary parts ω I even beyond the Nariai limit (Figs. 30 and 31). The best-fit curves for the numerical results in Fig. 30 near the Nariai limit are ω I(1) /T H ≈ 3.1559, 3.0757, 3.013 and ω I(2) /T H ≈ 3.1555, 3.0887, 2.996 for U 1 and U 2 cases when β = 1, 1/2, 1/3, respectively, while the analytic result is ω I /T H = π from (108) for the lowest mode, n = 0. This result may indicate the accuracy of our WKB approach itself. This is contrary to the real parts ω R which depart from numerical results beyond the Nariai limit (Figs. 32 and 33). But this would be partly due to the non-constant or scale-dependent nature of the real part of ω/κ + in (108) for the EBI-dS case, which may be compared with the RN-dS case in [29] or the Sch-dS case in [28,32].  Before finishing this section, we note that the Nariai limit is also associated with the large-l limit,  Fig. 30. The orange lines denote the exact solutions near the Nariai limit at the bottom and the blue lines represent the best-fit curves near the limit. 2 and the linear dependence ω R ≈ σl with σ = κ + (r ++ − r + )/2r 2 + as observed in Figs. 19, 21, and 23. which may be compared with the RN case in [31] or the Sch-dS case in [33]. On the other hand, the asymptotic approach of ω I to a limiting value corresponds to ω I = κ + (n + 1/2) in (108) 12 .
We also remark that Figs. 34 and 35 show, as r + → 0, the divergent ω as with (ā,b,c,d) ≈ (1.0239, 2.7715, 0.9834, 0.4133), (1.0356, 2.9366, 1.0215, 1.0623) for U 1 , U 2 , respectively. This is consistent with the similar behavior in Fig. 16. It is interesting to note that the similar divergence can be also seen from the r + → 0 limit in the Nariai limit formula (108) with a quite close exponentā = 1 for ω I , while somewhat different exponent c = 3/2 for ω R , even though the point-like horizon limit r + → 0 would be far beyond the Nariai limit.

VI. DISCUSSION
We have studied, using the Schutz-Iyer-Will's WKB method, QNMs for the axial gravitational perturbations of electrically charged black holes in EBI gravity with a positive cosmological constant. We have found that all the axial perturbations are stable, i.e. "nonnegative" ω I , including the flat (Λ = 0) as well as dS (Λ > 0) cases, for all cases where the WKB method applies, correcting some errors in the literature [21].
It is also found that there are cases where the conventional WKB method does not apply, like the three-turning-points problems for the lower-l cases in Figs. 18 -23, and a more generalized formalism is necessary for studying their QNMs and stabilities [26]. Regarding the degenerate horizons, where the WKB formula can be marginally applicable, our results seem to be consistent with stability [19,22,34] 13 .
On the other hand, we may note in our case that there are actually three types of degenerate horizons with the vanishing Hawking temperature. The first is the Nariai-type horizons, where the (outer) black hole horizons, regardless of whether it is the RN-type or the Sch-type, coincide with the cosmological horizon. In this case, we find that the QNMs are completely "frozen", i.e., ω ≈ 0 and this would indicate the solution as the final state 12 Of course, these two behaviors may not be directly compared with those of Figs 19, 21, and 23 because of different proportional coefficients depending on whether it is near the Nariai limit or not. Actually, we have ω = −0.01763i + 0.12083 l, − 0.02278i + 0.10775 l, − 0.02089i + 0.08954 l near the Nariai limit, while ω = −0.03055i + 0.1092 l, − 0.03349i + 0.1006 l, − 0.03345i + 0.0910 l for the cases of Figs 19, 21, and 23, which correspond to ǫ = 0.99091, 0.70872, 0.57611, respectively. It is rather surprising that we have rough agreements already even beyond the Nariai limit. 13 This may be contrasted with the "horizon instability" of axisymmetric extremal horizons [35]. It would be a challenging problem to study the horizon instability of axisymmetric extremal horizons within WKB approaches or in a more generalized approach [26].