Including mode mixing in a higher-multipole model for gravitational waveforms from nonspinning black-hole binaries

TL;DR

This paper enhances gravitational waveform models for nonspinning black-hole binaries by including and accurately modeling mode mixing effects in subdominant modes, achieving high faithfulness for improved detection and parameter estimation.

Contribution

It introduces a method to incorporate mode mixing into analytical waveform models for subdominant modes, improving their accuracy and faithfulness.

Findings

Mode mixing significantly affects subdominant gravitational wave modes.

The new models achieve over 99.6% faithfulness.

Inclusion of mode mixing improves waveform accuracy for large mass ratios.

Abstract

As gravitational-wave (GW) observations of binary black holes are becoming a precision tool for physics and astronomy, several subdominant effects in the GW signals need to be accurately modeled. Previous studies have shown that neglecting subdominant modes in the GW templates causes an unacceptable loss in detection efficiency and large systematic errors in the estimated parameters for binaries with large mass ratios. Our recent work [Mehta et al., Phys. Rev. D 96, 124010 (2017)] constructed a phenomenological gravitational waveform family for nonspinning black-hole binaries that includes subdominant spherical harmonic modes $(\ell = 2, m = \pm 1)$, $(\ell = 3, m = \pm 3)$, and $(\ell = 4, m = \pm 4)$ in addition to the dominant quadrupole mode, $(\ell = 2, m=\pm 2)$. In this article, we construct analytical models for the ($\ell = 3, m = \pm 2$) and ($\ell = 4, m = \pm 3$) modes and include them in the existing waveform family. Accurate modeling of these modes is complicated by the mixing of multiple spheroidal harmonic modes. We develop a method for accurately modeling the effect of mode mixing, thus producing an analytical waveform family that has faithfulness >99.6%.

Tables

Table 1. Table 1: Summary of the parameters of the NR waveforms used in this paper: q := m 1 / m 2 assign 𝑞 subscript 𝑚 1 subscript 𝑚 2 q:=m_{1}/m_{2} is the mass ratio of the binary, M ​ ω orb 𝑀 subscript 𝜔 orb M\omega_{\mathrm{orb}} is the orbital frequency after the junk radiation and e 𝑒 e is the residual eccentricity. The waveforms listed under the title Fitting are used to produce the analytical fits described in Sec.  II.3 while those listed under the title Verification are used for assessing the faithfulness of the analytical model in Sec.  II.5 .

Simulation ID	$q$	$M{\omega }_{\mathrm{orb}}$	$e$	# orbits
Fitting
SXS:BBH:0198	$1.20$	$0.015$	$2.0\times {10}^{-4}$	$20.7$
SXS:BBH:0201	$2.32$	$0.016$	$1.4\times {10}^{-4}$	$20.0$
SXS:BBH:0200	$3.27$	$0.017$	$4.1\times {10}^{-4}$	$20.1$
SXS:BBH:0182	$4.00$	$0.020$	$6.8\times {10}^{-5}$	$15.6$
SXS:BBH:0297	$6.50$	$0.021$	$5.9\times {10}^{-5}$	$19.7$
SXS:BBH:0063	$8.00$	$0.019$	$2.8\times {10}^{-4}$	$25.8$
SXS:BBH:0301	$9.00$	$0.023$	$5.7\times {10}^{-5}$	$18.9$
SXS:BBH:0185	$9.99$	$0.021$	$2.9\times {10}^{-4}$	$24.9$
Verification
SXS:BBH:0066	$1.00$	$0.012$	$6.4\times {10}^{-5}$	$28.1$
SXS:BBH:0184	$2.00$	$0.018$	$7.6\times {10}^{-5}$	$15.6$
SXS:BBH:0183	$3.00$	$0.019$	$6.3\times {10}^{-5}$	$15.6$
SXS:BBH:0182	$4.00$	$0.020$	$6.8\times {10}^{-5}$	$15.6$
SXS:BBH:0187	$5.04$	$0.019$	$5.0\times {10}^{-5}$	$19.2$
SXS:BBH:0181	$6.00$	$0.017$	$7.9\times {10}^{-5}$	$26.5$
SXS:BBH:0298	$7.00$	$0.021$	$4.0\times {10}^{-4}$	$19.7$
SXS:BBH:0063	$8.00$	$0.019$	$2.8\times {10}^{-4}$	$25.8$
SXS:BBH:0301	$9.00$	$0.023$	$5.7\times {10}^{-5}$	$18.9$
SXS:BBH:0185	$9.99$	$0.021$	$2.9\times {10}^{-4}$	$24.9$

Equations

$$\mathbf{h}(t;\iota ,{\phi }_0)=\sum _{\ell \ge 2}\sum _{|m|\le \ell }{Y}_{\ell m}(\iota ,{\phi }_0){}^{\text{Y}}{\mathbf{h}}_{\ell m}(t).$$

$$\mathbf{h}(t;\iota ,{\phi }_0)=\sum _{{\ell }^{\prime }\ge 2}\sum _{|m|\le {\ell }^{\prime }}\sum _{n\ge 0}{S}_{{\ell }^{\prime }mn}(\iota ,{\phi }_0){}^{\text{S}}{\mathbf{h}}_{{\ell }^{\prime }mn}(t).$$

$$S_{{\ell }^{\prime }mn}=\sum _{\ell \ge |m|}{\mu }_{m\ell {\ell }^{\prime }n}^{\ast }{Y}_{\ell m},$$

