A matrix method for quasinormal modes: Schwarzschild black holes in asymptotically flat and (anti-) de Sitter spacetimes

A black hole, long considered to be a physical as well as mathematical curiosity, is derived in general relativity as a generic prediction. Through gravitational collapse, a stellar-mass black hole can be formed at the end of the life cycle of a very massive star, when its gravity overcomes the neutron degeneracy pressure. A crucial feature of a black hole is the existence of the event horizon, a boundary in spacetime beyond which events cannot affect an outside observer. Despite its invisible interior, however, the properties of a black hole can be inferred through its interaction with other matter. By quantum field theory in curved spacetime, it is shown that the event horizons emit Hawking radiation, with the same spectrum of black body radiation at a temperature determined by its mass, charge and angular momentum [2]. The latter completes the formulation of black hole thermodynamics [3], which describes the properties of a black hole in analogy to those of thermodynamics by relating mass to energy, area to entropy, and surface gravity to temperature. Quasinormal modes (QNMs) arise as the temporal oscillations owing to perturbations in black hole spacetime [5]. Owing to the energy loss through flux conservation, these modes are not normal. Consequently, when writing the oscillation in an exponential form, $\exp \left(-\text{i}\omega t\right)$, the frequency of the modes is a complex number. The real part, ${{\omega}_{R}}$, represents the actual temporal oscillation; and the imaginary part, ${{\omega}_{I}}$, indicates the decay rate. Therefore, these modes are commonly referred to as quasinormal. The stability of the black hole spacetime guarantees that all small perturbation modes are damped. Usually, QNMs can be conditionally divided into three stages. The first stage involves a short period of the initial outburst of radiation, which is sensitively dependent on the initial conditions. The second stage consists of a long period dominated by the quasinormal oscillations, where the amplitude of the oscillation decays exponentially in time. This stage is characterized by only a few parameters of the black holes, such as their mass, angular momentum, and charge. The last stage takes place when the QNMs are suppressed by power-law or exponential late-time tails. The properties of QNMs have been investigated in the context of the AdS/CFT correspondence [6, 7]. As a matter of fact, practically every stellar object oscillates, and oscillations produced by very compact stellar objects and their detection are of vital importance in physics and astrophysics. In 2015, the first observation of gravitational waves from a binary black hole merger was reported [4]. The observation provides direct evidence of the last remaining unproven prediction of general relativity and reconfirms its prediction of space-time distortion on the cosmic scale.

Mathematically, the QNMs are governed by the linearized equations of general relativity constraining perturbations around a black hole solution. The resulting master field equation is a linear second order partial differential equation. Due to the difficulty in finding exact solutions to most problems of interest, various approximate methods have been proposed [8]. If the inverse potential, which can be viewed as a potential well, furnishes a well-defined bound state problem, the QNMs can be evaluated by solving the associated Schrödinger equation. In particular, when a smooth potential well can be approximated by the Pöschl-Teller potential, QNM frequencies can be obtained through the known bound states [9]. For general potential function, approaches such as continued fraction method [10], Horowitz and Hubeny (HH) method [6], asymptotic iteration method [11] can be utilized. A common feature of the above methods is that the corresponding master field equation is obtained by representing the wave function with power series. Higher precision is therefore achieved by considering higher order expansions. A semi-analytic technique to obtain the low-lying QNMs is based on a matching of the asymptotic WKB solutions at spatial infinity and on the event horizon [12]. The WKB formula has been extended to the sixth order [13]. Further generalization to a higher order, however, is not straightforward. Finite difference method is developed to numerically integrate the master field equation [14], and the temporal evolution of the perturbation can be obtained.

In this work, by discretizing the linear partial different equation [1], we transfer the master field equation as well as its boundary conditions into a homogeneous matrix equation. In our approach, the master field equation is presented in terms of linear equations describing N discretized points where the wave function is expanded up to Nth order for each of these points. This leads to a non-standard eigenvalue problem and can be solved numerically for the quasinormal frequencies. The present paper is organized as follows. In the next section, we briefly review how to reformulate the master field equation in terms of a matrix equation of non-standard eigenvalue problem. In sections 3 to 5, we investigate the quasinormal modes of Schwarzschild, Schwarzschild de Sitter and Schwarzschild anti- de Sitter black hole spacetime respectively. The precision and efficiency of the present approach are studied by comparing to the results obtained by other methods. Discussions and speculations are given in the last section.

