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

Abstract

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

1. Introduction

The remnant of merging black holes (BHs) is a perturbed compact object (a Kerr BH [1], or even something more exotic) characterized by the set of complex frequencies known as quasinormal modes (QNMs), and the gravitational radiation from this remnant is called the ringdown phase of gravitational waves (GWs). Ringdown GWs from binary BH mergers are now routinely detected by ground-based GW detectors, such as Advanced LIGO and Virgo [2,3]. An accurate measurement of QNMs encoded in the ringdown signal, therefore, offers various GR tests regarding the compact object, for example, to disclose the remnant property (the ergo region of Kerr geometry [4,5,6], ringdown test of general relativity (GR) [7], etc.), and to verify GR itself in the strong-field regime; we refer readers to Ref. [8] for a review of the BH QNMs, and Ref. [9] for that of the ringdown GWs.

The theoretical investigation of QNMs has a long history since early 1970, including Vishveshwara [10], Press [11], Teukolsky and Press [12], and Chandrasekhar and Detweiler [13], and Detweiler [14]. In the 1980s, Detweiler [15] basically completed the analysis of QNMs for Kerr BHs. Leaver [16] gave a standard method to calculate the QNM frequencies very accurately (also see, e.g., Refs. [17,18,19], for more modern techniques to compute QNMs, motivated by the recent development in the BH perturbation theory or high-energy physics).

At the same time, the work on the data analysis of ringdown GWs also goes back to the 1980s. Echeverria [20] and Finn [21] showed that one can indeed extract information of the mass and spin of the remnant Kerr BH from the ringdown signals. Flanagan and Hughes [22,23] evaluated the signal-to-noise ratio (SNR) and parameter estimation errors for practical GW detectors. Creighton [24] analyzed data from the Caltech 40-m prototype interferometer with a single-filter search. Arnaud et al. [25] and Nakano et al. [26,27] developed effective search methods of ringdown GWs. Tsunesada et al. [28,29,30] discussed the detection efficiency, event selection, and parameter estimation of TAMA 300 (‘prototype’ KAGRA) data [31] in the matched filtering analysis. The LIGO-Virgo data in 2005 and the LIGO-Virgo data in 2005–2010 were analyzed in Refs. [32,33], respectively (also see Refs. [34,35,36,37,38] for the earlier related works).

Among that orderly development, the first direct GW detection of merging binary BHs: GW150914 [39] has particularly burst the activity of QNMs and the ringdown GW data analysis. In Ref. [40], three tests of GR involving the ringdown phase were attempted for GW150914: consistency test between the remnant BH parameters evaluated by the inspiral and post-inspiral signals; consistency check between the remnant BH parameters predicted from an inspiral-merger-ringdown waveform in GR and evaluated by the least-damped QNM; and constraining deviations from GR with a parametrized waveform model. This analysis was based on the leading order, longest-lived QNM, i.e., the $(\ell =2,m=2,n=0)$ mode; here, $(\ell ,m)$ are angular indices, and n shows the overtone index. (From this point onward, we refer to $n=0$ mode as the fundamental tone, while we refer to $n\ge 1$ mode as the overtone.), and it was a mainstream strategy for the ringdown GW data analysis (see, e.g., Ref. [41], in the context of multiband GW observation with B-DECIGO [42,43,44,45]).

Nevertheless, Berti et al. [46] showed that the ringdown analysis with only $(\ell =2,$$m=2,n=0)$ mode can lose $10\%$ of potential LIGO events, and the recent development has indeed revealed the importance of including QNM overtones in the ringdown GW analysis. Notably, Giesler et al. [47] demonstrated that one can have a larger SNR for the ringdown phase and better parameter estimations of the remnant BH with enough QNM tones up to $n=7$ (also see Refs. [48,49,50] for earlier work). The study was followed up with more detailed theoretical investigation [51,52,53,54,55,56,57,58] to consider the QNM fits including overtones to the numerical relativity (NR) waveforms of binary BH mergers [59,60,61] after the time of peak amplitude (of some waveform quantity). The observational impacts of including QNM overtones in the ringdown GW analysis are examined for GW150914 [62] and for GW190521 [63], but the strong contribution of overtones and higher harmonics was not found in the ringdown signal. Because GW190521 with a total mass of $150M_{\odot }$ is the heaviest binary BH merger observed to date and its signal are dominated by the merger and ringdown phases, three ringdown models were particularly considered in this case [64]: a single damped sinusoid, a set of the $(\ell =2,m=2)$ QNM fundamental tone and overtones, i.e., $n=0,1,2$, and a set of all fundamental tone of the QNMs up to the $\ell =4$ harmonic mode by taking the QNM starting time into account (also see Ref. [65] about testing GR with ringdown GWs for binary BH events in GWTC-2 [2]).

To data, the contribution of QNM overtones to the ringdown has been largely investigated relying on the NR waveforms as a reference signal. Although NR waveforms ‘exactly’ describe the ringdown phase of the BH mergers, still, the relative error in the NR strain data would be typically limited at ∼${10}^{-4}$ around the time of peak amplitude (see Figure 2 of Ref. [47]), according to a publicly available NR catalogue of “SXS Gravitational Waveform Database” [66]; also see other NR catalogues including “CCRG@RIT Catalog of Numerical Simulations” [67], “Georgia tech catalog of gravitational waveforms” [68], and “SACRA Gravitational Waveform Data Bank” [69].

Interestingly, it is found that the dominant $(\ell =2,m=2)$ mode after the time of the peak amplitude could be fitted by only the QNM fundamental tone and overtones in the linearized perturbation around the remnant BH spacetime, making use of precise NR waveforms [70,71,72]. This facts implies that the second and higher order QNMs would not be excited enough, at least in the $(\ell =2,m=2)$ mode (we note, however, that the second order QNMs due to the self-coupling of the first-order, e.g., $(\ell =2,m=2)$ mode was already identified in the $(\ell =4,m=4)$ mode of the NR waveforms [73]; also see, e.g., Refs. [74,75,76], for earlier works on the second order QNMs). This motivates us to work with the close-limit approximation to the binary BH merges, proposed in the seminal work by Price and Pullin [77], and Abrahams and Prices [78].

