1 Introduction

Black holes (BHs), a generic prediction of Einstein’s General Relativity (GR) [1], are way more than just mathematical objects. After Hawking’s seminal work [2, 3] in which it was shown that BHs emit radiation from their horizon, these objects have attracted a lot of attention over the last decades, as they comprise an excellent laboratory to study and understand several aspects of gravitational theories. In particular, the study of BHs has received considerable attention in the context of scale dependent theories (see e.g. [4,5,6,7,8]), where certain deviations from the classical solution appear.

On the other hand, how a system responds to small perturbations has always been an important issue in Physics. The work of [9] marked the birth of BH perturbations, and it was later extended by [10,11,12,13,14]. The state-of-the art in BH perturbations is summarized in Chandrasekhar’s monograph [15]. When BHs are perturbed, the geometry of spacetime undergoes dumbed oscillations. Quasinormal modes (QNMs), with a non-vanishing imaginary part, carry unique information about the few BH parameters, since they do not depend on the initial conditions. After the LIGO direct detections of gravitational waves [16,17,18], that offer us the strongest evidence so far that BHs exist and merge, QNMs of black holes are more relevant than ever. By observing the quasinormal spectrum, that is frequencies and damping rates, we can determine the mass, angular momentum and charge of the BH, and even falsify the theoretical paradigm of the no-hair conjecture [19]. QNMs of BHs have been extensively studied in the literature. For a review on the subject see [20], and for a more recent one [21].

Relativistic scattering of waves has been traditionally studied in asymptotically flat spacetimes without a cosmological constant, such as Schwarzschild [22], Kerr [23] and Reissner–Norstrom [24] BHs, see e.g. [25,26,27,28,29,30,31]. However, due to inflation [32], the current cosmic acceleration [33, 34] and the AdS/CFT correspondence [35, 36], asymptotically non-flat spacetimes with a positive or negative cosmological constant have also been studied over the years [37,38,39,40,41,42,43,44,45,46,47]. In [48, 49], however, the authors have found black hole solutions in three and four dimensions that are neither asymptotically flat nor asymptotically (anti) de Sitter. In those works the model is described by the Einstein–Born–Infeld dilaton action. Originally the Born–Infeld non-linear electrodynamics was introduced in the 30’s in order to obtain a finite self-energy of point-like charges [50]. During the last decades this type of action reappears in the open sector of superstring theories [51, 52] as it describes the dynamics of D-branes [53, 54]. Furthermore, in the closed sector of all superstring theories at the massless level the graviton is accompanied by the dilaton that determines the string coupling constant. Since superstring theory is so far the only consistent theory of quantum gravity, it is more than natural to study the QNM of black hole solutions obtained in the framework of Einstein–Born–Infeld dilaton models.

Computing the QNM frequencies in an analytical way is possible in a few cases only [55,56,57,58,59,60,61], while in most of the cases some numerical scheme [62,63,64,65] or semi-analytical methods are employed, such as the well-known from standard quantum mechanics WKB method used extensively in the literature [66,67,68,69,70,71,72,73,74,75,76]. In the present work we obtain an exact analytical expression for the quasinormal spectrum of a four-dimensional Einstein–Born–Infeld dilaton spacetime, which is well motivated since it contains ingredients found in superstring theory. For a neutral BH, such as the Schwarzschild one, the scalar, vector and tensor perturbations can be studied separately. If, however, the BH is electrically charged then electromagnetic and gravitational perturbations are coupled and must be studied simultaneously [69]. Therefore, in this work we take the first step to study the stability against scalar perturbations by perturbing the BH with a probe scalar field \(\varPhi \), not to be confused with the dilaton \(\phi \) (see the discussion below), hoping to be able to address the problem of the coupled electromagnetic-gravitational perturbation in a future work.

Our work is organized as follows: After this introduction, we present the model and the BH solution in the next section. In Sect. 3 we discuss scalar perturbations where we present the effective potential of the Schrödinger-like equation, while in the fourth section we solve the radial equation in terms of hypergeometric functions. In Sect. 5 we obtain an exact expression for the quasinormal modes, and in Sect. 6 we compare our solution with numerical results obtained with a recently developed method [62]. Finally, we conclude our work in the last section. We use natural units such that \(c = \hbar = 1\) and metric signature \((-, +, +, +)\).