Recently, we proposed a non-grid-based interpolation scheme which can be used to solve the eigenvalue problem [1]. A key step of the method is to formally discretize the unknown eigenfunction in order to transform a differential equation as well as the boundary conditions into a homogeneous matrix equation. Based on the information about N scattered data point, Taylor series are carried out for the unknown eigenfunction up to Nth order for each discretized point. Then the resulting homogeneous system of linear algebraic equations is solved for the eigenvalue. Here, we briefly describe the discretization procedure. For a univariate function f(x), one applies the Taylor expansion of a function to N  −  1 discrete points in a small vicinity of another point. Without loss of generality, let us expand the function about x₂ to ${{x}_{1}},{{x}_{3}},{{x}_{4}},\cdots,{{x}_{N}}$, and therefore obtains N  −  1 linear relations between function values and their derivatives as follows

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta \mathcal{F}=MD, \end{array}\end{eqnarray} \tag{ 2.1 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta \mathcal{F}={{\left(\,f\left({{x}_{1}}\right)-f\left({{x}_{2}}\right),f\left({{x}_{3}}\right)-f\left({{x}_{2}}\right),\cdots,f\left({{x}_{j}}\right)-f\left({{x}_{2}}\right),\cdots,f\left({{x}_{N}}\right)-f\left({{x}_{2}}\right)\right)}^{T}}, \end{array}\end{eqnarray} \tag{ 2.2 }$$

$$\begin{eqnarray}M=\left(\begin{array}{c} {{x}_{1}}-{{x}_{2}} & \frac{{{\left({{x}_{1}}-{{x}_{2}}\right)}^{2}}}{2} & \cdots & \frac{{{\left({{x}_{1}}-{{x}_{2}}\right)}^{i}}}{i!} & \cdots & \frac{{{\left({{x}_{1}}-{{x}_{2}}\right)}^{N-1}}}{(N-1)!} \\ {{x}_{3}}-{{x}_{2}} & \frac{{{\left({{x}_{3}}-{{x}_{2}}\right)}^{2}}}{2} & \cdots & \frac{{{\left({{x}_{3}}-{{x}_{2}}\right)}^{i}}}{i!} & \cdots & \frac{{{\left({{x}_{3}}-{{x}_{2}}\right)}^{N-1}}}{(N-1)!} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ {{x}_{j}}-{{x}_{2}} & \frac{{{\left({{x}_{j}}-{{x}_{2}}\right)}^{2}}}{2} & \cdots & \frac{{{\left({{x}_{j}}-{{x}_{2}}\right)}^{i}}}{i!} & \cdots & \frac{{{\left({{x}_{j}}-{{x}_{2}}\right)}^{N-1}}}{(N-1)!} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ {{x}_{N}}-{{x}_{2}} & \frac{{{\left({{x}_{N}}-{{x}_{2}}\right)}^{2}}}{2} & \cdots & \frac{{{\left({{x}_{N}}-{{x}_{2}}\right)}^{i}}}{i!} & \cdots & \frac{{{\left({{x}_{N}}-{{x}_{2}}\right)}^{N-1}}}{(N-1)!} \end{array}\right),\end{eqnarray} \tag{ 2.3 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} D={{\left(\,{{f}^{\prime}}\left({{x}_{2}}\right),{{f}^{\prime \prime}}\left({{x}_{2}}\right),\cdots,{{f}^{(k)}}\left({{x}_{2}}\right),\cdots,{{f}^{(N)}}\left({{x}_{2}}\right)\right)}^{T}}. \end{array}\end{eqnarray} \tag{ 2.4 }$$

Now, the above equation implies that all the derivatives at x  =  x₂ can be expressed in terms of the function values by using the Cramer's rule. In particular, we have

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{f}^{\prime}}\left({{x}_{2}}\right)=\det \left({{M}_{1}}\right)/\det (M), \\ {{f}^{\prime \prime}}\left({{x}_{2}}\right)=\det \left({{M}_{2}}\right)/\det (M), \end{array}\end{eqnarray} \tag{ 2.5 }$$

