Estimation of starting times of quasinormal modes in ringdown gravitational waves with the Hilbert-Huang transform

TL;DR

This paper introduces a new method to accurately estimate the start time, frequency, and amplitude of quasinormal modes in gravitational wave signals from black hole mergers, validated on simulations and real data.

Contribution

It presents a novel technique for identifying the QNM start time and parameters, improving analysis of gravitational wave signals from black hole mergers.

Findings

QNM frequencies match theoretical values within 1% for simulated data.

The method reveals a correlation between QNM start time and black hole spin.

Parameters estimated from GW150914 are consistent with LIGO's results.

Abstract

It is known that a quasinormal mode (QNM) of a remnant black hole dominates a ringdown gravitational wave (GW) in a binary black hole (BBH) merger. To study properties of the QNMs, it is important to determine the time when the QNMs appear in a GW signal as well as to calculate its frequency and amplitude. In this paper, we propose a new method of estimating the starting time of the QNM and calculating the QNM frequency and amplitude of BBH GWs. We apply it to simulated merger waveforms by numerical relativity and the observed data of GW150914. The results show that the obtained QNM frequencies and time evolutions of amplitudes are consistent with the theoretical values within 1% accuracy for pure waveforms free from detector noise. In addition, it is revealed that there is a correlation between the starting time of the QNM and the spin of the remnant black hole. In the analysis of GW150914, we show that the parameters of the remnant black hole estimated through our method are consistent with those given by LIGO and a reasonable starting time of the QNM is determined.

Figures

Tables

Table 1. Table 1: Relative errors of ω R subscript 𝜔 R \omega_{\mathrm{R}} and ω I subscript 𝜔 I \omega_{\mathrm{I}} estimated for noise-free simulated waveforms.

Relative errors [%]	${\omega }_R$	${\omega }_I$
average value	0.32	0.70
worst case	2.19	4.42

Table 2. Table 2: Initial parameters of BBHs of waveforms 0305, 0262 and 0230. Optimal SNRs defined by Eq. ( 16 ) are also shown. In this table, χ 1 ​ z subscript 𝜒 1 𝑧 \chi_{1z} and χ 2 ​ z subscript 𝜒 2 𝑧 \chi_{2z} respectively denote the initial spin components of the primary and the secondary BHs aligned with the orbital angular momentum of the BBH. For these three waveforms, either of BHs does not have the component of the initial spin perpendicular to the orbital angular momentum.

ID	0305	0262	0230
mass ratio $q$	$1.22$	$3$	$1$
initial spin ${\chi }_{1z}$	$0.1$	$-0.6$	$0.8$
initial spin ${\chi }_{2z}$	$-0.09$	$0.0$	$0.8$
SNR $\rho$	$127.2$	$83.07$	$159.6$

Table 3. Table 3: Parameters of source BBHs including t 0 subscript 𝑡 0 t_{0} estimated by applying our method to noise-free waveforms and the estimated values of t 0 subscript 𝑡 0 t_{0} , ω R subscript 𝜔 R \omega_{\mathrm{R}} and ω I subscript 𝜔 I \omega_{\mathrm{I}} for the simulated observed data. From the first to the last column, the waveform ID, the initial mass ratio m 1 / m 2 subscript 𝑚 1 subscript 𝑚 2 m_{1}/m_{2} , the z 𝑧 z -components of the initial dimensionless spins of the primary and secondary BHs, χ 1 ​ z subscript 𝜒 1 𝑧 \chi_{1z} and χ 2 ​ z subscript 𝜒 2 𝑧 \chi_{2z} , respectively, the ratio of the remnant mass M rem subscript 𝑀 rem M_{\mathrm{rem}} to the initial total mass M t subscript 𝑀 t M_{\mathrm{t}} , the remnant dimensionless spin χ rem subscript 𝜒 rem \chi_{\mathrm{rem}} , the real and imaginary parts of the theoretical value of the QNM frequency of the remnant BH, ω R subscript 𝜔 R \omega_{\mathrm{R}} and ω I subscript 𝜔 I \omega_{\mathrm{I}} , respectively, the estimated value of the starting time of the QNM t 0 subscript 𝑡 0 t_{0} for each of noise-free waveforms, the optimal SNR defined by Eq. ( 16 ) for the simulated observed data, assuming z = 0.09 𝑧 0.09 z=0.09 , M t = 60 ​ M ⊙ subscript 𝑀 t 60 subscript 𝑀 direct-product M_{\mathrm{t}}=60M_{\odot} and Advanced LIGO’s design sensitivity, and the estimated values of t ^ 0 subscript ^ 𝑡 0 \hat{t}_{0} , ω ^ R subscript ^ 𝜔 R \hat{\omega}_{\mathrm{R}} , | ω ^ I | subscript ^ 𝜔 I |\hat{\omega}_{\mathrm{I}}| are listed. For each estimated value, the median value is shown with the region in which 68.2% of the estimated values are included.