2 The model and the BH background

Our starting point is the model considered in [49] described by the action

$$\begin{aligned} S[g_{\mu \nu },A_{\mu \nu },\phi ]&= \int \mathrm {d}^4 x \sqrt{-g} \Bigl [ R-2(\nabla \phi )^2-V(\phi ) \ \nonumber \\&\quad +4 \gamma e^{-2 \kappa \phi } (1-\sqrt{1+Y}) \Bigl ], \end{aligned}$$
(1)

where

$$\begin{aligned} Y = \frac{F_{\mu \nu } F^{\mu \nu }}{2 \gamma }, \end{aligned}$$
(2)

and where R is the Ricci scalar, g is the determinant of the metric tensor \(g_{\mu \nu }\), \(F_{\mu \nu }\) is the electromagnetic field strength, \(\phi \) is the dilaton with a self-interaction potential \(V(\phi )=2 \varLambda e^{-2 \kappa \phi }\), \(\gamma \) is the Born–Infeld parameter, and \(\kappa \) is the dilaton coupling constant. Assuming static spherically symmetric solutions, the line element of the metric is found to be [49]

$$\begin{aligned} \mathrm {d} s^2=-f(r) \mathrm {d}t^2 + f(r)^{-1} \mathrm {d}r^2 + e^{2 \kappa \phi } (\mathrm {d} \theta ^2 + \sin ^2 \theta \, \mathrm {d} \varphi ^2), \end{aligned}$$
(3)

while the dilaton is given by [49]

$$\begin{aligned} \phi (r) = \frac{\kappa }{1+\kappa ^2} \ln (br-c), \end{aligned}$$
(4)

where b, c are constants of integration. In the following we set for convenience and without loss of generality \(b=1\) and \(c=0\). Since the model is string inspired, in the following we shall consider the string coupling case \(\kappa =1\). Then the line element takes the form

$$\begin{aligned} \mathrm {d}s^2&=- \left( \frac{r}{L} - r_0 \right) \mathrm {d}t^2 + \left( \frac{r}{L} - r_0 \right) ^{-1} \mathrm {d}r^2\nonumber \\&\quad + r (\mathrm {d} \theta ^2 + \sin ^2 \theta \, \mathrm {d} \varphi ^2), \end{aligned}$$
(5)

where the constant \(r_0\) is related to the mass of the black hole, \(r_0=4M\) [49], while L is given by

$$\begin{aligned} L^{-1} = 2 (1-\varLambda -2 H), \end{aligned}$$
(6)

where the constant H is given by [49]

$$\begin{aligned} H = -\gamma + \sqrt{\gamma (Q^2+\gamma )}, \end{aligned}$$
(7)

and the charge Q of the black hole is given by [49]

$$\begin{aligned} Q^2 = \frac{1+\sqrt{1+16 \gamma ^2}}{8 \gamma }. \end{aligned}$$
(8)

There is a single event horizon \(r_H=L r_0\), and therefore the metric function can be written down equivalently as \(f(r)=(r-r_H)/L\). Overall, the model is characterized by 3 free parameters, namely \(\gamma , \varLambda , M\). The horizon depends on all of them while while L does not depend on the mass of the black hole.

3 Scalar perturbations

In this section we study the propagation of a probe minimally coupled massless scalar field \(\varPhi (t,r,\theta , \varphi )\) in a given gravitational background of the form

$$\begin{aligned} \mathrm {d}s^2 = -f(r) \mathrm {d}t^2 + f(r)^{-1} \mathrm {d}r^2 + r (\mathrm {d} \theta ^2 + \sin ^2 \theta \, \mathrm {d} \varphi ^2), \end{aligned}$$
(9)

with a known metric function \(f(r)=(r-r_H)/L\). The starting point is the well-known wave equation

$$\begin{aligned} \frac{1}{\sqrt{-g}} \partial _\mu (\sqrt{-g} g^{\mu \nu } \partial _\nu ) \varPhi = 0, \end{aligned}$$
(10)

which is a partial differential equation for the scalar field. Next we seek solutions where the time and angular dependence are known as follows

$$\begin{aligned} \varPhi (t,r,\theta , \varphi ) = e^{-i \omega t} R(r) Y_l^m(\theta , \varphi ), \end{aligned}$$
(11)