where M_i is the matrix formed by replacing the ith column of M by the column vector $\Delta \mathcal{F}$. Now, by permuting the N points, ${{x}_{1}},{{x}_{2}},\cdots,{{x}_{N}}$, we are able to rewrite all the derivatives at the above N points as linear combinations of the function values at those points. Substituting the derivatives into the eigenequation, one obtains N equations with $f\left({{x}_{1}}\right),\cdots,f\left({{x}_{N}}\right)$ as its variables. It was shown [1] that the boundary conditions can be implemented by properly replacing some of the above equations. Usually, the equations which are closer to the boundary of the problem are chosen to be replaced, since those equations are the least precise ones. For instance, in the case of asymptotically flat Schwarzschild spacetime below, we choose to replace the first and the last line in the matrix equation and implement the boundary condition by replacing equation (3.11) with equation (3.13).

Now we apply the above method to investigate the master field equation of for QNM. For simplicity, here we only investigate the scalar perturbation in black hole spacetime. According to the action of the massless scalar field with minimal coupling in curved four dimensional spacetime:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} S={\int}^{}{{\text{d}}^{4}}x\sqrt{-g}\mathcal{L}={\int}^{}{{\text{d}}^{4}}x\sqrt{-g}\left({{\partial}_{\mu}} \Phi {{\partial}^{\mu}} \Phi \right), \end{array}\end{eqnarray} \tag{ 2.6 }$$

the equation of motion for the massless scalar field reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{g}^{\mu \nu}}{{\nabla}_{\mu}}{{\nabla}_{\nu}} \Phi =0. \end{array}\end{eqnarray} \tag{ 2.7 }$$

Consider the following static spherical metric

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{d}{{s}^{2}}=-F(r)\text{d}{{t}^{2}}+\frac{\text{d}{{r}^{2}}}{F(r)}+{{r}^{2}}\left(\text{d}{{\theta}^{2}}+{{\sin}^{2}}\theta \text{d}{{\varphi}^{2}}\right), \end{array}\end{eqnarray} \tag{ 2.8 }$$

and rewriting the scalar field by using the separation of variables $\Phi =\frac{\phi (r)}{r}Y(\theta ){{\text{e}}^{-\text{i}\omega t+\text{i}m\varphi}}$, we obtain the following well-known Schrödinger-type equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{\text{d}}^{2}}\phi}{\text{d}r_{\ast}^{2}}+\left[{{\omega}^{2}}-V(r)\right]\phi =0 \end{array}\end{eqnarray} \tag{ 2.9 }$$

where $V(r)=F(r)\left(\frac{{{F}^{\prime}}(r)}{r}+\frac{L(L+1)}{{{r}^{2}}}\right)$ is the effective potential, and ${{r}_{\ast}}={\int}^{}\frac{\text{d}r}{F(r)}$ is tortoise coordinate. As discussed below, the boundary conditions in asymptotically flat, de Sitter and anti-de Sitter spacetimes are different. For the interpolation in equation (2.5) to be valid, appropriate coordinate transformation shall be introduced, which will be discussed in detail in the following sections.

In Schwarzschild spacetime, one has

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} F(r)=1-\frac{2M}{r}, \end{array}\end{eqnarray} \tag{ 3.1 }$$

and r_h  =  2M corresponds to the event horizon of the black hole. The potential vanishes on the horizon F(r_h)  =  0 and at infinity $r\to \infty$, therefore, the wave function has the asymptotic solution $\phi (r)\sim {{\text{e}}^{\pm \text{i}\omega {\int}^{}\frac{\text{d}r}{F(r)}}}$, where  ±  correspond to the wave travelling in positive and negative direction respectively. Since the wave function must be an ingoing wave on the horizon and an outgoing wave at infinity, the boundary conditions of equation (2.9) read

$$\begin{eqnarray}\begin{array}{*{35}{l}} \phi (+\infty ) & \sim & {{\text{e}}^{\text{i}\omega {{r}_{\ast}}}}, \\ \phi \left({{r}_{h}}\right) & \sim & {{\text{e}}^{-\text{i}\omega {{r}_{\ast}}}}, \end{array}\end{eqnarray} \tag{ 3.2 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{r}_{\ast}}=r+{{r}_{h}}\ln \left(\frac{r}{{{r}_{h}}}-1\right) \end{array}\end{eqnarray} \tag{ 3.3 }$$

is the tortoise coordinate. We study the QNM only in the region $r\in \left[{{r}_{h}},\infty \right)$. By taking into account the above boundary conditions, we first make use of the coordinate transformation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} x=1-\frac{{{r}_{h}}}{r}, \end{array}\end{eqnarray} \tag{ 3.4 }$$