ID	$m_1/m_2$	${\chi }_{1z}$	${\chi }_{2z}$	$M_{\mathrm{rem}}/M_t$	${\chi }_{\mathrm{rem}}$	${\omega }_R$/rad s-1	$\left|{\omega }_I\right|$/rad s-1	$t_0$/ms	${\rho }_{\mathrm{IMR}}$	${\hat{t}}_0$/ms	${\hat{\omega }}_R$/rad s-1	$\left|{\hat{\omega }}_I\right|$/rad s-1
0209	1	$-0.90$	$-0.50$	0.965	0.46	1584	301	4.88	107.7	${2.2}_{-1.0}^{+1.2}$	${1496}_{-56.3}^{+47.7}$	${254}_{-27.7}^{+28.5}$
0210	1	$-0.90$	$+0.00$	0.961	0.54	1668	297	4.88	114.6	${2.2}_{-1.7}^{+1.7}$	${1578}_{-97.7}^{+53.4}$	${247}_{-42.4}^{+41.2}$
0211	1	$-0.90$	$+0.90$	0.951	0.68	1860	288	4.88	127.7	${2.2}_{-1.5}^{+1.7}$	${1746}_{-120.8}^{+58.5}$	${241}_{-44.6}^{+38.8}$
0212	1	$-0.80$	$-0.80$	0.967	0.43	1553	302	3.91	105.0	${2.4}_{-2.2}^{+1.7}$	${1474}_{-105.0}^{+54.2}$	${254}_{-65.4}^{+44.5}$
0218	1	$-0.50$	$+0.50$	0.952	0.69	1862	288	9.03	128.7	${2.2}_{-1.5}^{+1.7}$	${1746}_{-103.6}^{+59.2}$	${262}_{-51.5}^{+41.1}$
0221	1	$-0.40$	$+0.80$	0.945	0.74	1974	281	5.86	135.7	${2.4}_{-1.7}^{+1.2}$	${1850}_{-117.6}^{+53.2}$	${240}_{-38.3}^{+26.2}$
0223	1	$+0.30$	$+0.00$	0.947	0.73	1946	283	4.39	133.9	${2.2}_{-1.2}^{+1.7}$	${1820}_{-133.2}^{+57.5}$	${253}_{-47.8}^{+42.8}$
0230	1	$+0.80$	$+0.80$	0.911	0.91	2515	235	8.54	159.6	${2.9}_{-1.0}^{+1.2}$	${2300}_{-143.4}^{+76.9}$	${208}_{-32.1}^{+23.4}$
0231	1	$+0.90$	$+0.00$	0.935	0.82	2149	268	6.84	145.4	${2.4}_{-1.2}^{+1.5}$	${1991}_{-127.0}^{+68.7}$	${223}_{-33.1}^{+33.2}$
0232	1	$+0.90$	$+0.50$	0.920	0.88	2390	247	7.08	156.2	${2.9}_{-1.5}^{+1.2}$	${2203}_{-167.5}^{+71.9}$	${213}_{-43.9}^{+27.0}$
0259	2.5	$+0.00$	$+0.00$	0.967	0.58	1699	293	5.37	105.9	${2.0}_{-1.5}^{+1.7}$	${1593}_{-91.8}^{+64.8}$	${234}_{-43.7}^{+40.1}$
0262	3	$-0.60$	$+0.00$	0.978	0.27	1426	303	4.64	83.1	${3.2}_{-2.0}^{+1.5}$	${1352}_{-50.0}^{+57.9}$	${252}_{-43.8}^{+43.7}$
0263	3	$-0.60$	$+0.60$	0.977	0.29	1440	303	5.13	86.1	${2.7}_{-1.7}^{+1.7}$	${1365}_{-46.1}^{+59.1}$	${245}_{-39.7}^{+44.5}$
0264	3	$-0.60$	$-0.60$	0.979	0.25	1412	303	4.64	80.8	${3.2}_{-2.4}^{+1.5}$	${1345}_{-77.5}^{+53.0}$	${253}_{-63.2}^{+45.2}$
0265	3	$-0.60$	$-0.40$	0.979	0.26	1416	303	3.42	80.6	${2.7}_{-2.7}^{+1.7}$	${1342}_{-86.7}^{+56.7}$	${244}_{-63.8}^{+47.2}$
0274	3	$-0.23$	$+0.85$	0.973	0.47	1576	298	4.39	95.4	${2.7}_{-2.2}^{+2.0}$	${1490}_{-83.7}^{+59.6}$	${234}_{-51.5}^{+52.0}$
0276	3	$+0.30$	$+0.00$	0.972	0.53	1637	295	4.64	96.3	${2.4}_{-1.7}^{+1.5}$	${1539}_{-60.0}^{+59.3}$	${235}_{-37.5}^{+44.7}$
0283	3	$+0.30$	$+0.30$	0.965	0.68	1828	285	5.62	107.8	${2.7}_{-1.5}^{+1.2}$	${1710}_{-85.3}^{+59.5}$	${230}_{-40.4}^{+36.3}$
0287	3	$+0.60$	$-0.60$	0.960	0.78	2012	270	5.86	111.7	${2.2}_{-1.5}^{+1.5}$	${1852}_{-136.5}^{+76.5}$	${214}_{-43.9}^{+35.7}$
0289	3	$+0.60$	$+0.00$	0.957	0.80	2054	267	6.35	115.8	${2.2}_{-1.2}^{+1.5}$	${1897}_{-139.5}^{+70.4}$	${216}_{-42.2}^{+35.3}$
0292	3	$+0.73$	$-0.85$	0.956	0.83	2128	259	6.59	114.5	${2.2}_{-1.7}^{+1.5}$	${1955}_{-171.8}^{+73.2}$	${208}_{-42.7}^{+38.1}$
0293	3	$+0.85$	$+0.85$	0.936	0.91	2478	224	7.81	132.1	${2.4}_{-1.2}^{+2.4}$	${2239}_{-216.5}^{+106.0}$	${183}_{-49.4}^{+37.8}$
0297	6.5	$+0.00$	$+0.00$	0.987	0.35	1467	298	3.42	64.6	${2.7}_{-2.4}^{+1.7}$	${1377}_{-82.3}^{+71.6}$	${218}_{-57.4}^{+60.1}$
0303	10	$+0.00$	$+0.00$	0.992	0.26	1400	299	4.15	49.4	${2.2}_{-2.2}^{+2.0}$	${1294}_{-77.6}^{+77.0}$	${200}_{-55.4}^{+68.3}$
0305	1.22	$+0.10$	$-0.09$	0.952	0.69	1871	287	5.62	127.2	${2.4}_{-1.5}^{+1.2}$	${1759}_{-92.0}^{+50.1}$	${246}_{-34.6}^{+30.3}$
0317	3.33	$+0.30$	$-0.02$	0.964	0.75	1937	275	6.10	107.2	${2.0}_{-1.2}^{+1.5}$	${1803}_{-95.0}^{+62.5}$	${218}_{-34.4}^{+37.0}$

