Spectral Problems for Quasinormal Modes of Black Holes

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\'e approximants, etc. 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.


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 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 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

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: 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,h 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, there exist no bound states. This point is crucially important in 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 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 where R is the reflection coefficient. In the region x ∼ −∞, only the transmitted wave exists: where T is the transmission coefficient. It is more convenient to change the overall normalization as follows: where A = 1/R and B = T/R. The resonant states are defined by the absence of the incoming wave: 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 parameterh. If settingh = i h, 3 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

the eigen-equation (1) is related to (2) with
where E bound n (h) is the bound state energy for the Schrödinger equation (1), and E resonance n ( h) 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).

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. 4 Concerning boundary conditions, one can see the following observation. If the wave function behaves as e − √ −Ex/h in x → ∞, then it satisfies the bound state condition in x → ∞. After the analytic continuation, the wave function behaves as e i √ Ex/ h , 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. Also this relation is 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.
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: 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: Then the Schrödinger equation is rewritten as −g 2 d 2 dq 2 + e 2q − 2e q ψ(q) = εψ(q).
We change the variable z = (2/g)e q . The differential equation (10) becomes We further perform the transform 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: zy (z) + (α + 1 − z)y (z) + νy(z) = 0, where 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 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: This is easily solved, and we finally obtain 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: where (a) k is the Pochhammer symbol. We have used instead of = because the right hand side contains formal divergent series, i.e., 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 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.
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: Therefore the number of the allowed bound states, N bound , is given by 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 The eigenvalues and the eigenfunctions for the inverted potential V are given by 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 states 7 while the number of the bound states is finite. To construct highly excited resonant states, we need unphysical "false bound states."

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.

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 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 h = −i h. This corresponds to including n = −1, −2, −3, . . . . where We assume that Q(x) does not have any singular points for x ∈ R. We consider two independent solutions satisfying conditions: 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 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 One can easily check that w E (x) satisfies the following non-linear differential equation: The two solutions are conversely reconstructed by Now, we would like to construct a solution that satisfies the decaying boundary condition at x = −∞. Let us consider the following function: This is a solution to (23) because the function ψ(x) is written as a superposition of the two solutions ψ 1 (x) and ψ 2 (x).
Also, this function behaves as ψ(x) ∼ C/(xw 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 Let us define 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 satisfies 9 This is a kind of quantization conditions. 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, is also a counting function. We show in Figure 4 the behaviors of N Milne (E) and N exact (E) forh = 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 case ofh = 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 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.

Wronskian method
We also review the well-known Wronskian method. The Wronskian is useful to evaluate connection coefficients among local solutions. Let y p1 (x) and y p2 (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 a1 (x) = C 11 y b1 (x) + C 12 y b2 (x), The connection coefficients C ij are evaluated by Wronskians. Let us assume that y a1 (x) and y b1 (x) satisfy boundary conditions we are interested in. Then these two boundary conditions require This determines the eigenvalues of the boundary value problem. In other words, C 12 = 0 means that the two solutions y a1 (x) and y b1 (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 = ∞: 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): 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 (A143), the Bore sum of S ∞ (z) is analytically given by 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 (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).
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 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 Borel ∞1 (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.

Perturbation theory
In this subsection, we consider perturbative expansions of eigenvalues in the quantum parameterh. Since the Schrödinger equation (1) is the form of singular perturbation inh, 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 = √h q. The Schrödinger equation (1) is then written as It is more convenient to define new variables by The Schrödinger equation is now written as In our analysis, the potential generically has the following form: v(gq) Taking the limit g → 0 with q kept finite, we can expand the potential functions as v 0 (gq) where Ω = v 0 (0) and v m,j = v (j) m (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 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: 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 + 3l, where we use a freedom of N n (g) so that A n n,l = δ l0 .
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 where we set A k n,l = 0 for k > n + 3k. This determines A k n,l for k = n + 3k, n + 3k − 1, . . . , n + 1 uniquely. Next, the energy correction is obtained by (53) Finally, we fix the remaining coefficients by 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 Let us apply it to the Morse potential. We start with the rescaled eigen-equation (10). After changing variables q → √ gq and E = (ε + 1)/(2g), this is rewritten as It turns out that the energy eigenvalue has only the first order perturbative correction: 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 (58) 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 resum it by the Borel summation method (or Padé approximants). This point is important to obtain correct numerical values.

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.

Standard WKB
Let us start with the original WKB method. We use (23). Since the equation takes the form of the singular perturbation inh, we put the ansatz: 11 where p(x) satisfies the non-linear differential equation: 12 11 If one puts the ansatz ψ(x) = exp[ ī h x p(x )dx ], one obtains the Riccati equation for p(x). It is solved order by order inh. 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) We can solve (60) perturbatively inh: At the lowest order, we have two branches of the solution: 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 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: The non-linear equation (60) now becomes where We solve (65) perturbatively, 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 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 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 Therefore the derivative termDU 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 orderh 12 (k = 6), we have 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 onh explicitly. For example, in the next section, we encounter the situation such as Even in this case, we can still use the above formulae. The equation (60) is now written as Then the equation (65) becomes 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: Let us define The solution satisfying the boundary condition at x → −∞ is given by where we ignore irrelevant constants. Using the connection formula in [19], we can continue this solution to the region in x → +∞. The result is where In the region x → +∞, the first term satisfies the boundary condition. We conclude that the bound state boundary condition requires 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: 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 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 or equivalently Two remarks should be mentioned in order. As for usual perturbative series, the left hand side is a formal divergent series inh. 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 complexh-plane. Therefore, in general, the left hand side may have non-perturbative corrections inh. 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 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 where a and b are two turning points, and with some phase shift δ that is not relevant in our analysis. See [19]. As seen in subsection 2.1, the resonant states require 1/R = 0, and semiclassically it yields Therefore we obtain the Bohr-Sommerfeld condition for the resonant energy: 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 criteria 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 16 The case of n < 0 corresponds to another analytic continuationh = −i h in our framework.
It is straightforward to derive the quantum corrected condition: or equivalently 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 These are real as long as −1 ≤ ε < 0. The classically allowed integral is exactly given by The Bohr-Sommerfeld condition turns out to give the exact result in this case: The computation of the resonant states for the inverted potential is almost the same: It would be interesting to check that the quantum corrections really do not contribute to the all-order Bohr-Sommerfeld quantization condition.

Uniform WKB
The uniform WKB method has a little bit 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: where D ν (z) is the parabolic cylinder function that satisfies the Weber equation: This differential equation is essentially same as the Schödinger equation for the harmonic oscillator. The other solution is given by D −ν−1 (iz). Plugging the uniform ansatz into the Schrödinger equation (1), one obtains 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 At the lowest order, we have It is formally solved by 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: Then the solution u 0 (x) has the series representation At the first order, we have 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 Similarly at the second order, the regularity of u 2 at x = 0 as well as u 0 and u 1 determines 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 subsection 2.4. It turns out that both are in agreement as expected. See Eq. (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.

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 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 subsection 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.

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/(2M) rather than the radial variable r itself. The master equation then takes the form where 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 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 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.

Wronskian method
We first apply the Wronskian method to this problem. For this purpose, we use the z-variable. Transforming the wave function by 18 then the new function y(z) satisfies the differential equation where 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 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 Plugging it into the confluent Heun equation, we obtain the following three-term recursion: where 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) The leading asymptotic behavior determines a and b as 32 of 81 The QNM boundary condition at infinity is satisfied by the former case. Therefore we seek the formal solution with the form The coefficients c k again satisfy the following three-term recursion: where 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 subsection 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 (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 (z 0 ) have the same numerical accuracy. 19 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  ], the accuracy estimation is about O(10 −30 ). Therefore we can get 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 way for comparison. For instance, the numerical result (125) shows about 30-digit agreement with Leaver's result as expected.

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, Of course, this is not the unique way to introduceh. Ultimately, one can freely introduce it by hand, and seth = 1 at last. Such a freedom are related to a choice of "base functions" in [33]. Here we choose a very particular base function, which naturally relatesh to the multipole index . It has an advantage that we do not have to introduce an additional freedom of parameters. Then, we obtain where 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 limith → 0. We should note that = 0 is singular in our picture. 20 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 We finally obtain the Schrödinger-like equation for the z-variable: where Note that the domain we are interested in is z > 1, not the whole real line. 20 Note that the = 0 mode exists only for s = 0.
We apply the WKB method to the differential equation (131). This is essentially the same form as (23), but Q depends onh 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 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 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. 21 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: 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 the result, we immediately find an ODE for Ω 0 (E), 21 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.
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, 22 .
The general solution is given by 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, 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, 22 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.

Our requirement for the integration constants is to reproduce this expansion from (138).
With the help of Mathematica, we find the following exact values: 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 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: where the superscript (0) means the leading Bohr-Sommerfeld approximation. This value is compared to the QNM frequency for the s = 1 perturbation, The leading approximation already gives 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 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 where¯ We seek a good differential relation for the integrand Z 1 (z)p 0 (z). We find the following one where D 1 is a second order differential operator, 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, 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 We confirmed this statement up to k = 6. Now we can evaluate the Bohr-Sommerfeld condition with the quantum corrections. We want to solve 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, 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 The more detailed result is shown in Table 1.

Perturbative method
Next, we use perturbation theory reviewed in subsection 2.4. Fortunately we have a natural quantum parameterh 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-h expansion of the QNM frequency. As we have already mentioned, the same result is also obtained by the uniform WKB method in subsection 2.5.
For this purpose, the tortoise variable (112) is useful. We rewrite (128) as where The potential now depends onh, 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 (157) Recall that we are considering the resonant problem for the QNMs. Strictly speaking, to use perturbation theory in subsection 2.4, we need to analytically continue the Planck parameterh → −ih in order to invert the potential. See [20] for detail. However, the result in subsection 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 where g = √h and x = gq, v(gq) The perturbative series of can be computed by the method reviewed in subsection 2.4. The angular frequency Ω is given by 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)) where we regard A k 0,l = 0 for k > 3l and for k < 0. This determines A k 0,l for k = 3l, 3l − 1, . . . , 1 in this order. Then the energy correction is obtained by Recall that we set A 0 0,l = δ l0 . For l = 1, the recurrence relation (162) gives , Then we find 0,1 = 0. For l = 2, we have the following non-zero components: 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 Plugging Ω = −4 √ 2i/27 and s = 2, we finally get 23 Recall thath = 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]. Also, the result can be compared with an expansion in [55] 24 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ω pert s=2, =2,n=0 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 23 For Ω = 4 √ 2i/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. 24 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.
(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].

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 where 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 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 withh = 2/[ ( + 1)]. We can use the general formulation in subsection 2.4. Note that to get 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.

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 subsection 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 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-dimension. 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 where f (r) = 1 − r H /r, and the z-direction is compactified. The master equation reduces to the Schrödinger type equation, where x is the tortoise variable, defined by The integration constant c is fixed as follows.
The potential is explicitly given by 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 . Also 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 where E = ω 2 . Milne's counting function is then obtained by 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 This equation can be solved analytically, There is the critical value k c ≈ 0.876 at which the counting function changes its sign.
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: 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 where the integration over (−∞, x min ) can be computed exactly, Our final expression is thus given by 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.

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 in future developments on 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.
Although we did not review in this article, there are also new perspetives 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]. 25 The ideas are based on the common mathematical structure in the spectral problems. These new perspectives provide us powerful analytic treatments on the spectral problems. However, a nice physical interpretation of the relation is still obscure.
where I, J, K run over 1, 2, 3, the operator L ξ (I) denotes the Lie derivative along ξ µ (I) , and I JK is the totally anti-symmetric symbol with 123 = 1. Defining the Casimir operator L 2 as the relations 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 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 holds. In particular, the Killing equations L ξ g µν = 0 become 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 where L is the separation constant and x i denote the coordinates except for λ, we obtain a relation The solution of this equation can be written by 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.
(A20) (2) , the relations hold for different modes m. Using the commutation relations [L ξ (t) , ] = [L ξ (3) , ] = [L 2 , ] = 0, we can also show that Φ satisfies The regular solution of the above equations can be written by 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. 26 Appendix A.

