Bernstein spectral method for quasinormal modes of a generic black hole spacetime and application to instability of dilaton-de Sitter solution

Proper oscillation frequencies of black holes, quasinormal modes Kokkotas:1999bd ; Berti:2009kk ; Konoplya:2011qq usually cannot be found analytically and a number of numerical or semianalytical methods are used. Simple and quick and, in a sense, automatic methods, such as the WKB one Mashhoon:1982im ; Schutz:1985zz ; Iyer:1986np ; Konoplya:2003ii ; Konoplya:2003ii ; Matyjasek:2017psv ; Konoplya:2019hlu unfortunately do not allow one to find quasinormal modes in all required cases, for example, when the real oscillation frequency is much smaller than the damping rate or zero Konoplya:2019hlu . The latter case includes important situations: instability which is represented by exponentially growing modes without oscillations when one is constrained by spherically symmetric solutions Konoplya:2008yy . Unstable modes can, in principle be detected via the Frobenius method Leaver:1985ax , yet it is based on trial and error procedure when searching roots of algebraic equations and there is a risk that some of the roots could be simply missed.

The method which is free of such deficiency is the time-domain integration of the wave equation at a fixed spacial coordinate Gundlach:1993tp . This method includes contribution of all overtones and the instability is simply shown as an exponential growth of the absolute values of the wave function at late times. However, higher overtones cannot be easily extracted from the time-domain profile within this approach and one is usually constrained by detecting two-three first overtones with reasonable accuracy. Moreover, even detecting instability may be challenging, because the instability may start at very late times after a long period of damped oscillations, as it occurs, for example, for black holes with higher curvature corrections Konoplya:2008ix . The latter would require time-consuming numerical integration of the wave equations until very late times.

The method which is efficient for finding purely imaginary, that is, nonoscillatory, quasinormal modes is a spectral method based on the Bernstein polynomial nonorthogonal basis functions Fortuna:2020obg . As was shown in Fortuna:2020obg while reproducing only several first overtones during a short computing time, it perfectly detects the algebraically special mode of the Schwarzschild black hole with high accuracy, which occurs at $n=8$, where $n$ is the overtone number. Another spectral method based on the Chebyshev polynomials was used to obtain a family of the purely imaginary modes of the Schwarzschild-de Sitter black holes Jansen:2017oag . The spectral method was also successfully employed to detect the instabilities of the rotating black holes in AdS Monteiro:2009ke ; Dias:2010 ; Dias:2013 and higher-dimensional rotating black holes Dias:2009 and black strings Dias:2022mde (see also Dias:2015nua for review). This makes us expect that the Bernstein spectral method will be also useful for finding the instability region.

Here, using the initial code of Fortuna:2020obg as a basis, we develop it to be even more user-friendly and solve two main problems related to it:

• •

  The Bernstein spectral method produces a great number of roots of the matrix equation, not all of which are true quasinormal frequencies, so that only varying the size of the matrix, one can see which roots survive and distinguish them as quasinormal modes. We incorporated this procedure of the comparison of matrices into the code.

• •

  The code now can be automatically applied (without further adjustment) to asymptotically flat, de Sitter or anti-de Sitter spacetimes.

If the wave equation depends on a number of parameters, the problem of finding the instability region via the time-domain integration is extremely time consuming, because the integration must be performed plenty of times to cover all the possible values of the parameters. Thus, a problem of this kind could be an excellent playground for the Bernstein spectral method.

Such an example could provide the instability of a charged scalar field in the background of the charged asymptotically de Sitter black hole Zhu:2014sya ; Konoplya:2014lha . The region of instability in this case depends on the charges of the field and black hole, the black hole mass and the cosmological constant. While the instability of the charged scalar field in the Reissner-Nordström-de Sitter background was studied in Zhu:2014sya ; Konoplya:2014lha ; Dias:2018ufh , no such analysis was done for the dilatonic analog of this configuration. Indeed, quasinormal modes of a neutral scalar field around the dilatonic black hole at nonzero cosmological constant was considered for the first time in Fernando:2016ftj with the help of the WKB approximation, while the charged scalar filed of the asymptotically flat case was considered in Konoplya:2002ky . Thus, to the best of our knowledge, no quasinormal modes of the charged scalar field in the background of the asymptotically de Sitter dilatonic black hole were considered so far.

