{'fraction_non_alphanumeric': 0.04086908473773533, 'fraction_numerical': 0.025930457433792478, 'mean_word_length': 4.621186910704382, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quasi Normal Modes of Black Holes and Detection in Ringdown Process

Quasi-normal modes (QNMs) of a black hole (BH) are the eigen modes describing the dissipative oscillation of various fields in that spacetime, which can be intrinsically produced by the linear perturbation theory. With the discovery of the first gravitational waves (GWs) event, GW150914, a new window into the universe has been opened, allowing for the detection of QNMs associated to the ringdown process, which will enable more accurate measurements of the BHs parameters as well as further testing of general relativity. This article discusses the linear perturbation theory of BHs and provides review of several QNMs calculation methods including the newly developed methods. We will also focus on the connection between QNMs and the detection of GWs as well as some recent advancements in this area.


INTRODUCTION
Black holes (BHs); Einstein (1915) and gravitational waves (GWs); Einstein (1916) are proposed in the 20th century as the significant ingredients predicted by general relativity (GR). In the next several decades, people are attempt to find out and explore their existent evidences and characteristics; Askar et al. (2019). Since 2015, several GWs events, such as GW150914; Abbott et al. (2016) and GW170817A; Abbott et al. (2017), from the binary black holes (BBHs) of the stellar mass and the binary Neutron stars (BNSs) respectively, have been detected; GWO (2022); Abbott et al. (2019Abbott et al. ( , 2021a, while the optical observations for the shadow from M87 also provided further indirect evidence for supermassive BHs (SMBHs); Collaboration et al. (2019).
BHs own a natural distinguishable feature -an event horizon. This surface, as a one-way causal boundary, separates the communication of information in a classical level and brings a huge obstacle for us to observe the interior; Schwarzschild (1916). Therefore, the BHs, which are the vacuum solutions can only determined by several parameters without the complicated equations of state like stellar or neutron stars; Kokkotas and Schmidt (1999); Nollert (1999), such as Schwarzschild and Kerr BHs, are Schmidt (1999); Berti et al. (2009); Konoplya and Zhidenko (2011). These articles provide a comprehensive exploration of this topic from the aspects of theory, method, detection, etc. However, with the development of the past decade, there have been significant advancements in all aspects: Initially, QNMs equations in more different spacetime for tensor cases associated to GWs were explored, while the reconstruction of metric was developed. In the meanwhile, further methods for calculating QNMs or related eigen functions have been proposed with several available programs. Furthermore, with the successful detection of GWs events; Abbott et al. (2016), the analyses for the detection of GWs during the ringdown process are now under discussion. We will focus on these advancements and review these aspects as follows: We review the QNMs produced by linear perturbation theory in Sec.2 and the method to calculate the QNMs in Sec.3. At the end of the article in Sec.4, we review the detection advancements based on the detected data from LIGO; Aasi et al. (2015), VIRGO;Acernese et al. (2014) and KAGRA; Akutsu et al. (2021). Without further explication, we'll use the units = c = G = 1.

LINEAR PERTURBATION THEORY AND MASTER EQUATIONS OF QNMS
In general, the perturbations of a BH can result from either an additional field injected into spacetime (such as a particle with the mass m M BH falling into a BH; Davis et al. (1972)) or directly from the perturbing metric (such as the ringdown process at the end of a binary BH merger event). Within general relativity (as well as several other gravity theories), the linear perturbation theory requests us to focus on the first order of perturbation and ignore the reaction to the background. Fields or gravitational radiation will propagate in spacetime in the form of damping oscillations, the characteristics of which are typically governed by a set of radial Schrödinger-like equations in frequency domain with the corresponding angular equations There are two methods to study the linear perturbations in the background spacetime: Firstly, one can parameterize the perturbations as the variation of the coefficients of metric directly and insert them into the Einstein equation or Maxwell equation in the curved spacetime. Or, we can also study them in the form of Newman-Penrose (N-P) formalism; Newman and Penrose (1962) via the N-P equations. Both of the two methods in several cases has been summarizes in the monograph from Chandrasekhar; Chandrasekhar and Thorne (1985).
In this section, we give a general review of the Schrödinger-like equations that govern a BH's quasi-normal modes, specifically for Schwarzschild and Kerr BHs. The master equations governing the propagation of fields or gravitational radiation will be our starting point for discussion.
• Scalar Field in Background Spacetime (Scalar Perturbations). The motion of a massless scalar field Φ in the background spacetime can be obtained from Klein-Gordon equation: where ∇ µ is the covariant derivative. The equation mentioned above may be formally rewritten as: with g the determinant of the background metric g µν .

Frontiers
• Maxwell Field in Background Spacetime (Vector Perturbations). In this case, the Maxwell equations govern the propagation of a massless vector field A µ in background spacetime: which can be rewritten in a more explicit form in the curved background spacetime as: • Linear Gravitational Perturbation in Background Spacetime (Tensor Perturbations). Within a gravity theory, "gravitational perturbation" refers to the perturbation of spacetime itself. For metric perturbations, the metric can be expressed explicitly as: whereg µν is the metric of background spacetime, δg µν is the linear perturbation term, the terms of two order perturbations g 2 µν and higher orders are disregarded due of their little impact in comparison to δg µν . And the governing equations are well known Einstein equations provided by: After applying the separation of variables to a decoupled master equation, the radial Schrödinger-like equations can be obtained in the form of: where ω is the eigen frequency coming from the separation of t under the Fourier transform, and V eff is the effective potential determined by the background metric with in general the asymptotic boundary behaviors at infinity as (for example, we show the effective potential in Schwarzschild case in Figure.1): leading to the asymptotic boundary behaviors of eigen function Ψ at infinity determined by the solutions of: The general solutions near boundaries can be written in the combination of e +iωx and e −iωx describing outgoing and incoming waves respectively. And the different combination of the both solutions usually results in three different directions: • Ψ → e −iωx (incoming) at −∞ and Ψ → e +iωx (outgoing) at +∞. This will result in the most fundamental field of quasi-normal modes (QNMs), which is the subject of this article. In this case the imaginary component of ω is often negative due to the stability; Regge and Wheeler (1957).

Frontiers
• Ψ → e −iωx (incoming) at −∞ and Ψ → e −kx at x → +∞. This will give rise to a new discipline referred to as quasi-bound states (QBS) after considering the mass of fields. k ≡ m 2 p − ω 2 where m p denotes the mass of the massive scalar perturbation. The QBS's often employed to explore superradiant instability Brito et al. (2020), which is not involved in this article.  (2019), which will not be taken in our consideration.
For scalar case, there is just one component equation, which is inherently decoupled. However, it is challenging to obtain a decoupled equation in other cases since a Maxwell field is governed by six coupled component equations (Eq.(4)) whereas gravitational perturbations are governed by ten (Eq. (6)). Actually, avenue to the decoupled equation must take the symmetry of background spacetime corresponding to the gauge-invariant variables into account with expressing the master equation in terms of them; Kodama et al. (2000), or N-P formalism may be also helpful; Newman and Penrose (1962). We will discuss how to address this problem in some situations, properly speaking in the Schwarzschild and Kerr cases. At the end of this section, we also summarize several publications including the QNMs equations for tensor perturbation case.

Perturbations in Schwarzschild Spacetime
With the effect of a scalar field in vacuum, we start our discussion in Schwarzschild spacetime. The static spherical BH solution is the well known Schwarzschild metric given by: where M is the mass of BH. As is mentioned above, Eq.
(2) is an inherently decoupled equation with Φ now the function of (t, r, θ, φ). Similar to the technique taken to solve the hydrogen atom problem in quantum mechanics, by substituting the background metric provided by Eq.(10) into the master equation and applying the separation of variables: the QNMs governing equation can be obtained in time domain: where the effective potential is given by: where ( + 1) is the angular separation constant arising from the separation of angular part and the angular equations are: with γ = diag(1, sin 2 θ) the metric of the unit sphere surface S 2 (θ, φ) and Y m (θ, φ) the scalar spherical harmonics.
Suppose that the solution of the perturbation equation with the time dependence R(t, r) = e iωt R(ω, r) corresponding to the Fourier transformation in the standard procedure of the normal modes analysis; Nollert (1999); Kokkotas and Schmidt (1999), and substituting it into Eq.(7) yields the radial QNMs equations: with r * the tortoise coordinate defined by: The gravitational perturbation in Schwarzschild spacetime is then taken into consideration, beginning with the symmetry of the Schwarzschild geometry corresponding to gauge-invariant variables; Thompson et al. (2017); Nagar and Rezzolla (2005); Martel and Poisson (2005); Sarbach and Tiglio (2001). Because of the static spherical symmetry of the background manifold M 4 (t, r, θ, φ), it can be regard as the product of a Lorentzian 2-dimension manifold M 2 (t, r) and a 2-dimension unit sphere surface manifold S 2 (θ, φ) with the metric γ = diag(1, sin 2 θ) as is mentioned in the scalar case. By taking advantage of this, the metric perturbations h µν ≡ δg µν can be decomposed in multipoles known as odd-parity or even-parity depending on their transformation features under parity. The definitions of both multipoles are as follows: under a parity transformation (θ, φ) → (π − θ, π + φ), the odd (or axial) parity part transforms as (−1) +1 while the even (or polar) parity transforms as (−1) +1 . Thus, the metric perturbations h µν in Eq.(5) can be expressed as: Similar to the scalar case, the following step involves the separation of the variables. To do this, it is necessary to introduce the vector spherical harmonics; Nollert (1999); Dewitt and Dewitt (1973); Edmonds (1996) and the tensor spherical harmonics; Thompson et al. (2017); Nagar and Rezzolla (2005); Thorne (1980); Zerilli (1970a); Mathews (1962); Regge and Wheeler (1957) as the angular dependence. And the perturbation of the odd parity could be described as: while that of the even parity is: Frontiers with h 0 , h 1 , H 0 , H 1 , H 2 and K the parameterization coefficients of the perturbation metric as for the function of r, and the Fourier transformation has been employed. The above formalism is derived from the work of Vishveshwara (1970) by using spherical symmetry resulting to m = 0 under the Regge-Wheeler gauge; Kokkotas and Schmidt (1999); Nollert (1999); Nagar and Rezzolla (2005); Regge and Wheeler (1957); Thompson et al. (2017). After substituting the parameterized metric formalism into the Einstein equation provided in Eq.(6), ten coupled two-order differential equations governing gravitational perturbations will be obtained, with three for odd parity and seven for even; Berti et al. (2009).
In order to obtain a decoupled master equation, specific parameterization coefficient combinations must be introduced. However finding the specific combination is challenging, fortunately that of odd parity was first found in Regge and Wheeler (1957) with some mistakes, and rectified by Edelstein and Vishveshwara (1970). Then the decoupled QNMs equation in the form of Eq.(7) can be obtained: where r * is the defined by Eq.(16) and the effective potential is given by: which is identical to the scalar case in Eq.(13) except of the value of s. The QNMs governing equation above for the odd parity called Regge-Wheeler equation.
Meanwhile there is also a decoupled equation for even parity with the effective potential given by: with Λ = 1 2 ( − 1)( + 2) and the corresponding governing equation for even parity are called Zerilli equation originally derived by Zerilli (1970b,a), with the corrected version can be found in the Appendix A of Sago et al. (2003).
After that, by introducing gauge-invariant variables, which were initially proposed by Moncrief (1974), a set of normative and efficient procedures for gravitational perturbations in Schwarzschild spacetime was constructed from the following research can be found in Thompson et al. (2017); Nagar and Rezzolla (2005); Martel and Poisson (2005); Sarbach and Tiglio (2001); Sengupta (1979, 1980). The research for both of the multipoles above yields a significant property called isospectral first discovered by Chandrasekhar in his book; Chandrasekhar and Thorne (1985) with some discussion can be found in Appendix A of Berti et al. (2009) and the recent research in Jaramillo et al. (2022), indicating that various multipoles may generate the same characteristic spectrum. That implies it is sufficient to analyze either of the situation for simplification.
Along the same avenue as before, the Maxwell field in Schwarzschild spacetime can be considered by expressing the Maxwell equations Eq.(4) into the vector harmonics and decoupled into the QNMs equation in the form of Eq. (7): where r * has the same definition as Eq.(16) and the effective potential is given by: For a concise summary, the QNMs equation in Schwarzschild spacetime is provided by: with r * the tortoise coordinate defined the same as Eq. (16): who maps r from the region (2M, +∞) to the region (−∞, +∞) with 2M the horizon of Schwarzschild, and the effective potentials are: with the value of s corresponding to the perturbation types: Now, we can investigate the boundary behaviors of the eigen function by solving the equations as Eq. (9), and the boundary condition for QNMs are: with e +iωr * and e −iωr * denoting outgoing and incoming waves respectively.
There is a thorough discussion of the Schwarzschild spacetime in Dewitt and Dewitt (1973) with two further reviews; Nollert (1999); Kokkotas and Schmidt (1999). Methods for calculating equations with QNMs are discussed in Sec.3 and the reconstruction of the metric from the eigen function R in Eq.(25) can be found in Berti et al. (2009). The identical result can be obtained by using N-P formalism; Newman and Penrose (1962) (see details in Chandrasekhar (1975Chandrasekhar ( , 1984 and Chandrasekhar and Thorne (1985) with the relationship between both of the multipoles for tensor perturbations can also be found).

Perturbations in Kerr Spacetime
Then, we discuss the perturbations of a static rotating axisymmetric BH as characterized generally by the Kerr solution; Kerr (1963) in terms of the Boyer-Lindquist coordinate; Boyer and Lindquist (1967): where Σ ≡ r 2 + a 2 cos 2 θ, ∆ ≡ r 2 − 2M r + a 2 and a = J M is the parameter describing the rotating property with J the angular momentum.
As with the difficulty of finding specific combinations of metric coefficients in Schwarzschild case, it is challenging to deal with relevant problems by using the metric perturbations methods. The metric in Kerr spacetime is determined by two parameters M and a in contrast with Schwarzschild case only M remaining, which further complicates the problem. Separation of the dependence of t and φ is obviously available due to the symmetry from stationarity and axisymmetry respectively; Teukolsky (2015), meanwhile the discovery of the separability of r and θ in scalar case brought hope for dealing with this problem; Carter (1968).
However, it is also difficult to use metric perturbations methods to obtain a decoupled equation for Kerr case, fortunately the development of another method based on the N-P formalism; Newman and Penrose (1962) has proven its advantages for dealing with this problem in Schwarzschild case; Price (1972); Bardeen and Press (1973), thereby providing a superior method for analyzing such problems in Kerr case.
In N-P formalism, one chooses four normalized orthogonal null vectors l, n, m, m * as the basis of a tetrad with the first two of those real and the remaining two being complex and conjugated with each other. The components of those in Boyer-Lindquist coordinate are given by; Teukolsky (1973); Chandrasekhar and Thorne (1985): Following that, one may express the master equations and field quantities on N-P tetrad. For vector perturbations with the master equations given by Eq.(4), electromagnetic field tensor F µν can be donated in three independent complex scalar quantities φ 0 , φ 1 , φ 2 , where φ 0 and φ 2 describing the perturbations of a Maxwell field are defined as: while the vacuum Maxwell equations Eq.(4) can be denoted in four equations.

Frontiers
Similarly, the equations representing tensor perturbations in kerr spacetime Eq.(6) become 18 N-P equations provided by Ricci identities and 8 complex equations derived from Bianchi identities. Meanwhile, Weyl tensors C µνσλ turn to five N-P quantities Ψ 0 , Ψ 1 , Ψ 2 , Ψ 3 , Ψ 4 , and Ricci tensors are transformed into ten (including Ricci scalar), where Ψ 0 and Ψ 4 describing the incoming and outgoing gravitational radiation respectively are defined as: By applying the separation of variables to the field with spin s: with the specific form of the fields associated to different spin s in Eq.(34) called Teukolsky function shown in Table 1, one may obtain the decoupled equations for r and θ respectively given by: and with K ≡ r 2 + a 2 ω − am, s λ m ≡ s A m + a 2 ω 2 − 2amω and s A m the eigen value determined by the angular part equation Eq.(36) produced from the separation of the dependence of θ.
The above decoupled equations are known as the Teukolsky equations and were first proposed and discussed by Teukolsky (1972Teukolsky ( , 1973; Press and Teukolsky (1973), whose reasoning process can be seen in Teukolsky (2015), where he stressed the significance of a N-P variableρ written in Boyer-Lindquist coordinate as:ρ whose real and imaginary parts represent the divergence and curl of the outgoing principal null respectively. And the derivation process can be found from the Teukolsky's original essays as mentioned above, or from Chandrasekhar and Thorne (1985).
The separation of the dependence of θ leads to the angular part equations Eq.(36) with the eigen function s S m called spin-weighted spheroidal harmonics (SWSH) determined by the value of s, , m and aω.
For aω = 0 and s = 0, it reduces to Schwarschild case with s A m = ( + 1) and s S m becoming the scalar spherical harmonics as defined in Eq.(14). When aω = 0 and s = 0, the eigen functions turn to spin-weighted spherical harmonics; Goldberg et al. (1967) with the eigen value s A m = ( + 1) − s(s + 1). However, there is still no analytical solution for SWSH, therefore the determination of eigen values inevitably becomes a numerical problem; Press and Teukolsky (1973); Leaver (1985); Seidel (1989); Berti et al. (2006b) with the fully asymptotically behavior analysis can be found in Hod (2015).
The radial Teukolsky equations Eq.(35) governing the QNMs in Kerr spacetime does not seem to have the same form as Eq. (7). However, under the transformation given by Detweiler (1977) with the tensor case corrected in the Appendix B of Maggio et al. (2021), those can be transformed into the form of Eq. (7). Along the same approach as in the Schwarzschild case, we can derive the following boundary asymptotic behavior; Teukolsky and Press (1974) at spatial infinity: and near horizon: with k = ω − mω + , ω + = a 2M r + and r + the large root of ∆ = 0 corresponds to the event horizon. In the case of QNMs of the BHs, the boundary condition should be chosen that R m behaves incoming near horizon and outgoing at spatial infinity. Meanwhile, r * is the tortoise coordinate in Kerr spacetime defined as: It is still challenging to reconstruct the metric coefficients, due to the difficulties to derive the QNMs equations in Kerr spacetime through another way, namely the metric perturbation method, who similarly prevents a direct relationship between the metric coefficients and the eigen functions of the QNMs equations. The only remaining option is to attempt to reconstruct the metric using Weyl tensors with Ψ 0 and Ψ 4 associated with the eigen functions of Teukolsky equations Eq.(34) for s = ±2. However, because of the value of spin weight s = ±2, information on = 0 and = 1 associated with the perturbations of mass and angular momentum respectively are lost which must be provided by the rest of the N-P equations. Chandrasekhar attempted the construction and gave a set of methods in Chandrasekhar and Thorne (1985), but it was too complex for application in actual research.
Another method called CCK procedure is based on a key result first proposed by Chrzanowski (1975) and developed by Wald (1978); Stewart (1979); Kegeles and Cohen (1979) where they reconstructed the metric perturbation h µν from a spin-2 scalar Hertz potential with the adoption of radiation gauge; Barack and Ori (2001). The first case of metric reconstruction for the nonvacuum situation is given by Ori (2003) with some other applications can be found in Yunes and González (2006); Sano and Tagoshi (2014); Merlin et al. (2016) and a relatively thorough overview of the procedure can be found in van De Meent (2017); Toomani et al. (2021). Furthermore, a recent research seeks to expand the method to the general Lorentzian gauge in order to address the singularity problem in nonvacuum situations; Dolan et al. (2022). Meanwhile, Loutrel et al. (2021) attempt to develop a new method to reconstruct the first-order metric perturbation just from the solution of the first-order Teukolsky equation, without the requirement for Hertz potentials.
It is worth mentioned that the solutions of the homogeneous QNMs equations will exponentially diverge near infinity for example in Schwarzschild and Kerr cases as shown in Fig.2, which brings more difficulties in calculation. Fortunately, another equivalent form to Teukolsky equation Eq.(35) which is more friendly to numeral calculation was developed in Sasaki and Nakamura (1982).

Other Cases
We summarize the publications where the formalism of the QNMs governing equations for gravitational perturbation can be found.
• The QNMs governing equation of the axial gravitational perturbation in general spherical symmetric spacetime can be found in Zhang et al. (2021), while those of polar parity can be found in Liu et al. (2022).
• The perturbation equations for tensor case was discussed in Kerr-Newman-de Sitter spacetime ; Suzuki et al. (1998).
• Extension to the anti-de Sitter spacetime for Schwarzschild case can be found in Cardoso and Lemos (2001), while those of Kerr can be found in Tattersall (2018) under the slow rotation limit.
• The discussion for the linear tensor perturbations under several modified gravity theories can be found in Moulin et al. (2019).

METHODS FOR CALCULATING QNMS
In the previous section, we reviewed linear perturbation theory and the related QNMs equations. We have turned the problem of the perturbations in the curved spacetime into a set of Schrodinger-like equations Eq. (7) with the corresponding boundary conditions, specifically incoming at the horizon and outgoing at spatial infinity. In this section, we will discuss how to solve these equations and obtain the accurate eigenvalues or QNMs.
This seems to be a straightforward eigen value problem: directly integrate from one boundary to the other and use the shooting method (like Chandrasekhar and Detweiler (1975)) to obtain the appropriate eigen values. However, when one does so, the exponentially diverging asymptotic behaviors on the boundaries cause the numerical error to increase exponentially, which contradicts the requirement of the common shooting method to increase the value of r * (or r) as large as possible to match the asymptotic behavior at infinity, making it difficult to find accurate QNMs directly using numerical integration. In other words, singularities at the horizon and spatial infinity bring the difficulties of direct integration significantly by using shooting method. Fig.2a and Fig.2b depict the boundary behaviors of Schwarzschild case for example to illustrate the asymptotic behaviors of exponential growth at the boundary.
The analysis of the challenge of numerically of locating accurate QNMs producing from asymptotic boundary behavior can be found in Nollert and Schmidt (1992) who also introduced Green function methods to deal with the inhomogeneous equations resulting in a series of studies in nonvacuum case, with a thorough discussion can be found in Poisson et al. (2011).
In fact, the employment of analytical or semi-analytical methods to supplement purely numerical methods may simplify and improve the processing of related problems. We will illustrate the general idea of solving this problem by discussing two well known methods: WKB approximation methods and continue fraction methods. Meanwhile, some other methods are listed at the end of this section.

WKB Approximation Methods
WKB (also known as JWKB) approximation methods were fist proposed by Jeffreys (1925) with a general method of employing approximate solutions to solve linear second order differential equations including the Schrödinger equation, and developed by Wentzel (1926); Kramers (1926); Brillouin (1926) with the treatment of the turning points to address specific problems in quantum mechanics. Meanwhile, the fundamental concepts of the WKB method are usually summarized in almost every quantum mechanics literature; Fröman and Fröman (1965); Hall (2013).
Since its first application in the perturbations problem of a BH; Schutz and Will (1985), this method has undergone constant development, and it remains one of the most effective methods for exploring related problems. For convenience, we rewrite the QNMs governing equation Eq.(7) in another form: where Q(x) = ω 2 −V eff and is a small parameter to track the order of WKB approximation first introduced by Iyer and Will (1987) during his research for third order WKB method. After setting = 1, Eq.(41) returns to its original form Eq. (7).
The WKB approximation retains high precision only in the so-called classically allowed region defined by Q(x) > 0. Considering that Q(x) (or V eff ) is usually unimodal, Q(x) ∼ 0 produces two turning points and divides the whole integration domain into three regions as shown in Fig.3. In regions I and III, the WKB approximation is introduced by assuming the solution in the form of the asymptotic series expansion of as: by substituting the ansatz Eq.(42) into Eq.(41) and equating the same powers of , the specific form of S n can be solved order by order. For example, the fundamental and the first order solutions can be solved as: and with the sign of Eq.(43) determined by the asymptotic behavior taken at both of the boundaries. Under the consideration of only the fundamental solution with Ψ ∼ e S 0 , the boundary behavior of Ψ ∼ e ±iωx corresponds to S 0 ∼ ±iωx. Thus, after introducing the four solutions Ψ I + , Ψ I − , Ψ III + and Ψ III − to denote the corresponding signs in regions I and III respectively with the boundary behaviors in region I (spatial infinity) of Fig.3 as: and in region III (horizon): where the above "in" and "out" represent the waves incident from region I (or region III) to region II and the waves emitted from region II to region I (or region III) respectively (not the incoming and outgoing waves). We then obtain the general solutions in regions I and III given as: And the amplitudes in region I are associated with those in region III through the linear matrix: where M 11 , M 12 , M 21 and M 22 are the coefficients determined by the matching of WKB solutions Eq. (47) in regions I and III with the solution in region II respectively.
The determination of the elements of the matrix in Eq.(48) needs to consider the solution in region II by approximating Q(x) in the form of Taylor series at the peak of Q(x) as: with x 0 the point of the maximum of Q(x), Q 0 = Q(x 0 ) and Q 0 the second derivative with respect to x at the point x = x 0 . And the above Taylor expansion approximation is valid under the assumption that |x − x 0 | is a small value, or to be exact considering the scope of region II given by: with a small value which also gives the validity of the approximation. After that Eq.(41) can be rewritten in the form of parabolic cylinder equation; Bender (1978); Olver et al. (2010): with the substitution as; Iyer and Will (1987): And the general solution of this equation can be denoted as the linear combination of the parabolic cylinder functions as: Under the asymptotic behavior of the parabolic cylinder functions, the solutions become: and where x 1 is the smaller turning point and x 2 is the bigger one as shown in Fig.3. Now that we have calculated the asymptotic solutions close to both of the turning points on both sides, we must match them. Specifically, around the bigger turning point that serves as the dividing point between regions I and II, we must match the coefficients of Eq. (55) For QNMs case of a BH, the condition of the normal modes limits the coefficients in region I with Z I in = 0, while the BHs indicate that there is no wave reflected from the horizon with Z III in = 0. By applying the above conditions to Eq.(57), we obtain the limitation: and ν must be an integer corresponding to the overtone number n. After considering Eq.(53), we obtain the QNMs under the first order of WKB approximation determined by: where V 0 is the peak of the effective potential V eff and the signs of n denote the real part of ω as: n = 0, 1, 2, . . . , Re ω > 0 −1, −2, . . . , Re ω < 0 For a concise summary, to explore such problems using the WKB methods, the following stages are usually taken: • The whole integration domain is divided into several regions by the turning points determined by Q(x) ∼ 0, as the unimodal potential shown in Fig.3.
• By applying WKB approximation Eq.(42) to the QNMs governing equations Eq.(41) in the regions I and III determined by Q(x) > 0, we then obtain general solution in these regions as Eq.(47).
• By approximating Q(x) in region II using Taylor expansion and rewriting the equation into analytical parabolic cylinder equation, we obtain the general solutions in the form of the linear combination of parabolic cylinder functions Eq.(54) with the asymptotic behaviors near the turning points as Eq. (55) and Eq.(56).
• By matching the corresponding coefficients and eliminating the A and B near different turning points, we obtain the matrix in Eq.(57).
• By considering the substitution Eq.(53) and applying the corresponding coefficients according to the specific physical problem such as Z I in = Z III in = 0, we then obtain the first WKB order estimate values of QNMs determined by Eq.(60).
The higher WKB approximation methods lead to the same form of Eq.(57), with the modified expression for R in Eq.(58) still only based on ν. In the meanwhile, the higher order Taylor expansion series result in the different substitution of ν in Eq.(53) that leads to QNMs determined by: where Λ i are the functions of the values of the effective potential and the derivatives (up to the i-th order) at the maximum of the effective potential. The explicit modified terms Λ 2 , Λ 3 of third WKB order approximation methods can be found in Iyer and Will (1987) with the calculation of QNMs; Iyer (1987); Kokkotas and Schutz (1988); Seidel and Iyer (1990). Meanwhile, Λ 4 , Λ 5 and Λ 6 for sixth WKB order approximation can be found in Konoplya (2003Konoplya ( , 2004. And those of the thirteenth WKB order approximation were provided by Matyjasek and Opala (2017), with introducing Padé approximation instead of the Taylor series, leading to more accurate results than those of the sixth WKB order in several cases; Konoplya et al. (2019).

Continued Fraction Methods (Leaver's Methods)
The study of continued fraction in mathematics goes back hundreds of years, however, it was not introduced in the eigen value problems until 1934 by Jaffé (1934) where the bound state of the hydrogen molecule ion; Hylleraas (1931) was studied and he obtained a solution with the proof of convergence, while the same discovery was made by Baber and Hassé (1935) independently. Meanwhile, the early related works were reviewed in Leaver (1986b).
With the observation; Leaver (1986b) that the Teukolsky equations are the subclass of spheroidal wave equations arising during the process of Jaffé (1934); Baber and Hassé (1935), Leaver first introduced the continue fraction into the linear perturbation problems and calculated the QNMs in Schwarzschild and Kerr spacetime in Leaver (1985). After decades of development, this method is one of the most effective ways for estimating QNMs and can provide almost the most accurate value of QNMs. We will illustrate the general thought in Kerr spacetime with the unit c = G = 2M = 1 the same as Leaver (1985): Instead of transforming the equations into the form of Eq.(7) in the tortoise coordinate, we often discuss the Teukolsky equations Eq.(35) directly in r coordinates with the boundary conditions for QNMs of a Kerr BH given as: with σ + = ωr + −am √ 1−4a 2 . Following the Leaver's approach, we assume the expression of the Teukolsky functions being finite at the regular singular points or the boundaries as: By substituting the above ansatz into the radial Teukolsky equations Eq.(35) and equating the coefficients of each orders to zero, the expression coefficients satisfy the following three-term recursion relation: and α r k a r k+1 + β r k a r k + γ r k a r k−1 = 0, k = 1, 2 . . .
where α r k , β r k and γ r k are the recursion coefficients functions of k, ω and a, s, m, A m as for the parameters of the QNMs governing equations Eq.(35), with the specific formalism can be found in Leaver (1985). Then, it turns to the problem of dealing with the three-term recursion relation whose properties explored by Gautschi (1967). Eq.(66) leads to the continued fraction which determines the values of QNMs ω as: where the continued fraction R k can be regarded as the function of ω for given a, s, m and A m with the boundary conditions as k → ∞ and k = 0. The analysis of the convergence of the expansion coefficients as k → ∞ indicates; Gautschi (1967): Meanwhile, boundary condition at k = 0 are given from Eq.(65). By substituting it into the continued fraction for R 0 , we obtained the characteristic equation determining the QNMs as: with the equivalent formalism by inverting an arbitrary number of times k given as: For given a, s, m, A m , and by setting k = k c in a large cutoff value with R kc = 1 due to the boundary condition from Eq.(68), QNMs ω now become the roots of Eq.(69) or Eq.(70) and can be obtained through a numeral method.
However, the determination of A m may be a nontrivial problem as is mentioned in Sec.2.2. Following the same approach, the angular Teukolsky equations Eq.(36) can be solved by supposing the series solution for the angular eigen functions as: with u = cos θ. And it turns to the similar three-term recursion relation as: where the form of the corresponding coefficients α θ k , β θ k and γ θ k can be found in Leaver (1985). And the eigen value A m for given a, s and m can be obtained from the above method.
As the overtone value n increases, the convergence of continued fraction worsens as well Starinets (2002) which leads to the calculation for higher overtone requiring larger cutoff value k c with more computing power. Based on this difficulty, a method applicable to higher overtone QNMs was generalized by expand continued fraction R k in the series of 1 √ k ; Nollert (1993); Zhidenko (2006).
Meanwhile, the continued fraction from the Frobenius series can be found in Konoplya and Zhidenko (2011). And the application of this method in several cases can be found in Leaver (1985); Onozawa (1997); Berti et al. (2004) for Kerr BHs, in Leaver (1990) for Reissner-Nordström BHs, in Berti and Kokkotas (2005) for Kerr-Newman BHs.

Other Methods
Due to the complexity of the gravity theories and the resulting spacetime geometries, in many cases, we must deal with the corresponding perturbation problems case by case. We summarize some publications that employ additional methods and asymptotic formalism.
• Shooting Method and Its Extension. For this kind of eigenvalue problem, the obvious option is to integrate directly and use the shooting method. The most straightforward strategy is to integrate directly from horizon to a large cutoff value with equating the coefficient of the outgoing wave to zero; Press and Teukolsky (1973). The effective extension of this method can be found in Chandrasekhar and Detweiler (1975), where he applied Taylor series expansion at the horizon and spatial infinity respectively, before integrating to a specific intermediate point and matching the both solutions by equating the Wronskian of them to zero, while the details and extension in a matrix formalism can be found in Molina et al. (2010). An interesting example using this method can be found in Mai et al. (2022) where the unstable QNMs in a special gravity theory were found.
• The 'Phase-amplitude' Method This method try to deal directly with singularities that appear with r → +∞ by choosing a specific integral curve in the complex plane to make numerical integration methods possible; Fröman et al. (1992), while Andersson (1992) use this method to calculate QNMs in Schwarzschild case.
• Exact Solutions for Special Potentials Even though it is challenging to obtain exact solutions for the actual BHs QNMs equations, there are always Schrodinger-like equations with special potentials associated to analytic solutions. The researches in Blome and Mashhoon (1984) and Ferrari and Mashhoon (1984) analyzed the connection between QNMs and bound states of the inverted effective potential while used inverted Pöschl-Teller potential to approximate the actual BHs effective potentials in Schwarzschild, Kerr and Reissner-Nordström cases. Some of the other useful potentials can be found in Boonserm and Visser (2011).
• Exact Solutions of Heun Equation The Heun equation is a generalization of the hypergeometric equation while the radial and angular perturbation equations can be rewritten in the form of confluent Heun equations ; Arscott et al. (1995); Baumann et al. (2019), with related researches in Fiziev (2006) for Schwarzschild case and in Borissov and Fiziev (2009);Fiziev (2009) (2007); Cho et al. (2022), thus the applicability to perturbation problems seems evident. The details of this method can be found in Mino et al. (1997); Sasaki and Tagoshi (2003).
The examples for Schwarzschild case can be found in Fujita (2012) with 22PN expansion and for Kerr in Fujita (2015) with 11PN expansion. However, such accuracy is still not enough for GWs detection; Sago et al. (2016).
• Asymptotic Iteration Method (AIM) This method is based on a mathematical theorem resulting in a equivalently condition for the linear homogeneous second order ordinary differential equations, including the QNMs equations. By dealing with the above equivalently condition in numeral methods, the accurate QNMs can be found. The details are reviewed in Ciftci et al. (2003); Cho et al. (2012) with the calculation of QNMs can be found in Mamani et al. (2022) for Schwarzschild and the available Julia package in .
• The Pseudo-spectral Method In this method, the continuous independent variables (radial coordinate for QNMs equations) are replaced by a discrete set of points called the grid, thus the eigen function can be approximated by a series of cardinal functions corresponding to the grid; Jansen (2017). Then, the coefficients of each order of eigen functions can be expanded as the series of ω which results in a matrix governing the eigenvalue problem. The calculation of QNMs can be found in Mamani et al. (2022) for Schwarzschild and the details with available Mathematica package in Jansen (2017), • MST Method This method is based on the formalism developed by Mano, Suzuki and Takasugi;Mano et al. (1996). The homogeneous radial Teukolsky solutions in Eq.(35) are expanded in the series of hypergeometric functions near horizon and Coulomb wave functions at spatial infinity, resulting in the three term recurrence relation for their expansion coefficients respectively, which is similar to Leaver's method; Leaver (1986b). However, the both three term recurrence relations among the expansion coefficients are the same which makes the analytically match possible; Fujita and Tagoshi (2004). The details can be found in Fujita and Tagoshi (2004); Fujita et al. (2009) with the program implemented in BHP (2022) and an example can be found in Piovano et al. (2020).

BLACK HOLE SPECTROSCOPY AND DETECTION ADVANCEMENTS
In the previous sections, we reviewed linear perturbation theory leading to the related QNMs equations in Sec.2 and the calculation of QNMs in Sce.3. Then, we wonder what the perturbed BHs at a great distance may appear like on our detectors, accurately speaking, they are the GWs waveform during the ringdown stage, which leads to the detection of BHs spectroscopy; Berti et al. (2006aBerti et al. ( , 2007a. One of the most common states of a perturbed BH is the remnant of the final ringdown stage of a binary BH merger event, which can be regarded as the perturbation of a rotating black problem described by the Teukolsky equations Eq.(36) and Eq.(35). Meanwhile the charges of a BH have been proven to have no detectable effect on the ringdown waveform; Carullo et al. (2022). Due to the difficulties of reconstructing the metric in terms of Teukolsky functions as is mentioned at the end of Sec.2.2, it is difficult to obtain the GWs waveform by directly applying TT gauge; Flanagan and Hughes (2005). However, the asymptotic behavior of Ψ 4 at infinity is naturally associated with the both polarization modes of the outgoing GWs; Teukolsky (1973): where h + = h θθ and h × = h θφ are the both polarization modes and the dots on the top denote the derivative with respect to time t. With the relation between Ψ 4 and the Teukolsky functions Eq.(34) provided in Table  1, the GWs waveform at infinity can be written in; Berti et al. (2006a): where the overtone index n is now introduced to denote the eigen functions with different QNMs ω, ι is the angle between the angular momentum vector and the line-of-sight vector while ϕ 0 is the phase angle based on source frame; Varma et al. (2014). It should be noted that this formalism of reconstruction into the linear superposition of different modes is incompleteness, but the numerical simulation shows the applicability of this formalism at intermediate times; Berti et al. (2006a); Szpak (2004); Beyer (1999Beyer ( , 2001. Then, by separating the real and imaginary parts of QNMs ω = ω mn + i τ mn and substituting the outgoing boundary behavior given from Eq.(38) as R mn = r 3 Z out mω e −iωr with Z out mω = M A mn e iφ mn , the polarization amplitudes can be obtained as: where A +,× mn and φ +,× mn are the amplitude and the original phase respectively in general regarded as the free parameters or determined by the previous stage; ; Taracchini et al. (2012). Meanwhile, ω mn = 2πf mn is the QNMs' real part with f mn the frequency of the oscillation, and τ mn is damping time given from the values of QNMs; Berti et al. (2009) In general, the ringdown waveform is dominated by the mode with = |m| = 2 for Kerr case, while the other multipoles are subdominant; Berti et al. (2007b); Buonanno et al. (2007). And the QNMs for given ( , m) are sorted by the damping time τ mn , where the fundamental mode noted by n = 0 has the longest damping time with the integer overtone index n > 0 labeling the other modes with shorter damping time.
And the mode with ( , m, n) = (2, 2, 0) or noted as (2, 2, 0) mode is usually the fundamental mode for Kerr case. Meanwhile the measurement of a detector is given by: with the pattern functions given as: F × (θ S , φ S , ψ S ) = 1 2 1 + cos 2 θ S cos 2φ S sin 2ψ S + cos θ S sin 2φ S cos 2ψ S .
Where θ S and φ S denote the polar and azimuth angles of the source in the sky based on detector frame, while ψ S is the azimuth angles of the angular momentum vector based on radiation frame; Varma et al. (2014).
Besides the parameters of the source directly determining the waveform generally including the mass M and the spin parameter a, another critical parameter associated to the detectability and measurability is the signal-to-noise ratio (SNR) ρ defined asFinn (1992); Flanagan and Hughes (1998a,b): where S h (f ) is the noise power spectral density (PSD) or sensitivity curve of different detectors; Aasi et al. In general, SNR is the threshold to determine if a signal has been detected. When the SNR exceeds a certain threshold, such as ρ > 2.5 in Cabero et al. (2020), we consider the signal have been detected. Another method of applying the Bayesian model to examine the likelihood of detecting QNMs will yield the Bayes factor defined as: where B AB > 3.2 denotes "substantial" support for H A over H B , B AB > 10 denotes "strong" support and B AB > 100 is "decisive"; Kass and Raftery (1995); Cabero et al. (2020).
Before the first detected GW event GW150914; Abbott et al. (2016), there have been several studies to predict the range of measurable sources, such as the results of Flanagan and Hughes (1998a): As Earth-based GWs detectors continuously probe GWs, data-driven searches for the existence or the accurate detection of QNMs (or overtone) are taking a new direction. Meanwhile, the sources of some detected GWs events are just within the range predicted by Eq.(85); GWO (2022); Abbott et al. (2019Abbott et al. ( , 2021a, and we summarize some recent advancements in this subject: • Detection of Fundamental Mode and Higher Overtone. Using a waveform model with ringdown process has shown obvious advancements in the parameter estimates of mass M and spin a contrast with that without ringdown process; Baibhav et al. (2018); Isi et al. (2019). However, due to the short damping time of the QNMs, it is difficult to identify if the QNMs included in the data and when they happened; London et al. ( , 2016. Therefore, the standard QNMs tests only take the fundamental (2, 2, 0) mode into consideration. However, the fundamental mode alone is not enough to estimate the accurate values of mass M and spin a because of lack of the information when the ringdown process happens. The research from Giesler et al. (2019) tried to consider a model including overtones up to n = 7 and claimed that the spacetime can be well described as a linear perturbed BH directly after the peak. This means that extending the ringdown process directly after the peak by introducing higher overtones may be feasible, with the great significance for how to connect the ringdown waveform after the previous waveform and has led to a series of researches in Bhagwat et al. (2020) Meanwhile the revisiting for the first detected GW GW150914 brought a contradict problem. As is pointed in Cotesta et al. (2022), both of the researches in Bustillo et al. (2021); Abbott et al. (2021e) provided a weak evidence in favor of "overtone model" with log 10 -Bayes factor ∼ 0.6, contradicting with the research in Isi et al. (2019) who claimed at least one of overtones detected with 3.6σ confidence. As this problem is under further exploration, a startling new point appears that the overtones already detected may be noise-dominated because of the low Bayes factors! And some of the related researches can be found in Cotesta et al. (2022); Isi and Farr (2022).
In general, introducing overtones into waveform model indeed brings better effects. However, how to take them into a waveform model and if the overtones detected are noise-dominated may require further exploration and more accurate detection.
• Detection of Higher Angular Modes. In general, the (2, 2, 0) mode is indeed the dominant mode, while the sub-dominant mode is sometimes not the (2, 2, 1) mode but the modes with > 2 or m > 2 called higher angular modes. And some related publications are Capano et al. (2022) • Detection of Nonlinear QNMs. We have already known that under the inclusion of overtones, the waveform can be well-modelled. However, a binary BH merger is a highly nonlinear system, and we don't know if this nonlinear behavior propagates to infinity and contribute to the waveform on our detectors. In Sec.2, the QNMs are produced under the linear perturbations. When we take the second or higher order terms in Eq.(5) into consideration, the nonlinear QNMs can be calculated with the details in Brizuela et al. (2009) for Schwarzschild case and in Loutrel et al. (2021); Ripley et al. (2021) for Kerr. In the meanwhile, the current detectability research revealed that the detection of nonlinear QNMs requires more precise detectors, which may be available in the next generation; Cheung et al. (2022); Mitman et al. (2022).

DISCUSSION
We review the QNMs produced by linear perturbation theory in Sec.2, where we discuss the difficulties of reconstruction of metric and summarize some publications including the formalism of QNMs equations for tensor perturbations. The method to calculate the QNMs are reviewed in Sec.3 including the newly developed methods in the past few decades. At the end of the article in Sec.4, we review the detection advancements and highlight the current difficulties in detection of overtones.   . This is the schematic diagram of WKB method where two turning points x 1 and x 2 are determined by ω 2 ∼ V eff , which make the whole integration domain divided into three regions.