Table 4. Table 4: Properties of filters we applied to the observed data of GW150914. These are digital filters based on the Butterworth infinite impulse response filter of order 4. In this table, f central subscript 𝑓 central f_{\mathrm{central}} , Δ ​ f pass Δ subscript 𝑓 pass \Delta f_{\mathrm{pass}} and Δ ​ f stop Δ subscript 𝑓 stop \Delta f_{\mathrm{stop}} denote the central frequency, pass bandwidth and stop bandwidth of a notch filter, respectively.

	Properties of notch filters [Hz]
	$f_{\mathrm{central}}$	$\Delta f_{\mathrm{pass}}$	$\Delta f_{\mathrm{stop}}$
Notch 1	36	1	0.1
Notch 2	40	1	0.1
Notch 3	60	1	0.1
Notch 4	331	10	1
	Cutoff frequencies [Hz]
		Lower	Upper
Bandstop		330.8	671.6
Bandpass		20.0	320.0

Table 5. Table 5: Estimated parameters of GW150914 by our method. The dimensionless spin parameter χ rem subscript 𝜒 rem \chi_{\mathrm{rem}} and the mass M rem subscript 𝑀 rem M_{\mathrm{rem}} of the remnant BH are estimated from the QNM frequency ω = ω R + i ​ ω I 𝜔 subscript 𝜔 R i subscript 𝜔 I \omega=\omega_{\mathrm{R}}+\mathrm{i}\omega_{\mathrm{I}} by using the conversion table  Berti  et al. ( 2006 ); Ber , assuming the redshift z = 0.09 𝑧 0.09 z=0.09 and the QNM mode of l = m = 2 𝑙 𝑚 2 l=m=2 and n = 0 𝑛 0 n=0 . The error ranges of ω R subscript 𝜔 R \omega_{\mathrm{R}} and ω I subscript 𝜔 I \omega_{\mathrm{I}} are evaluated from only the standard deviation of the fitting of the IF and IA. The measures of errors of χ rem subscript 𝜒 rem \chi_{\mathrm{rem}} and M rem subscript 𝑀 rem M_{\mathrm{rem}} corresponding the errors of ω R subscript 𝜔 R \omega_{\mathrm{R}} and ω I subscript 𝜔 I \omega_{\mathrm{I}} are also shown.

${\hat{t}}_0$/ms	${\hat{\omega }}_R$/rad s-1	$\left|{\hat{\omega }}_I\right|$/rad s-1	${\hat{\chi }}_{\mathrm{rem}}$	${\hat{M}}_{\mathrm{rem}}$/$M_{\odot }$
$5.223$	$1670\pm 9.96$	$211.2\pm 3.13$	${0.8151}_{-0.096}^{+0.092}$	${71.65}_{-0.75}^{+0.77}$

Equations

$$h_{\mathrm{QNM}}(t)$$

$$\displaystyle=A_{0}\mathrm{e}^{\omega_{\mathrm{I}}(t-t_{0})}\cos\bigl{(}\omega_{\mathrm{R}}(t-t_{0})+\varphi_{0}\bigr{)},$$

$$s(t)$$

$$=\sum _{i=1}^{N_{\mathrm{IMF}}}{a}_i(t)\cos ({\varphi }_i(t))+r(t),$$

$$f_i(t)=\frac{1}{2\pi }\frac{\mathrm{d}{\varphi }_i(t)}{\mathrm{d}t}.$$

$$a_j(t)$$

$${\varphi }_j(t)$$

$$\mathrm{RMSE}(n_0,N)=\sqrt{\frac{1}{N}\sum _{n=n_0}^{n_0+N-1}{\left(\ln a_1[n]-bt[n]-c\right)}^2},$$

$${\hat{n}}_0(N)=\underset{n_0}{\mathrm{argmin}}\mathrm{RMSE}(n_0,N).$$

$$\hat{N}=\underset{N}{\mathrm{argmin}}\left[\mathrm{Err}(N_{\min },N)+\mathrm{Err}(N+1,N_{\max })\right],$$

$$\mathrm{Err}(N_1,N_2)=\underset{a,b}{\min }\sqrt{\frac{{\sum }_{N=N_1}^{N_2}{\left(e(N)-(aN+b)\right)}^2}{N_2-N_1}},$$

$$e(N)=\mathrm{RMSE}({\hat{n}}_0(N),N).$$

$$\ln a_1(t)={\omega }_{\mathrm{I}}t-{\omega }_{\mathrm{I}}{t}_0+\ln A_0,$$

$${\hat{\omega }}_{\mathrm{I}}$$