Notice, that the dilaton-de Sitter black hole is interesting on its own, because of the holographic duality relating quantum gravity on the de Sitter space and the conformal field theory Witten:2001kn . Then, the quasinormal modes of black holes are important not only from the point of view of observation of gravitational waves, but also as poles of the corresponding Green functions in the dual field theory, allowing to describe the relaxation of the corresponding quantum fields at finite temperature Son:2007vk . In the low energy limit of string theory, the Einstein action is modified, among other terms, by a scalar, dilaton, field.

Thus, our paper is two-fold: on the one hand it represents the improved Bernstein polynomial procedure for finding quasinormal modes of a broad class of black holes and, on the other hand, the code is used here for investigating the instability of the asymptotically de Sitter dilatonic black holes.

The paper is organized as follows. In section II we provide the basic information about the dilatonic black hole in the presence of cosmological constant. Section III is devoted to the discussion of the nature of instability – superradiation. Section IV briefly reviews the time-domain integration which we also used as a complementary method for checking the obtained results. The Bernstein spectral method is discussed in section V. In section VI we discuss the obtained numerical data on the instability region. Finally, in the conclusion we summarize the obtained results and mention some open problems.

The well-known four-dimensional and asymptotically flat charged black-hole solutions with the dilaton field were obtained by Gibbons and Maeda Gibbons:1987ps and, independently, by Garfinkle et al. Garfinkle:1990qj . However, as was shown in Poletti:1994ff , no such extension is possible in the presence of cosmological constant with exponential coupling. Instead, the dilation potential can be represented is the sum of the cosmological constant (with some dimensionless coefficient) and two Liouville type terms Gao:2004tu . Here we will study this form of the dilaton potential.

The action for the dilaton gravity can be written as,

where $R$ is the scalar curvature, $F_{\mu\nu}$ is the Maxwell’s field strength and $\varphi$ is the dilation field. The potential for the dilation field is given by $V(\varphi)$, which has the following form

The metric of the dilaton-de Sitter black hole is Gao:2004tu

where

Here, $\Lambda$ is the cosmological constant, $M$ is the mass and, $D$ is the dilation charge. Notice that when $\Lambda=0$, the black hole solution in Eq. (3) becomes the Gibbons-Maeda-Garfinkle-Horowitz-Strominger (GMGHS) black hole Gibbons:1987ps ; Garfinkle:1990qj . When $D=0$, the space-time reduces the Schwarzschild-de Sitter solution. Some basic properties of this black hole solution were studied in Fernando:2016ftj ; Benakli:2021fvv ; Dehyadegari:2020tau .

The dilation field $\varphi$, dilation charge $D$, and electric field $F_{01}$, for the above solution are given by the following relations,

where $\varphi_{0}$ is the dilation field at $r\to\infty$ and $Q$ is the electric charge of the black hole.

The two largest roots of the equation $f(r)=0$ are the event $r_{+}$ and cosmological $r_{c}$ horizons. For large mass $M$ there are no horizons and the spacetime is a naked singularity. Since the metric has a singularity at $r=2D$, the space-time becomes a black hole only if the radius of the event horizon $r_{+}>2D$.

We shall further compare the dilaton-de Sitter black hole with the Reissner-Nordström-de Sitter black hole, $\varphi=0$,

We designate the Cauchy, event, and cosmological horizons as $r_{-}$, $r_{+}$, and $r_{c}$ respectively ($r_{-}<r_{+}<r_{c}$). Then, the metric function $f(r)$ can be written as,

where $r_{n}=r_{-}+r_{+}+r_{c}$ for the Reissner-Nordström-de Sitter black hole (6) and $r_{n}=0$ for the dilaton-de Sitter black hole (4).