The close-limit approximation describes the last stage of a binary BH coalescence as a perturbation of a single (remnant) BH, demonstrating a remarkably good agreement with full NR simulation [79,80]; other recent applications of close-limit approximations are, for example, for the initial data of the NR simulation [81] and for binary BH mergers beyond GR [82]. A particularly relevant aspect of the closed-limit approximation here is that the approximated system can be evolved highly accurately in the framework of the linear BH perturbation theory. A range of methods for numerically evolving the linearized Einstein equation is described in literature, and the numerical accuracy that we can achieve in this work is at the level of ${10}^{-7}$ in the GW strain data at the time of the peak amplitude and ${10}^{-13}$ in the late time of simulations (see Figure 2). Therefore, our close-limit waveforms in the BH perturbation approach allow us the more depth study regarding to the QNM fit than using NR waveform with numerically accuracy available for now.

This paper is organized as follows. In Section 2, we briefly review the close-limit approximation of the post-Newtonian (PN) metric [83], which we use to prepare the initial data for the BH perturbation calculation. In the calculation presented in this work, for simplicity, we focus on the head-on collision of non-spinning binary BHs. Therefore, the remnant BH can be considered as a (non-spinning) Schwarzschild BH. Our numerical method in the BH perturbation approach is presented in Section 3, where we carefully check the numerical accuracy. In Section 4, we propose a modified QNM fitting formula by introducing an orthonormal set of mode function to analyze the ringdown waveforms. In Section 5, we use the QNM fitting formulae in Section 4 and apply them to the GW data obtained in Section 3. We will find the late-time power-law tail [84] in the fit residuals, as well as a convergent behavior for the modified QNM fitting formula. In Section 6, we summarize our analysis, and discuss some applications. A brief analysis of the dominant $(\ell =2,m=2)$ mode of NR waveforms for non-precessing binary BHs and the fitting with the QNMs overtones is given in Appendix A and Appendix B. Here, we find some problematic behavior in NR waveforms for the cases with a highly spinning remnant BH, and this is discussed in Appendix C.

Throughout this paper, we adopt the negative metric signature $(-+++)$, geometrized units with $c=G=1$, except in Section 2. The coordinates, $x^{\alpha }=\{ct,r,\theta ,\phi \}$, are used as the Schwarzschild coordinates for Schwarzschild spacetime and the Boyer-Lindquist coordinates for Kerr spacetime.

2. Post-Newtonian Initial Conditions in the Close Limit: Head-On Collisions

We wish to describe the evolution of non-spinning binary BHs in its last stage, making use of the close-limit approximation starting from a given initial condition. For this purpose, we adopt the close-limit form of the post-Newtonian metric for two non-spinning BHs (or point particles) of masses $m_A(A=1,2)$ derived by Le Tiec and Blanchet [83]. This close-limit, PN metric is particularly convenient to examine the ringdown phase of asymmetric mass-ratio binaries. To set the stage for our analysis, we review their essential discussion and results here.

We assume that the binary system admits two different (dimensionless) small parameters. The first one is the close-limit parameter [77,78,80,85,86,87] in which the binary separation $r_{12}$ is sufficiently small compared to the distance r from any field point to a reference source point (e.g., the center of mass of the system):

$${ϵ}_{\mathrm{CL}}\sim \frac{r_{12}}{r}\ll 1.$$

(1)

The other is the post-Newtonian (PN) parameter [88,89,90,91] where the typical value of the relative orbital velocity $v_{12}$ is small enough compared to the speed of light (or that of the separation between two masses $r_{12}$ is sufficiently large):

$${ϵ}_{\mathrm{PN}}\sim \frac{v_{12}^2}{c^2}\sim \frac{GM}{c^2{r}_{12}}\ll 1.$$

(2)

Here, $M=m_1+m_2$ is the total mass of the binary system, and we shall say that a term relative to $O\left({ϵ}_{\mathrm{PN}}^N\right)$ is Nth PN order. We note that the PN expansion is valid when $r\ll r_{12}/\sqrt{GM/\left(c^2{r}_{12}\right)}$ that defines the (so-named) near-zone.

The crux of this approach is that the metric of the binary system is formally expanded in powers of both in ${ϵ}_{\mathrm{CL}}$ and ${ϵ}_{\mathrm{PN}}$ (though the close-limit approximation would imply $r_{12}\gtrsim GM/c^2$, i.e., ${ϵ}_{\mathrm{PN}}\lesssim 1$ [84]; also see Refs. [83,92,93] for an elaboration of this double expansion). This allows us to recast the (late inspiral phase of) PN binary spacetime to a PN vacuum perturbation of the single Schwarzschild spacetime. Specifically, we work with the 2PN near-zone metric $g_{\mu \nu }^{\mathrm{PN}}$ of two point masses [94,95,96,97], and then further expand it in the power of ${ϵ}_{\mathrm{CL}}$. The resultant metric can be manipulated to be identified with

$$g_{\mu \nu }^{\mathrm{PN}}=g_{\mu \nu }^{\mathrm{Schw}.}+h_{\mu \nu },$$

(3)

where $g_{\mu \nu }^{\mathrm{Schw}.}$ is the Schwarzschild metric with mass M, and $h_{\mu \nu }$ is the perturbation that takes the schematic form

$$h_{\mu \nu }=G\sum _{n=0}^2\sum _{k\ge 0}{h}_{\mu \nu }^{(n,k)}{ϵ}_{\mathrm{PN}}^n{ϵ}_{\mathrm{CL}}^{k+1}+O\left(G^2\right).$$

(4)

