The matrix method for black hole quasinormal modes

The black hole is one of the most exotic and intriguing subjects in all of theoretical physics. On the experimental side, black hole merger is a magnificent astrophysical phenomenon, it releases an enormous amount of energy in the form of gravitational radiation, larger even than that carried by the light from all stars in the entire visible Universe combined. The gravitational waves traverse the Universe, carrying with them the information about the merger. The first detection of gravitational waves in 2016 was heralded as a remarkable breakthrough in fundamental physics [1, 2]. The observation was announced by LIGO and Virgo collaborations, and it matches the predictions of gravitational waves that emanated from the inward spiral and merger of a pair of black holes [3-5]. Subsequently, further measurements confirmed the gravitational waves from compact-object binaries, including black holes and neutron stars [6-8]. These observations inaugurate a revolutionary era of astronomy, as astrophysicists and cosmologists are now able to make precise observations based on gravitational waves, in addition to electromagnetic radiation. Indeed, while black hole candidates have been identified through electromagnetic signals from nearby matter [9-12], gravitational wave measurements provide a direct test of general relativity, as well as unique access to the properties of spacetime, in the strong field regime. The final fate of a black hole binary consists of three distinct stages: inspiral, merger, and ringdown [13]. In the first stage of the inspiral, the orbit of the black hole binary shrinks gradually due to the emission of gravitational radiation. The inspiral dynamics can be calculated by using the post-Newtonian approximation, which results from a systematic expansion of the full Einstein equations. As the two black holes further evolve and spiral inward, they eventually reach the merger stage, which relates to the strong field, the dynamical regime of general relativity. During the merger stage, the two black holes coalesce into a single, highly distorted remnant black hole, which is surrounded by the combined event horizon. The process is no longer quasi-adiabatic and incredibly sophisticated, where most analytic or semi-analytic methods break down, and numerical relativity is usually employed to calculate the dynamical properties of the system. Gravitational wave emission reaches a peak during the second stage. The remnant black hole eventually settles down into a more dormant state, where distortions from the spherical shape rapidly decrease while emitting gravitational waves. This last stage is known as “ringdown”, in analogy with the way the ringing of a bell is suppressed in time after it is struck. The form of the gravitational wave is characterized by exponentially damped sinusoids. Various techniques in the black hole perturbation theory can be used to tackle the problem. The topic of the present paper are the black hole quasinormal modes, which are closely related to the last stage of the black hole merger process.

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. Revealedbytemporal evolution of small perturbations, the dynamic process of a stable black hole metric also involves three stages. These include the initial outburst, triggered by the source of the perturbation, the quasi-normal oscillations, and late-time asymptotic tail. Usually, the amplitude of the oscillations decays exponentially with time. Otherwise, if small oscillations increase, this indicates that the black hole metric is not stable. The frequency of a quasinormal mode is usually complex, and is known as quasinormal frequency. The real part of the quasinormal frequency defines the frequency of oscillations. The imaginary part, on the other hand, is negative for a stable metric, and it plays the role of suppressing the amplitude of oscillations.

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.

The primary feature of the matrix method [45-47] , which makes it distinct from the others, lies in the fact that the series expansion of the wavefunction is carried out in the vicinity of a series of points, where the latter are not necessarily distributed evenly in the domain of the wavefunction. This feature provides flexibility for achieving a reasonable degree of precision, without sacrificing efficiency. Roughly speaking, in terms of the above grid point expansion, one discretizes the master equation of the perturbation field, and subsequently, rewrites it together with the boundary conditions in the form of a matrix equation. The latter can then be handled readily by various algorithms for solving a system of linear equations.

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:

where $ r_*=\int{{\rm d}r/f(r)} $ is the tortoise coordinate, $ \omega $ and V(r) are the quasinormal frequency and effective potential, respectively. Here, f (r) is determined by the metric, while the black hole event horizon, r_h , and the cosmological horizon, r_c , are often determined by $ f(r_{\rm h})=0 $ and $ f(r_{\rm c})=0 $ , respectively. We proceed by noting that the above form of the master equation can generally be applied to a wide variety of metrics. However, it is not applicable neither to the case of rotating black holes, nor to that given in terms of Eddington–Finkelstein coordinates. The specific form of the master equation in these cases will be discussed with the corresponding topics in the following sections. Usually, the asymptotic behavior of the effective potential at the boundary is relatively simple, and, therefore, in the present work, we will just assume that the effective potential vanishes at the horizon and approaches a constant at infinity $ V(\infty)\equiv V_\infty=\rm{const} $ .

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

