Next Article in Journal
Nontrivial Isometric Embeddings for Flat Spaces
Previous Article in Journal
Near-AdS2 Spectroscopy: Classifying the Spectrum of Operators and Interactions in N=2 4D Supergravity
Previous Article in Special Issue
Ringing of the Regular Black Hole with Asymptotically Minkowski Core
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

Spectral Problems for Quasinormal Modes of Black Holes

Department of Physics, Rikkyo University, Toshima, Tokyo 171-8501, Japan
*
Author to whom correspondence should be addressed.
Submission received: 31 October 2021 / Revised: 26 November 2021 / Accepted: 28 November 2021 / Published: 4 December 2021

Abstract

:
This is an unconventional review article on spectral problems in black hole perturbation theory. Our purpose is to explain how to apply various known techniques in quantum mechanics to such spectral problems. The article includes analytical/numerical treatments, semiclassical perturbation theory, the (uniform) WKB method and useful mathematical tools: Borel summations, Padé approximants, and so forth. The article is not comprehensive, but rather looks into a few examples from various points of view. The techniques in this article are widely applicable to many other examples.

1. Introduction

The main motivation to write up this article is the following. We collect various traditional approaches to one-dimensional eigenvalue problems in quantum mechanics. Some of them are well explained in the usual textbooks, but some are not. We would like to demonstrate that these approaches are widely applicable to spectral problems in black hole perturbation theory. We also present a unified manner to treat bound states and resonant states together. The latter problem is particularly important in the analysis of quasinormal modes (QNMs) of black holes.
We guess that some of the readers are familiar with general relativity but may not be so familiar with quantum mechanics. Conversely, some may be familiar with quantum mechanics but not with general relativity. The article is pedagogical and designed for both of these people.
There have already been many excellent review articles [1,2,3,4,5,6] on black hole perturbation theory. For differentiation from these articles, the contents of the current article are intentionally quite biased. Although we review the linear perturbation of the four-dimensional Schwarzschild spacetime in Appendix A, the reader who wants to learn black hole perturbation theory in more general cases should see other good reviews [1,2,3,4,5,6] or textbooks [7,8,9,10] and references therein. In the main text, we rather concentrate our attention on practical aspects of computational schemes on eigenvalue problems associated with black hole perturbation theory.
Let us briefly comment on physical significances of QNMs. They are defined as poles of Green’s function in initial value problems of the linear perturbation around black hole solutions, and describe exponentially damped oscillation as the late time behavior of the perturbation1 [11,12,13,14,15]. Interestingly, they also describe the last stage of a gravitational waveform in a black hole merger, which was firstly observed by LIGO–Virgo [16], while the black hole merger process should be treated as a fully non-linear problem. A recent argument [17] shows that the gravitational waveform of the binary black hole merger can be fit by the superposition of QNMs even before the peak amplitude. This suggests that the linear perturbation treatment around the final state black hole is a good approximation even just after the black hole merger. In general relativity, because of the uniqueness of the Kerr black holes, the QNM spectra are completely characterized by their mass and spin. If we can check it by observation, it will be a good test for general relativity. In other words, we have a possibility to test theories beyond general relativity through the observation of QNMs [18]. For this purpose, the QNM spectra are particularly important. The computational technique introduced in the present article can be also applied for such problems.
Mathematically, the spectral problems in this article are connection problems of local solutions to ordinary differential equations at two different spatial points, that is, two-point boundary value problems. Since these problems require global information on the solutions, they are not solved analytically in general. This fact makes spectral problems rich and interesting. We would like to investigate the spectral problems as analytic as possible.
The organization of the article is as follows. In the next section, we review several classic approaches to spectral problems in quantum mechanics. Some of them are less familiar. It includes analytical/numerical treatments, perturbation theory and the WKB method. All of them have direct applications to black hole physics. In most cases, spectral problems are not exactly solvable. It is desirable to present widely applicable ways to various problems. To understand this section, the reader needs basics on second order ordinary differential equations, asymptotic analysis and Padé approximants, which are summarized in Appendices Appendix B and Appendix C. In Section 3, we discuss how to apply these techniques to spectral problems in black hole physics. We consider two spectral problems: quasinormal modes and mode (in)stability of black holes by seeing concrete examples. More basics on black hole perturbation theory are explained in Appendix A.

2. A Quantum Mechanical Approach to Spectral Problems

In this section, we review various ways for spectral problems in quantum mechanics. In the next section, we will discuss actual applications to black hole physics.

2.1. Bound States and Resonant States

We begin by reviewing bound states and resonant states in quantum mechanics. Let us consider the one-dimensional time-independent Schrödinger equation of the form:
2 d 2 d x 2 + V ( x ) ψ ( x ) = E ψ ( x ) .
This type of the eigen-equations appears in many places in theoretical physics. Different contexts often provide us different perspectives to the same problem. In many cases, just appears as a parameter of the equation. A typical problem in quantum mechanics is to compute energy eigenvalues of bound states for a given potential. Another interesting problem is a potential scattering. As we see just below, this problem is related to resonant states. These two problems turn out to be interrelated to each other. The former problem is related to mode stability of black holes, and the latter has a direct connection with quasinormal modes.
As is well-known, the bound states require the normalizability of the wave function. Therefore the potential needs to have a global minimum in a considered region. However, this does not mean that a potential with a minimum always has a bound state. Unless a well is deep enough, no bound states exist. This point is crucially important in the analysis of (in)stability problems of black holes.
The resonant states are less familiar. Let us consider a scattering problem for a potential with a wall, shown in Figure 1.2
The Schrödinger equation in this setup is:
˜ 2 d 2 d x 2 + V ˜ ( x ) ψ ˜ ( x ) = E ˜ ψ ˜ ( x ) .
There is an incoming wave from x = + , and it scatters with the potential wall. As a result, there are a reflected and a transmitted waves. In the region x + , the wave function is written as:
ψ ˜ ( x ) ψ ˜ in ( x ) + R ψ ˜ r ( x ) , x + ,
where R is the reflection coefficient. In the region x , only the transmitted wave exists:
ψ ˜ ( x ) T ψ ˜ t ( x ) , x ,
where T is the transmission coefficient. It is more convenient to change the overall normalization as follows:
B ψ ˜ t ( x ) x ψ ˜ ( x ) x + A ψ ˜ in ( x ) + ψ ˜ r ( x ) ,
where A = 1 / R and B = T / R . The resonant states are defined by the absence of the incoming wave:
A = 1 R = 0 .
Clearly, it happens for characteristic energies that are located on singularities of the reflection coefficient (or equivalently the transmission coefficient). Hence, this phenomenon is called a resonance. In general, these characteristic eigenvalues are discrete but complex-valued [19]. Note that the boundary condition for the resonant states is the same as that for the quasinormal modes of black holes. Computing the resonant energies is the main issue in this article.
Next, we see a relation between bound states and resonant states. In this article, we particularly focus on two types of potentials V ( x ) and their sign-flipped ones V ˜ ( x ) = V ( x ) , as shown in Figure 2. Since the potential V ( x ) has a well in these cases, it admits bound states. Then its inverted potential V ˜ ( x ) has resonant states. The energy eigenvalues of these two states are simply related by an analytic continuation of the quantum parameter . If setting = i ˜ ,3 the eigen-Equation (1) is related to (2) with E ˜ = E and ψ ˜ ( x ) = ψ ( x ) | ˜ .4 Therefore physics in the system (2) may be described by that in the system (1). In particular, we expect the following exact spectral relation:5
E ˜ n resonance ( ˜ ) = E n bound ( = i ˜ ) , n = 0 , 1 , 2 , ,
where E n bound ( ) is the bound state energy for the Schrödinger Equation (1), and E ˜ n resonance ( ˜ ) is the resonant energy for its sign-flipped one (2). We will test this equality in a few examples in detail. We can extract the resonant energy from the bound state energy. This idea has been applied to the computation of the QNM frequencies in [20,21] (see [22,23,24] for earlier works).

2.2. Lessons from an Exactly Solvable Model

In the context of black hole perturbation theory, the Pöschl–Teller potential V ˜ PT ( x ) = V 0 / ( 2 cosh 2 β x ) is usually taken as an exactly solvable example [22]. This is reasonable because its shape is similar to effective potentials of master equations appearing in black hole perturbations, as in the right panel of Figure 2. For differentiation, we take another less-familiar example: the Morse potential. The Morse potential is also exactly solvable. Its asymptotic behavior is quite different from the potentials in black hole perturbations. See the left panel of Figure 2. However, from the viewpoint of the singularity structure of ordinary differential equations (see Appendix B), the Morse potential is more similar to the spectral problem for asymptotically flat black holes. The Pöschl–Teller potential is rather similar to the spectral problem for asymptotically (anti-)de Sitter black holes. Since the Pöschl–Teller potential has been reviewed in detail in [5], we do not discuss it here.
Let us consider the following very special potential:
V Morse ( x ) = V 0 ( e 2 β x 2 e β x ) , < x < ,
where we assume V 0 > 0 and β > 0 . Its shape is shown in the left panel in Figure 2. It has the global minimum at x = 0 . A simple reason of the solvability of the Morse potential is a shape invariance [25]. Another explanation is that the Schrödinger equation is mapped to the well-known equation as we will see below. To make expressions simpler, we introduce the following rescaled parameters:
q : = β x , ε : = E V 0 , g : = β V 0 .
Then the Schrödinger equation is rewritten as:
g 2 d 2 d q 2 + e 2 q 2 e q ψ ( q ) = ε ψ ( q ) .
We change the variable z = ( 2 / g ) e q . The differential Equation (10) becomes:
z d d z z d d z + z 4 z 4 g ψ ( z ) = ε g 2 ψ ( z ) .
We further perform the transform:
ψ ( z ) = e z / 2 z κ g y ( z ) ,
where κ : = ε . Since we are interested in the bound states, we implicitly assume 1 < ε < 0 and κ > 0 . However, in the construction of the solutions, we do not need this assumption. Then the differential equation reduces to the generalized (or associated) Laguerre equation:
z y ( z ) + ( α + 1 z ) y ( z ) + ν y ( z ) = 0 ,
where
α = 2 κ g , ν = 1 κ g 1 2 .
The generalized Laguerre equation has the regular singular point at z = 0 and the irregular singular point at z = . It is essentially equivalent to the confluent hypergeometric equation. The (characteristic) exponents at z = 0 are 0 and α . The former solution corresponds to the generalized Laguerre function L ν α ( z ) . Some of basics on ordinary differential equations are explained in Appendix B.
Let us see the eigenvalues of the bound states. In the limit z 0 (i.e., q ), since the wave function behaves as ψ ( z ) z κ g or ψ ( z ) z κ g , the normalizability requires the regularity of y ( z ) at z = 0 . This means that we have to choose y ( z ) = L ν α ( z ) . Therefore the analytic solution satisfying the boundary condition at z = 0 is given by
ψ ( z ) = e z / 2 z κ g L ν α ( z ) .
This solution of course does not satisfy the boundary condition at z = for arbitrary E because L ν α ( z ) exponentially grows in z . The boundary condition at infinity is satisfied if and only if L ν α ( z ) is a polynomial. This requires ν to be a non-negative integer. Therefore, we obtain an exact quantization condition for the bound states:
1 ε n g 1 2 = n , n = 0 , 1 , 2 ,
This is easily solved, and we finally obtain:
ε n = 1 g n + 1 2 2 , n = 0 , 1 , 2 .
The same result is obtained by the asymptotic expansion of L ν α ( z ) . It is well-known that the generalized Laguerre function has the following asymptotic expansion:
L ν α ( z ) ( z ) ν Γ ( ν + 1 ) k = 0 ( 1 ) k ( ν ) k ( α ν ) k k ! z k sin ( π ν ) Γ ( α + ν + 1 ) z α ν 1 e z π k = 0 ( ν + 1 ) k ( α + ν + 1 ) k k ! z k ( z ) ,
where ( a ) k is the Pochhammer symbol. We have used ≃ instead of = because the right hand side contains formal divergent series, that is, the radius of convergence is just zero. We have to treat it in a careful manner. We return to this issue in the next subsection and in Appendix B.3. Since the second line in this asymptotic expansion is a source of the exponential grow, the boundary condition at infinity requires the absence of it: sin ( π ν ) = 0 . Then, the right hand side in the first line has finite terms, and it reduces to a polynomial.
Let us count the number of the bound states. The positivity κ n = ε n > 0 for the bound states leads to the upper bound to the quantum number n:
n < 1 g 1 2 .
Therefore, the number of the allowed bound states, N bound , is given by:
N bound = 1 g 1 2 ,
where p means the least integer greater than or equal to p. For instance, if g 2 , there are no bound states although the potential still has the well. We show N bound against the parameter g in the left panel of Figure 3.
Let us recall the exact spectrum (17). At first glance, it seems to give an infinite number of the eigenvalues for n = 0 , 1 , 2 , . This is not the case. As seen above, for relatively large n, κ n becomes imaginary. How do we exclude such n’s by looking at only the spectrum (17)? In the representation (17), this is understood as two branches of the square root of ε n . We show in in the right panel of Figure 3 the spectrum ε n against the parameter g. For a given negative energy ranging 1 < ε < 0 , we have two ε n with quantum number n corresponding to different g. The left branch (solid line) gives the physical bound state energy because it is continuously connected with the classical point g = 0 . The right branch (dashed line) is an unphysical mode that does not satisfy the bound state boundary condition.6 For a fixed value of g, highly excited states are mostly on the unphysical branch.
Things are different for the resonant states. As mentioned in the introductory section, the resonant states are related to the bound states by the analytic continuation: g = i g ˜ (or g = i g ˜ ). In this continuation, the parameters become:
z = 2 g e q = 2 i g ˜ e q , α = 2 κ i g ˜ , ν = 1 κ i g ˜ 1 2 .
The eigenvalues and the eigenfunctions for the inverted potential V ˜ are given by:
ε ˜ n : = ε n = 1 i g ˜ n + 1 2 2 , ψ ˜ n ( q ) = e i e q / g ˜ 2 i g ˜ e q i κ n / g ˜ L n α 2 i g ˜ e q .
One can check that this eigenfunction always satisfies the boundary condition of the resonant states for all n = 0 , 1 , 2 , . We conclude that there are always an infinite number of the resonant states7 while the number of the bound states is finite. To construct highly excited resonant states, we need unphysical “false bound states”.

2.3. Numerical Methods

There are many ways to evaluate the eigen-energies numerically. Since our purpose is not to cover all of these methods, we present only two particular methods that are useful in our later analysis.

2.3.1. Milne’s Method

Here we revisit a very old result by Milne [26] because it is useful to numerically count the number of bound states. We first rewrite the Schrödinger Equation (1) as:
2 d 2 d x 2 + Q ( x ) ψ ( x ) = 0 ,
where
Q ( x ) = E V ( x ) .
We assume that Q ( x ) does not have any singular points for x R . We consider two independent solutions satisfying conditions:
ψ 1 ( x 0 ) = 1 , ψ 1 ( x 0 ) = 0 , ψ 2 ( x 0 ) = 0 , ψ 2 ( x 0 ) = 1 ,
where x 0 is a certain point on the real axis. It is well-known that the Wronskian of these two solutions does not depend on x, and it always equals to unity. Let us define
w E ( x ) : = ψ 1 ( x ) 2 + ψ 2 ( x ) 2 ,
where we use the subscript E to emphasize the energy dependence of w. The fact that the Wronskian of ψ 1 ( x ) and ψ 2 ( x ) is the unity guarantees that w E ( x ) never vanishes. It satisfies
w E ( x 0 ) = 1 , w E ( x 0 ) = 0 .
One can easily check that w E ( x ) satisfies the following non-linear differential equation:
2 d 2 d x 2 + Q ( x ) w E ( x ) = 2 w E ( x ) 3 .
The two solutions are conversely reconstructed by:
ψ 1 ( x ) = w E ( x ) cos x 0 x d x w E ( x ) 2 , ψ 2 ( x ) = w E ( x ) sin x 0 x d x w E ( x ) 2 .
Now, we would like to construct a solution that satisfies the decaying boundary condition at x = . Let us consider the following function:
ψ ( x ) = w E ( x ) sin x d x w E ( x ) 2 .
This is a solution to (23) because the function ψ ( x ) is written as a superposition of the two solutions ψ 1 ( x ) and ψ 2 ( x ) .
ψ ( x ) = ψ 1 ( x ) sin x 0 d x w E ( x ) 2 + ψ 2 ( x ) cos x 0 d x w E ( x ) 2 .
This function also behaves as ψ ( x ) C / ( x w E ( x ) ) in x , and actually satisfies the boundary condition at x = . Therefore, this is what we want. Then, the decaying boundary condition at x = requires:
sin d x w E ( x ) 2 = 0 .
Let us define
N Milne ( E ) : = 1 π d x w E ( x ) 2 1 .
The important fact found in [26] is that this function is a non-decreasing function of E.8 It turns out that for the bound state energy E n , the function satisfies9
N Milne ( E n ) = n , n = 0 , 1 , 2 , .
This is a kind of quantization condition. Therefore, N Milne ( E ) can be regarded as a counting function of the bound states. If m 1 N Milne ( E ) < m , then there are m bound states below the energy E.
Now we apply this simple method to the Morse potential (8). In the actual computation, we use (10). To fit the notation in this subsection, we equivalently fix the parameters V 0 = β = 1 . There is no difficulty to solve (28) with (27) numerically for given E. This is just an initial value problem of the second order ordinary differential equation. We can easily evaluate N Milne ( E ) . Solving N Milne ( E ) = n , we obtain the bound state energy at level n. In the computation, we set x 0 = 0 . It is interesting to notice that we have another exact quantization condition (16). Since the left hand side in (16) is obviously an increasing function,
N exact ( E ) : = 1 E 1 2
is also a counting function.
We show in Figure 4 the behaviors of N Milne ( E ) and N exact ( E ) for = 1 / 5 . Interestingly, these two functions are quite different, but N Milne ( E ) indeed gives the correct bound state spectra. The total number of the bound states is evaluated by the values of the counting functions at E = 0 because V Morse 0 in x . In the case of = 1 / 5 , since N Milne ( 0 ) 4.1 or N exact ( 0 ) = 4.5 , we have the five bound states. Of course it is consistent with the exact result (20). In other words, we always have
N bound = N exact ( 0 ) = N Milne ( 0 ) .
We can use this relation in the mode stability analysis of black holes. On the other hand, it seems hard to apply this method to resonant state problems. To compute their eigenvalues numerically, one needs other approaches.

2.3.2. Wronskian Method

We also review the well-known Wronskian method. The Wronskian is useful to evaluate connection coefficients among local solutions. Let y p 1 ( x ) and y p 2 ( x ) be two independent solutions around x = p . Let us consider a two-point boundary value problem between x = a and x = b . The four solutions should be related by the connection coefficients:
y a 1 ( x ) = C 11 y b 1 ( x ) + C 12 y b 2 ( x ) , y a 2 ( x ) = C 21 y b 1 ( x ) + C 22 y b 2 ( x ) .
The connection coefficients C i j are evaluated by Wronskians. Let us assume that y a 1 ( x ) and y b 1 ( x ) satisfy boundary conditions we are interested in. Then these two boundary conditions require:
C 12 = W [ y a 1 , y b 1 ] W [ y b 2 , y b 1 ] = 0 W [ y a 1 , y b 1 ] = 0 .
This determines the eigenvalues of the boundary value problem. In other words, C 12 = 0 means that the two solutions y a 1 ( x ) and y b 1 ( x ) are linearly dependent, and their Wronskian must vanish. This method works not only for bound state problems but also for resonant state problems.
Let us apply it to the Morse potential. It is more convenient to use the z-variable. We consider the boundary value problem for the differential Equation (13). As already seen, the regular solution at z = 0 is given by y 01 ( z ) = L ν α ( z ) . The construction of the local solutions at infinity is more involved because z = is the irregular singular point. See Appendix B.2 on formal series solutions at an irregular singular point. Here we take a shortcut. We use the asymptotic expansion (18). This tells us the all-order expressions of the two formal series solutions at z = :
y 1 formal ( z ) = z ν k = 0 ( 1 ) k ( ν ) k ( α ν ) k k ! z k , y 2 formal ( z ) = z α ν 1 e z k = 0 ( ν + 1 ) k ( α + ν + 1 ) k k ! z k .
One can check that these are actually solutions to (13). The former solution satisfies the boundary condition at infinity. An important caution is that the sums on the right hand side are divergent series. They do not converge for any values of z and generic values of α and ν . This is why we call these solutions formal. To see it in detail, let us consider the truncated sum in the first equation in (39):
S m ( z ) : = k = 0 m ( 1 ) k ( ν ) k ( α ν ) k k ! z k .
We plot some values of S m ( z ) for ν = 1 / 3 , α = 3 / 2 and z = 1 , 2 , 3 in Figure 5. Obviously S m diverges when m gets large. The dashed lines are the values of the Borel summation (see Appendix B.2) of the infinite sum S ( z ) . By setting p = 1 in (A142), the Bore sum of S ( z ) is analytically given by:
S Borel ( z ) = 0 d ζ e ζ 2 F 1 ( ν , α ν ; 1 ; ζ / z ) ,
where 2 F 1 ( a , b ; c ; z ) is the Gauss hypergeometric function. This integral is convergent for any z ( 0 , ) , and reproduces the asymptotic series S ( z ) . Up to a certain order, S m ( z ) gets close to the Borel sum S Borel ( z ) . This maximal order is called an optimal order. The optimal order m opt ( z ) depends on z. Its simple estimation method is explained in [30]. In the present case, m opt ( 1 ) = 3 , m opt ( 2 ) = 4 and m opt ( 3 ) = 6 . See the right panel in Figure 5. Beyond the optimal order, the approximation of S Borel ( z ) by the finite sum S m ( z ) gets worse. This fact sometimes causes a confusion in numerics. If a series expansion is a divergent series, one has to choose a truncation order very carefully. A facile calculation makes an approximation worse.10 The Borel summation avoids this problem.
Now we apply the Borel summation to the formal power series (39). The Borel summation of the first equation in (39) is given by:
y 1 Borel ( z ) = z ν S Borel ( z ) = z ν 0 d ζ e ζ 2 F 1 ( ν , α ν ; 1 ; ζ / z ) .
By construction, it has the correct asymptotic expansion in 1 / z . This solution is analytic in a Stokes sector π < arg z < π , but the Stokes phenomenon happens along arg z = ± π . See Appendix B.3 for the Stokes phenomenon. In our analysis here, the Stokes phenomenon is not important because we are interested in arg z = 0 . We look for the zeros of the Wronskian for the two analytic solutions y 01 ( z ) and y 1 Borel ( z ) . One can check that the Wronskian for z > 0 vanishes when ν is a non-negative integer.
Of course, in most examples, we cannot construct analytic solutions at boundary points. We have to use truncated series solutions to evaluate their Wronskians. If a boundary point is an ordinary point or a regular singular point of a differential equation, we can immediately apply Padé approximants to the truncated series solutions. It provides us an approximate solution of the analytical one in the complex domain. For irregular singular points, formal series solutions are usually divergent. For divergent series, Padé approximants usually still works (see Appendix C.2), but sometimes, due to the Stokes phenomenon, we need the Borel–Padé summation method or numerical solutions.