We note that this expression is truncated at $O\left(G\right)$ so that $h_{\mu \nu }$ consistently satisfies the linearized Einstein equation to the 2PN order. The general explicit expressions for the coefficients $h_{\mu \nu }^{(n,k)}$ are too unwieldy to be presented here; the expressions up to $k=2$ in the (slightly altered) Cartesian harmonic coordinates are listed in Equation (2.12) of Ref. [83], for example.

Because the Schwarzschild spacetime is spherically symmetric, the angular dependence of $h_{\mu \nu }$ in Equation (4) can be conveniently separated to decompose it in a given tensor spherical harmonics. Using a standard Schwarzschild coordinates, this mode decomposition takes a form (e.g., Section 4.2 of Ref. [98])

$$h_{\mu \nu }=\sum _{i=1}^{10}\sum _{\ell ,m}{\overline{h}}^{\left(i\right),\ell m}(t,r)Y_{\mu \nu }^{\left(i\right),\ell m}(\theta ,\phi ),$$

(5)

where $Y_{\mu \nu }^{\left(i\right),\ell m}$ is a basis of tensor spherical harmonics. A set of 10 time-radial functions ${\overline{h}}^{\left(i\right),\ell m}$ further decouples into two subsets: seven for even-parity perturbations, and three for odd-parity ones, depending on whether they are “even” or “odd” under the transformation $(\theta ,\phi )\to (\pi -\theta ,\pi +\phi )$.

Model Problem: Head-On Collisions of Two Non-Spinning Black Holes

We now specialize the close-limit, 2PN metric to a head-on collision of two non-spinning BHs. The head-on collision is not an astrophysically relevant scenario, but it provides a useful and simple example to explore the ringdown radiation with a good numerical accuracy in the remaining sections.

We adopt the notations of Le Tiec and Blanchet (Section 3 of Ref. [83]) and work with the Regge-Wheeler basis of tensor Spherical harmonics and in the Regge-Wheeler gauge to simplify the expressions [99] (also see Appendix A of Ref. [100]). In the case of the head-on collision, the odd-parity perturbations are identically zero, and the non-trivial components of even-parity perturbations (in the Regge-Wheeler gauge) are:

$$\begin{array}{cc}\hfill h_{00}& =\left(1-\frac{2M}{r}\right)\sum _{\ell ,m}{\tilde{H}}_0^{\ell m}{Y}_{\ell m},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill h_{rr}& ={\left(1-\frac{2M}{r}\right)}^{-1}\sum _{\ell ,m}{\tilde{H}}_2^{\ell m}{Y}_{\ell m},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill h_{\phi \phi }& =h_{\theta \theta }{sin}^2\theta =r^2\sum _{\ell ,m}{\tilde{K}}^{\ell m}{sin}^2\theta Y_{\ell m}.\hfill \end{array}$$

(8)

We then read off ${\tilde{H}}_0^{\ell m},{\tilde{H}}_2^{\ell m}$ and ${\tilde{K}}^{\ell m}$ from the close-limit expansion to the 2PN metric in Equation (4). To model the ringdown radiation yet in the simple setup, we shall consider only the $\ell =2$ and 3 perturbations that satisfy the linearized Einstein equations at the first order in the (symmetric) mass ratio $\nu \equiv m_1{m}_2/M^2$ up to terms $O(G^2,c^{{}^6},r_{12}^4)$. This restriction gives

$${\tilde{H}}_0^{\ell m}\equiv {\tilde{H}}^{\ell m}={\tilde{H}}_2^{\ell m}.$$

(9)

With this relation, the explicit expressions of $\ell =2$ components on the time symmetric initial surface are (Equation (3.6) of Ref. [83]; $r_{12}$ is the binary’s separation, and we assume that the relative orbital velocity $v_{12}$, orbital angular velocity ${\omega }_{12}$, and the orbital phase $\beta$ are all zero):

$$\begin{array}{cc}\hfill {\tilde{H}}_{\mathrm{CL}}^{2,0}=& -2\sqrt{\frac{\pi }{5}}\nu \frac{M{r}_{12}^2}{r^3}={\tilde{K}}_{\mathrm{CL}}^{2,0},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\tilde{H}}_{\mathrm{CL}}^{2,\pm 2}=& -\sqrt{\frac{3}{2}}{\tilde{H}}_{\mathrm{CL}}^{2,0},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\tilde{K}}_{\mathrm{CL}}^{2,\pm 2}=& -\sqrt{\frac{3}{2}}{\tilde{K}}_{\mathrm{CL}}^{2,0},\hfill \end{array}$$

(10)

and those of $\ell =3$ components are:

$$\begin{array}{cc}\hfill {\tilde{H}}_{\mathrm{CL}}^{3,\pm 1}=& \mp \sqrt{\frac{3\pi }{7}}\nu \frac{\delta M{r}_{12}^3}{r^4}={\tilde{K}}_{\mathrm{CL}}^{3,\pm 1},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\tilde{H}}_{\mathrm{CL}}^{3,\pm 3}=& -\sqrt{\frac{5}{3}}{\tilde{H}}_{\mathrm{CL}}^{3,\pm 1},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\tilde{K}}_{\mathrm{CL}}^{3,\pm 3}=& -\sqrt{\frac{5}{3}}{\tilde{K}}_{\mathrm{CL}}^{3,\pm 1},\hfill \end{array}$$

(11)

with $\delta M\equiv m_1-m_2$.

3. Numerical Method

In this section, we shall use the close-limit, 2PN metric in Equations (10) and (11) (in the Regge-Wheeler gauge) as initial data and numerically evolve it with the time-domain (TD) implementation of the Regge-Wheeler-Zerilli equations, to obtain the associated ringdown waveform.

Before proceeding, we note that Equation (10) in the equal-mass limit $(\nu =1/4)$ essentially agrees with the Brill-Lindquist geometry for the head-on collision [101], after a certain coordinate adjustment to the initial separation $r_{12}$. This will provide reassurance that the physics in our TD numerical waveform should be largely independent of our specific choice of the PN initial data (see Section 5 of Ref. [83] for further elaboration).