and rewrite the scalar wave function as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \phi ={{\text{e}}^{\frac{\text{i}\omega {{r}_{h}}}{1-x}}}{{(1-x)}^{-\text{i}\omega {{r}_{h}}}}{{x}^{-\text{i}\omega {{r}_{h}}}}R(x). \end{array}\end{eqnarray} \tag{ 3.5 }$$

In this case, the boundary conditions become R(0)  =  R₀ and R(1)  =  R₁, where R₀ and R₁ are indeterminate constants.

The boundary conditions can be further simplified by introducing

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \chi (x)=x(1-x)R(x), \end{array}\end{eqnarray} \tag{ 3.6 }$$

so that

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \chi (1)=\chi (0)=0. \end{array}\end{eqnarray} \tag{ 3.7 }$$

As will be seen below, the boundary condition in equation (3.7) guarantees that the resulting matrix equation is homogeneous. The corresponding field equation now becomes

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tau}_{0}}(x){{\chi}^{\prime \prime}}(x)+{{\lambda}_{0}}(x){{\chi}^{\prime}}(x)+{{s}_{0}}(x)\chi (x)=0 \end{array}\end{eqnarray} \tag{ 3.8 }$$

where

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\tau}_{0}}(x) & = & {{x}^{2}}{{(1-x)}^{2}}{{A}_{2}}(x), \\ {{\lambda}_{0}}(x) & = & x(x-1)\left[x(x-1){{A}_{1}}(x)+2(1-2x){{A}_{2}}(x)\right], \\ {{s}_{0}}(x) & = & 2{{A}_{2}}(x)+x(x-1)\left[x(x-1){{A}_{0}}(x)+(1-2x){{A}_{1}}(x)+6{{A}_{2}}(x)\right] \end{array}\end{eqnarray} \tag{ 3.9 }$$

with

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{A}_{2}}(x) & = & -x{{(1-x)}^{2}}, \\ {{A}_{1}}(x) & = & 4\text{i}M\omega \left(2{{x}^{2}}-4x+1\right)-(1-3x)(1-x), \\ {{A}_{0}}(x) & = & (1-x)\left(1-8\text{i}M\omega \right)+16{{M}^{2}}{{\omega}^{2}}(x-2)+L(L+1). \end{array}\end{eqnarray} \tag{ 3.10 }$$

Now, we discretize the interval $x\in [0,1]$ by introducing N evenly distributed points with x₁  =  0 and x_N  =  1. By equation (2.5), one may rewrite the above partial different equation in a matrix form:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\mathcal{M}}\mathcal{F}=0, \end{array}\end{eqnarray} \tag{ 3.11 }$$

where $\mathcal{F}=\left(\,{{f}_{1}},{{f}_{2}},\cdots,{{f}_{i}},\cdots,{{f}_{N}}\right){{}^{T}}$ with ${{f}_{i}}=\chi \left({{x}_{i}}\right)$, and the matrix $\bar{\mathcal{M}}$ is a function of the quasinormal frequency, ω. The boundary conditions ${{f}_{1}}={{f}_{N}}=0$ can be implemented by defining