$${\hat{\omega }}_{\mathrm{R}}$$

$${\hat{\omega }}_{\mathrm{R}}^{src.}$$

$$\rho =2{\left({\int }_{f_{\min }}^{f_{\max }}\frac{|\tilde{h}(f){|}^2}{S_{\mathrm{n}}(f)}\mathrm{d}f\right)}^{1/2}.$$

Full text

I Introduction

Advanced LIGO detected two events of gravitational waves (GWs) in 2015, both of which are generated by binary black hole (BBH) mergers. The waveforms of GWs from BBHs consist of three phases, namely, the inspiral, merger and ringdown phases. The inspiral phase corresponds to an orbital motion of the two black holes (BHs) and can be described by the post-Newtonian approximation of general relativity (GR). The merger phase is described as a full nonlinear evolution of Einstein equations and derived by numerical relativity simulations. The ringdown phase corresponds to quasinormal modes (QNMs) of the remnant BH. Since the QNMs are intrinsic in the BH metric and can be seen only with GWs, analyses of the QNMs with GWs enable us to test GR at a strong and dynamical gravitational field.

A QNM is expressed as a damped sinusoid

[TABLE]

where $\omega=\omega_{\mathrm{R}}+\mathrm{i}\omega_{\mathrm{I}}$ is the QNM frequency, and $A_{0}$, $t_{0}$, $\varphi_{0}$ are the initial amplitude, time, and phase, respectively. Nakano et al. Nakano et al. (2015) have proposed a method to test GR using observed QNM GWs. While the waveform in the ringdown phase consists of the sum of an infinite number of QNMs, they focused on the dominant, slowest-damped mode [$(l,m,n)=(2,2,0)$] and revealed that the QNM frequency can be estimated with an accuracy at the $5\sigma$ level for a typical case of the mass and spin of the remnant black hole if the signal-to-noise ratio (SNR) of observed GWs is sufficiently high Kinugawa et al. (2016). It means that we can conclude GR is inconsistent with the event if the QNM frequency estimated from it takes a value in the region forbidden by GR. This method can determine how high SNR for the QNM signal is required to test GR for each mass of BHs. We need to estimate the starting time of the QNM in the ringdown phase to precisely estimate the QNM frequency and to calculate the SNR for the QNM from observed GWs from BBHs. According to Nakano et al. Nakano et al. (2015), the roughly estimated rate of detections of BBHs with SNR $>50$ by the second generation GW detectors, such as Advanced LIGO Abbott et al. (2016a); Aasi et al. (2015), Advanced Virgo Abbott et al. (2016a); Acernese et al. (2015), and KAGRA Kuroda et al. (1999); Aso et al. (2013), is 2 yr*-1*. To make the best of such chances to test GR, establishing a method of estimating the starting time is important. However, it is difficult to phenomenologically define the starting time, because the transition from the merger to ringdown phase still has not been well understood. In a theoretical approach, effects of the determination of the starting time for residuals in the fitting of QNM signals and for extractions of parameters have been investigated by using only QNM signals Dorband et al. (2006) and by using full numerical-relativity waveforms of BBH Berti et al. (2007a, b). In this paper, we propose a method to estimate the starting time in a signal-processing approach by using the Hilbert-Huang transform (HHT).

Since the Hilbert-Huang transform is not limited by the time-frequency uncertainty relation, unlike the short-time Fourier transform (STFT) and the wavelet transform (WT), it provides high time-frequency resolution, and therefore enables us to investigate phenomena that have even rapid changes in frequency as well as no or slow changes, while the STFT and the WT can detect only the latter. Application of the HHT to GW data analysis has been examined with various ways Camp et al. (2007); Stroeer et al. (2009); Takahashi et al. (2013); Kaneyama et al. (2016); Sakai et al. (2017).

The HHT consists of two steps: a mode decomposition step and a spectral analysis step. At the mode decomposition step, input data will be decomposed into some intrinsic mode functions (IMFs). Then, at the spectral analysis step, the time evolution of the amplitude and the phase of each IMF are calculated by means of Hilbert transform. This spectral analysis is called Hilbert spectral analysis (HSA). Hence, we can express input data $s(t)$ as

[TABLE]

where $N_{\mathrm{IMF}}$ is the number of IMFs, and $c_{i}(t)$, $a_{i}(t)$, $\phi_{i}(t)$ are the $i$th IMF, the instantaneous amplitude (IA), and the instantaneous phase (IP) respectively, while $r(t)$ is the residual or trend, which is the nonoscillatory part mode of the data. Instantaneous frequency (IF) $f_{i}(t)$ can be defined as

[TABLE]

The original HHT (proposed in Huang et al. (1998)) uses the empirical mode decomposition (EMD) as a decomposition method. It is known, however, that the EMD has many drawbacks, such as mode-mixing Dätig and Schlurmann (2004); Rilling and Flandrin (2008), mode-splitting Yeh et al. (2010), and lack of mathematical foundation Colominas et al. (2015), and therefore many improved methods of the EMD have been being proposed Meignen and Perrier (2007); Hong et al. (2009); Huang and Kunoth (2013); Colominas et al. (2015); Yeh et al. (2010); Colominas et al. (2014). We will adopt the ensemble EMD (EEMD), which is proposed by Z. Wu and N. E. Huang Wu and Huang (2009) to overcome the mode-mixing problem. The detail of the EEMD is also written in Kaneyama et al. (2016).