3.1. Initial Data and Gravitational Waveform in Terms of Master Functions

In the practical implementation, it is more preferable to use gauge-invariant master functions that can be constructed from ${\tilde{H}}_{\mathrm{CL}}^{\ell m}$ and ${\tilde{K}}_{\mathrm{CL}}^{\ell m}$ in the Regge-Wheeler gauge. For even-parity perturbations, we use (the Regge-Wheeler-gauge expression for) the Zerilli-Moncrief function given by (Equation (5.1) of Ref. [83]; also see Equation (126) of Ref. [98])

$${\Psi }^{\ell m}=\frac{r}{2({\lambda }_{\ell }+1)}\left[{\tilde{K}}^{\ell m}+\frac{r-2M}{{\lambda }_{\ell }r+3M}\left({\tilde{H}}^{\ell m}-r{\partial }_r{\tilde{K}}^{\ell m}\right)\right],$$

(12)

where ${\lambda }_{\ell }\equiv (\ell -1)(\ell +2)/2$. Substituting Equations (10) and (11) into this, the $\ell =2$ and 3 multipolar coefficients for the head-on collisions are then translated to

$$\begin{array}{cc}\hfill {\Psi }^{2,0}& =\frac{1}{3}\sqrt{\frac{\pi }{5}}\nu M\frac{r_{12}^2}{r^2}\frac{5M-6r}{3M+2r},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\Psi }^{3,\pm 1}& =\pm \frac{1}{4}\sqrt{\frac{\pi }{21}}\nu \delta M\frac{r_{12}^3}{r^3}\frac{7M-10r}{3M+5r},\hfill \end{array}$$

(14)

with ${\Psi }^{2,\pm 2}=-\sqrt{3/2}{\Psi }^{2,0}$ and ${\Psi }^{3,\pm 3}=-\sqrt{5/3}{\Psi }^{3,\pm 1}$.

Furthermore, the Zerilli-Moncrief function in Equation (12) is directly related to the GW strain

$$h_{+}-i{h}_{\times }=\sum _{\ell ,m}\frac{\sqrt{(\ell -1)\ell (\ell +1)(\ell +2)}}{r}{\Psi }^{\ell m}(t,r){}_{-2}{Y}_{\ell m}(\theta ,\phi ),$$

(15)

where ${}_{-2}{Y}_{\ell m}$ is the spin-weighted spherical harmonics with $s=-2$ (note that the contribution from the odd-parity perturbation identically vanishes in the case of the head-on collision). To avoid the dependence of the polarization on the observer’s location, we focus on ${\Psi }^{\ell m}$ in the remaining sections.

3.2. Time Domain Integration of the Zerilli Equation

In our model problem of the head-on collision, the Zerilli-Moncrief function ${\Psi }^{\ell m}$ satisfies the homogeneous Zerilli equation. We write the equation in the double null coordinates, in practice,

$$\left[\frac{{\partial }^2}{\partial u\partial v}+\frac{1}{4}{V}_{\ell }^{\mathrm{Z}}\left(r\right)\right]{\Psi }^{\ell m}(u,v)=0,$$

(16)

where $u=t-r^{*}$, and $v=t+r^{*}$, with the tortoise coordinate $r^{*}=r+2Mln[r/\left(2M\right)-1]$, and

$$V_{\ell }^{\mathrm{Z}}=\left(1-\frac{2M}{r}\right)\left[\frac{\ell (\ell +1)}{r^2}-\frac{6M}{r^3}\frac{{\lambda }_{\ell }({\lambda }_{\ell }+2)r^2+3M(r-M)}{{({\lambda }_{\ell }r+3M)}^2}\right]$$

(17)

is the Zerilli potential.

To evolve Equation (16), we use the finite-difference scheme developed by Lousto and Price [102]. We introduce an uniform null grid in the $1+1$ numerical domain as the left panel of Figure 1. Consider a grid cell with the size of $h\times h$ in the domain. Integrating Equation (16) over the cell gives

$${\Psi }_1-{\Psi }_2-{\Psi }_3+{\Psi }_4+\frac{h^2}{8}{V}_{\ell }^{\mathrm{Z}}\left(r_c\right)({\Psi }_2+{\Psi }_3)+O\left(h^4\right)=0,$$

(18)

where $r_c$ is the value of r at the center of the cell, and ${\Psi }_i$ for $i=1,2,3,4$ correspond to the values of $\Psi$ (we have omitted the indices $(\ell ,m)$ here) at the vertices of the cell in the right panel of Figure 1. By using the above equation, we can obtain ${\Psi }_1$ from ${\Psi }_2$, ${\Psi }_3$ and ${\Psi }_4$ with the local error term of $O\left(h^4\right)$,

$${\Psi }_1={\Psi }_2+{\Psi }_3-{\Psi }_4-\frac{h^2}{8}{V}_{\ell }^{\mathrm{Z}}\left(r_c\right)({\Psi }_2+{\Psi }_3)+O\left(h^4\right).$$

(19)

The local error is accumulated through the integration over v and u to give a global error of $O\left(h^2\right)$.

The finite-difference scheme on the initial surface should be derived separately. To implement the time-symmetric initial conditions in Equations (13) and (14), we also assume the time symmetry on the initial surface; therefore, ${\partial }_t{\Psi |}_{t=t_0}=0$, which leads to

$${\Psi }_1={\Psi }_4+O\left(h^3\right),$$

(20)

if the vertices 2 and 3 are on the surface. From this relation, we find that the finite-difference scheme on the initial surface is

$${\Psi }_1=\frac{{\Psi }_2+{\Psi }_3}{2}-\frac{h^2}{16}{V}_{\ell }^{\mathrm{Z}}\left(r_c\right)({\Psi }_2+{\Psi }_3)+O\left(h^3\right).$$

(21)

While the local error term in this case is $O\left(h^3\right)$, it eventually gives the global error of $O\left(h^2\right)$ because the error is accumulated only over the initial surface. Hence, the convergence of our numerical solution through the finite difference scheme in Equations (19) and (21) is quadratic with respect to h (see Figure 2 below).