2.4. Perturbation Theory

In this subsection, we consider perturbative expansions of eigenvalues in the quantum parameter . Since the Schrödinger Equation (1) is the form of singular perturbation in , we cannot naïvely apply the ordinary method of Rayleigh–Schrödinger perturbation theory. Normally singular perturbations are treated by the WKB method that is reviewed in the next subsection. Here, we use a more unconventional treatment. We will show that results from this method agree with those from the uniform WKB method in the next subsection.
We first rescale the variable by x = q . The Schrödinger Equation (1) is then written as
1 2 d 2 d q 2 + V ( q ) 2 ψ ( q ) = E 2 ψ ( q ) .
It is more convenient to define new variables by:
g = , v ( x ) = V ( x ) V ( 0 ) 2 , ϵ = E V ( 0 ) 2 .
The Schrödinger equation is now written as:
1 2 ψ ( q ) + v ( g q ) g 2 ψ ( q ) = ϵ ψ ( q ) .
In our analysis, the potential generically has the following form:
v ( g q ) g 2 = 1 g 2 m = 0 M g m v m ( g q ) .
Taking the limit g 0 with q kept finite, we can expand the potential functions as:
v 0 ( g q ) g 2 = 1 2 Ω 2 q 2 + j = 3 g j 2 v 0 , j q j , g m 2 v m ( g q ) = j = 1 g j + m 2 v m , j q j , m 1 ,
where Ω = v 0 ( 0 ) and v m , j = v m ( j ) ( 0 ) / j ! . We have used v ( 0 ) = 0 . Therefore (45) can be regarded as a perturbation of the harmonic oscillator. In this picture, we zoom in the potential minimum. Near the minimum the potential is a deformation of the harmonic potential.
Now we can apply perturbation theory. The eigenvalue at the n-th energy level is perturbatively expanded as:
ϵ n = l = 0 g l ϵ n , l , ϵ n , 0 = Ω n + 1 2 .
The usual textbook method is quite complicated to compute high-order corrections for the potential (47). Instead we use a smart way by Bender and Wu [31] and a generalization by Sulejmanpasic and Ünsal [32]. The idea is very simple. We put the following ansatz of the normalizable eigenfunction:
ψ n ( q ) = N n ( g ) e Ω q 2 / 2 l = 0 g l F n , l ( q ) .
In the perturbative computation, we can fix the overall factor N n ( g ) as we want. The zeroth order term F n , 0 ( q ) is of course the n-th Hermite polynomial H n ( Ω q ) for the unperturbed harmonic oscillator. The l-th order correction F n , l ( q ) is a polynomial of at most degree n + 3 l ,
F n , l ( q ) = k = 0 n + 3 l A n , l k q k ,
where we use a freedom of N n ( g ) so that
A n , l n = δ l 0 .
This ansatz is very powerful. In fact, the Schrödinger equation determines all the coefficients in the polynomial F n , k ( q ) as well as ϵ n , k recursively. This fact was first observed in the quartic anharmonic oscillator [31], and it was recently generalized in [32] to arbitrary potentials with a locally harmonic minimum.
We borrow the result in [32]. First, we have relations:
A n , l k = 1 2 Ω ( k n ) [ ( k + 2 ) ( k + 1 ) A n , l k + 2 + j = 1 l 1 2 ϵ n , j A n , l j k 2 j = 1 l v 0 , j A n , l j k j 2 2 j = 1 l m = 1 M v m , j A n , l j + 2 m k j ] ( n + 1 k n + 3 k ) ,
where we set A n , l k = 0 for k > n + 3 k . This determines A n , l k for k = n + 3 k , n + 3 k 1 , , n + 1 uniquely. Next, the energy correction is obtained by:
ϵ n , l = ( n + 2 ) ( n + 1 ) 2 A n , l n + 2 j = 1 l 1 ϵ n , j A n , l j n + j = 1 l v 0 , j A n , l j n j 2 + j = 1 l m = 1 M v m , j A n , l j + 2 m n j .
Finally, we fix the remaining coefficients by:
A n , l k = 1 2 Ω ( k n ) [ ( k + 2 ) ( k + 1 ) A n , l k + 2 + j = 1 l 2 ϵ n , j A n , l j k 2 j = 1 l v 0 , j A n , l j k j 2 2 j = 1 l m = 1 M v m , j A n , l j + 2 m k j ] ( 0 k n 1 ) .
These relations are easily solved by symbolic computation systems. In general, there are no odd order corrections in the energy: ϵ n , l = 0 for all odd l. The perturbative series of the original energy is then given by:
E n = V ( 0 ) + 2 ϵ n = V ( 0 ) + 2 Ω n + 1 2 + 2 l = 1 l ϵ n , 2 l .
Let us apply it to the Morse potential. We start with the rescaled eigen-Equation (10). After changing variables q g q and E = ( ε + 1 ) / ( 2 g ) , this is rewritten as:
1 2 d 2 d q 2 + 1 2 q 2 + k = 1 g k / 2 2 k + 2 2 2 ( k + 2 ) ! q k + 2 ψ ( q ) = E ψ ( q ) .
It turns out that the energy eigenvalue has only the first order perturbative correction:
E n = n + 1 2 g 2 n + 1 2 2 + 0 g 2 + 0 g 3 + 0 g 4 + .
Of course, this is a reflection of the hidden supersymmetric structure in the Morse potential. This cancellation is directly confirmed by the Mathematica package in [32] up to any desired orders in g. The result is consistent with the exact result (17). However, we should note that the wave function receives an infinite number of corrections. For instance, the ground state eigenfunction is given by:
ψ 0 ( q ) = e q 2 / 2 ( 1 g 1 / 2 6 q ( q 2 + 3 ) + g 72 q 2 ( q 4 + 3 q 2 + 9 ) g 3 / 2 6480 q 3 ( 5 q 6 + 54 q 2 + 135 ) + g 2 155520 q 4 ( 5 q 8 30 q 6 + 81 q 4 + 162 q 2 + 405 ) + ) .
It is a good exercise to compare it with the exact eigenfunction.
A merit of the perturbative series is that one can easily analytically continue the Planck parameter to the complex domain. Therefore, the resonant eigenvalues are obtained in this way. In the case of the Morse potential, the perturbative series accidentally stops at the first order, but this is not the case for most models. Moreover the perturbative series (48) is not convergent in general. Therefore we need to resume it by the Borel summation method (or Padé approximants). This point is important to obtain correct numerical values.

2.5. The WKB Method

The WKB method is a very powerful tool to investigate global properties of the wave function. Typically, it is used to analyze quantum tunneling effects. We review two WKB methods.

2.5.1. Standard WKB

Let us start with the original WKB method. We use (23). Since the equation takes the form of the singular perturbation in , we put the ansatz:11
ψ ( x ) = 1 p ( x ) exp i x p ( x ) d x ,
where p ( x ) satisfies the non-linear differential equation:12
p ( x ) 2 Q ( x ) + 2 p ( x ) 2 p ( x ) 3 p ( x ) 2 4 p ( x ) 2 = 0 .
We can solve (60) perturbatively in :
p ( x ) = k = 0 2 k p k ( x ) .
At the lowest order, we have two branches of the solution:
p 0 ± ( x ) = ± Q ( x ) = ± E V ( x ) .
For each branch, the quantum corrections are uniquely fixed. Hence these two branches leads to the two solutions of the Schrödinger equation. In general, these two are independent. Once one of them is constructed, the other is easily obtained. Therefore, we abbreviate the subscripts ±. The first two corrections in p ( x ) are given by:
p 1 = p 0 4 p 0 2 + 3 p 0 2 8 p 0 3 , p 2 = p 0 16 p 0 4 5 p 0 p 0 8 p 0 5 13 p 0 2 32 p 0 5 + 99 p 0 2 p 0 32 p 0 6 297 p 0 4 128 p 0 7 .
The higher order corrections rapidly get so lengthy.
For practical computations, a sophisticated way in [33] is very useful.13 Let us introduce the following notations:
q ( x ) : = p ( x ) p 0 ( x ) , D ^ : = 1 p 0 ( x ) d d x .
The non-linear Equation (60) now becomes:
1 q ( x ) 2 + 2 ϵ 0 ( x ) + q ( x ) 1 / 2 D ^ 2 q ( x ) 1 / 2 = 0 ,
where
ϵ 0 ( x ) : = p 0 ( x ) 3 / 2 d 2 d x 2 p 0 ( x ) 1 / 2 = 5 Q ( x ) 2 16 Q ( x ) 3 Q ( x ) 4 Q ( x ) 2 .
We solve (65) perturbatively,
q ( x ) = 1 + k = 1 2 k Y k ( x ) .
By construction we have p k ( x ) = Y k ( x ) p 0 ( x ) . The non-linear Equation (65) leads to a recurrence relation for Y k ( x ) [33]. This function can be separated as:
Y k ( x ) = Z k ( x ) + D ^ U k ( x ) .
Note that this separation is not unique. We require that Z k ( x ) gets as “simple” as possible, as shown below. It has the following advantage. In the WKB approach, a contour integral around two turning points.14
C ( a , b ) p k ( x ) d x ,
is particularly important. Here the contour C ( a , b ) is a closed curve encircling only the two turning points x = a and x = b . Using the above notation, this integral is written as:
C ( a , b ) p k ( x ) d x = C ( a , b ) Y k ( x ) p 0 ( x ) d x = C ( a , b ) Z k ( x ) p 0 ( x ) + d d x U k ( x ) d x = C ( a , b ) Z k ( x ) p 0 ( x ) d x .
Therefore, the derivative term D ^ U k ( x ) does not contribute to the contour integral. This fact drastically reduces the computational cost. For our purpose in spectral problems, it is sufficient to compute Z k . These are much simpler than p k . Up to order 12 ( k = 6 ), we have:
Z 1 = 1 2 ϵ 0 , Z 2 = 1 8 ϵ 0 2 , Z 3 = 1 32 ( 2 ϵ 0 3 ϵ 1 2 ) , Z 4 = 1 128 ( 5 ϵ 0 4 10 ϵ 0 ϵ 1 2 + ϵ 2 2 ) , Z 5 = 1 512 ( 14 ϵ 0 5 70 ϵ 0 2 ϵ 1 2 + 14 ϵ 0 ϵ 2 2 ϵ 3 2 ) , Z 6 = 1 2048 ( 42 ϵ 0 6 420 ϵ 0 3 ϵ 1 2 35 ϵ 1 4 + 126 ϵ 0 2 ϵ 2 2 + 20 ϵ 2 3 18 ϵ 0 ϵ 3 2 + ϵ 4 2 ) ,
where ϵ k : = D ^ k ϵ 0 . We perform the separation (68) so that Z k does not contain higher derivative functions ϵ k . This requirement uniquely fixes the separation. The explicit forms of U k are found in [33,35]. They are not important in our analysis.
Sometimes, the function Q ( x ) depends on explicitly. For example, in the next section, we encounter the situation such as:
Q ( x ) = Q 0 ( x ) + 2 Q 1 ( x ) .
Even in this case, we can still use the above formulae. Equation (60) is now written as:
p ( x ) 2 Q 0 ( x ) + 2 Q 1 ( x ) + p ( x ) 2 p ( x ) 3 p ( x ) 2 4 p ( x ) 2 = 0 .
Then the Equation (65) becomes:
1 q ( x ) 2 + 2 ϵ 0 ( x ) + Q 1 ( x ) p 0 ( x ) 2 + q ( x ) 1 / 2 D ^ 2 q ( x ) 1 / 2 = 0 ,
where p 0 ( x ) = Q 0 ( x ) . If we define:
ϵ ¯ 0 ( x ) : = ϵ 0 ( x ) + Q 1 ( x ) p 0 ( x ) 2 ,
then the formulae (71) still hold by replacing ϵ k ϵ ¯ k = D ^ k ϵ ¯ 0 .
In general, the turning points in the contour integral are complex, and we have to consider the WKB solutions in the complex domain. The WKB method still works in the complex domain [36,37]. In the resonant spectral problem, such a complex analysis is particularly important. Since the lowest function p 0 ( x ) is a multivalued function, we have to choose branch cuts carefully in numerical calculations.
Let us see the computations of the bound state energy and the resonant energy in the WKB method. To find these eigenvalues, we need to solve the connection problems at the semiclassical level. We start with the connection problem for the bound states shown in the left panel in Figure 6. There are two real turning points a < b . We use the lowest order WKB solution:
ψ ( x ) 1 p 0 ( x ) exp i x p 0 ( x ) d x .
Let us define:
w ( x 1 , x 2 ) : = 1 x 1 x 2 E V ( x ) d x , σ ( x 1 , x 2 ) : = 1 x 1 x 2 V ( x ) E d x .
The solution satisfying the boundary condition at x is given by:
ψ ( x ) 1 p 0 ( x ) e σ ( x , a ) ,
where we ignore irrelevant constants. Using the connection formula in [19], we can continue this solution to the region in x + . The result is:
1 p 0 ( x ) e σ ( x , a ) x ψ ( x ) x + cos ϕ p 0 ( x ) e σ ( b , x ) + 2 sin ϕ p 0 ( x ) e σ ( b , x ) ,
where
ϕ = π 2 w ( a , b ) .
In the region x + , the first term satisfies the boundary condition. We conclude that the bound state boundary condition requires:
sin ϕ = 0 w ( a , b ) = π n + 1 2 , n = 0 , 1 , 2 , .
where we have used w ( a , b ) > 0 for a < b . This equation is well-known as the Bohr–Sommerfeld quantization condition.
It seems hard to extend this argument when the quantum corrections are taken into account. Instead we rewrite the quantization condition as follows:
C ( a , b ) p 0 ( x ) d x = 2 π n + 1 2 ,
where p 0 ( x ) = Q ( x ) should be understood as a multi-valued complex function that defines a Riemann surface. Geometrically the left hand side gives an area surrounded by the curve p 2 + V ( x ) = E in the phase space ( x , p ) . This quantization condition is simply derived by the following argument. We assume that there are no singular points of the Schrödinger equation inside the closed curve C ( a , b ) . Then the wave function should be single-valued when x goes around C ( a , b ) . After encircling it, the WKB solution (76) apparently becomes
1 p 0 ( x ) exp i x p 0 ( x ) d x 1 p 0 ( x ) exp i C ( a , b ) p 0 ( x ) d x + i x p 0 ( x ) d x ,
where the minus sign comes from the multi-valuedness of p 0 ( x ) .15 The single-valued condition requires the quantization condition (82). This argument is easily extended to the formal WKB solution (59) [38]. The result is given by:
k = 0 2 k C ( a , b ) p k ( x ) d x 2 π n + 1 2 ,
or equivalently
k = 0 2 k C ( a , b ) Z k ( x ) p 0 ( x ) d x 2 π n + 1 2 .
Two remarks should be mentioned in order. As for usual perturbative series, the left hand side is a formal divergent series in . Therefore the equality in (84) and (85) should be understood in the asymptotic sense. The Borel resummation of the left hand side causes the Stokes phenomenon in the complex -plane. Therefore, in general, the left hand side may have non-perturbative corrections in . Its analysis is quite complicated. This is beyond the scope of this article. We refer to [34,37,39,40] on this issue. The other remark is that if there is a singular point of the Schrödinger equation inside C ( a , b ) , one has to consider a monodromy (see Appendix B.1). Then the required condition is a matching condition of two monodromies of the analytic solution and of the WKB solution. The monodromy of the WKB solution is given by (83), and our task is to know the monodromy of the true solution.
The argument of the resonant states is in parallel. See the right panel in Figure 6. Let us define:
w ˜ ( x 1 , x 2 ) : = 1 ˜ x 1 x 2 E ˜ V ˜ ( x ) d x , σ ˜ ( x 1 , x 2 ) : = 1 ˜ x 1 x 2 V ˜ ( x ) E ˜ d x .
In the case of the resonance, we start with the transmitted wave in x , and connect it to the waves in x + . The connection formula leads to:
T 1 p ˜ 0 ( x ) e i w ˜ ( x , a ˜ ) x ψ ˜ ( x ) x + 1 p ˜ 0 ( x ) e i w ˜ ( b ˜ , x ) + R 1 p ˜ 0 ( x ) e i w ˜ ( b ˜ , x ) ,
where a ˜ and b ˜ are two turning points, and
R = e i δ 1 + e 2 σ ˜ ( a ˜ , b ˜ ) , T = e σ ˜ ( a ˜ , b ˜ ) i δ 1 + e 2 σ ˜ ( a ˜ , b ˜ ) ,
with some phase shift δ that is not relevant in our analysis. See [19]. As seen in Section 2.1, the resonant states require 1 / R = 0 , and semiclassically, it yields:
1 + e 2 σ ˜ ( a ˜ , b ˜ ) = 0 .
Therefore, we obtain the Bohr–Sommerfeld condition for the resonant energy:
σ ˜ ( a ˜ , b ˜ ) = π i n + 1 2 , n = 0 , ± 1 , ± 2 , .
For the resonant states, σ ˜ ( a ˜ , b ˜ ) must be complex-valued. There are no reasons to exclude n < 0 . However, the spectra for n = m and for n = m 1 ( m = 0 , 1 , 2 , ) are symmetric.16 This is checked by taking the complex conjugate to (90). We can restrict on n 0 without loss of generality.
If two turning points a ˜ and b ˜ are both real, σ ˜ ( a ˜ , b ˜ ) is necessarily real. Therefore, the resonance condition is never satisfied. This implies that the resonant energy should be complex-valued. We have to carefully choose the correct complex turning points when we solve the quantization condition (90). A simple criterion is to trace continuous variations from real turning points to complex ones. To extend to the quantum corrected result, it is useful to rewrite the quantization condition as:
C ( a ˜ , b ˜ ) p ˜ 0 ( x ) d x = 2 π ˜ n + 1 2 .
It is straightforward to derive the quantum corrected condition:
k = 0 ˜ 2 k C ( a ˜ , b ˜ ) p ˜ k ( x ) d x 2 π ˜ n + 1 2 ,
or equivalently
k = 0 ˜ 2 k C ( a ˜ , b ˜ ) Z ˜ k ( x ) p ˜ 0 ( x ) d x 2 π ˜ n + 1 2 .
We have seen that the (all-order) Bohr–Sommerfeld quantization conditions for the bound states and the resonant states are almost same. These are related by the analytic continuation. This is a reason why we can expect the spectral relation (7).
For the Morse model (10), we compute the Bohr–Sommerfeld integral. The two turning points are given by:
a = log ( 1 1 + ε ) , b = log ( 1 + 1 + ε ) .
These are real as long as 1 ε < 0 . The classically allowed integral is exactly given by:
w ( a , b ) = 1 g a b ε e 2 q + 2 e q d q = π g ( 1 ε ) .
The Bohr–Sommerfeld condition turns out to give the exact result in this case:
π g ( 1 ε n ) = π n + 1 2 ε n = 1 g n + 1 2 2 .
The computation of the resonant states for the inverted potential is almost the same:
σ ˜ ( a ˜ , b ˜ ) = 1 g ˜ a ˜ b ˜ e 2 q + 2 e q ε ˜ d q = π g ˜ ( 1 ε ˜ ) ,
where17 a ˜ = log ( 1 1 ε ˜ ) , b ˜ = log ( 1 + 1 ε ˜ ) and
π g ˜ ( 1 ε ˜ n ) = π i n + 1 2 ε ˜ n = 1 i g ˜ n + 1 2 2 .
It would be interesting to check that the quantum corrections really do not contribute to the all-order Bohr–Sommerfeld quantization condition.