In this paper, we propose a method of estimating the starting time of the QNM in the ringdown phase of BBH GW, and report the results of evaluation for the efficiency of the method. Following Nakano et al. Nakano et al. (2015), we will focus on the dominant, slowest-damped mode [$(l,m,n)=(2,2,0)$]. It is realized by applying a bandpass filter to each waveform before the HHT and choosing appropriate parameters of the EMD, as described in Sec. II.1 below. The basic concepts and algorithm of the method are described in Sec. II, and the setups and the results of simulation for the evaluation are shown in Sec. III. In Sec. IV, the result of application of the method to the observed data of GW150914 is described. Section V is devoted to a summary.

Throughout this paper, discrete sequences are represented with brackets, such as $x[n]$, and continuous functions are represented with parentheses, such as $s(t)$. The Fourier transform of a continuous function $s(t)$ is represented as $\tilde{s}(f)$.

We perform calculation of the HHT using KAGRA Algorithmic Library (KAGALI) Oohara et al. (2017).

II Strategy of our method

We assume it is already known that the data contains a GW signal from a BBH merger.

If the QNM is perfectly extracted in the $j$th IMF, it follows from Eqs. (1) and (3) that the IA and IP of the IMF are give by

[TABLE]

In reality, the IMF contains also other modes before the QNM starts, and noise components become dominant after the QNM is sufficiently damped. Equations (5) and (6) do not hold in the merger phase and the noise-dominant segment. If the IA fits Eq. (5) most properly during a significant period from a certain time, it can be defined as the starting time of the QNM, and we call the segment a QNM-dominant segment.

The general flowchart of our proposed method is shown in Fig. 1. The method consists of three steps: (A) applying digital filters and the HHT, (B) determining a QNM-dominant segment, (C) and fitting IA and IP of the IMF of interest with the theoretical forms.

II.1 Applying digital filters and the HHT

Raw strain data contain detector noise generated from various origins, such as seismic noise, quantum shot noise, suspension resonance, and so on. Some of them are relatively strong and narrow band noises, called line noises. As shown in Ref. Abbott et al. (2016b), many line noises exist at the Advanced LIGO’s first observing run (O1).

In this step, to extract a sensitive band and to attenuate line noises from raw data, some digital filters are applied to the data. For extraction of the sensitive band, one bandpass filter will be applied, and for attenuation of each line noise, a corresponding notch filter will be applied. After the data processed by the filters, the HHT is applied in order to extract the QNM signal.

We must determine some parameters of the filters and the HHT. The lower cutoff frequency of the bandpass filter should be equal to the lower limit of the properly calibrated band. For other parameters, however, since there is not a universal principle applicable in various situations, they should be determined based on each targeted data. In our method, we attempt to extract the QNM of $(l,m,n)=(2,2,0)$ into the IMF1, which contains the highest-frequency components, in order to prevent the other IMFs from absorbing some components of the target mode. Wu and Huang Wu and Huang (2004) revealed that IMF1 sometimes contains a small but unignorable fraction of components in frequency bands of other IMFs. Therefore, if we attempt to extract a target mode into the second or higher IMF, some fraction of the component may be absorbed into the IMF1. To extract the QNM into the IMF1, we choose the parameters that maximize the signal energy of the IMF1 in the segment containing the merger phase and the beginning of the ringdown phase. The segment will be represented as $[t_{\mathrm{m}},\,t_{\mathrm{m}}+T_{\mathrm{m}}]$, if $t_{\mathrm{m}}$ denotes the time when the amplitude of the GWs is at its maximum and $T_{\mathrm{m}}$ denotes a certain duration that covers the merger phase and the beginning of the ringdown phase. In this paper, we determined to fix $T_{\mathrm{m}}$ at 10 ms for all input data after several checks. We have confirmed that this value of $T_{\mathrm{m}}$ does not cause any trouble for our tests, but if we consider larger total mass than $\sim 60M_{\odot}$, $T_{\mathrm{m}}$ should take a greater value. We choose the upper cutoff frequency of the bandpass filter to make the dominant $(l,m)=(2,2)$ mode be the highest frequency components of the filtered data. The nondominant higher modes, such as $(l,m)=(3,3)$ mode, will be eliminated in this step.

II.2 Determining QNM-dominant segment

After the QNM signal is extracted into the IMF1, we will estimate the starting time of the QNM $\hat{t}_{0}$ and the duration of QNM-dominant segment $\hat{T}$. Figure 3 shows the flowchart of determining the segment.

The IA in the QNM-dominant segment decreases monotonically, but monotonically decreasing segments are not always the QNM-dominant segment. This means that we need to determine the QNM-dominant segment in these segments. To begin with, we search for the longest segment where the IA decreases monotonically. We assume that it extends from $t[n_{\mathrm{b}}]$ to $t[n_{\mathrm{e}}]$, or in the interval $[n_{\mathrm{b}},n_{\mathrm{e}}]$. And then we make a linear regression of the logarithm of IA, $\ln a_{1}[n]$, on a linear function $bt[n]+c$ for every possible subsegment $[n_{0},\,n_{0}+N]$, where $N_{\mathrm{min}}\leq N\leq n_{\mathrm{e}}-n_{\mathrm{b}}$ and $n_{\mathrm{b}}\leq n_{0}\leq n_{\mathrm{e}}-N$, with $N_{\mathrm{min}}$ being a predetermined minimum size of the fitting interval. There is no definite principle of determining $N_{\mathrm{min}}$, but we set it to 5 here. Defining the root mean squared error (RMSE) of the fitting as

[TABLE]

