Fundamental Tone and Overtones of Quasinormal Modes in Ringdown Gravitational Waves: A Detailed Study in Black Hole Perturbation

Ringdown gravitational waves of compact object binaries observed by ground-based gravitational-wave detectors encapsulate rich information to understand remnant objects after the merger and to test general relativity in the strong field. In this work, we investigate the ringdown gravitational waves in detail to better understand their property, assuming that the remnant objects are black holes. For this purpose, we perform numerical simulations of post-merger phase of binary black holes by using the black hole perturbation scheme with the initial data given under the close-limit approximation, and generate data of ringdown gravitational waves with smaller numerical errors than that associated with currently available numerical relativity simulations. Based on the analysis of the data, we propose an orthonormalization of the quasinormal mode functions describing the fundamental tone and overtones to model ringdown gravitational waves. Finally, through some demonstrations of the proposed model, we briefly discuss the prospects for ringdown gravitational-wave data analysis including the overtones of quasinormal modes.


Introduction
The remnant of merging black holes (BHs) is a perturbed compact object (a Kerr BH [1], or even something more exotic) characterized by the set of complex frequencies known as quasinormal modes (QNMs), and the gravitational radiation from this remnant is called the ringdown phase of gravitational waves (GWs). Ringdown GWs from binary BH mergers are now routinely detected by ground based GW detectors such as Advanced LIGO and Virgo [2,3]. An accurate measurement of QNMs encoded in the ringdown signal therefore offers various GR tests regarding the compact object, for example, to disclose the remnant property (the ergo region of Kerr geometry [4][5][6], ringdown test of general relativity (GR) [7] etc.), and to verify GR itself in the strong-field regime; we refer readers to Ref. [8] for a review of the BH QNMs, and Ref. [9] for that of the ringdown GWs.
The theoretical investigation of QNMs has a long history since early 1970, including Vishveshwara [10], Press [11], Teukolsky and Press [12], and Chandrasekhar and Detweiler [13], and Detweiler [14]. In 1980s, Detweiler [15] basically completed the analysis of QNMs for Kerr BHs. Leaver [16] gave a standard method to calculate the QNM frequencies very accurately (see also, e.g., Refs. [17][18][19] for more modern techniques the remnant BH spacetime, making use of precise NR waveforms [70][71][72]. This facts implies that the second and higher order QNMs would not be excited enough at least in the ( = 2, m = 2) mode (we note, however, that the second order QNMs due to the self-coupling of the first-order, e.g., ( = 2, m = 2) mode was already identified in the ( = 4, m = 4) mode of the NR waveforms [73]; see also, e.g., Refs. [74][75][76] for earlier works on the second order QNMs). This motivates us to work with the close-limit approximation to the binary BH merges, proposed in the seminal work by Price and Pullin [77], and Abrahams and Prices [78].
The close-limit approximation describes the last stage of a binary BH coalescence as a perturbation of a single (remnant) BH, demonstrating a remarkably good agreement with full NR simulation [79,80]; other recent applications of close-limit approximations are, for example, for the initial data of the NR simulation [81] and for binary BH mergers beyond GR [82]. A particularly relevant aspect of the closed-limit approximation here is that the approximated system can be evolved highly accurately in the framework of the linear BH perturbation theory. A range of methods for numerically evolving the linearized Einstein equation is described in literature, and the numerical accuracy that we can achieve in this work is at the level of 10 −7 in the GW strain data at the time of the peak amplitude and 10 −13 in the late time of simulations (see Figure 2). Therefore, our close-limit waveforms in the BH perturbation approach allow us the more depth study regarding to the QNM fit than using NR waveform with numerically accuracy available for now. This paper is organized as follows. In Section 2, we briefly review the close-limit approximation of the post-Newtonian (PN) metric [83], which we use to prepare the initial data for the BH perturbation calculation. In the calculation presented in this work, for simplicity, we focus on the head-on collision of nonspinning binary BHs. Therefore, the remnant BH can be considered as a (nonspinning) Schwarzschild BH. Our numerical method in the BH perturbation approach is presented in Section 3 where we carefully check the numerical accuracy. In Section 4, we propose a modified QNM fitting formula by introducing an orthonormal set of mode function to analyze the ringdown waveforms. In Section 5, we use the QNM fitting formulae in Section 4 and apply them to the GW data obtained in Section 3. We will find the late-time power-law tail [84] in the fit residuals, and a convergent behavior for the modified QNM fitting formula. In Section 6, we summarize our analysis, and discuss some application. A brief analysis of the dominant ( = 2, m = 2) mode of NR waveforms for nonprecessing binary BHs and the fitting with the QNMs overtones is given in Appendices A and B. Here, we find some problematic behavior in NR waveforms for the cases with a highly spinning remnant BH, and this is discussed in Appendix C.
Throughout this paper, we adopt the negative metric signature (− + ++), geometrized units with c = G = 1, except in Section 2. The coordinates, x α = {ct, r, θ, ϕ}, are used as the Schwarzschild coordinates for Schwarzschild spacetime and the Boyer-Lindquist coordinates for Kerr spacetime.

Post-Newtonian initial conditions in the close limit: head-on collisions
We wish to describe the evolution of nonspinning binary BHs in its last stage, making use of the close-limit approximation starting from a given initial condition. For this purpose, we adopt the close-limit form of the post-Newtonian metric for two nonspinning BHs (or point particles) of masses m A (A = 1, 2) derived by Le Tiec and Blanchet [83]. This close-limit, PN metric is particularly convenient to examine the ringdown phase of asymmetric mass-ratio binaries. To set the stage for our analysis, we review their essential discussion and results here.
We assume that the binary system admits two different (dimensionless) small parameters. The first one is the close-limit parameter [77,78,80,[85][86][87] in which the binary separation r 12 is sufficiently small compared to the distance r from any field point to a reference source point (e.g., the center of mass of the system): The other is the post-Newtonian (PN) parameter [88][89][90][91] where the typical value of the relative orbital velocity v 12 is small enough compared to the speed of light (or that of the separation between two masses r 12 is sufficiently large): Here, M = m 1 + m 2 is the total mass of the binary system, and we shall say that a term relative to O( N PN ) is Nth PN order. We note that the PN expansion is valid when r r 12 / GM/(c 2 r 12 ) that defines the (so named) near-zone.
The crux of this approach is that the metric of the binary system is formally expanded in powers of both in CL and PN (though the close-limit approximation would imply r 12 GM/c 2 i.e., PN 1 [84]; see also Refs. [83,92,93] for an elaboration of this double expansion). This allows us to recast the (late inspiral phase of) PN binary spacetime to a PN vacuum perturbation of the single Schwarzschild spacetime. Specifically, we work with the 2PN near-zone metric g PN µν of two point masses [94][95][96][97], and then further expand it in the power of CL . The resultant metric can be manipulated to be identified with where g Schw.

µν
is the Schwarzschild metric with mass M and h µν is the perturbation that takes the schematic form We note that this expression is truncated at O(G) so that h µν consistently satisfies the linearized Einstein equation to the 2PN order. The general explicit expressions for the coefficients h (n, k) µν are too unwieldy to be presented here; the expressions up to k = 2 in the (slightly altered) Cartesian harmonic coordinates are listed in Eq. (2.12) of Ref. [83], for example.
Because the Schwarzschild spacetime is spherically symmetric, the angular dependence of h µν in Eq. (4) can be conveniently separated to decompose it in a given tensor spherical harmonics. Using a standard Schwarzschild coordinates, this mode decomposition takes a form (e.g., Section 4.2 of Ref. [98]) where Y (i), m µν is a basis of tensor spherical harmonics. A set of 10 time-radial functionsh (i), m further decouples into two subsets: seven for even-parity perturbations while three for odd-parity ones, depending on whether they are "even" or "odd" under the transformation (θ, ϕ) → (π − θ, π + ϕ).

Model problem: head-on collisions of two nonspinning black holes
We now specialize the close-limit, 2PN metric to a head-on collision of two nonspinning BHs. The head-on collision is not an astrophysically relevant scenario, but it provides a useful and simple example to explore the ringdown radiation with a good numerical accuracy in the remaining sections.
We adopt the notations of Le Tiec and Blanchet (Section 3 of Ref. [83]), and work with the Regge-Wheeler basis of tensor Spherical harmonics and in the Regge-Wheeler gauge to simplify the expressions [99] (see also Appendix A of Ref. [100]) . In the case of the head-on collision, the odd-parity perturbations are identically zero, and the non-trivial components of even-parity perturbations (in the Regge-Wheeler gauge) are We then read offH m 0 ,H m 2 andK m from the close-limit expansion to the 2PN metric in Eq. (4). To model the ringdown radiation yet in the simple setup, we shall consider only the = 2 and 3 perturbations that satisfy the linearized Einstein equations at the first order in the (symmetric) mass ratio ν ≡ m 1 m 2 /M 2 up to terms O(G 2 , c 6 , r 4 12 ). This restriction gives With this relation, the explicit expressions of = 2 components on the time symmetric initial surface are (Eq. (3.6) of Ref. [83]; r 12 is the binary's separation and we assume that the relative orbital velocity v 12 , orbital angular velocity ω 12 and the orbital phase β are all zero) and those of = 3 components areH with δM ≡ m 1 − m 2 .

Numerical method
In this section, we shall use the close-limit, 2PN metric in Eqs. (10) and (11) (in the Regge-Wheeler gauge) as initial data, and numerically evolve it with the time-domain (TD) implementation of the Regge-Wheeler-Zerilli equations, to obtain the associated ringdown waveform.
Before proceeding, we note that Eq. (10) in the equal-mass limit (ν = 1/4) essentially agrees with the Brill-Lindquist geometry for the head-on collision [101], after a certain coordinate adjustment to the initial separation r 12 . This will provide reassurance that the physics in our TD numerical waveform should be largely independent of our specific choice of the PN initial data (see Section 5 of Ref. [83] for further elaboration).

Initial data and gravitational waveform in terms of master functions
In the practical implementation, it is more preferable to use gauge-invariant master functions that can be constructed fromH m CL andK m CL in the Regge-Wheeler gauge. For even-parity perturbations, we use (the Regge-Wheeler-gauge expression for) the Zerilli-Moncrief function given by (Eq. (5.1) of Ref. [83]; see also Eq. (126) of Ref. [98]) where λ ≡ ( − 1)( + 2)/2. Substituting Eqs. (10) and (11) into this, the = 2 and 3 multipolar coefficients for the head-on collisions are then translated to with . Furthermore, the Zerilli-Moncrief function in Eq. (12) is directly related to the GW strain where −2 Y m is the spin-weighted spherical harmonics with s = −2 (note that the contribution from the odd-parity perturbation identically vanishes in the case of the head-on collision). To avoid the dependence of the polarization on the observer's location, we focus on Ψ m in the remaining sections.

Time domain integration of the Zerilli equation
In our model problem of the head-on collision, the Zerilli-Moncrief function Ψ m satisfies the homogeneous Zerilli equation. We write the equation in the double null coordinates, in practice, where u = t − r * and v = t + r * with the tortoise coordinate r * = r + 2M ln[r/(2M) − 1], and is the Zerilli potential. To evolve Eq. (16), we use the finite-difference scheme developed by Lousto and Price [102]. We introduce an uniform null grid in the 1 + 1 numerical domain as the left panel of Figure 1. Consider a grid cell with the size of h × h in the domain. Integrating Eq. (16) over the cell gives . We impose the close-limit, PN initial condition in Eqs. (13) and (14)  where r c is the value of r at the center of the cell, and Ψ i for i = 1, 2, 3, 4 correspond to the values of Ψ (we have omitted the indices ( , m) here) at the vertices of the cell in the right panel of Figure 1. By using the above equation, we can obtain Ψ 1 from Ψ 2 , Ψ 3 and Ψ 4 with the local error term of O(h 4 ), The local error is accumulated through the integration over v and u to give a global error of O(h 2 ). The finite-difference scheme on the initial surface should be derived separately. To implement the time-symmetric initial conditions in Eqs. (13) and (14), we also assume the time symmetry on the initial surface, and therefore ∂ t Ψ| t=t 0 = 0, which leads to if the vertices 2 and 3 are on the surface. From this relation, we find that the finite-difference scheme on the initial surface is While the local error term in this case is O(h 3 ), it eventually gives the global error of O(h 2 ) because the error is accumulated only over the initial surface. Hence, the convergence of our numerical solution through the finite difference scheme in Eqs. (19) and (21) is quadratic with respect to h (see Figure 2 below).

Simulation parameters
In this paper, we calculate the Zerilli-Moncrief function with the close-limit, 2PN initial data for head-on collisions of two nonspinning BHs in the following two cases: • Case A: an equal mass collision with the initial separation of r 12 = 3.5M.
• Case B: an asymmetric mass collision with the mass ratio of ν = 0.2 and the initial separation of r 12 = 4.4M.
These system parameters are chosen so that they agree with those employed in Ref. [83]. In the case A, we focus only on ( = 2, m = 0) mode because ( = 2, m = ±2) modes are just proportional to this mode. Similarly, in the case B, we look at only ( = 3, m = 1) modes because ( = 3, m = −1) differs by only an overall minus sign, and ( = 3, m = ±3) modes are proportional to that mode; recall Eqs. (13) and (14).

Code validation
In our numerical simulations, the initial data given in Section 3.1 with the parameters shown in Section 3.3 is set on the region of −1000 ≤ r * /M ≤ 3000 at t = t i = 0 and the Zerilli equation in Eq. (16) are integrated within the future domain of dependence of the initial surface (the upper triangular domain in the left panel of Figure 1). The produced gravitational waveform is extracted on  Absolute difference

Figure 2.
Absolute differences between the results with different resolutions, h/M = (0.08, 0.04, 0.02, 0.01, 0.005). For example, Ψ h=0.04 − Ψ h=0.08 (the blue solid curve) denotes the absolute differences between ( = 2, m = 0) mode of the Zerilli-Moncrief function for the case A in the h/M = 0.08 and 0.04 simulations. The plot demonstrates the second order convergence of our time domain code. The numerical accuracy with the highest resolution is roughly estimated as 10 −7 at the time of the peak amplitude u peak and 10 −13 in the late time of (u − u peak )/M 200.
To check the convergence of our TD code, we consider the absolute differences between the Zerilli-Moncrief function Ψ (we have omitted the indices ( , m) here) computed with different resolutions, where h 1 and h 2 denote adjacent resolutions. The result for the ( = 2, m = 0) Zerilli-Moncrief function for the case A is shown in Figure 2. We can find that the difference becomes a quarter when the resolution gets a half. This shows the O(h 2 ) convergence of the data, as expected. The numerical accuracy with the highest resolution, h = 0.005M, is roughly estimated as 10 −7 at the time of the peak amplitude. In the late time of (u − u peak )/M 200, the numerical accuracy becomes 10 −13 .

Modeling of ringdown waveforms
In this section, we introduce waveform models to describe the ringdown data in our close-limit calculations, as a superposition of QNMs associated with the remnant Schwarzschild BH. The method to compute QNMs are well developed in literature (see, e.g., Ref. [103]). We use the accurate numerical data of QNM frequencies of the remnant BH provided by Berti [103] (for = 2) and those generated with the Black Hole Perturbation Club (B.H.P.C.) code [104] We should note that the zero-spin remnant BH here is a limitation of our close-limit approximation to the binary nonspinning BHs. When we treat general binary BH mergers in the full NR simulations with an astrophysically realistic initial data, the remnant BHs after merger are Kerr with a non-zero spin in general. For example, a nonspinning equal-mass binary BH of the quasicircular merger can produce a remnant Kerr BH with the nondimensional spin ∼ 0.67 [105].

Standard quasinormal-mode fitting formula
We write the ringdown waveforms as a sum of the fundamental tone and overtones of QNMs decomposed into spin-weighted spherical harmonics with angular indices ( , m) and spin s = −2, with where we truncate the summation with a finite overtone index of n = N although there are infinite overtones. Here, Q m n is a complex amplitude, ω QNM mn is the Schwarzschild QNM frequency of the ( , m, n) mode, and t 0 is the 'starting time' before which we do not include the model to fit our numerical waveforms. As mentioned later, we choose the time of the peak amplitude of the wave, t peak , as t 0 .
With an appropriate choice of the initial phase, our numerical waveforms of the head-on collisions are given as purely real functions of the retarded time, u. Therefore, instead of the complex representation of the waveform in Eq. (24), we adopt the following real expression for the Zerilli-Moncrief function in practice, making use of the relation in Eq. (15), with where C m k (k = 0, · · · , 2N + 1) are real constants, u 0 is the starting time in terms of the retarded time, ω mn and τ mn are derived from the real and imaginary parts of the QNM frequency, respectively, The prefactor in the mode functions φ m k (u) of Eqs. (26) and (27) is chosen so that they are certainly normalized (see the next subsection). We note that the model in Eq. (25) is a linear function with respect to the fitting coefficients C m k . This means that one can obtain the best fit values of C m k through linear least squares.

Modified quasinormal-mode fitting formula with an orthonormal set of mode functions
The mode functions φ m k (u) defined in the previous subsection are linearly independent but not orthogonal each other. For convenience of our analysis, we introduce an 'orthonormal set' of mode functions by using the Gram-Schmidt procedure. For this purpose, we first define the inner product between two real functions of f (u) and g(u) given in u ≥ u 0 For this definition, we can find that φ m k in Eqs. (26) and (27) is normalized as (φ m k , φ m k ) = 1. Using the inner product in Eq. (30) and the standard Gram-Schmidt procedure, we obtain the orthonormal set of mode functionsφ m k . Their explicit listing is . . .
However, we note that a direct numerical implementation of above construction does not work well because of the accumulation of round-off errors. In practice, we use the modified Gram-Schmidt algorithm [106] to obtain accurately the orthonormal set of mode functions. From the above orthonormal set of mode functionsφ m k , we can extract the fitting coefficients easily byC from the numerically generated Zerilli-Moncrief function Ψ m num , and construct the modified model waveform asΨ up to a finite overtone index of n = N.

Results
We now present the demonstration of the QNM fits including overtones given in Eqs. (25) and (37) with the close-limit TD waveform for head-on collisions explained in Section 3.       Figure 3 (the blue solid lines). We find that the residuals are all larger than our numerical uncertainty (the green dotted lines), and shows the power-law behavior in the late time (u − u peak > 100M, the orange dashed lines). This corresponds to the late-time tail in the evolution of a perturbation of the gravitational field around a remnant Schwarzschild BH, studied in various literature; see, e.g., Refs. [107][108][109] for the Schwarzschild background spacetime, and, e.g., Refs. [110,111] for the Kerr background spacetime. The power-law tail in a sufficiently late time is also investigated using a sophisticated numerical method with compactification along hyperboloidal surfaces [112]. Each of QNM fit residuals deviates from the power-law line when u becomes larger because the null infinity approximation v u is broken (see Section IV in Ref. [113]). In Figure 4, we also find that the QNM fit residuals in the early-time (u − u peak 100M) do not improve even if one includes more overtones into the model waveform in Eq. (37). Although we do not have a mathematically rigorous proof, we speculate that this limitation may arise from the QNM fit including overtones (that describes the exponential decay of the perturbation) to the power-law tail component of the perturbation, based on a simple numerical experiment detailed below. First, we prepare a mock data of the tail component of the perturbation as

Quasinormal-mode fitting and residuals
where D is an arbitrary constant, and we define a window function by with u peak < u 1 < u 2 . The window function is introduced to remove the divergent behavior of the late-time tail at u = u peak . Next, we fit this mock data using the QNM model including overtones in Eq. (37) in the same manner as we do for the TD data. Figure 5 shows the QNM fit including overtones to the mock tail data and corresponding fit residuals. We see that the residuals in the early-time ((u − u peak )/M 100) show the similar behaviour in Figure 4. This result supports our expectation about the early-time residual in the TD data.

Convergence of the fitting coefficients
Because the model waveform in Eqs. (25) and (37) involves many overtones k (or N) it is instructive to check their convergence in terms of fitting coefficients on k, to ask how higher overtones should be included in practice when we wish to model the TD data from the peak time at u = u peak .
In the left panel of Figure 6 where the TD data is for the ( = 2, m = 0) mode in the case A, we show the fitting coefficients for the original QNM basis C m k (the orange filled triangles) and those for the orthonormal set of mode functionsC m k (the blue filled circles) as a function of overtones up to k ≤ 21 (i.e., N = 10). It is found that |C m k | decreases roughly in power law of k, while |C m k | does not. Due to the benefit of the orthonormalization of mode functions in Eq. (37), it can be conveniently used to assess the fraction of power in each mode Based on the fact that one more mode is added to the orthonormal set when k increases by two, we introduce an estimator to characterize the match (or overlap) between the TD data and the QNM fit model including overtones at a given n The right panel of Figure 6 shows ∆ n as a function of mode n. We see that the inclusion of the first overtone (n = 1) can improve the match by 10%, and that the QNM fits up to n = 5 will suppress the mismatch less than 0.1%.
In the top and bottom panels of Figure 7, we show the same figures as for Figure 6 but with the ( = 2, m = 0) and ( = 3, m = 1) modes in the case B, respectively. We see the similar convergent behavior inC m k here, and again find that (n ≥ 5) higher overtones contributes to the match only at the level of less than 0.1%. These results highlight the potential benefit to use the orthonormal set of mode functions in the data analysis of ringdown GWs.

Summary and Discussion
It is now widely recognized the importance to have a ringdown waveform model including higher overtones of QNMs (see, e.g., Ref. [47]). With the overtones, one can set the starting time of the ringdown GWs at the time of the peak amplitude much earlier than a time (∼ 10M -20M after the peak time) to obtain unbiased remnant BH parameters (see, e.g., Figure 5 for GW150914 in Ref. [40]). This allows the increased SNR of observed ringdown signals as well as better parameter estimation of the remnant BH.
In this paper, we examined in detail the ringdown GWs of binary BHs with QNM fits including overtones in Eqs. (25) and (37), using the accurate close-limit waveform in the BH perturbation theory as a test bed. Our analysis is restricted to the narrow case of head-on collisions of two nonspinning BHs  (based on the 2PN initial data and linear BH perturbation theory), but it suffices to highlight some of the key features of the full problem, and we found three main results: (i) we reconfirm the importance of QNMs overtones to fit the ringdown waveforms after the time of the peak amplitude. This agrees with the previous findings in literature [47,72]; (ii) the small contributions of late-time tail exist in the fit residuals at the level of O(10 −5 ) or below (see Figure 4); (iii) the fitting coefficients decay with overtones when one uses the orthonormal set of mode functions in Eq. (36) (see Figures 6 and 7).
A natural next step of our study is to assess the merit of the orthonormal set of mode functions in the context of GW data analysis, accounting for the detector's noise. A standard technique for the analysis of ringdown GWs is the frequency-domain approach (see, e.g., Ref. [114] for various methods to extract the ringdown GWs), and we can directly apply the orthonormal set of mode function presented in Section 4.2 with the (so named) noise-weighted inner product in the matched filtering method (see, e.g., Ref. [115]). The implementation of this method to the GW data analysis and the analysis of real data from GW detectors will be a future work.
Another interesting future work would be to add BH's spin to our close-limit waveform. Although it is known for any field points that the 2.5PN near-zone metric with spin effects [95,96], the numerical computation of either metric or curvature perturbations in Kerr spacetime is particularly challenging in the time domain [116][117][118], and would need much coding efforts. It is probably more viable, instead, to perturbatively include the spin effect into our non-spinning close-limit waveform (e.g., [119]). It may be also possible to calculate the close-limit waveform in Kerr   technique (and taking the extreme mass-ratio limit). Because, in either case, the extension is involved, we leave them for future work.
From the theoretical point of view, it will be useful to include information on excitation coefficients and factors for the fundamental tone and overtones of QNMs derived from the BH perturbation approach (See, e.g., Refs. [120][121][122][123]). Also, although we focused only on the ( = 2, m = 0) and ( = 3, m = 1) modes in the bulk of this paper (and will focus only on the ( = 2, m = 2) mode in Appendix A), the inclusion of the other subdominant (higher harmonic) modes will be beneficial. In particular, the ( = 4, m = 4) mode of NR waveforms for binary BHs is especially interesting because this mode includes the second order modes composed by the non-linear, self-coupling of first-order ( = 2, m = 2) × ( = 2, m = 2) mode, and it has already been found in NR simulations [73].
Nevertheless, astute readers may ask what extent the results based on the close-limit approximation (to the head-on collisions of two nonspinning BHs) are consistent with the full NR waveforms in more astrophysically relevant scenarios (e.g., coalescences observed by LIGO/Virgo). We conclude our paper to briefly answer this question, analyzing the NR waveform of SXS:BBH:0305 in the Simulating eXtreme Spacetimes (SXS) catalogue [66] with the same approach that we took in Section 5; Appendices A and B also provide supplemented examples of the analysis for NR waveforms. Figure 8 shows the ( = 2, m = 2) spheroidal mode of the waveform SXS:BBH:0305 [66,124,125] (see Eq. (A1) for the construction of the spheroidal modes from the spherical ones). Here, we set the time of the peak amplitude as the starting time of the QNM fit (see the orange dashed line in the left panel of Figure 8), and we use the QNM fitting formula introduced in Eq. (37). Unlike the TD data generated within the close-limit approximation and the BH perturbation approach, we cannot confirm the power-law tail certainty due to some constant fit residuals of ∼ 10 −5 after t − t peak 80M; recall that the tail contribution in Figure 4 is at the level of O(10 −5 ), which is as large as the remaining residual in the NR waveform. We expect that future NR simulation with higher numerical accuracy (or, e.g., the close-limit approximation to bound-orbit BBH mergers) will provide a more robust answers of the need of the tail contribution in the ringdown waveform modeling, and eventually in the measurement of remnant properties of BHs. Figure 9 displays the fitting coefficients for the original QNM basis C m k and those for the orthonormalized QNM basisC m k introduced in Eq. (36), again focusing on the the ( = 2, m = 2) spheroidal mode of SXS:BBH:0305. We see that the trend of C m k andC m k are similar to that of the close-limit, TD data displayed in Figures 6 and 7 (see also Appendix B). Namely,C m k converges as increasing the overtone k while C m k (from the least square fit by using Eq. (25)) does not. The convergence property ofC m k supports the argument that the dominant ( = 2, m = 2) mode of NR waveforms after the time of the peak amplitude can be described only by the QNMs in the linearized BH perturbation (see Section 1 for the other confirmations).

Author Contributions:
The authors contribute equally to this paper.

Acknowledgments:
We would like to thank Ryuichi Fujita and Takahiro Tanaka for useful discussion. S. I. is grateful to Leor Barack for his continuous encouragement.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Analysis of some numerical-relativity waveforms in time domain
In this appendix and Appendix B, we focus on the ( = 2, m = 2) (spheroidal harmonic) mode of the ringdown signal (after the time of the peak amplitude) provided by NR waveforms for nonprecessing, spinning binary BH mergers, and analyze them following the same approach as in Section 5. Our goal here is to examine if we find i) the late-time, power-law tail (this appendix) and ii) the convergence of the orthonormalized QNM fits (Appendix B) in more astrophysically relevant NR waveforms. We note that there are also various works on the mismatch and parameter estimation errors based on the ringdown portion of the NR waveforms and their QNM fits with overtones [47,[51][52][53][54]57], memory (mainly the m = 0 mode) [55] and mirror (negative m) modes [56,58].
We find that all the SXS waveform data in the late time have the same constant residue of O(10 −5 ) as for SXS:BBH:0305, except for BBH:1477 that has the shorter time stretch between the peak and the end of data including the original data than other data set. We cannot identify the late-time tails in the QNM fit  Table A1 and their QNM fit residuals. Each plot corresponds to that in the right panel of Figure Table A1 and their QNM fit residuals. residuals (at least) in these NR waveforms, and some improvement in accuracy of NR simulation may be needed in order to reveal the tails.
We also find the slower damping fit residuals than the ( = 2, m = 2) mode itself in the late-time part of SXS:BBH:0178 and SXS:BBH:1124 (with highly spinning remnant BHs, χ rem ∼ 0.95). This is rather unexpected result because any QNMs (even including the second and higher order perturbations) have a faster damping time than the fundamental tone ( = 2, m = 2, n = 0) in the high-spin range; see Figure A5 in Appendix C.
Before concluding this appendix, we should note the possible second-order contribution, i.e., the self-coupling of the two first-order QNMs computed from the linear perturbation theory, to our analysis. For example, the second-order, ( = 2, m = 2) × ( = 2, m = 2) mode is found in NR simulations [73]. This mode is corresponding to the second order ( = 4, m = 4) mode, and we did not consider it here. Also, in principle, the ( = 2, m = 2) mode include the second-order contribution like ( = 2, m = 1) × ( = 2, m = 1) mode, but we cannot find these second-order contributions in our analysis. This is consistent with the previous work in Ref. [73].

Appendix B. Analysis of the fitting coefficients of some numerical-relativity data
As a by-product of the analysis in Appendix A, we obtain the fitting coefficients (in terms of the original QNM basis and the orthonormal set of mode functions) for the NR data listed in Table A1, and can check the convergence in the same manner as shown in Figure 9. We summarize the results in Figure A3 for the SXS data and Figure A4 for the RIT data. From these results, we expect that the contribution from the overtones to the ringdown waveform becomes more significant with larger remnant spins. The further study will be required to clarify the relation between the parameters of binary BHs and the significance of the overtones (A comprehensive study of the dependence of the mode excitations on the source parameters has been done in the extreme mass-ratio limit [122,123]).

Appendix C. Frequencies of Kerr quasinormal modes
Although many plots for Kerr QNM frequencies are available in literature, e.g., in Refs. [8,133,134], we reproduce a QNM figure here for m = modes using publicly available data provided in "Ringdown" by Berti [103] (see also "Kerr Quasinormal Modes: s = −2, n = 0-7" by Cook [135]). They give the lowest damping rate in each mode, and help to check the presence of the second order perturbation in our analysis. Below the indices ( , m) refer to the spin-weighted spheroidal basis with s = −2. Figure A5 shows the Kerr QNM frequencies (the fundamental tone and overtones up to n = 7) for the nondimensional spin 0 ≤ χ ≤ 0.999, but it suffices to consider the QNM frequencies up to the filled circles (χ = 0.95) because the NR simulations used for our study have the remnant BH spin χ 0.95 after merger. In particular, we see that the ( = 2, m = 2, n = 0) mode has the lowest damping rate (recall Eq. (29)). Importantly, the first order ( = 2, m = 2, n = 0) mode remains to be the longest-lived mode even when considering the second and higher order QNMs as we establish below.
Suppose an expansion of the metric of the form, µν + (higher-order terms) , where g µν is a background (Kerr) spacetime, and h (1) µν is the first order perturbation here. The second order perturbation h (2) µν formally satisfies the linearized Einstein equation where G (1) µν is the linearized Einstein tensor (in a given gauge), and G (2) µν [h (1) , h (1) ] is the term quadratic in h (1) in the expansion of the Einstein tensor. Although the real part of second-order QNM frequencies of h (2) µν can take various values due to the self-mode coupling of h (1) µν through G (2) µν , Figure A5 indicates that for χ 0.95 the imaginary part of any second-order QNM frequencies ω (2) I is bounded at [74,75] where ω (1) I, 220 is the imaginary part of the first order ( = 2, m = 2, n = 0) QNM frequency. This in turn implies that the damping time of any second order QNMs satisfies Therefore, the first order ( = 2, m = 2, n = 0) QNM has the slowest damping time, at least to the second order in the expansion of the metric and for 0 ≤ χ ≤ 0.95 of the remnant BH (whose imaginary part of the first order QNMs is not too close to zero).  Here ω R and ω I mean the real and imaginary parts of the QNM frequency, respectively. We plot the values for 2 ≤ ≤ 7 in the range of the nondimensional spin of 0 ≤ χ ≤ 0.999. The filled circles show the QNM values at χ = 0.95. For χ ≤ 0.95, the n = 0 mode always has the highest real frequency and the slowest damping time. There is an interesting feature in the fifth overtone (n = 5) of the m = = 2 mode. To make these plots, we use the QNM data provided in Ref. [103].