A charged, massive scalar field $\phi$ in curved space-time obeys the Klein-Gordon equation

where $e$ and $\mu$ are, respectively, the charge and mass of the field and

is the electromagnetic 4-potential of the black hole.

After the standard separation of angular variables with the help of spherical harmonics and introduction a new wave function $\Psi$, the above equation of motion can be reduced to the following form:

where

is the tortoise coordinate, the dilaton field is

and the effective potential has the form

The above effective potential $V(r)$ has the form of the positive definite potential barrier with a single maximum for $\ell\geq 1$, but it has a negative gap for $\ell=0$ case (see Fig. 1).

Here we will we consider the scattering problem for a charged scalar field in the background of the dilatonic de Sitter black hole. Owing to the $\partial_{t}$ Killing vector for stationary space-times, the following ansatz $\Psi(r,t)=\exp[-i\omega t]\Psi(r)$ transforms (9) to the Schrödinger-like form

Due to the following symmetry of (12):

it is sufficient to consider only the range $eQ>0$

We will use the usual scattering boundary conditions, implying that the wave coming from the cosmological horizon will be partially reflected by the potential barrier and partially come back to the cosmological horizon, while at the event horizon of a classical black hole purely incoming wave is always required:

Since the Wronskian of the complex conjugated solutions is constant, one can derive the following relation

Amplification of the incident wave, that is, a superradiance, occurs when

When the cosmological constant $\Lambda$ vanishes, the superradiant condition Eq. (17) is reduced to the one for asymptotically flat space-times formulated by Bekenstein Bekenstein:1973mi ; Hod:2012-2013 . Notice, that $\omega$ is the real oscillation frequency in (17).

In order to understand how the condition of superradiance is related to the instability we will consider, generally, complex frequencies and impose the quasinormal boundary conditions,

Following Konoplya:2014lha we will show that the real part of $\omega$ satisfying the superradiance inequality (17) is the necessary, but not sufficient, condition for the instability. For this purpose we multiply (12) by the complex conjugated $\Psi^{\star}$ and integrate the first term by parts,

The right-hand side is real, since the effective potential $V(r)$ is real. Because $\Phi(r)$ is a monotonically decreasing function, taking imaginary part of both sides, we find that if either $Re(\omega)\geq\Phi(r_{+})\geq\Phi(r)$ or $Re(\omega)\leq\Phi(r_{c})\leq\Phi(r)$ is satisfied, then we have $Im(\omega)<0$. This proves that the instability can take place only provided

For $eQe^{2\varphi_{0}}$ larger than the threshold of instability the real part of the dominant mode also satisfies Eq. (22), being damped. This is qualitatively different from the higher-dimensional asymptotically anti-de Sitter black holes, for which the necessary condition is also the sufficient one, when the black hole is charged or rotating Kodama:2009rq ; Uchikata:2011zz ; Li:2012rx ; Wang:2014eha .

One should notice that the asymptotic value of the dilaton field $\phi_{0}$ only changes the units of the dilaton charge $D$ and electric charge of the test field $e$. Therefore, without loss of generality we shall further consider $\phi_{0}=0$.

When the event horizon radius is close to the de Sitter radius, that is the black hole occupies almost the whole de Sitter space, the spectral problem can be treated analytically Cardoso:2003sw ; Molina:2003ff ; Churilova:2021nnc . In this extreme limit we expand the value of the surface gravity $\kappa$ in terms of a small difference between the radii of the event and cosmological horizons,

Then we obtain

and the tortoise coordinate takes the following simple form:

Therefore, one can find a closed form for $r$ in terms of $r_{\star}$,

Substituting Eq. (25) into the above expression for $f(r)$ we have

Using Eq. (25) in (11) we see that the potential approaches the Pöschl-Teller potential in the near-extreme limit,

where the nonnegative constant $A$ is defined as follows,

and