where $ \bar{\omega} $ is defined as $ \bar{\omega}\equiv \sqrt{\omega^2-V_\infty} $ , so that $ \bar{\omega}=\omega $ for $ V_\infty=0 $ . If the form of the effective potential is not sufficiently simple, one usually has to resort to numerical methods. Eqs.(1) and (2) furnish an eigenvalue problem, where the boundary conditions are defined at $ r_*\rightarrow\pm\infty $ in terms of the tortoise coordinate. We note that the above discussion is valid for asymptotically flat and dS spacetimes, while the boundary conditions for asymptotically AdS spacetime, and the angular part of the master equation, will be addressed below.

Similarly to the other methods based on series expansion, we first change the domain of the radial variable to a finite range by introducing $ z=z(r) $ . For convenience, we choose the new domain to be $ z\in[0,1] $ , and the boundaries correspond to the points z=0 and z=1, respectively. However, since the wave function is also transformed, the boundary conditions in Eq.(2) will also be affected. In this regard, we rewrite the wave function in the form $ \Phi=A(z)R(z) $ , where $ A(z(r)) $ carries the asymptotic behavior of the wave function at the boundary given in Eq.(2). As a result, R(z) is expected to be well behaved and, in particular, to be regular at the boundary. To be specific, the master equation reads

Here, $ {\cal H}(\omega,r) $ is a linear operator, determined by Eq.(1), which depends on $ \omega $ and z. C₀ and C₁ are constants, since the asymptotic part of the wave function has been factorized. It is convenient to further introduce

where $ z(1-z) $ can be replaced by other smooth functions which are zero only at z=0 and z=1. We note that, for the case of AdS spacetime, since $ R(r\rightarrow\infty)\rightarrow 0 $ , the transformation Eq.(5) can be replaced by $ F(z)=R(z)(1-z) $ where $ z=r_{\rm h}/r $ .

Subsequently, the master equation can be rewritten with F(z)

with the boundary conditions in an even simpler form

Here, the linear operator $ {\cal G}(\omega,z) $ is derived from the form of $ {\cal H}(\omega,z) $ and Eq.(5). We note that, as a numerical trick, Eq.(5) is not an absolutely necessary procedure. In practice, the resultant wave function no longer possesses a boundary condition involving an undetermined constant. Furthermore, when C₀ and C₁ turn out to be quite distinct in their values, the factor $ {z(1-z)} $ helps to suppress the function values at the boundary points, as well as in their vicinity. Intuitively, the only scenario when the procedure does not work as intended is when the form of the wave function becomes more fluctuating due to the factor $ {z(1-z)} $ for some particular reason, and the subsequent Taylor expansion becomes less efficient. But, as we do not have a good estimate of the form of the wave function at this point, the introduction of Eq.(5) helps to maintain the resultant equation simple, at the least. The existing numerical results are in most cases found to be consistent with the above speculations.

As the next step, we descretize the master equation Eq.(6) by introducing N interpolation points z_i with $ i=1, 2, \cdots, N $ , in the interval $ z\in[0,1] $ . Subsequently, the wave function F(z) is also represented by descrete values $ f_i=F(z_i) $ on the grid. As a result, one discretizes the differential equation Eq.(6) and rewrites it in terms of f_i. In the present algorithm, various derivatives of the wave function at z_i are obtained by inverting a $ N\times N $ matrix [45], so that the information about the function values in the entire interval is utilized. Besides, owing to the fact that the above inverting process is formal, it is carried out in practice before the real numerical calculations, and, therefore, it does not have any negative impact on the efficiency of the method. The resulting matrix equation can be formally written as

where $ \bar{\cal M} $ represents the descretized version of the operator $ {\cal G}(\omega,z) $ , and is a $ N\times N $ matrix which is still a function of the quasinormal frequency $ \omega $ , and