$$\begin{eqnarray}{{\mathcal{M}}_{k,i}}=\left\{\begin{array}{c} {{\delta}_{k,i}}, & k=1~\text{or}~N, \\ {{{\bar{\mathcal{M}}}}_{k,i}}, & k=2,3,\cdots,N-1, \end{array}\right. \end{eqnarray} \tag{ 3.12 }$$

and replacing equation (3.11) with

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{M}\mathcal{F}=0. \end{array}\end{eqnarray} \tag{ 3.13 }$$

Equation (3.13) furnishes a non-standard eigenvalue problem, obtained by discretizing the master equation for massless scalar field equation (2.8) in Schwarzschild black hole spacetime. As a homogeneous matrix equation, for eigenvalues $\omega ={{\omega}_{0}}$, the determinant

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \det \left(\mathcal{M}\left({{\omega}_{0}}\right)\right)=0. \end{array}\end{eqnarray} \tag{ 3.14 }$$

Equation (3.14) is the desired algebraic equation for the quasinormal frequencies, which can be solved numerically by using, for example, Mathematica. In table 1, we show the calculated values of the quasinormal frequencies, which are compared to those obtained by sixth order WKB method. It is inferred from the results that the present method is consistent with the WKB method.

In order to show that the present method gives convergent therefore reliable results, we show in figure 1 the calculated frequencies, as well as relative errors as functions of the number of interpolation points N. It is found that the results indeed converge well at big N. We note that the algebraic equation (3.14), by construction, usually has a finite number of roots. Therefore, one has to identify the correct eigenvalues from other by-products of the numerical solution. In order to seek out the relevant eigenvalue corresponding to the quasinormal mode, one may first study a simpler case by taking a specific limit of the metric parameter where the corresponding quasinormal frequency is already known. After pinning down the relevant quasinormal frequency, one can then vary continuously the parameters to restore the general case in question. In practice, the relevant eigenvalue of equation (3.14) is likely to be the one with the smallest imaginary part.

Zoom In Zoom Out Reset image size

In Schwarzschild de Sitter spacetime, one has

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} F(r)=1-\frac{2M}{r}- \Lambda {{r}^{2}}= \Lambda \left(1-\frac{{{r}_{h}}}{r}\right)\left({{r}_{c}}-r\right)\left(r+{{r}_{c}}+{{r}_{h}}\right). \end{array}\end{eqnarray} \tag{ 4.1 }$$

Here $M=\frac{{{r}_{h}}{{r}_{c}}\left({{r}_{h}}+{{r}_{c}}\right)}{2\left(r_{h}^{2}+{{r}_{h}}{{r}_{c}}+r_{c}^{2}\right)}$ and $\Lambda =\left(r_{h}^{2}+{{r}_{h}}{{r}_{c}}+r_{c}^{2}\right){{}^{-1}}$ are the black hole mass and the cosmological constant, where r  =  r_h and r  =  r_c represent the event horizon and cosmological horizon respectively. The boundary conditions, owing to the existence of the two horizons, become

$$\begin{eqnarray}\begin{array}{*{35}{l}} \phi \left({{r}_{c}}\right) & \sim & {{\text{e}}^{\text{i}\omega {{r}_{\ast}}}}, \\ \phi \left({{r}_{h}}\right) & \sim & {{\text{e}}^{-\text{i}\omega {{r}_{\ast}}}}. \end{array}\end{eqnarray} \tag{ 4.2 }$$

We study QNM in the radial interval ${{r}_{h}}\leqslant r\leqslant {{r}_{c}}$, and the the tortoise coordinate reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{r}_{\ast}}={{\eta}_{h}}\ln \left(r-{{r}_{h}}\right)+{{\eta}_{c}}\ln \left({{r}_{c}}-r\right)+{{\eta}_{i}}\ln \left(r+{{r}_{h}}+{{r}_{c}}\right), \end{array}\end{eqnarray} \tag{ 4.3 }$$

where

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\eta}_{h}} & = & \frac{{{r}_{h}}}{ \Lambda }{{\left({{r}_{h}}-{{r}_{c}}\right)}^{-1}}{{\left(2{{r}_{h}}+{{r}_{c}}\right)}^{-1}}, \\ {{\eta}_{c}} & = & \frac{{{r}_{c}}}{ \Lambda }{{\left(r_{h}^{2}+{{r}_{h}}{{r}_{c}}-2r_{c}^{2}\right)}^{-1}}, \\ {{\eta}_{i}} & = & \frac{{{r}_{h}}+{{r}_{c}}}{ \Lambda }{{\left(2{{r}_{h}}+{{r}_{c}}\right)}^{-1}}{{\left(2{{r}_{h}}+{{r}_{c}}\right)}^{-1}}. \end{array}\end{eqnarray} \tag{ 4.4 }$$

In order to transfer the radial interval into [0, 1], we introduce the coordinate transformation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} y=\frac{r-{{r}_{h}}}{{{r}_{c}}-{{r}_{h}}}, \end{array}\end{eqnarray} \tag{ 4.5 }$$

which gives y  =  0 at r  =  r_h and y  =  1 at r  =  r_c. In accordance with the boundary conditions, we rewrite the scalar field as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \phi ={{(1-y)}^{\text{i}\omega {{\eta}_{c}}}}{{y}^{-\text{i}\omega {{\eta}_{h}}}}G(\,y). \end{array}\end{eqnarray} \tag{ 4.6 }$$