3.3. Simulation Parameters

In this paper, we calculate the Zerilli-Moncrief function with the close-limit, 2PN initial data for head-on collisions of two non-spinning BHs in the following two cases:

• Case A: an equal mass collision with the initial separation of $r_{12}=3.5M$.
• Case B: an asymmetric mass collision with the mass ratio of $\nu =0.2$ and the initial separation of $r_{12}=4.4M$.

These system parameters are chosen so that they agree with those employed in Ref. [83]. In case A, we focus only on $(\ell =2,m=0)$ mode because $(\ell =2,m=\pm 2)$ modes are just proportional to this mode. Similarly, in case B, we look at only $(\ell =3,m=1)$ modes because $(\ell =3,m=-1)$ differs by only an overall minus sign, and $(\ell =3,m=\pm 3)$ modes are proportional to that mode; recall Equations (13) and (14).

3.4. Code Validation

In our numerical simulations, the initial data given in Section 3.1 with the parameters shown in Section 3.3 is set on the region of $-1000\le r^{*}/M\le 3000$ at $t=t_i=0$, and the Zerilli equation in Equation (16) is integrated within the future domain of dependence of the initial surface (the upper triangular domain in the left panel of Figure 1). The produced gravitational waveform is extracted on $v=v_f=3000M$ as a function of the retarded time $u=t-r^{*}$. We take the resolution as $h/M=(0.08,0.04,0.02,0.01,0.005)$ with the number of grid-point N = (50,000, 100,000, 200,000, 400,000, 800,000), respectively. From this point onward, we may set $M=1$.

To check the convergence of our TD code, we consider the absolute differences between the Zerilli-Moncrief function $\Psi$ (we have omitted the indices $(\ell ,m)$ here) computed with different resolutions,

$$\Delta \Psi =\left|{\Psi }_{h_2}-{\Psi }_{h_1}\right|,$$

(22)

where $h_1$ and $h_2$ denote adjacent resolutions. The result for the $(\ell =2,m=0)$ Zerilli-Moncrief function for case A is shown in Figure 2. We can find that the difference becomes a quarter when the resolution gets a half. This shows the $O\left(h^2\right)$ convergence of the data, as expected. The numerical accuracy with the highest resolution, $h=0.005M$, is roughly estimated as ${10}^{-7}$ at the time of the peak amplitude. In the late time of $(u-u_{\mathrm{peak}})/M\gtrsim 200$, the numerical accuracy becomes ${10}^{-13}$.

4. Modeling of Ringdown Waveforms

In this section, we introduce waveform models to describe the ringdown data in our close-limit calculations, as a superposition of QNMs associated with the remnant Schwarzschild BH. The method to compute QNMs are well developed in literature (see, e.g., Ref. [103]). We use the accurate numerical data of QNM frequencies of the remnant BH provided by Berti [103] (for $\ell =2$) and those generated with the Black Hole Perturbation Club (B.H.P.C.) code [104] (for $\ell =3$).

We should note that the zero-spin remnant BH here is a limitation of our close-limit approximation to the binary non-spinning BHs. When we treat general binary BH mergers in the full NR simulations with an astrophysically realistic initial data, the remnant BHs after merger are Kerr with a non-zero spin in general. For example, a non-spinning equal-mass binary BH of the quasicircular merger can produce a remnant Kerr BH with the nondimensional spin ∼$0.67$ [105].

4.1. Standard Quasinormal-Mode Fitting Formula

We write the ringdown waveforms as a sum of the fundamental tone and overtones of QNMs decomposed into spin-weighted spherical harmonics with angular indices $(\ell ,m)$ and spin $s=-2$, ${}_{-2}{Y}_{\ell m}(\theta ,\phi )$:

$$h=h_{+}-i{h}_{\times }=\frac{1}{r}\sum _{\ell ,m}{h}_N^{\ell m}\left(t\right){}_{-2}{Y}_{\ell m}(\theta ,\phi )(t\ge t_0),$$

(23)

with

$$h_N^{\ell m}\left(t\right)=\sum _{n=0}^N{Q}_n^{\ell m}{e}^{-i{\omega }_{\ell mn}^{\mathrm{QNM}}(t-t_0)},$$

(24)

where we truncate the summation with a finite overtone index of $n=N$, although there are infinite overtones. Here, $Q_n^{\ell m}$ is a complex amplitude, ${\omega }_{\ell mn}^{\mathrm{QNM}}$ is the Schwarzschild QNM frequency of the $(\ell ,m,n)$ mode, and $t_0$ is the ’starting time’ before which we do not include the model to fit our numerical waveforms. As mentioned later, we choose the time of the peak amplitude of the wave, $t_{\mathrm{peak}}$, as $t_0$.

With an appropriate choice of the initial phase, our numerical waveforms of the head-on collisions are given as purely real functions of the retarded time, u. Therefore, instead of the complex representation of the waveform in Equation (24), we adopt the following real expression for the Zerilli-Moncrief function in practice, making use of the relation in Equation (15),

$${\Psi }_N^{\ell m}\left(u\right)=\sum _{n=0}^N\left[C_{2n}^{\ell m}{\varphi }_{2n}^{\ell m}\left(u\right)+C_{2n+1}^{\ell m}{\varphi }_{2n+1}^{\ell m}\left(u\right)\right],$$

(25)

with

$$\begin{array}{cc}\hfill {\varphi }_{2n}^{\ell m}\left(u\right)& =\sqrt{\frac{4({\omega }_{\ell mn}^2{\tau }_{\ell mn}^2+1)}{{\omega }_{\ell mn}^2{\tau }_{\ell mn}^3}}{e}^{-(u-u_0)/{\tau }_{\ell mn}}sin\left[{\omega }_{\ell mn}(u-u_0)\right],\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\varphi }_{2n+1}^{\ell m}\left(u\right)& =\sqrt{\frac{4({\omega }_{\ell mn}^2{\tau }_{\ell mn}^2+1)}{({\omega }_{\ell mn}^2{\tau }_{\ell mn}^2+2){\tau }_{\ell mn}}}{e}^{-(u-u_0)/{\tau }_{\ell mn}}cos\left[{\omega }_{\ell mn}(u-u_0)\right],\hfill \end{array}$$