$$\mathbf{h}(t;\iota ,{\phi }_0)=\sum _{{\ell }^{\prime }\ge 2}\sum _{|m|\le {\ell }^{\prime }}\sum _{n\ge 0}\sum _{\ell \ge |m|}{\mu }_{m\ell {\ell }^{\prime }n}^{\ast }{Y}_{\ell m}(\iota ,\phi ){}^{\text{S}}{\mathbf{h}}_{{\ell }^{\prime }mn}(t).$$

$${}^{\text{Y}}{\mathbf{h}}_{\ell m}(t)=\sum _{{\ell }^{\prime }\ge |m|}\sum _{n\ge 0}{}^{\text{S}}{\mathbf{h}}_{{\ell }^{\prime }mn}(t){\mu }_{m\ell {\ell }^{\prime }n}^{\ast }.$$

$${}^{\text{Y}}{\mathbf{h}}_{\ell m}(t)\simeq \sum _{{\ell }^{\prime }\le \ell }{}^{\text{S}}{\mathbf{h}}_{{\ell }^{\prime }m0}(t){\mu }_{m\ell {\ell }^{\prime }0}^{\ast }.$$

$${}^{\text{S}}{\mathbf{h}}_{320}(t)$$

$${}^{\text{S}}{\mathbf{h}}_{430}(t)$$

$$\begin{array}{rl}{}^{\text{S}}{\tilde{\mathbf{h}}}_{320}^{\mathrm{R}}(f)& \simeq ({\alpha }_1{\mu }_{2330}^{\mathrm{R}}-{\alpha }_2{\mu }_{2330}^{\mathrm{I}})/|{\mu }_{2330}{|}^2,\\ {}^{\text{S}}{\tilde{\mathbf{h}}}_{320}^{\mathrm{I}}(f)& \simeq ({\alpha }_2{\mu }_{2330}^{\mathrm{R}}+{\alpha }_1{\mu }_{2330}^{\mathrm{I}})/|{\mu }_{2330}{|}^2,\end{array}$$

$$\begin{split}\alpha_{1}&:={}^{\text{Y}}\tilde{\mathbf{h}}_{32}^{\mathrm{R}}(f)-\bigg{(}{}^{\text{Y}}\tilde{\mathbf{h}}_{22}^{\mathrm{R}}(f)\rho_{2320}^{\mathrm{R}}+{}^{\text{Y}}\tilde{\mathbf{h}}_{22}^{\mathrm{I}}(f)\rho_{2320}^{\mathrm{I}}\bigg{)},\\ \alpha_{2}&:={}^{\text{Y}}\tilde{\mathbf{h}}_{32}^{\mathrm{I}}(f)+\bigg{(}{}^{\text{Y}}\tilde{\mathbf{h}}_{22}^{\mathrm{R}}(f)\rho_{2320}^{\mathrm{I}}-{}^{\text{Y}}\tilde{\mathbf{h}}_{22}^{\mathrm{I}}(f)\rho_{2320}^{\mathrm{R}}\bigg{)}.\end{split}$$

$$A_{\ell m}^{\text{U}}(f)$$

$${\Psi }_{\ell m}^{\text{U}}(f)$$