This implies that G(0)  =  G₀ and G(1)  =  G₁, where G₀ and G₁ are indeterminate constants. Again, in order to transfer the boundary conditions into the desired form, we further introduce

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \zeta (\,y)=y(1-y)G(\,y), \end{array}\end{eqnarray} \tag{ 4.7 }$$

and transfer the field equation to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tau}_{0}}(\,y){{\zeta}^{\prime \prime}}(\,y)+{{\lambda}_{0}}(\,y){{\zeta}^{\prime}}(\,y)+{{s}_{0}}(\,y)\zeta (\,y)=0 \end{array}\end{eqnarray} \tag{ 4.8 }$$

with

$$\begin{eqnarray}&&{{\tau}_{0}}(\,y)= \Lambda {{\left({{r}_{c}}-{{r}_{h}}\right)}^{2}}y(1-y)\frac{(\,y-2){{r}_{h}}-(1+y){{r}_{c}}}{(\,y-1){{r}_{h}}-y{{r}_{c}}},\\ &&{{\lambda}_{0}}(\,y)=\frac{(\,y-1)y\tau _{0}^{\prime}(\,y)+2\text{i}{{\tau}_{0}}(\,y)\left(\,y\omega {{\eta}_{c}}+{{\eta}_{h}}(\omega -y\omega )+\text{i}(2y-1)\right)}{(\,y-1)y},\\ &&{{s}_{0}}(\,y)=\frac{{{\tau}_{0}}(\,y)}{{{(\,y-1)}^{2}}{{y}^{2}}}\left\{\frac{\text{i}(\,y-1)y\tau _{0}^{\prime}(\,y)\left(\,y\omega {{\eta}_{c}}-(\,y-1)\omega {{\eta}_{h}}+2\text{i}y-\text{i}\right)}{{{\tau}_{0}}(\,y)}\right. \\ &&+\,\frac{{{(\,y-1)}^{2}}{{y}^{2}}\left({{r}_{c}}-{{r}_{h}}\right){{}^{2}}\left({{\omega}^{2}}-V(\,y)\right)}{{{\tau}_{0}}{{(\,y)}^{2}}}+6(\,y-1)y+2\\ &&\left. +\,\omega \left[(\,y-1){{\eta}_{h}}\left(2y\omega {{\eta}_{c}}+5\text{i}y-3\text{i}\right)-y{{\eta}_{c}}\left(\,y\omega {{\eta}_{c}}+5\text{i}y-2\text{i}\right)-{{(\,y-1)}^{2}}\omega \eta _{h}^{2}\right]\right\},\\ &&V(\,y)={{\tau}_{0}}(\,y)\left(\frac{\frac{\tau _{0}^{\prime}(\,y)}{{{r}_{c}}-{{r}_{h}}}}{{{r}_{h}}+\left({{r}_{c}}-{{r}_{h}}\right)y}+\frac{L(L+1)}{{{\left({{r}_{h}}+\left({{r}_{c}}-{{r}_{h}}\right)y\right)}^{2}}}\right).\end{eqnarray} \tag{ 4.9 }$$

The boundary conditions now read $\zeta (0)=\zeta (1)=0$. Now we are in the position to utilize the same numerical procedure to discretize the wave function in the interval $0\leqslant y\leqslant 1$, and solve for the quasinormal frequencies. The obtained the quasinormal frequencies are presented in table 2 compared to those obtained by the WKB method. It is found that the results from the present method are consistent with those from the WKB method.

Finally, we study the Schwarzschild Anti-de Sitter spacetime. We have

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} F(r)=1-\frac{2M}{r}+ \Lambda {{r}^{2}}= \Lambda \left(1-\frac{{{r}_{h}}}{r}\right)\left({{r}^{2}}+{{r}_{h}}r+r_{h}^{2}+\frac{1}{ \Lambda }\right). \end{array}\end{eqnarray} \tag{ 5.1 }$$

Here $M=\frac{1}{{{r}_{h}}}\left(1+ \Lambda r_{h}^{2}\right)$ is the mass of the black hole. Equation (5.1) implies that the effective potential V(r) defined below equation (2.9) diverges at infinity [6], and therefore the wave function vanishes at infinity. Following [6], one utilizes the coordinate transformation v  =  t  +  r_*. The resulting black hole metric reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{d}{{s}^{2}}=-F(r)\text{d}{{v}^{2}}+2\text{d}v\text{d}r+{{r}^{2}}\left(\text{d}{{\theta}^{2}}+{{\sin}^{2}}\theta \text{d}{{\varphi}^{2}}\right). \end{array}\end{eqnarray} \tag{ 5.2 }$$