(27)

where $C_k^{\ell m}$ ($k=0,\cdots ,2N+1$) are real constants, $u_0$ is the starting time in terms of the retarded time, and ${\omega }_{\ell mn}$ and ${\tau }_{\ell mn}$ are derived from the real and imaginary parts of the QNM frequency, respectively,

$$\begin{array}{c}\hfill {\omega }_{\ell mn}=\Re ({\omega }_{\ell mn}^{\mathrm{QNM}}),\end{array}$$

(28)

$$\begin{array}{c}\hfill {\tau }_{\ell mn}=-{\displaystyle \frac{1}{\Im ({\omega }_{\ell mn}^{\mathrm{QNM}})}}.\end{array}$$

(29)

The prefactor in the mode functions ${\varphi }_k^{\ell m}\left(u\right)$ of Equations (26) and (27) is chosen so that they are certainly normalized (see the next subsection). We note that the model in Equation (25) is a linear function with respect to the fitting coefficients $C_k^{\ell m}$. This means that one can obtain the best fit values of $C_k^{\ell m}$ through linear least squares.

4.2. Modified Quasinormal-Mode Fitting Formula with an Orthonormal Set of Mode Functions

The mode functions ${\varphi }_k^{\ell m}\left(u\right)$ defined in the previous subsection are linearly independent but not orthogonal each other. For convenience of our analysis, we introduce an ‘orthonormal set’ of mode functions by using the Gram-Schmidt procedure. For this purpose, we first define the inner product between two real functions of $f\left(u\right)$ and $g\left(u\right)$ given in $u\ge u_0$

$$(f,g)\equiv {\int }_{u_0}^{\infty }f\left(u\right)g\left(u\right)du.$$

(30)

For this definition, we can find that ${\varphi }_k^{\ell m}$ in Equations (26) and (27) is normalized as $({\varphi }_k^{\ell m},{\varphi }_k^{\ell m})=1$.

Using the inner product in Equation (30) and the standard Gram–Schmidt procedure, we obtain the orthonormal set of mode functions ${\tilde{\varphi }}_k^{\ell m}$. Their explicit listing is:

$$\begin{array}{ccc}\hfill {\tilde{\varphi }}_0^{\ell m}& \equiv & {\varphi }_0^{\ell m},\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill {\tilde{\varphi }}_1^{\ell m}& \equiv & \frac{1}{\sqrt{1-{({\tilde{\varphi }}_0^{\ell m},{\varphi }_1^{\ell m})}^2}}\left[{\varphi }_1^{\ell m}-({\tilde{\varphi }}_0^{\ell m},{\varphi }_1^{\ell m}){\tilde{\varphi }}_0^{\ell m}\right],\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill {\tilde{\varphi }}_2^{\ell m}& \equiv & \frac{1}{\sqrt{1-{({\tilde{\varphi }}_1^{\ell m},{\varphi }_2^{\ell m})}^2-{({\tilde{\varphi }}_0^{\ell m},{\varphi }_1^{\ell m})}^2}}\hfill \end{array}$$

$$\begin{array}{ccc}& & \times \left[{\varphi }_2^{\ell m}-({\tilde{\varphi }}_1^{\ell m},{\varphi }_2^{\ell m}){\tilde{\varphi }}_1^{\ell m}-({\tilde{\varphi }}_0^{\ell m},{\varphi }_2^{\ell m}){\tilde{\varphi }}_0^{\ell m}\right],\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill {\tilde{\varphi }}_3^{\ell m}& \equiv & \frac{1}{\sqrt{1-{({\tilde{\varphi }}_2^{\ell m},{\varphi }_3^{\ell m})}^2-{({\tilde{\varphi }}_1^{\ell m},{\varphi }_3^{\ell m})}^2-{({\tilde{\varphi }}_0^{\ell m},{\varphi }_3^{\ell m})}^2}}\hfill \end{array}$$

$$\begin{array}{ccc}& & \times \left[{\varphi }_3^{\ell m}-({\tilde{\varphi }}_2^{\ell m},{\varphi }_3^{\ell m}){\tilde{\varphi }}_2^{\ell m}-({\tilde{\varphi }}_1^{\ell m},{\varphi }_3^{\ell m}){\tilde{\varphi }}_1^{\ell m}-({\tilde{\varphi }}_0^{\ell m},{\varphi }_3^{\ell m}){\tilde{\varphi }}_0^{\ell m}\right],\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill ⋮& & \end{array}$$

$$\begin{array}{ccc}\hfill {\tilde{\varphi }}_k^{\ell m}& \equiv & {\left[1-\sum _{p=0}^{k-1}{({\tilde{\varphi }}_p^{\ell m},{\varphi }_k^{\ell m})}^2\right]}^{-1/2}\left[{\varphi }_k^{\ell m}-\sum _{p=0}^{k-1}({\tilde{\varphi }}_p^{\ell m},{\varphi }_k^{\ell m}){\tilde{\varphi }}_p^{\ell m}\right].\hfill \end{array}$$

(35)

However, we note that a direct numerical implementation of above construction does not work well because of the accumulation of round-off errors. In practice, we use the modified Gram-Schmidt algorithm [106] to obtain accurately the orthonormal set of mode functions.

From the above orthonormal set of mode functions ${\tilde{\varphi }}_k^{\ell m}$, we can extract the fitting coefficients easily by

$${\tilde{C}}_k^{\ell m}=({\tilde{\varphi }}_k^{\ell m},{\Psi }_{\mathrm{num}}^{\ell m}),$$

(36)