We notice that $V(r_{\star})={\cal O}\left(\kappa\right)^{2}$. Therefore, expanding the frequency in the series with respect to $\kappa$, we find that

where $\omega_{\kappa}$ is the eigenvalue of the following equation:

which is written in terms of the new coordinate

Using the analytic formula for the Pöschl-Teller potential eigenvalues, we find, for small $eQ$,

implying that the scalar field in the near-extreme limit is stable for sufficiently small $eQ$, unless $A=0$, that is,

When $\ell=0$ and $\mu=0$, then $A=0$ and the imaginary part of the dominant quasinormal mode ($n=0$),

which cannot be computed within the simple approach presented in this section. However, Eq. (32) allows us to conclude that whatever small mass $\mu>0$ stabilizes the scalar field, at least for sufficiently small values of its charge $e$.

For the Reissner-Nordström-de Sitter black hole the instability was observed only when $\ell=0$ and $\mu=0$ Zhu:2014sya ; Konoplya:2014lha . Therefore, when studying the parametric region of the scalar-field instability in the background of the dilaton-de Sitter black hole, we will represent here mostly the spherically symmetric ($\ell=0$) perturbations of the massless scalar field ($\mu=0$).

Here in order to integrate the wave equation (9) and analyze the spectrum of the perturbation, we shall use, first of all, the time-domain integration, which includes contribution from all modes. The discretization scheme which we shall use was proposed in Abdalla:2010nq . Defining $\Psi(r_{\star},t)=\Psi(j\Delta r_{\star},i\Delta t)=\Psi_{j,i}$, $V(r(r_{\star}))=V(j\Delta r_{\star})=V_{j}$ and $\Phi(r(r_{\star}))=\Phi(j\Delta r_{\star})=\Phi_{j}$, we can write down (9) as

The initial Gaussian wave-package has the form

Then, the evolution of $\Psi$ can be described by the following expression

Following Abdalla:2010nq , we choose the parameters $a=0$ and $b=\sqrt{10}$ in the Gaussian wave package and use

making sure that $\Delta t$ is small enough to achieve the required precision of the profile.

In order to calculate $Re(\omega)$ we used the Prony method of fitting the time-domain profile data by superposition of damping exponents Berti:2007dg

We consider a late time period, which starts at $t_{0}$ and ends at $t=N\Delta t+t_{0}$, where $N$ is an integer and $N\geq 2p-1$. Then the formula (34) is valid for each value from the profile data:

The Prony method allows us to find $z_{i}$ in terms of the known $x_{n}$ and, since $\Delta t$ is also known, to calculate the quasinormal frequencies $\omega_{i}$.

First, we will discuss the Bernstein spectral method for asymptotically de Sitter spacetimes, when the purely outgoing wave is imposed at the de Sitter horizon. Following Fortuna:2020obg , we introduce the compact coordinate, which is defined as follows:

and define the function $\psi(u)$, which is regular for $0\leq u\leq 1$ when $\omega$ is a quasinormal mode,

where $\Omega_{c}(\omega)$ and $\Omega_{+}(\omega)$ are obtained from the characteristic equations at the singular points, $u=0$ and $u=1$, of the wavelike equation (12).

In order to satisfy the quasinormal boundary conditions (20), we fix the values of $\Omega_{c}(\omega)$ and $\Omega_{+}(\omega)$, such that

We represent $\psi(u)$ as a sum

where

are the Bernstein polynomials.

Substituting (37) into (12) and using a Chebyshev collocation grid of $N+1$ points,

we obtain a set of linear equations with respect to $C_{k}$, which has nontrivial solutions iff the corresponding coefficient matrix is singular. Since the elements of the coefficient matrix are polynomials (of degree 2) of $\omega$, the problem is reduced to the eigenvalue problem of a matrix pencil (of order 2) with respect to $\omega$, which can be solved numerically. Once the eigenvalue problem is solved, one can calculate the corresponding coefficients $C_{k}$ and explicitly determine the polynomial (39), which approximates the solution to the wave equation (12).