$$A_{\ell m}^{\text{U, mod}}(f)=\{\begin{array}{ll}{A}_{\ell m}^{\mathrm{IM}}(f),& f<f_{\ell m}^{\mathrm{A}},\\ A_{\ell m}^{\mathrm{RD}}(f),& f\ge f_{\ell m}^{\mathrm{A}},\end{array}$$

$$A_{\ell m}^{\mathrm{IM}}(f)=A_{\ell m}^{\mathrm{PN}}(f)\left[1+\sum _{k=0}^{k=1}\left({\alpha }_{k,\ell m}+{\alpha }_{k,\ell m}^{\text{L}}\ln v_f\right)v_f^{k+8}\right],$$

$$A_{\ell m}^{\mathrm{RD}}(f)=w_{\ell m}{e}^{-{\lambda }_{\ell m}}\left|\text{mathpzc}{B}_{\ell m}(f)\right|,$$

$$\text{mathpzc}{B}_{\ell m}(f)=\frac{{\sigma }_{\ell m}-\mathrm{i}f}{f_{\ell m}^2+({\sigma }_{\ell m}-\mathrm{i}f{)}^2}.$$

$${\Psi }_{\ell m}^{\text{U, mod}}(f)=\{\begin{array}{ll}{\Psi }_{\ell m}^{\mathrm{IM}}(f),& f<f_{\ell m}^{\mathrm{P}},\\ {\Psi }_{\ell m}^{\mathrm{RD}}(f),& f\ge f_{\ell m}^{\mathrm{P}},\end{array}$$

$${\Psi }_{\ell m}^{\mathrm{IM}}(f)={\Psi }_{\ell m}^{\mathrm{PN}}(f)+\sum _{k=0}^{k=3}({\beta }_{k,\ell m}+{\beta }_{k,\ell m}^{\mathrm{L}}\ln v_f+{\beta }_{k,\ell m}^{L2}{\ln }^2{v}_f)v_f^{k+8},$$

$${\Psi }_{\ell m}^{\mathrm{RD}}(f)=2\pi f{t}_{\ell m}^{\mathrm{P}}+{\varphi }_{\ell m}^{\mathrm{P}}+\arctan \text{mathpzc}{B}_{\ell m}(f),$$

$$\begin{array}{rl}{\alpha }_{i,\ell m}& =a_{i,\ell m}^{\alpha }+b_{i,\ell m}^{\alpha }\eta +c_{i,\ell m}^{\alpha }{\eta }^2,\\ {\alpha }_{i,\ell m}^{\mathrm{L}}& =a_{i,\ell m}^{\alpha ,\mathrm{L}}+b_{i,\ell m}^{\alpha ,\mathrm{L}}\eta +c_{i,\ell m}^{\alpha ,\mathrm{L}}{\eta }^2,\\ {\beta }_{k,\ell m}& =a_{k,\ell m}^{\beta }+b_{k,\ell m}^{\beta }\eta +c_{k,\ell m}^{\beta }{\eta }^2,\\ {\beta }_{k,\ell m}^{\mathrm{L}}& =a_{k,\ell m}^{\beta ,\mathrm{L}}+b_{k,\ell m}^{\beta ,\mathrm{L}}\eta +c_{k,\ell m}^{\beta ,\mathrm{L}}{\eta }^2,\\ {\beta }_{0,\ell m}^{L2}& =a_{0,\ell m}^{\beta ,L2}+b_{0,\ell m}^{\beta ,L2}\eta +c_{0,\ell m}^{\beta ,L2}{\eta }^2,\\ {\lambda }_{\ell m}& =a_{\ell m}^{\lambda }+b_{\ell m}^{\lambda }\eta +c_{\ell m}^{\lambda }{\eta }^2,\\ f_{\ell m}^X& =(a_{\ell m}^X+b_{\ell m}^X\eta +c_{\ell m}^X{\eta }^2)/M,\end{array}$$

$${}^{\text{Y}}{\tilde{\mathbf{h}}}_{\ell m}^{\text{mod}}(f):=A_{\ell m}^{\text{U, mod}}(f)e^{\mathrm{i}{\Psi }_{\ell m}^{\text{U, mod}}(f)},$$

$$\begin{array}{rl}{}^{\text{M}}{\tilde{\mathbf{h}}}_{32}^{\mathrm{R},\text{ mod}}(f)& ={}^{\text{Y}}{\tilde{\mathbf{h}}}_{22}^{\mathrm{R},\text{ mod}}(f){\rho }_{2320}^{\mathrm{R}}-{}^{\text{Y}}{\tilde{\mathbf{h}}}_{22}^{\mathrm{I},\text{ mod}}(f){\rho }_{2320}^{\mathrm{I}}\\ & +{}^{\text{U}}{\tilde{\mathbf{h}}}_{320}^{\mathrm{R},\text{ mod}}(f){\mu }_{2330}^{\mathrm{R}}-{}^{\text{U}}{\tilde{\mathbf{h}}}_{320}^{\mathrm{I},\text{ mod}}(f){\mu }_{2330}^{\mathrm{I}},\end{array}$$

$${}^{\text{U}}{\tilde{\mathbf{h}}}_{320}^{\text{mod}}(f):={\epsilon }_{32}{}^{\text{Y}}{\tilde{\mathbf{h}}}_{\ell m}^{\text{mod}}(f)e^{\mathrm{i}{\delta }_{32}\pi }.$$

$${\epsilon }_{\ell m}$$

$${\delta }_{\ell m}$$

$$A_{\ell m}(f)$$

$${\Psi }_{\ell m}(f)$$

Full text

I Introduction

The detection of compact binary coalescences is now commonplace for the Advanced LIGO Aasi et al. (2015) and Advanced Virgo Acernese et al. (2015) detectors, and they have now produced their first catalogue of such detections Abbott et al. (2018). The accurate extraction of the parameters and hence the science output from these events depends on the accurate modeling of the gravitational waves (GWs) from such sources. Detection of GWs from compact binaries primarily relies on the method of matched filtering, which requires hundreds of thousands of signal templates to be compared against the data (e.g., Dal Canton and Harry (2017); Mukherjee et al. (2018)). Inference of source parameters from observed signals also relies on comparing the data with theoretical waveform templates Abbott et al. (2016). Although numerical relativity (NR) provides the most accurate template waveforms, the large computational cost and sparse parameter space coverage of the NR simulations make the direct implementation of NR waveforms in GW data analysis challenging (see Abbott et al. (2016); Lange et al. (2017) for some recent work in this direction). Over the past decade, there has thus been a considerable effort devoted to developing quick-to-evaluate accurate waveform models for the detection and parameter estimation of GWs from the inspiral, merger, and ringdown of binary black holes, e.g., Buonanno and Damour (1999, 2000); Pan et al. (2011, 2011, 2014); Damour and Nagar (2009); Damour et al. (2009, 2013); Ajith et al. (2008, 2011); Santamaria et al. (2010); Hannam et al. (2014); Husa et al. (2016); Khan et al. (2016); Taracchini et al. (2014); Babak et al. (2017); Bohé et al. (2017); Nagar et al. (2018); Khan et al. (2018); McWilliams (2019).

Most of these (semi) analytical waveforms contain only the dominant multipoles (quadrupole) of the gravitational radiation, though higher multipole models are now starting to be developed Pan et al. (2011); Blackman et al. (2015, 2017); London et al. (2018); Cotesta et al. (2018); Varma et al. (2019); Nagar et al. (2019); Varma et al. (2019). Studies show that neglecting the higher modes can result in a considerable reduction in the sensitivity of searches for high-mass, higher-mass-ratio binary black holes Varma and Ajith (2017); Calderón Bustillo et al. (2016); Varma et al. (2014); Harry et al. (2018). Neglecting these modes can also lead to systematic biases in the parameter estimation of LIGO events from binaries with large mass ratios or high inclination angles, thus biasing our inference of the astrophysical properties of the sources Varma and Ajith (2017); Varma et al. (2014); Abbott et al. (2017). The inclusion of higher multipoles is also expected to provide several other advantages, such as improvements in the precision of parameters extracted from the data Sintes and Vecchio (2000); Van Den Broeck and Sengupta (2007); Arun et al. (2007); Trias and Sintes (2008); Arun et al. (2009); Graff et al. (2015); Lange et al. (2017); O’Shaughnessy et al. (2017); London et al. (2018); Kumar et al. (2019) and in the accuracy of various observational tests of GR Pang et al. (2018), the detection of GW memory Lasky et al. (2016); Talbot et al. (2018), etc.

In this paper, we extend our previous higher multipole waveform family for nonspinning binary black holes Mehta et al. (2017) to include some additional subdominant spherical harmonic modes ($\ell=3,m=2$ and $\ell=4,m=3$) which have a more complicated behavior in their post-merger part of the waveforms due to an effect known as mode mixing.111Since we are considering nonprecessing binaries, we will only refer to the modes with positive $m$ explicitly, as the modes with negative $m$ can be obtained from the positive $m$ modes by symmetry. Additionally, we will refer to the modes using just $\ell m$ for the remainder of the paper. As a consequence, they posses some unusual bumps in the post-merger part which make it difficult to model them accurately. The dominant cause for mode mixing is the mismatch between angular basis functions used in NR (i.e., spin $-2$ weighted spherical harmonics) and in Kerr black hole perturbation theory (i.e., spin $-2$ weighted spheroidal harmonics) Kelly and Baker (2013).222See Buonanno et al. (2007) for an initial study of this effect and London et al. (2014) for a further study of mode mixing in numerical relativity waveforms. As GW observations are entering a regime of precision astronomy, such as precision tests of general relativity (GR) Abbott et al. (2019); Ghosh et al. (2018); Meidam et al. (2018); Dhanpal et al. (2019), modeling of such subtle effects in the waveforms becomes important.

We introduce a method for approximately extracting the unmixed spheroidal harmonic modes from the spherical harmonic modes (the quasi-normal modes separate in the spheroidal basis). We then model these unmixed modes using suitable phenomenological functions motivated by black hole perturbation theory and finally reintroduce the mixing to obtain the model for the spherical harmonic modes. The resulting waveform model is highly faithful (faithfulness $>99.6\%$) and fast to evaluate.

We describe our method for removing the mode mixing approximately and modeling the resulting “unmixed” modes in Sec. II, where we also describe how we add the mode mixing back in to the final model and test its accuracy by computing matches with hybrid waveforms. We summarize and conclude in Sec. III. Additionally, we list the waveforms used for calibration and validation in Appendix A and give some additional plots in Appendix B. Throughout the paper we denote the binary’s total mass by $M$. We also denote the real and imaginary parts of quantities by a superscript R and I, respectively.

II An improved waveform model for subdominant modes

II.1 Mixing of spherical and spheroidal harmonic modes

The two polarizations $h_{+}(t)$ and $h_{\times}(t)$ of GWs can be expressed as a complex waveform $\mathbf{h}(t):=h_{+}(t)-\mathrm{i}\,h_{\times}(t)$. It is convenient to expand this in terms of the spin $-2$ weighted spherical harmonics so that the radiation along any direction $(\iota,\varphi_{0})$ in the source frame can be expressed as

[TABLE]

The spherical harmonic modes ${}^{\text{Y}}\mathbf{h}_{\ell m}(t)=A_{\ell m}(t)\,e^{\mathrm{i}\,\phi_{\ell m}(t)}$ are purely functions of the intrinsic parameters of the system (such as the masses and spins of the binary), while all the angular dependence is captured by the spherical harmonic basis functions $Y_{\ell m}(\iota,\varphi_{0})$. Here, by convention, the polar angle $\iota$ is measured with respect to the orbital angular momentum of the binary. The leading contribution to $\mathbf{h}(t;\iota,\varphi_{0})$ comes from the quadrupolar $22$ modes.

We construct spherical harmonic modes of hybrid waveforms for different modes using the method described in Varma et al. (2014); Mehta et al. (2017). Figure 1 shows the amplitude (solid lines in left panel) and instantaneous frequency (solid lines in right panel) of the second time derivative of different spherical harmonic modes of the hybrid waveforms with mass ratio $q=4$.333We consider the second time derivative of $\mathbf{h}$ (i.e., the Weyl scalar $\psi_{4}$) here instead of $\mathbf{h}$ itself in order to give a cleaner illustration. If we make the same plot using $\mathbf{h}$, we find additional oscillations, even in modes that are not expected to have significant mode mixing. These oscillations appear to be due primarily to additional constant and linear terms in $\mathbf{h}$ that are removed by taking the time derivatives. Taking a single time derivative of $\mathbf{h}$ (i.e., considering the Bondi news) removes most of the oscillations, but taking a second time derivative removes some remaining oscillations. Since we are concerned with removing the mode mixing in the frequency domain, where the time derivatives correspond to a multiplicative factor, there is nothing lost in illustrating the mode mixing removal in the time domain using $\psi_{4}$. We compute the second time derivative by second-order accurate finite differencing. We note that the $22$, $33$, $44$, and $21$ modes, for which an analytical phenomenological model was presented in Mehta et al. (2017), have smoothly varying amplitude and frequency. On the other hand, the $32$ and $43$ modes have some bumps in the post-merger regime $(t>0)$. The unusual behavior of these modes is attributed to what is known as mode-mixing, where multiple spheroidal harmonic modes are getting mixed in one spherical harmonic mode. The prime cause of the mode mixing is the mismatch between the angular basis that is used in NR simulations to extract waveforms (spherical harmonics) and the one that is used to separate the Teukolsky equations in Kerr black hole perturbation theory (spheroidal harmonics) Kelly and Baker (2013).

The mixing of multiple spheroidal harmonic modes creates multiple frequencies in the ringdown waveform that makes it hard to model them using simple analytical functions. Figure 2 shows an example of the Fourier domain amplitude $|{}^{\text{Y}}\tilde{\mathbf{h}}^{\text{R}}_{\ell m}(f)|$ of different hybrid modes — note the non-monotonic behavior seen in the higher frequencies of the $32$ and $43$ modes. These modes were thus left out in the phenomenological model presented in Mehta et al. (2017). In this paper, we present a phenomenological model for the $32$ and $43$ modes. Our approach is to subtract the effect of mode mixing from these modes which allows us to model these “unmixed” modes using methods that were used earlier, and then reintroduce the effects of mode mixing to obtain the final model.

II.2 Removal of mode mixing from the 32 and 43 modes

The binary merger produces a perturbed black hole which settles into a stationary Kerr black hole. Teukolsky’s solution for GWs from a perturbed Kerr black hole has a natural decomposition in spin $-2$ weighted spheroidal harmonics, $S_{\ell mn}\equiv S_{\ell m}(a_{f}\omega_{\ell mn})$ associated with quasinormal mode (QNM) frequencies $\omega_{\ell mn}$, where $M_{f}a_{f}$ is the spin angular momentum of the final black hole (of mass $M_{f}$). See, e.g., Berti et al. (2006) for information about the properties of these functions. Thus, GW polarizations from the ringdown can be written as

[TABLE]

Here the overtone index $n$ measures the magnitude of the imaginary part of the quasinormal mode frequencies $\omega_{\ell mn}$. Note that the spheroidal harmonic basis functions $S_{\ell^{\prime}mn}$ can be expressed in terms of (spin $-2$ weighted) spherical harmonics $Y_{\ell m}$ as

[TABLE]

where $\mu_{m\ell\ell^{\prime}n}$ are mixing coefficients which can be computed simply using the fits provided by Berti and Klein Berti and Klein (2014) (there are more complicated fits given in London and Fauchon-Jones (2018)) and the star denotes the complex conjugate.444We actually substitute $\mu_{m\ell\ell^{\prime}n}\to(-1)^{\ell+\ell^{\prime}}\mu_{m\ell\ell^{\prime}n}$, where the prefactor corrects for the difference in the sign convention for spin-weighted spherical harmonics that we use—the same convention as Ajith et al. (2007), which is also the one used in the SpEC code Boy —and the one used by Berti and Klein. There is an additional factor of $(-1)^{m}$ that we neglect, as it is fixed for each mode we consider (including its mixed modes). By inserting this expansion in Eq. (2), we have

[TABLE]

Comparing this with Eq. (1), we get

[TABLE]

Thus, spherical harmonic modes of the hybrid waveforms can be written in terms of the spheroidal harmonic modes. From inspection of the different spherical harmonic modes of the NR data, we get an understanding of the relative amplitudes of these modes (see, e.g., Fig. 1 in both Pan et al. (2011) and Cotesta et al. (2018)). We thus make the following approximations when removing the mode mixing:

• •

The amplitudes of the higher spheroidal overtones are negligible because their damping times are factors of $\gtrsim 3$ smaller than those of the leading overtone $n=0$. Hence we will only consider mixing from the leading overtone.

• •

For $\ell=m$ spherical modes, the mixing contribution from any mode except the $\ell\ell 0$ spheroidal mode is negligible.

• •

For a general $\ell m$ spherical mode, contribution from spheroidal modes with $\ell^{\prime}>\ell$ is negligible, since the higher mode amplitudes are much smaller than the $\ell m0$ spheroidal mode, and they are also multiplied by the mixing coefficient which is already small.

As a result of these approximations, a particular $\ell m$ spherical mode will have contribution from spheroidal modes $\ell^{\prime}m0$ with $\ell^{\prime}\leq\ell$ (and the obvious restriction of $\ell^{\prime}\geq|m|$). We thus have

[TABLE]

To determine the spheroidal modes ${}^{\text{S}}\mathbf{h}_{\ell^{\prime}m0}(t)$ from the spherical modes ${}^{\text{Y}}\mathbf{h}_{\ell m}(t)$, we observe that it is a perfectly determined system of coupled equations when we consider different $\ell m$ spherical modes. To be specific, we compute the following spheroidal modes:555The same procedure also works for modes with $m\leq\ell-2$ that have three or more spheroidal modes mixed into the spherical mode in our approximation, e.g., the 42 mode studied in Tiwari (2016). However, this mode has a small enough amplitude that we do not include it in the present study.

[TABLE]

These spheroidal harmonic modes for a binary with $q=4$ are shown as dashed lines in Fig. 1 (as discussed there, we plot the second time derivatives to give a cleaner illustration). It can be seen that the amplitude oscillations seen in the spherical modes (solid lines) are largely absent in the spheroidal modes (dashed lines). In addition, the instantaneous frequency (right panel) of the spheroidal modes approaches the corresponding quasi-normal-mode frequency.

We can also convert Eqs. (7) into the frequency domain, so that we can remove the mode mixing from the frequency domain waveforms. Here we want to compute the Fourier transforms of the real and imaginary parts separately, since in this nonprecessing case we can focus on just modeling the real part, and the imaginary part can be obtained from the real part by a phase shift of $\pi/2$. However, we give the expression for the imaginary part as well, for completeness. A straightforward calculation, i.e., taking the real and imaginary parts of Eq. (7a) and expressing them in the frequency domain, gives us the following form for the $32$ mode:

[TABLE]

where

[TABLE]

Here $\rho_{2320}:=\mu_{2320}/\mu_{2220}$ and ${}^{\text{Y}}\tilde{\mathbf{h}}^{\mathrm{R}}_{\ell m}(f)$, ${}^{\text{Y}}\tilde{\mathbf{h}}^{\mathrm{I}}_{\ell m}(f)$ are the Fourier transforms of the real and imaginary parts of ${}^{\text{Y}}\mathbf{h}_{\ell m}(t)$, respectively. The expressions for the $43$ mode are analogous.

The amplitude $|{}^{\text{Y}}\tilde{\mathbf{h}}_{\ell m}^{\mathrm{R}}(f)|$ in the Fourier domain is shown for the $32$ and $43$ modes in Fig. 3 (lighter shades). There are clearly two features in the $32$ mode at close to the QNM frequencies of the $320$ and $220$ modes and similarly for the $43$ mode. Now the “unmixed” modes are constructed as follows:

[TABLE]

where $A^{\text{U}}_{\ell m}(f)$ and $\Psi^{\text{U}}_{\ell m}(f)$ represent the amplitude and phase of the unmixed modes respectively, while $f^{\text{mix}}_{\ell m}$ is a transition frequency. The parameters $w_{\ell m}^{\text{U}}$ and $\phi_{\ell m}^{\text{U}}$ are determined by demanding the continuity of the amplitude and phase at $f^{\text{mix}}_{\ell m}$, respectively.

To determine $f^{\text{mix}}_{\ell m}$, we note that the bump in the amplitude of a certain $\ell m$ spherical harmonic mode due to the mixing of the $(\ell-1)m$ mode always appears at frequencies slightly below the $\ell m$ mode’s dominant QNM frequency $f_{\ell m}^{\text{QNM}}$. When $f^{\text{mix}}_{\ell m}$ is allowed to be a free parameter, it becomes degenerate with the model parameters [Eq. (18)] and thus makes the model fail, i.e., the parameters appearing in Eq. (18) do not have a simple dependence on $\eta$. We find that fixing $f^{\text{mix}}_{\ell m}=0.9f_{\ell m}^{\text{QNM}}$ gives good agreement of the model parameters with quadratic functions of $\eta$. In Fig. 3, we also plot the unmixed modes (dashed lines). The bumps in the amplitudes of the spherical harmonic modes due to mode mixing are significantly suppressed in the unmixed modes.

II.3 Construction of the analytical waveform model

To construct models for the amplitude $A^{\text{U}}_{\ell m}(f)$ and phase $\Psi^{\text{U}}_{\ell m}(f)$, $\ell m\in\{32,43\}$, we follow exactly the same procedure as Sec. @slowromancapii@ B of Mehta et al. (2017). We calibrate the model to the same hybrid waveforms used to calibrate the model in Mehta et al. (2017). These hybrid waveforms are constructed by matching NR waveforms from the SXS Gravitational Waveform Database SXS ; Mroué et al. (2013); Blackman et al. (2015); Chu et al. (2016), listed in Appendix A, with post-Newtonian (PN)/effective-one-body waveforms, using the procedure described in Varma et al. (2014); Mehta et al. (2017). The imaginary part of the unmixed mode (in the time domain) is related to the real part by a phase shift of $\pi/2$, due to the symmetry of nonprecessing binaries. Hence, we only model the Fourier transform of the real part. The amplitude model is thus

[TABLE]

where $f^{\rm A}_{\ell m}$ denotes the transition frequency from the inspiral-merger part of the waveform to the ringdown in the amplitude. The inspiral-merger part is modelled as

[TABLE]

where $v_{f}=(2\pi Mf/m)^{1/3}$ and $A^{\mathrm{PN}}_{\ell m}(f)$ is the Padé resummed version of the Fourier domain 3PN amplitude of the 32 and 43 modes. The Fourier domain amplitude is obtained using the stationary phase approximation as in Van Den Broeck and Sengupta (2007), starting from the time-domain PN results in Blanchet et al. (2008). We use $P^{0}_{4}$ and $P^{0}_{3}$ Padé approximants for the $32$ and $43$ modes, respectively, similar to our treatment of the other modes in Mehta et al. (2017).

The modeling of $A_{\ell m}^{\mathrm{RD}}(f)$ exactly follows Eq. (2.10) of Mehta et al. (2017), i.e.,

[TABLE]

where $\mathpzc{B}_{\ell m}(f)$ is the Fourier transform of a damped sinusoid:

[TABLE]

The frequencies $f_{\ell m}$ and $\sigma_{\ell m}$ are the real and imaginary parts of the $\ell m0$ quasi-normal mode frequency of a Kerr black hole $\Omega_{\ell m0}=2\pi\,(f_{\ell m}+\mathrm{i}\,\sigma_{\ell m})$, determined from the mass and spin of the final black hole using the fits from Ber ; Berti et al. (2006, 2009). The phenomenological parameter $\lambda_{\ell m}$ in Eq. (13) is determined from fits to numerical Fourier transforms of the hybrid waveforms, while $w_{\ell m}$ is a normalization constant to make the amplitudes continuous at the merger-ringdown matching frequency $f^{\mathrm{A}}_{\ell m}$. The mass and spin of the final black hole are computed from the masses of the initial black holes using fitting formulae calibrated to NR simulations, given in Pan et al. (2011).

Similarly, for the phase model we have

[TABLE]

[TABLE]

where $\Psi^{\rm PN}_{\ell m}(f)$ is the PN phasing of the $\ell m$ mode and $f^{\rm P}_{\ell m}$ denotes the transition frequency from the inspiral-merger part of the waveform to the ringdown in the phase. The ringdown part of the phase is modelled as

[TABLE]

where $t^{\mathrm{P}}_{\ell m}$ and $\phi^{\mathrm{P}}_{\ell m}$ are computed by matching two phases ($\Psi^{\rm IM}_{\ell m}$ and $\Psi^{\rm RD}_{\ell m}$) and their first derivative at the matching frequency $f^{\mathrm{P}}_{\ell m}$.

Now, the phenomenological parameters appearing in the analytical models (for the $32$ and $43$ modes) are represented as quadratic functions of the symmetric mass ratio $\eta$:

[TABLE]

where the index $i$ runs from 0 to 1 and $k$ runs from 0 to 3, while $X\in\{A,P\}$. We also refit the phase of the $22$ mode using the smaller number of coefficients given in (16); the fit in Mehta et al. (2017) has the same form, except that the sum extends up to $k=4$ instead of $k=3$. We use this refit since it improves the $22$ mode’s overlap with high mass ratio hybrid waveforms. Figure 4 provides a comparison of the amplitudes and phases of the unmixed modes in the Fourier domain with the analytical fits given by Eqs. (11) and (15). Figure 5 shows the values of the phenomenological parameters estimated from the hybrid waveforms, as well as the fits given in Eq. (18).666We find that quadratic polynomials in $\eta$ provide sufficiently accurate fits in terms of mismatches. Hence we do not consider higher order fits, even though there appears to be some substructure that would require a higher-order polynomial to fit (Fig. 5). It is possible that some of the structure seen in the $32$ and $43$ modes’ coefficients for $\eta$ close to $0.25$ is related to the fact that the mode mixing removal does not work as well for $q<3$, as discussed below.

II.4 Adding the mode mixing contribution into unmixed modes

Having constructed analytical models for the amplitude and phase of unmixed modes, we need to add the mode mixing contribution back into this model in order to get the analytical model for the amplitude and phase of the spherical harmonic modes ${}^{\text{Y}}\tilde{\mathbf{h}}_{\ell m}^{\mathrm{R}}(f)$. This is done as follows: We denote the (Fourier domain) model waveform by

[TABLE]

where $A^{\text{U, mod}}_{\ell m}(f)$ and $\Psi^{\text{U, mod}}_{\ell m}(f)$ are given by Eqs. (11) and (15). We then write

[TABLE]

where

[TABLE]

The expressions for the $43$ mode are analogous. Here we have introduced two free parameters, $\varepsilon_{\ell m}$ and $\delta_{\ell m}$, corresponding to the amplitude ratio and phase difference at $f_{\ell m}^{\text{mix}}$. We fit these parameters by minimizing the mismatch of ${}^{\text{M}}\tilde{\mathbf{h}}_{32}^{\mathrm{R}}(f)$ with the corresponding hybrid mode. They are represented as cubic functions of the symmetric mass ratio; we find similar functional behavior when we compute the amplitude ratio and phase difference between mixed (spherical) and spheroidal modes [Eq. (8)] at $f_{\ell m}^{\text{mix}}$. Specifically,

[TABLE]

where $\ell m\in\{32,43\}$. The fits for parameters $\varepsilon_{\ell m}$ and $\delta_{\ell m}$ are shown in the Fig. 5.

The amplitude and phase of the complete model for ${}^{\text{Y}}\tilde{\mathbf{h}}_{\ell m}^{\mathrm{R}}(f)$ for the mixed modes ($32$ and $43$) are finally constructed as follows:

[TABLE]

The parameters $w_{\ell m}^{\text{M}}$ and $\phi_{\ell m}^{\text{M}}$ ensure the continuity of amplitude and phase at $f_{\ell m}^{\text{mix}}$, respectively. We compare the results of the final model for the spherical harmonics with the hybrids in Fig. 6.

So far, we have used $q=4$ for all our illustrations. We chose this mass ratio to give a clean illustration in a case where the higher modes are relatively prominent and the mode mixing is still fairly large. (The mode mixing decreases as the mass ratio increases for nonspinning binary black holes, as the final spin decreases with increasing mass ratio.) We find that the mode mixing removal is less effective for smaller mass ratios, for reasons that we do not understand. Nevertheless, we still find that our model provides an accurate representation of the waveforms in these cases, as is shown by the match calculations below. We give illustrations of the mode mixing removal and the accuracy of the model for $q=2.32$ in Appendix B.

II.5 Assessing the accuracy of the analytical model

We assess the accuracy of our model by computing mismatches with the same set of $10$ hybrid waveforms used to validate the model in Mehta et al. (2017) (which only share $4$ waveforms—primarily high mass-ratio ones—with the set of $8$ waveforms that are used to construct the model; see Appendix A). The overlaps are computed assuming the design power spectrum of Advanced LIGO (in the “high-power, zero-detuning” configuration adl ),777This noise curve has recently been updated slightly with newer predictions for the thermal noise upd . We use the older version. assuming a low-frequency cutoff of 20 Hz, for a range of total masses.

Figure 7 shows the comparison of our waveform model in the time domain against the hybrid waveforms for two cases, firstly when the model waveforms contain only four modes, i.e., $22$, $33$, $44$, and $21$, and secondly when it also includes the $32$ and $43$ modes in addition to the four modes mentioned before. The hybrid waveforms contain all modes with $\ell\leq 4$, except for the $m=0$ modes, which are small and not well-resolved in the NR simulations. We see that the inclusion of the two additional modes improves the agreement between the hybrid and phenomenological waveforms. Additionally, comparing with Fig. 3 in Mehta et al. (2017), we see the improvement in the face-on case due to the refit of the $22$ mode, as well as the inclusion of the $32$ mode; the $43$ mode does not contribute for a face-on binary.

While Fig. 7 only shows qualitative agreement between the phenomenological and hybrid waveforms, Fig. 8 shows the mismatch plots. The top panel plots show the mismatch (unfaithfulness) between the hybrid waveforms and the case where the model waveforms contain only four modes for various inclination angles. The bottom panel plots show the mismatch after including the $32$ and $43$ modes in the model waveforms for the same inclination angles. We see that for high mass ratio waveforms, the maximum mismatch reduces from $1\%$ to $0.2\%$ for the highest inclination angle $\iota=\pi/2$. However, mismatches are even lower ($\sim 0.05\%$) for other inclination angles. The lower mass ratio cases are almost unaffected, though they show a little improvement. This is expected, because the contribution of higher modes is significant for high mass ratio and inclination angles.

We also show the improvement in the accuracy of the model for the 22 mode alone for higher mass ratios in Fig. 9. This comes from a refit of this mode’s phase. Unfortunately, this improvement for higher mass ratios comes at the cost of a somewhat larger mismatch for mass ratios of $1$ and $2$. Future work will consider improvements to the structure of the model to improve the mismatch for small mass ratios.

III Summary and Conclusions

In this paper, we extend our analytical frequency-domain phenomenological higher-mode model for gravitational waveforms from nonspinning binary black holes Mehta et al. (2017) to include two additional subdominant modes, namely the $32$ and $43$ modes, in addition to the $22$, $33$, $44$, and $21$ modes. This waveform family now has a faithfulness of $>99.6\%$ for binaries up to a mass ratio of $10$ and a total mass of $200M_{\odot}$, using the design sensitivity Advanced LIGO noise curve and a low-frequency cutoff of $20$ Hz. The two additional modes that we model in this paper, i.e., $32$ and $43$, exhibit the effects of mode-mixing, i.e., having multiple spheroidal harmonic ringdown modes mixed into a single spherical harmonic mode. This leads to bumps in the ringdown part of the waveform. We have introduced a simple way of approximately extracting the unmixed (spheroidal harmonic) modes. We then model these unmixed modes using the method used for the other modes in Mehta et al. (2017). We then reinstate the mode mixing using the models for the unmixed modes to obtain the final model for the spherical harmonic modes. We also refit our model for the dominant $22$ mode to improve its accuracy for large mass ratios. This is the first analytical inspiral-merger-ringdown waveform family that models the mode mixing effect.

We note that efforts are underway to construct effective-one-body/phenomenological models for binary black holes with non-precessing spins London et al. (2018); Cotesta et al. (2018). The simple prescription that we use to model mode mixing effects may be used in these models as well. Future models aiming for a highly accurate description of non-quadrupole modes may also need to consider other sources of mode mixing, e.g., the mode-mixing due to boosts and displacements from the origin, discussed in Boyle (2016). Additionally, as discussed in Appendix B, more accurate determination of the waveforms at infinity in numerical simulations, e.g., through Cauchy-characteristic extraction Reisswig et al. (2009); Taylor et al. (2013); Chu et al. (2016); Handmer et al. (2016), will likely be necessary input for precise models. As GW observations are becoming precision probes of physics and astrophysics, accuracy requirements on GW templates can only grow.

Acknowledgements.

We are grateful to the SXS collaboration for making a public catalog of numerical-relativity waveforms. We also thank Emanuele Berti and Michael Boyle for useful discussions and clarifications, Frank Ohme for helpful comments, and Mark Scheel for providing a Cauchy-characteristic extraction SXS waveform for comparison. A. K. M., P. A., and V. V. acknowledge support from the Indo-US Centre for the Exploration of Extreme Gravity funded by the Indo-US Science and Technology Forum (IUSSTF/JC-029/2016). N. K. J.-M. acknowledges support from the AIRBUS Group Corporate Foundation through a chair in “Mathematics of Complex Systems” at the International Centre for Theoretical Sciences (ICTS) and from STFC Consolidator Grant No. ST/L000636/1. Also, this work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 690904. P. A.’s research was supported by the Max Planck Society through a Max Planck Partner Group at ICTS-TIFR, and by the Canadian Institute for Advanced Research through the CIFAR Azrieli Global Scholars program. V. V.’s research was supported by the Sherman Fairchild Foundation, and NSF grants PHY–170212 and PHY–1708213 at Caltech. Computations were performed at the ICTS cluster Alice. This document has LIGO preprint number LIGO-P1800203-v6.

Appendix A Numerical relativity waveforms

In Table 1, we list the SXS waveforms that are used to fit the coefficients appearing in Sec. II.3 and to assess the accuracy of our model as described in Sec. II.5.

Appendix B Mode mixing removal for lower mass ratios

We give analogs of Figs. 1 and 6 for a mass ratio of $q=2.32$ in Figs. 10 and 11. The first of these figures illustrates that the mode mixing removal is still effective in improving the agreement of the instantaneous frequency with the expected QNM frequency, and in reducing the amplitude oscillations of the 32 mode. However, the mode mixing removal is less effective for $q\lesssim 3$, for reasons we do not fully understand. Nevertheless, the second figure shows that the final model for the mixed modes still agrees well with the Fourier transform of the hybrid.

The numerical noise we find in the instantaneous frequency of the modes is reduced when considering the NR waveform with no extrapolation to infinity. Experimentation with the equal-mass nonspinning Cauchy-characteristic extraction SXS waveform from Taylor et al. (2013) finds that this does not suffer from the numerical noise that is present in the instantaneous frequencies of the analogous finite radius or extrapolated equal-mass nonspinning SXS waveform modes.