$\hat{n}_{0}(N)$, which denotes the optimal value of $n_{0}$ for each $N$, is determined by

[TABLE]

Finally, $\hat{N}$, which is the optimal value of $N$, is determined as the transition point of a slope of $N$ – $\mathrm{RMSE}(\hat{n}_{0}(N),N)$ plot, based on following two hypotheses:

• •

If the segment $[\hat{n}_{0}(N),\,\hat{n}_{0}(N)+N]$ is a part of the QNM-dominant segment, RMSE gradually increases with $N$, because only detector noise contributes to fitting error.

• •

If the segment $[\hat{n}_{0}(N),\,\hat{n}_{0}(N)+N]$ contains other modes, a merger phase or noise-dominant segment, RMSE rapidly increases with $N$, because fitting basis is not proper for the modes.

Specifically, $\hat{N}$ is determined as

[TABLE]

Figure 4 illustrates how to determine $\hat{N}$. It is plotted for the same data with Fig. 1010(a), which will be described in Sec. III later. In this case, $\hat{N}$ is determined to be $22$.

II.3 Fitting IA and IP in the segment

Since taking the logarithm of Eq. (5) gives

[TABLE]

we can estimate $\omega_{\mathrm{I}}$ and $A_{0}$ if $\ln a_{1}[n]$ is successfully fitted onto a linear regression function $bt[n]+c$ in a segment $[n_{0},\,n_{0}+N]$, such that

[TABLE]

Also, $\omega_{\mathrm{R}}$ and $\varphi_{0}$ can be estimated by fitting $\phi_{1}[n]$ with $dt[n]+e$ as

[TABLE]

Since GWs are redshifted, the QNM frequencies estimated here should be shifted by the factor $1/(1+z)$, where $z$ is the cosmological redshift parameter of the source. For tests of GR, frequencies should be evaluated at the source frame and compared with those estimated with the theory. The values of frequencies given for the rest of the paper are what are evaluated at the source frame,

[TABLE]

Hereinafter, they will be denoted simply by $\omega_{\mathrm{R}}$ and $\omega_{\mathrm{I}}$ without the superscript “src”.

III Simulation

We conduct some simulations in this section to evaluate the efficiency of our method. We make use of waveforms obtained through numerical relativity simulation of BBH mergers by the SXS (Simulating eXtreme Spacetimes) project Chu et al. (2016); SXS . First, we make analysis of pure GWs free from detector noise and compare the obtained values of the QNM frequencies to the theoretical values calculated with the parameters given by the SXS. And then, we make “simulated observation data” by injecting these waveforms to simulated detector noise and apply the procedure to them in order to simulate realistic analysis with observed strain data.

The SXS project provides gravitational waveforms from BBH mergers in their web site SXS . After the detection of GW150914, 96 waveforms newly presented to serve to validate and improve aligned-spin waveform models for GW science Chu et al. (2016). Each waveform is characterized by the mass ratio, $q=m_{1}/m_{2}$, and the initial dimensionless spins, $\vec{\chi}_{1}$ and $\vec{\chi}_{2}$, where $m_{1}$ and $\vec{\chi}_{1}$ are the mass and the initial spin of the primary BH of the BBH, respectively, while $m_{2}$ and $\vec{\chi}_{2}$ are those of the secondary BH. The waveforms can be classified into several categories: Neither is spinning, only one is spinning, and both are spinning. The direction of the spins may be aligned or antialigned with the orbital angular momentum. The magnitude of the spins may be equal or different if both are spinning. We select some waveforms from each category such that the largest mass ratio case, the largest spin magnitude case, and so on are included. As a result, we selected a total of 24 waveforms out of the 96 waveforms. Moreover, we add two waveforms, which best fit GW150914 and GW151226. The distribution of remnant BH’s masses and spins is shown in Fig. 5.

According to Kinugawa et al. Kinugawa et al. (2014), the mass distribution of population III BHs has a peak at $30M_{\odot}$, and for this reason, we set the total mass of a BBH to $60M_{\odot}$. We assume that the BBH is located at $z=0.09$, as GW150914 and GW151226 are, and adopt a flat $\Lambda$CDM cosmology with the cosmological parameters observed by the Planck project Ade et al. (2016), where Hubble parameter $H_{0}=67.9$ $\mathrm{km\;s^{-1}\;Mpc^{-1}}$ and matter density parameter $\varOmega_{\mathrm{m}}=0.306$.

In order to make simulated observation data, we generate Gaussian noises based on Advanced LIGO’s design sensitivity, the zero-detuned high-power sensitivity curve aLI , and add simulated waveforms to them.

To begin with, we apply our method to noise-free gravitational waveforms. Figure 6 illustrates how the method works. The blue solid line shows the waveform analyzed, the red dashed line and the green dot-dashed line are the IA and the IF calculated with the Hilbert spectral analysis, respectively. The QNM-dominant segment we determined is represented by the area filled with the color gray, although only the horizontal time span of the area is meaningful. The upper and lower orange dotted lines show the exponentially decaying curve fitted to the IA and its sign-inversion, respectively. It is apparent that the IF becomes invariant with time in the QNM-dominant segment as expected from Eq. (1), while the IF oscillates in the latter part of this segment. The oscillation is caused by the fact that the amplitude becomes very small, and therefore numerical errors grow significant.