In order to exclude the spurious eigenvalues, which appear due to finiteness of the polynomial basis in (39), we compare both the eigenfrequencies and corresponding approximating polynomials for different values of $N$. First, from each set of the solutions we take the eigenvalues that differ less than the required accuracy. Then, following Konoplya:2022xid , for each pair of the corresponding eigenfunction, $\psi^{(1)}$ and $\psi^{(2)}$, we calculate

where $\alpha$ is the angle between the vectors $\psi^{(1)}$ and $\psi^{(2)}$ in the $L^{2}$-space. If all values of $\alpha$ are sufficiently small, we conclude that the obtained eigenvalues $\omega$ approximate the quasinormal frequencies, and better approximations correspond to larger $N$. The error estimation can be done by calculating the difference between the approximate eigenvalues of $\omega$, corresponding to different values of $N$.

Notice that, unlike Fortuna:2020obg , we compare the polynomial approximations for the eigenfunctions without a normalization, so that the eigenfunctions, obtained for different $N$, have different complex constant prefactor.

It is a valuable feature of the Bernstein method that it can be similarly applied to asymptotically flat ($\Lambda=0$) and AdS black holes ($\Lambda<0$), which correspond to a qualitatively different boundary condition at infinity. In this case the compact coordinate is introduced as follows

Notice that although for $r_{c}\to\infty$ the quasinormal spectrum of the asymptotically de Sitter black hole approaches the one of the flat black hole and Eq. (36) leads to the correct expression for the compact coordinate (40). The prefactor in (37) and the equation for the regular function $\psi(u)$ cannot be obtained by taking this limit. The reason is that the singular point $u=0$ becomes irregular for the asymptotically flat spacetime. The regular function $\psi(u)$ is defined as follows:

where $\Omega_{c}(\omega)$, $\alpha(\omega)$, and $\Omega_{+}(\omega)$ are again determined by solving the characteristic equations, and the quasinormal boundary conditions imply that $\Omega_{c}(\omega)$ and $\Omega_{+}(\omega)$ satisfy Eq. (38).¹¹1Similarly, the singular point $u=1$ becomes irregular for the extreme charge of the Reissner-Nordström-(A)(dS) black hole.

For the AdS black hole, the regular function $\psi(u)$ is introduced as

so that the function $\Psi(u)$ obeys the Dirichlet boundary condition at spatial infinity ($u=0$).

After obtaining of the differential equation for the regular function $\psi(u)$, one applies the expansion (39) for dS, flat, and AdS case, and finally applies the spectral method Jansen:2017oag ; Fortuna:2020obg .

We publicly share the Mathematica® package allowing one to obtain the equation for the regular function $\psi(u)$, defined by Eqs. (37) (de Sitter), (41) (flat), or (42) (anti-de Sitter), solve the eigenvalue problem for the finite Bernstein polynomial basis, and compare different approximations package . Indeed, the region of instability shown in Fig. 2 is visually indistinguishable when plotted via the data from the time-domain integration or Bernstein spectral method.

We calculated quasinormal modes with the help of the two above methods: time-domain integration and the Bernstein spectral method. Both methods very well agree on the region of instability, though the time-domain integration is much more (at least two orders) time consuming. Indeed, the curves in Fig. 2 are visually indistinguishable when constructed on the data from either of the two methods. The larger is $r_{+}$ the larger is the stability region, that is, at the smaller $eQ$ the configuration is stabilized. The regime of tiny $eQ$ for the weakly charged Reissner-Nordström-de Sitter black hole is unstable unless $r_{+}$ is not larger than some critical value (see Fig. 6 in Konoplya:2014lha ). The dilaton charge $D$ changes the situation: at a fixed $r_{+}$ and $r_{c}$ the system is stabilized at a larger $eQ$. In a similar fashion with the Reissner-Nordström-de Sitter solution Konoplya:2014lha ) larger values of $r_{+}$ make the region of stability larger.