Explicit calculation
We write the form of the Klein-Gordon field Φ as with where we assumed the time dependence as e −iωt and omitted to write indices , m in Φ (0) . Because Eq. (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 to 27 where the effective potential V (0) is given by Thus, we obtain the master equation We comment that if we consider the massive Klein-Gordon equation ( − m 2 )Φ = 0, the effective potential is changed as In this subsection, i, j denote the tensor indices on S 2 , and we raise or lower tensor indices by γ ij . We denote by∇ i the covariant derivative with respect to γ ij . Any regular vector field v i on S 2 can be uniquely decomposed into the scalar part and the vector part as The proof is simple. Taking the divergence of Eq (A32), we obtain an equation∇ i v i = ∇ i∇ i S. Solving this Poisson equation on S 2 ,∇ i S is uniquely determined. 29 Then, V i is given by In a similar way, any regular symmetric tensor field t ij on S 2 can be uniquely decomposed as 30 To prove this, first, taking the trace of Eq. (A34), we obtain t L = t i i /2. Next, acting an operator∇ i∇j on Eq. (A34), after some calculation, we obtain (A36) 29 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. 30 While one may think that there might exitst a tensor part which satisfies t t ij with trace free and divergence free conditions t t i i =∇ i t t ij = 0, such a term does not exist for S 2 . If we assume that a tensor part t t ij exists, we can consider a perturbed metricγ ij = γ ij + t t ij , where γ ij is the metric of a unit two sphere. Using the conditions t t i i =∇ i t t ij = 0, we can show that the Ricci scalar ofγ ij becomes 2 + O( 2 ). This implies thatγ ij is the metric of a unit two sphere at this order, and the metric can be transformed into the same form as Eq. (A31) by a gauge transformation at O( ). Thus, t t ij should be t t ij = 2∇ (i ζ j) with some vector field ζ j . If we write ζ i =∇ i S + V i with∇ i V i = 0, the conditions t T i i =∇ i t T ij = 0 implies∇ i S = V i = 0, and then we conclude t t ij = 0.
Solving Eq. (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)γ ijˆ )t T is uniquely determined. 31 Finally, taking a divergence of Eq. (A34), we obtain an equation The solution of Eq. (A37) with∇ i t v i = 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 v j) .

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 m i which can express any divergence free regular vector field on S 2 . 32 The vector spherical harmonics V m i satisfies the equations 33 where = 1, 2, · · · , and m = 0, ±1, ±2, · · · , ± . We can explicitly construct V m i by using the spherical harmonics Y m as 31 While there are homogeneous solutions in t T for = 0 and 1 modes, these do not affect (∇ i∇j − (1/2)γ ijˆ )t T . 32 The proof of the completeness can be seen in [84]. 33 The second equation can also be written asˆ V m whereˆ ij is the Levi-Civita tensor of γ ij . 34 Using the vector spherical harmonics, any divergence free vector V i on S 2 can be written by where d m are constants. Then, a vector field v i in Eq. (A32) becomes Next, we consider the case of the vacuum Einstein equations. 35 We denote by h µν the small deviation from the Schwarzschild metric where is a small parameter. The total metric at O( ) becomes (A47) 34 If we use the differential form,ˆ ij becomesˆ = sin θdθ ∧ dφ. 35 See [87] for the case with source terms.
We need to solve the vacuum Einstein equations G µν = 0 or equivalently R µν = 0 at O( ) where ∇ µ denotes the covariant derivative with respect to the Schwarzschild metric, and X denotes the operator s.t.Xh µν gives the linearized Einstein tensor. The perturbed metric h µν can be written by where A, B run t, r, and i, j run θ, φ. We can regard h AB , h Ai , h ij 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 m where h (+) m µν is the even parity perturbation and h We should note that if a symmetric tensor h µν satisfies the above equations, h µν should take the form of Eqs.
hold andXh m µν should take the same form of Eqs. (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 hold, where we used Pˆ ij = −ˆ ij . 36 Using these relations, h Because the background metric g (Sch) µν is invariant under the parity transformation andX has the same symmetry as that of g (Sch) µν , the commutation relation [P,X] = 0 holds, and then holds. We writeXh where C ± are constants, and Z and this impliesX 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.
for the even parity modes, and 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 Due to the spherical symmetry of the background spacetime, the master equations for each mode only depend on but not m. 37 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 Eqs. (A82)-(A85). In general, this type of gauge choice is called the complete gauge fixing. 38 We should note that the perturbed metric with the 37 Because the ladder operators L ± = −iL ξ (1) ± L ξ (2) are constructed from the Killing vector on S 2 , Eq. (A21) shows that different m modes can be generated by the infinitesimal coordinate transformation associated with the Killing vectors on S 2 . 38 We assume that the Einstein equations can be decouposed into some modes which can be separately duscussed by each mode, e.g., 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∇ (µ ζ ν) complete gauge fixing does not admit a pure gauge. 39 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 Eqs. (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 settingh µν = 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 Eqs (A72)-(A78) with Eqs. (A82)-(A85) are gauge invariant quantities.
f d dr f dΦ 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 Eqs. (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 where Thus, we only need to analyze the Regge-Wheeler equation in Eqs. (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 where We expand the gauge field A µ as with 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 the gauge field A µ transforms as 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α 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. 42 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. 43 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.
then the rank r is given by 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 (A133) into (13). The leading behavior in z → ∞ also determines the parameters a and b as 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, 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 (A145) has discontinuities. Therefore, y Borel ∞1 (z) is analytic only in each domain of 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 (A145) 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 (A146) 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 Borel ∞1 (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 (A146) 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 formal ∞2 (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 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:  Figure A1. Branch cuts of y Borel ∞1 (z) (Left) and y Borel ∞2 (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.
Note that y Borel ∞1 (z) does not have a discontinuity along the positive real axis. Interestingly the combinations y Borel ∞2 (z ± i0) ± i 2 y Borel ∞1 (z) have no discontinuity along the positive real axis. It turns out that we have the following identities: The Airy function Bi(z) has the Stokes lines along |arg z| = 0, 4π/3, 8π/3, . . . . The relation (A151) is understood as the Stoke phenomenon of Bi(z) along the Stokes line arg z = 0. We show the branch cut structure of y Borel ∞1 (z) and y Borel ∞2 (z) (or equivalently the Stokes lines of Ai(z) and Bi(z)) in Figure A1.

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 its Padé approximants of degree [M/N] are defined by where the M + N + 1 coefficients {a m ; b n } are determined uniquely by the condition 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 (A153) as the [M/N] Padé approximant. Also 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.  Figure A2. The Padé approximants of (A155) have good extrapolations to z > 1 where the Taylor series approximation necessarily breaks down. 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. 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 (A155). In the left panel of Figure A3, we show zeros and poles of the [80/80] Padé approximant of (A155). 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:

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: Obviously it is a divergent series. Its Borel summation is exactly given by where we have taken p = 1 in (A143) as usual. We compute the Padé approximant of the Stieltjes series (A157). For comparison, we also evaluate the truncated sum  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 (A158) 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.  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 (A144), and compute its Padé approximant. A big difference is that the infinite series (A144) is convergent. We can evaluate the Borel summation (A143) by replacing B p (ζ/z) by its Padé approximant. . It implies the branch cut on the negative real axis. This is consistent with the discontinuity of the Borel sum (A158). 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.