2.5.2. Uniform WKB

The uniform WKB method has a little bit of a different flavor. In the standard WKB method, it is not so easy (but possible) to derive the semiclassical perturbative expansion of the energy eigenvalue itself. In perturbation theory, this is easy, but it is rather hard (but possible) to see the global structure of the solutions because we zoom in the potential well/wall. The uniform WKB method has both advantages of these two approaches in a sense. This is particularly powerful in the analysis of non-perturbative instanton corrections to the spectrum. See [41,42] for instance. In this article, we just show a perturbative aspect to mention a relation to previous works in black hole physics.
The starting point is the following uniform WKB ansatz of the solution:
ψ ( x ) = 1 u ( x ) D ν u ( x ) ,
where D ν ( z ) is the parabolic cylinder function that satisfies the Weber equation:
y ( z ) + ν + 1 2 z 2 4 y ( z ) = 0 .
This differential equation is essentially same as the Schrödinger equation for the harmonic oscillator. The other solution is given by D ν 1 ( i z ) . Plugging the uniform ansatz into the Schrödinger Equation (1), one obtains:
1 4 u 2 u 2 + E V ν + 1 2 u 2 + 2 3 u 2 4 u 2 u 2 u = 0 .
For simplicity, we assume that the potential has the (global) minimum V min = 0 at x = 0 . This is easily realized by constant shifts of x and E. The equation is solved perturbatively:
u ( x ) = l = 0 l u l ( x ) , E = l = 1 l E l .
At the lowest order, we have:
u 0 ( x ) u 0 ( x ) = ± 2 V ( x ) .
It is formally solved by:
u 0 ( x ) 2 = 4 0 x V ( x ) d x ,
where we have chosen the positive branch. We have also required the condition u 0 ( 0 ) = 0 to fix the integration constant. The reason is as follows. If u 0 ( 0 ) 0 , then u 0 ( 0 ) = 0 because of V ( 0 ) = 0 and (103). In this case, the wave function is divergent at x = 0 . We exclude this unphysical case. As in the standard WKB, two branches of the square root in (104) lead to two independent solutions. We can choose the plus sign without loss of generality.
Let us expand the potential around x = 0 :
V ( x ) = k = 2 V k x k .
Then the solution u 0 ( x ) has the series representation:
u 0 ( x ) = 2 V 2 1 / 4 x + V 3 3 2 V 2 3 / 4 x 2 + 13 V 3 2 + 36 V 2 V 4 144 2 V 2 7 / 4 x 3 + O ( x 4 ) , x 0 .
At the first order, we have:
E 1 ν + 1 2 u 0 2 + 1 2 u 0 u 0 ( u 0 u 1 ) = 0 .
It is not easy to solve this equation for u 1 analytically. However, our goal is to get the spectrum. This is done by the regularity condition of u 1 at x = 0 . The result is:
E 1 = 2 V 2 ν + 1 2 .
Similarly at the second order, the regularity of u 2 at x = 0 as well as u 0 and u 1 determines:
E 2 = 7 + 60 α 2 32 V 3 2 V 2 2 + 3 ( 1 + 4 α 2 ) 8 V 4 V 2 ,
where α : = ν + 1 / 2 . In this way, we can fix the perturbative coefficients order by order. There is no difficulty to push the computation for higher orders.
So far, we do not require the normalizability of the wave function. It is realized by the special condition ν = n ( n = 0 , 1 , 2 , ) because in this case the parabolic cylinder function D ν ( z ) becomes normalizable. The perturbative coefficients of the energy eigenvalue should be compared with the result from the perturbative method in Section 2.4. It turns out that both are in agreement as expected. See Equation (2.7) in [20].
A very similar approach is found in the context of black hole perturbation theory [43,44,45]. There, the connection problem near the potential peak was solved by using the parabolic cylinder function. In black hole perturbation theory, the method in [43,44,45] is conventionally called the WKB approach. In this article, we refer to it as the uniform WKB approach to distinguish it from the standard WKB reviewed in the previous subsection. At the perturbative level, the uniform WKB method leads to the same result in perturbation theory. Its true value is revealed in non-perturbative analysis, which is beyond the scope of this article.

3. Applications to Black Hole Perturbation Theory

Now we proceed to applications to two kinds of problems in black hole perturbation theory. First, we review computations of the QNM frequencies of Schwarzschild black holes. They are well-studied and a very good playground for the applications. Next, we look at a mode stability problem of black strings in five-dimension. It is well-known that they show the Gregory–Laflamme instability. We show that the method in Section 2.3 allows us to evaluate a critical value of a parameter in this phase transition numerically. This method also judges a mode (in)stability of a given black hole easily.

3.1. Quasinormal Modes of Schwarzschild Black Holes

For notational simplicity, we abbreviate tildes in this subsection though we consider a resonant state problem. It will cause no confusions. Let us see the QNM problem of the four-dimensional Schwarzschild spacetime. We will present basics on linear perturbation theory of this geometry in Appendix A. There are also several excellent review articles on the QNMs [2,3,4,5,6]. The odd parity gravitational perturbation finally leads to the master equation with the effective potential (A91). As we will discuss in Appendix A, it is well-known that the odd parity and the even parity perturbations have the same QNM spectra. Since the former master equation is simpler than the latter one, we can concentrate our attention to the odd parity sector.
For later convenience, we use a dimensionless variable z = r / r H = r / ( 2 M ) rather than the radial variable r itself. The master equation then takes the form:
f ( z ) d d z f ( z ) d d z + ( 2 M ω ) 2 V s ( z ) ϕ ( z ) = 0 ,
where
f ( z ) = 1 1 z , V s ( z ) = f ( z ) ( + 1 ) z 2 + ( 1 s 2 ) 1 z 3 .
Here, s = 0 , 1 , 2 is the spin weight of the perturbing field. To see a direct connection to the Schrödinger equation, it is useful to introduce the so-called tortoise coordinate by:
d x d z = 1 f ( z ) x = z + log z 1 + c ,
where we fix the integration constant c so that the potential V s has a peak at x = 0 . This choice is not essential in the following analysis. In the tortoise coordinate x, the potential V s has the shape shown as the dashed line in the right panel of Figure 2. The analytic form in terms of x is complicated. The resonant boundary condition is then given by:
ϕ e 2 i M ω x ( z 1 ) 2 i M ω ( x ) , e + 2 i M ω x e + 2 i M ω z ( x + ) .
As we will see below, sometimes it is more convenient to use the original variable z rather than the tortoise variable x.
We apply previous methods to compute quasinormal mode frequencies (i.e., resonant energies), satisfying this boundary condition. Our purpose is to present widely applicable ways. The computation here is just an example. For instance, one can easily extend it to other spherically symmetric geometries such as the Reissner-Nördstrom solution. Although the Kerr solution is not spherically symmetric, the approach in this article is probably applicable even in this case. To do so, we note that an alternative master equation, which is isospectral to Teukolsky’s master equation [46], for the Kerr spacetime was recently conjectured in [47]. This master equation is expected to be useful for studying the QNM spectra.

3.1.1. Wronskian Method

We first apply the Wronskian method to this problem. For this purpose, we use the z-variable. Transforming the wave function by18
ϕ ( z ) = z 1 + s ( z 1 ) 2 i M ω e 2 i M ω z y ( z ) ,
then the new function y ( z ) satisfies the differential equation:
y ( z ) + γ z + δ z 1 + ϵ y ( z ) + α z q z ( z 1 ) y ( z ) = 0 ,
where
q = ( + 1 ) s ( s + 1 ) + 4 i M ω ( 1 + 2 s ) , α = ( 4 M ω ) 2 + 4 i M ω ( 1 + s ) , γ = 1 + 2 s , δ = 1 4 i M ω , ϵ = 4 i M ω .
This differential equation is well-known as the confluent Heun equation [48]. The relation between the master Equation (110) and the confluent Heun equation was discussed in great detail in [49]. We construct formal series solutions at z = 1 (event horizon) and z = (spacial infinity). The local solutions at z = 1 can be expressed by the local confluent Heun solution in [48], but we rather construct the series solutions directly for extensions to other cases.
The point z = 1 is a regular singular point of (115). Its exponents are ρ 1 = 0 and ρ 2 = 1 δ . Recalling the transformation (114), the QNM boundary condition at z = 1 is satisfied by the solution with ρ 1 = 0 . Therefore the Frobenius series solution we want is given by:
y 11 ( z ) = k = 0 a k ( z 1 ) k .
Plugging it into the confluent Heun equation, we obtain the following three-term recursion:
α k 1 a k + 1 + β k 1 a k + γ k 1 a k 1 = 0 ,
where
α k 1 = ( k + 1 ) ( k + δ ) , β k 1 = k 2 + k ( γ + δ + ϵ 1 ) + α q , γ k 1 = ( k 1 ) ϵ + α .
Using this recursion, we can easily evaluate the coefficients a k to very high orders. The series solution is convergent for | z 1 | < 1 in general. In the practical computation, we have to truncate the sum at a certain finite order. The standard method to reconstruct an (approximate) analytic solution from a truncated convergent series solution is Padé approximants. We review the Padé approximants in Appendix C.
Let us proceed to the solutions at infinity. The confluent Heun equation has the irregular singular point at z = . The Poincaré rank of this singular point is just r = 1 , and the asymptotic series solutions take the form (see Appendix B.2):
y formal ( z ) = e a z z b k = 0 c k z k .
The leading asymptotic behavior determines a and b as:
( a , b ) = 0 , α ϵ , ϵ , γ δ + α ϵ .
The QNM boundary condition at infinity is satisfied by the former case. Therefore, we seek the formal solution with the form:
y 1 formal ( z ) = z α / ϵ k = 0 c k z k .
The coefficients c k again satisfy the following three-term recursion:
α k c k + 1 + β k c k + γ k c k 1 = 0 ,
where
α k = ( k + 1 ) ϵ , β k = k 2 + k γ + δ ϵ 1 2 α ϵ + q α + α ( γ + δ 1 ) ϵ α 2 ϵ 2 , γ k = k 1 + α ϵ k γ + α ϵ .
A big difference from the solution near the horizon z = 1 is that this formal series is not convergent for any z. We have to truncate or resum it properly as in the Morse potential. As shown in Section 2.3, one way is the Borel summation. In practice, however, the computational cost is much saved by using the Padé approximants. We discuss an application of the Padé approximants to divergent series in Appendix C.2.
In summary, to evaluate the Wronskian practically, we use the Padé approximants y 11 [ M 1 / N 1 ] ( z ) and y 1 [ M / N ] ( z ) of the formal series solutions (117) and (122) by regarding z 1 and 1 / z as expansion parameters, respectively. We have a freedom to choose z 0 , at which we search zeros of the Wronskian. We fix it so that both the approximants y 11 [ M 1 / N 1 ] ( z 0 ) and y 1 [ M / N ] ( z 0 ) have the same numerical accuracy. Such accuracy is roughly estimated by | 1 y [ ( M 1 ) / ( N + 1 ) ] ( z 0 ) / y [ M / N ] ( z 0 ) | . Therefore our requirement for z 0 is:
| 1 y 11 [ ( M 1 1 ) / ( N 1 + 1 ) ] ( z 0 ) / y 11 [ M 1 / N 1 ] ( z 0 ) | | 1 y 1 [ ( M 1 ) / ( N + 1 ) ] ( z 0 ) / y 1 [ M / N ] ( z 0 ) | .
Let us show an explicit result. We set ( s , ) = ( 2 , 2 ) , and computed the series solutions (117) and (122) up to ( z 1 ) 120 and 1 / z 120 respectively for a given numerical value of M ω . By using the recurrence relations, they are very quickly done. We used the Padé approximants y 11 [ 60 / 60 ] ( z ) and y 1 [ 60 / 60 ] ( z ) . Setting z 0 = 8 , we searched zeros of the Wronskian W [ y 11 [ 60 / 60 ] , y 1 [ 60 / 60 ] ] by Newton’s method for an initial value M ω = 0.4 0.07 i . Then we finally get:
M ω s = 2 , = 2 QNM 0.373671684418041835793492 0.0889623156889356982804609 i .
Comparing it with a zero of W [ y 11 [ 59 / 61 ] , y 1 [ 59 / 61 ] ] , the accuracy estimation is about O ( 10 30 ) . Therefore, we can obtain quite good numerical values of the QNM eigenvalues. However, in this approach, it is not easy to identify the overtone number of the obtained QNMs. This identification is usually fixed by comparing with other approaches like the perturbative method or the (uniform) WKB method. In these methods, we can easily know the overtone number of the QNM, as will be seen below.
The Wronskian method presented here is similar to the well-known continued fraction method by Leaver [50]. We do not explain Leaver’s method in detail because it is discussed in many places. An advantage of the Wronskian method is that it still works even when we do not have an explicit recurrence relation. We just need two series (or numerical) solutions satisfying proper boundary conditions. This can be done without information on the recurrence relation in principle. Moreover, to use Leaver’s continued fraction method, one has to treat a recurrence relation carefully if it is not a three-term relation [51,52]. In this sense, the application range of the Wronskian method is wider than that of Leaver’s method. However, since Leaver’s method seems to be the best numerical method to obtain the precise QNM frequencies for the Schwarzschild spacetime, we can use it as a reference for comparison. For instance, the numerical result (125) shows about 30-digit agreement with Leaver’s result as expected.

3.1.2. WKB Method

We proceed to another method. In the WKB method, we have two options which variables x or z we use. It turns out that the z-variable looks more useful technically. The similar approach is found in [35].
To see a quantum mechanical aspect more clearly, we rewrite the master Equation (110) as follows. First, we divide it by ( + 1 ) / 2 ,
2 ( + 1 ) f ( z ) d d z f ( z ) d d z + 2 ( + 1 ) ( 2 M ω ) 2 f ( z ) 2 z 2 + 2 ( + 1 ) 1 s 2 z 3 ϕ ( z ) = 0 .
We regard the coefficient in front of the derivative term as the square of a Planck parameter,
: = 2 ( + 1 ) 1 / 2 .
Of course, this is not a unique way to introduce . Ultimately, one can freely introduce it by hand, and set = 1 at last. Such freedom is related to a choice of “base functions” in [33]. Here, we choose a very particular base function, which naturally relates to the multipole index . It has the advantage that we do not have to introduce an additional freedom of parameters.
Then, we obtain
2 f ( z ) d d z f ( z ) d d z + E f ( z ) 2 z 2 + 2 1 s 2 z 3 ϕ ( z ) = 0 ,
where
E : = 2 ( + 1 ) ( 2 M ω ) 2 .
In terms of the tortoise variable x, this is nothing but the Schrödinger equation. The potential receives the quantum correction in general. In our picture, the master equation for general spins s is regarded as a quantum deformation from that for the electromagnetic perturbation s = 1 . Note that this picture is just a computational technique.
The role of the multipole number is now mapped to (the inverse of) the Planck parameter. The eikonal limit corresponds to the classical limit 0 . We should note that = 0 is singular in our picture.19 We have to consider this case separately. The treatment in this particular case is, however, straightforward by introducing a formal Planck parameter by hand.
Next, we rewrite the differential equation as the so-called normal form. To do so, we transform the dependent variable by:
ψ ( z ) = f ( z ) ϕ ( z ) .
We finally obtain the Schrödinger-like equation for the z-variable:
2 d 2 d z 2 + Q 0 ( z ) + 2 Q 1 ( z ) ψ ( z ) = 0 ,
where
Q 0 ( z ) = 2 2 z + E z 3 z ( z 1 ) 2 , Q 1 ( z ) = 1 4 s 2 + 4 s 2 z 4 z 2 ( z 1 ) 2 .
Note that the domain we are interested in is z > 1 , not the whole real line.
We apply the WKB method to the differential Equation (131). This is essentially the same form as (23), but Q depends on explicitly. Even in this case, we can apply the WKB method, as was discussed in the previous section. The all-order Bohr–Sommerfeld quantization condition is given by (93). We start with:
p 0 ( z ) = Q 0 ( z ) = 2 2 z + E z 3 ( z 1 ) z .
Three zeros of Q 0 ( z ) are turning points. For 0 < E < 8 / 27 , two of them are in the regime z > 1 . Let us denote these two turning points by a ( E ) and b ( E ) with a ( E ) < b ( E ) , as shown in Figure 6. In the computation of the quasinormal mode frequencies, we need to analytically continue these turning points to the complex domain.
We want to compute the contour integral:
Ω 0 ( E ) : = C ( a , b ) p 0 ( z ) d z = C ( a ( E ) , b ( E ) ) 2 2 z + E z 3 ( z 1 ) z d z .
At first glance, it looks hard to perform the integral analytically. In principle, it is evaluated numerically by analytic continuations of a ( E ) and b ( E ) for complex values of E. In the numerical computation, we have to choose directions of branch cuts carefully.20 In our case, however, it is evaluated analytically. We need no integration formulae.
A basic strategy to obtain an analytic result is to look for a differential equation that Ω 0 ( E ) satisfies. To do so, we notice that the known function p 0 ( z ) satisfies the following partial differential equation:
8 E 2 ( 8 27 E ) 3 E 3 + 4 E ( 40 189 E ) 2 E 2 + 2 ( 16 123 E ) E + 15 p 0 = z 8 z ( 15 + 20 z 3 z 2 ( 2 + 3 E ) z 3 + E z 4 ) ( 2 2 z + E z 3 ) 3 / 2 .
Since the integration contour C ( a , b ) is a closed path that does not encircle the branch point z = 0 nor the other turning point, the contour integral of the right hand side trivially vanishes. As a result, we immediately find an ODE for Ω 0 ( E ) ,
8 E 2 ( 8 27 E ) d 3 d E 3 + 4 E ( 40 189 E ) d 2 d E 2 + 2 ( 16 123 E ) d d E + 15 Ω 0 ( E ) = 0 .
Such ODEs for contour integrals are known as Picard–Fuchs equations [53]. To find the Picard–Fuchs Equation (136), the partial differential Equation (135) for p 0 is crucial. This happens very non-trivially, and we do not have a systematic algorithm to find it. An approach in [53] is a hint to do this for general cases. Here we have put ansatz of differential operators appropriately.
Once we find the Picard–Fuchs Equation (136), it is not difficult to solve it. The standard technique is to construct series solutions at singular points, as reviewed in Appendix B. In the current case, we find analytic solutions in terms of generalized hypergeometric functions by using a new variable.21
ξ = 8 27 E .
The general solution is given by:
Ω 0 ( E ) = C 1 ξ + C 2 ξ 1 / 6 3 F 2 1 6 , 1 6 , 2 3 ; 1 3 , 5 3 ; ξ + C 3 ξ 5 / 6 3 F 2 5 6 , 5 6 , 4 3 ; 5 3 , 7 3 ; ξ .
There are three integration constants, which should be fixed by conditions of Ω 0 ( E ) with the integral form (134). This is obtained by considering the limit E 8 / 27 or ξ 1 . In this limit, the two turning points a and b meet together at z = 3 / 2 . The integrand p 0 ( z ) has the expansion around ξ = 1 ,
p 0 ( z ) = 2 z 3 3 ( z 1 ) 2 ( z + 3 ) 3 z 2 z 2 3 ( z 1 ) ( 2 z 3 ) 2 z 3 ( z + 3 ) ( ξ 1 ) + 2 z 2 ( z 3 9 z + 9 ) ( z 1 ) ( z + 3 ) ( 2 z 3 ) 3 2 z 3 ( z + 3 ) ( ξ 1 ) 2 + O ( ( ξ 1 ) 3 ) .
The turning points a and b are now shrinking in all orders in ξ 1 . Then, the integral around the turning points a and b is simply evaluated by the residue theorem:
Ω 0 ( E ) = 2 π i ( ξ 1 ) + 49 π i 36 2 ( ξ 1 ) 2 + O ( ( ξ 1 ) 3 ) , ξ 1 .
Our requirement for the integration constants is to reproduce this expansion from (138). With the help of Mathematica, we find the following exact values:
C 1 = 8 π i 3 2 3 , C 2 = i 2 3 2 Γ 2 ( 1 / 6 ) Γ ( 1 / 3 ) , C 3 = 3 i 8 3 2 Γ 2 ( 5 / 6 ) Γ ( 2 / 3 ) .
Once we find the exact form of Ω 0 ( E ) , it is a easy task to continue it to the complex E-plane. Our leading Bohr–Sommerfeld quantization condition is:
Ω 0 ( E ) = 2 π n + 1 2 , E = 2 ( 2 M ω ) 2 , = 2 ( + 1 ) 1 / 2 .
Note that there is no spin dependence at the leading order while the multipole dependence is included through the Planck parameter. This is expected in the eikonal limit. The spin dependence appears as “quantum corrections” to the potential. As mentioned above, at the leading order, the quantization condition should give an approximate value for s = 1 . The integer n characterizes quasinormal modes. It corresponds to the overtone number.
By solving the quantization condition (142) for = 1 and n = 0 , we find the following values of the frequencies:
M ω = 1 , n = 0 ( 0 ) 0.2675 0.09695 i ,
where the superscript ( 0 ) means the leading Bohr–Sommerfeld approximation. This value is compared to the QNM frequency for the s = 1 perturbation,
M ω s = 1 , = 1 , n = 0 QNM 0.248263264178 0.0924877179529 i .
The leading approximation already provides a relatively good numerical value. Of course, the semiclassical approximation should get better for the eikonal limit 1 . The = 1 case is the most “severe” situation.
The approximation is expected to be improved by including the quantum corrections as in (93). To check it, we need:
Ω k ( E ) = C ( a , b ) Z k ( z ) p 0 ( z ) d z ,
where Z k are given by (71) but we have to replace ϵ 0 by ϵ ¯ 0 in (75) due to the quantum correction Q 1 ( z ) to the potential.
It turns out that these integrals are also evaluated analytically in the Schwarzschild case. The idea is similar to the Picard–Fuchs equation. We look for certain differential operators acting on Ω 0 ( E ) . Let us see the computation for k = 1 in detail. We want to evaluate
Ω 1 ( E ) = C ( a , b ) Z 1 ( z ) p 0 ( z ) d z = 1 2 C ( a , b ) ϵ ¯ 0 ( z ) p 0 ( z ) d z ,
where
ϵ ¯ 0 ( z ) = p 0 ( z ) 3 / 2 d 2 d z 2 p 0 ( z ) 1 / 2 + Q 1 ( z ) p 0 ( z ) 2 .
We seek a good differential relation for the integrand Z 1 ( z ) p 0 ( z ) . We find the following one:
Z 1 ( z ) p 0 ( z ) = D 1 p 0 ( z ) + z a 1 , 0 + a 1 , 1 z + + a 1 , 5 z 5 z ( 2 2 z + E z 3 ) 3 / 2 ,
where D 1 is a second order differential operator,
D 1 = 1 12 ( 5 24 s 2 ) 9 8 ( 1 6 s 2 ) E 2 E 2 + 1 12 ( 1 12 s 2 ) 3 4 ( 1 6 s 2 ) E E + 1 16 + 3 32 ( 1 6 s 2 ) E ,
and a 1 , i are some coefficients. We do not know a general algorithm to find (148), but once we put such an ansatz correctly by trials and errors, we can fix all the coefficients uniquely. With the similar argument for the Picard–Fuchs equation, the relation (148) immediately leads to the contour integral relation:
Ω 1 ( E ) = D 1 Ω 0 ( E ) .
Since we know Ω 0 ( E ) exactly, we obtain an analytic form of Ω 1 ( E ) without any integration!
This nice property seems to exist in the higher corrections. We observe that there exist second order differential operators such as
Ω k ( E ) = D k Ω 0 ( E ) .
We confirmed this statement up to k = 6 .
Now we can evaluate the Bohr–Sommerfeld condition with the quantum corrections. We want to solve:
k = 0 2 k Ω k ( E ) 2 π n + 1 2 .
Of course, in the practical computation, we have to truncate the sum on the left hand side at a certain finite order. We set ( s , ) = ( 2 , 2 ) , for example. We solve the 12th order ( k = 6 ) quantization condition for n = 0 , and obtain the numerical eigenvalue,
M ω s = 2 , = 2 , n = 0 ( 12 ) 0.3736290 0.08894092 i .
This is indeed close to the correct value (125). The approximation is improved by using Padé approximants for the asymptotic series (152). If we use the [ 6 / 6 ] Padé approximant of the left hand side, we get:
M ω s = 2 , = 2 , n = 0 [ 6 / 6 ] 0.3736770 0.08896645 i .
The more detailed result is shown in Table 1.

