Extraction of multiple quasinormal modes QNMs from ringdown gravitational waves emitted from a binary black hole coalescence is a touchstone to test whether a remnant black hole is described by the Kerr spacetime in general relativity. However it is not straightforward to check the consistency between the ringdown signal and the QNM frequencies predicted by the linear perturbation theory. While the longestlived mode can be extracted in a stable manner the higher overtones damp more quickly and hence the fitting of overtones tends to end up with the overfit. To improve the extraction of overtones we propose an iterative procedure consisting of fitting and subtraction of the longestlived mode of the ringdown waveform in the time domain. Through the analyses of the mock waveform and numerical relativity waveform we clarify that the iterative procedure allows us to extract the overtones in a more stable manner.
['gravitational waves', 'black hole', 'quasi-normal modes', 'General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena']
Iterative extraction of overtones from black hole ringdown
{'Iterative extraction of overtones from black hole ringdown': 'Kazuto Takahashi 1, ∗ and Hayato Motohashi 2, †  \n1 Department of Electrical and Electronic Engineering, Tokyo Institute of Technology,  \n2-12-1 Ookayama, Meguro-ku, Tokyo, 152-8552, Japan  \n2 Division of Liberal Arts, Kogakuin University, 2665-1 Nakano-machi, Hachioji, Tokyo, 192-0015, Japan  \nExtraction of multiple quasinormal modes from ringdown gravitational waves emitted from a binary black hole coalescence is a touchstone to test whether a remnant black hole is described by the Kerr spacetime in general relativity. However, it is not straightforward to check the consistency between the ringdown signal and the quasinormal mode frequencies predicted by the linear perturbation theory. While the longest-lived mode can be extracted in a stable manner, the higher overtones damp more quickly and hence the fitting of overtones tends to end up with the overfit. To improve the extraction of overtones, we propose an iterative procedure consisting of fitting and subtraction of the longest-lived mode of the ringdown waveform in the time domain. Through the analyses of the mock waveform and numerical relativity waveform, we clarify that the iterative procedure allows us to extract the overtones in a more stable manner.', 'I. INTRODUCTION': 'At the final stage of a binary black hole coalescence, where a remnant black hole is formed, emitted gravitational waves have a peak strain amplitude followed by damped oscillations. Ringdown is such regime of the damped oscillations of gravitational waves emitted during the remnant black hole going down to a stationary state. Black hole linear perturbation theory predicts that the ringdown gravitational waves are described by a superposition of exponentially damped sinusoids characterized by complex frequencies known as the quasinormal modes (QNMs) [1-4]. The QNMs are discrete complex frequencies ω ℓmn , which are, in addition to the two angular indices ( ℓ, m ), labeled by overtone index n according to the damping time -1 / Im( ω ℓmn ). The longest-lived mode n = 0 is called the fundamental mode, and modes with n ≥ 0 are n -th overtones and have shorter damping time. Since the QNM frequencies of a Kerr black hole are fully determined by mass and spin [5], the measurement of the frequency and damping time of the fundamental mode would provide an estimation of the mass and spin of the remnant black hole [6, 7]. Further, an independent measurement of multiple modes would allow us to test whether the remnant black hole is consistent with the Kerr spacetime, known as the black hole spectroscopy [8, 9].  \nIt is known that the QNM spectrum is unstable with respect to a small change of the effective potential [1016]. Such small perturbation would be caused by nonlinear effects in dynamical formation of a deformed black hole and/or surrounding material. However, the ringdown waveform in the time domain is affected only at late time, which would be actually subject to noises in the observational signals, and hence it is well described by the unperturbed QNMs [10, 11, 17-19]. In this sense,  \nthe black hole spectroscopy program is robust and still works. The apparent discrepancy is due to the fact that the frequency domain analysis and time domain analysis are not necessarily equivalent when the full data is not available. Thus, given that we cannot observe the signal for a very long time and with sufficiently high precision, it is important to analyze both the frequency domain and the time domain in a complementary manner.  \nSince the merger of binary black holes is highly nonlinear process, it is in principle nontrivial if the ringdown gravitational waves can be described by a superposition of QNM damped sinusoids obtained by the linear perturbation theory. Nevertheless, from the firstly observed ringdown gravitational waves GW150914 [20], the fundamental mode of ( ℓ, m ) = (2 , 2) mode is robustly extracted [21]. However, the extraction of the overtones is not straightforward. The possibility of detection of the overtones from the GW150914 is still under debate [2231]. The ringdown analyses of GW190521 has also been focused recently [23, 32-36]. Since the observational ringdown data is inevitably subject to noises, it would be worthwhile to revisit a simpler and ideal setup such as mock waveform and/or numerical relativity waveform.  \nActually, there has also been an active discussion on the possibility of the extraction of the overtones from numerical relativity waveform. It was shown in [37] that the fit starting from the peak of the strain amplitude improves as one includes more overtones in the fitting function. However, the significance of the overtones in numerical relativity waveform was questioned in [38, 39] and the possibility of the overfit was pointed out. It was clarified that the fit of the overtones is sensitive to the start time of the data interval used for the fit. Further, it was shown that one can improve the mismatch even with overtones of frequencies irrelevant to the Kerr QNMs corresponding to the remnant spin. The difficulty of extraction of overtones is also pointed out in [40] with an analysis of simulation of head-on collision of two nonspinning black holes. An interesting approach was considered in [41], where the fundamental mode and lower overtones are filtered out to reveal subdominant effects  \nin the ringdown such as nonlinear effects and sphericalspheroidal mode mixing [42]. Such an approach would be useful to clarify the stability of the extraction of the overtones, but the frequency domain filter not only removes the mode of interest but also modifies the contribution from other modes. It would be worthwhile to consider a time domain approach to remove the dominant mode to help the extraction of subdominant modes.  \nThe aim of the present paper is to scrutinize the extraction of the overtones from the ringdown gravitational waves. In contrast to the conventional fitting method, where one extracts the coefficients of all the modes of interest from a single fitting, we consider an iterative procedure of fitting and subtraction to extract the longest-lived mode at each step. Starting from analyses of mock waveform, we shall see that we can fit the longest-lived mode in the most stable manner compared to the subdominant modes, and when we subtract a damped sinusoid of the longest-lived mode from the time domain waveform, this nature can be taken over to the next-longest-lived mode. Hence, this procedure improves the stability of the fit of the subdominant modes. We clarify that, by making use of iterative fitting and subtractions of the longest-lived mode, we can improve the stability of the overtones compared to the conventional fit.  \nThe rest of the paper is organized as follows. In § II, we review the fitting algorithm used throughout the present paper. We investigate the iterative fitting of mock waveforms in § III and § IV. First, in § III, we deal with a mock waveform composed of a superposition of pure damped sinusoids. Second, in § IV, we add a mock waveform a small constant mimicking noise and/or tail in the ringdown signal. For these toy examples, we shall see how we can improve the fit of overtones by iteratively extracting the longest-lived mode. In § V, we apply the iterative fitting method to numerical relativity waveform in the Simulating eXtreme Space-times (SXS) catalog [43, 44] and show an improvement of the extraction of overtones. § VI is devoted to conclusion and discussion.  \nThroughout the present paper, we work in the natural unit where c = G = 1, and employ the total binary mass M before merger to normalize physical quantities such as time and frequency. We mainly focus on ℓ = m = 2 mode, so we often omit the subscript ℓ and m for ω ℓmn etc. when no confusion occurs. Numerical calculations are performed with machine precision but when we present specific numbers, values are rounded to four significant digits unless otherwise noted.', 'II. FITTING ALGORITHM': "In this section, we briefly review the fitting algorithm known as the eigenvalue method [45], which we exploit throughout the present paper. While the optimization algorithms such as curve fit in scipy or NonLinearModelFit in Mathematica are widely used in the literature, it is known that for these algorithms the  \nresults depend on the initial guess for the fitting parameters. Such ambiguities are not ideal when we discuss the possibility of overfit. The eigenvalue method does not require the initial guess and hence provides a robust result.  \nLet Ψ( t ) be a given ringdown waveform of gravitational wave strain which we would like to fit by a fitting function. Regarding Ψ( t ), we shall consider mock waveform in § III and § IV, and numerical relativity waveform in § V. Regarding the fitting function, we use the following function:  \nψ fit [ K,N ] ( t ) = N ∑ n = K C n ψ n ( t ) , (1)  \nwhich is a superposition of the damped sinusoids ψ n ( t ) = e -iω n t with complex frequencies ω n . In the present paper, we perform the frequency-fixed fitting, i.e., we determine complex coefficient for each given QNM mode by the fitting algorithm. We thus set ω n the QNM frequencies corresponding to the given waveform Ψ( t ). Specifically, we use the fiducial value frequencies for the analysis of the mock waveform, and the Kerr QNM frequencies associated with the remnant black hole for the analysis of the numerical relativity waveform. The fitting parameters are complex coefficients C n . For the following, we shall also use the notation C n = A n e iϕ n with amplitude A n and phase ϕ n being real parameters.  \nTo quantify the goodness of fit, let us define the overlap ρ between the ringdown waveform Ψ and the fitting function ψ fit by  \nρ 2 = |⟨ ψ fit | Ψ ⟩| 2 ⟨ Ψ | Ψ ⟩⟨ ψ fit | ψ fit ⟩ , (2)  \nwhere  \n⟨ ψ | ϕ ⟩ = ∫ t e t i ψ ∗ ϕ dt. (3)  \nHere, ψ ∗ is the complex conjugate of ψ , t i is the start time of the fit, and t e is the end time of the fit. By introducing the notation B n = ⟨ ψ n | Ψ ⟩ and D nm = ⟨ ψ n | ψ m ⟩ = ( e i ( ω ∗ n -ω m ) t e -e i ( ω ∗ n -ω m ) t i ) / [ i ( ω ∗ n -ω m )], we can rewrite Eq. (2) as  \nρ 2 = | ∑ n C ∗ n B n | 2 ⟨ Ψ | Ψ ⟩ ∑ n,m C ∗ n D nm C m . (4)  \nWe determine the coefficient C n by maximizing the overlap. Requiring a derivative of ρ 2 with respect to C n to vanish, we obtain  \nC n = ∑ m ( D -1 ) nm B m . (5)  \nHere, D -1 is the inverse of a matrix D , whose argument is D nm . The maximum ρ 2 max of the overlap is then given by  \nρ 2 max = ∑ n,m B ∗ n ( D -1 ) nm B m ⟨ Ψ | Ψ ⟩ . (6)  \nWe can quantify the goodness of fit by using the mismatch defined by  \nM = 1 -ρ max . (7)  \nSince the ringdown waveform is given as a discrete data of ( t a , Ψ( t a )), we need to evaluate the integrals B n and ⟨ Ψ | Ψ ⟩ in Eq. (6) numerically. In doing so, we exploit the Romberg integration T 2 , 2 [46]. Although the Simpson integration T 1 , 1 is commonly used for the numerical integration, we find that it does not provide a sufficient accuracy for our purpose. Namely, the mismatch calculated by the Simpson's rule sometimes becomes negative. Therefore, following the Romberg's method with the Richardson extrapolation, we improve the accuracy of numerical integration by employing T 2 , 2 , which we use for the numerical integration throughout the present paper.", 'III. PURE DAMPED SINUSOIDS': 'Before considering the extraction of QNMs from numerical relativity ringdown waveform, we study mock ringdown waveform. For the following, we shall study the fitting of two kinds of mock waveforms. In this section, we consider a mock waveform composed of a superposition of fiducial damped sinusoids. In § IV, we shall add the mock waveform a constant to mimic the noise and/or tail of the ringdown signal.  \nLet us consider a mock waveform of a superposition of fiducial damped sinusoids:  \nΨ [0 , 7] ( t ) = 7 ∑ n =0 a n e -i ( ω n t -θ n ) . (8)  \nAs a representative mock waveform, we consider a waveform mimicking the GW150914-like numerical relativity waveform SXS:BBH:0305 in the Simulating eXtreme Space-times (SXS) catalog [43, 44]. We set ω n the QNM frequencies for ( ℓ, m ) = (2 , 2) mode of Kerr black hole with the spin parameter 0 . 6921, which corresponds to the remnant dimensionless spin of SXS:BBH:0305. The QNM frequencies are listed in Table VIII in Appendix A. The n = 0 mode is known as the fundamental mode, and n ≥ 1 is the n -th overtone. The damping time -1 / Im( ω n ) decreases as the overtone index n increases. We set the fiducial value for a n and θ n as listed in Table I, which we determine by the best-fit value obtained by the conventional fit (see Table III). Namely, we fit the gravitational strain Ψ SXS = ( h + -ih × ) r/M of ( ℓ, m ) = (2 , 2) mode of SXS:BBH:0305 by the fitting algorithm described in § II with a fitting function ψ fit [0 , 7] . Here, h + and h × are two orthogonal polarization modes of gravitational waves, and we normalize the strain by multiplying r/M , where r is the distance between the binary and observer and M is the total binary mass before merger. We generate a discrete data of ( t a , Ψ( t a ))  \nFIG. 1. Mock waveform Ψ [0 , 7] (blue) given in Eq. (8) and Ψ [1 , 7] (orange) after the subtraction of the fundamental mode.  \n<!-- image -->  \n-  \nusing the same sequence of sampling time t a of the simulation SXS:BBH:0305 with the resolution Lev6 and the extrapolation order N ext = 2. Typical time data spacing is about ∆ t a /M = O (10 -1 ) for the ringdown phase.  \nTABLE I. Fiducial values used for the mock waveforms Ψ [0 , 7] in (8) and Ψ ( c ) [0 , 7] in (11).  \nThe mock waveform Ψ [0 , 7] given in (8) is dominated by the fundamental mode at late time, whereas the contribution of the overtones mainly affect the early time waveform. Figure 1 depicts | Re(Ψ [0 , 7] ( t -t peak )) | (blue) and | Re(Ψ [1 , 7] ( t -t peak )) | (orange) as a function of time t -t peak , where t peak is the peak time at which the strain amplitude | Ψ SXS | = √ h 2 + + h 2 × of the numerical relativity waveform has a maximum. To compare with the fit of the simulated ringdown waveform which we shall analyze in § V, we write the time t -t peak measured from the peak time. The late-time waveform of | Re(Ψ [0 , 7] ) | is dominated by a single damped sinusoid, which is the fundamental mode a 0 e -i [ ω 0 ( t -t peak ) -θ 0 ] . On the other hand, the late-time waveform of | Re(Ψ [1 , 7] ) | is dominated by the first overtone a 1 e -i [ ω 1 ( t -t peak ) -θ 1 ] , whose damping time is shorter than the fundamental mode. Since Ψ [0 , 7] is purely a superposition of damped sinusoids, for  \nFigures 3 and 4 respectively show the amplitudes A n and phases ϕ n for the coefficients in the fitting function (1) for each N , where the dotted lines indicate the fiducial values. While the mismatch M monotonically improves as we increase the number of overtones included in the fitting function, as demonstrated in [38, 40], the t i -dependence of the fit sheds light on different aspects of the fitting. It helps us to understand overfit, stability, and correctness of the fit of overtones.  \n<!-- image -->  \n-  \nFIG. 2. Mismatch M between the mock waveform Ψ [0 , 7] and the fitting function ψ fit [0 ,N ] .  \nt < t peak the waveform diverges and does not have physical meaning.', 'A. Conventional fit': 'First, we apply the conventional fitting method to the mock waveform. Namely, we fit the mock waveform Ψ [0 , 7] ( t -t peak ) by using the fitting algorithm described in § II and the fitting function ψ fit [0 ,N ] ( t -t peak ), and obtain the best-fit values C n for all modes n = 0 , · · · , N from a single fitting. In doing so, we need to specify the region of the data used for the fit. Let us denote the time interval for the region for the fit by [ t i , t e ]. It is important to clarify how much the fitting is sensitive to the choice of the start time of the fit t i and the end time of the fit t e . As we can see from Fig. 1, the damped sinusoids continue to the end of the plotted region. We confirm that the result is insensitive to the choice of the end time t e of the fit so long as we set t e sufficiently late, say ( t e -t peak ) /M ≳ 100. As a representative value, we set ( t e -t peak ) /M = 150.  \nOn the other hand, as we shall see below, the goodness of the fit depends on the choice of the start time t i of the fit. Figure 2 shows the mismatch M between the mock waveform (8) and the fitting function (1) as a function of the start time of the fit t i . For the fitting function (1), we consider the sum from the fundamental mode K = 0 to the N -th overtone. For t i -t peak > 0, the slope of M are different for each N . The larger N , the earlier start time of the fit t i at which the mismatch M reduces. This trend is consistent with [37]. This is because the overtones contribute to the waveform at earlier times. In particular, the lower overtones improve the fit significantly even when they are subdominant. For the case with N = 7, since the fitting function is given by a superposition of the the same damped sinusoids used  \nin the mock waveform, we expect the mismatch becomes tiny. Indeed, in Fig. 2, the mismatch for the case with N = 7 remains O (10 -8 ) with fluctuations, and this behavior is expected to be due to errors in the numerical integration.  \nWe can see from Fig. 3 that for the early time, the overtones in the fitting function (1) can fit the overtones in the mock waveform (8) correctly, matching the amplitude A n to the fiducial value. However, for the late time, the overtones becomes much smaller than the envelope of the waveform, and the fitting becomes more difficult. As we move the start time of the fit t i later, the amplitudes for the higher overtones obtained from the fit diverge, and this behavior can be interpreted as the overfit. Since the divergence of A n compensates the damping of the overtone ψ n , the contribution of the overtone in the fitting function remains roughly constant even for the fit starting from the late time. Such behavior is not consistent with the expected decay of the contribution of the overtones in the mock waveform. Therefore, if we set t i late time, the higher overtones in the fitting function cannot extract the overtones in the mock waveform correctly, but overfit the mock waveform.  \nWe can see the same trend in the phase ϕ n in Fig. 4. The region of the plateau consistent with the fiducial values in Fig. 4 matches to the one in Fig. 3. There are some places where the value of ϕ n suddenly jumps by π , such as n = 1 in N = 6. This is due to the fact that the coefficients C n = A n e iϕ n actually pass near the origin in the complex plane. Indeed, we can observe the corresponding sharp dips in Fig. 3.  \nFurther, Figs. 3 and 4 allow us to check the stability and correctness of the fit. If there is a flat region, or plateau, where the variations of A n and ϕ n are small with respect to t i , we can say that the fit is stable. In addition, if the values of A n and ϕ n at the plateau are consistent with the fiducial values a n and θ n , we can say that the fit extracts the correct amplitude of the mode. We can see that as we increase the number of overtones included in the fitting function, the fitting of the fundamental mode and lower overtones becomes more stable. Specifically, the plateau of n = 0 and n = 1 is extended as we increase N , and ( A 0 , ϕ 0 ) and ( A 1 , ϕ 1 ) are consistent with the fiducial values. However, if one includes too many overtones in the fitting function, the fitting of the lower overtones are actually destabilized, as we shall discuss below.  \nThe N = 7 case corresponds to the fit of the mock waveform (8) by the fitting function (1) with the same  \nFIG. 4. Phase ϕ n for the fit of the mock waveform Ψ [0 , 7] by the fitting function ψ fit [0 ,N ] .  \n<!-- image -->  \n-  \n-  \n-  \n-  \nFIG. 3. Amplitude A n for the fit of the mock waveform Ψ [0 , 7] by the fitting function ψ fit [0 ,N ] .  \n<!-- image -->  \n-  \n-  \n-  \n-  \nnumber of the overtones. Naively, one may expect that the fitting works well for this case since the mismatch remains indeed small as we saw in Fig. 2. However, in Figs. 3 and 4, while the fitting reproduces the fiducial value of ( A n , ϕ n ) around t i = t peak , the plateau is not so long for higher overtones. This is consistent with the  \nresults obtained in [38]. As we stressed in § II, in our case the fitting algorithm does not require the initial guess, and there are no ambiguities. Nevertheless, higher overtones cannot be extracted at late time. This may be due to the fact that the higher modes damp faster and become negligible at the later stage. Therefore, even if we  \nfit the mock waveform of n = 0 to 7 modes by the fitting function composed of the same modes, the fit does not reproduce the correct value for all region of t i , especially for higher overtones.  \nIt is also worthwhile to note that the fit of lower overtones are actually destabilized as we increase the number of overtones included in the fitting function. For instance, in Fig. 3, the time when n = 2 curve starts to diverge is ( t i -t peak ) /M ≃ 40 for N = 3 , 4, but it becomes earlier for N ≥ 5. The fit of the first overtone n = 1 is also destabilized in N = 6 and 7. Although a seemingly counterintuitive result, this observation indicates that increasing the number of modes used in the fitting function does not always result in a better fitting.  \nOn the other hand, at least in the linear perturbation theory, the ringdown waveform includes an infinite number of QNM damped sinusoids. Thus, we cannot include the same number of modes in the fitting function of a superposition of a finite number of modes. So long as we deal with the fitting function with the finite number modes, we always fit the ringdown waveform with fewer number modes. To consider a similar situation in the analysis of the mock waveform (8) up to the seventh overtone, it would be reasonable to focus on the case where the fitting function (1) includes the fewer number modes. Specifically, we focus on the results with N ≤ 5.  \nTo extract coefficients consistent with the fiducial values, we need to choose a suitable start time of the fit. It is common to adopt the start time of the fit as the time when the decrease of the mismatch saturates. However, as we saw above, the mismatch itself is not sufficient to check if the overtones are overfitted. We would like to develop an alternative criterion to choose a start time of the fit taking into account the stability of the fit.  \nTo visualize the accuracy of the fit, let us consider the relative error between C n obtained by the fit and the fiducial value C fid n = a n e iθ n , i.e.,  \n∣ ∣ ∣ ∣ C n C fid n -1 ∣ ∣ ∣ ∣ . (9)  \nIn Fig. 5, we present the relative error as a function of the start time of the fit t i . As N increases, the relative error has a minimum at earlier time. For the mock waveform analysis, we can identify the most favorable start time of the fit as a start time of the fit t i with which the relative error from the fiducial value is minimum. However, we are ultimately interested in the fitting of the ringdown signal where the fiducial value is unknown. We thus would like to consider an indicator that helps us to extract such start time of the fit that gives the best-fit value without using the knowledge of the fiducial value.  \nAs such an indicator, it is useful to quantify the flatness of the plateau of the coefficient C n = A n e iϕ n . Let us introduce the rate γ n of change of the coefficient C n with respect to the start time of the fit t i as  \nγ n = ∣ ∣ ∣ ∣ 1 C n dC n dt i ∣ ∣ ∣ ∣ . (10)  \nA time t i when γ n takes a minimum value indicates that the rate of change of C n is smallest. We expect that we can use such a t i to extract the best-fit values close to the fiducial values. In practice, since the data is discrete, the derivative in (10) is understood as a finite difference ∆ C n / ∆ t i .  \nFigure 6 shows γ n as a function of the start time of the fit t i . If γ n persists small for some range of t i , we can say that the A n and ϕ n have a plateau there. As a concrete example, let us define the plateau region as a region that satisfies the criterion γ n < 10 -2 . For instance, for the case of N = 6, based on the criterion, we can say from Fig. 6 that there exists plateau up to fourth overtone, whereas there is no plateau in the fifth and sixth overtone. Indeed, we confirm that there is a correspondence between the plateau in Figs. 3 and 4 and the region where γ n is small in Fig. 6. As we increase N , we can see that γ n becomes smaller at earlier start time of the fit t i , which is consistent with the behavior of the plateau in A n .  \nFrom Figs. 3-6, we can infer that it is reasonable to adopt the best-fit value of C n by the fitting starting from t i where γ 1 is minimum. We see that there is a correspondence between t i where γ 1 is minimum in Fig. 6 and t i when the relative error of C n ( n = 0 , · · · N ) have the minimum in Fig. 5. We also confirm that there exist the plateau consistent with the fiducial value in Figs. 3 and 4 around t i where γ 1 is minimum in Fig. 6. The t i when γ n ≥ 2 takes the minimum is almost the same as the t i when γ 1 takes the minimum, and there are little changes in C n around those times.  \nNote that the start time of the fit t i when γ n is minimum is more or less the same, except the longest-lived mode, which is n = 0 in the present case. From Fig. 6, we see that γ 0 takes the minimum at later start time of the fit t i compared to γ n ≥ 1 . In such late times, we expect that only the fundamental mode contributes to the fit. Indeed, for instance, as we can see in the panel N = 3 of Fig. 6, γ 0 is minimum at ( t i -t peak ) /M ≃ 40, but there is no corresponding plateau in A n in Fig. 3 and the relative error of C n is not minimum in Fig. 5. Thus, we adopt γ 1 as a more appropriate indicator than γ 0 to determine the best-fit value of C n , including C 0 .  \nComparing the panels for N = 0 , · · · , 7 in Fig. 6, we note that the minimum of γ n for each n decreases as we increase N , but more or less saturates when we superpose three or four modes. Therefore, we expect that, in order to extract C 0 close to the fiducial values, it would be reasonable to use the fitting function superposing four modes, and to read off the best-fit values by the fitting starting from t i where γ 1 is minimum.  \nAlso, note that γ n for n ≥ 1 asymptotically approaches to almost constant for the late start time of the fit. This originates from the divergence of A n due to the overfit, which simply compensates the exponential growth of ψ n = e -iω n t . The compensation implies C n ∼ e iω n t i , for which γ n = | d ln C n /dt i | remains an almost constant. Indeed, the asymptotic value of γ n for n ≥ 1 roughly co-  \nFIG. 6. Rate of change γ n for the fit of the mock waveform Ψ [0 , 7] by the fitting function ψ fit [0 ,N ] .  \n<!-- image -->  \n-  \n-  \n-  \n-  \nFIG. 5. Relative error of C n for the fit of the mock waveform Ψ [0 , 7] by the fitting function ψ fit [0 ,N ] .  \n<!-- image -->  \n-  \n-  \n-  \n-  \nincides with the damping rate | Im( ω n ) | . Therefore, we can use γ n as an indicator for both of the stable fit and overfit.  \nIn summary, we conclude that we can improve the conventional fit and extract the best-fit value of C n close to the fiducial values by the fitting starting from t i when  \nγ n of the next longest-lived mode takes the minimum value. In particular, the fit of the longest-lived mode is most stable and gives the best-fit value close to the fiducial value. To extract the longest-lived mode and lower overtones in a stable manner from the mock waveform Ψ [0 , 7] , we should superpose four damped sinusoids as the  \nFigures 8 and 9 respectively show the amplitude A n and phase ϕ n for the fit of the mock waveform Ψ [1 , 7] by the fitting function ψ fit [1 ,N ] as a function of the start time of the fit t i . Compared to the case of the mock waveform Ψ [0 , 7] in Figs. 3 and 4, we can see the same qualitative behavior. In particular, comparing the N = 3 panels in Figs. 3 and 4 and the N = 4 panels in Figs. 8 and 9, A 2 and ϕ 2 start to diverge at ( t i -t peak ) /M ≃ 40 in Figs. 3 and 4, but in Figs. 8 and 9, A 2 and ϕ 2 do not diverge up to ( t i -t peak ) /M ≃ 50 and the plateau is extended. The same improvement applies to the plateau of A 3 and ϕ 3 . We see that the subtraction process stabilizes the fit of overtone coefficients.  \n<!-- image -->  \n-  \nFIG. 7. Mismatch M between the mock waveform Ψ [1 , 7] and the fitting function ψ fit [1 ,N ] .  \nfitting function. Superposition of five or more modes do not improve the relative error significantly. Rather, increasing the number of modes used in the fitting function sometimes destabilizes the fitting of the lower overtones. We expect that this setup captures the qualitative behavior of a more realistic situation, where one tries to fit the actual ringdown signal composed of a superposition of an infinite number of QNMs, whose fiducial values are unknown, by a superposition of the finite number of damped oscillations.', 'B. Iterative fit and subtractions': 'So far, we consider the conventional way of fitting, i.e., the fit of the mock waveform Ψ [0 , 7] given in Eq. (8) by the fitting function (1). In particular, we can extract the longest-lived mode close to the fiducial value in a stable manner. However, we found that the plateau for the higher overtone with respect to the start time of the fit is short and hence the fitting of overtones is still not robust. This trend is consistent with the previous work [38]. Usually, the fundamental mode n = 0 is the longest-lived mode, and hence the extraction of overtones that have shorter damping time is challenging.  \nSuppose that the first overtone n = 1 was the longestlived mode of the mock waveform. Then, it is natural to expect that we would be able to extract the first overtone stably. Therefore, if we subtract the best-fit fundamental mode from the original waveform, the first overtone becomes the longest-lived mode and we expect that we can extract the first overtone more stably. This is the main idea of the iterative extraction method which we propose in the present paper.  \nIn this section, from the above point of view, we consider the fitting of a mock waveform, where the fundamental mode is subtracted, and examine how the ex-  \ntraction of overtones can be improved. In practice, we cannot extract the longest-lived mode without numerical errors. Before performing a more realistic subtraction of the best-fit longest-lived mode in the subsequent sections, in this section, we assume an ideal subtraction and consider a mock waveform Ψ [ K, 7] composed of a superposition of the modes from n = K to 7 with K ≥ 1.  \nFirst, let us consider a situation that we subtract the fundamental mode from the waveform Ψ [0 , 7] . Assuming an ideal subtraction without errors, we consider the mock waveform Ψ [1 , 7] without the fundamental mode, which is shown by the orange curve in Fig. 1.  \nWe then perform the fit of the mock waveform Ψ [1 , 7] by the fitting function ψ fit [1 ,N ] . Figures 7 shows the mismatch M between Ψ [1 , 7] and ψ fit [1 ,N ] as a function of the start time of the fit t i . Compared to the case of the mock waveform Ψ [0 , 7] in Fig. 2, we can see some common features. Specifically, the mismatch does not reduce until we move the start time of the fit t i sufficiently late if we fit the mock waveform by a single damped sinusoid, and it reduces for the earlier t i if we increase the number of overtones included in the fitting function. As expected, the first overtone in Fig. 7 takes over typical features of the longest-lived mode in Fig. 2.  \nLet us check the relative error of C n and the rate of change γ n , which are shown in Figs. 10 and 11 as a function of t i , respectively. Again, we can see qualitatively same behavior as in Figs. 5 and 6 for the case of the mock waveform Ψ [0 , 7] . In particular, γ 1 in Fig. 11 now has a global minimum, which is located at a different t i of the minimum for other modes γ n ≥ 2 . This precisely matches the behavior of γ 0 in Fig. 6. Further, we see that our indicator for the best-fit values works well: We use t i when γ n for the next-longest-lived mode, which is γ 2 for the present case, is minimum. It matches the start time of the fit t i when the relative error for C n is minimum.  \nTherefore, as we expected, if the longest-lived mode is completely subtracted, the original next-longest-lived mode takes over the role of the longest-lived mode. Since we can extract the longest-lived mode most stably, this process improves the stability of the extraction of the original next-longest-lived mode. Indeed, the plateau of A n and ϕ n is extended, and the minimum of the relative error of C n is improved. Comparing the fit by the four  \nFIG. 9. Phase ϕ n for the fit of the mock waveform Ψ [1 , 7] by the fitting function ψ fit [1 ,N ] .  \n<!-- image -->  \n-  \n-  \n-  \nFIG. 8. Amplitude A n for the fit of the mock waveform Ψ [1 , 7] by the fitting function ψ fit [1 ,N ] .  \n<!-- image -->  \n-  \n-  \n-  \nlowest modes, i.e., the fit of Ψ [0 , 7] by ψ fit [0 , 3] and the fit of Ψ [1 , 7] by ψ fit [1 , 4] , the relative error of C 1 at the best-fit value improves from 6 . 8 × 10 -6 to 4 . 1 × 10 -7 .  \nWe expect that, by iterating the subtraction of the longest-lived mode, we can improve the fit of higher overtones step by step. Namely, the first subtraction of n = 0  \nmode makes n = 1 mode the longest-lived mode, and the second subtraction of n = 1 mode makes n = 2 mode the longest-lived mode, and so on. We then adopt the bestfit value when the mode of interest becomes the longestlived mode. To examine the efficiency of the fit based on the above procedure, assuming the ideal subtraction,  \n<!-- image -->  \n-  \n-  \n-  \nFIG. 10. Relative error of C n for the fit of the mock waveform Ψ [1 , 7] by the fitting function ψ fit [1 ,N ] .FIG. 11. Rate of change γ n for the fit of the mock waveform Ψ [1 , 7] by the fitting function ψ fit [1 ,N ] .  \n<!-- image -->  \nwe investigate the fit of Ψ [ K, 7] by ψ fit [ K,N ] subsequently for K = 0 , · · · , 4. We present the results in Fig. 12 to see how the fit of n = 1 , · · · , 4 modes is improved by the subtraction procedure. The first row of Fig. 12 shows the amplitude A 1 , the phase ϕ 1 , and relative error of C 1 as a function of the start time of the fit t i for the  \niterative fit of the mock waveform Ψ [ K, 7] by the fitting function ψ fit [ K,K +3] . Specifically, in the first row of Fig. 12, the dot-dashed curve corresponds to the fit of Ψ [0 , 7] by ψ fit [0 , 3] , and the solid curve corresponds to the fit of Ψ [1 , 7] by ψ fit [1 , 4] . We can see that as we subtract the longestlived mode, the fitting improves; the plateau in A 1 and  \n1  \nA  \n2  \nA  \n4  \nA  \n10  \n10  \n10  \n10  \n10  \n10  \n10  \n10  \n3  \n2  \n1  \n0  \n3  \n2  \n1  \n0  \nFiducial value  \n[  \nK,N  \n] = [0  \n,  \n3]  \n[  \nK,N  \n] = [1  \n,  \n4]  \n0  \n(  \nti  \nFiducial value  \n[  \nK,N  \n] = [0  \n,  \n3]  \n[  \nK,N  \n] = [1  \n,  \n4]  \n[  \nK,N  \n] = [2  \n,  \n5]  \n0  \n20  \n40  \n60  \n-  \n(  \nti  \nt  \nt  \npeak  \npeak  \n)  \n/ M  \n)  \n/ M  \nFIG. 12. The improvement of the amplitude A n (left column), phase ϕ n (middle column), and relative error of C n (right column) from n = 1 mode (first row) to n = 4 mode (fourth row) for the iterative fit of the mock waveform Ψ [ K, 7] by the fitting function ψ fit [ K,N ] with N = K +3.  \n<!-- image -->  \n10  \n10  \n10  \n10  \n10  \n4  \n3  \n2  \n1  \n0  \n-  \nFiducial value  \n[  \nK,N  \n] = [1  \n,  \n4]  \n[  \nK,N  \n] = [2  \n,  \n5]  \n[  \nK,N  \n] = [3  \n,  \n6]  \n[  \nK,N  \n] = [4  \n,  \n7]  \n0  \n20  \n40  \n60  \n-  \n(  \nti  \nt  \npeak  \n)  \n/ M  \n(  \nti  \n(  \nti  \n(  \nti  \n(  \nti  \n0  \n20  \n40  \n60  \n0  \n20  \n40  \n60  \n0  \n20  \n40  \n60  \n0  \n20  \n40  \n60  \n-  \nt  \nt  \npeak  \npeak  \nt  \nt  \n)  \n/ M  \n)  \n/ M  \npeak  \npeak  \n)  \n/ M  \n)  \n/ M  \n(  \nti  \n(  \nti  \n(  \nti  \n(  \nti  \n0  \n20  \n40  \n60  \n0  \n20  \n40  \n60  \n0  \n20  \n40  \n60  \n0  \n20  \n40  \n60  \n-  \nt  \nt  \nt  \nt  \npeak  \npeak  \npeak  \npeak  \n)  \n/ M  \n)  \n/ M  \n)  \n/ M  \n)  \n/ M  \nϕ 1 becomes longer and the relative error from the fiducial value becomes smaller. This improvement due to the iterative fit and subtraction procedure is more manifest for n = 2 , 3, and 4, which are presented in the second, third, and fourth row of Fig. 12, respectively. For instance, let us focus on the case of n = 3 in the third row of Fig. 12. It may be difficult to extract the correct value from the fit of Ψ [0 , 7] by ψ fit [0 , 3] shown by the dotted curve as the plateau is marginal and the fitting is not so sta-  \n. However, iterating the (ideal) subtractions, we can obtain a sufficiently long plateau for the fit of Ψ [3 , 7] by ψ fit [3 , 6] shown by the solid curve, which is also consistent with the fiducial value C fid 3 .  \nLet us compare the amplitude A n and phase ϕ n obtained by the conventional fit and iterative fit of the mock waveform Ψ [0 , 7] . Here we provide two ways of comparison. First, in Figs. 13 and 14, we compare the results of the conventional fit with the fitting function ψ fit [0 , 7] and  \n-  \n-  \n-  \n-  \n-  \n-  \n20  \n40  \n60  \n-  \nπ  \n/  \n1  \nφ  \nπ  \n/  \n2  \nφ  \nπ  \n/  \n3  \nφ  \nπ  \n/  \n4  \nφ  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \n0  \n.  \n2  \n0  \n.  \n4  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n6  \n0  \n.  \n8  \n1  \n.  \n0  \n1  \n.  \n2  \n1  \n.  \n4  \n1  \n.  \n6  \n0  \n.  \n75  \n0  \n.  \n50  \n0  \n.  \n25  \n0  \n.  \n00  \n0  \n.  \n25  \n0  \n.  \n50  \n0  \n.  \n25  \n0  \n.  \n50  \n0  \n.  \n75  \n1  \n.  \n00  \n1  \n.  \n25  \n1  \n.  \n50  \n|  \n1  \n-  \nfiducial1  \n/C  \n1  \nC  \n10  \n10  \n10  \n-  \n10  \n-  \n|  \n10  \n-  \n10  \n10  \n|  \n1  \n-  \nfiducial2  \n/C  \n2  \nC  \n10  \n10  \n-  \n|  \n10  \n-  \n10  \n10  \n10  \n|  \n1  \n-  \nfiducial3  \n/C  \n3  \nC  \n|  \n10  \n10  \n-  \n10  \n-  \n10  \n10  \n10  \n|  \n1  \n-  \nfiducial4  \n/C  \n4  \nC  \n|  \n10  \n-  \n3  \n1  \n1  \n3  \n5  \n4  \n2  \n0  \n2  \n4  \n7  \n5  \n3  \n1  \n1  \n3  \n7  \n4  \n1  \n2  \nFIG. 14. Comparison of the phase ϕ n obtained by fitting the mock waveform Ψ [0 , 7] by the conventional method (top) with the fitting function ψ fit [0 , 7] and the iterative method (bottom) with the fitting function ψ fit [ K,N ] with K = 0 , · · · , 4 and N = K +3. The dotted lines represent the fiducial values.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 13. Comparison of the amplitude A n obtained by fitting the mock waveform Ψ [0 , 7] by the conventional method (top) with the fitting function ψ fit [0 , 7] and the iterative method (bottom) with the fitting function ψ fit [ K,N ] with K = 0 , · · · , 4 and N = K +3. The dotted lines represent the fiducial values.  \nthe results of the iterative fit with the fitting function ψ fit [ K,N ] with K = 0 , · · · , 4 and N = K + 3. This is because, as we stressed above, the longest-lived mode is more stably extracted when we fit it by a superposition of several modes. As a representative number of the superposition, we adopt the four modes. Consequently, while both methods use n = 0 , · · · , 7 overtones for the fitting function, we obtain C n with n = 0 , · · · , 7 for the conventional method and C n with n = 0 , · · · , 4 for the iterative method. The top panels show the result of the conventional fitting, which is the same as the panels for N = 7 in Figs. 3 and 4. We see that while the fundamental mode can be extracted in a stable manner, the overtones have shorter plateau and are sensitive to the choice of the start time of the fit. The bottom panels show the result of the fitting based on the iterative fitting method. Namely, we present the result of the subsequent fit of the mock waveform Ψ [ n, 7] by the fitting function ψ fit [ n,n +3] for  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nn = 0 , · · · , 4, where the n -th mode is extracted when it is the longest-lived mode. We see that for both of A n and ϕ n , the plateau consistent with the fiducial value is extended for the overtones. Thus, the iterative fitting method improves the stability of the extraction of the higher overtones.  \nAs an alternative way of comparison of the conventional and iterative fitting methods, in Fig. 15, we compare the results of the conventional fit with the fitting function ψ fit [0 , 4] and the results of the iterative fit with the fitting function ψ fit [ K, 4] with K = 0 , · · · , 4. For both methods, we use n = 0 , · · · , 4 overtones for the fitting function, and obtain the same number of C n . We see that the iterative fitting method provides a longer plateau in A n and ϕ n for overtones, and hence improves the stability of their extraction.  \nIn summary, through the analysis of a simple mock waveform of a superposition of damped sinusoids, we see  \nFIG. 16. Mock waveform Ψ ( c ) [0 , 7] (blue) with a constant given in Eq. (11) and the waveform Ψ ( c ) [0 , 7] -C iter 0 ψ 0 (orange) after the first subtraction, which is m = 1 case of Eq. (12).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 15. Comparison of the amplitude A n (top) and phase ϕ n (bottom) obtained by fitting the mock waveform Ψ [0 , 7] by the conventional method (dashed) with the fitting function ψ fit [0 , 4] and the iterative method (solid) with the fitting function ψ fit [ K, 4] with K = 0 , · · · , 4. The dotted lines represent the fiducial values.  \nthat the longest-lived mode can be extracted in a stable manner, while a stable extraction of the overtones is in general difficult. We can improve the fit of the overtones by subtracting the longest-lived mode. After the subtraction, the next-longest-lived mode in the original waveform plays the role of the longest-lived mode, and its fit is stabilized. By iteratively subtracting the longest-lived modes, we can realize a stable extraction of the overtones. The subtraction of the longest-lived mode does not affect the contribution from other modes, and we can directly extract the fiducial value from the mock waveform. However, we should be bear in mind that we treated an ideal mock waveform without noise nor powerlaw tail and assumed an ideal subtraction without errors. In § IV, we consider another mock waveform taking into account these points.  \n<!-- image -->  \n-', 'IV. DAMPED SINUSOIDS WITH CONSTANT': "In § III, we analyze the mock waveform of a superposition of pure damped sinusoids assuming the ideal subtraction. By subtracting the longest-lived mode, the next-longest-lived mode can be extracted in a more stable manner. By iteratively subtracting the longest-lived mode from the waveform, we see that we can extract higher overtones such as n = 3 mode.  \nIn this section, we consider a bit more realistic analysis with increasing complexity. We consider a mock waveform with a constant added and perform the actual subtraction of the longest-lived mode by using the best-fit values. We shall also discuss how the stability of the fit will change for the data with different sampling rates at the end of this section. We shall see that the iterative fitting method is still efficient for waveform data with low sampling rates.  \nWe add a constant to the mock waveform to mimic a more realistic situation. For the numerical relativity waveform, which is the main focus on the present paper, it is known that the waveform data mostly approaches a small constant at the late time due to numerical errors. The observed ringdown waveform is also contaminated by noises, which may not be removed completely. Even if noises were completely removed, we would not be able to observe damped sinusoids very long time, since the ringdown waveform would be dominated by the powerlaw tail at the late time. Thus, it would be inevitable that these effects hide the damped sinusoids at the late time below a certain order of threshold.  \nTaking into account the above facts, we consider the  \nFIG. 18. Mismatch M (left top), amplitude A n (right top), phase ϕ n (left middle), relative error of C n (right middle), and rate of change γ n (bottom) for the fit of the mock waveform Ψ ( c ) [0 , 7] -C iter 0 ψ 0 by the fitting function ψ fit [1 ,N ] with N = 4.  \n<!-- image -->  \n-  \n<!-- image -->  \n<!-- image -->  \n-  \n<!-- image -->  \n<!-- image -->  \n-  \nFIG. 17. Mismatch M (left top), amplitude A n (right top), phase ϕ n (left middle), relative error of C n (right middle), and rate of change γ n (bottom) for the fit of the mock waveform Ψ ( c ) [0 , 7] by the fitting function ψ fit [0 ,N ] with N = 3.  \nfollowing mock waveform for t ≥ t peak :  \nΨ ( c ) [0 , 7] ( t ) = 7 ∑ n =0 a n e -i ( ω n t -θ n ) + c, (11)  \nwhere we set a n and θ n the same fiducial value as in (8), which is listed in Table I. Here, we introduced a small complex constant c = -(0 . 8 + 1 . 8 i ) × 10 -5 as a 'noise', whose value is the same as the constant appearing in the late time waveform of the numerical simulation SXS:BBH:0305. For the region t < t peak , to make the mock waveform closer to those of the numerical simulation, we use the waveform of SXS:BBH:0305.  \nWhile the actual noise and tail are not constant, the qualitative behavior that the damped sinusoids become subdominant at the late time is well captured by this mock waveform. Actually, the conventional fit of a mock waveform with/without tail and numerical relativity waveform is extensively studied in [38], and we shall see that our results of the conventional fit of Ψ ( c ) [0 , 7] are  \n-  \n-  \nπ  \n/  \nn  \nφ  \n10  \n-  \n10  \n-  \nM  \n10  \n-  \n0  \n.  \n5  \n0  \n.  \n0  \n0  \n.  \n5  \n1  \n.  \n0  \n1  \n.  \n5  \n2  \n.  \n0  \n-  \n-  \n-  \n-  \n2  \n4  \n6  \n-  \nN  \n= 1  \nN  \n= 2  \nN  \n= 3  \nN  \n= 4  \nN  \n= 5  \nN  \n= 6  \nN  \n= 7  \n40  \n-  \n20  \n0  \n20  \n40  \n(  \nti  \n(  \nti  \npeak  \n-  \nt  \nN  \n= 4  \n-  \nt  \n0  \n20  \n40  \npeak  \n)  \n/ M  \n10  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n0  \n1  \n2  \n3  \n4  \nn  \nγ  \n)  \n/ M  \n0  \n20  \n40  \npeak  \n(  \nti  \nN  \n= 4  \n-  \nt  \n)  \n/ M  \nqualitatively consistent with the previous results of the fit of the mock waveform with tail.  \nWe present the mock waveform Ψ ( c ) [0 , 7] ( t -t peak ) by a blue curve in Fig. 16. The damped sinusoids begin to be hidden for ( t -t peak ) /M ≳ 125 as they reach to the order of the added constant ≃ 10 -5 . Compared to the case of the mock waveform Ψ [0 , 7] without noise in Fig. 1, due to the constant, we have less data available for the fitting, and hence expect that the efficiency of the fit becomes worse.  \nWe perform the fitting of the mock waveform Ψ ( c ) [0 , 7] ( t -t peak ) by the fitting function ψ fit [0 ,N ] ( t -t peak ), and present the results in Fig. 17. For the amplitude, relative error, and the rate of change in Fig. 17, we display only N = 3 as the representative case with the four modes in the fitting function as we discussed in § III. The all cases N = 0 , · · · , 7 are supplemented in Figs. 43 in Appendix B.  \nThe first panel of Fig. 17 shows the mismatch. Unlike the fit of the mock waveform without constant shown in  \n|  \n1  \n-  \nfid  \nn  \n/C  \nn  \n10  \n-  \nC  \n|  \n10  \n-  \n0  \n20  \n40  \npeak  \n)  \n/ M  \n(  \nti  \n<!-- image -->  \n10  \n10  \n3  \n1  \n1  \n3  \n-  \nN  \n= 4  \n-  \nt  \nFig. 2, we see that the mismatch in Fig. 17 is bounded from below due to the existence of the noise in the mock waveform. The lower bound increases for the fitting starting from later time. This behavior is typical for the fit of the waveform with noise or tail, originating from the larger ratio of the noise or tail to the damping waveform at later time. In the case of the mock waveform without noise in Fig. 2, we obtain smaller value of mismatch at earlier t i reaching the order of the numerical errors. Such behavior of the mismatch does not occur for a more realistic case with noises, and the fitting becomes more difficult in general.  \nThe second and third panels of Fig. 17 show the amplitude A n and phase ϕ n as a function of the start time of the fit t i for the fit of the mock waveform Ψ ( c ) [0 , 7] by the fitting function ψ fit [0 ,N ] . Compared to the fit of Ψ [0 , 7] without noise in Figs. 3 and 4, the plateau is shortened due to the existence of noise. For the N = 3 case, the length of the plateau for A 2 and ϕ 2 differs by about a factor of 2 with and without noise. The plateau is also shortened for other N cases shown in Fig. 43. For instance, for the N = 1 case, A 1 and ϕ 1 remains almost constant until ( t i -t peak ) /M ≃ 50 in Figs. 3 and 4 but ( t i -t peak ) /M ≃ 40 in Fig. 43. The divergence tends to begin earlier with higher overtones. As pointed out in [38], the introduction of noise shortens the plateau and reduces the stability of the fit.  \nWe present the relative error of C n in the fourth panel of Fig. 17 for N = 3 and Fig. 43 for the all cases N = 0 , · · · , 7. For the most modes, the behavior of the relative error of C n is similar to that of Fig. 5. For instance, t i at which the relative error reaches a minimum is roughly the same regardless of N , and the larger N is, the earlier t i at which the relative error reaches a minimum becomes in the early ringdown period. However, the differences show up in the fundamental mode for the fitting starting from late time. For the N = 2 case in Fig. 5 for the mock waveform without noise, the relative error in the fundamental mode has the minimum at t i that were significantly different from the minima in the other modes. In contrast, in Fig. 43 for the mock waveform with noise, the relative error has minima at similar t i in all modes. This originates from the fact that the waveform in the late time is contaminated by the noise. The tendency of the lower bound for the fitting starting from late time t i is also consistent with the behavior of the mismatch in Fig. 17.  \nThe bottom panels of Fig. 17 and Fig. 43 show the rate of change γ n , where we see several differences from Fig. 6 for the case without noise. The minima of γ n appear more sharply when the noise is present than when the noise is absent. Also, a lower bound shows up at the late time, which is similar to the mismatch. Further, there is no exceptional behavior of the minimum of γ 0 , and all γ n take the minima in the same start time of the fit.  \nEven with the noise in the mock waveform, we confirm that the rate of change γ n for the next-longest-lived mode works well as an indicator of the optimal start time of  \nthe fit with the minimum relative error. In Fig. 17, we see that there exists a clear correspondence between the minimum of the relative error of C n and the minimum of the rate of change γ n for the next-longest-lived mode. The correspondence also holds for other N cases shown in Fig. 43. Compared to the amplitude A n and phase ϕ n in Fig. 43, we see that in each N , the minimum of γ n is located just before the highest overtone A n and ϕ n begin to diverge. It would be reasonable to consider that the minimum of γ n is located before the contamination of the fit by noise. Thus, it is still efficient to use the minimum of γ n as an indicator for the survey of the best-fit value even when the noise is present. Incidentally, there exist relatively sharp peaks just after γ n reaches a minimum, which are caused by the passage of C n near the origin on the complex plane, as confirmed by the dips in A n and the sudden changes in ϕ n .  \nNext, we perform the subtraction of the longest-lived mode. While in § III B we assumed the ideal subtraction and used Ψ [1 , 7] , in this section we subtract the longestlived mode from the mock waveform by using the bestfit value of C n . Namely, the waveform after the m -th subtraction reads  \nΨ ( c ) [0 , 7] -m -1 ∑ n =0 C iter n ψ n . (12)  \nHere, C iter n is the best-fit value of C n obtained by the fitting of the waveform after the n -th subtraction, where n -th overtone becomes the effective longest-lived mode.  \nThe waveform after the first subtraction is shown by orange curve in Fig. 16, where we can see that the damped sinusoid corresponding to the first overtone. Furthermore, we can clearly see that the constant mimicking the noise and tail becomes dominant at the late time. Compared to the mock waveform Ψ [1 , 7] in § III assuming the absence of the noise/tail and the ideal subtraction, it is clear that the noise shortens the region where we can observe the damped sinusoids. We cannot fit the first overtone for the region where it is subdominant compared to the noise, and hence we have less data available for the fit. Specifically, for the original mock waveform Ψ ( c ) [0 , 7] before the subtraction we can observe the damped sinusoids up to ( t -t peak ) /M ≃ 140, whereas for the waveform Ψ ( c ) [0 , 7] -C iter 0 ψ 0 we can only see the damped sinusoids up to ( t -t peak ) /M ≃ 50. This is because the first overtone appeared after the subtraction has shorter damping time than the fundamental mode and hence becomes comparable to the noise more quickly. Thus, for the waveform with the noise, iterating the subtraction of the longest-lived mode, we have less data available for the fit of overtones. Consequently, for the iterative fitting method, we gradually decrease the end time of the fit t e at each step. Specifically, we set ( t e -t peak ) /M = 50 for the fit of the waveform Ψ ( c ) [0 , 7] -C iter 0 ψ 0 .  \nTo the mock waveform Ψ ( c ) [0 , 7] -C iter 0 ψ 0 obtained after the first subtraction, we apply the fitting algorithm  \nFIG. 19. The improvement of the amplitude A n (left column), phase ϕ n (middle column), and relative error of C n (right column) from n = 1 mode (first row) to n = 4 mode (fourth row) for the iterative fit of the mock waveform given in (12) by the fitting function ψ fit [ K,N ] with N = K +3.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nwith the fitting function ψ fit [1 ,N ] , and present the results in Fig. 18. Again, for the amplitude, phase, relative error, and the rate of change, we display N = 4 case in Fig. 18 and the all cases N = 1 , · · · , 7 in Fig. 44 in Appendix B.  \nThe first panel of Fig. 18 shows the mismatch. Compared to the fit of the original waveform before the subtraction of the fundamental mode in Fig. 17, there are no significant difference for the range 0 < ( t i -t peak ) /M ≲ 20. However, the tilt of the lower bound of the mis-  \ntch becomes larger around ( t i -t peak ) /M ≃ 20. This is caused by the shortening of the interval available for the fit due to the noise appeared after the subtraction of the fundamental mode. The position of the start time of the fit t i where the mismatch takes a minimum value is almost unchanged.  \nWe present the amplitude A n and phase ϕ n in the second and third panels of Fig. 18 for N = 4 and Fig. 44 for the all cases N = 1 , · · · , 7. Compared to the fit of  \nFIG. 21. Comparison of the phase ϕ n obtained by fitting the mock waveform Ψ ( c ) [0 , 7] by the conventional method (top) with the fitting function ψ fit [0 , 7] and the iterative method (bottom) with the fitting function ψ fit [ K,N ] with K = 0 , · · · , 4 and N = K +3. The dotted lines represent the fiducial values.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 20. Comparison of the amplitude A n obtained by fitting the mock waveform Ψ ( c ) [0 , 7] by the conventional method (top) with the fitting function ψ fit [0 , 7] and the iterative method (bottom) with the fitting function ψ fit [ K,N ] with K = 0 , · · · , 4 and N = K +3. The dotted lines represent the fiducial values.  \nthe original mock waveform before the subtraction of the fundamental mode, we see that the plateau of the first overtone is slightly extended. As for the second overtone, the plateau is even more extended. However, the improvement is not as large as that in the case of the waveform without noise, i.e., from Figs. 3 and 4 to Figs. 8 and 9. This is because at the late time we are actually trying to fit the constant noise, which is dominant over the fiducial damped sinusoids, by the fitting function of a superposition of damped sinusoids. Such inappropriate fitting inevitably leads to exponential divergence of the coefficients determined by the damping rate of each mode. While we observed this phenomenon for the fit of the higher overtones of the mock waveform without noise, in the case of the mock waveform with noise, it also occurs for the fit of the longest-lived mode, which prevents a significant extension of the plateau.  \nWe present the relative error of C n and the rate of  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nchange γ n in the fourth and fifth panels of Fig. 18 for N = 4 and Fig. 44 for the all cases N = 1 , · · · , 7, respectively. As mentioned above, we determine the best-fit value by using the start time of the fit where γ n for the next-longest-lived mode takes the minimum. Comparing the fit of the mock waveform before and after the subtraction of the fundamental mode, the relative errors of the best-fit value are more or less the same for the longestlived mode, the next-longest-lived mode, and so on. In other words, focusing on the same n -th mode, the relative error reduces by virtue of the subtraction. Also, we can see that the plateau is extended. For n = 1 mode, the region where γ 1 < 10 -2 is clearly extended from Fig. 17 to Fig. 18. For the relative error of C 1 , the dip shows up in Fig. 18 at late time, but the relative errors of other modes do not take the minimum there. We also observed such peculiar dip in the analysis of the mock waveform without noise in § III. Even with the noise, our criteria  \nFIG. 22. Comparison of the amplitude A n (top) and phase ϕ n (bottom) obtained by fitting the mock waveform Ψ ( c ) [0 , 7] by the conventional method (dashed) with the fitting function ψ fit [0 , 4] and the iterative method (solid) with the fitting function ψ fit [ K, 4] with K = 0 , · · · , 4. The dotted lines represent the fiducial values.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nstill work and we can determine the optimal start time of the fit by the minimum of γ n for the next-longest-lived mode.  \nHaving confirmed that our strategy works well, we iterate the subtraction of the longest-lived mode and obtained the best-fit value for the coefficient C n for n = 0 , · · · , 4. We highlight the improvement of the fit of each mode by virtue of the subtraction in Fig. 19, which should be contrasted with Fig. 12 for the case assuming the ideal subtraction and the absence of noise. We can see that the plateau is extended in Fig. 19, but the improvement is less significant than Fig. 12. Clearly, this is due to the existence of the constant, which places the upper bound for the plateau.  \nSuch behavior of the plateau can also be observed in Figs. 20-22, where we compare the conventional and iterative fit in the two ways described in § III B. Although the  \niterative fit provided very long plateau for the case assuming the ideal subtraction and the absence of noise, in a more realistic case of the mock waveform with constant and performing numerical subtraction, the extension of the plateau by the iterative fit in Figs. 20-22 is limited compared to Fig. 13-15. Having said that, we can see that the iterative fitting method helps us to improve the stability of the fit of the overtones for n ≥ 2 compared to the conventional fit. Apparently, while the plateau is short, the conventional fit in Fig. 20 seems to succeed the extraction of the fiducial amplitude for the overtones. This is because in this case we are fitting a superposition of the seven damped sinusoids by a superposition of the same number of damped sinusoids. Indeed, in Fig. 22 and also in Figs. 43 and 44 in Appendix B, there is no plateau for overtones in the conventional method when the number of modes included in the fitting function is fewer. As we discussed, in a realistic case, we need to fit a superposition of infinite number of damped sinusoids by a superposition of finite number of damped sinusoids. Also, the nonlinear effects become dominant at early time, and prevent the stable extraction of overtones. Hence, the overtones in the conventional fit are not likely to possess even short plateau. We shall see in § V that this is indeed the case for the fit of the numerical relativity waveform (see e.g., Figs. 29 and 31).  \nIn Fig. 22 for the alternative way of comparison, the plateau is also extended by the iterative fit. However, from Fig. 22 we note that while the plateau in the phase ϕ n is consistent with the fiducial value, the plateau in the amplitude A n is off from the fiducial value for higher overtones, with the difference in A 4 reaching approximately a factor of 2. On the other hand, in Fig. 20, A 4 extracted by the iterative method is close to the fiducial value. This is because we restrict ourselves to use the fitting function consisting of the modes up to fourth overtone in Fig. 22, whereas we use the modes up to seventh overtone in Fig. 20.  \nFinally, let us discuss the stability of the fit with respect to the sampling rate of the waveform data. So far, we use a discrete data ( t a , Ψ ( c ) [0 , 7] ( t a )) of the mock waveform with the same sampling times as BBH:SXS:0305, which we denote waveform data 0. The time data spacing is not exactly constant but remains ∆ t a /M ≈ O (10 -1 ) for the ringdown phase. In general, even with the same waveform, sampling rate affects the fitting, i.e., different sampling rates can yield different fitting results [31]. Based on the waveform data 0, ( t a , Ψ ( c ) [0 , 7] ( t a )) with a = 0 , · · · , N samp , we generate the following three sets of data with reduced sampling rates:  \n- 1. Waveform data 1: data with the half number of sampling points ⌊ N samp / 2 ⌋ , using even-numbered samples of data 0.\n- 2. Waveform data 2: data with the quarter number of sampling points ⌊ N samp / 4 ⌋ , using even-numbered samples of data 1.  \nn  \nA  \nn  \nA  \n10  \n10  \n10  \n10  \n10  \n-  \n3  \n2  \n1  \n0  \n1  \n-  \n10  \n10  \n-  \n-  \n5  \n0  \n5  \n10  \n15  \n20  \n-  \n(  \nt  \ni  \nt  \npeak  \n)  \n/ M  \nFIG. 23. Same plots as in Fig. 20 for the mock waveform data Ψ ( c ) [0 , 7] with reduced number of sampling points and its interpolation. Conventional fit (left column) and iterative fit (right column) for data 1 (top row) with half sampling points, data 2 (middle row) with quarter sampling points, and data 3 (bottom row) generated by the interpolation of the data 2.  \n<!-- image -->  \n10  \n10  \n10  \n10  \n10  \n-  \n3  \n2  \n1  \n0  \n1  \n-  \n-  \n5  \n0  \n5  \n10  \n15  \n20  \n-  \n(  \nt  \ni  \nt  \npeak  \n)  \n/ M  \n(  \nt  \n(  \nt  \n(  \nt  \ni  \ni  \ni  \n-  \n-  \n-  \nn  \n= 0  \nn  \n= 1  \nn  \n= 2  \nn  \n= 3  \nn  \n= 4  \n-  \nn  \n= 0  \nn  \n= 1  \nn  \n= 2  \nn  \n= 3  \nn  \n= 4  \n-  \nn  \n= 0  \nn  \n= 1  \nn  \n= 2  \nn  \n= 3  \nn  \n= 4  \n-  \n5  \n0  \n5  \n10  \n15  \n20  \nt  \npeak  \n)  \n/ M  \n)  \n/ M  \n)  \n/ M  \n- 3. Waveform data 3: data with the number of sampling points 2 ×⌊ N samp / 4 ⌋-1 doubled by the cubic spline interpolation of data 2.  \nIn Figs. 23 and 24, we show the amplitude A n and  \nphase ϕ n obtained by the conventional (left column) and iterative fitting (right column) between the above three different sampling rates of the mock waveform data. Compared to the results for the data 0 in Figs. 20 and  \n5  \n0  \n5  \n10  \n15  \n20  \npeak  \n5  \n0  \n5  \n10  \n15  \n20  \nt  \nt  \npeak  \nn  \n= 0  \nn  \n= 1  \nn  \n= 2  \nn  \n= 3  \nn  \n= 4  \nn  \n= 5  \nn  \n= 6  \nn  \n= 7  \nn  \n= 0  \nn  \n= 1  \nn  \n= 2  \nn  \n= 3  \nn  \n= 4  \nn  \n= 5  \nn  \n= 6  \nn  \n= 7  \nn  \nA  \nn  \nA  \nn  \nA  \n10  \n10  \n10  \n10  \n10  \n-  \n3  \n2  \n1  \n0  \n1  \n-  \n10  \n10  \n10  \n10  \n10  \n-  \n3  \n2  \n1  \n0  \n1  \n-  \n10  \n10  \n10  \n10  \n10  \n-  \n3  \n2  \n1  \n0  \n1  \n-  \n10  \n10  \n10  \nFIG. 24. Same plots as in Fig. 21 for the mock waveform data Ψ ( c ) [0 , 7] with reduced number of sampling points and its interpolation. Conventional fit (left column) and iterative fit (right column) for data 1 (top row) with half sampling points, data 2 (middle row) with quarter sampling points, and data 3 (bottom row) generated by the interpolation of the data 2.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n21, the fitting becomes unstable as the sampling rate decreases. Further, these figures highlight the efficiency of the iterative fitting. For the data 1 (top rows), the plateau in the conventional fit disappears, whereas the  \niterative fit still keep almost same structure as the result for the data 0. Compared to the conventional fitting, the iterative fitting is able to extract the fiducial values in a stable manner even for data with low sampling rates.  \nOf course, as the number of sampling points is reduced, the iterative fitting eventually break down as well. For the data 2 (middle rows), both fitting methods fail to extract the fiducial values. One could encounter in such a situation for the analysis of the low resolution data, so we check the fitting of the data 3, generated by the interpolation of the data 2. From the bottom rows of Fig. 23 and 24, we see that the fitting of the interpolated data can recover the result of data 1.  \nThe lessons from the mock waveform analysis in this section is that the constant, which mimics noise and/or tail in observed ringdown signal, is dominant for the late time and hence affects the fit if the start time of the fit is late. The iterative fitting method helps us to extract higher overtones, although, due to the noise, the region where the meaningful data is available becomes shorter and shorter as we iterate the subtractions. The noise places the upper bound for the plateau, and the improvement of the stability of the fit by the iterative subtractions is limited compared to the case without the noise. Up to which overtone a sufficiently long plateau can be obtained depends on the magnitude of the noise and the damping time of the overtone. While we dealt with the mock waveform with a constant in this section, these lessons are robust even for the random noises and/or power-law tail. Further, we confirm that the iterative fitting method is robust to the data with low sampling rates, compared to the conventional fit.", 'V. NUMERICAL RELATIVITY WAVEFORM': 'From the analyses of the mock waveform in § III and IV, we figure out that we can improve the stability of the fit of the overtones by iteratively subtracting the longestlived mode step by step from the ringdown waveform. It is because the longest-lived mode can be most stably fitted, and the next-longest-lived mode takes over the same feature after the subtraction.  \nIn this section, we apply the iterative fitting method to numerical relativity waveform in the SXS catalog listed in Table II and explore the fit of the ( ℓ, m ) = (2 , 2) mode of these waveforms. The Kerr QNM frequencies corresponding to each value of the remnant dimensionless spin are listed in Table VIII-XI in Appendix A. We describe the estimation of the numerical errors in Appendix C.  \nFirst, we analyze SXS:BBH:0305, which is well-known for its similarity to GW150914 [20], the first observed gravitational waves from BBH. Highlighting differences from the analysis in § III and IV, the analysis of SXS:BBH:0305 in this section helps us to understand how the nonlinearity and/or spherical-spheroidal mixing affect the results of fitting. Since previous work have mainly analyzed SXS:BBH:0305, the analysis of SXS:BBH:0305 is ideal for comparing and examining the validity of our iterative fitting method.  \nIn addition to the SXS:BBH:0305, we also investigate high spin case SXS:BBH:0158 and low spin case  \nTABLE II. Remnant dimensionless spin χ rem , large mass ratio q , and resolution level of the SXS numerical simulations [43, 44] examined in § V.  \nSXS:BBH:0156 for comparison. We choose them as the equal mass binaries having the highest/lowest remnant dimensionless spin among the SXS simulations with high resolutions. Further, we also pick up another low spin simulation, SXS:BBH:1108, which is the simulation with the parameters similar to the observed event GW190814 [47]. While the numerical errors in these simulations are larger than in SXS:BBH:0305, it is interesting to investigate how the efficiency of the ringdown fitting depends on the black hole spin.  \nIt is known that most of the Kerr QNM frequencies, except for the fifth overtone, have higher frequency and longer damping time for higher spin Kerr black holes, and they become eventually degenerate as they approach the accumulation point at the extremal limit [48-52]. On the one hand, slowly damping and fast oscillating ringdown signal would be more easily observed. On the other hand, damped sinusoids with almost degenerate frequencies and damping time would be hard to be resolved. In contrast, the opposite applies to the QNMs for lower spin black holes. They damp faster and oscillate more slowly, which would restrict the range of the time domain data available for the fitting. However, the differences between the QNM frequencies are more distinct for lower spin black holes, and may be more easily distinguished. The fitting analysis of the waveforms with various values of black hole spin would clarify these points.  \nWe denote the best-fit waveform by using the iterative fitting method as  \nψ iter [0 ,N ] = N ∑ n =0 C iter n ψ n . (13)  \nHere, ψ n = e -iω n t is the damped sinusoid with QNM frequencies ω n of ( ℓ, m ) = (2 , 2) mode of gravitational waves emitted from a Kerr black hole with the remnant dimensionless spin of each SXS simulations, and C iter n is the best-fit value coefficient obtained by each step of the iterative procedure. Namely, first, we fit the original SXS waveform Ψ SXS = ( h + -ih × ) r/M by a fitting function ψ fit [0 , 7] given in (1) with the fitting algorithm described in § II, and obtain the best-fit value C iter 0 for the fundamental mode by setting the start time of the fit as the time when the rate of change γ 1 significantly reduces. We then subtract the fundamental mode and consider Ψ SXS -ψ iter [0 , 0] , where the original first overtone is  \nFIG. 26. Mismatch M (left top), amplitude A n (right top), phase ϕ n (left bottom), and rate of change γ n (right bottom) for the fit of SXS:BBH:0305 waveform Ψ SXS -c by the fitting function ψ fit [0 ,N ] .  \n<!-- image -->  \n-  \nFIG. 25. Numerical relativity waveform SXS:BBH:0305, waveform raw data (black dashed) and waveforms after the subtraction of the numerical constant (blue) with numerical errors (gray, see Appendix C). We iteratively fit and subtract the longest-lived mode. Waveform after the subtraction of the best-fit fundamental mode (orange) and then the best-fit first overtone (green) are shown.  \nnow the longest-lived mode. Next, we fit the waveform Ψ SXS -ψ iter [0 , 0] by a fitting function the ψ fit [1 , 7] , and obtain the best-fit value C iter 1 for the first overtone. With C iter 1 , we can subtract the first overtone and consider Ψ SXS -ψ iter [0 , 1] . By iterating this procedure, we obtain the best-fit values C iter n and construct the fitted waveform ψ iter [0 ,N ] step by step.  \nTo further improve the stability of the fit, we also subtract a numerical constant c which is dominant in the late time waveform. Therefore, we insert a preliminary step for the identification and subtraction of c . After that we apply the subtraction procedure to the waveform Ψ SXS -c .', 'A. SXS:BBH:0305 - GW150914-like simulation': 'First, let us investigate the GW150914-like waveform SXS:BBH:0305. Through the iterative fitting method, we consider several waveforms step by step. Among them, we show first three waveforms in Fig. 25. The black curve depicts the raw data waveform Ψ SXS of SXS:BBH:0305 from the SXS catalog. As we can see, this waveform includes the effect of a numerical constant at late time, after ( t -t peak ) /M ≈ 130. We identify this constant value by averaging the waveform for ( t -t peak ) /M ≥ 120 after the subtraction of the fitted fundamental mode. We obtain c = -(0 . 8 + 1 . 8 i ) × 10 -5 , which is consistent with the value used in the literature [38, 53]. As we discussed in § IV, such a constant reduces the accuracy of the fit. Therefore, we consider the waveform Ψ SXS -c , in which the numerical constant is subtracted from the original  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nwaveform. The subtraction of the numerical constant was also performed in [37] and known to provide a result similar to the one with the mapping to the super rest frame [54-58].  \nWe confirm that the mismatch between Ψ SXS and ψ fit [0 ,N ] without the subtraction of the numerical constant is consistent with the results in [45], where the lower bound of the mismatch is slightly lifted, and that the difference originates from the constant [58].  \nThe waveform after the subtraction of the constant is shown by blue curve in Fig. 25. We can see an improvement that the damped sinusoid continues in the waveform after ( t -t peak ) /M = 130. We then extract the coefficient for fundamental mode C iter 0 by fitting the waveform Ψ SXS -c by the fitting function ψ fit [0 , 3] of a superposition of four QNM modes as discussed in § III. We choose the start time of the fit as the time when γ 1 takes the minimum. We denote the fitted fundamental mode as ψ iter [0 , 0] = C iter 0 ψ 0 , and we shall subtract it from blue curve to obtain the orange curve in Fig. 25, which is Ψ SXS -c -ψ iter [0 , 0] . Gray curve in Fig. 25 is the conservative estimation of the numerical errors (see Appendix C for details).  \nBefore discussing the iterative subtraction of the longest-lived mode, let us take a closer look of the fit of Ψ SXS -c by the fitting function ψ fit [0 , 3] . Figure 26 shows the mismatch M , amplitude A n , and rate of change γ n as a function of the start time of the fit t i . It should be compared with Fig. 17, but the panel for the relative error of C n is omitted here because there is no known  \nfiducial values for the numerical relativity waveform.  \nThe first panel of Fig. 26 shows the mismatch M between the SXS:BBH:0305 waveform Ψ SXS -c and the fitting function ψ fit [0 ,N ] , which is consistent with the results in [37]. Each curve in Fig. 26 represents what we used as the fitting function ψ fit [0 ,N ] . The larger the number of N , the smaller the mismatch M with earlier t i , which is the same trend as the mock waveform which we discussed in § III and § IV. Compared to Fig. 2 for the mock waveform without constant, the mismatch M in Fig. 26 has a lower bound in each curve.  \nThe second panel of Fig. 26 shows that the amplitude A n as a function of the start time of the fit t i , whose behavior is in agreement with [38]. We also present the phase ϕ n in the third panel of Fig. 26. Unlike the first panel, each curve represents the n -th overtone fitted by the fitting function ψ fit [0 , 3] . As expected, the fit of the fundamental mode is robust and not sensitive to the choice of the start time of the fit, which is also consistent with what we learned from the mock waveform analysis. On the other hand, for the overtones n ≥ 1, A n and ϕ n diverge, similar to Fig. 17 for the mock waveform with constant. Even though we subtracted the numerical constant observed in the late time of the SXS waveform, the fit is not as stable as the fit of the mock waveform analysis without a constant. We can see that the start time of the fit t i of the onset of divergence is earlier than the case of Fig. 17. This may be due to the spherical-spheroidal mode mixing in the SXS waveform, which plays a similar role to random noise, reducing the accuracy of the fit.  \nThe fourth panel of Fig. 26 shows the rate of change γ n . The each curve represents the n -th overtone same as the second panel. In parallel to A n and ϕ n , the behavior of γ n is also similar to Fig. 17 for the mock waveform with constant. γ 0 remains small even for late t i , which suggests that the plateau continues up to late start time of the fit and the fit is stable. On the other hand, for n ≥ 1, γ n approaches to a constant, which roughly coincides with | Im( ω n ) | , signaling the overfit.  \nAs we discussed in § III and § IV, we use the minimum of γ n as a criterion to extract the best-fit value of C n . In each step of the iterative procedure, to obtain the best-fit value of the longest-lived mode, we adopt the start time of the fit t i , when γ n for the next-longest-lived mode has the minimum. The reason why we focus on the minimum of γ n for the next-longest-lived mode rather than the of the longest-lived mode was that the latter can be isolated from the minima for other γ n , as in the case of the mock waveform without constant. Of course, it might be the case that γ n for the next-longest-lived mode would have the minimum isolated form other minima. Ultimately one needs to check the behavior of γ n as a function of t i to decide the optimal start time of the fit. In the case of the SXS:BBH:0305 waveform, we can see that the all γ n takes the minimum value at ( t i -t peak ) /M ≈ 11, which we regard as the optimal start time of the fit. We then adopt C 0 with the optimal start time of the fit as the best-fit value C iter 0 .  \nFIG. 27. Mismatch M (left top), amplitude A n (right top), phase ϕ n (left bottom), and rate of change γ n (right bottom) for the fit of SXS:BBH:0305 waveform Ψ SXS -c -ψ iter [0 , 0] after the subtraction of the fundamental mode by the fitting function ψ fit [1 ,N ] .  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nAfter we obtain the best-fit value C iter 0 , we subtract the fundamental mode from the blue curve to obtain the orange curve in Fig. 25, which is Ψ SXS -c -ψ iter [0 , 0] . Up to ( t -t peak ) /M ≈ 50, we see that the orange curve can be mainly described by a single damped sinusoid, which is the first overtone, but after that, we can see a peculiar behavior like the beat, which we discuss in detail in the next section.  \nIn Fig. 27, we present the result of the fit of Ψ SXS -c -ψ iter [0 , 0] waveform colored by orange in Fig 25 by the fitting function ψ fit [1 ,N ] . Similar to Fig. 26, we show mismatch M , the amplitude A n , and the rate of change γ n .  \nBy subtracting the fundamental mode, the mismatch shown in the first panel of Fig. 27 changed significantly from Fig. 26. In Fig. 26, the mismatch is bounded from below but the lower bound remain almost constant with respect to t i . However, in Fig. 27, the lower bound of the mismatch increases as t i increases. This trend is close to the case of the mock waveform with noise rather than the mock waveform without noise. It implies that even after subtracting the late time constant, the waveform still contains an effective noise such as the residual fundamental mode and spherical-spheroidal mixing.  \nThe amplitude A n and the rate of change γ n for each overtone in the fitting function ψ fit [1 , 4] are shown in the second and fourth panels of Fig. 27, respectively. By virtue of the subtraction, the stability of the first overtone n = 1 is slightly improved from Fig. 26. For instance, while the minimum value of γ 1 is about 10 -3 and is almost the  \nFIG. 28. The improvement of the amplitude A n (left column), phase ϕ n (right column) from n = 1 mode (first row) to n = 7 mode (seventh row) for the iterative fit of SXS:BBH:0305 by the fitting function ψ fit [ K,N ] with K = 0 , · · · , 7 and N = min( K +3 , 7).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \nγ  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nsame as Fig. 26, the region where γ n ≲ 10 -2 is longer due to the subtraction of the fundamental mode. The location of the plateau tends to move to earlier t i after subtraction. As expected, the improvement is not as manifest as in the mock analysis without constant.  \nFrom the fitting of Ψ SXS -c -ψ iter [0 , 0] , we extract the best-fit value C iter 1 . We then subtract the best-fit first overtone from the waveform, after which we obtain Ψ SXS -c -ψ iter [1 , 0] , which is shown by the green curve in Fig 25. By fitting this waveform, next we extract the best-fit value C iter 2 for the second overtone. We can continue the iterative fitting so long as the available data interval is longer than the time scale of the effective longestlived mode.  \nThus, starting from the original SXS:BBH:0305 waveform, we iteratively fit and subtract the longest-lived mode step by step and obtain the best-fit values C iter n , focusing on ( ℓ, m ) = (2 , 2) QNMs. This procedure accumulates the slight improvements of the stability of the fit of the overtones. The results are summarized in Figs. 28-31. Figure 28 shows the results of the fit by the fitting function ψ fit [ K,N ] with K = 0 , · · · , 7 and N = min( K + 3 , 7). The amplitude A n , phase ϕ n , and rate of change γ n are presented in each row for n = 1 , · · · , 7. For n = 1 , · · · , 4 modes, we use a superposition of four QNMs as a fitting function, since the superposition of the four modes is efficient to extract the overtones as we examined above. Here, as a reference, we also provide the result for n = 5 , 6 , 7 modes, which are fitted by fewer number of modes. The type of each curve indicates the fit at each step of the iterative fitting procedure. The solid curve represents the fit when the mode of interest is the longest-lived mode, from which we extract the best-fit value C iter n . The dot-dashed, dashed, and dotted curves represents the fit at the intermediate steps when the mode of interest is the second, third, and fourth longest-lived mode, respectively. Regarding the amplitude A n , we see that the plateau is extended and the fit becomes more stable as the mode of interest approaches the longest-lived mode. In particular, looking at the n = 2 mode, we find that the plateau of A 2 is almost nonexistent in the first fit where it is the third longest-lived mode, whereas the plateau extends roughly from t i ≃ 0 to 10 for the final fit where the n = 2 mode is the longest-lived mode. This is also the case for the phase ϕ n , which is also becoming more stable by iterative subtractions. The phase ϕ n sometimes fluctuates widely, but the iterative fitting method reduces the behavior. The γ n in the third column of Fig. 28 allows us to quantitatively confirm that the plateau is extended. As a demonstration, let us focus on the range of t i /M when γ n < 10 -1 is satisfied. For the n = 2 mode, in the first fit, such a range is only about 5. After iterating subtractions, the range is extended to 15 in the final fit, which is three times longer than without the subtraction procedure. For the other modes, we can also see that the extension of the plateau can be evaluated qualitatively by A n , ϕ n and quantitatively by γ n . The iterative fit-  \ntly improves the stable extraction of the overtones. This also applies to n = 5 , 6 , 7 modes which are fitted by a fewer number of modes.  \nIn Fig. 29, we highlight the difference between the conventional fit and the iterative fit. The dashed curves represent the conventional fit, which are obtained by a single fit of the waveform Ψ SXS -c with the start time of the fit t i = t peak and the fitting function ψ fit [0 , 7] including up to seventh overtone, same as [37]. The behavior of the dashed curves are consistent with Fig. 6 in [38]. While the fundamental mode and the first overtone are fitted in a stable manner, the fit of the overtones n ≥ 2 are sensitive to the choice of the start time of the fit, suggesting that the fit is not stable. The stability of the overtones is improved in the iterative fit, as shown by the solid curves in Fig. 29. From Fig. 29 and the γ n in Fig. 28, we note that the plateau begins from earlier t i for higher overtones, which makes sense given that the higher overtones are responsible to early time ringdown.  \nIt might appear from Fig. 29 that the higher overtones up to n = 7 may be extracted in a stable manner. However, we should interpret the results carefully by checking the waveform at each iterative subtraction. Figure 30 represents SXS:BBH:0305 waveform (black solid), the extracted longest-lived mode at each iterative step (colored solid) by fitting the waveform obtained by the previous step (dashed), and the waveform after the subtraction of the extracted mode (dotted). For lower overtones, we can observe several oscillations before the waveform is contaminated by numerical errors. On the other hand, for higher overtones, since they decay more quickly, the waveform is contaminated by numerical errors before exhibiting several oscillations. Therefore, the range of data of waveform over the numerical errors which we can use for the fitting becomes shorter. Here, as a rough criterion, let us adopt the number of oscillations between the peak time and the time when numerical errors start to become dominant. We see that for n ≥ 5, the waveform does not demonstrate a single period of damped sinusoid between the time range. Of course, the peak time itself does not necessarily mean the starting time of the ringdown, so a more physical criterion needs to be established. Nevertheless, we should bear in mind that even if the plateau is extended for n ≥ 5 in Fig. 29, it is necessary to carefully determine up to which overtones can be extracted without overfit. Below, we focus on the modes up to n = 4 for the iterative fitting method.  \nIn Fig. 31, we compare the results in the alternative way, where we use the fitting function consisting of the modes up to the fourth overtone for both conventional and iterative fit. Again, the plateau is extended by the iterative fit compared to the conventional fit. The improvement is also consistent with what we expected from the improvement in Figs. 20-22 for the mock waveform analysis. As we can see from Figs. 29 and 31, the difference between the extracted values in the two ways becomes larger for higher overtones, which is also in parallel to Figs. 20-22 for the case of mock waveform with noise.  \nFIG. 29. Comparison of the amplitude A n (top row) and phase ϕ n (bottom row) obtained by fitting SXS:BBH:0305 by the conventional method (dashed) with the fitting function ψ fit [0 , 7] and the iterative method (solid) with the fitting function ψ fit [ K,N ] with K = 0 , · · · , 7 and N = min( K +3 , 7).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nBased on the comparison of Figs. 20-22 and the fiducial value, we expect that the plateau in Fig. 29 gives a closer estimation of the underlying coefficient of each mode contained in the waveform.  \nTherefore, for the following, we mainly focus on the former way of comparison, i.e., Fig. 29 with the modes up to n = 4. While there exists the above subtleties, we confirm that the lessons from the mock waveform analyses indeed work for the numerical relativity waveform and the subtractions helps us to extract the overtones in a more stable manner.  \nThe best-fit values obtained by the conventional fit and iterative fit of the waveform Ψ SXS -c in Fig. 29 are listed in Table III. For the conventional fit, we fix t i = t peak and ( t e -t peak ) /M = 150 and obtain the amplitude and phase by fitting Ψ SXS -c by the fitting function ψ fit [0 , 7] . While there is a slight difference, e.g., the fitted waveform is Ψ SXS in [37] but Ψ SXS -c for our analysis, other setup of the fitting is the same and the resultant values of A n and ϕ n are consistent with [37]. For t i and t e we show  \nonly two digits as we adopt the closest value in a given grid of SXS waveform.  \nOn the other hand, for the iterative fit, we iteratively adopt the fitting function ψ fit [ K,N ] with K = 0 , · · · , 7 and N = min( K +3 , 7) as in Fig. 29. We optimize the range [ t i , t e ] of the data for the fit as follows: When extracting n -th overtone, we determine t i as the fit start time when γ m with m = min( n +1 , 7) is minimum as we described above and t e as the time around which the corresponding waveform in Fig. 30 is comparable to the numerical errors, though we confirm that the result is not so sensitive to t e even if we extend it into the regime where the numerical errors are dominant. The choice of t i is more important, and we expect that t i is earlier for higher overtones as higher overtones are dominant at earlier times. Indeed, from Table III we confirm the trend that t i for the higher overtones are getting earlier. This trend stops when t i < t peak for n ≥ 4.  \nAs we explained above, for n ≥ 5, the iterative method is less efficient since the subtracted waveform is contam-  \n-  \nTABLE III. Best-fit values for the amplitude A n and phase ϕ n by fitting the SXS:BBH:0305 waveform subtracting the numerical constant with the conventional and iterative methods. For the conventional fit, we set the fit start time ( t i -t peak ) /M = 0 and the fit end time ( t e -t peak ) /M = 150.  \ninated by numerical errors before exhibiting a single period of oscillation after the peak time. Therefore, we consider that the extracted values A n and ϕ n for n ≥ 5 may be less meaningful, and present them in parentheses in Table III as a reference value. Still, a criterion up to which overtones can be extracted without overfit is nontrivial. In addition to the number of oscillations, another simple criterion would be whether t i for the best fit is earlier than some threshold value. If we adopt the peak time as a threshold, n = 4 would not be extracted. However, while t i for the best fit decreases and becomes even before the peak time, note that the plateau in Fig. 29 extends after the peak time and hence the values evaluated with t i after the peak time are the same order. For instance, ( A 4 , ϕ 4 /π ) = (21 . 46 , -0 . 8107) for t i = t peak whereas ( A 4 , ϕ 4 /π ) = (21 . 99 , -0 . 8115) for ( t i -t peak ) /M = -2. Note that this criterion is relevant to the debate on the start time of ringdown, which varies from before to ∼ 10 M after the peak time, making it difficult to draw conclusions at this stage.  \nLet us focus on the amplitudes and phases up to n = 4. From Table III, we can see that the difference between the two fitting methods becomes larger for higher overtones. Specifically, for the fundamental mode, the relative difference between conventional and iterative method in C n is about 10 -3 , but it reaches about 20% for n = 3, and 40% for n = 4. From the analysis of the mock waveform shown in Fig. 20, where the fiducial values are stably and correctly extracted by the iterative fit, we expect that the best-fit values obtained by the iterative fit capture the contribution of each mode more appropriately.', 'B. Spherical-spheroidal mixing in SXS:BBH:0305': 'As we saw above, in Fig. 27, the subtraction of the lower modes reveals the existence of the beat. To determine whether the beat originates from the numerical errors, we evaluate the numerical errors by comparing the numerical relativity waveform obtained with different resolution and extrapolation order, which are shown  \nas green and magenta thin dotted curves, respectively (see Appendix C). The envelope of the all error curves is shown as a gray curve, which we adopt as the conservative estimation of the numerical errors. This is because for the ringdown analysis, higher order extrapolation tends to lead an overfit, and hence a mild extrapolation order is recommended [44]. We note that the rough amplitude of the first beat in the waveform Ψ SXS -c -ψ iter [0 , 0] around ( t -t peak ) /M ≈ 60 is larger than the errors estimated by the next-highest resolution. Further, the beat shows up more clearly above the estimated numerical errors in the waveform Ψ SXS -c -ψ iter [1 , 0] after the subtraction of the first overtone, as shown by the green curve in Fig. 25. Therefore, it would be reasonable to attribute the beat to the effect originating from the spherical-spheroidal mixing rather than the numerical errors. Indeed, in the case of the mock waveform, after the subtraction of the fundamental mode we see the first overtone and such behavior does not show up. The appearance of the beat is peculiar to the numerical relativity waveform.  \nActually, the QNM frequency of ( ℓ, m, n ) = (3 , 2 , 0) mode is close to the (2 , 2 , 0) fundamental mode (see Table VIII). Considering the possibility that the fundamental mode is not completely subtracted, it is expected that the beat between the residual of (2 , 2 , 0) mode and (3 , 2 , 0) mode via the spherical-spheroidal mixing occurs. Indeed, the beat frequency is given by M ( ω 220 -ω 320 ) = -0 . 2268 + 0 . 003404 i , which is consistent with the time scale of the beat appearing in Fig. 25.  \nLet us discuss the impact of the spherical-spheroidal mixing in the fit of the ( ℓ, m ) = (2 , 2) mode of the waveform SXS:BBH:0305. By taking into account the dominant mixing of ( ℓ, m, n ) = (3 , 2 , 0) mode, we consider the following fitting function  \nψ fit [ K,N ]+320 = N ∑ n = K C 22 n ψ 22 n + C 320 ψ 320 . (14)  \nIn this subsection, we recover the subscript ( ℓ, m ). In the iterative fitting procedure, we subtract the best-fit  \nFIG. 30. SXS:BBH:0305 waveform (black solid), the extracted longest-lived mode at each iterative step (colored solid) by fitting the waveform obtained by the previous step (dashed), and the waveform after the subtraction of the extracted mode (dotted).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 32. Numerical relativity waveform SXS:BBH:0305. Waveforms after the subtraction of the numerical constant (blue), the best-fit values of ( ℓ, m, n ) = (2 , 2 , 0) mode (light brown), (3 , 2 , 0) mode (orange), and then (2 , 2 , 1) mode (green).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 31. Comparison of the amplitude A n (top) and phase ϕ n (bottom) obtained by fitting SXS:BBH:0305 by the conventional method (dashed) with the fitting function ψ fit [0 , 4] and the iterative method (solid) with the fitting function ψ fit [ K, 4] with K = 0 , · · · , 4.  \nlongest-lived mode in each step. Therefore, following the ordering | Im ω 220 | < | Im ω 320 | < | Im ω 221 | < · · · , we subtract (2 , 2 , 0), (3 , 2 , 0), (2 , 2 , 1), and so on. For the each step of the iteration, we denote the sum of the best-fit effective longest-lived modes as ψ iter [0 , 0] , ψ iter [0 , 0]+320 , ψ iter [0 , 1]+320 and so on.  \nWe present the waveform of the ( ℓ, m ) = (2 , 2) mode of the SXS:BBH:0305 simulation after subtracting the numerical constant (blue) in Fig. 32, together with the waveform after the subtraction of ( ℓ, m, n ) = (2 , 2 , 0) mode (light brown), (3 , 2 , 0) mode (orange), and then (2 , 2 , 1) mode (green) as well. After the subtraction of the best-fit value of the (3 , 2 , 0) mode, we can observe oscillations originating from (2 , 2 , 1) mode a little longer up to ( t i -t peak ) /M ∼ 60, and the beat disappears. Further, after the subtraction of the best-fit value of the (2 , 2 , 1) mode, the beat in the green curve in Fig. 25 disappears in Fig 32. Therefore, it is reasonable to attribute  \n<!-- image -->  \n-  \nthe origin of the beat to the spherical-spheroidal mixing.  \nWe present the comparison of the results obtained by the conventional fit and iterative fit when taking into account the (3 , 2 , 0) mode in Fig. 33, and the best-fit values in Table IV. Naively, taking the spherical-spheroidal mixing into account, the stability of the fit of overtones may be improved. While the contribution of the (3 , 2 , 0) mode is as small as A 320 = O (10 -3 ), it affects the fit of the higher overtones of the (2 , 2 , n ) modes. Compared to Table III without the (3 , 2 , 0) mode, while there is relatively small difference for n ≤ 2, the (2 , 2 , n ) modes with n ≥ 3 differ significantly. As another difference, t i for the best fit is shifted to later time. Note that the validity of n = 3 , 4 are marginal in the sense that the fitting interval is shorter than the lower modes and hence the fit of higher overtones is still challenging as we discussed above.  \nOur analysis shows that the iterative fitting method is efficient to isolate the spherical-spheroidal mixing. The stability of the fit of the lower modes is improved when the spherical-spheroidal mixing compared to the conventional fit. However, after the inclusion of the mixing, the fit of the higher overtones becomes more sensitive and the best-fit values deviate significantly. Unfortunately, the analysis of the mixing is limited by the numerical errors. While our analysis suggests the efficiency of the iterative fitting method to study the spherical-spheroidal mixing, a more detailed analysis with reduced errors would be necessary to obtain a rigorous understanding.  \n-  \n-  \n-  \n-  \nTABLE IV. Best-fit values for the amplitude A n and phase ϕ n by fitting the ( ℓ, m ) = (2 , 2) mode of SXS:BBH:0305 waveform subtracting the numerical constant with the conventional and iterative methods, taking into account the spherical-spheroidal mixing with ( ℓ, m, n ) = (3 , 2 , 0) mode. For the conventional fit, we set ( t i -t peak ) /M = 0 and ( t e -t peak ) /M = 150.', 'C. SXS:BBH:0158 - high spin simulation': 'Next, we explore the fit of the waveform SXS:BBH:0158, whose remnant dimensionless spin χ rem = 0 . 9450, as a representative example of the high spin simulation. As we mentioned above, as we increase the black hole spin, the QNMs damp more slowly and oscillate more rapidly, and most QNM frequencies become eventually degenerate. By analyzing the high spin simulation, we shall see how these factors affect the extraction of overtones.  \nWe present the waveform of SXS:BBH:0158 in Fig. 34. It is clear that the raw data waveform shown by the black curve is contaminated by a numerical constant at late time, whose value is about c = (7 . 4 -4 . 8 i ) × 10 -5 . While the waveform data of SXS:BBH:0158 are available up to ( t -t peak ) /M = 400, we did not show the full waveform since the late time waveform simply remains the constant and does not have physical meaning. For the following we use the waveform Ψ SXS -c shown by the blue curve.  \nAs expected, we can observe more oscillations compared to Fig. 25, which would make the fitting easier than the case of SXS:BBH:0305. We present the result of the fitting in Fig. 45 in Appendix B, which should be contrasted with Figs. 26 and 27. As we did in the SXS:BBH:0305 analysis, we fit the waveform Ψ SXS -c and obtain the best-fit value C iter 0 by using the fit starting from the time where γ 1 takes the minimum. We then subtract the fundamental mode from the blue curve to obtain the orange curve in Fig. 34, the latter of which is Ψ SXS -c -ψ iter [0 , 0] . Up to ( t -t peak ) /M ≈ 70, we see that the orange curve can be mainly described by a single damped sinusoid, which is the first overtone. The damping time seems changing for ( t -t peak ) /M ≳ 70 but we expect that the waveform is subject to the numerical errors there.  \nWe iterate the fit and subtraction and collect the best-fit value for the longest-lived mode at each step as C iter n . The results are shown in Fig. 35. We can see  \nthat the plateau of A n and ϕ n is longer than the case of SXS:BBH:0305. Even before the subtraction of the longest-lived mode, we can observe a mild plateau up to n = 2, while there is no plateau on the second overtone in SXS:BBH:0305 before the subtractions. The stability of the fit vary greatly depending on the ringdown waveform.  \nWe can evaluate the stability or the length of the plateau quantitatively by the plot of γ n in Fig. 35. The region where γ n remains small is longer than the case of SXS:BBH:0305. Specifically, for the fit of Ψ SXS -c , the range of t i where γ 1 < 10 -1 is about -5 < ( t i -t peak ) /M < 20 for SXS:BBH:0305, but it is about 10 < ( t i -t peak ) /M < 45 for SXS:BBH:0158, i.e., it is extended about ∆ t/M ≃ 10. Furthermore, the region of t i where γ 2 < 10 -1 is about 5 < ( t i -t peak ) /M < 10 for SXS:BBH:0305, while it is about 10 < ( t i -t peak ) /M < 30 for SXS:BBH:0158, which is about four times larger than the range for SXS:BBH:0305. Therefore, as expected, we find that the fit tends to be more stable for higher spin black holes, as the QNMs have slower damping time and higher real frequency.  \nThe iterative fitting method further improves the stability of the fit, and the plateau of A n and ϕ n is extended. From the middle column in Fig. 35, we see that the stability of ϕ n improves significantly in all panels. Especially, the improvement of ϕ 3 (red) and ϕ 4 (purple) are manifest. While A 3 and A 4 are also flattened, they are still tilted. We expect that additional subtractions would make them sufficiently flat, for which we need the waveform with smaller numerical errors.  \nWe represent a comparison of the A n and ϕ n obtained by the conventional fit and the iterative fit in Fig. 36. It shows that the iterative fit provides more stable fit than the conventional fit for all overtones up to the fourth overtone. On the other hand, we can see that some values of A n obtained by each fit are significantly different between the conventional fit and our method, especially for the higher overtone. It can also be confirmed from Ta-  \nFIG. 34. Numerical relativity waveform SXS:BBH:0158. Raw data (black dashed) and waveforms after the subtraction of the numerical constant (blue) and the best-fit fundamental mode (orange).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 33. Comparison of the conventional fit (dashed) and iterative fit (solid) of SXS:BBH:0305 waveform including ( ℓ, m, n ) = (3 , 2 , 0) mode.  \nble V, where we compare the best-fit values. *1 It should be contrasted to the case of SXS:BBH:0305, where the conventional fit and the iterative fit provide the best-fit values of the same order. The two fitting methods do not necessarily provide a similar best-fit values. The consistency is not a priori guaranteed, but depends on each specific waveform. From the point of the stability of the fit, it would be plausible to expect that the best-fit values obtained from the iterative fit provides a more reasonable extraction. Note also that the best-fit values of overtones  \n<!-- image -->  \n-  \nup to n = 4 are extracted at relatively large t i . Therefore, if the numerical errors are smaller, one may be able to extract overtones n ≥ 5.', 'D. SXS:BBH:0156 - low spin simulation': 'Let us investigate the fit of the waveform of the SXS:BBH:0156 with χ rem = 0 . 3757 shown in Fig. 37 as an example of the low spin simulation. Again, while the waveform data of SXS:BBH:0156 are available up to ( t -t peak ) /M = 400, we did not show the full waveform since at late time it remains roughly constant c = (1 -0 . 3 i ) × 10 -4 and fluctuate around it. As expected, compared to SXS:BBH:0305 and SXS:BBH:0156, we can observe a fewer number of oscillations before the contamination by the constant, which originates from the short damping time of the QNM frequencies of low spin Kerr black holes. Since we do not have sufficient number of sampling points, we use the cubic spline interpolation to increase the number of sampling points by two times for the regime ( t -t peak ) /M ≥ -60. For this regime, the raw data spacing is almost constant ∆ t/M ≈ 0 . 5, so after the interpolation, we have ∆ t/M ≈ 0 . 25.  \nSince the interval of available data is shorter, we expect that the fit of the low spin simulation is more challenging. Nevertheless, from Figs. 38 and 39, we see that relatively clear plateaus show up. Actually, while the range of the plateau is short, A n remains almost constant, in contrast to the fact that for the high spin SXS:BBH:0158 waveform, A n shows a continuous mild growth as a function of t i .  \nAfter the subtraction of the third overtone, the waveform exhibits about half period of oscillation between the peak time and the time when the waveform becomes  \nFIG. 35. The improvement of the amplitude A n (left column), phase ϕ n (right column) from n = 1 mode (first row) to n = 4 mode (fourth row) for the iterative fit of SXS:BBH:0158 by the fitting function ψ fit [ K,N ] with N = K +3.  \n<!-- image -->  \ncomparable to numerical errors. A half period would be more or less a threshold of the necessary range of data to extract damped sinusoid. Further, the time of the strain peak may not exactly coincide with the time of the onset of the QNM oscillations. They may be triggered even before the peak time. Taking into account the waveform before the peak time, we can include more than half period of oscillation to the fitting analysis. Indeed, we see that the plateau in Fig. 39 is extended before the strain peak time.  \nIn Table VI, We compare the best-fit values obtained by the conventional fit and iterative fit. We see that in this case the two best-fit values are the same order. Note that extracted amplitudes and phases slightly change when the number of sampling points changes.  \nTABLE V. Best-fit values for the amplitude A n and phase ϕ n by fitting the SXS:BBH:0158 waveform subtracting the numerical constant with the conventional and iterative methods. For the conventional fit, we set ( t i -t peak ) /M = 10 and ( t e -t peak ) /M = 220.  \n×  \n-  \n-  \n-  \n-  \nTABLE VI. Best-fit values for the amplitude A n and phase ϕ n by fitting the SXS:BBH:0156 waveform subtracting the numerical constant with the conventional and iterative methods. For the conventional fit, we set ( t i -t peak ) /M = 0 and ( t e -t peak ) /M = 150.  \n-  \n-  \n-  \n-  \n-', 'E. SXS:BBH:1108 - GW190814-like simulation': 'Finally, we consider the SXS:BBH:1108 waveform, which is related to the GW190814. GW190814 [47] is gravitational wave signal from a coalescence of 23 M ⊙ black hole with a 2 . 6 M ⊙ compact object, the latter of which is either a light black hole or a heavy neutron star. The source has the most unequal mass ratio among the gravitational waves observed so far, and the dimensionless spin of the primary black hole is constrained to ≤ 0 . 07. Assuming the secondary component is a black hole, SXS:BBH:1108 serves as a waveform consistent with GW190814.  \nCompared to the three numerical relativity waveforms considered so far, the SXS:BBH:1108 waveform shown in Fig. 40 is contaminated by large numerical errors. The waveform roughly approaches a constant c = ( -3 . 1 -8 . 8 i ) × 10 -6 , but it is below the estimated numerical errors. Also, the waveform damps quickly since the SXS:BBH:1108 is the low spin simulation with χ rem = 0 . 2772. However, in the case of SXS:BBH:1108 the sampling points are not insufficient, so we do not use the spline interpolation.  \nThe results of the fitting analysis are shown in Figs. 41 and 42 and Table VII. Their qualitative behavior is the  \nsame as the results for the SXS:BBH:0156. Our results imply that the extraction of the QNM depends mainly on the spin of the remnant black hole, rather than the mass ratio of the binary before merger.', 'VI. CONCLUSIONS AND DISCUSSION': 'The extraction of QNMs from black hole ringdown is a key to understand gravity at dynamical and strong field regime. For the extraction of the first overtone from the observed ringdown gravitational waves, the conflicting results have been reported mainly due to the low signal-tonoise ratio of the observational data. Even for the simulated waveform by numerical relativity calculations, the extraction of overtones is not straightforward, since the higher overtones more quickly damp and become subdominant compared to lower modes, and hence it tends to end up with the overfit. In the present paper, we proposed an improved way of extracting the overtone QNMs. In contrast to the conventional fitting method, where one extracts the coefficients of all the QNMs of interest from a single fitting, we considered an iterative procedure of fitting and subtraction to extract the longest-lived QNM at each step. The iterative fitting method is summarized  \n-  \n-  \n-  \n-  \nTABLE VII. Best-fit values for the amplitude A n and phase ϕ n by fitting the SXS:BBH:1108 waveform subtracting the numerical constant with the conventional and iterative methods. For the conventional fit, we set ( t i -t peak ) /M = 0 and ( t e -t peak ) /M = 150.  \nFIG. 37. Numerical relativity waveform SXS:BBH:0156. Raw data (black dashed) and waveforms after the subtraction of the numerical constant (blue) and the best-fit fundamental mode (orange).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 36. Comparison of the conventional fit (dashed) and iterative fit (solid) of SXS:BBH:0158 waveform.', 'as follows:': "- 1. Fit a given waveform by a fitting function composed of a superposition of the four longest-lived damped sinusoids with coefficients C n as fitting parameters.  \n<!-- image -->  \n-  \nSet the end time of the fit t e around the time when the signal becomes comparable to noises.  \n- 2. Check if there exists a certain range of the start time of the fit t i when the rate of change γ n = | 1 C n dC n dt i | remains small compared to | Im( ω n ) | to guarantee the stable fit rather than the overfit.\n- 3. Read off the best-fit value of the coefficient C n for the longest-lived mode by a fit starting from a time when the decrease of the rate of change γ n for the next-longest-lived mode saturates (not necessarily exactly the minimum since it tends to oscillate).\n- 4. Subtract the best-fit longest-lived damped sinusoid from the waveform.\n- 5. Iterate the above steps so long as the available data interval is longer than the time scale of the effective longest-lived mode.  \nFIG. 38. The improvement of the amplitude A n (left column), phase ϕ n (right column) from n = 1 mode (first row) to n = 4 mode (fourth row) for the iterative fit of SXS:BBH:0156 by the fitting function ψ fit [ K,N ] with N = K +3.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \nγ  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nThis method enables us to reveal the subdominant contribution of the higher overtones by iteratively subtracting the dominant contribution of the lower modes. Through the analyses of the mock waveform and numerical relativity waveform, we clarified that by collecting the best-fit values of the longest-lived mode at each step we can improve the stability of the fit of overtones. We confirmed that the coefficients for the higher overtones are sensitive to the start time of the fit, and they begin to diverge exponentially for the fit starting from late time.  \nThe divergence of the coefficient simply compensates the exponential damping, which can be clearly seen from the fact that the rate of change γ n roughly coincides with the damping rate | Im( ω n ) | , signaling the overfit. Therefore, we should adopt the start time of the fit where γ n remains smaller than | Im( ω n ) | , for which it is important to include higher overtones for the fitting function to improve the stability of the fitting of the lower overtones. Indeed, we found that as we increase the number of overtones, the plateau of the fit of the lower overtone coeffi-  \nFIG. 40. Numerical relativity waveform SXS:BBH:1108. Raw data (black dashed) and waveforms after the subtraction of the numerical constant (blue) and the best-fit fundamental mode (orange).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 39. Comparison of the conventional fit (dashed) and iterative fit (solid) of SXS:BBH:0156 waveform.  \ncients are extended. This is because the higher overtones are responsible for describing the early time waveform, and hence with the higher overtones included, the lower overtones are appropriately fitted at early time. However, inclusion of too many modes in the fitting function tends to destabilize the fitting of the lower overtones. This is why we adopted a fitting function composed of four longest-lived damped sinusoids.  \nA caveat for the above step 1 is that it is still unclear whether it is always optimal to adopt the four longestlived modes. We adopted the four modes based on the mock waveform related to SXS:BBH:0305, whose dimensionless remnant spin ∼ 0 . 7 is more or less typical value for the remnant black hole. In this sense, the four modes should be regarded as a representative number, and it would be safe to identify the optimal number for each waveform by checking the stability of the extracted coefficients.  \nAlso, in the step 3, one may think that it would be more natural to use γ n for the longest-lived mode as an indicator of the stability of the fit. The reason why we  \n<!-- image -->  \n-  \nfocus on the γ n for the next-longest-lived mode is based on our analysis of the mock waveform of the pure damped sinusoids, where the minimum of γ n for the longest-lived mode differs from the minimum of γ n for the other modes. However, for waveforms with noises, the mimima of γ n tend to match each other. In such cases, there is no reason to focus on γ n for the next-longest-lived mode, and we are simply adopting the start time of the fit t i where γ n are minimized.  \nWe examined how the fit of the SXS numerical relativity waveform can be improved by the iterative extractions. Iterating the subtractions helps to reveal the subdominant contribution including the spherical-spheroidal mixing and improve the stability of the fit of the overtones. We clarified that for the GW150914-like numerical relativity waveform SXS:BBH:0305, the overtones with n ≥ 2 can be extracted in a more stable manner, compared to the conventional fit where only n = 0 , 1 can be stably extracted. In addition to SXS:BBH:0305, we also investigated high spin simulation SXS:BBH:0158 and low spin simulations SXS:BBH:0156 and SXS:BBH:1108. As expected, for the high spin case, we have more oscillating signals. Even before the subtraction of the longest-lived mode, we can observe a mild plateau for the fit of overtones up to n = 2, while there is no plateau on the second overtone in SXS:BBH:0305 before the subtraction. We found that the best-fit values for SXS:BBH:0158 obtained by the iterative fit differ significantly from the conventional fit, especially for the higher overtones. On the other hand, for the low spin case, we confirmed that the subtraction still improves the fit of overtones, while the range of the plateau is shorter than the high spin case.  \nWe studied how the fit is affected by noises, which eventually hides the ringdown signal at late time. While  \n1  \nA  \n3  \nA  \n4  \nA  \n10  \n10  \n-  \n10  \n10  \n10  \n10  \n-  \n10  \n10  \n10  \n10  \n-  \n10  \n-  \n0  \n1  \n-  \n10  \n0  \n10  \n20  \npeak  \n[  \nK,N  \n] = [0  \n,  \n3]  \n[  \nK,N  \n] = [1  \n,  \n4]  \n(  \nt  \ni  \n-  \nt  \n)  \n/ M  \nFIG. 41. The improvement of the amplitude A n (left column), phase ϕ n (right column) from n = 1 mode (first row) to n = 4 mode (fourth row) for the iterative fit of SXS:BBH:1108 by the fitting function ψ fit [ K,N ] with N = K +3.  \n<!-- image -->  \n4  \n2  \n0  \n2  \n5  \n3  \n1  \n1  \n3  \n[  \nK,N  \n] = [0  \n,  \n3]  \n[  \nK,N  \n] = [1  \n,  \n4]  \n[  \nK,N  \n] = [2  \n,  \n5]  \n[  \nK,N  \n] = [3  \n,  \n6]  \n[  \nK,N  \n] = [1  \n,  \n4]  \n[  \nK,N  \n] = [2  \n,  \n5]  \n[  \nK,N  \n] = [3  \n,  \n6]  \n[  \nK,N  \n] = [4  \n,  \n7]  \n-  \n5  \n0  \n5  \n10  \n15  \n20  \n(  \nt  \ni  \ni  \n-  \n-  \n-  \nt  \nt  \npeak  \n)  \n/ M  \n)  \n/ M  \n(  \nti  \n(  \nti  \n(  \nti  \n(  \nti  \n-  \n-  \n-  \n-  \nt  \nt  \nt  \nt  \n10  \n0  \n10  \n20  \n-  \n-  \n-  \n-  \n10  \n0  \n10  \n20  \n10  \n0  \n10  \n20  \n10  \n0  \n10  \n20  \npeak  \npeak  \npeak  \npeak  \n)  \n/ M  \n)  \n/ M  \n)  \n/ M  \n)  \n/ M  \n5  \n0  \n5  \n10  \n15  \n20  \n5  \n0  \n5  \n10  \n15  \n20  \n5  \n0  \n5  \n10  \n15  \n20  \n(  \nt  \ni  \ni  \ni  \ni  \n-  \n-  \n-  \n-  \nt  \nt  \nt  \nt  \npeak  \n10  \n0  \n10  \n20  \n-  \n-  \n-  \n-  \n10  \n10  \n10  \npeak  \n)  \n/ M  \n)  \n/ M  \n)  \n/ M  \n)  \n/ M  \nthe stability of the fit is reduced by noises at late time, the iterative fit is efficient to extract overtones at early time. Up to which overtone a sufficiently stably extracted depends on the magnitude of the noise and the damping time of the overtone. A rough criterion would be given by the number of oscillations between the peak time and the time when the noise becomes dominant. For overtones whose number of oscillations is more than one, we expect its stable extraction.  \nWe also confirmed that the iterative fitting method is  \nrobust to the data with low sampling rates, compared to the conventional fit. These natures are advantageous for the extraction of the QNMs from observed data. In particular, a recent analysis of the observed ringdown gravitational waves shows that the conventional analysis is sensitive to the sampling rate [31]. Iterative fitting method would be more advantageous to apply to the observed data as it provides robust fitting results even for data with low sampling rates, compared to conventional fitting method.  \n(  \nt  \npeak  \n(  \nt  \npeak  \n-  \n5  \n0  \n5  \n10  \n15  \n20  \n(  \nt  \npeak  \n-  \n-  \n10  \n10  \nπ  \n/  \n1  \nφ  \nπ  \n/  \n2  \nφ  \nπ  \n/  \n3  \nφ  \nπ  \n/  \n4  \nφ  \n-  \n-  \n-  \n-  \n-  \n-  \n-  \n1  \n.  \n5  \n1  \n.  \n0  \n0  \n.  \n5  \n0  \n.  \n0  \n0  \n.  \n5  \n0  \n.  \n5  \n0  \n.  \n0  \n0  \n.  \n5  \n1  \n.  \n0  \n1  \n.  \n5  \n1  \n.  \n5  \n1  \n.  \n0  \n0  \n.  \n5  \n0  \n.  \n0  \n0  \n.  \n5  \n1  \n.  \n0  \n0  \n.  \n5  \n0  \n.  \n0  \n0  \n.  \n5  \n1  \n.  \n0  \n1  \nγ  \n2  \nγ  \n3  \nγ  \n4  \nγ  \n10  \n10  \n10  \n-  \n10  \n-  \n10  \n10  \n10  \n-  \n10  \n-  \n10  \n10  \n-  \n10  \n-  \n10  \n10  \n10  \n-  \n10  \n-  \n1  \n0  \n1  \n2  \n1  \n0  \n1  \n2  \n0  \n1  \n2  \n1  \n0  \n1  \n2  \n-  \n-  \n-  \n(  \nt  \nFIG. 42. Comparison of the conventional fit (dashed) and iterative fit (solid) of SXS:BBH:1108 waveform.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nWe expect that the iterative fitting still works to the case of the QNM spectral instability under small perturbation of the effective potential. Such a perturbation could drastically alter the QNM frequencies, but it modifies the ringdown signal only at the late time, so it is not fatal for the black hole spectroscopy program [18, 19]. Since the iterative fitting method allows us a more direct way to extract the QNMs from the ringdown signal in the time domain, it is plausible that the method is robust against the spectral instability.  \nThe fitting of the ringdown by QNMs is a nontrivial task and various approaches would be possible. An approach closely related to ours is to consider a filter in the frequency domain, which removes the QNM poles of interest on the complex frequency plane [28, 29, 41]. The filter also affects to the contribution from other modes, and hence the amplitudes and phases obtained by the fit of the filtered waveform are different from the original ones. While one can reconvert them and estimate the original values, it might be the case that the modification associated with the filtering affects the fit of the  \nsubdominant modes. For instance, if the amplitude of a subdominant mode becomes significantly small due to the filtering, it may be the case that the fit of the subdominant mode would become difficult or be destabilized. Since the amplitudes and phases are related to the excitation coefficient of QNMs [59-64], it would be ideal if we could extract them without additional modifications by the mode cleaning process. The iterative fitting method in the time domain which we consider in the present paper does not affect other modes, and allows us to extract subdominant modes in a stable manner.  \nAnother approach of a similar spirit would be the 'bootstrap procedure' suggested in [38], where one first identifies longer-lived modes by a theory-agnostic fitting, and then fixes their frequencies when searching for shorter-lived modes. Very recently, this fitting strategy is further explored in [65], where after checking the stability of the amplitudes and phases of QNMs, the unstable modes are iteratively removed from the fitting function so that the all the modes are stably fitted. Our strategy also adopts the iterative procedure but is complementary in the sense that we iteratively extract and subtract the most stable longest-lived mode one-by-one from the waveform and then consider the fitting of subdominant modes.  \nWhile we mainly focus on the extraction of overtones in the present paper, the extraction of subdominant modes with the iterative fitting method deserves further study as we clarified that the higher overtone coefficients are changed depending on whether we include the sphericalspheroidal mixing in the fitting function. Nonlinear effects predicted by second-order perturbation theory [6674] would also be an interesting target [41, 58, 75-80]. Further, while we focused on the frequency-fixed fitting of numerical relativity waveform in the present paper, we expect that our fitting method may also be efficient for a theory-agnostic fitting and/or observed ringdown signal. For instance, we expect to improve the theoryagnostic fitting by performing, e.g., the Prony method first to identify the QNM frequency and coefficient of the longest-lived mode, subtracting it from the ringdown data, and iterating the procedure to extract higher overtones. It would also be interesting to apply our fitting method to analyze the overtones in the ringdown gravitational waves emitted from a particle plunging into a black hole [64, 81-84]. We leave these topics for future work.", 'ACKNOWLEDGMENTS': 'H.M. was supported by Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research (KAKENHI) Grant No. JP22K03639.  \nTABLE VIII. Kerr QNM frequencies corresponding to SXS:BBH:0305 for ( ℓ, m ) = (2 , 2) and (3 , 2).TABLE IX. Kerr QNM frequencies corresponding to SXS:BBH:0158 for ( ℓ, m ) = (2 , 2).', 'Appendix A: Kerr QNM frequencies': 'In this appendix, we provide the QNM frequencies corresponding to the SXS simulations used in the main text. For SXS:BBH:0305, we provide the QNM frequencies for ( ℓ, m ) = (2 , 2) and (3 , 2) modes in Table VIII. For other simulations, the QNM frequencies of ( ℓ, m ) = (2 , 2) mode are listed in Tables IX-XI, respectively. All QNM frequencies are generated by using the qnm package [85] with machine precision, and here we present values rounded to four digits as a reference.  \nNote that, since we use the normalization by the total binary mass M , to obtain the QNM frequencies describing the ringdown signal after merger, we multiply a factor M/M rem to the Kerr QNM frequencies associated with the dimensionless remnant spin, where M rem is the remnant mass after merger. As a specific example, let us consider SXS:BBH:0305, whose dimensionless remnant spin is a/M rem = 0 . 6921. The fundamental mode of the ( ℓ, m ) = (2 , 2) mode of the Kerr QNM frequencies for a/M rem = 0 . 6921 is given by M rem ω = 0 . 5291 -0 . 08109 i . Using the normalization by the total binary mass M before merger, the QNM frequency is given by Mω = (0 . 5291 -0 . 08109 i ) M/M rem . Since M rem /M = 0 . 9520 for SXS:BBH:0305, we obtain Mω = 0 . 5558 -0 . 08517 i . The QNM frequencies in Tables VIII-XI are obtained in this way.  \nTABLE X. Kerr QNM frequencies corresponding to SXS:BBH:0156 for ( ℓ, m ) = (2 , 2).TABLE XI. Kerr QNM frequencies corresponding to  \nSXS:BBH:1108 for ( ℓ, m ) = (2 , 2).', 'Appendix B: Supplemental plots': 'In this appendix, we supplement additional figures which are not included in the main text. Figure 43 shows the results of the fit of the mock waveform with constant Ψ ( c ) [0 , 7] considered in § IV, and Fig. 44 shows those after the subtraction of the fundamental mode. Comparing them to Figs. 3-6 and Figs. 8-11 for the case of the mock waveform without constant, we see that the constant reduces the stability of the fit.  \nFigures 45-47 show the results of the fit of SXS:BBH:0158, SXS:BBH:0156, and SXS:BBH:1108, respectively, with Ψ SXS -c (left column), and Ψ SXS -c -ψ iter [0 , 0] (right column) after the subtraction of the fundamental mode. By comparing with the results for SXS:BBH:0158 and the ones for SXS:BBH:0305 in Figs. 26 and 27, we can examine how the stability of the fit is changed for the case of high spin. Specifically, it can be seen that the start time of the fit t i at which the mismatch value reaches the lower bound is generally later. In particular, comparing the fit of Ψ SXS -c with N = 7, for SXS:BBH:0305 the mismatch is quite small and takes the minimum around ( t i -t peak ) /M = 0, whereas for SXS:BBH:0158, the mismatch has not yet fallen completely around ( t i -t peak ) /M = 0 but around ( t i -t peak ) /M = 10. To make the mismatch for SXS:BBH:0158 take the minimum at the fit starting from the peak time, it would be necessary to superpose more  \novertones in the fitting function. However, while the mismatch would be reduced, such higher overtones cannot be extracted in a stable manner but would end up with overfit. As argued in the main text, the plateau of A n and ϕ n is extended than that in SXS:BBH:0305. Regarding the rate of change, for the fundamental mode (blue), γ 0 remain the order of 10 -4 for both cases, whereas for the overtones, the region where γ 1 remains small is extended. The left column of Fig. 45 show the result of the fit after the subtraction of the fundamental mode, where we can observe a similar change from Fig. 26 to 27.  \nOn the other hand, from the results for the low spin cases depicted in Figs. 46 and 47, we see that that the trend is opposite to the case of the high spin of SXS:BBH:0158. For instance, the mismatch is already reduced and reaches the lower bound at t peak with a fewer number of overtones. The plateau is short, and the fit of the low spin event is more challenging.', 'Appendix C: Estimation of numerical errors': "In this appendix we explain how we estimate the numerical errors used in the main text. When using a SXS simulation, we adopt the highest resolution level and the extrapolation order N ext = 2 as the base waveform Ψ (0) . To estimate the numerical errors, we compare the base waveform with the following two kinds of waveforms: Ψ (1) with the next-highest resolution level and N ext = 2, and Ψ (2) with the highest resolution level and N ext = 3. However, we should bear in mind that the higher extrapolation order is not necessarily good especially for the analysis of the ringdown [44], so the comparison with the waveform Ψ (2) is treated as a conservative estimation of the numerical errors.  \nTo estimate the numerical errors between Ψ (0) and Ψ ( j ) with j = 1 , 2, we closely follow the procedure described in Appendix B in [44]. First, we truncate the first 500 M of the waveform data to avoid the initial transients. Then, we define the time t begin and t end as the earliest and latest time that is covered by the both waveforms of interest. We would like to compare the two waveforms in this time range, but we should note that in general a mismatch between two waveforms could originate merely from an overall coordinate rotation and/or an overall time shift. Therefore, we choose an coordinate system, where the mismatch between the two waveforms is minimum. First, we choose a z -axis so that the orbital angular momentum L = µ r × v , where µ is the reduced mass and r and v are relative distance and velocity, lies along the + z direction at a reference time t = t begin +1000 M . We consider the waveforms observed from a particular direction and apply the Planck-taper window function given in Eq. (7) in [86] with t 1 = t begin , t 4 = t end , t 2 being the 10th zerocrossing time of Re (Ψ 22 ) after t = t 1 , and t 3 being 50 M after the peak of the waveform. We then calculate a mismatch between Ψ (0) ( t, φ ) and Ψ ( j ) ( t + δt ( j ) , φ + δφ ( j ) ), and optimize δt ( j ) and δφ ( j ) so that the mismatch is min-  \ne sampling points are different for the two waveforms in general, to evaluate the mismatch, we first interpolate the data, and then perform numerical integration. After optimization, we adopt the difference | Ψ (0) -Ψ ( j ) | as estimation of the numerical errors.  \nWe perform the estimation of numerical errors for various directions of observation. Instead of the latitudelongitude lattice, we adopt a spherical Fibonacci lattice, which is known to provide an almost isotropic distributions [87]. In the spherical coordinates, the k -th direction ( k = -N, · · · , N ) is given by  \nφ = 2 πk τ , θ = π 2 -arcsin ( 2 k 2 N +1 ) (C1)  \nwhere τ = (1+ √ 5) / 2 is the golden ratio. We adopt N = 10 throughout the present paper, i.e, we estimate the numerical errors for the observations from almost evenly distributed 21 directions.  \nAfter generating the spherical Fibonacci lattice, we implement the above process of the choice of the z -axis and the rotation of δφ ( j ) by multiplying the rotation matrices as follows. Let r k be the locations of spherical Fibonacci lattice on the unit sphere, and R a ( α ) be the rotation matrix which rotates a vector by an angle α around the a -axis. First, we define  \nr (0) k = R z ( φ (0) L ) R y ( θ (0) L ) r k , (C2)  \nwhere θ (0) L and φ (0) L be the polar and azimuthal angles of the orbital angular momentum L (0) , respectively. r (0) k yields the spherical Fibonacci lattice with the z -axis in the direction of the orbital angular momentum. We use the polar and azimuthal angles of r (0) k to calculate the strain Ψ (0) . Second, we define  \nr ( j ) k = R z ( φ ( j ) L ) R y ( θ ( j ) L ) R z ( δφ ( j ) ) r k , (C3)  \nwhere θ ( j ) L and φ ( j ) L be the polar and azimuthal angles of the orbital angular momentum L ( j ) , respectively. Compared to r (0) k , r ( j ) k yields the lattice rotated by the angle δφ ( j ) around the direction of the orbital angular momentum. We use the polar and azimuthal angles of r ( j ) k and the time shifted by δt ( j ) to calculate the strain Ψ ( j ) . We then optimize δφ ( j ) and δt ( j ) to minimize the mismatch between Ψ (0) and Ψ ( j ) , after which we evaluate the difference | Ψ (0) -Ψ ( j ) | as the estimation of the numerical errors.  \nHaving said that, since for the four simulations we considered in the main text, the orbital angular momentum is along the + z direction, we did not multiply the rotation matrices R z ( φ (0) L ) R y ( θ (0) L ) and R z ( φ ( j ) L ) R y ( θ ( j ) L ) as they are close to identity matrix.  \nIn Figs. 25, 34, 37, 40, green and magenta thin dotted curves represent the errors estimated by | Ψ (0) -Ψ (1) | and | Ψ (0) -Ψ (2) | , respectively. As expected, the errors estimated by using the higher-order extrapolation is larger  \nFIG. 43. Amplitude A n , phase ϕ n , relative error of C n , and rate of change γ n for the fit of the mock waveform Ψ ( c ) [0 , 7] by the fitting function ψ fit [0 ,N ] .  \n<!-- image -->  \nthan the ones estimated by using the next-highest resolution. Note that the error may be small for observations  \n- [1] H.-P. Nollert, Quasinormal modes: the characteristic 'sound' of black holes and neutron stars , Class. Quant. Grav. 16 (1999) R159.\n- [2] K. D. Kokkotas and B. G. Schmidt, Quasinormal modes of stars and black holes , Living Rev. Rel. 2 (1999) 2 [ arXiv:gr-qc/9909058 ].\n- [3] E. Berti, V. Cardoso and A. O. Starinets, Quasinormal modes of black holes and black branes , Class. Quant. Grav. 26 (2009) 163001 [ arXiv:0905.2975 ].\n- [4] R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory , Rev. Mod. Phys. 83 (2011) 793 [ arXiv:1102.4014 ].\n- [5] S. L. Detweiler, Resonant oscillations of a rapidly rotating black hole , Proc. Roy. Soc. Lond. A 352 (1977) 381.\n- [6] F. Echeverria, Gravitational Wave Measurements of the Mass and Angular Momentum of a Black Hole , Phys. Rev. D 40 (1989) 3194.\n- [7] L. S. Finn, Detection, measurement and gravitational radiation , Phys. Rev. D 46 (1992) 5236 [ arXiv:gr-qc/9209010 ].\n- [8] O. Dreyer, B. J. Kelly, B. Krishnan, L. S. Finn, D. Garrison and R. Lopez-Aleman, Black hole spectroscopy: Testing general relativity through gravitational wave observations , Class. Quant. Grav. 21 (2004) 787 [ arXiv:gr-qc/0309007 ].\n- [9] E. Berti, V. Cardoso and C. M. Will, On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA , Phys. Rev. D 73 (2006) 064030 [ arXiv:gr-qc/0512160 ].\n- [10] H.-P. Nollert, About the significance of quasinormal modes of black holes , Phys. Rev. D 53 (1996) 4397 [ arXiv:gr-qc/9602032 ].\n- [11] E. Barausse, V. Cardoso and P. Pani, Can environmental effects spoil precision gravitational-wave astrophysics? , Phys. Rev. D 89 (2014) 104059 [ arXiv:1404.7149 ].\n- [12] R. G. Daghigh, M. D. Green and J. C. Morey, Significance of Black Hole Quasinormal Modes: A Closer Look , Phys. Rev. D 101 (2020) 104009 [ arXiv:2002.07251 ].\n- [13] J. L. Jaramillo, R. Panosso Macedo and L. Al Sheikh, Pseudospectrum and Black Hole Quasinormal Mode Instability , Phys. Rev. X 11 (2021) 031003 [ arXiv:2004.06434 ].\n- [14] W.-L. Qian, K. Lin, C.-Y. Shao, B. Wang and R.-H. Yue, Asymptotical quasinormal mode spectrum for piecewise approximate effective potential , Phys. Rev. D 103 (2021) 024019 [ arXiv:2009.11627 ].\n- [15] M. H.-Y. Cheung, K. Destounis, R. P. Macedo, E. Berti and V. Cardoso, Destabilizing the Fundamental Mode of Black Holes: The Elephant and the Flea , Phys. Rev. Lett. 128 (2022) 111103 [ arXiv:2111.05415 ].\n- [16] J. L. Jaramillo, R. Panosso Macedo and L. A. Sheikh, Gravitational Wave Signatures of Black Hole Quasinormal Mode Instability , Phys. Rev. Lett. 128  \nfrom particular directions. The envelope of the all error curves is shown as a gray curve, which we adopt as the conservative estimation of the numerical errors.", '(2022) 211102 [ arXiv:2105.03451 ].': "- [17] H.-P. Nollert and R. H. Price, Quantifying excitations of quasinormal mode systems , J. Math. Phys. 40 (1999) 980 [ arXiv:gr-qc/9810074 ].\n- [18] E. Berti, V. Cardoso, M. H.-Y. Cheung, F. Di Filippo, F. Duque, P. Martens et al., Stability of the fundamental quasinormal mode in time-domain observations against small perturbations , Phys. Rev. D 106 (2022) 084011 [ arXiv:2205.08547 ].\n- [19] K. Kyutoku, H. Motohashi and T. Tanaka, Quasinormal modes of Schwarzschild black holes on the real axis , Phys. Rev. D 107 (2023) 044012 [ arXiv:2206.00671 ].\n- [20] LIGO Scientific, Virgo collaboration, B. Abbott et al., Observation of Gravitational Waves from a Binary Black Hole Merger , Phys. Rev. Lett. 116 (2016) 061102 [ arXiv:1602.03837 ].\n- [21] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., Tests of general relativity with GW150914 , Phys. Rev. Lett. 116 (2016) 221101 [ arXiv:1602.03841 ].\n- [22] M. Isi, M. Giesler, W. M. Farr, M. A. Scheel and S. A. Teukolsky, Testing the no-hair theorem with GW150914 , Phys. Rev. Lett. 123 (2019) 111102 [ arXiv:1905.00869 ].\n- [23] LIGO Scientific, Virgo collaboration, R. Abbott et al., Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog , Phys. Rev. D 103 (2021) 122002 [ arXiv:2010.14529 ].\n- [24] J. Calder'on Bustillo, P. D. Lasky and E. Thrane, Black-hole spectroscopy, the no-hair theorem, and GW150914: Kerr versus Occam , Phys. Rev. D 103 (2021) 024041 [ arXiv:2010.01857 ].\n- [25] R. Cotesta, G. Carullo, E. Berti and V. Cardoso, Analysis of Ringdown Overtones in GW150914 , Phys. Rev. Lett. 129 (2022) 111102 [ arXiv:2201.00822 ].\n- [26] M. Isi and W. M. Farr, Revisiting the ringdown of GW150914 , arXiv:2202.02941 .\n- [27] E. Finch and C. J. Moore, Searching for a ringdown overtone in GW150914 , Phys. Rev. D 106 (2022) 043005 [ arXiv:2205.07809 ].\n- [28] S. Ma, L. Sun and Y. Chen, Black Hole Spectroscopy by Mode Cleaning , Phys. Rev. Lett. 130 (2023) 141401 [ arXiv:2301.06705 ].\n- [29] S. Ma, L. Sun and Y. Chen, Using rational filters to uncover the first ringdown overtone in GW150914 , Phys. Rev. D 107 (2023) 084010 [ arXiv:2301.06639 ].\n- [30] M. Crisostomi, K. Dey, E. Barausse and R. Trotta, Neural posterior estimation with guaranteed exact coverage: The ringdown of GW150914 , Phys. Rev. D 108 (2023) 044029 [ arXiv:2305.18528 ].\n- [31] Y.-F. Wang, C. D. Capano, J. Abedi, S. Kastha, B. Krishnan, A. B. Nielsen et al., A frequency-domain perspective on GW150914 ringdown overtone , arXiv:2310.19645 .\n- [32] LIGO Scientific, Virgo collaboration, R. Abbott et al., GW190521: A Binary Black Hole Merger with a  \nFIG. 44. Amplitude A n , phase ϕ n , relative error of C n , and rate of change γ n for the fit of the mock waveform Ψ ( c ) [0 , 7] -C iter 0 ψ 0 by the fitting function ψ fit [1 ,N ] .  \n<!-- image -->  \nFIG. 46. Mismatch M (first row), amplitude A n (second row), phase ϕ n (third row), and rate of change γ n (fourth row) for the fit of SXS:BBH:0156 waveform Ψ SXS -c by the fitting function ψ fit [0 ,N ] (left column), and Ψ SXS -c -ψ iter [0 , 0] after the subtraction of the fundamental mode by the fitting function ψ fit [1 ,N ] (right column).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFIG. 45. Mismatch M (first row), amplitude A n (second row), phase ϕ n (third row), and rate of change γ n (fourth row) for the fit of SXS:BBH:0158 waveform Ψ SXS -c by the fitting function ψ fit [0 ,N ] (left column), and Ψ SXS -c -ψ iter [0 , 0] after the subtraction of the fundamental mode by the fitting function ψ fit [1 ,N ] (right column).  \nTotal Mass of 150 M ⊙ , Phys. Rev. Lett. 125 (2020) 101102 [ arXiv:2009.01075 ].  \n[33] LIGO Scientific, Virgo collaboration, R. Abbott et al., Properties and Astrophysical Implications of the 150 M ⊙ Binary Black Hole Merger GW190521 , Astrophys. J. Lett. 900 (2020) L13 [ arXiv:2009.01190 ].  \n[34] C. D. Capano, M. Cabero, J. Westerweck, J. Abedi, S. Kastha, A. H. Nitz et al., Multimode Quasinormal  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nSpectrum from a Perturbed Black Hole , Phys. Rev. Lett. 131 (2023) 221402 [ arXiv:2105.05238 ].  \n[35] C. D. Capano, J. Abedi, S. Kastha, A. H. Nitz, J. Westerweck, Y.-F. Wang et al., Statistical validation of the detection of a sub-dominant quasi-normal mode in GW190521 , arXiv:2209.00640 .  \n- [36] H. Siegel, M. Isi and W. M. Farr, Ringdown of GW190521: Hints of multiple quasinormal modes with a  \nFIG. 47. Mismatch M (first row), amplitude A n (second row), phase ϕ n (third row), and rate of change γ n (fourth row) for the fit of SXS:BBH:1108 waveform Ψ SXS -c by the fitting function ψ fit [0 ,N ] (left column), and Ψ SXS -c -ψ iter [0 , 0] after the subtraction of the fundamental mode by the fitting function ψ fit [1 ,N ] (right column).  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nprecessional interpretation , Phys. Rev. D 108 (2023) 064008 [ arXiv:2307.11975 ].  \n- [37] M. Giesler, M. Isi, M. A. Scheel and S. Teukolsky, Black Hole Ringdown: The Importance of Overtones , Phys. Rev. X 9 (2019) 041060 [ arXiv:1903.08284 ].\n- [38] V. Baibhav, M. H.-Y. Cheung, E. Berti, V. Cardoso, G. Carullo, R. Cotesta et al., Agnostic black hole spectroscopy: Quasinormal mode content of numerical  \n- relativity waveforms and limits of validity of linear perturbation theory , Phys. Rev. D 108 (2023) 104020 [ arXiv:2302.03050 ].\n- [39] P. J. Nee, S. H. Volkel and H. P. Pfeiffer, Role of black hole quasinormal mode overtones for ringdown analysis , Phys. Rev. D 108 (2023) 044032 [ arXiv:2302.06634 ].\n- [40] P. Mourier, X. Jim'enez Forteza, D. Pook-Kolb, B. Krishnan and E. Schnetter, Quasinormal modes and their overtones at the common horizon in a binary black hole merger , Phys. Rev. D 103 (2021) 044054 [ arXiv:2010.15186 ].\n- [41] S. Ma, K. Mitman, L. Sun, N. Deppe, F. H'ebert, L. E. Kidder et al., Quasinormal-mode filters: A new approach to analyze the gravitational-wave ringdown of binary black-hole mergers , Phys. Rev. D 106 (2022) 084036 [ arXiv:2207.10870 ].\n- [42] E. Berti and A. Klein, Mixing of spherical and spheroidal modes in perturbed Kerr black holes , Phys. Rev. D 90 (2014) 064012 [ arXiv:1408.1860 ].\n- [43] A. H. Mroue et al., Catalog of 174 Binary Black Hole Simulations for Gravitational Wave Astronomy , Phys. Rev. Lett. 111 (2013) 241104 [ arXiv:1304.6077 ].\n- [44] M. Boyle et al., The SXS Collaboration catalog of binary black hole simulations , Class. Quant. Grav. 36 (2019) 195006 [ arXiv:1904.04831 ].\n- [45] G. B. Cook, Aspects of multimode Kerr ringdown fitting , Phys. Rev. D 102 (2020) 024027 [ arXiv:2004.08347 ].\n- [46] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis. Springer, 2002, 10.1007/978-0-387-21738-3.\n- [47] LIGO Scientific, Virgo collaboration, R. Abbott et al., GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object , Astrophys. J. Lett. 896 (2020) L44 [ arXiv:2006.12611 ].\n- [48] S. L. Detweiler, BLACK HOLES AND GRAVITATIONAL WAVES. III. THE RESONANT FREQUENCIES OF ROTATING HOLES , Astrophys. J. 239 (1980) 292.\n- [49] V. Cardoso, A Note on the resonant frequencies of rapidly rotating black holes , Phys. Rev. D 70 (2004) 127502 [ arXiv:gr-qc/0411048 ].\n- [50] S. Hod, Slow relaxation of rapidly rotating black holes , Phys. Rev. D 78 (2008) 084035 [ arXiv:0811.3806 ].\n- [51] H. Yang, F. Zhang, A. Zimmerman, D. A. Nichols, E. Berti and Y. Chen, Branching of quasinormal modes for nearly extremal Kerr black holes , Phys. Rev. D 87 (2013) 041502 [ arXiv:1212.3271 ].\n- [52] H. Yang, A. Zimmerman, A. Zengino˘glu, F. Zhang, E. Berti and Y. Chen, Quasinormal modes of nearly extremal Kerr spacetimes: spectrum bifurcation and power-law ringdown , Phys. Rev. D 88 (2013) 044047 [ arXiv:1307.8086 ].\n- [53] N. Sago, S. Isoyama and H. Nakano, Fundamental Tone and Overtones of Quasinormal Modes in Ringdown Gravitational Waves: A Detailed Study in Black Hole Perturbation , Universe 7 (2021) 357 [ arXiv:2108.13017 ].\n- [54] O. M. Moreschi, Supercenter of Mass System at Future Null Infinity , Class. Quant. Grav. 5 (1988) 423.\n- [55] O. M. Moreschi and S. Dain, Rest frame system for asymptotically flat space-times , J. Math. Phys. 39 (1998) 6631 [ arXiv:gr-qc/0203075 ].\n- [56] S. Dain and O. M. Moreschi, General existence proof for rest frame systems in asymptotically flat space-time ,\n- Class. Quant. Grav. 17 (2000) 3663 [ arXiv:gr-qc/0203048 ].\n- [57] K. Mitman et al., Fixing the BMS frame of numerical relativity waveforms , Phys. Rev. D 104 (2021) 024051 [ arXiv:2105.02300 ].\n- [58] L. Maga˜na Zertuche et al., High precision ringdown modeling: Multimode fits and BMS frames , Phys. Rev. D 105 (2022) 104015 [ arXiv:2110.15922 ].\n- [59] E. W. Leaver, Spectral decomposition of the perturbation response of the Schwarzschild geometry , Phys. Rev. D 34 (1986) 384.\n- [60] N. Andersson, Excitation of Schwarzschild black hole quasinormal modes , Phys. Rev. D 51 (1995) 353.\n- [61] N. Andersson, Evolving test fields in a black hole geometry , Phys. Rev. D 55 (1997) 468 [ arXiv:gr-qc/9607064 ].\n- [62] K. Glampedakis and N. Andersson, Quick and dirty methods for studying black hole resonances , Class. Quant. Grav. 20 (2003) 3441 [ arXiv:gr-qc/0304030 ].\n- [63] E. Berti and V. Cardoso, Quasinormal ringing of Kerr black holes. I. The Excitation factors , Phys. Rev. D 74 (2006) 104020 [ arXiv:gr-qc/0605118 ].\n- [64] Z. Zhang, E. Berti and V. Cardoso, Quasinormal ringing of Kerr black holes. II. Excitation by particles falling radially with arbitrary energy , Phys. Rev. D 88 (2013) 044018 [ arXiv:1305.4306 ].\n- [65] M. H.-Y. Cheung, E. Berti, V. Baibhav and R. Cotesta, Extracting linear and nonlinear quasinormal modes from black hole merger simulations , Phys. Rev. D 109 (2024) 044069 [ arXiv:2310.04489 ].\n- [66] R. J. Gleiser, C. O. Nicasio, R. H. Price and J. Pullin, Second order perturbations of a Schwarzschild black hole , Class. Quant. Grav. 13 (1996) L117 [ arXiv:gr-qc/9510049 ].\n- [67] K. Ioka and H. Nakano, Second and higher-order quasi-normal modes in binary black hole mergers , Phys. Rev. D 76 (2007) 061503 [ arXiv:0704.3467 ].\n- [68] H. Nakano and K. Ioka, Second Order Quasi-Normal Mode of the Schwarzschild Black Hole , Phys. Rev. D 76 (2007) 084007 [ arXiv:0708.0450 ].\n- [69] S. Okuzumi, K. Ioka and M.-a. Sakagami, Possible Discovery of Nonlinear Tail and Quasinormal Modes in Black Hole Ringdown , Phys. Rev. D 77 (2008) 124018 [ arXiv:0803.0501 ].\n- [70] D. Brizuela, J. M. Martin-Garcia and M. Tiglio, A Complete gauge-invariant formalism for arbitrary second-order perturbations of a Schwarzschild black hole , Phys. Rev. D 80 (2009) 024021 [ arXiv:0903.1134 ].\n- [71] E. Pazos, D. Brizuela, J. M. Martin-Garcia and M. Tiglio, Mode coupling of Schwarzschild perturbations: Ringdown frequencies , Phys. Rev. D 82 (2010) 104028 [ arXiv:1009.4665 ].\n- [72] N. Loutrel, J. L. Ripley, E. Giorgi and F. Pretorius, Second Order Perturbations of Kerr Black Holes: Reconstruction of the Metric , Phys. Rev. D 103 (2021) 104017 [ arXiv:2008.11770 ].\n- [73] J. L. Ripley, N. Loutrel, E. Giorgi and F. Pretorius,  \n- Numerical computation of second order vacuum perturbations of Kerr black holes , Phys. Rev. D 103 (2021) 104018 [ arXiv:2010.00162 ].\n- [74] M. Lagos and L. Hui, Generation and propagation of nonlinear quasinormal modes of a Schwarzschild black hole , Phys. Rev. D 107 (2023) 044040 [ arXiv:2208.07379 ].\n- [75] L. London, D. Shoemaker and J. Healy, Modeling ringdown: Beyond the fundamental quasinormal modes , Phys. Rev. D 90 (2014) 124032 [ arXiv:1404.3197 ].\n- [76] L. Sberna, P. Bosch, W. E. East, S. R. Green and L. Lehner, Nonlinear effects in the black hole ringdown: Absorption-induced mode excitation , Phys. Rev. D 105 (2022) 064046 [ arXiv:2112.11168 ].\n- [77] M. H.-Y. Cheung et al., Nonlinear Effects in Black Hole Ringdown , Phys. Rev. Lett. 130 (2023) 081401 [ arXiv:2208.07374 ].\n- [78] K. Mitman et al., Nonlinearities in Black Hole Ringdowns , Phys. Rev. Lett. 130 (2023) 081402 [ arXiv:2208.07380 ].\n- [79] J. Redondo-Yuste, G. Carullo, J. L. Ripley, E. Berti and V. Cardoso, Spin dependence of black hole ringdown nonlinearities , Phys. Rev. D 109 (2024) L101503 [ arXiv:2308.14796 ].\n- [80] B. Bucciotti, A. Kuntz, F. Serra and E. Trincherini, Nonlinear quasi-normal modes: uniform approximation , JHEP 12 (2023) 048 [ arXiv:2309.08501 ].\n- [81] S. A. Hughes, A. Apte, G. Khanna and H. Lim, Learning about black hole binaries from their ringdown spectra , Phys. Rev. Lett. 123 (2019) 161101 [ arXiv:1901.05900 ].\n- [82] H. Lim, G. Khanna, A. Apte and S. A. Hughes, Exciting black hole modes via misaligned coalescences: II. The mode content of late-time coalescence waveforms , Phys. Rev. D 100 (2019) 084032 [ arXiv:1901.05902 ].\n- [83] N. Oshita, Ease of excitation of black hole ringing: Quantifying the importance of overtones by the excitation factors , Phys. Rev. D 104 (2021) 124032 [ arXiv:2109.09757 ].\n- [84] H. Lim, S. A. Hughes and G. Khanna, Measuring quasinormal mode amplitudes with misaligned binary black hole ringdowns , Phys. Rev. D 105 (2022) 124030 [ arXiv:2204.06007 ].\n- [85] L. C. Stein, qnm: A Python package for calculating Kerr quasinormal modes, separation constants, and spherical-spheroidal mixing coefficients , J. Open Source Softw. 4 (2019) 1683 [ arXiv:1908.10377 ].\n- [86] D. J. A. McKechan, C. Robinson and B. S. Sathyaprakash, A tapering window for time-domain templates and simulated signals in the detection of gravitational waves from coalescing compact binaries , Class. Quant. Grav. 27 (2010) 084020 [ arXiv:1003.2939 ].\n- [87] A. Gonz'alez, Measurement of Areas on a Sphere Using Fibonacci and Latitude-Longitude Lattices , Math. Geosci. 42 (2010) 49 [ arXiv:0912.4540 ]."}