The data given in the tables shows that the fundamental quasinormal modes computed by the Bernstein spectral method are in excellent agrement with those extracted from the time-domain profiles. A tiny difference between them, must be interpreted in favor of the Bernstein spectral method, because in order to find quasinormal frequency from the time-domain profile, one needs not only to integrate the weave equation with enormous accuracy, but also accurately guess the beginning and end of the ringdown period, in order to avoid the mixture of the initial outburst and asymptotic tails.

From the Tables 1 and 2 one can see that the region of instability indeed can be determined with the help of the Bernstein spectral method, which allows us to calculate the dominant modes (either with positive or negative imaginary part) with good accuracy. In Table 3 we show that the method is also efficient for the accurate calculation of the overtones. This can be seen by comparing the obtained frequencies with those extracted from the time-domain profiles for larger multipole number, $\ell\geq 1$, when the ringing occurs for sufficiently long time, allowing us to determine the dominant modes with higher accuracy.

With the help of the spectral method one can easily check that the instability occurs for the spherically symmetric perturbations ($\ell=0$) only. Using the Mathematica® package package , it is possible to find the accurate value of the scalar-field instability threshold (critical scalar-field charge) for any given black-hole parameters, $r_{+}$, $r_{c}$, and $D$.

Quasinormal modes of gravitational and other spin field perturbations of the Schwarzschild-de Sitter black holes are known to consist from the two branches: the Schwarzschild-like modes deformed by the nonzero value of the cosmological constant and the purely imaginary (i.e., nonoscillatory) modes of the empty de Sitter spacetime Lopez-Ortega:2006aal corrected by the nonzero mass of the black hole. When the black hole is small in comparison with the cosmological scale, the de Sitter modes become the least damping ones, thereby, dominating in a signal at late times and showing itself as exponential tails Konoplya:2022xid . The de Sitter branch cannot be reproduced with the most frequently used WKB method Konoplya:2022gjp and requires usage of accurate methods based on the converging procedures, such as the Frobenius expansion or Bernstein spectral method.

In Konoplya:2022kld it was shown that the above result is general for small black holes in an arbitrary metric theory of gravity immersed in the de Sitter spacetime: the dominant quasinormal modes of the uncharged massless fields are purely imaginary, approaching the de Sitter modes according to the following universal law

where $s=0$ for the scalar field and $\tilde{n}=0,1,2\ldots$ enumerates the purely imaginary modes.

In Table 4 we compare the above analytic formula (43) and the numerical values obtained with the help of the Bernstein spectral method. We can see that the Bernstein method is very efficient not only for finding purely imaginary algebraically special modes, but also for purely imaginary modes satisfying the quasinormal boundary condition.

Here we represented the improved user-friendly Mathematica® procedure for finding quasinormal modes of a wide class of black holes with the help of the Bernstein spectral method. The code we are sharing is automatic at a great extent and can be used for a wide class of wave equations allowing for asymptotically flat, de Sitter, or anti-de Sitter boundary conditions. By applying the method to the charged scalar field perturbations around a charged dilatonic asymptotically de Sitter black hole, we show that the method is very efficient when studying the regime of instability. We have checked our results with the help of the (less economic) time-domain integration method and achieved excellent agreement between both methods. We have shown that the instability for a dilaton-de Sitter black hole has superradiant nature, in a similar way with the Reissner-Nordström-de Sitter case.

The difficulty when using the Bernstein spectral method is related to the great number of roots of the matrices which are the bigger, the more overtones one wish to find, or the higher accuracy of the particular modes are required. Then, the question is which roots represent quasinormal frequencies and which are the numerical artifacts. Our code includes the comparison for matrices of various sizes, allowing to discard the numerical artifacts. The code is automatic and does not require essential modifications when applying it to this or that black-hole perturbation equations.

The Bernstein spectral method could be a relatively quick way to detect instability of black holes or find purely imaginary modes. In principle, it could be extended to a system of chained wave equations and at least to some axially symmetric black holes. We believe that future publications will study these questions.