3.1.3. Perturbative Method

Next, we use perturbation theory, reviewed in Section 2.4. Fortunately we have a natural quantum parameter in terms of the multipole index . We derive the perturbative expansions of the QNM frequencies in this quantum parameter. Compared with the WKB method, this method has an advantage that directly gives the small- expansion of the QNM frequency. As we have already mentioned, the same result is also obtained by the uniform WKB method in Section 2.5.
For this purpose, the tortoise variable (112) is useful. We rewrite (128) as:
2 d 2 d x 2 + V ( x ) ϕ ( x ) = E ϕ ( x ) ,
where
V ( x ) = V 0 ( x ) + 2 V 1 ( x ) = f ( z ) 2 z 2 + 2 1 s 2 z 3 .
The potential now depends on , but we can use the results in the previous section. V 0 ( x ) has the maximum at z = 3 / 2 . We fix the integration constant c in (112) so that z = 3 / 2 is mapped to x = 0 . Hence we have c = 3 / 2 + log 2 . We then expand the potential V 0 ( x ) and V 1 ( x ) around x = 0 .
V 0 ( x ) = 8 27 32 x 2 729 + 128 x 3 19683 + 256 x 4 59049 256 x 5 177147 2048 x 6 7971615 + O ( x 7 ) , V 1 ( x ) = ( 1 s 2 ) 8 81 16 x 729 32 x 2 2187 + 128 x 3 19683 + 448 x 4 531441 8192 x 5 7971615 + 17408 x 6 215233605 + O ( x 7 ) .
Recall that we are considering the resonant problem for the QNMs. Strictly speaking, to use perturbation theory in Section 2.4, we need to analytically continue the Planck parameter i in order to invert the potential. See [20] for detail. However, the result in Section 2.4 can be formally used if we instead consider the analytic continuation of Ω . Let us see it in detail.
The Equation (155) can be written as:
1 2 ψ ( q ) + v ( g q ) g 2 ψ ( q ) = ϵ ψ ( q ) ,
where g = and x = g q ,
v ( g q ) g 2 = v 0 ( g q ) g 2 + g 2 v 4 ( g q ) , v 0 ( x ) = V 0 ( x ) V 0 ( 0 ) 2 , v 4 ( x ) = V 1 ( x ) V 1 ( 0 ) 2 ,
and
ϵ = E V ( 0 ) 2 .
The perturbative series of ϵ can be computed by the method reviewed in Section 2.4. The angular frequency Ω is given by:
Ω 2 = 2 v 0 ( 0 ) = 32 729 ,
which is now purely imaginary due to the resonant problem. However we can regard Ω as a formal parameter during the computation.
Here, we focus on the lowest overtone n = 0 for simplicity. The perturbative corrections are then computed as follows. For l = 1 , 2 , , we start with the recurrence relation (see (52)),
A 0 , l k = 1 2 Ω k [ ( k + 2 ) ( k + 1 ) A 0 , l k + 2 + j = 1 l 1 2 ϵ 0 , j A 0 , l j k 2 j = 1 l v 0 , j A 0 , l j k j 2 2 j = 1 l v 4 , j A 0 , l j 2 k j ] ( 1 k 3 l ) ,
where we regard A 0 , l k = 0 for k > 3 l and for k < 0 . This determines A 0 , l k for k = 3 l , 3 l 1 , , 1 in this order. Then the energy correction is obtained by:
ϵ 0 , l = A 0 , l 2 .
Recall that we set A 0 , l 0 = δ l 0 .
For l = 1 , the recurrence relation (162) gives:
A 0 , 1 3 = 64 59049 Ω , A 0 , 1 2 = 0 , A 0 , 1 1 = 64 19683 Ω 2 , A 0 , 1 0 = 0 .
Then we find ϵ 0 , 1 = 0 .
For l = 2 , we have the following non-zero components:
A 0 , 2 6 = 2048 3486784401 Ω 2 , A 0 , 2 4 = 32 59049 Ω + 5632 1162261467 Ω 3 , A 0 , 2 2 = 32 19683 Ω 2 + 5632 387420489 Ω 4 ,
and
ϵ 0 , 2 = 32 19683 Ω 2 5632 387420489 Ω 4 .
Up to this order, the results do not depend on s, but the spin-dependence appears for l 3 . In this way, we can easily fix ϵ 0 , l up to quite high orders.
The perturbative series of the fundamental QNM frequency is given by:
( M ω n = 0 pert ) 2 = ( + 1 ) 8 V ( 0 ) + Ω + 2 l = 1 l ϵ 0 , 2 l .
Plugging Ω = 4 2 i / 27 and s = 2 , we finally get22
( M ω s = 2 , n = 0 pert ) 2 = ( + 1 ) 8 [ 8 27 4 2 i 27 281 2 729 + 6163 i 3 26244 2 + 969319 4 17006112 + 60687295 i 5 2448880128 2 + 21306926297 6 528958107648 4107036251635 i 7 114254951251968 2 + 172267964131073 8 24679069470425088 + O ( 9 ) ] .
Recall that = 2 / [ ( + 1 ) ] . The higher order corrections are easily and quickly computed by using Mathematica.
Note that this result is closely related to the uniform WKB method in [43,44,45,54]. The result can also be compared with an expansion in [55]23 by re-expanding (168) in terms of L 1 : = ( + 1 / 2 ) 1 . As we have mentioned many times, the semiclassical perturbative series (168) is in general an asymptotic divergent series, and we have to truncate the sum at the optimal order or to resum it by the Borel(-Padé) summation [20,21] or the Padé approximants [56,57,58].
For = 2 , we show the numerical values of M ω s = 2 , = 2 , n = 0 pert for lower orders in Table 1. Of course, as in the WKB analysis, the perturbative approximation is better for larger .
An advantage of this method is that it is widely applicable to many potentials. In fact, to compute the QNM frequency we need only the information near the potential peak (i.e., the Taylor series around it). The global information is encoded in high-order corrections in principle. Moreover, the overtone number n of the QNMs naturally appears as the quantum number of the bound states.
The application range of perturbation theory is quite wide. One can consider other types of perturbative expansions of the QNM spectra. For instance, some parameters of black holes can be regarded as deformation parameters from simpler geometries (e.g., the Schwarzschild spacetime). The Reissner-Nördstrom spacetime is regarded as a deformation of the Schwarzschild one by the electric charge, and the Kerr case is also regarded as a deformation by the angular momentum. Then one can consider the perturbative expansions of the QNM spectra for such deformation parameters [47,59,60,61,62].

3.1.4. Comments on Even Parity Perturbation

So far, we have looked into the analysis in the odd parity sector as well as in the scalar and the electromagnetic perturbations. As discussed in Appendix A, the even parity gravitational perturbation of the Schwarzschild spacetime leads to the Zerilli Equation (A105). Since the QNM spectra of the Zerilli equation is exactly the same as those of the Regge–Wheeler equation, it is sufficient to study either of them. However we should comment that our formulation is also applicable to the Zerilli equation. We briefly see it.
In our convention in this section, the Zerilli equation is written as:
f ( z ) d d z f ( z ) d d z + ( 2 M ω ) 2 V + ( z ) ϕ ( z ) = 0 ,
where
f ( z ) = 1 1 z , V + ( z ) = f ( z ) λ 2 ( λ + 2 ) z 3 + 3 λ 2 z 2 + 9 λ z + 9 z 3 ( λ z + 3 ) 2 ,
with λ = 2 + 2 . This ODE has regular singular points at z = 0 , 1 , 3 / λ and an irregular singular point at z = with Poincaré rank 1. Then, the local series solutions near z = 1 and z = can be constructed straightforwardly. This means that the Wronskian method works in this case.
The semiclassical analyses are also applied. To see them, we rewrite the Zerilli potential as:
V + ( z ) = f ( z ) ( + 1 ) z 2 + δ V + ( z ) , δ V + ( z ) = 3 ( λ ( λ + 4 ) z 2 + 6 z 3 ) z 3 ( λ z + 3 ) 2 .
In the eikonal limit , we can see δ V + ( z ) O ( 1 ) , and then this part can be regarded as the “quantum” correction, as in the Regge–Wheeler potential. We can play the same game. In particular, the perturbative correction is computed by expanding the potential as:
V + ( z ) = f ( z ) 2 2 z 2 2 3 z 3 4 3 ( 2 z 3 ) z 4 6 3 ( 4 z 2 15 z + 12 ) 2 z 5 + O ( 8 ) ,
with = 2 / [ ( + 1 ) ] . We can use the general formulation in Section 2.4. Note that, to obtain a perturbative correction to the eigenvalue, we can truncate the sum in the potential (172) to a certain order. It is a good exercise to check that the perturbative correction to the QNM frequency for this potential agrees with (168) for the Regge–Wheeler potential. This is nothing but the isospectrality of the Regge–Wheeler/Zerilli potential. See Appendix A.4.7.

3.2. Mode Stability Analysis: Gregory–Laflamme Instability for Black Strings

In this subsection, we discuss the mode stability problem in black hole perturbation theory. If the linear perturbation of a spacetime geometry has an exponentially growing mode, then it is unstable. In the language of master equations, there exists a bound state solution with a negative energy ( E = ω 2 < 0 ) that corresponds to an unstable mode. As we reviewed in Section 2.3, the number of bound states of the Schrödinger type equation is captured by counting functions. We can use them to judge whether the linear perturbation of a given geometry is stable or not.
Here we use Milne’s counting function (33) for this purpose. The number of bound states below an energy E are computed by N Milne ( E ) . Therefore the mode stability condition is simply presented by:
N Milne ( 0 ) < 0 .
If this condition is satisfied, then there is no bound state for E < 0 .
Let us consider a concrete example. We discuss the mode (in)stability of black strings in five-dimensions. It is well-known that this black string shows the so-called Gregory–Laflamme instability. We show that Milne’s counting function actually explains the critical value of this transition. The metric of the black string solution in five-dimension is given by
d s 2 = f ( r ) d t 2 + d r 2 f ( r ) + r 2 ( d θ 2 + sin 2 θ d ϕ 2 ) + d z 2 ,
where f ( r ) = 1 r H / r , and the z-direction is compactified. The master equation reduces to the Schrödinger type equation,
d 2 d x 2 + ω 2 V ( x ) ψ ( x ) = 0 ,
where x is the tortoise variable, defined by:
x = r + r H log r r H 1 + c .
The integration constant c is fixed as follows.
The potential is explicitly given by:
V ( x ) = f ( r ) r 3 k 6 r 9 + 6 k 4 r 7 9 k 4 r 6 r H 12 k 2 r 4 r H + 9 k 2 r 3 r H 2 + r H 3 ( k 2 r 3 + r H ) 2 ,
where k is the wave number along the z-direction. For k 2 < 2 + 3 , the potential has a global minimum in r > r H . V 0 in x , and V k 2 in x . We focus on the case 0 < k 2 < 2 + 3 , and fix the integration constant c so that the potential V ( x ) has a minimum at x = 0 . The typical shape of V ( x ) is shown in the left panel of Figure 7.
We first solve the differential equation:
d 2 d x 2 + E V ( x ) w E ( x ) = 1 w E ( x ) 3 , w E ( 0 ) = 1 , w E ( 1 ) = 0 ,
where E = ω 2 . Milne’s counting function is then obtained by:
N Milne ( E ) = 1 π d x w E ( x ) 2 1 .
Of course, the non-linear differential Equation (178) cannot be solved analytically. We numerically solve it for a finite segment [ x min , x max ] . We finally want to evaluate N Milne ( 0 ) . We have to care about its numerical evaluation. In the region x > x max , the potential is greater than zero-energy. See the left of Figure 7. Since 1 / w 0 ( x ) 2 very rapidly decays in x > x max , we can ignore the contribution to N Milne ( 0 ) from this region as long as V ( ) = k 2 > 0 . On the other hand, in x < x min , the potential asymptotically goes to zero. The contribution from this region is not exponentially small for E = 0 . We need to evaluate this contribution approximately. If | x min | is large, we can ignore the potential. Therefore the approximate equation at E = 0 in x < x min is
w ¯ 0 ( x ) = 1 w ¯ 0 ( x ) 3 .
This equation can be solved analytically,
w ¯ 0 ( x ) 2 = C 1 + ( x + C 2 ) 2 C 1 ,
where C 1 and C 2 are integration constants. We have to connect it, at x = x min , with the (numerical) solution in x min x x max smoothly. This determines the integration constants. The result is the following:
C 1 = w 0 ( x min ) 2 1 + w 0 ( x min ) 2 w 0 ( x min ) 2 , C 2 = x min + w 0 ( x min ) 3 w 0 ( x min ) 1 + w 0 ( x min ) 2 w 0 ( x min ) 2 ,
where w 0 ( x min ) and w 0 ( x min ) are evaluated by solving (178) at E = 0 numerically. Therefore the counting function at E = 0 is approximately evaluated by
N Milne ( 0 ) = 1 π d x w 0 ( x ) 2 1 1 π x min x max d x w 0 ( x ) 2 + 1 π x min d x w ¯ 0 ( x ) 2 1 ,
where the integration over ( , x min ) can be computed exactly,
x min d x w ¯ 0 ( x ) 2 = x min d x C 1 C 1 2 + ( x + C 2 ) 2 = arctan [ w 0 ( x min ) w 0 ( x min ) ] + π 2 .
Our final expression is thus given by:
N Milne ( 0 ) 1 π x min x max d x w 0 ( x ) 2 + 1 π arctan [ w 0 ( x min ) w 0 ( x min ) ] 1 2 .
The remaining integral is evaluated by solving (178) with E = 0 numerically for x min x x max .
Now we apply this formalism to the black string potential (177). Technically, it is not hard to solve (178) numerically for a given potential because it is an initial value problem. In the right panel of Figure 7, we show N Milne ( 0 ) against the parameter k. We set x min = 100 and x max = 100 . It is clear to see that there is a value at which the sign of N Milne ( 0 ) changes. This is nothing but the Gregory–Laflamme instability. Our estimation of the critical value k c is:
k c 0.8762 .
For k < k c , there exists a negative bound state mode, and we conclude that the black string solution is unstable.
We note that there is another efficient way to discuss the mode stability of black holes [63,64,65,66]. The approach in this article is effectively useful to see the (in)stability quickly. On the other hand, the method in [63,64,65,66] is useful to show the (in)stability more rigorously.

4. Concluding Remarks

In this article, we reviewed various techniques in quantum mechanics, and explained how these techniques are applied to spectral problems appearing in black hole perturbation theory. Each of them has its own advantage, and it would be important to see the same problem from various perspectives. We hope that the techniques in this article will be useful for future developments of black hole perturbation theory.
Our main interest is QNM frequencies, which is particularly important in recent gravitational wave observations. As spectral problems, they are understood as resonant eigenvalues for Schrödinger-type wave equations. We showed various methods to compute the QNM frequencies with high precision both numerically and analytically.
We also proposed a simple criteria to see the mode (in)stability of the linear perturbation of a given geometry. The problem is mapped to an existence condition of negative energy bound states. We can easily write it down by using Milne’s very old idea in [26]. As an application, we explicitly showed the Gregory–Laflamme instability of 5D black string solutions.
In theories beyond general relativity, one often encounters coupled master equations. There are some approaches to analyze such systems [60,66,67,68,69,70,71,72,73]. It would be interesting to develop the techniques in this article for the coupled master equations.
Although we did not review them in this article, there are also new perspectives to study the black hole spectral problems from four-dimensional supersymmetric gauge theories [47,74,75] and from two-dimensional conformal filed theories [76,77,78,79,80].24 The ideas are based on the common mathematical structure in the spectral problems. These new perspectives provide us powerful analytic treatments for the spectral problems. However, a nice physical interpretation of the relation is still obscure.

Author Contributions

These authors (Y.H. and M.K.) contribute equally to this paper. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by JSPS KAKENHI Grant Numbers JP18K03657 (YH) and JP20H04746 (MK).

Data Availability Statement

Not applicable.

Acknowledgments

We are grateful to Hiroyuki Nakano for reading the draft carefully.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Derivation of the Master Equations around the Four-Dimensional Schwarzschild Spacetime

In this appendix, we review a derivation of the master equations around the Schwarzschild spacetime [82,83]. For simplicity, we assume that the spacetime is four dimensional and it has the global spherical symmetry. The basic technique used in this section can also be applied for the various cases such as modified gravity theories. See also reviews and textbooks [1,3,7,8,9]. For further general cases, which treat higher dimensions and the maximally symmetric angular part, are discussed, for example, in [84,85].