Figures 7 represents the relative errors of the QNM frequencies between what we obtain and the theoretical values, which are calculated for the spin and the mass of the remnant BH by using the table of QNM frequencies against parameters of the BH released by Berti Berti et al. (2006); Ber . The values of the spin and the mass are presented in the metadata of each simulation. It is shown that the obtained frequencies agree with theoretical values sufficiently with some exceptions, while discrepancy in $\omega_{\mathrm{I}}$ is somewhat larger than in $\omega_{\mathrm{R}}$. Some of the SXS waveforms in the ringdown phase look deformed, especially in these exceptional cases. If numerical errors in a waveform are comparable to the small amplitude of the QNM, it will be hard to estimate $t_{0}$ and therefore the discrepancy will be large. Estimation of $\omega_{\mathrm{I}}$ is likely to be affected more sensitively than $\omega_{\mathrm{R}}$ by the estimation error in $t_{0}$. The average and the highest values of the relative errors over the waveforms are shown in Table 1. It means that we can estimate the frequency and the damping rate of the QNM with sufficient accuracy if the QNM signal is exactly extracted from observed data.

To be more realistic, we make 1000 sets of “simulated observation data” for each waveform with simulated noise of Advanced LIGO’s design sensitivity and apply our method to each set of them. Some results of the analysis for three representative waveforms are shown in Fig. 8, and some parameters of these three waveforms are listed in Table. 2, where an optimal SNR, $\rho$, for a waveform $h(t)$ and a detector sensitivity $S_{\mathrm{n}}(f)$ is given by

[TABLE]

We set $f_{\mathrm{min}}$ to $40$ Hz and $f_{\mathrm{max}}$ to $2048$ Hz, following Ref. Allen et al. (2012). Figures 88(a)–88(i) are the histograms of parameters $t_{0}$, $\omega_{\mathrm{R}}$, $\omega_{\mathrm{I}}$, obtained through the analysis, and 88(j)–88(l) show the distributions of the estimated QNM frequencies. The black solid line in each of these figures represents the Schwarzschild limit, and the area above this line is the prohibited area by GR. The bin widths of the histograms are 0.5 s for $t_{0}$, 100 rad/s for $\omega_{\mathrm{R}}$, and 25 rad/s for $\omega_{\mathrm{I}}$. The gray-colored vertical line in each histogram shows the positions along the horizontal axis for the value obtained through the analysis for the corresponding noise-free waveform, which is almost the same value with the theoretical value for $\omega_{\mathrm{R}}$ and $\omega_{\mathrm{I}}$. The centers of the gray circle in each distribution plot indicates the value obtained for the corresponding noise-free waveform.

The source parameters including the estimated starting times of the QNM for noise-free waveforms and the values obtained through our analysis for all 26 waveforms are listed in Table 3.

We confirmed, for noise-free waveforms, that we can reasonably determine the QNM-dominant segment as shown by Fig. 6, and the QNM frequencies obtained in our analysis differ from theoretical values by only 1% or less, although there are larger discrepancies in some cases, as shown by Fig. 7 and Table 1. In Fig. 9, we plot the spin of the remnant BH against the starting time of the QNM. Apparently there is a correlation between the starting time and the remnant spin; correlation coefficient $r$ is 0.764. It results from the fact that a large value of the spin causes the merger phase to be long, and hence the starting time of the ringdown phase is delayed. It is consistent with the analysis by Zhang et al. Zhang et al. (2013).

However, the noise added to the waveform increases the estimation error of $t_{0}$ as shown by Figs. 88(a)–88(c) and Table 3. It tends to result in the smaller value of $t_{0}$, regardless of the SNR of a whole waveform and other parameters. A primary source of this error is imperfectness of the mode decomposition. In Fig. 10, a comparison is drawn between two examples of the decomposition and fitting with the same waveform, 0262, but a different time series of noise. The analysis with the data in Fig. 1010(a) yields the value of $t_{0}$ better than that in Fig. 1010(b), namely, 4.64 ms and 3.17 ms, respectively, while the value evaluated for the noise-free waveform is 4.64 ms. We can observe that the amplitude of IMF2 in Fig. 1010(b) is slightly higher than that in Fig. 1010(a). It means that the signal of the QNM is split into IMF1 and IMF2 in Fig. 1010(b), although that is almost perfectly extracted into IMF1 in Fig. 1010(a). This is caused by a drawback called “mode-splitting” of mode decomposition methods such as the EMD.

IV Applying to GW150914

In this section, we apply our method to the strain data of GW150914 from LIGO Hanford Collaboration (2016) to estimate its ability in analysis of real observed data. The spectral strain sensitivity of the data is shown in Fig. 1111(a). Strong spectral lines are seen in the data. To reduce large oscillation outside the most sensitive frequency band and to attenuate the strong spectral lines inside it, we applied digital filters based on the Butterworth infinite impulse response filter of order 4. The properties of the filters are listed in Table 4, and the spectral strain sensitivity of the filtered data is shown in Fig. 1111(b).