is a column vector composed of f_i. The boundary condition Eq.(7) implies that

We make use of the above condition to replace the first and last line of the matrix $ \bar{\cal M} $ . The reason behind this choice is to select a line or column in the original equation which makes use of the least amount of information about of grid points. Subsequently, Eq.(8) becomes

where the element $ {\cal M}_{k,i} $ of the matrix $ {\cal M} $ is defined by

The matrix equation indicates that $ {\cal F} $ is the eigenvector of $ {\cal M} $ , which implies

Eq.(13) is a non-linear algebraic equation for the quasinormal frequency $ \omega $ , which can be solved by any standard nonlinear equation solver.

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 $ {\cal M}(\omega) $ will be “almost” diagonal, since only the elements $ {\cal M}_{k,i} $ with $ i=k, k\pm 1 $ are non-vanishing. Subsequently, the resulting wave function can be obtained via an iterative method, similar to the continued fraction method, or the HH method. However, as mentioned above, a sparse matrix indicates that the corresponding discretization scheme does not fully make use of the information about the entire interval $ z\in[0,1] $ . On the other hand, in several different methods, such as the continued fraction method and the HH method, the series expansion of the wave function is only carried out in a single radial position, for instance, the horizon. From a different point of view, the resulting algebraic equation can in fact also be rephrased in the form of a matrix equation, Eq.(11). Owing to the fact that the iterative relation only involves a few expansion coefficients, the corresponding matrix is also rather sparse. Therefore, for a given N, the matrix method is related to a Taylor expansion of the order N, which subsequently provides better precision. Moreover, one may even seek out an optimized grid distribution, which provides a higher resolution for the region where the variation of the wavefunction is more dramatic, by inserting more grid points. In this regard, the present algorithm can be viewed as an improvement over many existing approaches.

The following section is devoted to a more detailed discussion of various applications.

While the general strategy for the evaluation of the black hole quasinormal modes presented in the previous section is concise and straightforward, many particular aspects, such as those related to the boundary conditions and specific characteristics of the master equations, might bring up complications for particular problems. In the present section, we discuss in detail several specific cases for the applications of the matrix method. In particular, our discussion addresses the boundary conditions for the black hole metrics in asymptotically flat, dS, and AdS spacetimes, as well as the associated ansatz for the method of separation of variables. The background metrics concern static, rotational, and extreme black holes. The algorithm presented in the following does not require the existence of an analytic form of the tortoise coordinate. We also comment on the case of a numerical black hole metric, and that where the master equation is a system of coupled ordinary differential equations.

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 r_c between the observer and infinity, which is also a Cauchy one-way membrane. An observer staying between r_h and r_c will not receive any signals from the inside of event horizon, nor from the outside of cosmological horizon. Therefore, the physically valid domain is $ r_{\rm h}\leqslant r \leqslant r_{\rm c} $ for the dS spacetime, while it is $ r_{\rm h} \leqslant r \leqslant \infty $ for the asymptotically flat and AdS spacetimes.

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 $ f'(r_{\rm h})=0 $ , which implies a Taylor expansion of f (r) around the horizon that does not contain the linear term. Subsequently, as shown below, the leading contribution in the tortoise coordinate becomes different.

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 $ \omega $ and $ \lambda $ , corresponding to the radial and angular parts of the wave functions. The latter restores to a simple case when the parameter of the black hole related to angular asymmetry vanishes. In particular, if the black hole metric reduces to a static spherical one, one finds $ \lambda=\ell(\ell+1) $ , namely, the angular momentum quantum number. The angular part of the master equation becomes decoupled, and the angular wavefunction reduces to spherical harmonics. On the other hand, if the angular part of the master equation is coupled to the radial part, its boundary condition is the usual periodic boundary condition.

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 $ r_*(r) $ , and as a result it is not straightforward to obtain the boundary condition in an analytic form. The second difficulty is apparently even more severe, and it concerns the numerical black hole metric. Fortunately, as we show below, the matrix method is quite versatile in handling these subtle situations. In fact, it is relatively intuitive to understand the reason. The critical procedure of the matrix method is the discretization of continuous functions into discrete values on a grid. Naturally, as a result, the implementation of the method does not rely on analytic expressions of the tortoise coordinate, as shown in the following subsection. Moreover, the method neither relies on the analytic form of the black hole metric itself. It is rather straightforward to slightly alter the procedure to adopt it to the numerical black hole metric. This is not the case with the other approaches, such as the HH method, where the numerical difficulty is enhanced by the requirement of precision values for the high order derivatives of the metric. Last but not least, for the difficulties with the numerical tortoise coordinate, one does not require an explicit form of r_* but its asymptotical behavior, and, therefore, one can merely expand the form of r_* near the boundary and extract the desired 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 $ f(r)=\hat{f}_1(r-r_{\rm h})+\displaystyle\frac{\hat{f}_2}{2}(r-r_{\rm h})^2+{\cal O}(r-r_{\rm h}) $ (where $ \hat{f}_i\equiv f^{(i)}(r_{\rm h}) $ and $ {\cal O} $ represents terms of order higher than $ (r-r_{\rm h})^2 $ ), we have