Appendix A.1. Four-Dimensional Schwarzschild Spacetime and Killing Vectors

The metric of the Schwarzschild spacetime is given by
d s 2 = g μ ν ( Sch ) d x μ d x ν = f d t 2 + d r 2 f + r 2 ( d θ 2 + sin 2 θ d ϕ 2 ) ,
with f = 1 r H / r . This spacetime admits four linearly independent Killing vectors
ξ ( t ) μ x μ = t ,
ξ ( 1 ) μ x μ = sin ϕ θ cos θ cos ϕ sin θ ϕ ,
ξ ( 2 ) μ x μ = cos ϕ θ cos θ sin ϕ sin θ ϕ ,
ξ ( 3 ) μ x μ = ϕ ,
where those vectors satisfy the Killing equations μ ξ μ + ν ξ μ = 0 . The generator of the spherical symmetry ξ ( I ) μ with I = 1 , 2 , 3 form the S O ( 3 ) algebraic relations
[ L ξ ( I ) , L ξ ( J ) ] = ϵ I J K L ξ ( K ) ,
where I , J , K run over 1 , 2 , 3 , the operator L ξ ( I ) denotes the Lie derivative along ξ ( I ) μ , and ϵ I J K is the totally anti-symmetric symbol with ϵ 123 = 1 . Defining the Casimir operator L 2 as
L 2 : = I = 1 3 L ξ ( I ) L ξ ( I ) ,
the relations
[ L 2 , L ξ ( I ) ] = 0 ,
hold for I = 1 , 2 , 3 .

Appendix A.2. Massless Klein–Gordon Equation

First, we consider the massless Klein–Gordon equation on the Schwarzschild spacetime
Φ = 0 ,
where = μ μ . Before showing the explicit calculation, we discuss the general relation among the separation of variables and Killing vectors.

Appendix A.2.1. Killing Vector and Separation of Variable

When a spacetime admits a Killing vector ξ μ , we can choose the coordinate system so that ξ μ becomes one of the coordinate basis. We denote the corresponding coordinate by λ , then
ξ μ x μ = λ ,
holds. In particular, the Killing equations L ξ g μ ν = 0 become
g μ ν λ = 0 ,
in this coordinate system and the relation [ / λ , ] = 0 holds because the d’Alembertian, which acts on scalar fields, = ( 1 / | g | ) μ ( | g | g μ ν ν ) does not depend on λ . If we assume the form of the Klein–Gordon field Φ as
Φ = e i L λ Φ ˜ ( x i ) ,
where L is the separation constant and x i denote the coordinates except for λ , we obtain a relation
λ Φ = Φ λ = i L e i L λ Φ ˜ ( x i ) = i L Φ .
The solution of this equation can be written by
Φ = e i L λ F ( x i ) ,
where F is a function of x i . Thus, we conclude that the Klein–Gordon equation reduces to the equation F ( x i ) = 0 which only contains x i , while the explicit form of the function F ( x i ) is not known.
If the spacetime admits two or more commutative Killing vectors, we can choose coordinate system so that these Killing vectors are coordinate bases, and the same discussion as single Killing vector case holds by using the method of separation of variables with respect to the corresponding coordinates.

Appendix A.2.2. Separation of Variable Using SO(3) Symmetry

While the S O ( 3 ) Killing vectors ξ ( I ) μ are not commutative, we can consider the simultaneous eigenfunction of one of the Killing vectors and L 2 . We denote the simultaneous eigenfunction with respect to L ξ ( t ) , L ξ ( 3 ) and the Casimir operator L 2 by e i ω t e i m ϕ S ( θ ) where ω and m are constants, and we assume the eigenvalue of L 2 is ( + 1 ) . The operator L 2 satisfies
[ L 2 , ] = 0 ,
because L 2 = I = 1 3 L ξ ( I ) L ξ ( I ) and [ L ξ ( I ) , ] = 0 . Assuming the form of the Klein–Gordon field Φ as
Φ = e i ω t e i m ϕ S ( θ ) Φ ˜ ( r ) ,
the relations
L ξ ( t ) Φ = i ω Φ ,
L ξ ( 3 ) Φ = i m Φ ,
L 2 Φ = ( + 1 ) Φ ,
hold. Because the operator L 2 is the Laplacian on S 2 , L 2 Φ = θ 2 Φ + cot θ θ Φ + csc 2 θ ϕ 2 Φ , the regular solutions to the above equations exist only when = 0 , 1 , 2 , and m = 0 , ± 1 , ± 2 , , ± from the Sturm-Liouville theory, then e i m ϕ S ( θ ) is the spherical harmonics Y m
e i m ϕ S ( θ ) = Y m .
Defining L ± = i L ξ ( 1 ) ± L ξ ( 2 ) , the relations
L ± Y m = ( m ) ( ± m + 1 ) Y m ± 1 ,
hold for different modes m. Using the commutation relations [ L ξ ( t ) , ] = [ L ξ ( 3 ) , ] = [ L 2 , ] = 0 , we can also show that Φ satisfies
L ξ ( t ) Φ = i ω Φ ,
L ξ ( 3 ) Φ = i m Φ ,
L 2 Φ = ( + 1 ) Φ .
The regular solution of the above equations can be written by
Φ = e i ω t Y m F ( r ) ,
where F ( r ) is a function of r. This implies that the Klein–Gordon equation reduces to an ordinary differential equation F ( r ) = 0 which only depends on r.25

Appendix A.2.3. Explicit Calculation

We write the form of the Klein–Gordon field Φ as
Φ = , m Φ m ,
with
Φ m = e i ω t Y m ( θ , ϕ ) Φ ( 0 ) ( r ) r ,
where we assumed the time dependence as e i ω t and omitted to write indices , m in Φ ( 0 ) . Because Equation (A25) holds for each mode , m , and the spherical harmonics Y m with different indices are linearly independent, we can separately discuss different , m modes. The Klein–Gordon equation for each mode reduces to26
Φ m = e i ω t Y m r f f d d r f d Φ ( 0 ) d r + ω 2 V ( 0 ) Φ ( 0 ) = 0 ,
where the effective potential V ( 0 ) is given by
V ( 0 ) = f ( + 1 ) r 2 + r H r 3 .
Thus, we obtain the master equation
f d d r f d Φ ( 0 ) d r + ω 2 V ( 0 ) Φ ( 0 ) = 0 .
We comment that if we consider the massive Klein–Gordon equation ( m 2 ) Φ = 0 , the effective potential is changed as V ( 0 ) V ( 0 ) + f m 2 .

Appendix A.3. Vector Harmonics on S2

Before discussing the case of the Einstein equations, we introduce some mathematical arguments.

Appendix A.3.1. Helmholtz–Hodge Decomposition

We introduce the Helmholtz–Hodge decomposition of vector and symmetric two tensor fields on S 2 .27 Let γ i j be the metric of a unit two sphere S 2
γ i j d x i d x j = d θ 2 + sin 2 θ d ϕ 2 .
In this subsection, i , j denote the tensor indices on S 2 , and we raise or lower tensor indices by γ i j . We denote by ^ i the covariant derivative with respect to γ i j .
Any regular vector field v i on S 2 can be uniquely decomposed into the scalar part and the vector part as
v i = ^ i S + V i ,
with
^ i V i = 0 .
The proof is simple. Taking the divergence of Equation (A32), we obtain an equation ^ i v i = ^ i ^ i S . Solving this Poisson equation on S 2 , ^ i S is uniquely determined.28 Then, V i is given by V i = v i ^ i S which clearly satisfies ^ i V i = 0 .
In a similar way, any regular symmetric tensor field t i j on S 2 can be uniquely decomposed as29
t i j = ^ i ^ j 1 2 γ i j ^ t T + t L γ i j + 2 ^ ( i t j ) v ,
with
^ i t i v = 0 .
To prove this, first, taking the trace of Equation (A34), we obtain t L = t i i / 2 . Next, acting an operator ^ i ^ j on Equation (A34), after some calculation, we obtain
^ ^ t T + 2 ^ t T = 2 ^ i ^ j ( t i j t L γ i j ) .
Solving Equation (A36) with respect to ^ t T , an equation ^ t T = A with a regular function A is obtained, and then we obtain t T by solving this equation. Then, the form of ( ^ i ^ j ( 1 / 2 ) γ i j ^ ) t T is uniquely determined.30 Finally, taking a divergence of Equation (A34), we obtain an equation
^ i t i j ^ i ^ j 1 2 γ i j ^ t T t L γ i j = ^ t j v + t j v .
The solution of Equation (A37) with ^ i t i v = 0 is uniquely determined except for the homogenous solutions. Because the homogenous solutions correspond to the Killing vectors on S 2 , they do not affect the form of ^ ( i t j ) v .

Appendix A.3.2. Vector Harmonics

It is well known that any regular function on S 2 can be written by the spherical harmonics Y m . Here, we introduce the vector spherical harmonics V i m which can express any divergence free regular vector field on S 2 .31 The vector spherical harmonics V i m satisfies the equations32
L ξ ( 3 ) V i m = i m V i m ,
L 2 V i m = ( + 1 ) V i m ,
with
^ i V i m = 0 ,
where = 1 , 2 , , and m = 0 , ± 1 , ± 2 , , ± . We can explicitly construct V i m by using the spherical harmonics Y m as
V i m = ϵ ^ i j ^ j Y m ,
where ϵ ^ i j is the Levi–Civita tensor of γ i j .33
Using the vector spherical harmonics, any divergence free vector V i on S 2 can be written by
V i = m d m V i m ,
where d m are constants. Then, a vector field v i in Equation (A32) becomes
v i = ^ i S + V i = , m v i m = , m c m ^ i Y m + d m V i m .
We note that each , m component v i m satisfies
L ξ ( 3 ) v i m = i m v i m ,
L 2 v i m = ( + 1 ) v i m .
On the other hand, if v i satisfies Equations (A44) and (A45), v i should take the form of
v i = v i m = c m ^ i Y m + d m V i m .
This can be shown as follows. Substituting Equation (A43) into Equations (A44) and (A45), we obtain equations ¯ , m ¯ ( m ¯ m ) v i ¯ m ¯ = 0 and ¯ , m ¯ ( ¯ ( ¯ + 1 ) ( + 1 ) ) v i ¯ m ¯ = 0 . Then v i ¯ m ¯ with ¯ and m ¯ m should vanish because of the linear independence of Y m and V i m .

Appendix A.4. Einstein Equations

Appendix A.4.1. Linear Metric Perturbation

Next, we consider the case of the vacuum Einstein equations.34
We denote by ϵ h μ ν the small deviation from the Schwarzschild metric where ϵ is a small parameter. The total metric at O ( ϵ ) becomes
g μ ν = g μ ν ( Sch ) + ϵ h μ ν .
We need to solve the vacuum Einstein equations G μ ν = 0 or equivalently R μ ν = 0 at O ( ϵ )
G μ ν = ϵ 1 2 μ ν h α α 1 2 α α h μ ν + α ( μ h ν ) α + 1 2 g μ ν ( α α h β β α β h α β ) = : ϵ X ^ h μ ν = 0 ,
where μ denotes the covariant derivative with respect to the Schwarzschild metric, and X ^ denotes the operator s.t. X ^ h μ ν gives the linearized Einstein tensor.
The perturbed metric h μ ν can be written by
h μ ν d x μ d x ν = h A B d x A d x B + h A i d x A d x i + h i j d x i d x j ,
where A , B run t , r , and i , j run θ , ϕ . We can regard h A B , h A i , h i j as scalar, vector and symmetric tensor on S 2 , respectively. Then, using the Helmholtz–Hodge decomposition, perturbed metric h μ ν can be expressed by Y m and V i m
h μ ν d x μ d x ν = , m h μ ν m d x μ d x ν = , m h μ ν ( + ) m + h μ ν ( ) m d x μ d x ν ,
where h μ ν ( + ) m is the even parity perturbation
h μ ν ( + ) m d x μ d x ν = e i ω t H 0 d t 2 + 2 H 1 d t d r + H 2 d r 2 Y m 2 r e i ω t H t Ω d t + H r Ω d r ^ i Y m d x i ( + 1 ) + 2 r 2 e i ω t K γ i j Y m + K 2 ( + 1 ) i ^ ^ j Y m 1 2 γ i j ^ Y m d x i d x j ,
and h μ ν ( ) m is the odd parity perturbation
h μ ν ( ) m d x μ d x ν = 2 e i ω t ( h 0 d t + h 1 d r ) V i m d x i 2 r 2 e i ω t h Ω ^ ( i V j ) m d x i d x j ( + 1 ) 1 .
We assumed the time dependence of h μ ν as e i ω t . There are seven unknown functions H 0 , H 1 , H 2 , H t Ω , H r Ω , K , K 2 for the even parity perturbation, and three unknown functions h 0 , h 1 , h Ω for the odd parity perturbation. Note that we omitted to write indices , m in these unknown functions.
Similar to the case of the Klein–Gordon equation, h μ ν m satisfies relations
L ξ ( t ) h μ ν m = i ω h μ ν m ,
L ξ ( 3 ) h μ ν m = i m h μ ν m ,
L 2 h μ ν m = ( + 1 ) h μ ν m .
We should note that if a symmetric tensor h μ ν satisfies the above equations, h μ ν should take the form of Equations (A51) and (A52) from a similar discussion to obtain Equation (A46) for the vector case.
Because the operator X ^ in Equation (A48) depends only on the covariant derivative of the background metric g μ ν ( Sch ) , X ^ and the Lie derivative with respect to the Killing vectors are commute, and the commutation relations [ X ^ , L ξ ( t ) ] = [ X ^ , L ξ ( 3 ) ] = [ X ^ , L 2 ] = 0 hold. Thus, the relations
L ξ ( t ) X ^ h μ ν m = i ω X ^ h μ ν m ,
L ξ ( 3 ) X ^ h μ ν m = i m X ^ h μ ν m ,
L 2 X ^ h μ ν m = ( + 1 ) X ^ h μ ν m ,
hold and X ^ h μ ν m should take the same form of Equations (A51) and (A52) but coefficients are different functions. This implies that we can separately discuss different , m , ω modes and we need to solve ordinary differential equations of r for each , m , ω mode.

Appendix A.4.2. Parity Transformation

In this subsection, we show that we can separately discuss the even and odd parity perturbations. We define the parity transformation operator P which is generated by the coordinate transformation
( θ , ϕ ) ( π θ , π + ϕ ) .
Then, the relations
P Y m = ( 1 ) Y m ,
P V i m = ( 1 ) + 1 V i m ,
hold, where we used P ϵ ^ i j = ϵ ^ i j .35 Using these relations, h μ ν ( ± ) m satisfy
P h μ ν ( + ) m = ( 1 ) h μ ν ( + ) m ,
P h μ ν ( ) m = ( 1 ) + 1 h μ ν ( ) m .
Because the background metric g μ ν ( Sch ) is invariant under the parity transformation and X ^ has the same symmetry as that of g μ ν ( Sch ) , the commutation relation [ P , X ^ ] = 0 holds, and then
P X ^ h μ ν ( + ) m = ( 1 ) X ^ h μ ν ( + ) m ,
holds. We write X ^ h μ ν ( + ) m as
X ^ h μ ν ( + ) m = C + Z μ ν ( + ) + C Z μ ν ( ) ,
where C ± are constants, and Z μ ν ( + ) and Z μ ν ( ) take the same form as Equations (A51) and (A52), respectively. On the other hand, using Equation (A64), we obtain
( 1 ) X ^ h μ ν ( + ) m = P X ^ h μ ν ( + ) m = ( 1 ) C + Z μ ν ( + ) + ( 1 ) + 1 C Z μ ν ( ) ,
and this implies
X ^ h μ ν ( + ) m = C + Z μ ν ( + ) C Z μ ν ( ) .
Equations (A65) and (A67) show C Z μ ν ( ) = 0 , and we conclude that X ^ h μ ν ( + ) m does not contain the odd parity components. In a similar way, we can see that X ^ h μ ν ( ) m does not contain the even parity components. Thus, even and odd parity perturbations can be separately discussed because they are not coupled in the Einstein equations.

Appendix A.4.3. Gauge Transformation

The total line element at O ( ϵ ) is given by
d s 2 = ( g μ ν ( Sch ) + ϵ h μ ν ) d x μ d x ν .
If we introduce a coordinate transformation
x μ x μ ϵ ζ μ ,
the line element changes
d s 2 g μ ν ( Sch ) d x μ d x ν + ϵ 2 ( μ ζ ν ) + h μ ν d x μ d x ν ,
where μ is the covariant derivative with respect to the Schwarzschild metric. Thus, there is an gauge ambiguity due to the term 2 ( μ ζ ν ) at O ( ϵ ) . Writing ζ μ as
ζ μ d x μ = , m e i ω t ( ζ t d t + ζ r d r ) Y m e i ω t ζ S ^ i Y m d x i ( + 1 ) + e i ω t ζ V V i m d x i ,
the perturbed metric transforms
H 0 H 0 2 i ω ζ t f f ζ r ,
H 1 H 1 i ω ζ r f f ζ t + ζ t ,
H 2 H 2 + f f ζ r + 2 ζ r ,
H t Ω H t Ω ( + 1 ) r ζ t i ω r ζ S ,
H r Ω H r Ω ( + 1 ) r ζ r 2 r 2 ζ S + 1 r ζ S ,
K K + f r ζ r + ( + 1 ) 2 r 2 ζ S ,
K 2 K 2 ( + 1 ) r 2 ζ S ,
for the even parity modes, and
h 0 h 0 i ω ζ V ,
h 1 h 1 2 r ζ V + ζ V ,
h Ω h Ω ( + 1 ) 1 2 r 2 ζ V ,
for the odd parity modes, where ′ denotes the radial derivative.
Hereafter, we focus on 2 modes because = 0 , 1 modes do not contain dynamical degrees of freedom. Then, we can choose H t Ω = H r Ω = K 2 = h Ω = 0 gauge, which is called the Regge–Wheeler gauge, by choosing ζ μ as
ζ t = r H t Ω ( + 1 ) i r 2 ω K 2 ( + 1 ) ,
ζ r = r H r Ω ( + 1 ) + r 2 K 2 ( + 1 ) ,
ζ S = r 2 K 2 ( + 1 ) ,
ζ V = 2 r 2 h Ω ( + 1 ) 1 .
Due to the spherical symmetry of the background spacetime, the master equations for each mode only depend on but not m.36 Thus, we only need to consider m = 0 cases. If we need the explicit metric form of m 0 , it can be obtained by acting the ladder operators of spherical harmonics on m = 0 modes.

Appendix A.4.4. Remarks on Complete Gauge Fixing and Gauge Invariant Quantities

In the gauge choice H t Ω = H r Ω = K 2 = h Ω = 0 , the gauge functions ζ μ are algebraically determined in Equations (A82)–(A85). In general, this type of gauge choice is called the complete gauge fixing.37 We should note that the perturbed metric with the complete gauge fixing does not admit a pure gauge.38 This implies that non-trivial h μ ν with the complete gauge fixing always corresponds to physically meaningful perturbations.
Another remark is that we can construct gauge invariant quantities. Let h μ ν be an arbitrary perturbed metric whose form is h μ ν = h ¯ μ ν + ( μ ζ ν ) ( 1 ) with some vector field ζ μ ( 1 ) . We consider a gauge transformation h μ ν h μ ν + 2 ( μ ζ ν ) ( 2 ) with a vector field ζ μ ( 2 ) . We assume that we can completely fix the gauge by choosing ζ μ ( 2 ) . If we substitute the form of ζ μ ( 2 ) for the complete gauge fixing, which are Equations (A82)–(A85) for the Schwarzschild case, to the perturbed metric h μ ν , all components of the perturbed metric do not contain ζ μ ( 1 ) . This is because if ζ μ ( 1 ) is contained, we can construct the perturbed metric h μ ν with pure gauge by setting h ¯ μ ν = 0 , and this contradicts with the complete gauge fixing. This implies that all components of the perturbed metric with the complete gauge fixing corresponds to the gauge invariant quantities. For the Schwarzschild case, the non-trivial components of the right hand side of Equations (A72)–(A78) with Equations (A82)–(A85) are gauge invariant quantities.

Appendix A.4.5. Odd Parity Perturbation

For the odd parity perturbation, we need to consider the perturbed metric with
h μ ν ( ) 0 = 2 e i ω t sin θ Y 0 θ ( h 0 d t + h 1 d r ) .
The ( θ , ϕ ) component of Einstein equations gives39
h 0 = i f ω ( f h 1 ) .
The ( r , ϕ ) component also gives the second order differential equation of h 1
f 2 h 1 f ( 2 r 5 r H ) r 2 h 1 + 2 r 2 6 r r H + 5 r H 2 ( + 1 ) f r 2 r 4 h 1 + ω 2 h 1 = 0 .
Introducing the master variable
Φ ( 2 ) = r r H r 2 h 1 ,
we obtain the master equation for the odd parity perturbation
f d d r f d Φ ( 2 ) d r + ( ω 2 V ( 2 ) ) Φ ( 2 ) = 0 ,
where V ( 2 ) is the effective potential
V ( 2 ) = f ( + 1 ) r 2 3 r H r 3 .
Equation (A90) is the Regge–Wheeler equation [82].40