The scalar field equation becomes

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} F(r)\frac{{{\text{d}}^{2}}\phi (r)}{\text{d}{{r}^{2}}}+\left[\frac{\text{d}F(r)}{\text{d}r}-2\text{i}\omega \right]\frac{\text{d}\phi (r)}{\text{d}r}-U(r)\phi (r)=0, \end{array}\end{eqnarray} \tag{ 5.3 }$$

where $U(r)=\frac{1}{r}\frac{\text{d}F(r)}{\text{d}r}+\frac{L(L+1)}{{{r}^{2}}}$. We introduce the coordinate transformation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} z={{r}_{h}}/r, \end{array}\end{eqnarray} \tag{ 5.4 }$$

to transfer the radial coordinate into the interval $z\in [0,1]$, with z  =  1 at r  =  r_h and z  =  0 at $r\to \infty$. On the event horizon, the boundary condition reads $\phi {{|}_{z=1}}={{\phi}_{0}}$, where ${{\phi}_{0}}$ is a constant. By further introducing

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \varrho (z)=(z-1)\phi (z), \end{array}\end{eqnarray} \tag{ 5.5 }$$

one obtains the desired boundary condition $\varrho (0)=\varrho (1)=0$ for the field equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tau}_{0}}(z){{\varrho}^{\prime \prime}}(z)+{{\lambda}_{0}}(z){{\varrho}^{\prime}}(z)+{{s}_{0}}(z)\varrho (z)=0, \end{array}\end{eqnarray} \tag{ 5.6 }$$

with

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\tau}_{0}}(z) & = & r_{h}^{2}\left(\frac{ \Lambda }{{{z}^{2}}}- \Lambda z\right)-z+1, \\ {{\lambda}_{0}}(\,y) & = & \frac{2\text{i}\omega {{r}_{h}}+ \Lambda \left(-{{z}^{2}}+2z+2\right)r_{h}^{2}+(2-z)z}{{{z}^{2}}}, \\ {{s}_{0}}(z) & = & \frac{r_{h}^{2}\left[\Lambda {{z}^{2}}\left({{z}^{2}}-2z-2\right)-(z-1)U\left(\frac{r0}{z}\right)\right]-2\text{i}\omega {{z}^{2}}{{r}_{h}}+(z-2){{z}^{3}}}{(z-1){{z}^{4}}}. \end{array}\end{eqnarray} \tag{ 5.7 }$$

Again, the eigenequation can be obtained following the same procedures as before. The calculated quasinormal frequencies are shown in table 3 and compared with the results obtained by the HH method. It is inferred from the results that the present method is in accordance with the HH method.

Table 3. The quasinormal frequencies in asymptotically Anti de Sitter black hole spacetime obtained by the present method. The interpolation makes use of 22 points. It is compared to those obtained by HH method. Both calculations consider L  =  0 and n  =  0.

$\left(\Lambda,{{r}_{h}}\right)$	ω (HH method [6])	ω (present method)
{1, 100}	184.9534  −  266.3856i	184.956  −  266.385i
{1, 50}	92.4937  −  133.1933i	92.4949  −  133.193i
{1, 10}	18.6070  −  26.6418i	18.6073  −  26.6417i
{1, 5}	9.4711  −  13.3255i	9.471 29  −  13.3255i
{1, 1}	2.7982  −  2.6712i	2.797 78  −  2.670 47i
{1, 0.8}	2.5878  −  2.1304i	2.586 24  −  2.128 76i
{1, 0.6}	2.4316  −  1.5797i	2.425 92  −  1.572 12i
{1, 0.4}	2.3629  −  1.0064i	2.381 52  −  0.938 149i

In this work, we proposed a new interpolation scheme to discretize the master field equation for the scalar quasinormal modes. It is shown that the method can be applied to different black hole spacetimes. By appropriately introducing coordinate transformations, the resulting homogeneous matrix equation possesses very similar characteristics. And therefore, the quasinormal frequencies can be obtained by the same numerical solver for algebraic equations.