Figure 12 illustrates the results of the analysis, which shows the time series data and the IA of the IMF1 obtained through the HHT analysis, the estimated QNM-dominant segment, and the exponentially decaying curve best-fitted to the IA. It is similar to Fig. 1010(a), suggesting that the level of the accuracy in parameter estimates for real data is same as that for simulated data. The evaluated values of parameters are listed in Table 5. The error ranges of $\omega_{\mathrm{R}}$ and $\omega_{\mathrm{I}}$ shown here are evaluated from only the standard deviation of the fitting of the IF and IA. The corresponding frequency and damping time are $f=\omega_{\mathrm{R}}/2\pi=266\pm 1.6$ Hz and $\tau=1/|\omega_{\mathrm{I}}|=4.73\pm 0.07$ ms, respectively. Assuming that the redshift $z=0.09$, and hence luminosity distance $d_{\mathrm{L}}=410$ Mpc, following the parameter estimation in Ref. Abbott et al. (2016c), and that the QNM we extracted is of $(l,m,n)=(2,2,0)$, the dimensionless spin parameter $\chi_{\mathrm{rem}}$ and the mass $M_{\mathrm{rem}}$ of the remnant BH are evaluated from the QNM frequency $\omega$ by using the conversion table Berti et al. (2006); Ber . They are shown in Table 5, too, as well as measures of statistical errors corresponding to those of $\omega_{\mathrm{R}}$ and $\omega_{\mathrm{I}}$. These values are apparently larger than those given by Ref. Abbott et al. (2016d). The reason is that, as shown in Fig. 8, the systematic error due to noise in the observed data is large, and both $\omega_{\mathrm{R}}$ and $\omega_{\mathrm{I}}$ are tend to be underestimated. So that the larger values of $\chi_{\mathrm{rem}}$ and $M_{\mathrm{rem}}$ are given.

V Summary and Discussion

We proposed a method of estimating the starting time of the QNM in the ringdown phase of GWs from BBH merger, and made its evaluation. The determination of the starting time of the QNM is essential for the test of GR proposed by Nakano et al. Nakano et al. (2015). We consider it as the commencing time of the segment in which the IA calculated with the HHT is well-fitted with an exponentially decaying form as it is expected for the QNM. Also, we can estimate the QNM frequency of a remnant BH with this method. The accuracy of the method was evaluated by applying it to noise-free gravitational waveforms free from detector noise provided by SXS group, simulated observation data, and then the real observed data of GW150914. The simulated observation data are made by injecting waveforms into Gaussian noise based on Advanced LIGO’s design sensitivity.

The results of analyzing noise-free waveforms show that the recovered QNM frequencies are consistent with theoretical values with accuracy of 0.32 % and 0.70 %, on average, for the real and imaginary parts, respectively. This means that we can recover parameters of the remnant BH correctly ONLY if the QNM signal is perfectly extracted. In future, third generation GW detectors, such as the Einstein Telescope Punturo et al. (2010) and B-DECIGO Nakamura et al. (2016) will observe BBH signals with enormously high SNR. With these detectors, our method can work properly and contribute to some tests of GR using QNMs. We find also that starting times of the QNM have a correlation with spins of remnant BHs. It is attributed to the fact that the merger phase remain longer, and hence the starting time of the QNM becomes later as the spin of the remnant BH is larger.

At present, the efficiency of our method for observed data is not so high. It is mainly caused by the fact that the EMD is still in developing stage. Noise in observed data will often give rise to the mode-splitting in the EMD used in the HHT. In this case, the peak value of the IMF is reduced, and therefore the signal will be observed to decay more slowly. It also causes the amplitude of the IMF to fit from earlier to an exponentially decaying curve. It results in the smaller value of the imaginary part of the QNM frequency and the smaller value of the starting time of the QNM. The deformation of the IMF due to noise usually leads to the smaller value of the IF through the Hilbert spectral analysis, and thus the real part of the QNM will be underestimated.

Even though several improvements of the EMD have been proposed, the mode-splitting problem has not yet been resolved. If the mode-splitting is controlled more effectively and thus the signal of the QNM is extracted into a single IMF, the parameters will be estimated more accurately through our method. One promising solution is to take sparsity in the frequency domain into account in the mode decomposition. This idea is based on the fact that the instantaneous frequency of the QNM is constant. It should be noted, however, that this solution can be applied to only some specific cases, in which the frequency of the mode of interest is known to be constant, and that it will not always resolve the mode-splitting problem.

We are now focusing on the only dominant, slowest-damped mode $(l,m,n)=(2,2,0)$, by filtering higher modes such as $(l,m)=(3,3)$ modes. According to London, et al. London et al. (2014), although $(l,m,n)=(2,2,1)$ and $(l,m,n)=(2,1,0)$ modes are in the same band of the target mode, their amplitudes are 1 order of magnitude lower than the target mode. We have not considered the contributions of these modes in this paper. For further work, we are planning to analyze mixtures of various modes. Improvement of mode decomposition will enhance the ability to decompose different modes into different IMFs properly, and therefore we will be able to investigate higher QNMs, as well as the fundamental mode, accurately.

Acknowledgement

This research made use of data obtained from the LIGO Open Science Center LOS , a service of LIGO Laboratory and the LIGO Scientific Collaboration, and made use of data generated by the Simulating eXtreme Spacetimes. LIGO is funded by the U.S. National Science Foundation. This work was in part supported by MEXT Grant-in-Aid for Scientific Research on Innovative Areas “New Developments in Astrophysics Through Multi-Messenger Observations of Gravitational Wave Sources” (Grants No. 24103005 and No. 24103006) and “Gravitational Wave Physics and Astronomy: Genesis” (Grants No. 17H06357 and No. 17H06358). This work was also supported in part by the joint research program of the Institute for Cosmic Ray Research, University of Tokyo (K.S., K.O., M.K. and H.T.), JSPS Grant-in-Aid for Scientific Research (C) (Grant No. 15K05071; K.O., Grant No. 16K05347; H.N., Grant No. 17K05437; H.T.) and JSPS Grant-in-Aid for Young Scientists (B) (Grant No. 26800129; H.T.). We would like to thank Nami Uchikata for carefully reading our manuscript.