which implies that the boundary condition at the horizon is given by

Case 2. For the extreme black hole, near the horizon, if $ f(r)\!=\!\displaystyle\frac{\hat{f}_2}{2}(r\!-\!r_{\rm h})^2\!+\!\displaystyle\frac{\hat{f}_3}{6}(r\!-\!r_{\rm h})^3\!+\!\displaystyle\frac{\hat{f}_4}{24}(r\!-\!r_{\rm h})^4\!+\!{\cal O}(r\!-\!r_{\rm h}) $ , where $ {\cal O} $ are terms of order higher than $ (r-r_{\rm h})^4 $ , we find

and the boundary condition at horizon is given by

Case 3. In the asymptotically flat spacetime, for $ r\to \infty $ , if $ f(r)=1+A r^\alpha+{\cal O} $ , where $ \alpha<0 $ and $ {\cal O} $ represents higher order terms less significant than $ r^{\alpha} $ , we have

Thus, the boundary condition at infinity is given by

and the second term could be significant when $ \alpha>-1 $ .

Case 4. In the asymptotically de Sitter spacetime, according to Case 1 and Case 2, one obtains similar results by replacing r_h with r_c. Therefore, following Eq.(15), near the non-extreme cosmological horizon, the boundary condition is

On the other hand, for the extreme cosmological horizon, we have

where $ \bar{f}_i\equiv f^{(i)}(r_{\rm c}) $ .

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

where the term $ r^{\frac{i\omega}{\hat{f}_1}} $ is introduced to cancel out the effect of the factor $ \left(r-r_{\rm h}\right)^{-\frac{i\omega}{\hat{f}_1}} $ at infinity, while it gives a rather smooth contribution at the horizon. Thus, it can be readily verified that the resultant expression, Eq.(22), indeed possesses the desired asymptotic forms at $ r\to r_{\rm h} $ and $ r\to\infty $ found in Eqs.(15) and (18). Likewise, for the extreme black hole, we have

As discussed above, for the black hole in asymptotically de Sitter spacetime, it is required that at the event horizon $ \hat{f}_1\not=0 $ and $ \hat{f}_1=0 $ for the non-extreme and extreme cases, respectively. Similarly, at the cosmological horizon, we have $ \bar{f}_1\not=0 $ and $ \bar{f}_1=0 $ , respectively for the non-extreme and extreme scenarios. Therefore, the boundary conditions in de Sitter spacetime are satisfied by assuming the following form for the wave functions

For the asymptotically anti-de Sitter spacetime, the Horowitz-Hubeny method usually employs the ingoing Eddington-Finkelstein coordinates $ (v,r) $ , with $ v=r_*+t $ . In what follows, we will also use this convention. The corresponding master equation reads

where $ U(r)=V(r)/f(r) $ , and since the potential is divergent at infinity, the boundary condition requires $ \Phi $ to be constant at the event horizon, while it must vanish at infinity.

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:

The boundary conditions are at $ u=\pm1 $ , where one defines $ W^2_\pm ={\cal W}(u=\pm 1,\omega,\lambda) $ . The asymptotic solutions at the boundary are given by

As it is required that the angular part of the wave function Y(u) always remains finite, it can be assumed to have the following form