On the one hand, we obtain the desired the boundary conditions through appropriate choice of coordinate, so that quasinormal modes for different black hole spacetimes are obtained through the same numerical scheme. The precision of the present method, on the other hand, can be easily improved by increasing the total number of discretization points N, which is a convenient feature. By taking the advantage of the efficiency of the existing matrix as well as algebraic equation solvers, such as Matlab and Mathematica, the present method is practical and efficient. In particular, we have deliberately transferred the radial variable into the interval [0, 1]. This is because the evaluation of equation (2.5) can be quite time-consuming for a high-rank matrix. However, once the radial interval is given, such calculations become independent of the specific form of metric. As a result, equation (2.5) can be carried out beforehand which in turn increases the efficiency of the present method.

In this work, we studied quasinormal frequencies. We note that the corresponding wave function can also be obtained easily by substituting the obtained frequency ${{\omega}_{A}}$ into $\mathcal{M}(\omega )$ and numerically evaluating the column matrix ${{\xi}_{A}}$ which satisfies $\mathcal{M}\left({{\omega}_{A}}\right){{\xi}_{A}}=0$. For instance, this can be achieved by using Eigensystem command of Mathematica to acquire the eigenvector ${{\xi}_{A}}$ corresponding to the null eigenvalue.

It is noting that the present method is particularly advantageous when applied to asymptotically AdS spacetime. This is because in this case, the derivation of equation (5.6) does not require the knowledge of an analytic form of the tortoise coordinates, which also applies to other metrics in AdS spacetime. In the asymptotically flat as well as dS spacetime, on the other hand, an analytic form of the tortoise coordinate usually provides considerable convenience to acquire the desired boundary conditions, in order that the problem can be transformed into a homogeneous matrix equation for the quasinormal modes. In fact, the above mathematical difficulties are also encountered for other approaches such as continued fraction and asymptotic iteration methods. Furthermore, we observe that the proposed method is quite general and can be employed to investigate more sophisticated and physically interesting cases. It possesses flexibility and therefore the potential to explore some black hole metrics where the applications of other traditional methods become less straightforward. As an example, the quasinormal modes of a rotational black hole is characterized by, besides the quasinormal frequency ω, a second eigenvalue λ whose physical content is associated with the angular quantum number L. Its numerical solution, therefore, involves finding the two eigenvalues, ω and λ, simultaneously. The continued fraction method is fit for the task, but its success relies on the derivation of a recurrence relation of the coefficients, which might not be obvious for some sophisticated metrics. A preliminary attempt [15] shows that the approach proposed in this work, on the other hand, can be applied straightforwardly in a more intuitive fashion. Another example is the quasinormal modes of massive Dirac field [16], where the present method is also expected to introduce significant convenience when handling the coupled equations of spinor components.

In practice, it is found that the most time-consuming part the calculation is to invert the matrix using equation (2.5). The efficiency is closely related to the existing algorithm for algebriac equation solvers, such as Mathematica and Matlab. However, this procedure is identical for most calculations once the grid points are fixed, therefore the efficiency of the method can be significantly improved for such similar problems. The HH method usually is carried out up to 40th order to achieve the precision presented in this work, while the precision of WKB method cannot be adjusted freely. The matrix method achieve the same precision by making use of 22 grid points. Furthermore, we observe that the proposed method is quite general and can be employed to investigate more sophisticated and physically interesting cases. It possesses flexibility and therefore the potential to explore some black hole metrics where the applications of other traditional methods become less straightforward. As an example, the quasinormal modes of a rotational black hole is characterized by, besides the quasinormal frequency ω, a second eigenvalue λ whose physical content is associated with the angular quantum number L. Its numerical solution, therefore, involves finding the two eigenvalues, ω, and λ, simultaneously. A preliminary attempt [15] shows that the approach proposed in this work can be applied straightforwardly in an intuitive fashion. Another example is the quasinormal modes of massive Dirac field [16] or massive vector field [17]. Here, the mathematical difficulties mostly stem from the coupled equations of spinor components. The solution of the problem was achieved, and a vital step was to appropriately introduce new variables and rewrite the master equations as decoupled equations (see [16–18]). For our present method, however, the coupled master equations can be handled directly which consequently lead to nothing but a bigger system of linear algebraic equations, the latter does not imply any further mathematical difficulty. These are worthy topics for further investigations.

We gratefully acknowledge the financial support from Brazilian funding agencies Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), as well as National Natural Science Foundation of China (NNSFC) under contract No.11573022 and 11375279.