from the numerically generated Zerilli-Moncrief function ${\Psi }_{\mathrm{num}}^{\ell m}$ and construct the modified model waveform as

$${\tilde{\Psi }}_N^{\ell m}\left(u\right)=\sum _{k=0}^{2N+1}{\tilde{C}}_k^{\ell m}{\tilde{\varphi }}_k^{\ell m}\left(u\right),$$

(37)

up to a finite overtone index of $n=N$.

5. Results

We now present the demonstration of the QNM fits including overtones given in Equations (25) and (37) with the close-limit TD waveform for head-on collisions explained in Section 3.

5.1. Quasinormal-Mode Fitting and Residuals

Figure 3 shows the QNM fit including overtones in Equation (37) to the $(\ell =2,m=0)$ mode of the Zerilli-Moncrief function ${\Psi }^{\ell m}$ in case A (top). Similarly, in case B, we show the QNM fits to the $(\ell =2,m=0)$ mode (bottom-left) and to the $(\ell =3,m=1)$ mode (bottom-right), respectively. In each panel, the QNM overtone in the model waveform is truncated at $N=7$, and the time of the peak amplitude, $u_{\mathrm{peak}}$, is chosen as the starting time, $u_0$.

Figure 4 provides the QNM fit residuals of Figure 3 (the blue solid lines). We find that the residuals are all larger than our numerical uncertainty (the green dotted lines), and shows the power-law behavior in the late time ($u-u_{\mathrm{peak}}>100M$, the orange dashed lines). This corresponds to the late-time tail in the evolution of a perturbation of the gravitational field around a remnant Schwarzschild BH, studied in various literature; see, e.g., Refs. [107,108,109], for the Schwarzschild background spacetime, and, e.g., Refs. [110,111], for the Kerr background spacetime. The power-law tail in a sufficiently late time is also investigated using a sophisticated numerical method with compactification along hyperboloidal surfaces [112]. Each of QNM fit residuals deviates from the power-law line when u becomes larger because the null infinity approximation $v\gg u$ is broken (see Section IV in Ref. [113]).

In Figure 4, we also find that the QNM fit residuals in the early-time ($u-u_{\mathrm{peak}}\lesssim 100M$) do not improve even if one includes more overtones into the model waveform in Equation (37). Although we do not have a mathematically rigorous proof, we speculate that this limitation may arise from the QNM fit including overtones (that describes the exponential decay of the perturbation) to the power-law tail component of the perturbation, based on a simple numerical experiment detailed below.

First, we prepare a mock data of the tail component of the perturbation as

$${\Psi }_{\mathrm{tail}}^{\ell m}\left(u\right)=D_{\ell }W\left(u\right){\left(\frac{u-u_{\mathrm{peak}}}{M}\right)}^{-\ell -1},$$

(38)

where $D_{\ell }$ is an arbitrary constant, and we define a window function by