with \(Y_l^m\) being the usual spherical harmonics. Using the above ansatz it is straightforward to obtain the radial equation, which is an ordinary differential equation

$$\begin{aligned} R'' + \left( \frac{h'}{h}+\frac{1}{r}\right) R' + \left( \frac{\omega ^2}{h^2}-\frac{l (l+1)}{r h}\right) R = 0, \end{aligned}$$
(12)

where the prime denotes differentiation with respect to radial distance r. Next, we recast the equation for the radial part into a Schrödinger-like equation of the form

$$\begin{aligned} \frac{\mathrm {d}^2 \psi }{\mathrm {d} x^2} + (\omega ^2 - V(x)) \psi = 0, \end{aligned}$$
(13)

by defining new variables, a dependent \(R \rightarrow \psi \) as well as an independent one \(r \rightarrow x\) as follows

$$\begin{aligned} R= & {} \frac{\psi }{\sqrt{r}}, \end{aligned}$$
(14)
$$\begin{aligned} x= & {} \int \frac{\mathrm {d} r}{f(r)}=L \ln \Bigl (\frac{r-r_H}{d}\Bigl ), \end{aligned}$$
(15)

with x being the so-called tortoise coordinate, and d is a constant of integration which will be taken as unity. Therefore, we obtain the expression for the effective potential

$$\begin{aligned} V(r) = h(r) \, \left( \frac{l (l+1)}{r}+\frac{h'(r)}{2 r}-\frac{h(r)}{4 r^2} \right) , \end{aligned}$$
(16)

which can be simplified to be

$$\begin{aligned} V(r) = V_0 - \frac{r_H l (l+1)}{L r} - \frac{r_H^2}{4 L^2 r^2}, \end{aligned}$$
(17)

where the constant term is given by \(V_0=(Ll(l+1)+1/4)/L^2\). The effective potential barrier versus the radial coordinate can be seen in Figs. 1 and 2 where we plot it for different \(\gamma \) and \(\varLambda \), respectively. Since it does not exhibit a maximum the WKB approximation is not applicable, and therefore we shall turn our attention to a different method for the numerical verification of our main analytical expression, see eq. (43) below.

Fig. 1
figure 1

Effective potential versus r for \(l = 0, M=2, \varLambda =0.1\) and \(\gamma =0.1\) (solid red line), \(\gamma =0.5\) (dashed blue line) and \(\gamma =2\) (dotted-dashed magenta line)

Full size image
Fig. 2
figure 2

Effective potential versus r for \(l = 0, M=2, \gamma =0.5\) and \(\varLambda =0.001\) (solid red line), \(\varLambda =0.01\) (dashed blue line) and \(\varLambda =0.1\) (dotted-dashed magenta line)

Full size image

Finally, the Schrödinger-like equation must be supplemented by appropriate boundary conditions at horizon and at infinity, which are the following [77]

$$\begin{aligned} \psi (x) \rightarrow \left\{ \begin{array}{lcl} A e^{i \omega x} &{} \text{ if } &{} x \rightarrow - \infty \\ &{} &{} \\ C_+ e^{i k_{\infty } x} + C_- e^{-i k_{\infty } x} &{} \text{ if } &{} x \rightarrow + \infty \end{array} \right. \end{aligned}$$
(18)

where \(A, C_+, C_-\) are arbitrary constants, and \(k_{\infty }\) depends on the value of the effective potential at infinity. Up to now, following the procedure just described one can compute the so-called greybody factors, which show the modification of the original spectrum of Hawking radiation due to the effective potential barrier, and where the frequency is real and takes continuous values. For an incomplete list see e.g. [25,26,27,28,29,30,31, 37,38,39,40,41,42,43,44,45,46,47] and references therein. Now, the QNMs are determined requiring that the first coefficient of the second condition vanishes, i.e. \(C_+ = 0\) [78]. The purely ingoing wave physically means that nothing can escape from the horizon, while the purely outgoing wave corresponds to the requirement that no radiation is incoming from infinity. We thus obtain an infinite set of discrete complex numbers \(\omega _n=\omega _R + \omega _I i\), which are precisely the QNM frequencies of the black hole. Given the time dependence of the probe scalar field \(\varPhi \sim e^{-i \omega t}\), it is clear that unstable modes correspond to \(\omega _I > 0\), while stable modes correspond to \(\omega _I < 0\). In addition, the real part of the mode \(\omega _R\) determines the period of the oscillation, \(T=2 \pi /\omega _R\), while the imaginary part \(|\omega _I|\) describes the decay of the fluctuation at a time scale \(t_D=1/|\omega _I|\).

Since at the horizon the effective potential vanishes, the general solution for the function \(\psi \) close to the horizon (where \(\omega ^2 \gg V(x)\)) is given by

$$\begin{aligned} \psi (x) = A_{+} e^{i \omega x} + A_{-} e^{-i \omega x}, \end{aligned}$$
(19)

while requiring purely ingoing solution we set \(A_{-}=0\) [44, 46], and thus the solution becomes

$$\begin{aligned} \psi (x) = A e^{i \omega x}. \end{aligned}$$
(20)

On the other hand, it is easy to check that at large r (or at large x, since when \(r \gg r_H\), \(r \simeq e^{x/L}\)) the potential tends to the constant \(V_0\), and therefore defining \(\varOmega \equiv \sqrt{\omega ^2-V_0}\) the solution for \(\psi \) is given by

$$\begin{aligned} \psi (x)&= D_{+} e^{i \varOmega x} + D_{-} e^{-i \varOmega x}. \end{aligned}$$
(21)

Therefore, the far-field solution expressed in the tortoise coordinate x takes the form of ingoing and outgoing plane waves provided that \(\omega ^2 > V_0\), while the QNMs are determined by requiring that \(D_+=0\).

4 Solution of the full radial equation in terms of hypergeometric functions

Next, we find an exact solution of the radial equation (12) in terms of hypergeometric functions by introducing \(z=1-r_H/r\). The new equation for z reads

$$\begin{aligned} z (1-z) R_{zz} + (1-z) R_z + \left( \frac{A}{z} + \frac{B}{-1+z} \right) R = 0, \end{aligned}$$
(22)

where \(A=(\omega L)^2, B=-(\omega L)^2+L l (l+1)\). To get rid of the poles we set

$$\begin{aligned} R = z^\alpha (1-z)^\beta F, \end{aligned}$$
(23)

where now F satisfies the following differential equation

$$\begin{aligned}&z (1-z) F_{zz} + [1+2 \alpha - (1+2 \alpha +2 \beta ) z] F_z\nonumber \\&\quad + \left( \frac{\bar{A}}{z} + \frac{\bar{B}}{-1+z} - C \right) F = 0, \end{aligned}$$
(24)

and the new constants are given by

$$\begin{aligned} \bar{A}= & {} A + \alpha ^2, \end{aligned}$$
(25)
$$\begin{aligned} \bar{B}= & {} B + \beta - \beta ^2, \end{aligned}$$
(26)
$$\begin{aligned} C= & {} (\alpha +\beta )^2. \end{aligned}$$
(27)

Demanding that \(\bar{A} = 0 = \bar{B}\) we obtain the Gauss’ hypergeometric equation

$$\begin{aligned} z (1-z) F_{zz} + [c-(1+a+b) z] F_z - ab F = 0, \end{aligned}$$
(28)

and we determine the parameters \(\alpha , \beta \) as follows

$$\begin{aligned} \alpha= & {} i \omega L, \end{aligned}$$
(29)
$$\begin{aligned} \beta= & {} \frac{1}{2} + i \sqrt{(\omega L)^2 - L l (l+1) -\frac{1}{4}} . \end{aligned}$$
(30)

Finally, the three parameters of Gauss’ equation are given by

$$\begin{aligned} c= & {} 1+2 \alpha , \end{aligned}$$
(31)
$$\begin{aligned} a= & {} \alpha + \beta , \end{aligned}$$
(32)
$$\begin{aligned} b= & {} \alpha + \beta . \end{aligned}$$
(33)

Note that the parameters a, b, c satisfy the condition \(c-a-b=1-2 \beta \). Therefore, the solution for the radial part is given by

$$\begin{aligned} R(z) = D z^\alpha (1-z)^\beta F(a,b;c;z), \end{aligned}$$
(34)

where D is an arbitrary coefficient, and the hypergeometric function can be expanded in a Taylor series as follows

$$\begin{aligned} F(a,b;c;z) = 1 + \frac{a b}{c} \,z + \cdots \end{aligned}$$
(35)

Note that the above solution for the choice of \(\alpha =i \omega L\) reproduces the purely ingoing solution at the horizon (20), as it can be seen from the fact that close to the horizon (\(z \rightarrow 0\)), the radial part becomes \(R(z) \simeq D z^\alpha \), and the parameter z can be written approximately \(z \simeq (r-r_H)/r_H = e^{x/L}/r_H\).

5 Exact quasinormal spectrum

To see how the radial part behaves in the far-field zone \(r \gg r_H\) (where \(z \rightarrow 1\)) we use the transformation [79]

$$\begin{aligned} \begin{aligned} F(a,b;c;z) = \ \frac{\varGamma (c) \varGamma (c-a-b)}{\varGamma (c-a) \varGamma (c-b)} \ \\&\times F(a,b;a+b-c+1;1-z) \ +\\ (1-z)^{c-a-b}&\frac{\varGamma (c) \varGamma (a+b-c)}{\varGamma (a) \varGamma (b)} \ \\&\times F(c-a,c-b;c-a-b+1;1-z), \end{aligned} \end{aligned}$$
(36)

and therefore the radial part as \(z \rightarrow 1\) reads

$$\begin{aligned} R(z \rightarrow 1)&=D (1-z)^\beta \frac{\varGamma (1+2 \alpha ) \varGamma (1-2 \beta )}{\varGamma (1+\alpha -\beta ) \varGamma (1+\alpha -\beta )} \nonumber \\&\quad +D (1-z)^{1-\beta } \frac{\varGamma (1+2 \alpha ) \varGamma (-1+2 \beta )}{\varGamma (\alpha +\beta ) \varGamma (\alpha +\beta )}. \end{aligned}$$
(37)

Since \(z=1-(r_H/r)\), the radial part R(r) for \(r \gg r_H\) can be written down as follows

$$\begin{aligned} R(r) \simeq D \frac{\varGamma (1+2 \alpha ) \varGamma (1-2 \beta )}{\varGamma (1+\alpha -\beta ) \varGamma (1+\alpha -\beta )} \left( \frac{r_H}{r}\right) ^{\beta } \nonumber \\ +D \frac{\varGamma (1+2 \alpha ) \varGamma (-1+2 \beta )}{\varGamma (\alpha +\beta ) \varGamma (\alpha +\beta )} \left( \frac{r_H}{r}\right) ^{1-\beta }. \end{aligned}$$
(38)

Upon defining new constants \(D_1, D_2\) as follows

$$\begin{aligned} D_1= & {} D \, \frac{\varGamma (1+2 \alpha ) \varGamma (1-2 \beta )}{\varGamma (1+\alpha -\beta ) \varGamma (1+\alpha -\beta )}, \end{aligned}$$
(39)
$$\begin{aligned} D_2= & {} D \, \frac{\varGamma (1+2 \alpha ) \varGamma (-1+2 \beta )}{\varGamma (\alpha +\beta ) \varGamma (\alpha +\beta )}, \end{aligned}$$
(40)

the function \(\psi \) that satisfies the Schrödinger-like equation takes the form

$$\begin{aligned} \psi \simeq D_1 \ r_H^\beta e^{-i \text {Im}(\beta ) \frac{x}{L}} + D_2 \ r_H^{1-\beta } e^{i \text {Im}(\beta ) \frac{x}{L}}. \end{aligned}$$
(41)

In the final step of the calculation we apply the boundary condition \(D_2=0\), as we mentioned when we discussed scalar perturbation in Sect. 3. We require that the Gamma function in the denominator has a pole, and therefore the QNMs are determined imposing the condition

$$\begin{aligned} \alpha + \beta = -n, \end{aligned}$$
(42)

with \(n=0,1,2,...\) being the overtone number. Using the previous expressions for \(\alpha \) and \(\beta \) we obtain the formula

figure a

which is our main result in the present work. We can immediately observe the following three features of the spectrum, namely a) the QNMs are purely imaginary, b) they do not depend on \(r_H\), so they depend on \(\gamma \) and \(\varLambda \) only, but not on the mass of the black hole, and c) all modes for \(n=0\) become \(\omega _{n=0}=-l (l+1) i\), and therefore they do not depend on any of the BH properties. In particular, the fundamental mode \(l=0,n=0\) is precisely zero, while for \(l > 0\) all modes corresponding to \(n=0\) are stable. In Table 1 we compare our exact values with the ones computed numerically, while in the Figs. 3 and 4 we show how the imaginary part of the frequencies change with \(\gamma \) and with \(\varLambda \) respectively for \(l=0, n=1\). Our figures show that the mode changes sign depending on the value of the cosmological constant \(\varLambda \) as well as the Born–Infeld parameter \(\gamma \) (or equivalently the electric charge Q). In particular, for low charge (large \(\gamma \)) the imaginary part is positive, whereas as the charge grows at a certain point the imaginary part becomes negative. Interestingly enough, a behaviour similar to that found in [80, 81] is observed, although the scalar field that perturbs the BH in the present work is not electrically charged. Contrary to these works, however, where it was found that all modes with \(l >0\) were stable, our results show that for any value of the angular momentum there is a certain value of the overtone number after which the modes become unstable.

Table 1 Scalar QNMs of Einstein–Born–Infeld dilaton black hole for various values of \(\gamma \) and \(\varLambda \). \(l,\,n\) are the angular momentum and overtone number, respectively, and \(L^{-1}=2(1-\varLambda -2H)\), see text. The values without the parenthesis are the exact QNMs, while the ones in the parenthesis are the numerical values
Full size table
Fig. 3
figure 3

Imaginary part of the QN modes versus \(\gamma \) for \(l = 0, n=1\) and \(\varLambda =0.2\) (solid red line), \(\varLambda =0.25\) (dashed blue line) and \(\varLambda =0.3\) (dotted-dashed magenta line)

Full size image
Fig. 4
figure 4

Imaginary part of the QN modes versus \(\varLambda \) for \(l = 0, n=1\) and \(\gamma =0.1\) (solid red line), \(\gamma =0.15\) (dashed blue line) and \(\gamma =0.2\) (dotted-dashed magenta line)

Full size image

6 Numerical results

We review a non grid-based interpolation scheme, proposed by Lin et al. [62]. This method makes use of data points in a small region of a query point to estimate its derivatives by employing Taylor expansion. The data points can be scattered, therefore they do not sit on a grid. A key step of the method is to discretize the unknown eigenfunction in order to transform a differential equation and its boundary conditions into a homogeneous matrix equation. Based on the information about N scattered data points, Taylor series are carried out for the unknown eigenfunction up to N-th order for each discretized point. The resulting homogeneous system of linear algebraic equations is solved for the eigenvalue. A huge advantage of this method is that the discretization of the wave function and its derivatives are made to be independent of any specific metric through coordinate transformation.

This method has been tested thoroughly for its accuracy and efficiency to various differential equation and eigenvalue problems in [62, 63], and [64]. The QNM results have been compared with WKB approximation [66, 67] (up to the sixth order), Horowit z-Hubeny method [82], and continued fraction method [83] achieving very good precision.

We have applied the present method to compute the scalar QNMs of a four-dimensional Einstein–Born–Infeld dilaton black hole, and our numerical results are summarized in Table 1.

We observe that the numerical results agree perfectly with our main result shown in eq. (43). We immediately see that for \(n=0\) the modes depend only on the angular momentum l, irrespectively of the choice of \(\gamma \) and \(\varLambda \), which is in agreement with \(\omega _{n=0}=-l (l+1) i\) obtained in the previous section. For \(l=0\), the aforementioned fundamental mode is exactly zero, while for \(l\ge 1\), as l increases, the modes decay faster while more stable overtones begin to appear. This is to be expected, since larger angular momentum increases the height of the potential stabilizing the system. Finally, we observe that as L decreases the stable overtones decay slower while the unstable ones grow faster.

7 Conclusions

To summarize, in this article we have studied the stability under scalar perturbations of (1+3)-dimensional Einstein–Born–Infeld dilaton spacetimes, and we have provided an exact analytical expression for the frequencies, which are found to be purely imaginary. We have confirmed our results computing the frequencies numerically using a recently developed non grid-based numerical scheme. In addition, an instability similar to that seen in charged scalar perturbations of the Reissner–Nordström black hole is observed, although in our work the scalar field that perturbs the BH is not electrically charged.