Appendix A.4.6. Even Parity Perturbation

For the even parity perturbation, we need to consider the perturbed metric with
h μ ν ( + ) 0 = e i ω t Y 0 ( H 0 d t 2 + 2 H 1 d t d r + H 2 d r 2 + 2 K r 2 ( d θ 2 + sin 2 d ϕ 2 ) ) .
The traceless part of the angular components of the Einstein equations shows the relation
H 0 = f 2 H 2 .
The ( t , θ ) component of the Einstein equations gives
H 2 = 2 K f + i ω f ( f H 1 ) .
Combining ( t , r ) , ( r , r ) and ( r , θ ) components of the Einstein equations, we obtain two equations
H 1 = r H ( 6 r H + r ( 2 + 3 λ ) ) 2 r 2 ( 3 r H + r λ ) f + 2 r 2 ω 2 ( 3 r H + r λ ) f H 1
+ 9 r H 2 + 4 r r H ( 2 + λ ) 4 r 2 λ r ( 3 r H + r λ ) f 2 + 4 r 3 ω 2 ( 3 r H + r λ ) f 2 ( i ω K ) , ( i ω K ) = ( 2 + λ ) ( 2 r H + r λ ) 4 r 2 ( 3 r H + r λ ) r ω 2 3 r H + r λ H 1
+ r H ( 6 r H + r ( 4 + λ ) ) 2 r 2 ( 3 r H + r λ ) f 2 r 2 ω 2 ( 3 r H + r λ ) f ( i ω K ) ,
with λ = 2 + 2 . Introducing new variables Φ + ( 2 ) and D +
H 1 = a 1 Φ + ( 2 ) + a 2 D + ,
ω K = a 3 Φ + ( 2 ) + a 4 D + .
with functions a 1 , a 2 , a 3 , a 4
a 1 = i f 2 i r λ r λ + 3 r H ,
a 2 = 2 i f ,
a 3 = λ ( λ + 2 ) r 2 + 3 r r H λ + 6 r H 2 2 r 2 ( r λ + 3 r H ) ,
a 4 = 1 ,
Φ + ( 2 ) and D + satisfy
f d Φ + ( 2 ) d r = D + ,
f d D + d r = ( V + ( 2 ) ω 2 ) Φ + ( 2 ) ,
where V + ( 2 ) is the effective potential
V + ( 2 ) = f λ 2 ( λ + 2 ) r 3 + 3 λ 2 r H r 2 + 9 λ r H 2 r + 9 r H 3 r 3 ( λ r + 3 r H ) 2 .
Thus, we obtain the master equation for the even parity perturbation
f d d r f d Φ + ( 2 ) d r + ( ω 2 V + ( 2 ) ) Φ + ( 2 ) = 0 ,
which is called the Zerilli equation [83].

Appendix A.4.7. Relation between the Regge–Wheeler and Zerilli Equations

While the effective potentials for the Regge–Wheeler and Zerilli equations in Equations (A91), (A90), (A105) and (A106) are different, it is known that these two equations have same spectra (see e.g., [5] for further discussion). In fact, the solutions of these two equations are connected by transformations
Φ ± ( 2 ) = 1 β ω 2 W Φ ( 2 ) + f d Φ ( 2 ) d r ,
where
β = 4 λ 2 ( λ + 1 ) 2 9 r H 2 ,
W = f 3 r H r ( 3 r H + 2 r λ ) 2 λ ( λ + 1 ) 3 r H .
Thus, we only need to analyze the Regge–Wheeler equation in Equations (A91) and (A90) to discuss the spectra of the system.

Appendix A.5. Maxwell Equations

In a similar way, we can also derive the master equations for the Maxwell equations
μ F μ ν = 0 ,
where F μ ν = μ A ν ν A μ . We expand the gauge field A μ as
A μ d x μ = , m A μ ( + ) m + A μ ( ) m d x μ ,
with
A μ ( + ) m d x μ = e i ω t ( a t d t + a r d r ) Y m e i ω t a S ^ i Y m d x i ( + 1 ) ,
A μ ( ) m d x μ = e i ω t a V V i m d x i .
Note that = 0 mode does not contain dynamical degrees of freedom, and then we focus on 1 . Because the Maxwell equations are invariant under the gauge transformation A μ A μ + μ Ψ with a scalar field Ψ , we can choose one of components of A μ ( + ) m as zero. Expanding Ψ as
Ψ = , m e i ω t ψ Y m ,
the gauge field A μ transforms as
a t a t i ω ψ ,
a r a r + r ψ ,
a S a S + ( + 1 ) ψ .
We can choose a S = 0 gauge and this completely fixes the gauge. We note that a V is gauge invariant. Defining two master variables
Φ + ( 1 ) = f a r i ω ,
Φ ( 1 ) = a V ,
after some calculations, we obtain the master equations
f d d r f d Φ ± ( 1 ) d r + ( ω 2 V ( 1 ) ) Φ ± ( 1 ) = 0 ,
with the effective potential
V ( 1 ) = f ( + 1 ) r 2 .

Appendix A.6. Master Equations for Different Spins

The master equations for the Klein-Gordon, Maxwell and Einstein equations are given by Equations (A29), (A30), (A90), (A91), (A105), (A106), (A120) and (A121). Because the Regge–Wheeler and Zerilli equations have same spectra as discussed in Appendix A.4.7, we can discuss the spectra of the master equations with different spins by studying a single equation in the form
f d d r f d ϕ d r + ( ω 2 V s ) ϕ = 0 ,
where V s is the effective potential
V s = f ( + 1 ) r 2 ( 1 s 2 ) r H r 3 ,
and s denotes the spin of the fields, that is, s = 0 , 1 , 2 for the Klein-Gordon, Maxwell and Einstein equations, respectively. The range of is given by s .

Appendix B. Some Basics on Second Order Ordinary Differential Equations

In this appendix, we review fundamental results on second order linear ordinary differential equations in the complex domain. We do not show any proofs on these results. We just refer the reader to [89,90,91] for rigorous treatments.
Let us consider a second order ordinary differential equation:
y ( z ) + p 1 ( z ) y ( z ) + p 2 ( z ) y ( z ) = 0 .
In this article, it is sufficient to assume that p 1 ( z ) and p 2 ( z ) are rational functions. For the Morse potential (8), by setting z = e β x , the differential equations are transformed into this form. See (13). The master Equation (110) is of this form from the beginning. Moreover most master equations including the Teukolsky equation [46] in black hole perturbation theory belong to this class.
The first step to analyze the differential Equation (A124) is to see singularities. We refer to singular points of p 1 ( z ) and p 2 ( z ) as singular points of the differential Equation (A124). In our assumption, all these singular points of (A124) are poles. In linear differential equations, we have the following nice property:
singular points of solutions to ( A 124 ) singular points of ( A 124 ) .
Note that the converse is not true in general. This statement is not true for non-linear differential equations because there appear movable singular points that depend on integration constants.
If both p 1 ( z ) and p 2 ( z ) are analytic at z = z 0 , this point is an ordinary point of (A124). If p 1 ( z ) has at most a simple pole at z = z 0 and p 2 ( z ) has at most a double pole at the same point, then z = z 0 is a regular singular point of (A124). Otherwise, singular points are called irregular singular points. The singularity at z = are discussed by the variable transformation z = 1 / z . It is particularly important to construct formal series solutions in the neighborhood of a singular point.

Appendix B.1. Regular Singular Points

Let us see series solutions at a regular singular point. The method is well-known as the Frobenius method. Without loss of generality, a regular singular point is located at z = 0 . From the condition of regular singularity, p 1 ( z ) and p 2 ( z ) have the following Laurent expansions
p 1 ( z ) = p 1 , 1 z + O ( 1 ) , p 2 ( z ) = p 2 , 2 z 2 + O ( z 1 ) , z 0 .
The indicial equation is then given by
ρ 2 + ( p 1 , 1 1 ) ρ + p 2 , 2 = 0 .
Let ρ i ( i = 1 , 2 ) be two roots of the indicial equation. These roots are called characteristic exponents.
If ρ 1 ρ 2 is not an integer, two independent series solutions are simply constructed by
y i ( z ) = z ρ i k = 0 a i k z k , i = 1 , 2 .
Let γ be a counterclockwise closed path encircling z = 0 but not encircling any other singular points. Then the analytic continuation of the two solutions (A127) along γ generates two solutions y i γ ( z ) and they are represented by
( y 1 γ ( z ) , y 2 γ ( z ) ) = ( y 1 ( z ) , y 2 ( z ) ) e 2 π i ρ 1 0 0 e 2 π i ρ 2 .
The matrix on the right hand side is referred to as a monodromy matrix.
If ρ 1 ρ 2 is an integer, the classification is slightly intricate. In this case, one of the solutions can be constructed in the same way above. If ρ 1 ρ 2 = m where m is a non-negative integer, the solution corresponding to the exponent ρ 1 is given by
y 1 ( z ) = z ρ 1 k = 0 a 1 k z k .
The second solution in general takes the form
y 2 ( z ) = α y 1 ( z ) log z + z ρ 2 k = 0 a 2 k z k ,
where α is a constant. The monodromy matrix of these solutions is given by
( y 1 γ ( z ) , y 2 γ ( z ) ) = ( y 1 ( z ) , y 2 ( z ) ) e 2 π i ρ 1 1 2 π i α 0 1 ,
where we have used e 2 π i ρ 2 = e 2 π i ρ 1 . In most cases, α is non-zero, and the logarithmic term exists in the second solution.41 In particular, α is never zero for m = 0 (or ρ 1 = ρ 2 ). However in exceptional cases, it is accidentally vanishing, and neither solutions have logarithmic singularities. In this exceptional case, the monodromy matrix is proportional to the unit matrix: e 2 π i ρ 1 diag ( 1 , 1 ) . Furthermore, if ρ 1 is an integer, then the monodromy matrix becomes trivial. In this very special case, the singular point is referred to as an apparent singular point.42
The important fact is that all the infinite sums appearing formal series solutions at a regular singular point are necessarily convergent. Thus the formal solutions naturally lead to analytic solutions. This is not the case for series solutions at an irregular singular point.

Appendix B.2. Irregular Singular Points

Constructions of local series solutions at an irregular singular point are much more involved. Discussing them with full generality is beyond the scope of this article. We just see a few examples. Let z = be an irregular singular point of (A124). We need asymptotic expansions to construct local solutions at the irregular singular point z = . For all the examples in this article, formal local solutions at z = have the following asymptotic expansion:
y formal ( z ) = e a z r z b k = 0 c k z r k ( c 0 = 1 ) ,
where r is referred to as the Poincaré rank.43 This asymptotic expansion is generally divergent, and the solution (A132) is purely formal. Note that the Poincaré rank r is easily computed from the differential Equation (A124) [48]. If p 1 ( z ) and p 2 ( z ) in (A124) behaves as
p 1 ( z ) = O ( z K 1 ) , p 2 ( z ) = O ( z K 2 ) , z ,
then the rank r is given by
r = 1 + max K 1 , K 2 2 .
If p 1 ( z ) is identically zero, we have r = 1 + K 2 / 2 .
We first revisit the generalized Laguerre Equation (13). Using the above prescription, we immediately get r = 1 . Then we plug the ansatz (A132) into (13). The leading behavior in z also determines the parameters a and b as
( a , b ) = ( 0 , ν ) , ( 1 , α ν 1 ) .
These provides two independent solutions. For each of them, one can fix the coefficients c k uniquely, and finally obtains the formal series solutions (39).
Next let us see Airy’s differential equation in detail,
y ( z ) z y ( z ) = 0 .
It is easy to check that it has the irregular singular point at z = . The Poincaré rank is 3 / 2 in this case. Plugging the ansatz (A132) into the Airy equation, one obtains
( a , b ) = 2 3 , 1 4 , 2 3 , 1 4 .
For ( a , b ) = ( 2 / 3 , 1 / 4 ) , the coefficients c k satisfy the following recurrence relation:
48 k c k + ( 6 k 1 ) ( 6 k 5 ) c k 1 = 0 , c 0 = 1 .
It is solved by
c k = 1 2 π k ! 3 4 k Γ k + 1 6 Γ k + 5 6 .
Obviously the other case has almost the same structure. Hence the two formal solutions are given by
y 1 formal ( z ) = e 2 3 z 3 / 2 z 1 / 4 k = 0 1 2 π k ! 3 4 k Γ k + 1 6 Γ k + 5 6 1 z 3 k / 2 , y 2 formal ( z ) = e 2 3 z 3 / 2 z 1 / 4 k = 0 1 2 π k ! 3 4 k Γ k + 1 6 Γ k + 5 6 1 z 3 k / 2 .
Since the coefficients factorially divergent c k k ! , the radius of convergence of these series is zero. In the next subsection, we discuss how to treat these formal solutions.

Appendix B.3. Borel Summations and Stokes Phenomena

We want to construct analytic solutions from formal series solutions in the previous subsection. To do so, we use the Borel summation method. Let us consider a formal power series
F formal ( z ) = k = 0 f k z k .
We often encounter that the coefficient f k diverges factorially, f k k ! ( k ). The (generalized) Borel summation is defined by
F p Borel ( z ) : = 0 d ζ e ζ ζ p 1 B p ζ z , p > 0 ,
where
B p ( ζ ) : = k = 0 f k Γ ( k + p ) ζ k .
Usually we set p unity, but it is convenient to keep it freely. If f k diverges factorially, the Borel transformed series (A143) has a finite radius of convergence. We can analytically continue it to the whole complex domain except for its singular points. Therefore the Borel summation method provides us a definite value from a formal divergent series. One can see the Borel summation (A142) is regarded as the Laplace transform. Then B p ( ζ ) is obtained by the inverse Laplace transform.
We apply the Borel summation to the formal series (A140). By setting p = 5 / 6 , we obtain a good analytic expression of the first solution:
y 1 Borel ( z ) = Γ ( 1 / 6 ) 2 2 / 3 π e 2 3 z 3 / 2 z 1 / 4 0 d ζ e ζ ζ 1 / 6 ( 4 + 3 ζ / z 3 / 2 ) 1 / 6 .
To derive this result, we have implicitly assumed that z is real and positive. In this slice, the integral in (A144) is always convergent. It turns out that as long as z > 0 , this solution is exactly related to the well-known Airy function:
y 1 Borel ( z ) = 2 π Ai ( z ) , z > 0 .
In other words, the asymptotic expansion of the Airy function Ai ( z ) is given by the first equation in (A140). If one is interested in the solution in the complex domain, then the Stokes phenomenon appears. We explain it briefly. It is clear that the integrand of (A144) has two singularities at ζ = 0 and at
ζ = 4 3 z 3 / 2 .
The latter causes the Stokes phenomenon. When arg z = ± 2 π / 3 , ± 2 π , ± 10 π / 3 , , the singularity is located on the positive real axis in the ζ -plane, and the integral is not defined as a result. Moreover, along these rays, the integral (A144) has discontinuities. Therefore, y 1 Borel ( z ) is analytic only in each domain of
2 π 3 ( 2 m 1 ) < arg z < 2 π 3 ( 2 m + 1 ) , m Z .
Let us recall that Airy’s differential equation has only the singular point at z = . Hence its solutions must be analytic in the whole complex plane. In fact, the two Airy functions Ai ( z ) and Bi ( z ) are both entire functions. The Borel resummed solution (A144) is understood as a building block to construct globally analytic solutions because it is not analytic in the complex plane. From these facts, we can extend the equality (A145) to the sector arg z < 2 π / 3 . Beyond this sectorial domain, the asymptotic expansion of Ai ( z ) discontinuously changes, and it is well-known as the Stokes phenomenon. Clearly this phenomenon is associated with the discontinuity of y 1 Borel ( z ) . The rays arg z = 2 π / 3 , 2 π , 10 π / 3 , are Stokes lines, and the sectorial domains separated by these lines are Stokes sectors. Outside the Stokes sector, the relation (A145) must be modified due to the Stokes phenomenon. Ai ( z ) receives an additional contribution. We should emphasize that Ai ( z ) itself does not have any discontinuities along the Stokes lines. The Borel resums of the asymptotic series solutions have them. This point is somewhat confusing. The additional contribution outside the Stokes sector is related to the Borel summation of the second solution y 2 formal ( z ) . As the result, the asymptotic series expansion of Ai ( z ) changes discontinuously.
The similar thing happens for the second solution. The Borel resum of the second solution is given by
y 2 Borel ( z ) = Γ ( 1 / 6 ) 2 2 / 3 π e 2 3 z 3 / 2 z 1 / 4 0 d ζ e ζ ζ 1 / 6 ( 4 3 ζ / z 3 / 2 ) 1 / 6 .
However, this is not defined for positive z. To evaluate it, we have to displace z in imaginary directions slightly. The positive real axis is just a Stokes line. There is a discontinuity along this line. It is not hard to find the following exact discontinuity:
y 2 Borel ( z + i 0 ) y 2 Borel ( z i 0 ) = i y 1 Borel ( z ) , z > 0 .
Note that y 1 Borel ( z ) does not have a discontinuity along the positive real axis. Interestingly the combinations y 2 Borel ( z ± i 0 ) ± i 2 y 1 Borel ( z ) have no discontinuity along the positive real axis. It turns out that we have the following identities:
y 2 Borel ( z + i 0 ) + i 2 y 1 Borel ( z ) = y 2 Borel ( z i 0 ) i 2 y 1 Borel ( z ) = π Bi ( z ) , z > 0 .
The Airy function Bi ( z ) has the Stokes lines along arg z = 0 , 4 π / 3 , 8 π / 3 , . The relation (A150) is understood as the Stoke phenomenon of Bi ( z ) along the Stokes line arg z = 0 . We show the branch cut structure of y 1 Borel ( z ) and y 2 Borel ( z ) (or equivalently the Stokes lines of Ai ( z ) and Bi ( z ) ) in Figure A1.
Figure A1. Branch cuts of y 1 Borel ( z ) (Left) and y 2 Borel ( z ) (Right) in the Airy equation. These correspond to the Stokes lines of Ai ( z ) and Bi ( z ) respectively. We start with the sheet 0, and cross the cuts in two opposite ways. The colored shaded domains are the Stokes sectors, and the dashed lines represent virtual cuts on other sheets.
Figure A1. Branch cuts of y 1 Borel ( z ) (Left) and y 2 Borel ( z ) (Right) in the Airy equation. These correspond to the Stokes lines of Ai ( z ) and Bi ( z ) respectively. We start with the sheet 0, and cross the cuts in two opposite ways. The colored shaded domains are the Stokes sectors, and the dashed lines represent virtual cuts on other sheets.
Universe 07 00476 g0a1

Appendix C. Padé Approximants

Padé approximants are rational functions constructed from a given power series. They are a very powerful and practical tool in numerics because they are not merely approximate functions but also extrapolating functions of the original power series beyond its radius of convergence. They often work even if the radius of convergence is zero. Moreover they have important analytic information on the original function. Padé approximants are discussed in great detail in [92]. We introduce very basic aspects on Padé approximants.
Let us define them more precisely. For a given formal power series
F formal ( z ) = k = 0 f k z k ,
its Padé approximants of degree [ M / N ] are defined by
F [ M / N ] ( z ) : = a 0 + a 1 z + + a M z M 1 + b 1 z + + b N z N ,
where the M + N + 1 coefficients { a m ; b n } are determined uniquely by the condition
F formal ( z ) F [ M / N ] ( z ) = O ( z M + N + 1 ) , z 0 .
Therefore to obtain F [ M / N ] ( z ) , we need the formal power series up to z M + N . We do not discuss an algorithm to determine these coefficients. It is easily done by recent symbolic computational systems in practice. We refer to (A152) as the [ M / N ] Padé approximant. Padé approximants with M = N are called diagonal Padé approximants. Obviously, for N = 0 , the [ M / 0 ] Padé approximant is just the (truncated) original formal power series itself.
Theoretical aspects of Padé approximants seem complicated [92]. For instance, it is unclear for us what is the best choice of M and N in general when M + N is fixed. Empirically, the diagonal Padé approximant seems to be the best. In this appendix, we rather focus on its practical aspects by seeing a few examples.

Appendix C.1. Convergent Series