In the above expressions for the wave function with separation of variables, Eqs.(22)-(28), we have implied that the functions R(r) and S(u) are finite and well behaved in the entire domain of r and u. That is to say, $ r_{\rm h}\leqslant r < \infty $ for the asymptotically flat and anti-de Sitter spacetimes, $ r_{\rm h}\leqslant r\leqslant r_{\rm c} $ for the asymptotically de Sitter spacetime, and $ -1\leqslant u\leqslant 1 $ for the angular part of the master equation. This is because the asymptotic form of the wave function at the boundary has already been factorized and, therefore, effectively extracted from the l.h.s. of Eqs.(22)-(28).

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 $ x=1-\displaystyle\frac{r_{\rm h}}{r} $ or $ z=r_{\rm h}/r $ , where r_h is the event horizon, so that at the event horizon x=0 and z=1, while x=1 and z=0 represent infinity. For the asymptotically de Sitter black hole spacetime, due to the existence of the cosmological horizon r_c, we choose $ y=\displaystyle\frac{r-r_{\rm h}}{r_{\rm c}-r_{\rm h}} $ , so that y=0 represents the event horizon while y=1 corresponds to the cosmological horizon. For the angular coordinate, we introduce $ \nu=\displaystyle\frac{1+u}{2} $ , so that one has $ \nu=0 $ and $ \nu=1 $ for u=−1 and u=1, respectively.

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.

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 $ nN\times nN $ , where N is the number of interpolation points introduced above. As an example, in the case of n=2, we have

where $ {\cal B} $ is also a $ N\times N $ matrix. In this case, one may rewrite the above equations in terms of a double-sized matrix equation, namely,

which can be readily solved as before. The only disadvantage is that the matrix on the l.h.s. of Eq.(30) might be a sparse matrix, which subsequently implies a loss of computational efficiency.

Alternatively, one can either rewrite Eq.(30) as

or as

Here, the expressions involve the evaluation of $ {\cal B}_a^{-1} $ or $ {\cal B}_b^{-1} $ , when the relevant inverse exists. To solve Eq.(31) or Eq.(32) , it is necessary to find the roots of a nonlinear equation in $ \omega $ , associated with the $ N\times N $ matrices. The second approach may potentially improve the efficiency.

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:

where “precs” is the desired precision of the function “TU”. If one assigns precs=100, the function output will be given in terms of a fraction truncated to a hundred significant digits.

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 $ N\to \infty $ . Nonetheless, to properly remove all unphysical roots is not a simple task. The situation becomes even more complicated when Eq.(13) is highly nonlinear. In this regard, the above mentioned Macro can be employed to alleviate the situation.

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.

where the quasinormal frequency of the (k+1)-th interation, $ \omega_{k+1} $ , is determined by the k-th result. The operator $ {\cal Y} $ is a functional of the determinant satisfying $ \omega={\cal Y}\left(\omega,\left|{\cal M}\right|\right) $ , where $ \omega $ is the root of Eq.(13). By choosing an appropriate form of $ {\cal Y} $ , one obtains the result for quasinormal frequency by repeatedly using Eq.(25). The primary challenge of this approach is to find a proper iteration relation for $ {\cal Y}(x,y) $ which efficiently leads to a convergent result.

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.

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

where $ \tilde{\omega}=\sqrt{\omega^2-V_\infty} $ for $ r\rightarrow\infty $ in the asymptotically flat spacetime, and $ \tilde{\omega}=\omega $ for $ r\rightarrow r_{\rm c} $ in the asymptotically dS spacetime. We note that, in a similar fashion, the master equation can also be derived in the outgoing Eddington-Finkelstein coordinates.

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 C₀ and C₁ 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.

In this paper, a comprehensive survey of the applications of the matrix method for black hole quasinormal modes was presented. It was shown that the proposed algorithm is independent of any particular analytical form of the tortoise coordinate. In fact, the method can even be employed to a numerical black hole metric. Moreover, it is argued that the present approach can readily be applied to the case where the master equation involves a system of coupled equations. Our discussion covered various types of black hole metrics, master equations, and the corresponding boundary conditions. The proposed method is reasonably accurate, efficient, as well as flexible for distinct physical scenarios. We look forward to carrying out further studies of several open questions addressed in this survey.