HTML
--> --> -->In general, a quasinormal mode is defined as an eigenmode of a dissipative system. In the context of general relativity, small perturbations of a black hole metric may lead to quasinormal modes [14-17]. Their importance is due to the fact that they are closely associated with the final fate of the black hole merger. Besides, quasinormal modes attracted much attention in recent years mainly due to the development of the holographic principle in the anti-de Sitter/conformal field theory (AdS/CFT) correspondence, which is widely recognized as a tool for exploring strongly coupled systems [18]. In other words, by studying the black hole perturbations in the bulk, one may extract important physical quantities, such as transport coefficients, e.g. viscosity, conductivity, and diffusion constants, of the dual system on the boundary.
The study of black hole quasinormal modes encompasses various types of perturbations, which include those of the scalar field, spinor field, vector field, or the metric itself. Revealed
From the mathematical point of view, the analysis of a quasinormal mode is related to a non-Hermitian eigenvalue problem for a system of coupled linear differential equations with associated boundary conditions. Apart of a few cases where analytic solutions for asymptotic quasinormal mode spectra can be obtained [19-21], one often resorts to semi-analytic techniques [22]. The P?shl-Theller potential approximation [23] replaces the effective potential in the master equation by a P?shl-Theller one, where the analytic eigenvalue is known. The WKB method [24, 25] is a semi-analytic approach proposed by Schutz et al, where the technique is applied specifically to the one-dimensional time-independent Schr?dinger-like equation for the quasi-eigenvalue problem. The 6th order formalism of the method was subsequently derived by Konoplya [26]. In practice, the application of the WKB method is rather straightforward, since it only requires the information about the background metric and the form of the effective potential. However, the uncertainty of the obtained result is somewhat undetermined beforehand, which is in part related to the specific shape of the effective potential. It is worth mentioning that, in the eikonal limit, the quasinormal modes of spherical black holes in asymptotically flat spacetime were found to be associated with the parameters of the circular null geodesic [27]. However, Konoplya and Stuchlik [28] showed analytically that the relationship is not valid for some particular cases, which has been further explored in [29-31]. In addition to the semi-analytic method, various numerical methods have also been proposed. The first type of numerical methods aims at evaluating the time-dependent wave equation directly for a given set of initial conditions. This is achieved by replacing the derivative with the finite difference [32-37]. Naturally, the resulting waveform contains the information about all three stages of temporal evolution of small perturbations. At the end of the calculation, one may further evaluate the quasinormal frequency by using Fourier analysis. A major disadvantage of the approach is that one cannot obtain the complete spectrum of quasinormal modes, due to the fact that usually only a few long-lived modes can be clearly identified. The second class of numerical schemes involves Fourier decomposition of the perturbed Einstein equations, while adopting appropriate boundary conditions. The latter assume that the resultant wave function corresponds to incoming waves on the horizon and outgoing at infinity. A series expansion of the wave function is usually carried out, and the method is typically iterative, namely, higher order results are derived from and improved using the lower order ones. Consequently, the precision of the method is controlled by the order of the expansion. These methods include the continued fraction method [38, 39], the Horowitz and Hubeny's (HH) method for the AdS black hole [40], and the asymptotic iteration method [41-44]. However, for the continued fraction method it is not so straightforward to derive the desired iterative relation for specific background metrics, and the HH method only applies to the asymptotically AdS spacetime. Recently, we proposed a matrix method [45-47] , where the spatial coordinate is discretized so that the differential equation, as well as its boundary conditions, are transformed into a homogeneous matrix equation. A vital feature of the method is that the eigenfunction is expanded in the vicinity of all grid points, and, therefore, the precision of the algorithm can be potentially improved. The method has been employed to study the quasinormal modes of the Schwarzschild black holes in asymptotically flat, (anti-)de Sitter (AdS/dS) spacetimes, as well as the Kerr and Kerr-Sen black hole metrics [46, 47]. In this survey, we further explore possible applications of the method to more general scenarios. We provide a detailed account of the procedure of how the algorithm is implemented. We argue that the implementation of the proposed method does not rely on the analytic form of the tortoise coordinate. Moreover, the method can be readily applied even to the case of the numerical black hole metric. Our discussion addresses black hole metrics of different asymptotic behavior, in particular where an analytic form of the tortoise coordinate is not available. Advantages, as well as improvements of the alternative approaches, are also addressed.
The paper is organized as follows. In the following section, we discuss the core algorithm of the matrix method. In Section III, we explore the applications of the method in various scenarios, including for different background metrics and the corresponding boundary conditions. In particular, we investigate the method of separation of variables for the master equation in asymptotically flat, dS, and AdS spacetimes, while considering different cases such as static, rotational, and extreme black holes. Moreover, the case where the master equation is a system of coupled ordinary differential equations is addressed. In addition, we provide a comprehensive discussion on the code implementation using Mathematica. Further discussions and concluding remarks are given in the last section of the paper.
Let us start by discussing the general strategy of applying the method to the standard scenario for quasinormal modes, where the master equation is of the Schr?dinger-type, given as follows:
$\frac{{{{\rm d}^2}}}{{{\rm d}r_*^2}}\Phi (r) + \left[ {{\omega ^2} - V(r)} \right]\Phi (r) = 0,$  | (2.1) |
In order to solve a differential equation, one also has to determine the associated boundary conditions. In the case of quasinormal modes, the boundary conditions are determined by the characteristics of the black hole at the horizon and infinity. As the apparent horizon of a black hole is a one-way membrane which prohibits the outgoing light rays to travel across it, the boundary conditions for the quasinormal modes are defined such that the solutions should be ingoing at the horizon and purely outgoing at infinity, namely
$\Phi \sim \left\{ {\begin{aligned} &{{{\rm e}^{{\rm i}\bar \omega {r_*}}}}&\!\!\!\!\!\!{r \to \infty }&\;\;{({\rm{for\ asymptotically \ flat \ spacetime}})}\\ &{{{\rm e}^{{\rm i}\omega {r_*}}}}&\!\!\!\!\!\!{r \to {r_{\rm c}}}&\;\;{({\rm{for \ asymptotically \ dS \ spacetime}})}\\ &{{{\rm e}^{ - {\rm i}\omega {r_*}}}}&\!\!\!\!\!\!{r \to {r_{\rm h}}}&\;\;{({\rm{for \ both \ spacetimes}})} \end{aligned}} \right.,$  | (2.2) |
Similarly to the other methods based on series expansion, we first change the domain of the radial variable to a finite range by introducing
${\cal H}(\omega ,z)R(z) = 0,$  | (2.3) |
${\rm with}~ R(z = 0) = {C_0}~{\rm{and}}~ R(z = 1) = {C_1}.$  | (2.4) |
$F(z) = R(z)z(1 - z),$  | (2.5) |
Subsequently, the master equation can be rewritten with F(z)
${\cal G}(\omega ,z)F(z) = 0,$  | (2.6) |
$F(z = 0) = F(z = 1) = 0.$  | (2.7) |
As the next step, we descretize the master equation Eq.(6) by introducing N interpolation points zi with
$ \bar{\cal M}{\cal F}=0, $  | (2.8) |
$ {\cal F}=(f_1,f_2,\cdots,f_i,\cdots,f_N)^T $  | (2.9) |
$ f_1=f_N=0. $  | (2.10) |
$ {\cal M}{\cal F}=0, $  | (2.11) |
${{\cal M}_{k,i}} = \left\{ {\begin{array}{*{20}{l}} {{\delta _{k,i}},}&{k = 1\;{\rm{or}}\;{{N}}}\\ {{{\bar {\cal M} }_{k,i}},}&{k = 2,3, \cdots ,N - 1} \end{array}} \right..$  | (2.12) |
$ \det\left({\cal M}(\omega)\right)\equiv\left|{\cal M}((\omega)\right|=0. $  | (2.13) |
In Refs. [46, 47], the above algorithm is employed to descretize the master equation, Eq.(6). Alternatively, one may also employ other numerical approaches for numerical differentiation, such as the differential quadrature method [48-50], and the Runge-Kutta method, among others. Regarding the finite difference formulas, if Newton's difference quotient is employed in Eq.(6), the resulting matrix
The following section is devoted to a more detailed discussion of various applications.
We first address the issue of different boundary conditions owing to different asymptotical properties of spacetime. Usually, spacetimes can be classified according to their asymptotical behavior at infinity, such as asymptotically flat, dS, and AdS spacetimes. It is known that the effective potential V(r) vanishes at infinity in the case of asymptotically flat spacetime. On the other hand, it diverges at infinity for asymptotically dS and AdS spacetimes. However, for the dS case, there exists a cosmological horizon rc between the observer and infinity, which is also a Cauchy one-way membrane. An observer staying between rh and rc will not receive any signals from the inside of event horizon, nor from the outside of cosmological horizon. Therefore, the physically valid domain is
For extreme black holes, the associated tortoise coordinate has a different feature near the black hole horizon than for non-extreme black holes. This is because, for an extreme black hole, one requires
For rotating black holes, an additional equation is implied, whose physical content is associated with the angular momentum quantum number. In this case, one usually gets a system of two coupled master equations with eigenvalues
In addition, one frequently encounters two difficult scenarios in realistic numerical studies. The first involves the difficulty in obtaining an analytic expression for the tortoise coordinate
2
3.1.Boundary conditions
In the following, we write down the explicit forms of the boundary conditions for five specific cases of asymptotic spacetime and its boundaries:Case 1. For the non-extreme black hole, near the horizon, if
$\begin{split} {r_*} = &\int {\displaystyle\frac{{{\rm d}r}}{{f(r)}}} = \int {\left( {\displaystyle\frac{1}{{{{\hat f}_1}(r - {r_{\rm h}})}} - \frac{{{{\hat f}_2}}}{{2\hat f_1^2}} + O} \right)} {\rm d}r\\ =& \displaystyle\frac{{\ln \left| {r - {r_{\rm h}}} \right|}}{{{{\hat f}_1}}} - \displaystyle\frac{{{{\hat f}_2}}}{{2\hat f_1^2}}(r - {r_{\rm h}}) + O, \end{split}$  | (3.1) |
$ \Phi\sim {\rm e}^{-{\rm i}\omega r_*}\sim B_a\equiv \left(r-r_{\rm h}\right)^{-\frac{i\omega}{\hat{f}_1}}. $  | (3.2) |
$\begin{split} {r_*} &= \int {\displaystyle\frac{{{\rm d}r}}{{f(r)}}} \\&= \int {\left( {\displaystyle\frac{2}{{{{\hat f}_2}{{(r - {r_{\rm h}})}^2}}} - \displaystyle\frac{{2{{\hat f}_3}}}{{3\hat f_2^2(r - {r_{\rm h}})}} + \displaystyle\frac{{4\hat f_3^2 - 3{{\hat f}_2}{{\hat f}_4}}}{{18\hat f_2^3}} + O} \right)} {\rm d}r\\ &= - \displaystyle\frac{2}{{{{\hat f}_2}(r - {r_{\rm h}})}} - \displaystyle\frac{{2{{\hat f}_3}\ln \left| {r - {r_{\rm h}}} \right|}}{{3\hat f_2^2}} + \displaystyle\frac{{4\hat f_3^2 - 3{{\hat f}_2}{{\hat f}_4}}}{{18\hat f_2^3}}(r - {r_{\rm h}}) + O, \end{split}$  | (3.3) |
$ \Phi\sim {\rm e}^{-{\rm i}\omega r_*}\sim B_b\equiv\left(r-r_{\rm h}\right)^{\frac{2i\omega \hat{f}_3}{3\hat{f}_2^2}}{\rm e}^{\frac{2i\omega}{\hat{f}_2(r-r_{\rm h})}}. $  | (3.4) |
${r_*} = \int {\frac{{{\rm d}r}}{{f(r)}}} = {\left. {\int {\left( {1 - A{r^\alpha } + O} \right)} {\rm d}r} \right|^{r \to \infty }} = r - \frac{A}{{1 + \alpha }}{r^{1 + \alpha }} + O.$  | (3.5) |
$ \Phi\sim {\rm e}^{{\rm i}\bar{\omega} r_*}\sim B_0\equiv {\rm e}^{{\rm i}\bar{\omega}\left(r-\frac{A}{1+\alpha}r^{1+\alpha}\right)},$  | (3.6) |
Case 4. In the asymptotically de Sitter spacetime, according to Case 1 and Case 2, one obtains similar results by replacing rh with rc. Therefore, following Eq.(15), near the non-extreme cosmological horizon, the boundary condition is
$ \Phi\sim {\rm e}^{i\omega r_*}\sim B_1\equiv\left(r_{\rm c}-r\right)^{\frac{i\omega}{\bar{f}_1}}. $  | (3.7) |
$ \Phi\sim {\rm e}^{{\rm i}\omega r_*}\sim B_2\equiv\left(r_{\rm c}-r\right)^{-\frac{2i\omega \bar{f}_3}{3\bar{f}_2^2}}{\rm e}^{\frac{2{\rm i}\omega}{\bar{f}_2(r_{\rm c}-r)}}, $  | (3.8) |
2
3.2.Separation of variables
By making use of the asymptotic behavior obtained above, we can write down the appropriate forms for the wave functions by the method of separation of variables.For the non-extreme black hole in asymptotically flat spacetime, we have
$ \Phi(r)={\rm e}^{{\rm i}\bar{\omega}\left(r-\frac{A}{1+\alpha}r^{1+\alpha}\right)}\left(r-r_{\rm h}\right)^{-\frac{i\omega}{\hat{f}_1}}r^{\frac{i\omega}{\hat{f}_1}}R(r), $  | (3.9) |
$\Phi (r) = {{\rm e}^{{\rm i}\bar \omega \left( {r - \frac{A}{{1 + \alpha }}{r^{1 + \alpha }}} \right)}}{\left( {r - {r_{\rm h}}} \right)^{\frac{{2i\omega {{\hat f}_3}}}{{3\hat f_2^2}}}}{{\rm e}^{\frac{{2{\rm i}\omega }}{{{{\hat f}_2}(r - {r_{\rm h}})}}}}{r^{ - \frac{{2i\omega {{\hat f}_3}}}{{3\hat f_2^2}}}}{{\rm e}^{ - \frac{{2{\rm i}\omega }}{{{{\hat f}_2}r}}}}R(r).$  | (3.10) |
$\Phi (r) = \left\{ {\begin{array}{*{20}{c}} {{B_a}{B_1}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{\not = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{\not = 0}}}\\ {{B_a}{B_2}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{\not = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{ = 0}}}\\ {{B_b}{B_1}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{ = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{\not = 0}}}\\ {{B_b}{B_2}R(r),}&{{\rm{when}}\;{{{\hat f}}_{\rm{1}}}{\rm{ = 0}}\;{\rm{and}}\;{{{\bar f}}_{\rm{1}}}{\rm{ = 0}}} \end{array}} \right..$  | (3.11) |
$ f(r)\frac{{\rm d}^2\Phi(r)}{{\rm d}r^2}+\left[\frac{{\rm d}f(r)}{{\rm d}r}-2i\omega\right]\frac{{\rm d}\Phi(r)}{{\rm d}r}-U(r)\Phi(r)=0, $  | (3.12) |
For the rotating black hole, one has to deal with the angular part of the master equation, which possesses the following form, similar to the radial part of the master equation:
$ (1-u^2)\frac{{\rm d}}{{\rm d}u}\left[(1-u^2)\frac{{\rm d}Y(u)}{{\rm d}u}\right]-{\cal W}(u,\omega,\lambda)Y(u)=0. $  | (3.13) |
$Y(u) \sim \left\{ {\begin{array}{*{20}{l}} {{{\left| {1 - u} \right|}^{ \pm \frac{{{W_ + }}}{2}}}}&{u \to 1}\\ {{{\left| {1 + u} \right|}^{ \pm \frac{{{W_ - }}}{2}}}}&{u \to - 1} \end{array}} \right..$  | (3.14) |
$Y(u) = {\left| {1 - u} \right|^{\left| {\frac{{{W_ + }}}{2}} \right|}}{\left| {1 + u} \right|^{\left| {\frac{{{W_ - }}}{2}} \right|}}S(u).$  | (3.15) |
Last but not least, we carry out another coordinate transformation to conveniently convert the domain of r and u into [0,1]. For the asymptotically flat and anti-de Sitter black hole spacetimes, the new coordinate is chosen as
By following the above procedure, the transformed master equation and its boundary conditions are given in Eq.(3) and Eq.(4). Subsequently, one may proceed with the calculations of the quasinormal frequency of the matrix equation as discussed in Section 2.
2
3.3.Master equation for a system of coupled equations
In the above discussion, we have assumed that the master equation can be expressed in terms of a system of decoupled equations, as we have treated the radial and angular parts of the master equation separately. In practice, such simplification might not always be possible. As an example, in the study of the quasinormal modes of a massive vector field, the resultant system of equations might be coupled. However, a system of coupled ordinary differential equations does not particularly pose a difficulty for the matrix method. This is because, in this case, one only needs to write the system of coupled equations in terms of a matrix equation of larger size. In other words, if a system of n equations cannot be decoupled, the complication formally manifests in the fact that one has to handle a matrix equation of the size $\begin{array}{l} {{\cal M}_a}{{\cal F}_a} = {{\cal B}_b}{{\cal F}_b},\quad {{\cal M}_b}{{\cal F}_b} = {{\cal B}_a}{{\cal F}_a}, \end{array}$  | (3.16) |
${{\cal M}_n}{{\cal F}_n} \equiv \left( {\begin{array}{*{20}{c}} {{{\cal M}_a}}&{ - {{\cal B}_b}}\\ { - {{\cal B}_a}}&{{{\cal M}_b}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\cal F}_a}}\\ {{{\cal F}_b}} \end{array}} \right) = 0,$  | (3.17) |
Alternatively, one can either rewrite Eq.(30) as
$ \left({\cal M}_a{\cal B}_a^{-1}{\cal M}_b-{\cal B}_b\right){\cal F}_b=0， $  | (3.18) |
$ \left({\cal M}_b{\cal B}_b^{-1}{\cal M}_a-{\cal B}_a\right){\cal F}_a=0. $  | (3.19) |
2
3.4.Code implementation
The code implementation has been carried out in Mathematica. It takes advantage of the efficiency of existing packages for matrix manipulation, as well as for solving nonlinear equations. The source codes for the Schwarzschild black holes in asymptotically flat, (anti-) de Sitter (AdS/dS) spacetimes, as well as for the Kerr and Kerr-Sen black hole metrics, have been released publicly [46, 47]. In this subsection, we discuss few technical aspects of the codes implemented in Mathematica.Firstly, Mathematica is not always very efficient in handling decimals, especially when decimals are further processed as arguments of internal functions, such as trigonometric or exponential functions. Therefore, in view of optimizing the computational time, one should avoid coding in decimals. As an example, one may use the following Macro to transform decimal output of a function “TU” into a fraction:
$ \begin{split} \rm{YOULI}[\rm{TU}\_] :=& \rm{Rationalize}[\rm{Chop}[\rm{Expand}\\ &[\rm{N}[\rm{TU},\rm{precs}]],10^{-\rm{precs}}], 10^{-\rm{precs}}], \end{split} $  | (3.20) |
In Mathematica, there are usually two methods of finding the roots of Eq.(13), which is essentially an algebraic equation. They are FindRoot and NSolve. The former is used to search for a root of a function starting from an initial point provided by the user, while the latter attempts to enumerate all existing roots of the function. For the first case, one has to find an appropriate initial guess for the quasinormal frequency, which can be achieved by the following approaches. (1) One can make use of the result of another less accurate method, such as P?shl-Theller potential approximation, as input for the FindRoot command. (2) When the black hole metric in question restores to a more straightforward case, for which the quasinormal frequency is already known, by continuously varying its parameters, one may obtain the desired result by gradually changing the parameters of the metric to the relevant values. We note that the above approaches may not work as intended for the case where the quasinormal frequency possesses a bifurcation in the parameter space. However, this problem exists in general for a variety of methods where the eigenvalue problem is expressed in terms of nonlinear equations. (3) One may also make use of the NSolve package to calculate the quasinormal frequencies. In principle, one may find quasinormal frequencies related to different quantum numbers by a single root-finding process. However, the challenge is how to pick out real physical roots in the numerical results. One should proceed with care, and compare the numerical values with those obtained from known cases of black hole metrics. Moreover, the following rule of thumb may provide some general guidance. For some cases, in particular those where the resultant algebraic equation, Eq.(13), is merely an N-th order polynomial, one can exhaust all roots via the NSolve command. Among these roots, some are not physical, and they appear because of a certain symmetry of the equations associated with a particular choice of interpolation points. These unphysical roots can be detected by varying the number of interpolation points, since the roots corresponding to the physical quasinormal frequency will always be present and approach a limit as
Apart from the above two methods, the matrix method can also be adapted to an iterative scheme. This third method is meaningful especially for situations where neither FindRoot nor NSolve is found to be efficient. In general, the iterative relation can be written as follows.
$ \omega_{k+1}={\cal Y}\left(\omega_k,\left|{\cal M}(\omega)\right|\right), $  | (3.21) |
In our code implementation in Refs.[46, 47], we have evenly divided the domain [0,1] ofz. In practice, it might be more favorable to utilize a nonuniform distribution which introduces more interpolation points in the region where the effective potential changes more radically. As a matter of fact, this is an issue that one encounters in other numerical schemes, such as the differential quadrature method and the Runge-Kutta method.
2
3.5.Additional considerations
There are a few more issues or topics which have not been explored in the existing studies but are worth mentioning.Apart for the asymptotically AdS spacetime, our analysis has been mostly carried out in the radial coordinate r. However, sometimes it can be more convenient to express the formalism in the ingoing Eddington–Finkelstein coordinates. In this case, the corresponding boundary condition has also to be modified
$ \Phi\sim \left\{\begin{array}{*{20}{l}} {\rm e}^{{\rm i}\left(\omega+\tilde{\omega}\right) r_*} & r\rightarrow\infty ~\,{\rm{or}}\,~ r\rightarrow r_{\rm c}\\ \Phi_0 & r\rightarrow r_{\rm h} \end{array}\right.. $  | (3.22) |
In practice, instead of carrying out the transformation Eq.(5) and rewriting the boundary condition as Eq.(7), the original boundary condition Eq.(4) might be directly employed in the calculations. As mentioned before, the primary motivation of introducing Eq.(5) is when C0 and C1 are very distinct in their magnitudes. It is understood that such a procedure is useful for the radial coordinate r , as discussed in the present work.
In reality, the black hole is not a static object. The topic of quasinormal modes for dynamic black holes is not only intriguing but also significant in practice. For a dynamical black hole, the event horizon, apparent horizon and infinite redshift surface are distinct concepts, and, therefore, the boundary conditions for quasinormal modes should be tackled with special care.
In certain modified gravity theories, the Lorentz invariance is broken. Owing to the tachyon, the concepts of the event horizon and apparent horizon have to be modified in these theories. Recent studies have proven the existence of the universal horizon [51, 52], which is defined as the boundary that even superluminal particles cannot escape to future infinity. As a part of the effort to find evidence of Lorentz invariance breaking, it is, therefore, interesting to investigate the quasinormal modes for the universal horizon in the framework of the matrix method.
The notion of quasinormal modes is not restricted to black holes. Quasinormal modes of compact stars have since long aroused considerable attention. The primary difference between the two phenomena, from the mathematical point of view, is the existence of the outgoing wave at the boundary of the star surface. This implies a more subtle condition for properly connecting the inside and outside wave functions at the boundary. Nonetheless, the matrix method can also be adapted to handle such problems appropriately.
Last but not least, we note that a quasinormal mode by itself is an eigenvalue of a secular equation. The matrix method is merely a discretized version of the eigenvalue problem in question. One of the alternatives to the eigenvalue problem is the shooting method, which solves the boundary value problem by reducing it to the initial value problem. In this context, one might also implement a specific version of the shooting method for the column vector, which is nothing else but a discretized eigenfunction.