Let us consider a simple function
F ( z ) = z + 1 z 2 + 1 = 1 + z 2 5 z 2 8 3 z 3 16 + 51 z 4 128 + 47 z 5 256 369 z 6 1024 + O ( z 7 ) .
We compare this exact function with its Padé approximants as well as its Taylor series expansion. Clearly the complex function F ( z ) has singularities at z = 1 , ± i , and the radius of convergence of its Taylor series is unity. Therefore the Taylor series approximation does not work for | z | > 1 . Even in this region its Padé approximants still work. Figure A2 shows graphs of F ( z ) , F [ 12 / 0 ] ( z ) , F [ 6 / 6 ] ( z ) , F [ 7 / 5 ] ( z ) and F [ 5 / 7 ] ( z ) , obtained by the same truncated sum up to z 12 . In this example, the diagonal [ 6 / 6 ] Padé approximant looks the “best” rational approximate function. Of course, the situations may change depending on M + N . It is not obvious for us whether the diagonal Padé approximant is always “best” or not. However, we can definitely say that the Padé approximants are better than the truncated power series because the latter cannot be used as an approximation beyond its radius of convergence while the former can.
Figure A2. The Padé approximants of (A154) have good extrapolations to z > 1 where the Taylor series approximation necessarily breaks down.
Figure A2. The Padé approximants of (A154) have good extrapolations to z > 1 where the Taylor series approximation necessarily breaks down.
Universe 07 00476 g0a2
Padé approximants also have information on the singularity structure of the original function [92]. At first glance, since all the singularities of Padé approximants are poles, it seems that they cannot describe other types of singularities such as branch points nor essential singularities. However, Padé approximants tell us even about these singularities. Let us see it for the example (A154). In the left panel of Figure A3, we show zeros and poles of the [ 80 / 80 ] Padé approximant of (A154). Clearly we can see the branch cut structure as a cluster of them. Interestingly the Padé approximant automatically chooses the direction of the branch cuts. They are stretched in the opposite direction to the origin.
We illustrate other examples:
G ( z ) = 1 z 2 + 1 exp z z + 1 = 1 + z z 2 z 3 3 + 2 z 4 3 + 2 z 5 15 14 z 6 45 + O ( z 7 ) , H ( z ) = 1 z 2 + 1 k = 0 z 2 k 1 = 1 + z z 2 2 + z 3 2 + 3 z 4 8 z 5 8 5 z 6 16 + O ( z 7 ) .
The function G ( z ) has an essential singularity at z = 1 , and H ( z ) has a natural boundary on | z | = 1 . These singularities are not obvious from the power series expansions at all. We can decipher them by the Padé approximants. We show the zeros and poles of G [ 80 / 80 ] ( z ) and H [ 80 / 80 ] ( z ) in the middle and right panels in Figure A3, respectively. It is interesting to observe that the branch cut structure in H ( z ) outside the natural boundary disappears in the Padé approximant. For more detail on singularities, see Section 2.2 in [92].
Figure A3. Padé approximants tell us singularity structures of complex functions. We show the distribution of zeros (light blue) and poles (light red) for F [ 80 / 80 ] ( z ) (left), G [ 80 / 80 ] ( z ) (middle) and H [ 80 / 80 ] ( z ) (right). The black dashed line is the convergence circle.
Figure A3. Padé approximants tell us singularity structures of complex functions. We show the distribution of zeros (light blue) and poles (light red) for F [ 80 / 80 ] ( z ) (left), G [ 80 / 80 ] ( z ) (middle) and H [ 80 / 80 ] ( z ) (right). The black dashed line is the convergence circle.
Universe 07 00476 g0a3

Appendix C.2. Divergent Series

Padé approximants for formal asymptotic series are involved but interesting. We show an illuminating example known as the Stieltjes series:
S formal ( z ) = k = 0 ( 1 ) k k ! z k .
Obviously it is a divergent series. Its Borel summation is exactly given by
S Borel ( z ) = 0 d ζ e ζ 1 + z ζ ,
where we have taken p = 1 in (A142) as usual. We compute the Padé approximant of the Stieltjes series (A156). For comparison, we also evaluate the truncated sum
S m ( z ) = k = 0 m ( 1 ) k k ! z k .
We show numerical values of S 2 m ( 1 / 10 ) , S [ m / m ] ( 1 / 10 ) and S [ m / m ] ( e 11 π i / 12 / 10 ) for 1 m 10 in Table A1. The Padé approximant rapidly converges to the Borel sum for z > 0 while the truncated sum has an optimal order. However, if z is close to the negative real axis, the convergence gets quite worse. This is because the Borel sum (A157) has a discontinuity along z < 0 . In fact, the Padé approximant implies the cut structure as shown in the left panel of Figure A4. This is perfectly consistent with the expectation from the Borel analysis. Finally, we show the convergence behavior of the Padé approximants in the right panel of Figure A4.
Table A1. Truncated sums and Padé approximants of the Stieltjes series are compared with the Borel sum. For z > 0 , the Padé approximant quickly converges to the numerical value of the Borel sum, but it is not the case if z is near the branch cut z < 0 of the Borel sum.
Table A1. Truncated sums and Padé approximants of the Stieltjes series are compared with the Borel sum. For z > 0 , the Padé approximant quickly converges to the numerical value of the Borel sum, but it is not the case if z is near the branch cut z < 0 of the Borel sum.
m S 2 m ( 1 / 10 ) S [ m / m ] ( 1 / 10 ) S [ m / m ] ( e 11 π i / 12 / 10 )
1 0.920000000000000 0.916666666666667 1.1171803947 0.0395971996 i
2 0.916400000000000 0.915662650602410 1.1194368945 0.0441279793 i
3 0.915920000000000 0.915634674922601 1.1185374181 0.0449593628 i
4 0.915819200000000 0.915633423180593 1.1186500893 0.0442145546 i
5 0.915819200000000 0.915633346051126 1.1185456350 0.0446695805 i
6 0.915899033600000 0.915633340032624 1.1186923372 0.0444312796 i
7 0.916148114432000 0.915633339468114 1.1185170114 0.0445036358 i
8 0.916932719052800 0.915633339406670 1.1186359811 0.0445494854 i
9 0.919778218477568 0.915633339399102 1.1186152844 0.0444638222 i
10 0.931942728518451 0.915633339398066 1.1185679982 0.0445073821 i
S Borel ( z ) 0.915633339397881 0.915633339397881 1.1186001330 0.0445042079 i
Figure A4. (Left) Zeros and poles of the Padé approximant S [ 80 / 80 ] ( z ) of the Stieltjes series (A156). It implies the branch cut on the negative real axis. This is consistent with the discontinuity of the Borel sum (A157). This branch cut is related to the Stokes phenomenon of the Borel summation. (Right) Convergence of the Padé approximants S [ m / m ] ( e i θ / 10 ) . As θ gets close to π , the convergence gets slower.
Figure A4. (Left) Zeros and poles of the Padé approximant S [ 80 / 80 ] ( z ) of the Stieltjes series (A156). It implies the branch cut on the negative real axis. This is consistent with the discontinuity of the Borel sum (A157). This branch cut is related to the Stokes phenomenon of the Borel summation. (Right) Convergence of the Padé approximants S [ m / m ] ( e i θ / 10 ) . As θ gets close to π , the convergence gets slower.
Universe 07 00476 g0a4
In summary, if one uses Padé approximants to divergent series, one has to see the branch cut structure of its Borel summation. If a parameter region one is interested in is far from these cuts, one can safely use the Padé approximants. If the parameter is close to the cuts, one needs a modification. A good resolution is to combine the Borel summation and the Padé approximants, which is sometimes referred to as the Borel–Padé summation. In this method, we first consider the Borel transform (A143), and compute its Padé approximant. A big difference is that the infinite series (A143) is convergent. We can evaluate the Borel summation (A142) by replacing B p ( ζ / z ) by its Padé approximant.

Notes

1
We note that there is also a branch cut contribution, which corresponds to the power law tails. We do not discuss this contribution in the present article.
2
In this section, we use ˜ for resonant state problems to distinguish them from bound state problems.
3
We have the freedom to set = i ˜ , and it leads to another branch of the resonance. Since these two branches are symmetric, it is sufficient to consider either of them.
4
Concerning boundary conditions, one can see the following observation. If the wave function behaves as e E x / in x , then it satisfies the bound state condition in x . After the analytic continuation, the wave function behaves as e i E ˜ x / ˜ , and satisfies the resonant boundary condition.
5
Strictly, there is a subtlety on the number of the allowed bound states. This point is discussed in the next section. This relation is also based on an assumption that the potential V ( x ) has bound states. In the case of the cubic potential for example, both V ( x ) and V ˜ ( x ) have only the resonant states, and the relation (7) should be modified for these eigenvalues.
6
In fact, this state satisfies the boundary condition such that there is no contamination of the decaying mode in | q | . That is, the solution purely grows in | q | .
7
Moreover the infinite number of the resonant states is “doubled” by the other analytic continuation = i ˜ . This corresponds to including n = 1 , 2 , 3 , .
8
We are not sure whether this is strictly proved or not. As far as we checked for many models, it indeed holds.
9
The function N Milne ( E ) in Equation (33) counts the number of nodes, that is, the number of zeros, of the wave function Equation (30). According to the nodal theorem [27,28,29], it corresponds to the number of the bound states with energy less than or equal to E.
10
However, up to the optimal order, the truncated sum gives a very good approximate value. This is why perturbation theory in physics is a successful approximation method
11
If one puts the ansatz ψ ( x ) = exp [ i x p ^ ( x ) d x ] , one obtains the Riccati equation for p ^ ( x ) . It is solved order by order in . The odd order part in the perturbative solution can be expressed by a derivative of the even order part, and it finally leads to (59). Hence p ( x ) has only the even order corrections as in (61). See [33,34] for a rigorous proof
12
We notice that this equation is formally equivalent to Milne’s non-linear Equation (28) if identifying p ( x ) = i / w E ( x ) 2 .
13
In the phase-integral formalism in [33], there is a freedom to choose a “base function” that depends on situations. We simply choose it as Q ( x ) itself. A change of the base function might improve the approximation. This technical issue is not a purpose of this review. See [33].
14
A turning point a is defined by Q ( a ) = E V ( a ) = 0
15
Since p 0 ( x ) has a square root branch cut, it comes back to the same value when going around the two turning points. However p 0 ( x ) changes its sign.
16
The case of n < 0 corresponds to another analytic continuation = i ˜ in our framework.
17
At this moment, we assume that a ˜ and b ˜ are real numbers such that a ˜ < b ˜ , but for the resonant eigenvalues these are actually complex. We easily perform an analytic continuation in this case. If we cannot evaluate the integral analytically, things are more delicate.
18
This transformation is not always necessary. We do it in order to map the problem into a known differential equation.
19
Note that the = 0 mode exists only for s = 0 .
20
For instance, the branch cut of z in Mathematica is along the negative real axis. This is not useful in evaluating the contour complex integral (134) numerically.
21
In the original variable E, the differential equation is also hypergeometric-type, but two of the solutions near E = 0 are expressed in terms of so-called Meijer’s G-functions
22
For Ω = 4 2 i / 27 , we get the complex conjugate to the right hand side in (168). Because of the physical requirement of the QNMs, we have to choose a branch of the square root of ( M ω n ) 2 so that the imaginary part of M ω n must be negative. Therefore the two possible values of Ω finally lead to two branches of M ω n whose real parts have opposite signs to each other.
23
In fact, the idea in [55] is very similar to the uniform WKB approach. The continuity condition in [55] just corresponds to the regularity condition in the uniform WKB method
24
These two approaches turn out to be interrelated due to the Alday–Gaiotto–Tachikawa relation [81]. However, the relation is quite non-trivial. See the introductory section in [74]
25
This conclusion holds even if we do not use the global spherical symmetry of the spacetime. Because the only non-trivial operator which commutes with all S O ( 3 ) Killing vectors is the Casimir operator L 2 , the derivatives with respect to the angular coordinates in □ should be proportional to L 2 . Thus, assuming Φ = e i ω t Φ ˜ ( r ) Y ( θ , ϕ ) where the function Y satisfies L 2 Y = λ Y with a separation constant λ , Φ should take the form of e i ω t Y F ( r ) .
26
We can easily calculate the explicit form of Φ m by using Mathematica.
27
Generalization to compact Riemannian (Einstein) manifolds can be seen in [86].
28
The homogeneous equation ^ i ^ i S = 0 has a constant solution, and this gives an ambiguity of the solution for = 0 mode. However, it does not affect ^ i S .
29
While one may think that there might exitst a tensor part which satisfies t i j t with trace free and divergence free conditions t i t i = ^ i t i j t = 0 , such a term does not exist for S 2 . If we assume that a tensor part t i j t exists, we can consider a perturbed metric γ ¯ i j = γ i j + ϵ t i j t , where γ i j is the metric of a unit two sphere. Using the conditions t i t i = ^ i t i j t = 0 , we can show that the Ricci scalar of γ ¯ i j becomes 2 + O ( ϵ 2 ) . This implies that γ ¯ i j is the metric of a unit two sphere at this order, and the metric can be transformed into the same form as Equation (A31) by a gauge transformation at O ( ϵ ) . Thus, t i j t should be t i j t = 2 ^ ( i ζ j ) with some vector field ζ j . If we write ζ i = ^ i S + V i with ^ i V i = 0 , the conditions t i T i = ^ i t i j T = 0 implies ^ i S = V i = 0 , and then we conclude t i j t = 0 .
30
While there are homogeneous solutions in t T for = 0 and 1 modes, these do not affect ( ^ i ^ j ( 1 / 2 ) γ i j ^ ) t T .
31
The proof of the completeness can be seen in [84].
32
The second equation can also be written as ^ V i m = ( ( + 1 ) + 1 ) V i m .
33
If we use the differential form, ϵ ^ i j becomes ϵ ^ = sin θ d θ d ϕ .
34
See [87] for the case with source terms.
35
We define the parity transformation of the tensor fields P W μ 1 μ 2 μ n as the components of the tensor after performing the coordinate transformation Equation (A59) to the whole tensor W μ 1 μ 2 μ n d x μ 1 d x μ 2 d x μ n which includes the coordinate basis, where ⊗ denotes the tensor product.
36
Because the ladder operators L ± = i L ξ ( 1 ) ± L ξ ( 2 ) are constructed from the Killing vector on S 2 , Equation (A21) shows that different m modes can be generated by the infinitesimal coordinate transformation associated with the Killing vectors on S 2 .
37
We assume that the Einstein equations can be decouposed into some modes which can be separately duscussed by each mode, for example, even and odd parity modes for the Schwarzschild case, and we focus on one of the modes where the generator of the gauge transformation ζ μ contains n functional degrees of freedom. At each mode, we can completely fix the gauge iff n independent components of h μ ν + 2 ( μ ζ ν ) can be set as zero for any h μ ν by algebraically choosing ζ μ . Then, ζ μ is expressed as linear combinations of n components of h μ ν , which will be set to be zero after the gauge choice, and their derivatives.
38
If a pure gauge is possible, the perturbed metric becomes h μ ν = 2 ( μ ζ ¯ ν ) with some non-trivial vector field ζ ¯ μ . From the assumption, we already completely fixed the gauge and n independent components of h μ ν vanish, where n is a functional degrees of freedom of the gauge transformation in the mode which we focus on. If we focus on such n components of the equation h μ ν = 2 ( μ ζ ¯ ν ) , we can solve them algebraically with respect to ζ ¯ μ . The vector ζ ¯ μ is expressed as linear combinations of four components of h μ ν , which are already set to be zero, and their derivatives. This implies ζ ¯ μ = 0 .
39
Similar to the case of the Klein-Gordon equation, we can compute the components of the linearized Einstein equations by using Mathematica.
40
See [88] for the other possible choices of the master variables and the effective potentials.
41
We note that two solutions with the form of Equation (A127) are linearly dependent when α is non-zero.
42
Note that since we can shift the exponents ρ i to ρ i + s by a transformation y ^ ( z ) = z s y ( z ) , the monodromy matrix e 2 π i ρ 1 diag ( 1 , 1 ) can be always reduced to diag ( 1 , 1 ) by y ^ ( z ) = z ρ 1 y ( z ) . In the new function y ^ ( z ) , the singular point z = 0 is now apparent.
43
The definition of the Poincaré rank seems different in the literature. In this article, we follow [48].