$$W\left(u\right)=\left\{\begin{array}{cc}0\hfill & (u<u_1),\hfill \\ {\displaystyle \frac{{(u-u_1)}^2(2u+u_1-3u_2)}{{(u_1-u_2)}^3}}\hfill & (u_1\le u\le u_2),\hfill \\ 1\hfill & (u_2<u),\hfill \end{array}\right.$$

(39)

with $u_{\mathrm{peak}}<u_1<u_2$. The window function is introduced to remove the divergent behavior of the late-time tail at $u=u_{\mathrm{peak}}$.

Next, we fit this mock data using the QNM model including overtones in Equation (37) in the same manner as we do for the TD data. Figure 5 shows the QNM fit including overtones to the mock tail data and corresponding fit residuals. We see that the residuals in the early-time ($(u-u_{\mathrm{peak}})/M\lesssim 100$) show the similar behaviour in Figure 4. This result supports our expectation about the early-time residual in the TD data.

5.2. Convergence of the Fitting Coefficients

Because the model waveform in Equations (25) and (37) involves many overtones k (or N), it is instructive to check their convergence in terms of fitting coefficients on k, to ask how higher overtones should be included in practice when we wish to model the TD data from the peak time at $u=u_{\mathrm{peak}}$.

In the left panel of Figure 6, where the TD data is for the $(\ell =2,m=0)$ mode in case A, we show the fitting coefficients for the original QNM basis $C_k^{\ell m}$ (the orange filled triangles) and those for the orthonormal set of mode functions ${\tilde{C}}_k^{\ell m}$ (the blue filled circles) as a function of overtones up to $k\le 21$ (i.e., $N=10$). It is found that $|{\tilde{C}}_k^{\ell m}|$ decreases roughly in power law of k, while $|C_k^{\ell m}|$ does not.

Due to the benefit of the orthonormalization of mode functions in Equation (37), it can be conveniently used to assess the fraction of power in each mode:

$${\rho }^2\equiv ({\tilde{\Psi }}_N^{\ell m},{\tilde{\Psi }}_N^{\ell m})=\sum _{k=0}^{2N+1}{|{\tilde{C}}_k^{\ell m}|}^2.$$

(40)

Based on the fact that one more mode is added to the orthonormal set when k increases by two, we introduce an estimator to characterize the match (or overlap) between the TD data and the QNM fit model including overtones at a given n:

$${\Delta }_n\equiv \frac{\sqrt{|{\tilde{C}}_{2n}^{\ell m}{|}^2+{|{\tilde{C}}_{2n+1}^{\ell m}|}^2}}{\rho }.$$

(41)

The right panel of Figure 6 shows ${\Delta }_n$ as a function of mode n. We see that the inclusion of the first overtone ($n=1$) can improve the match by $10\%$, and that the QNM fits up to $n=5$ will suppress the mismatch less than $0.1\%$.

In the top and bottom panels of Figure 7, we show the same figures as for Figure 6 but with the $(\ell =2,m=0)$ and $(\ell =3,m=1)$ modes in case B, respectively. We see the similar convergent behavior in ${\tilde{C}}_k^{\ell m}$ here, and again find that $(n\ge 5)$ higher overtones contributes to the match only at the level of less than $0.1\%$. These results highlight the potential benefit to use the orthonormal set of mode functions in the data analysis of ringdown GWs.

6. Summary and Discussion

Widely recognized now is the importance of having a ringdown waveform model including higher overtones of QNMs (see, e.g., Ref. [47]). With the overtones, one can set the starting time of the ringdown GWs at the time of the peak amplitude much earlier than a time (∼10–$20M$ after the peak time) to obtain unbiased remnant BH parameters (see, e.g., Figure 5 for GW150914 in Ref. [40]). This allows the increased SNR of observed ringdown signals, as well as better parameter estimation of the remnant BH.

In this paper, we examined in detail the ringdown GWs of binary BHs with QNM fits including overtones in Equations (25) and (37), using the accurate close-limit waveform in the BH perturbation theory as a test bed. Our analysis is restricted to the narrow case of head-on collisions of two non-spinning BHs (based on the 2PN initial data and linear BH perturbation theory), but it suffices to highlight some of the key features of the full problem, and we found three main results: (i) we reconfirm the importance of QNMs overtones to fit the ringdown waveforms after the time of the peak amplitude. This agrees with the previous findings in literature [47,72]; (ii) the small contributions of late-time tail exist in the fit residuals at the level of $O\left({10}^{-5}\right)$ or below (see Figure 4); and (iii) the fitting coefficients decay with overtones when one uses the orthonormal set of mode functions in Equation (36) (see Figure 6 and Figure 7).

A natural next step of our study is to assess the merit of the orthonormal set of mode functions in the context of GW data analysis, accounting for the detector’s noise. A standard technique for the analysis of ringdown GWs is the frequency-domain approach (see, e.g., Ref. [114] for various methods to extract the ringdown GWs), and we can directly apply the orthonormal set of mode function presented in Section 4.2 with the (so-named) noise-weighted inner product in the matched filtering method (see, e.g., Ref. [115]). The implementation of this method to the GW data analysis and the analysis of real data from GW detectors will be a future work.

Another interesting future work would be to add BH’s spin to our close-limit waveform. Although it is known for any field points that the $2.5$PN near-zone metric with spin effects [95,96], the numerical computation of either metric or curvature perturbations in Kerr spacetime is particularly challenging in the time domain [116,117,118] and would need much coding efforts. It is probably more viable, instead, to perturbatively include the spin effect into our non-spinning close-limit waveform (e.g., [119]). It may be also possible to calculate the close-limit waveform in Kerr case making use of the standard frequency-domain technique (and taking the extreme mass-ratio limit). Because, in either case, the extension is involved, we leave them for future work.

From the theoretical point of view, it will be useful to include information on excitation coefficients and factors for the fundamental tone and overtones of QNMs derived from the BH perturbation approach (see, e.g., Refs. [120,121,122,123]). In addition, although we focused only on the $(\ell =2,m=0)$ and $(\ell =3,m=1)$ modes in the bulk of this paper (and will focus only on the $(\ell =2,m=2)$ mode in Appendix A), the inclusion of the other subdominant (higher harmonic) modes will be beneficial. In particular, the $(\ell =4,m=4)$ mode of NR waveforms for binary BHs is especially interesting because this mode includes the second order modes composed by the non-linear, self-coupling of first-order $(\ell =2,$$m=2)\times (\ell =2,m=2)$ mode, and it has already been found in NR simulations [73].

Nevertheless, astute readers may ask what extent the results based on the close-limit approximation (to the head-on collisions of two non-spinning BHs) are consistent with the full NR waveforms in more astrophysically relevant scenarios (e.g., coalescences observed by LIGO/Virgo). We conclude our paper to briefly answer this question, analyzing the NR waveform of SXS:BBH:0305 in the Simulating eXtreme Spacetimes (SXS) catalogue [66] with the same approach that we took in Section 5; Appendix A and Appendix B also provide supplemented examples of the analysis for NR waveforms.

Figure 8 shows the $(\ell =2,m=2)$ spheroidal mode of the waveform SXS:BBH:0305 [66,124,125] (see Equation (A1) for the construction of the spheroidal modes from the spherical ones). Here, we set the time of the peak amplitude as the starting time of the QNM fit (see the orange dashed line in the left panel of Figure 8), and we use the QNM fitting formula introduced in Equation (37). Unlike the TD data generated within the close-limit approximation and the BH perturbation approach, we cannot confirm the power-law tail certainty due to some constant fit residuals of ∼${10}^{-5}$ after $t-t_{\mathrm{peak}}\gtrsim 80M$; recall that the tail contribution in Figure 4 is at the level of $\lesssim O\left({10}^{-5}\right)$, which is as large as the remaining residual in the NR waveform. We expect that future NR simulation with higher numerical accuracy (or, e.g., the close-limit approximation to bound-orbit BBH mergers) will provide a more robust answers of the need of the tail contribution in the ringdown waveform modeling, and eventually in the measurement of remnant properties of BHs.

Figure 9 displays the fitting coefficients for the original QNM basis $C_k^{\ell m}$ and those for the orthonormalized QNM basis ${\tilde{C}}_k^{\ell m}$ introduced in Equation (36), again focusing on the $(\ell =2,m=2)$ spheroidal mode of SXS:BBH:0305. We see that the trend of $C_k^{\ell m}$ and ${\tilde{C}}_k^{\ell m}$ are similar to that of the close-limit, TD data displayed in Figure 6 and Figure 7 (also see Appendix B). Namely, ${\tilde{C}}_k^{\ell m}$ converges as increasing the overtone k, while $C_k^{\ell m}$ (from the least square fit by using Equation (25)) does not. The convergence property of ${\tilde{C}}_k^{\ell m}$ supports the argument that the dominant $(\ell =2,m=2)$ mode of NR waveforms after the time of the peak amplitude can be described only by the QNMs in the linearized BH perturbation (see Section 1 for the other confirmations).