References

  1. Nakamura, T.; Oohara, K.; Kojima, Y. General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes. Prog. Theor. Phys. Suppl. 1987, 90, 1–218. [Google Scholar] [CrossRef]
  2. Kokkotas, K.D.; Schmidt, B.G. Quasi-Normal Modes of Stars and Black Holes. Living Rev. Rel. 1999, 2, 2. [Google Scholar] [CrossRef] [Green Version]
  3. Nollert, H.P. Quasinormal Modes: The Characteristic ‘sound’ of Black Holes and Neutron Stars. Class. Quantum Gravit. 1999, 16, R159–R216. [Google Scholar] [CrossRef]
  4. Ferrari, V.; Gualtieri, L. Quasi-Normal Modes and Gravitational Wave Astronomy. Gen. Relativ. Gravit. 2008, 40, 945–970. [Google Scholar] [CrossRef] [Green Version]
  5. Berti, E.; Cardoso, V.; Starinets, A.O. Quasinormal Modes of Black Holes and Black Branes. Class. Quantum Gravit. 2009, 26, 163001. [Google Scholar] [CrossRef]
  6. Konoplya, R.A.; Zhidenko, A. Quasinormal Modes of Black Holes: From Astrophysics to String Theory. Rev. Mod. Phys. 2011, 83, 793–836. [Google Scholar] [CrossRef] [Green Version]
  7. Chandrasekhar, S. The Mathematical Theory of Black Holes; Clarendon Press: Oxford, UK, 1998. [Google Scholar]
  8. Maggiore, M. Gravitational Waves: Volume 2: Astrophysics and Cosmology; Oxford University Press: Oxford, UK, 2018. [Google Scholar]
  9. Andersson, N. Gravitational-Wave Astronomy: Exploring the Dark Side of the Universe; Oxford University Press: Oxford, UK, 2020. [Google Scholar]
  10. Ferrari, V.; Gualtieri, L.; Pani, P. General Relativity and Its Applications: Black Holes, Compact Stars and Gravitational Waves; CRC Press: Boca Raton, FL, USA, 2020. [Google Scholar]
  11. Leaver, E.W. Spectral Decomposition of the Perturbation Response of the Schwarzschild Geometry. Phys. Rev. D 1986, 34, 384–408. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Nollert, H.P.; Schmidt, B.G. Quasinormal Modes of Schwarzschild Black Holes: Defined and Calculated via Laplace Transformation. Phys. Rev. D 1992, 45, 2617–2627. [Google Scholar] [CrossRef]
  13. Andersson, N. Excitation of Schwarzschild Black-Hole Quasinormal Modes. Phys. Rev. D 1995, 51, 353–363. [Google Scholar] [CrossRef]
  14. Andersson, N. Evolving Test Fields in a Black-Hole Geometry. Phys. Rev. D 1997, 55, 468–479. [Google Scholar] [CrossRef] [Green Version]
  15. Berti, E.; Cardoso, V. Quasinormal Ringing of Kerr Black Holes: The Excitation Factors. Phys. Rev. D 2006, 74, 104020. [Google Scholar] [CrossRef] [Green Version]
  16. Abbott, B.P.; LIGO Scientific and Virgo Collaborations; Abbott, R.; Abbott, T.D.; Abernathy, M.R.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; et al. Tests of General Relativity with GW150914. Phys. Rev. Lett. 2016, 116, 221101. [Google Scholar] [CrossRef]
  17. Giesler, M.; Isi, M.; Scheel, M.A.; Teukolsky, S.A. Black Hole Ringdown: The Importance of Overtones. Phys. Rev. X 2019, 9, 041060. [Google Scholar] [CrossRef] [Green Version]
  18. Barack, L.; Cardoso, V.; Nissanke, S.; Sotiriou, T.P.; Askar, A.; Belczynski, C.; Bertone, G.; Bon, E.; Blas, D.; Brito, R.; et al. Black Holes, Gravitational Waves and Fundamental Physics: A Roadmap. Class. Quantum Gravit. 2019, 36, 143001. [Google Scholar] [CrossRef] [Green Version]
  19. Konishi, K.; Paffuti, G. Quantum Mechanics: A New Introduction; Oxford University Press: Oxford, UK, 2009. [Google Scholar]
  20. Hatsuda, Y. Quasinormal Modes of Black Holes and Borel Summation. Phys. Rev. D 2020, 101, 024008. [Google Scholar] [CrossRef] [Green Version]
  21. Eniceicu, D.S.; Reece, M. Quasinormal Modes of Charged Fields in Reissner-Nordstrom Backgrounds by Borel-Pade Summation of Bender-Wu Series. Phys. Rev. D 2020, 102, 044015. [Google Scholar] [CrossRef]
  22. Ferrari, V.; Mashhoon, B. Oscillations of a Black Hole. Phys. Rev. Lett. 1984, 52, 1361–1364. [Google Scholar] [CrossRef]
  23. Ferrari, V.; Mashhoon, B. New Approach to the Quasinormal Modes of a Black Hole. Phys. Rev. D 1984, 30, 295–304. [Google Scholar] [CrossRef]
  24. Zaslavskii, O.B. Black-Hole Normal Modes and Quantum Anharmonic Oscillator. Phys. Rev. D 1991, 43, 605–608. [Google Scholar] [CrossRef] [PubMed]
  25. Cooper, F.; Khare, A.; Sukhatme, U. Supersymmetry and Quantum Mechanics. Phys. Rep. 1995, 251, 267–385. [Google Scholar] [CrossRef] [Green Version]
  26. Milne, W.E. The Numerical Determination of Characteristic Numbers. Phys. Rev. 1930, 35, 863–867. [Google Scholar] [CrossRef] [Green Version]
  27. Hartman, P. Ordinary Differential Equations, 2nd ed.; SIAM: Philadelphia, PA, USA, 1982. [Google Scholar]
  28. Courant, R.; Hilbert, D. Methods of Mathematical Physics; John Wiley & Sons: Hoboken, NJ, USA, 2008. [Google Scholar]
  29. Messiah, A. Quantum Mechanics; Courier Corporation: Washington, DC, USA, 2014. [Google Scholar]
  30. Mariño, M. Lectures on Non-Perturbative Effects in Large N Gauge Theories, Matrix Models and Strings. Fortsch. Phys. 2014, 62, 455–540. [Google Scholar] [CrossRef]
  31. Bender, C.M.; Wu, T.T. Anharmonic Oscillator. Phys. Rev. 1969, 184, 1231–1260. [Google Scholar] [CrossRef]
  32. Sulejmanpasic, T.; Ünsal, M. Aspects of Perturbation Theory in Quantum Mechanics: The BenderWu Mathematica Package. Comput. Phys. Commun. 2018, 228, 273–289. [Google Scholar] [CrossRef]
  33. Fröman, N.F.; Fröman, P.O. Phase-Integral Method: Allowing Nearlying Transition Points; Springer: New York, NY, USA, 1996. [Google Scholar]
  34. Kawai, T.; Takei, Y. Algebraic Analysis of Singular Perturbation Theory; American Mathematical Society: Providence, RI, USA, 2005. [Google Scholar]
  35. Fröman, N.; Fröman, P.O.; Andersson, N.; Hökback, A. Black-Hole Normal Modes: Phase-Integral Treatment. Phys. Rev. D 1992, 45, 2609–2616. [Google Scholar] [CrossRef] [PubMed]
  36. Landau, L.D.; Lifshitz, E.M. Quantum Mechanics: Non-Relativistic Theory; Elsevier: Amsterdam, The Netherlands, 1981. [Google Scholar]
  37. Voros, A. The Return of the Quartic Oscillator. The Complex WKB Method. Ann. IHP Phys. Théorique 1983, 39, 211–338. [Google Scholar]
  38. Dunham, J.L. The Wentzel-Brillouin-Kramers Method of Solving the Wave Equation. Phys. Rev. 1932, 41, 713–720. [Google Scholar] [CrossRef]
  39. Balian, R.; Parisi, G.; Voros, A. Discrepancies from Asymptotic Series and Their Relation to Complex Classical Trajectories. Phys. Rev. Lett. 1978, 41, 1141–1144. [Google Scholar] [CrossRef]
  40. Delabaere, E.; Pham, F. Resurgent Methods in Semi-Classical Asymptotics. Ann. IHP Phys. Théorique 1999, 71, 1–94. [Google Scholar]
  41. Álvarez, G. Langer–Cherry Derivation of the Multi-Instanton Expansion for the Symmetric Double Well. J. Math. Phys. 2004, 45, 3095–3108. [Google Scholar] [CrossRef]
  42. Dunne, G.V.; Ünsal, M. Uniform WKB, Multi-Instantons, and Resurgent Trans-Series. Phys. Rev. D 2014, 89, 105009. [Google Scholar] [CrossRef] [Green Version]
  43. Mashhoon, B. Quasi-Normal Modes of a Black Hole. In Proceedings of the Third Marcel Grossmann Meeting on General Relativity, Shanghai, China, 30 August–3 September 1982; Volume 18, pp. 599–608. [Google Scholar]
  44. Schutz, B.F.; Will, C.M. Black Hole Normal Modes—A Semianalytic Approach. Astrophys. J. 1985, 291, L33–L36. [Google Scholar] [CrossRef] [Green Version]
  45. Iyer, S.; Will, C.M. Black-Hole Normal Modes: A WKB Approach. I. Foundations and Application of a Higher-Order WKB Analysis of Potential-Barrier Scattering. Phys. Rev. D 1987, 35, 3621–3631. [Google Scholar] [CrossRef] [PubMed]
  46. Teukolsky, S.A. Rotating Black Holes: Separable Wave Equations for Gravitational and Electromagnetic Perturbations. Phys. Rev. Lett. 1972, 29, 1114–1118. [Google Scholar] [CrossRef] [Green Version]
  47. Hatsuda, Y. An Alternative to the Teukolsky Equation. Gen. Relativ. Gravit. 2021, 53, 93. [Google Scholar] [CrossRef]
  48. Ronveaux, A. Heun’s Differential Equations; Oxford University Press: Oxford, NY, USA, 1995. [Google Scholar]
  49. Fiziev, P.P. Exact Solutions of Regge–Wheeler Equation and Quasi-Normal Modes of Compact Objects. Class. Quantum Gravit. 2006, 23, 2447–2468. [Google Scholar] [CrossRef] [Green Version]
  50. Leaver, E.W. An Analytic Representation for the Quasi-Normal Modes of Kerr Black Holes. Proc. R. Soc. Lond. A 1985, 402, 285–298. [Google Scholar]
  51. Leaver, E.W. Quasinormal Modes of Reissner-Nordstrom Black Holes. Phys. Rev. D 1990, 41, 2986–2997. [Google Scholar] [CrossRef] [Green Version]
  52. Onozawa, H.; Mishima, T.; Okamura, T.; Ishihara, H. Quasinormal Modes of Maximally Charged Black Holes. Phys. Rev. D 1996, 53, 7033–7040. [Google Scholar] [CrossRef] [Green Version]
  53. Clemens, C.H. A Scrapbook of Complex Curve Theory; American Mathematical Society: Providence, RI, USA, 2002. [Google Scholar]
  54. Konoplya, R.A. Quasinormal Behavior of the D-Dimensional Schwarzschild Black Hole and the Higher Order WKB Approach. Phys. Rev. D 2003, 68, 024018. [Google Scholar] [CrossRef] [Green Version]
  55. Dolan, S.R.; Ottewill, A.C. On an Expansion Method for Black Hole Quasinormal Modes and Regge Poles. Class. Quantum Gravit. 2009, 26, 225003. [Google Scholar] [CrossRef]
  56. Matyjasek, J.; Opala, M. Quasinormal Modes of Black Holes: The Improved Semianalytic Approach. Phys. Rev. D 2017, 96, 024011. [Google Scholar] [CrossRef] [Green Version]
  57. Konoplya, R.A.; Zhidenko, A.; Zinhailo, A.F. Higher Order WKB Formula for Quasinormal Modes and Grey-Body Factors: Recipes for Quick and Accurate Calculations. Class. Quantum Gravit. 2019, 36, 155002. [Google Scholar] [CrossRef] [Green Version]
  58. Matyjasek, J.; Telecka, M. Quasinormal Modes of Black Holes. II. Padé Summation of the Higher-Order WKB Terms. Phys. Rev. D 2019, 100, 124006. [Google Scholar] [CrossRef] [Green Version]
  59. Cardoso, V.; Kimura, M.; Maselli, A.; Berti, E.; Macedo, C.F.B.; McManus, R. Parametrized Black Hole Quasinormal Ringdown: Decoupled Equations for Nonrotating Black Holes. Phys. Rev. D 2019, 99, 104077. [Google Scholar] [CrossRef] [Green Version]
  60. McManus, R.; Berti, E.; Macedo, C.F.B.; Kimura, M.; Maselli, A.; Cardoso, V. Parametrized Black Hole Quasinormal Ringdown. II. Coupled Equations and Quadratic Corrections for Nonrotating Black Holes. Phys. Rev. D 2019, 100, 044061. [Google Scholar] [CrossRef] [Green Version]
  61. Kimura, M. Note on the Parametrized Black Hole Quasinormal Ringdown Formalism. Phys. Rev. D 2020, 101, 064031. [Google Scholar] [CrossRef] [Green Version]
  62. Hatsuda, Y.; Kimura, M. Semi-Analytic Expressions for Quasinormal Modes of Slowly Rotating Kerr Black Holes. Phys. Rev. D 2020, 102, 044032. [Google Scholar] [CrossRef]
  63. Kodama, H.; Ishibashi, A. A Master Equation for Gravitational Perturbations of Maximally Symmetric Black Holes in Higher Dimensions. Prog. Theor. Phys. 2003, 110, 701–722. [Google Scholar] [CrossRef] [Green Version]
  64. Ishibashi, A.; Kodama, H. Stability of Higher-Dimensional Schwarzschild Black Holes. Prog. Theor. Phys. 2003, 110, 901–919. [Google Scholar] [CrossRef] [Green Version]
  65. Kimura, M. A Simple Test for the Stability of a Black Hole by S -Deformation. Class. Quantum Gravit. 2017, 34, 235007. [Google Scholar] [CrossRef]
  66. Kimura, M.; Tanaka, T. Robustness of the S -Deformation Method for Black Hole Stability Analysis. Class. Quantum Gravit. 2018, 35, 195008. [Google Scholar] [CrossRef] [Green Version]
  67. Amann, H.; Quittner, P. A Nodal Theorem for Coupled Systems of Schrödinger Equations and the Number of Bound States. J. Math. Phys. 1995, 36, 4553–4560. [Google Scholar] [CrossRef]
  68. Pani, P. Advanced Methods in Black-Hole Perturbation Theory. Int. J. Mod. Phys. A 2013, 28, 1340018. [Google Scholar] [CrossRef] [Green Version]
  69. Kimura, M. Stability Analysis of Schwarzschild Black Holes in Dynamical Chern-Simons Gravity. Phys. Rev. D 2018, 98, 024048. [Google Scholar] [CrossRef] [Green Version]
  70. Pierini, L.; Gualtieri, L. Quasinormal Modes of Rotating Black Holes in Einstein-Dilaton Gauss-Bonnet Gravity: The First Order in Rotation. Phys. Rev. D 2021, 103, 124017. [Google Scholar] [CrossRef]
  71. Wagle, P.; Yunes, N.; Silva, H.O. Quasinormal Modes of Slowly-Rotating Black Holes in Dynamical Chern-Simons Gravity. arXiv 2021, arXiv:2103.09913. [Google Scholar]
  72. Srivastava, M.; Chen, Y.; Shankaranarayanan, S. Analytical Computation of Quasinormal Modes of Slowly Rotating Black Holes in Dynamical Chern-Simons Gravity. Phys. Rev. D 2021, 104, 064034. [Google Scholar] [CrossRef]
  73. Cano, P.A.; Fransen, K.; Hertog, T.; Maenaut, S. Gravitational Ringing of Rotating Black Holes in Higher-Derivative Gravity. arXiv 2021, arXiv:2110.11378. [Google Scholar]
  74. Aminov, G.; Grassi, A.; Hatsuda, Y. Black Hole Quasinormal Modes and Seiberg-Witten Theory. arXiv 2020, arXiv:2006.06111. [Google Scholar]
  75. Bianchi, M.; Consoli, D.; Grillo, A.; Morales, J.F. QNMs of Branes, BHs and Fuzzballs from Quantum SW Geometries. arXiv 2021, arXiv:2105.04245. [Google Scholar]
  76. Novaes, F.; da Cunha, B.C. Isomonodromy, Painlevé Transcendents and Scattering off of Black Holes. J. High Energy Phys. 2014, 2014, 132. [Google Scholar] [CrossRef] [Green Version]
  77. da Cunha, B.C.; Novaes, F. Kerr Scattering Coefficients via Isomonodromy. J. High Energy Phys. 2015, 2015, 144. [Google Scholar] [CrossRef] [Green Version]
  78. Carneiro da Cunha, B.; Cavalcante, J.P. Confluent Conformal Blocks and the Teukolsky Master Equation. Phys. Rev. D 2020, 102, 105013. [Google Scholar] [CrossRef]
  79. Bonelli, G.; Iossa, C.; Lichtig, D.P.; Tanzini, A. Exact Solution of Kerr Black Hole Perturbations via CFT(2) and Instanton Counting. Greybody Factor, Quasinormal Modes and Love Numbers. arXiv 2021, arXiv:2105.04483. [Google Scholar]
  80. Carneiro da Cunha, B.; Cavalcante, J.P. Teukolsky Master Equation and Painlevé Transcendents: Numerics and Extremal Limit. Phys. Rev. D 2021, 104, 084051. [Google Scholar] [CrossRef]
  81. Alday, L.F.; Gaiotto, D.; Tachikawa, Y. Liouville Correlation Functions from Four-Dimensional Gauge Theories. Lett. Math. Phys. 2010, 91, 167–197. [Google Scholar] [CrossRef] [Green Version]
  82. Regge, T.; Wheeler, J.A. Stability of a Schwarzschild Singularity. Phys. Rev. 1957, 108, 1063–1069. [Google Scholar] [CrossRef]
  83. Zerilli, F.J. Gravitational Field of a Particle Falling in a Schwarzschild Geometry Analyzed in Tensor Harmonics. Phys. Rev. D 1970, 2, 2141–2160. [Google Scholar] [CrossRef]
  84. Ishibashi, A.; Kodama, H. Chapter 6 Perturbations and Stability of Static Black Holes in Higher Dimensions. Prog. Theor. Phys. Suppl. 2011, 189, 165–209. [Google Scholar] [CrossRef] [Green Version]
  85. Jansen, A.; Rostworowski, A.; Rutkowski, M. Master Equations and Stability of Einstein-Maxwell-Scalar Black Holes. J. High Energy Phys. 2019, 2019, 36. [Google Scholar] [CrossRef] [Green Version]
  86. Ishibashi, A.; Wald, R.M. Dynamics in Non-Globally-Hyperbolic Static Spacetimes: III. Anti-de Sitter Spacetime. Class. Quantum Gravit. 2004, 21, 2981–3013. [Google Scholar] [CrossRef] [Green Version]
  87. Martel, K.; Poisson, E. Gravitational Perturbations of the Schwarzschild Spacetime: A Practical Covariant and Gauge-Invariant Formalism. Phys. Rev. D 2005, 71, 104003. [Google Scholar] [CrossRef] [Green Version]
  88. Lenzi, M.; Sopuerta, C.F. Master Functions and Equations for Perturbations of Vacuum Spherically-Symmetric Spacetimes. Phys. Rev. D 2021, 104, 084053. [Google Scholar] [CrossRef]
  89. Coddington, E.A.; Levinson, N. Theory of Ordinary Differential Equations; McGraw-Hill: New York, NY, USA, 1955. [Google Scholar]
  90. Wasow, W. Asymptotic Expansions for Ordinary Differential Equations; Courier Dover Publications: New York, NY, USA, 2018. [Google Scholar]
  91. Delabaere, E.; Loday-Richaud, M.; Mitschi, C.; Sauzin, D. Divergent Series, Summability and Resurgence I-III; Springer International Publishing: Berlin/Heidelberg, Germany, 2016. [Google Scholar]
  92. Baker, G.A., Jr.; Graves-Morris, P. Pade Approximants: Encyclopedia of Mathematics and It’s Applications; Cambridge University Press: Cambridge, UK, 1996; Volume 59. [Google Scholar]
Figure 1. Scattering process by a potential wall. Resonant states are defined by the absence of the incoming wave. This occurs when 1 / R 0 (or roughly | T | | R | 1 ).
Figure 1. Scattering process by a potential wall. Resonant states are defined by the absence of the incoming wave. This occurs when 1 / R 0 (or roughly | T | | R | 1 ).
Universe 07 00476 g001
Figure 2. In this article, we mainly consider two types of potentials. These appear in the Morse potential (8) (left) and in an effective potential (111) in the linear perturbation of the Schwarzschild spacetime (right). Both cases can be treated in a unified manner.
Figure 2. In this article, we mainly consider two types of potentials. These appear in the Morse potential (8) (left) and in an effective potential (111) in the linear perturbation of the Schwarzschild spacetime (right). Both cases can be treated in a unified manner.
Universe 07 00476 g002
Figure 3. (Left) The number of the allowed bound states in the Morse potential. (Right) The exact spectrum ε n against the parameter g. The solid lines are the physical bound state energies while the dashed lines are on the unphysical branch. For 2 / 3 < g < 2 , the branch for n = 0 is physical but those for n 1 are all unphysical. The total number of the bound states is N bound = 1 , which is consistent with the left panel.
Figure 3. (Left) The number of the allowed bound states in the Morse potential. (Right) The exact spectrum ε n against the parameter g. The solid lines are the physical bound state energies while the dashed lines are on the unphysical branch. For 2 / 3 < g < 2 , the branch for n = 0 is physical but those for n 1 are all unphysical. The total number of the bound states is N bound = 1 , which is consistent with the left panel.
Universe 07 00476 g003
Figure 4. Two counting functions of the Morse potential. The (red) thick line is Milne’s one, and the (blue) thin line is N exact ( E ) .
Figure 4. Two counting functions of the Morse potential. The (red) thick line is Milne’s one, and the (blue) thin line is N exact ( E ) .
Universe 07 00476 g004
Figure 5. The (left) panel shows the divergence of the asymptotic sum S m ( z ) . We set ν = 1 / 3 , α = 3 / 2 and z = 1 , 2 , 3 . The dashed lines represent the values of the Borel sum S Borel ( z ) . The Borel sum gives the “true value” of the formal divergent sum S ( z ) . The (right) panel shows the optimal order of S m ( z ) . Beyond this order, S m ( z ) gets away from the “true value” S Borel ( z ) .
Figure 5. The (left) panel shows the divergence of the asymptotic sum S m ( z ) . We set ν = 1 / 3 , α = 3 / 2 and z = 1 , 2 , 3 . The dashed lines represent the values of the Borel sum S Borel ( z ) . The Borel sum gives the “true value” of the formal divergent sum S ( z ) . The (right) panel shows the optimal order of S m ( z ) . Beyond this order, S m ( z ) gets away from the “true value” S Borel ( z ) .
Universe 07 00476 g005
Figure 6. Two turning points in the WKB method.
Figure 6. Two turning points in the WKB method.
Universe 07 00476 g006
Figure 7. (Left) The shape of the potential (177) is shown. The parameter is fixed as k = 0.7 . Our interest is to know whether this potential admits a bound state or not. (Right) This figure show Milne’s counting function for the black string potential (177). There is the critical value k c 0.876 at which the counting function changes its sign.
Figure 7. (Left) The shape of the potential (177) is shown. The parameter is fixed as k = 0.7 . Our interest is to know whether this potential admits a bound state or not. (Right) This figure show Milne’s counting function for the black string potential (177). There is the critical value k c 0.876 at which the counting function changes its sign.
Universe 07 00476 g007
Table 1. The numerical fundamental QNM frequency for the gravitational perturbation s = 2 with = 2 . The numerical values in the WKB analysis are determined by solving (152) for n = 0 . The values in the perturbative analysis are obtained from (168). The same result as (168) is also obtained by the uniform WKB analysis. These approximate methods get better in the eikonal limit , and thus = 2 is the most severe.
Table 1. The numerical fundamental QNM frequency for the gravitational perturbation s = 2 with = 2 . The numerical values in the WKB analysis are determined by solving (152) for n = 0 . The values in the perturbative analysis are obtained from (168). The same result as (168) is also obtained by the uniform WKB analysis. These approximate methods get better in the eikonal limit , and thus = 2 is the most severe.
OrderWKBPerturbation/Uniform WKB
0 0.4686939 0.09646671 i 0.4714045
2 0.3824606 0.08924751 i 0.3748360 0.1210154 i
4 0.3745148 0.08847747 i 0.3723472 0.0896396 i
6 0.3734578 0.08859956 i 0.3735116 0.0882316 i
8 0.3734641 0.08878269 i 0.3737158 0.0887288 i
10 0.3735665 0.08889202 i 0.3736968 0.0889343 i
12 0.3736290 0.08894092 i 0.3736751 0.0889682 i
[ 6 / 6 ] Padé 0.3736770 0.08896645 i 0.3736670 0.0889626 i
Wronskian 0.3736717 0.08896232 i 0.3736717 0.08896232 i
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hatsuda, Y.; Kimura, M. Spectral Problems for Quasinormal Modes of Black Holes. Universe 2021, 7, 476. https://0-doi-org.brum.beds.ac.uk/10.3390/universe7120476

AMA Style

Hatsuda Y, Kimura M. Spectral Problems for Quasinormal Modes of Black Holes. Universe. 2021; 7(12):476. https://0-doi-org.brum.beds.ac.uk/10.3390/universe7120476

Chicago/Turabian Style

Hatsuda, Yasuyuki, and Masashi Kimura. 2021. "Spectral Problems for Quasinormal Modes of Black Holes" Universe 7, no. 12: 476. https://0-doi-org.brum.beds.ac.uk/10.3390/universe7120476

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop