# Time domain phenomenological model of gravitational wave subdominant
harmonics for quasi-circular non-precessing binary black hole coalescences

Héctor Estellés  Departament de Física, Universitat de les Illes Balears, IAC3
– IEEC, Crta. Valldemossa km 7.5, E-07122 Palma, Spain  Sascha Husa
Departament de Física, Universitat de les Illes Balears, IAC3 – IEEC, Crta.
Valldemossa km 7.5, E-07122 Palma, Spain  Marta Colleoni  Departament de
Física, Universitat de les Illes Balears, IAC3 – IEEC, Crta. Valldemossa km
7.5, E-07122 Palma, Spain  David Keitel  Departament de Física, Universitat de
les Illes Balears, IAC3 – IEEC, Crta. Valldemossa km 7.5, E-07122 Palma, Spain
Maite Mateu-Lucena  Departament de Física, Universitat de les Illes Balears,
IAC3 – IEEC, Crta. Valldemossa km 7.5, E-07122 Palma, Spain  Cecilio García-
Quirós  Departament de Física, Universitat de les Illes Balears, IAC3 – IEEC,
Crta. Valldemossa km 7.5, E-07122 Palma, Spain  Antoni Ramos-Buades  Max
Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg 1, D-14476 Potsdam, Germany  Departament de Física, Universitat de
les Illes Balears, IAC3 – IEEC, Crta. Valldemossa km 7.5, E-07122 Palma, Spain
Angela Borchers  Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-
Institut, Callinstr. 38, D-30167 Hannover, Germany  Departament de Física,
Universitat de les Illes Balears, IAC3 – IEEC, Crta. Valldemossa km 7.5,
E-07122 Palma, Spain

###### Abstract

In this work we present an extension of the time domain phenomenological model
IMRPhenomT for gravitational wave signals from binary black hole coalescences
to include subdominant harmonics, specifically the $(l=2,m=\pm 1)$,
$(l=3,m=\pm 3)$, $(l=4,m=\pm 4)$ and $(l=5,m=\pm 5)$ spherical harmonics. We
also improve our model for the dominant $(l=2,m=\pm 2)$ mode and discuss mode
mixing for the $(l=3,m=\pm 2)$ mode. The model is calibrated to numerical
relativity solutions of the full Einstein equations up to mass ratio 18, and
to numerical solutions of the Teukolsky equations for higher mass ratios. This
work complements the latest generation of traditional frequency domain
phenomenological models (IMRPhenomX), and provides new avenues to develop
computationally efficient models for gravitational wave signals from generic
compact binaries.

###### pacs:

Valid PACS appear here

##  I Introduction

In order to exploit the information encoded in gravitational wave (GW) signals
from compact binary systems, reliable and computationally efficient models of
such signals are required to perform Bayesian inference Veitch  _et al._
(2015); Ashton  _et al._ (2019) of the source properties. The development of
such waveform models has so far progressed in remarkable synchronicity with
the challenges posed by the improvement in sensitivity of the advanced LIGO
Aasi  _et al._ (2015) and Virgo Acernese  _et al._ (2014) detectors, and the
discovery of new types of compact binaries. A decade after the numerical
relativity (NR) breakthroughs to evolve black holes Pretorius (2005);
Campanelli  _et al._ (2006); Baker  _et al._ (2006), the third generation of
inspiral-merger-ringdown waveform models Husa  _et al._ (2016); Khan  _et al._
(2016); Hannam  _et al._ (2014); Bohé  _et al._ (2016); Pürrer (2016) could be
employed for parameter estimation of the first detected signal Abbott  _et
al._ (2016a), and an extensive study of systematic errors Abbott  _et al._
(2017) found the models to be sufficiently accurate for this relatively loud
signal, although improvements would soon be required for future events, in
particular to include sub-dominant harmonics. Models that include subdominant
harmonics were still not available to analyze the catalog of events from the
first two observation runs Abbott  _et al._ (2019), but within a year an
extensive study of the particularly massive GW170729 event Chatziioannou  _et
al._ (2019) with the first multi-mode models reported improved parameter
estimation results. Several models that include both precession and
subdominant harmonics were completed in time to study events in the third
observation run: clear evidence for such modes has been found for two events,
GW190412Abbott  _et al._ (2020a) and GW190814 Abbott  _et al._ (2020b), and
such models have been used to analyse events for the second catalog of GW
events LIGO Scientific Collaboration (2020).

Further improvements in waveform models will be required to keep up with the
increased sensitivity foreseen for the fourth and future observation runs of
the advanced detector network Abbott  _et al._ (2018), and to meet the goal of
being able to model generic signals. In a previous work Estellés  _et al._
(2020), we presented IMRPhenomTP, a phenomenological model in the time domain
for the dominant $l=2,m=\pm 2$ modes in the co-precessing frame, extended to
cover the $l=2$ subset of modes in the inertial frame of precessing systems
through the twisting-up procedure. This model forms the cornerstone of a new
family of time domain models, which are comparable in construction, accuracy
and computational efficiency with the IMRPhenomX family of frequency domain
models Pratten  _et al._ (2020a); García-Quirós  _et al._ (2020); Pratten  _et
al._ (2020b). The main motivation for such time domain models is twofold:
First, the time domain simplifies the separation of the inspiral, merger, and
ringdown regimes, and can thus simplify tests of general relativity, such as
inspiral-merger-ringdown consistency tests Ghosh  _et al._ (2016, 2017) or
parameterized tests Yunes and Pretorius (2009). Second, some effects may be
easier to model in the time domain, most notably the twisting-up approximation
that is often used to model precession Hannam  _et al._ (2014); Bohé  _et al._
(2016); Ramos-Buades  _et al._ (2020); Pratten  _et al._ (2020b) and in the
time domain does not require the commonly employed stationary phase
approximation Hannam  _et al._ (2014); Bohé  _et al._ (2016); Pratten  _et
al._ (2020b), or alternative techniques which still need to be developed
toward maturity Marsat and Baker (2018).

Here we present an extension of the core non-precessing model Estellés  _et
al._ (2020) that includes a subset of subdominant modes, in particular
$l=2,m=\pm 1$, $l=3,m=\pm 3$, $l=4,m=\pm 4$ and $l=5,m=\pm 5$. While this work
will not present results on precession, the procedures followed in Estellés
_et al._ (2020) carry over to the multimode extension discussed here, and in a
future publication the benefits of working in the time domain for
phenomenologically describing the merger-ringdown regions of precessing
waveforms will be analysed.

The paper is organized as follows. In Sec. II we discuss the model ansatz we
use, which is then calibrated across the parameter space to a set of numerical
waveforms in Sec. III. In Sec. IV we test the quality of the model across the
parameter space, and in Sec. V we test the performance of the model in
Bayesian parameter estimation examples. Finally we conclude in Sec. VI.

##  II Phenomenological model construction

Gravitational radiation is encoded in two independent polarisations $h_{+}$
and $h_{\times}$ that can be decomposed in a suitable basis of spin-weighted
spherical harmonics (SWSH) Goldberg  _et al._ (1967) (defined as in Wiaux  _et
al._ (2007)):

$\begin{split}h(t;\bm{\lambda};r,\theta,\phi)&=h_{+}-ih_{\times}\\\
&=\sum_{l}\sum_{-l\leq m\leq l}h_{lm}(t;\bm{\lambda};r)\
^{-2}Y_{lm}(\theta,\phi),\end{split}$ (1)

where the modes $h_{lm}(t;\bm{\lambda};r)$ express the dependence on time and
the intrinsic physical properties of the source, $\bm{\lambda}$, while the
relative orientation of the source with respect to the observer is encoded in
the SWSH basis functions ${}^{-2}Y_{lm}(\theta,\phi)$. This work will be
restricted to the modelling of binary black hole systems whose individual
spins are aligned or anti-aligned with the orbital angular momentum direction,
which is preserved in this situation. Then, the intrinsic degrees of freedom
are $\bm{\lambda}=\\{m_{1},m_{2},\chi_{1},\chi_{2}\\}$ where $m_{1,2}$ are the
component masses and $\chi_{1,2}=(\bm{S}_{1,2}/m_{1,2}^{2})\cdot\hat{\bm{L}}$
are the dimensionless spin components in the orbital angular momentum
direction $\hat{\bm{L}}$. For modelling purposes, the total mass of the system
$M=m_{1}+m_{2}$ is simply a scale parameter for vacuum spacetimes. We thus
factor it, out working in units where time is expressed in mass units, and
express the dependence on the masses through the symmetric mass ratio
$\eta=m_{1}m_{2}/M^{2}$.

Each mode is a complex function that can be decomposed into a real amplitude
and phase as

$h_{lm}(t)=H_{lm}(t)\exp\\{-i\phi_{lm}(t)\\},$ (2)

with $\dot{\phi}_{lm}(t)\equiv\omega_{lm}(t)>0$ for $m>0$ and
$\omega_{lm}(t)<0$ for $m<0$. The model provides a phenomenological
construction for the $m>0$ modes specified at the beginning, and the $m<0$
modes are computed using the symmetry relation:

$h_{lm}(t)=(-1)^{l}h_{l,-m}^{*},$ (3)

valid for non-precessing systems, where ∗ indicates complex conjugation. By
convention we shift the time coordinate such that $t=0$ occurs for the peak
amplitude of the $(2,2)$ mode.

###  II.1 Model overview

In a previous work Estellés  _et al._ (2020) a phenomenological model in the
time domain for the dominant $l=2$, $m=\pm 2$ modes was presented, IMRPhenomT.
In this work we provide the extension of the model to the modes $l=2,m=\pm 1$,
$l=3,m=\pm 3$, $l=4,m=\pm 4$ and $l=5,m=\pm 5$. The core structure of
IMRPhenomT are $C^{1}$ analytical expressions for the amplitude
$H_{22}(t)=|h_{22}(t)|$, the gravitational phase $\phi_{22}(t)$ and the wave
frequency $\omega_{22}(t)=\dot{\phi}_{22}(t)$ of the $l=2$, $m=2$ mode, in
terms of piece-wise functions defined on a physically motivated partition of
the evolution time of the system: inspiral, merger and ringdown regions.
Following the same strategy, the model presented here provides $C^{1}$
analytical expressions for the amplitudes $H_{lm}(t)=|h_{lm}(t)|$, the GW
phases $\phi_{lm}(t)$ and the wave frequencies
$\omega_{lm}(t)=\dot{\phi}_{lm}(t)$ for the included modes as piece-wise
functions:

$\phi_{lm}(t)=\begin{cases}\phi^{\text{insp}}_{lm}(t)&t\leq
t^{\omega}_{\text{cut},lm}\\\ \\\
\phi^{\text{merger}}_{lm}(t)&t^{\omega}_{\text{cut},lm}\leq t\leq
t^{\text{peak}}_{22}\\\ \\\ \phi^{\text{RD}}_{lm}(t)&t\geq
t^{\text{peak}}_{22},\end{cases}$ (4)

$\omega_{lm}(t)=\begin{cases}\omega^{\text{insp}}_{lm}(t)&t\leq
t^{\omega}_{\text{cut},lm}\\\ \\\
\omega^{\text{merger}}_{lm}(t)&t^{\omega}_{\text{cut},lm}\leq t\leq
t^{\text{peak}}_{22}\\\ \\\ \omega^{\text{RD}}_{lm}(t)&t\geq
t^{\text{peak}}_{22},\end{cases}$ (5)

$H_{lm}(t)=\begin{cases}H^{\text{insp}}_{lm}(t)&t\leq t^{H}_{\text{cut}}\\\
\\\ H^{\text{merger}}_{lm}(t)&t^{H}_{\text{cut}}\leq t\leq
t^{\text{peak}}_{lm}\\\ \\\ H^{\text{RD}}_{lm}(t)&t\geq
t^{\text{peak}}_{lm},\end{cases}$ (6)

where $t_{\text{cut}}$, the separation between the inspiral and merger
regimes, was selected in the previous work to correspond to the time of the
minimum energy circular orbit (MECO) Cabero  _et al._ (2017) for the phase and
frequency, and to half of this time for the amplitude. However, already in the
previous work it was realised that, while this time sets a natural ending for
the validity of the post-Newtonian (PN) constructions employed in the inspiral
regime, it complicates the extension of the model to high mass ratios since
the merger ansatz did not describe correctly the too elongated merger region
for these cases, and similarly for high spins where the merger region was too
short to be correctly reconstructed. For these reasons, we have now chosen a
different implementation of the separation times between regions, which will
be discussed in this section. We present details of the modelling of each of
these regions, highlighting the main differences with the version of
IMRPhenomT described in Estellés  _et al._ (2020).

###  II.2 Inspiral

For modelling the inspiral regime, we follow the approach of IMRPhenomXHM
García-Quirós  _et al._ (2020): We model the amplitude and the phase of the
$\ell=|m|=2$ mode by adding higher order pseudo-PN terms to known PN results,
and we describe the phase of the sub-dominant modes by augmenting the phase of
the dominant mode with PN terms.

####  II.2.1 Orbital frequency

In the inspiral regime, where the assumptions of the PN framework and the
energy balance condition are satisfied, PN analytical solutions for the
dynamical evolution of the orbital frequency
$\omega_{\mathrm{orb}}(t)\equiv\dot{\phi}_{\mathrm{orb}}(t)$ can be employed.
In the adiabatic approximation, we can consider the balance equations

$\dfrac{d\phi_{\mathrm{orb}}}{dt}-\dfrac{v^{3}}{M}=0,$ (7a)
$\dfrac{dv}{dt}+\dfrac{\mathcal{F}(v)}{ME^{\prime}(v)}=0,$ (7b)

where $v=\omega_{\mathrm{orb}}^{1/3}$ is the PN velocity,
$\phi_{\mathrm{orb}}(t)$ is the orbital phase of the binary, $\mathcal{F}(v)$
is the gravitational wave luminosity and $E(v)$ is the binding energy of the
system. Currently, the gravitational wave luminosity and the binding energy of
the system are known up to 3.5PN order for quasi-circular orbits Blanchet  _et
al._ (1995); Blanchet (1998); Blanchet  _et al._ (2004). In this work we use
the expressions for the energy fluxes and binding energy at 3.5PN order,
although we note that recently the 4.5PN coefficients for the energy flux have
also been calculated Marchand  _et al._ (2016). Expanding the quotient on
equation (7b), different families of analytical templates Buonanno  _et al._
(2009) can be found for the velocity derivative, including some analytical
approximations to the orbital phase and frequency, as in the TaylorT3
approximant:

$\omega_{\mathrm{orb}}(t)\approx\omega^{(T3)}_{n/2}(t)=\omega_{N}(t)\sum_{k=0}^{n}\hat{\omega}_{k}\theta^{k}(t),$
(8)

where $\theta(t)=[\eta(t_{0}-t)/(5M)]^{-1/8}$, $\omega_{N}=\theta^{3}/8$ and
the PN coefficients $\hat{\omega}_{k}$ are given in Appendix A1 of Estellés
_et al._ (2020).

Nevertheless, in order to improve the accuracy of the description up to the
end of the inspiral region, we need to expand the expression to higher order
adding 6 extra pseudo-PN orders with unknown coefficients:

$\omega^{\text{insp}}_{22}(t)=\omega^{(T3)}_{3.5}(t)+\omega_{N}(t)\sum_{k=8}^{13}\hat{c}_{k}\theta^{k}.$
(9)

In order to obtain the value of the coefficients across parameter space,
calibration with numerical simulations is required. However, instead of
performing parameter space fits of the unknown coefficients, following Husa
_et al._ (2016); Pratten  _et al._ (2020a) it is more convenient to fit the
value of the frequency at particular times, placing collocation points, and
then solve the coefficients through a linear system of equations:

$\omega_{N,i}\sum_{k=8}^{13}\hat{c}_{k}\theta_{i}^{k}=\lambda^{\omega}_{i}-\omega^{(T3)}_{3.5}(\theta_{i}).$
(10)

Expression (10) corresponds to a linear system of six equations that can be
solved to obtain the six unknown coefficients $\hat{c}_{k}$ in terms of known
frequency values $\lambda^{\omega}_{i}$ at particular $\theta_{i}$ values. The
frequency values $\lambda^{\omega}_{i}$ are fitted across parameter space from
numerical simulations.

Both the number of coefficients and the placement of the collocation points
differ from the original implementation presented in Estellés  _et al._
(2020). Instead of placing the collocation points in the late inspiral at
fixed times, we found that a better behaviour is achieved across parameter
space if they are placed at fixed $\theta(t)$ values instead. We have selected
the following placement:

$\theta_{i}=\\{0.33,\ 0.45,\ 0.55,\ 0.65,\ 0.75,\ 0.82\\}$ (11)

and their relation with time is illustrated in Fig. 1. With this placement,
the Newtonian term of TaylorT3, $\omega_{N}=\theta^{3}/8$, is the same for
each collocation point across parameter space. Furthermore, the number of
extra coefficients added to the TaylorT3 approximant has also increased with
respect to the previous implementation. This was decided in order to increase
the accuracy towards more unequal mass systems, especially when the primary
spin magnitude is high. As an illustrative example, in Fig. 2 we show the
$l=2$, $m=2$ frequency evolution of two configurations, both at mass ratio 8,
one for high positive and one for high negative primary spin magnitude,
computed from the SEOBNRv4 approximant Bohé  _et al._ (2017) and from the
TaylorT3 approximant. As can be seen in the figure, the four later collocation
points are needed for correcting the deviation of TaylorT3 with respect to the
correct solution represented by SEOBNRv4. The first two early collocation
points are needed for regularising the behaviour at low frequency after the
addition of the other extra terms, whose solution can sometimes lead to
deviations in the low frequency regime.

The data to be taken for each collocation point in order to produce the needed
parameter space fits is selected according to the actual time at which the
collocation point is placed in each case and the region of the parameter space
where the case is located. Inside the parameter space region covered by the NR
simulations, if the collocation point is placed inside the time interval of
the corresponding NR simulation, the value from the NR simulation is taken. If
the collocation point is placed at a time earlier than the minimum available
time for the corresponding NR simulation, then the result from a SEOBNRv4
waveform is taken. For mass ratios above the NR parameter space coverage, we
employ the data set of Teukolsky-EOB hybrids we have previously used in
García-Quirós  _et al._ (2020), in particular at mass ratios 80, 200 and 1000,
for different spin configurations.

Figure 1: Placement of inspiral orbital frequency collocation points. Dashed
vertical lines show the position of the different collocation points at fixed
$\theta$ values, and the solid coloured lines show the relation between time
and $\theta$ for different mass ratio configurations. Figure 2: Frequency
evolution of $l=2$, $m=2$ mode as a function of the $\theta$ parameter.
Purple: $q=8,\ \chi_{1}=0.8$ configuration. Orange: $q=8,\ \chi_{1}=-0.8$
configuration. Solid: from SEOBRNRv4 waveform. Dashed: 3.5PN TaylorT3
implementation. Dashed black vertical lines: collocation point placement.

####  II.2.2 Mode amplitude

For describing the inspiral amplitude of each mode, we employ the 3PN order
expressions from Blanchet  _et al._ (2008) with 2PN spin corrections from
Buonanno  _et al._ (2013) and 1.5PN contributions from Arun  _et al._ (2008).
For the amplitude of the $(2,2)$, we also include 3.5PN corrections from Faye
_et al._ (2012). The general form of the amplitude can be written as:

$H_{lm}^{PN}(t)=2\eta\sqrt{\dfrac{16}{5}}x(t)\sum_{k=0}^{7}\hat{h}^{lm}_{k}x(t)^{k/2}$
(12)

where
$x(t)=(\omega_{\mathrm{orb}}(t))^{2/3}=(\omega^{\mathrm{insp}}_{22}(t)/2)^{2/3}$
and the PN coefficients $\hat{h}^{lm}_{k}$ are in general complex functions of
$\\{\eta,\chi_{1},\chi_{2}\\}$. The original expression for the amplitude of
Eq. (12) for each mode can be found in Appendix A. However, it is more
convenient to rotate the complex amplitudes such that the dominant
contribution is given by their real part, since in this way the absolute value
of the amplitude can be approximated by the real part, while the smaller
imaginary part will contribute to the phase evolution as will be explained
below. The corresponding rotations applied to the expressions of Appendix A
are:

$\displaystyle H_{lm}^{{}^{\prime}PN}(t)$
$\displaystyle=\exp(i\Delta\phi_{lm})H_{lm}^{PN}(t),$ (13a)
$\displaystyle\Delta\phi_{22}$ $\displaystyle=0,$ (13b)
$\displaystyle\Delta\phi_{21}$ $\displaystyle=\pi/2,$ (13c)
$\displaystyle\Delta\phi_{33}$ $\displaystyle=-\pi/2,$ (13d)
$\displaystyle\Delta\phi_{44}$ $\displaystyle=\pi,$ (13e)
$\displaystyle\Delta\phi_{55}$ $\displaystyle=\pi/2.$ (13f)

Similar to what is done for the frequency, in order to improve the accuracy of
the description up to the end of the inspiral, we add three extra pseudo-PN
orders with unknown coefficients that will be obtained from calibration to
simulations:

$H_{lm}^{\mathrm{insp}}(t)={H_{lm}^{{}^{\prime}PN}}(t)+2\eta\sqrt{\dfrac{16}{5}}x(t)\sum_{k=8}^{10}\hat{d}^{lm}_{k}x(t)^{k/2},$
(14)

where now the extra coefficients $\hat{d}^{lm}_{k}$ are real functions of the
parameter space $\\{\eta,\chi_{1},\chi_{2}\\}$. As for the frequency
calibration, it is more convenient to perform parameter space fits of a set of
collocation points from the amplitude, and employ them to solve the value of
the coefficients for each case employing a linear system of equations:

$\displaystyle
2\eta\sqrt{\dfrac{16}{5}}\begin{bmatrix}x^{5}(t^{H}_{1})&x^{11/2}(t^{H}_{1})&x^{6}(t^{H}_{1})\\\
x^{5}(t^{H}_{2})&x^{11/2}(t^{H}_{2})&x^{6}(t^{H}_{2})\\\
x^{5}(t^{H}_{3})&x^{11/2}(t^{H}_{3})&x^{6}(t^{H}_{3})\end{bmatrix}\begin{bmatrix}\hat{d}^{lm}_{8}\\\
\hat{d}^{lm}_{9}\\\ \hat{d}^{lm}_{10}\end{bmatrix}$ (15)
$\displaystyle=\begin{bmatrix}\lambda^{H}_{1,lm}-\Re{H_{lm}^{{}^{\prime}PN}}(t^{H}_{1})\\\
\lambda^{H}_{2,lm}-\Re{H_{lm}^{{}^{\prime}PN}}(t^{H}_{2})\\\
\lambda^{H}_{3,lm}-\Re{H_{lm}^{{}^{\prime}PN}}(t^{H}_{3})\end{bmatrix},$

where the amplitude collocation points $(t^{H}_{i},\ \lambda^{H}_{i,lm})$ are
placed at fixed time positions in the late inspiral:

$t^{H}_{i}=\\{-2000M,-250M,-150M\\}.$ (16)

####  II.2.3 Phases of subdominant modes

For modelling the subdominant harmonic phases during the inspiral, a good
approximation is achieved by just rescaling the frequency (and the phase) of
the $(2,2)$ mode Baker  _et al._ (2008); Kelly  _et al._ (2011):

$\omega_{lm}^{\mathrm{insp}}(t)\approx\frac{m}{2}\omega_{22}^{\mathrm{insp}}(t),$
(17a)
$\phi_{lm}^{\mathrm{insp}}(t)\approx\frac{m}{2}\phi_{22}^{\mathrm{insp}}(t).$
(17b)

However, this approximation loses accuracy as the merger is approached,
accumulating some dephasing during the late inspiral. In the modelling of the
frequency domain model IMRPhenomXHM, it was realized that a good approximation
to correct the deviation of Eq. (17) is to employ the phase contribution from
the complex mode amplitude:

$\Lambda_{lm}(t)\equiv\arg{H_{lm}^{\mathrm{insp}}(t)}=\arctan{\Big{(}\dfrac{\Im{H_{lm}^{\mathrm{insp}}(t)}}{\Re{H_{lm}^{\mathrm{insp}}(t)}}\Big{)}}.$
(18)

For tracking also the correct phase relation between the modes satisfying the
specific tetrad convention, we need to rotate back by the global rotation of
Eq. (13) applied to the amplitude. With this, the inspiral phase of the
subdominant modes can be described with the approximation from Eq. (17),
corrected with the complex amplitude phase contribution of Eq. (18) and
rotated to the original mode orientation through the inverse global rotation
of eq. (13):

$\phi_{lm}^{\mathrm{insp}}(t)=(m/2)\phi_{22}^{\mathrm{insp}}(t)+\Lambda_{lm}(t)-\Delta\phi_{lm}.$
(19)

Consequently, the mode frequency $\omega_{lm}(t)=\dot{\phi}_{lm}(t)$ is given
by:

$\omega_{lm}^{\mathrm{insp}}(t)=(m/2)\omega_{22}^{\mathrm{insp}}(t)+\dfrac{d}{dt}\Lambda_{lm}(t).$
(20)

Since the inspiral amplitude is already a complex quantity, in the model
implementation there is no need to compute the phase contribution
$\Lambda_{lm}(t)$ if the modes are constructed as the complex amplitude times
the complex phase exponential:

$h^{\mathrm{insp}}_{lm}(t)=H_{lm}^{\mathrm{insp}}(t)\exp[(m/2)\phi_{22}^{\mathrm{insp}}(t)-\Delta\phi_{lm}].$
(21)

Phases and frequencies computed from the modes described by Eq. (21) will
satisfy Eq. (19). As an example of the relation between $\Lambda_{lm}(t)$ and
the deviation from approximation (17), in Fig. 3 the time derivative of
$\Lambda_{lm}(t)$ and the deviation of the frequency from the scaling
approximation are shown for the $l=2,m=1$ and $l=3,m=3$ modes of a NR
simulation and the surrogate model NRHybSur3dq8 for the same parameters, where
it can be observed that $d\Lambda_{lm}(t)/dt$ correctly reproduces the
behaviour.

Figure 3: Frequency contribution from complex inspiral amplitude. A weighted
frequency difference is shown for the modes $(l=2,m=1)$ and $(l=3,m=3)$ with
respect to the $(l=2,m=2)$ frequency for the SXS NR simulation SXS:BBH:0290
with parameters $q=3$, $\chi_{1}=0.6$, $\chi_{2}=0.4$. A difference of zero
would correspond to the approximation from Eq. (17) being exact. It can be
seen that the contribution $\Lambda_{lm}$ from the complex inspiral amplitude
matches the behaviour of the deviation as compared with the NR simulation and
to the model NRHybSur3dq8.

###  II.3 Ringdown

For describing the modes in the ringdown region, we employ the proposal
described in Damour and Nagar (2014), that we briefly summarise here. The
damped emission of the final black hole can be well approximated in terms of
linear perturbations of the Kerr solution, for which analytical approximations
are known in the framework of perturbation theory in terms of a linear
combination of damped quasinormal modes (QNMs) Kokkotas and Schmidt (1999):

$h^{RD}_{lm}(t)=\sum_{n=1}^{\infty}c_{nlm}\exp[i\sigma_{nlm}(t-t_{0})],$ (22)

where indices $n$, $l$ and $m$ refer to the energy level and to the spin-
weighted spheroidal mode, $c_{lnm}$ are amplitude coefficients and
$\sigma_{nlm}$ is the complex frequency of the mode level, from which the
asymptotic final frequency $\omega_{nlm}^{\text{RD}}$ and damping frequency
$\alpha_{nlm}$ can be obtained:

$\displaystyle\omega_{nlm}^{\text{RD}}$ $\displaystyle=\Re(\sigma_{nlm}),$
(23a) $\displaystyle\alpha_{nlm}$ $\displaystyle=\Im(\sigma_{nlm}).$ (23b)

The $n=1$ level, also called the ground state QNM, is the dominant QNM at late
times, since is the most long-lived. Defining a QNM-rescaled mode as

$\bar{h}_{lm}(t)=e^{i\sigma_{1lm}(t)}h_{lm}(t),$ (24)

the early ringdown behaviour can be modelled with the following
phenomenological expressions:

$\begin{split}\bar{\omega}_{lm}(t)&=\omega_{lm}(t)-\omega^{\mathrm{RD}}_{1lm}\\\
&=c_{1}\dfrac{c_{2}(c_{3}e^{-c_{2}t}+2c_{4}e^{-2c_{2}t})}{1+c_{3}e^{-c_{2}t}+c_{4}e^{-2c_{2}\
t}},\end{split}$ (25)

$\begin{split}|\bar{h}_{lm}(t)|&=e^{\alpha_{1lm}(t-t^{peak}_{lm})}|h|\\\
&=d_{1}\tanh[d_{2}(t-t^{\mathrm{peak}}_{lm})+d_{3}]+d_{4},\end{split}$ (26)

where the following coefficients are set by imposing that the amplitude
matches the value at the peak amplitude time $t^{peak}_{lm}$, the frequency
matches the value at $t=0$ and that the amplitude derivative at the peak is
zero:

$\displaystyle c_{1}$
$\displaystyle=\dfrac{1+c_{3}+c_{4}}{c_{2}(c_{3}+2c_{4})}\left(\omega_{1lm}^{\mathrm{RD}}-\omega_{lm}(t=t^{\mathrm{peak}}_{22})\right),$
(27a) $\displaystyle c_{2}$ $\displaystyle=(\alpha_{2lm}-\alpha_{1lm})/2,$
(27b) $\displaystyle d_{1}$
$\displaystyle=H_{lm}(t=t^{\mathrm{peak}}_{lm})\dfrac{\alpha_{1lm}\cosh^{2}(d_{3})}{d_{2}},$
(27c) $\displaystyle d_{2}$ $\displaystyle=(\alpha_{2lm}-\alpha_{1lm})/2,$
(27d) $\displaystyle d_{4}$
$\displaystyle=H_{lm}(t=t^{\mathrm{peak}}_{lm})-d_{1}\tanh(d_{3}),$ (27e)

and $c_{3}$, $c_{4}$, $d_{3}$ are free coefficients calibrated with NR data.

However, in the present work we choose a different treatment of the phase
coefficients to improve accuracy in regions where the relation (27b) is not
well satisfied, particularly for high unequal mass systems. In our model,
$c_{2}$ is treated as a free coefficient to calibrate across parameter space
for each mode while maintaining simplicity and avoiding complicated
interactions between the different coefficients fits and $c_{4}$ is set to 0
(from Eq. (25) it can be seen that setting $c_{4}=0$ transforms the ansatz
into a $\tanh$ function, as for the amplitude ansatz). While it can be argued
that for the same reason the amplitude condition (27d) should be relaxed, as
is done in other models that employ the same ringdown description as SEOBNRv4
and SEOBNRv4HM Cotesta  _et al._ (2018), we have found that maintaining the
original amplitude treatment gives a sufficiently accurate ansatz, so in
favour of simplicity and to not complicate the behaviour of the ansatz across
parameter space we decided to maintain relation (27d).

One caveat that has been found during the construction of this model is that
for the $l=5,m=5$ mode the amplitude ansatz (26) does not seem to correctly
reproduce the available NR data. This could be due to numerical issues, since
in general the $l=5,m=5$ is weak and difficult to extract in numerical
simulations. However, this behaviour has been found in a significant number of
cases, pointing to some physical explanation, that could be the interaction of
higher overtones or some mode mixing effect with other modes that complicates
the description in the spherical harmonic basis representation. Investigation
of this issue and possible improvements of the ringdown modelling is left for
future work.

In a similar spirit, we have investigated the modelling of spherical harmonic
modes with mode-mixing contributions. In particular, the phenomenological
model IMRPhenomXHM includes the $l=3$, $m=\pm 2$ modes, which as discussed in
García-Quirós  _et al._ (2020) present an important contribution from the
$l=2$, $m=2$ mode in the ringdown regime. We have followed the strategy of
García-Quirós  _et al._ (2020) to apply our ringdown ansatz to the $l=3$,
$m=\pm 2$ modes in a decomposition into spin-weighted spheroidal harmonics,
and then rotate back to the spin-weighted spherical harmonic modes. As is well
known, spin-weighted spheroidal spherical harmonics provide a natural basis
for ringdown in the Kerr spacetime of the remnant black hole. We thus perform
the phenomenological fits for the inspiral and merger in the spherical basis,
and for the ringdown in the spheroidal basis, using the free coefficients from
Eq. (27) in that basis, before rotating the resulting spheroidal modes back to
the spherical basis needed for constructing the polarisations.

Therefore, it is required to calibrate $H_{lm}(t^{\mathrm{peak}}_{lm})$ and
$\omega_{lm}(t^{\mathrm{peak}}_{22})$ for each mode across parameter space.
The values of $\omega_{1lm}^{\mathrm{RD}}$, $\alpha_{1lm}$ and $\alpha_{2lm}$
can be obtained from perturbation theory as functions of the final mass and
final spin of the remnant black hole. We have performed 1D fits of
$\omega_{1lm}^{\mathrm{RD}}$, $\alpha_{1lm}$ and $\alpha_{2lm}$ as a function
of the dimensionless final spin (i.e, the Kerr parameter) employing the
available numerical data from Berti  _et al._ (2009). In order to evaluate
these fits, a calibrated function for the final mass and final spin across
parameter space is needed. There are several NR fits in the literature Tichy
and Marronetti (2008); Rezzolla  _et al._ (2008); Buonanno  _et al._ (2008);
Barausse  _et al._ (2012); Healy  _et al._ (2014); Hofmann  _et al._ (2016);
Jiménez-Forteza  _et al._ (2017); Varma  _et al._ (2019a) for the properties
of the remnant black hole, and this model employs the results from Jiménez-
Forteza  _et al._ (2017) which are also employed in the phenomenological
models IMRPhenomXHM and IMRPhenomXAS. Finally, for the three free coefficients
$d_{3}$, $c_{3}$ and $c_{4}$ a calibration to numerical relativity is
performed across parameter space.

###  II.4 Merger

To describe the behaviour in the merger regions of the amplitude and frequency
of each mode, we present phenomenological expressions based on hyperbolic
functions, which reproduce well the NR behaviour:

$\omega^{\mathrm{merger}}_{lm}(t)=\omega_{1lm}^{\mathrm{RD}}\left(1-\bar{\omega}^{\mathrm{merger}}_{lm}(t)\right),$
(28)

where

$\bar{\omega}^{\mathrm{merger}}_{lm}(t)=\sum_{k=0}^{k=4}a_{k}\text{arcsinh}^{k}(\alpha_{1lm}t),$
(29)
$\begin{split}H_{lm}^{\mathrm{merger}}(t)=&b_{0}+b_{1}\tau_{lm}^{2}+b_{2}\text{sech}^{1/7}(2\alpha_{1lm}\tau_{lm})\\\
&+b_{3}\text{sech}(\alpha_{1lm}\tau_{lm}),\end{split}$ (30)

and $\alpha_{1lm}$ is the damping frequency of the ground state QNM of the
$(l,m)$ mode. $a_{i}$, $b_{i}$ are phenomenological coefficients to be
determined by imposing continuity and differentiability at the boundaries with
the inspiral and ringdown regions and by imposing the ansatz to match a
particular frequency or amplitude value at a specific time
$t_{cp}^{\mathrm{merger}}=-25M$, and $\tau_{lm}=t-t^{\mathrm{peak}}_{lm}$
where $t^{\mathrm{peak}}_{lm}$ is the peak amplitude time of each mode. As an
illustrative example, we show in Fig. 4 the normalized amplitude of the
different modes for a particular SXS NR simulation and the corresponding
solution from Eq. (30), aligning the peak of the modes at $t=0$.

Figure 4: Amplitude of the different modes normalized to the peak amplitude
and peak aligned at $t=0$, for the SXS NR simulation SXS:BBH:1460 with
parameters $q=8$, $\chi_{1}=0.124$, $\chi_{2}=0.109$. Solid: SXS NR, dashed:
IMRPhenomTHM.

In order to obtain a value for the coefficients for each case, a linear system
of equations is imposed for guaranteeing continuity and differentiability at
the region boundaries, and a collocation point is placed at a fixed time of
both for the amplitude and the frequency:

$\displaystyle\omega^{\mathrm{merger}}_{lm}(t^{\omega}_{\mathrm{cut}})$
$\displaystyle=\omega^{\mathrm{insp}}_{lm}(t^{\omega}_{\mathrm{cut}}),$ (31a)
$\displaystyle\omega^{\mathrm{merger}}_{lm}(t^{\mathrm{peak}}_{22})$
$\displaystyle=\omega^{\mathrm{RD}}_{lm}(t^{\mathrm{peak}}_{22}),$ (31b)
$\displaystyle\omega^{\mathrm{merger}}_{lm}(t^{\omega}_{\text{merger},lm})$
$\displaystyle=\lambda^{\omega}_{\text{merger},lm},$ (31c)
$\displaystyle\frac{d}{dt}\omega^{\mathrm{merger}}_{lm}(t^{\omega}_{\mathrm{cut}})$
$\displaystyle=\frac{d}{dt}\omega^{\mathrm{insp}}_{lm}(t^{\omega}_{\mathrm{cut}}),$
(31d)
$\displaystyle\frac{d}{dt}\omega^{\mathrm{merger}}_{lm}(t^{\mathrm{peak}}_{22})$
$\displaystyle=\frac{d}{dt}\omega^{\mathrm{RD}}_{lm}(t^{\mathrm{peak}}_{22}),$
(31e)

and

$\displaystyle H^{\mathrm{merger}}_{lm}(t^{H}_{\mathrm{cut}})$
$\displaystyle=H^{\mathrm{insp}}_{lm}(t^{H}_{\mathrm{cut}}),$ (32a)
$\displaystyle H^{\mathrm{merger}}_{lm}(t^{\mathrm{peak}}_{lm})$
$\displaystyle=H^{\mathrm{RD}}_{lm}(t^{\mathrm{peak}}_{lm}),$ (32b)
$\displaystyle H^{\mathrm{merger}}_{lm}(t^{H}_{\mathrm{merger}})$
$\displaystyle=\lambda^{H}_{\mathrm{merger},lm},$ (32c)
$\displaystyle\frac{d}{dt}H^{\mathrm{merger}}_{lm}(t^{H}_{\mathrm{cut}})$
$\displaystyle=\frac{d}{dt}H^{\mathrm{insp}}_{lm}(t^{H}_{\mathrm{cut}}).$
(32d)

The inspiral amplitude described in Eq. (14) is a complex quantity as we have
discussed, however Eq. (30) is a real function. In order to propagate the
phase contribution from the complex inspiral amplitude to the merger, we
evaluate the value of $\Lambda_{lm}(t)$ from Eq. (18) at the boundary time
between the inspiral and merger regions,
$\Lambda_{lm}^{\text{merger}}=\Lambda_{lm}(t^{H}_{\text{cut}})$, and we add
this contribution as a constant phase offset to the merger phase jointly with
the global amplitude rotation that contains the information of the relative
orientation between the modes:

$\phi_{lm}^{\text{merger}}(t)=\Lambda_{lm}^{\text{merger}}-\Delta\phi_{lm}+\int_{t^{\prime}=t^{\omega}_{\text{cut},lm}}^{t}dt^{\prime}\omega^{\text{merger}}_{lm}(t^{\prime}).$
(33)

This allows to keep the sign from the real part of the inspiral amplitude as
well as the relative mode orientation, and since Eq. (30) is a real quantity
it also allows to track possible sign changes, i.e discontinuities of $\pi$ in
the phase that could occur in the merger region. We have found during the
validation of the model that this phenomenon indeed happens in some region of
the parameter space for the $l=2,\ m=1$ mode, for mass ratios close to
$m_{1}/m_{2}=3$, high positive primary spin $\chi_{1}$ and high spin
difference $\chi_{1}-\chi_{2}$. Although we have not encountered this
situation in NR simulations, two different models with very different
modelling procedures agree in reproducing this phenomenon, the surrogate model
NRHybSur3dq8 and the Effective-One-Body model SEOBNRv4HM. A representative
case of this region of parameter space is shown in Fig. 5. While it can be
seen that agreement between the different models is quantitatively poor, the
qualitative behaviour is shared among them.

Figure 5: Example of transition to negative amplitude for the $l=2$, $m=1$
mode, for a case with parameters $q=2.894$, $\chi_{1}=0.871$,
$\chi_{2}=-0.849$. Top panel: absolute value of the amplitude for the
different models, where the point can be observed where the absolute
amplitudes vanish. Mid panel: frequency becomes negative at the points of
vanishing amplitude. Bottom panel: phase suffers a discontinuity of $\pi$ at
the times of vanishing amplitude. (For NRHybSur3dq8 this transition is
softened due to the model specific construction, which does not allow for
sudden discontinuities.)

##  III Calibration across parameter space

Once the phenomenological ansätze for describing the amplitude and frequency
of each mode have been specified, there is a set of quantities that needs to
be specified in order to reconstruct the modes. In summary, the following
quantities have to be specified for each mode (a total of 16 for the $l=2$,
$m=2$ mode and 10 for the subdominant modes):

  * •

Six inspiral frequency collocation points for the $l=2$, $m=2$ mode:

$\\{\lambda^{\omega}_{i}\\}_{i=1,...,6}$ (34)

  * •

Three inspiral amplitude collocation points for each mode:

$\\{\lambda^{H}_{lm,i}\\}_{i=1,2,3}$ (35)

  * •

Two merger collocation points for each mode:

$\\{\lambda^{\omega}_{\text{merger},lm},\ \lambda^{H}_{\text{merger,lm}}\\}$
(36)

  * •

Peak amplitude and frequency at $t^{\text{peak}}_{22}$ for each mode:

$\\{\omega_{lm}(t^{\text{peak}}_{22}),\ H_{lm}(t^{\text{peak}}_{lm})\\}$ (37)

  * •

Three ringdown coefficients for each mode:

$c^{lm}_{2},\ c^{lm}_{3},\ d^{lm}_{3}.$ (38)

###  III.1 Calibration procedure

The above quantities are calibrated to a dataset of NR simulations, numerical
Teukolsky waveforms and SEOBNRv4 waveforms (details on this dataset are given
in the next subsection) as functions of the intrinsic binary parameters
$\\{\eta,\chi_{1},\chi_{2}\\}$. As in the previous work Estellés  _et al._
(2020) and as it is done in the phenomenological models IMRPhenomXAS Pratten
_et al._ (2020a) and IMRPhenomXHM García-Quirós  _et al._ (2020), we employ
the hierarchical data-driven fitting approach Jiménez-Forteza  _et al._ (2017)
to calibrate the parameter space functions, which we briefly summarise here.

The approach consists of constructing a full parameter space ansatz through a
sequence of fits to lower dimensional subspaces, which are typically more
densely populated with numerical simulations. First, 1-dimensional fits are
performed for the dependence on the symmetric mass-ratio for the non-spinning
subset of the dataset, and the spin dependence of equal black holes. In order
to capture in a simpler way the full two-dimensional spin-dependence, it is
re-parameterized in terms of a dominant effective spin $\hat{S}$ and the spin
difference $\Delta\chi$:

$\displaystyle\hat{S}=\frac{m_{1}^{2}\chi_{1}^{2}+m_{2}^{2}\chi_{2}^{2}}{m_{1}^{2}+m_{2}^{2}},$
(39a) $\displaystyle\Delta\chi=\chi_{1}-\chi_{2},$ (39b)

which were employed in the construction of the final state fits of Jiménez-
Forteza  _et al._ (2017) and in some of the fits constructed for IMRPhenomXAS
and IMRPhenomXHM models. For the dominant effective spin effects, one performs
another 1-dimensional fit for the equal spin data at a particular mass-ratio,
which typically is set to $q=1$ for the fits of quantities of $m-$even modes,
and to $q=2$ for quantities of the $m-$odd modes, since the odd modes vanish
for $q=1$ equal spin cases due to symmetry. With both 1D fits, one constrains
a 2D fit $f_{2D}(\eta,\hat{S})$ over the equal spin subset. From this fit,
residuals with the unequal spin cases are computed, and these residuals are
fitted as a function $h(\eta,S,\Delta\chi)$ of spin difference, effective spin
and mass-ratio. The general form of a parameter space fit for a specific
quantity is:

$\text{fit}(\eta,\hat{S},\Delta\chi)=g(\eta)f_{2D}(\eta,\hat{S})+h(\eta,S,\Delta\chi),$
(40)

where $g(\eta)$ is a scale function of the symmetric mass ratio for taking
into account known effects at the symmetric mass ratio boundaries, like
vanishing amplitude when $\eta\rightarrow 0$ (then $g(\eta)=\eta$) or
vanishing amplitude for the odd modes when both black holes are identical
($g(\eta)=\eta\delta m$). Further explanation and details on this procedure
can be found in Pratten  _et al._ (2020a) and García-Quirós  _et al._ (2020),
since the procedure followed here is very similar to the procedure followed in
the PhenomX family of waveform models.

###  III.2 Calibration dataset

The data employed for the model calibration consist of 531 BBH nonprecessing
quasi-circular NR simulations performed with the SpEC code and released in the
SXS Collaboration Catalog of NR simulations Boyle  _et al._ (2019), 15
simulations from the BAM code and 61 cases obtained from numerical solutions
of the perturbative Teukolsky equation through the numerical code
$\it{Teukode}$ Harms  _et al._ (2014). The set of SpEC simulations offers a
good coverage of the parameter space region below $q=10$, though the number of
equal spin cases above $q=4$ is reduced. For this reason, we incorporate the
BAM simulations at equal spin points placed at $q=4$ and $q=8$. We also employ
BAM simulations at $q=18$ for different values of the primary spin, which help
to connect the comparable masses region with the high mass ratio region of the
Teukolsky waveforms, which are placed at $q=80$, $q=200$ and $q=1000$ for
different spin values. However, we have treated the $q=18$ as equal spin
cases, since the effective spin parameter employed is almost the same assuming
that the secondary spin has the same value as the primary spin and since this
helps the conditioning of the 2D fits in $\\{\eta,\hat{S}\\}$. Figure 6 shows
the parameter space coverage for $\eta$ and the effective spin $\hat{S}$.

For the $l=2$, $m=2$ inspiral frequency collocation points
$\\{\lambda^{\omega}_{i}\\}$ that are placed at earlier times than the
starting time of a particular NR simulation, the value has been taken from
waveforms produced with the SEOBNRv4 model at the same parameters as in the
corresponding NR simulation. This procedure is completely equivalent to
employing data from the EOB-NR hybridised waveforms that were employed in the
calibration of IMRPhenomXAS and IMRPhenomXHM, since the waveforms were
hybridised employing also the SEOBNRv4 model.

For the SpEC simulations, we have selected the highest resolution level
available from which we have used the rhOverM_Asymptotic_GeometricUnits_CoM.h5
file, which contains the RZW strain with the center of mass correction as
explained in Boyle  _et al._ (2019). From this file, we select the $N=3$
extrapolated level for the inspiral and merger collocation points of the
different modes, but we select the $N=2$ extrapolated level for the peak and
ringdown quantities, since we have found that in general this level contains
less noise for the peak and post-peak regions of the waveform modes.

For the BAM simulations and the Teukolsky waveforms, EOB-NR hybrid waveforms
for each mode are employed. These hybrids, also used in the calibration of
IMRPhenomXAS and IMRPhenomXHM, reproduce the $\psi_{4}$ Weyl scalar, and we
obtain the strain through fixed-frequency integration Reisswig and Pollney
(2011). The fixed-frequency integration technique is based on Fourier-
transforming the $\psi_{4}$ waveform and performing the double time
integration in the Fourier domain, which is equivalent to multiplying the
waveform by a factor $-1/(4\pi^{2}f^{2})$, where $f$ is the Fourier frequency.
However, to eliminate the effect of spectral leakage and removing spurious
non-linear drifts, $f$ is set to a value $f_{0}$ in the region outside the
required frequency range of the final time domain strain waveform following
Reisswig and Pollney (2011). As discussed in Reisswig and Pollney (2011), the
cutoff frequency $f_{0}$ needs to be adjusted to the waveform. Too low values
of $f_{0}$ lead to noisier time domain waveform, whereas too high values
typically cause a loss of amplitude at late times.

Figure 6: Parameter space coverage of the NR and Teukolsky dataset employed
for the calibration.

##  IV Performance of the phenomenological model

In order to test and validate the construction and calibration of the model,
as presented in the previous sections, here we will present a series of tests.
First, we compare the model with the public LVCNR Catalog of NR simulations
and with the surrogate model NRHybSur3dq8, which accurately reproduces and
interpolates a hybridised version of the SXS NR Catalog, employing the usual
mismatch quantity for assessment. Comparison with the phenomenological model
IMRPhenomXHM is also presented, which allows to assess the validity of the
model in more challenging regions of parameter space. And finally, a
computation of the gravitational wave recoil, with a discussion of the
accuracy implications of this quantity, will be compared to public NR fits for
this quantity and with the results of other multimode models.

Mismatch is a common quantity in the gravitational wave modelling community to
assess the disagreement between two gravitational waveforms. It is defined in
terms of a noise-weighted inner product called overlap:

$(h_{1}(t)|h_{2}(t))=4\Re\int_{f_{\text{min}}}^{f_{\text{max}}}df\dfrac{\tilde{h}_{1}(f)\tilde{h}^{*}_{2}(f)}{S_{n}(f)}$
(41)

where $\tilde{h}(f)$ is the Fourier transform of $h(t)$, ∗ denotes complex
conjugation and $S_{n}(f)$ is the power spectral density of a detector at a
particular frequency $f$. From it, the faithfulness or match is obtained by
maximising the normalised overlap over a global phase and time translation:

$\mathcal{M}\equiv\max_{\phi_{0},\Delta
t}\dfrac{(h_{1}(t)|h_{2}(t))}{\sqrt{|h_{1}(t)||h_{2}(t)|}},$ (42)

where the norm of a waveform $|h(t)|$ is defined in terms of the overlap
product. Mismatch is the deviation of the match from unity:

$\mathcal{MM}\equiv 1-\mathcal{M}.$ (43)

Two completely equivalent waveform will have zero mismatch, while two
completely different waveform will have a mismatch close to unity.

###  IV.1 Comparison with LVCNR Catalogue

A first performance test of the model is the comparison with the dataset of NR
simulations from the LVCNR Waveform Catalog Schmidt  _et al._ (2017), in
particular with 110 simulations from the SXS Collaboration SXS Collaboration
(2019a) produced by the SpEC code SXS Collaboration (2019b). NR simulations
are in general the most trustworthy source of information about the
gravitational wave emission in the nonlinear highly dynamical regime of the
binary merger, and provide information of this regime to the different
calibrated waveform model families. Therefore, an important test for a NR
calibrated waveform model is to show a good agreement with numerical
relativity.

Figure 7: Mismatch distributions of IMRPhenomTHM polarisations with the LVCNR
SXS Catalog. Histograms for the minimum, mean and maximum mismatch for each
case are shown for the three different inclinations employed in the
comparison. It can be seen that for face-on configuration and $\pi/3$
inclination, all cases are below $1\%$ mismatch, while for edge-on
configurations ($\pi/2$ inclination) a few cases are above $1\%$ mismatch.

Figure 8: Single mode mismatches of IMRPhenomTHM with LVCNR SXS Catalog as a
function of the total mass of the binary. Red: non-spinning configurations.
Blue: spinning configurations. It can be seen that for the weaker modes $l=4$,
$m=4$ and $l=5$ and $m=5$ a good portion of cases show poor agreement for all
mass configurations, due to NR waveform quality.

We compute mismatches (Eq. (43)) between polarisations obtained from the LVCNR
Catalog simulations and the model presented in this work, for three different
inclinations $\iota\in\\{0,\pi/3,\pi/2\\}$ and different total masses of the
binary system, from $M=20M_{\odot}$ to $M=320M_{\odot}$ in steps of
$40M_{\odot}$. Results can be seen in Fig. 7, where distributions of the
minimum, mean and maximum mismatch for each case are shown for the three
different inclinations. A degradation in the results can be observed with
increasing inclination, since the subdominant modes contribution is becoming
larger as the systems go from face-on configurations (zero inclination) to
edge-on configurations (inclination $\pi/2$) where the dominant $l=2,m=2$ mode
is turned off. Nevertheless, even for edge-on configurations, the maximum
mismatches of only a few cases are above $1\%$ mismatch.

Mismatches for individual modes are also computed, for the same total mass
list as in the previous comparison. (For individual modes, there is no need to
compute at different inclinations, since this only scales the waveform
amplitude for single mode content.) Results in terms of the total mass of the
system are shown in Fig. 8 and mismatch distribution of minimum, mean and
maximum mismatch for each mode is shown in Fig. 9. The parameter space
distribution of mismatches greater than $1\%$ is shown in Fig. 10. A
degradation of mismatch results can be observed when going from the dominant
mode $l=2$, $m=2$ to the weakest mode $l=5$, $m=5$, which is consistent with
the fact that for weaker modes, the quality of NR waveforms is worse since
they are more difficult to extract from the global complex strain.

Figure 9: Single mode mismatch distributions of IMRPhenomTHM with LVCNR SXS
Catalog. For each mode, the minimum, mean and maximum mismatch distributions
are shown. Figure 10: Parameter space distribution of mismatch greater than
$1\%$ for the single mode comparison with the LVCNR Catalog. The dominant
$l=2$, $m=2$ mode is included but not shown since no mismatch greater than
$1\%$ was obtained.

###  IV.2 Comparison with other multimode models

Although the dataset of NR simulations is the most trustworthy set of
waveforms in the highly dynamical regime of binary mergers, it is a finite
dataset, with less than a thousand simulations for non-precessing
configurations, and they are also limited in length due to computational cost.
As a more extensive test of the model, in this section we present comparisons
with a NR-EOB hybrid surrogate model called NRHybSur3dq8, which has been
constructed to infer waveforms of any length within the bounds of its training
region, $q\leq 9$, $|\chi_{i}|\leq 0.91$, which essentially corresponds to the
parameter space coverage of NR simulations. To test the model also outside the
training boundaries of NRHybSur3dq8, in more challenging regions of parameter
space, mismatch results with the Fourier domain phenomenological model
IMRPhenomXHM are also presented, which allow to test the model up to extremal
spin configurations, and for arbitrary mass ratio (although in the following
results we have restringed to mass ratios below $q=20$).

With each model, single mode mismatches between the coincident modes in both
models and full polarisations will be computed. Waveforms are generated for
random parameter choices between some parameter bounds, employing the function
SimInspiralFD from the LALSuite libraries LIGO Scientific Collaboration
(2020), which provides the complex Fourier domain polarisations needed for the
noise weighted product (41). Gaussian colored noise is modeled by the power
spectral density estimation for the zero detuning high power configuration of
the advanced LIGO detectors Barsotti  _et al._ (2018), selecting a sampling
rate of $8192$ Hz (corresponding to a maximum frequency of $4096$ Hz) and with
minimum frequencies for the match calculation and starting waveform generation
frequencies listed in Table 1.

Targeted model | NRHybSur3dq8 | IMRPhenomXHM
---|---|---
$M$ range | $[10,200]M_{\odot}$ | $[10,200]M_{\odot}$
$q$ range | $[1,9]$ | $[1,20]$
$\chi_{i}$ range | $[-0.91,0.91]$ | $[-1,1]$
$f_{\text{sampling}}$ |  $8192$ Hz
$f_{\text{min}}$ |  $20$ Hz
$f_{\text{start}}$ (polarisations) | 8 Hz | 10 Hz
$f_{\text{start}}$ (single modes) |  $(m/2)20$ Hz
$S_{n}(f)$ | aLIGOzerodetHighpower
Table 1: Setting for mismatch comparison.

Figure 11: Mismatch distributions for individual modes for IMRPhenomTHM with
the non-precessing state-of-the-art models NRHybSur3dq8 (top) and IMRPhenomXHM
(bottom). A correlation between model agreement and mode hierarchy can be seen
in the results, with the dominant $l=2$, $m=2$ mode the most similar between
the models.

In Fig. 11 we can see the distributions of mismatches for each mode with the
corresponding mode of the model NRHybSur3dq8 and IMRPhenomXHM. From both
comparisons, we can see a hierarchy in the agreement of the modes, with the
dominant $l=2$, $m=2$ mode having the best agreement with both models. The
$l=2$, $m=1$ and $l=3$, $m=3$ modes also show good agreement, with a peak of
their distributions close to $0.1\%$ mismatch. Finally, the weak $l=4$, $m=4$
and $l=5$, $m=5$ modes have the worst agreement, but still have the peak and
most of their mismatch distributions below $1\%$. This can be understood in
terms of the original NR data that has been employed in the construction of
the three models, where the dominant $l=2$, $m=2$ mode is typically very well
reproduced by simulations, but data quality degrades with increasing $l$, as
the modes are weaker and it is more difficult to cleanly extract them from the
simulations.

Figure 12: Mismatch distributions of IMRPhenomTHM with NRHybSur3dq8 and
IMRPhenomXHM for multimode polarisations. Each model contains its full
available mode content. Top: histogram of the distributions for all cases.
Bottom: parameter space locations of cases with higher mismatches ($>3\%$ for
NRHybSur3dq8 and $>10\%$ for IMRPhenomXHM).

In Fig. 12 we can see the distributions of mismatches for the full multimode
polarisations compared with both models. We see that the peak of both
distributions is close to $0.1\%$ mismatch, with the distribution of
mismatches against NRHybSur3dq8 being shifted to smaller values of mismatch
and therefore better agreement found between this model and NRHybSur3dq8.
Nevertheless, the parameter range for the comparison is not the same between
both models. As we have explained, the validity region of the surrogate model
NRHybSur3dq8 coincides with the parameter space coverage of the bulk of NR
simulations, and therefore calibrated models are expected to agree more in
that region, while the parameter range of the comparison with IMRPhenomXHM has
been extended to higher mass ratios (up to $q=20$) and extremal spins,
precisely to compare the model outside the NR coverage region. In the bottom
pannel of Fig. 12 we can see the distribution of cases with a mismatch greater
than $3\%$ for the NRHybSur3dq8 comparison and greater than the $10\%$ for the
IMRPhenomXHM comparison (in order to capture the tail of both distributions).
It can be seen that poor mismatch is typically located on the left boundary of
the symmetric mass ratio (high unequal mass systems) and for large positive
primary spin values, having a preference also towards edge-on configurations,
where the multimode structure of the model is more relevant.

###  IV.3 Gravitational recoil

While mismatch studies provide a good measure for quantifying agreement
between different waveforms and datasets across parameter space, some aspects
relevant for the modelling of a multimode waveform are harder to extract from
a single scalar quantity such as the mismatch. In particular, while individual
modes can be compared intrinsically, i.e comparing their frequency and
amplitude evolution, the relative orientation of the modes that compose the
multimode structure of the polarisations cannot be tested with single mode
mismatch studies, and in the polarisation mismatch studies this is typically
shadowed by the quality of agreement of the dominant mode.

As an additional test for the multimode structure of the model, we provide
comparisons for the gravitational recoil prediction. Gravitational recoil is
an effect of asymmetric black hole binaries due to the anisotropic radiation
of gravitational waves, which generates a net emission of angular momentum.
Only multimode models of gravitational waves can capture this effect, since a
single mode emission will not be anisotropic. Emitted angular momentum in the
$\hat{\bm{n}}$ direction is given by:

$\dfrac{d\bm{P}}{dt}=\lim_{r\rightarrow\infty}\dfrac{r^{2}}{16\pi}\int
d\Omega(\dot{h}_{+}+\dot{h}_{\times})\hat{\bm{n}}.$ (44)

With the spherical harmonic decomposition of the polarisations of Eq. (1), the
sky-sphere integral in the $\hat{\bm{n}}$ direction can be performed on the
spherical harmonic basis. In terms of the parallel and perpendicular
components of the emitted linear momentum with respect to the direction of the
orbital angular momentum $\hat{\bm{L}}$, for non-precessing systems the
parallel component is $P_{||}=0$ due to the equatorial symmetry of the system,
and the perpendicular component (which lies on the orbital plane) is given in
terms of the mode time derivatives:

$\begin{split}\dfrac{d\bm{P}_{\perp}}{dt}&=\lim_{r\rightarrow\infty}\dfrac{r^{2}}{8\pi}\sum_{l,m}\dot{h}_{l,m}(a_{l,m}\dot{h}^{*}_{l,m+1}\\\
&+b_{l,-m}\dot{h}^{*}_{l-1,m-1}-b_{l+1,m+1}\dot{h}^{*}_{l+1,m+1}),\end{split}$
(45)

where the factors $a_{l,m}$ and $b_{l,m}$ are given by:

$\displaystyle a_{l,m}$ $\displaystyle=\dfrac{\sqrt{(l-m)(l+m+1)}}{l(l+1)},$
(46a) $\displaystyle b_{l,m}$
$\displaystyle=\frac{1}{2l}\sqrt{\dfrac{(l-2)(l+2)(l+m)(l+m-1)}{(2l-1)(2l+1)}}.$
(46b)

Figure 13: Comparison of final gravitational recoil kick velocity for
different models and a NR fit of this quantity. Top: non-spinning
configurations. Middle: Unequal spin configurations. Bottom: Equal spin
configurations. Shaded regions correspond to the uncertainty of the NR fit.

The final recoil velocity of the remnant will, due to linear momentum
conservation, be the time integral over the whole binary history of the
emitted linear momentum in the opposite direction:

$\bm{v}_{f}=-\int dt\dfrac{d\bm{P}_{\perp}}{dt},$ (47)

with $\bm{v}_{f}$ being a complex number with the magnitude of the recoil kick
and its direction in the orbital plane.

A prediction for the recoil velocity is a test well suited for the multimode
structure of a model. Eq. (45) can be seen as a sum of pairs of interacting
modes with different mode number $m=m^{\prime}\pm 1$, and the integral of Eq.
(47) can be split into different contributions from each pair of modes. To
accurately predict the final recoil velocity, each of these contributing
integrals has to be accurate, which is a test for the intrinsic frequency and
amplitude of the modes; but also the total sum of the different terms has to
be accurate, and this is a test for the different phase relations between the
modes.

Here we compute the magnitude $v_{f}=|\bm{v}_{f}|$ of the recoil velocity for
three different models for different spin configurations as a function of the
symmetric mass ratio of the system, and we compare with an existing NR fit of
this quantity. Results for the three models and the Healy & Lousto NR fit
Healy and Lousto (2017) prediction can be seen in Fig. 13. A similar trend to
the other models and the NR fit can be seen to be reproduced by IMRPhenomTHM.
For the spinning configurations, most of the computed values are inside the
uncertainty interval of the NR fit, though for the non-spinning configurations
the recoil is in general a bit underestimated. It can be seen, however, that
in the underestimated region our model agrees quite well with NRHybSur3dq8.
This difference with respect to the NR fit and the Fourier domain model
IMRPhenomXHM could be due to the different mode content of the models.

##  V Performance in parameter estimation

One of the most important applications of gravitational wave models is the
parameter estimation for detected gravitational wave events, to infer the
physical properties of the source systems. In this section, we test the
performance of the model in parameter estimation for simulated injected
signals and for real gravitational wave detections, employing the Bayesian
inference suite specialised in gravitational wave parameter estimation Bilby
Ashton  _et al._ (2019), in a parallelizable pipeline suited for HPC
parallelBilby Smith  _et al._ (2020), which performs a stochastic sampling of
the likelihood function across parameter space employing the nested sampling
algorithm dynesty Speagle (2020).

###  V.1 NR Injection recovery

As a first test of the performance of the model in parameter estimation, we
test the parameter recovery for two synthetic signals where a NR simulation
from the LVCNR Catalog has been injected into Gaussian zero-mean noise. We
select the SXS non-precessing simulations SXS:BBH:0110, as it was done for the
validation of IMRPhenomXHM in García-Quirós  _et al._ (2020), and
SXS:BBH:0237. In Table 2 we list the intrinsic parameters of each simulation
as well as the extrinsic parameters chosen for the synthetic signal, followed
by the recovered value and its uncertainty. Parameter estimation has been
performed with parallelBilby and Dynesty settings of nlive=2048, walks=200,
nact=5 and maxmcmc=10000.

Parameters | SXS:BBH:0110 | SXS:BBH:0237
---|---|---
Injected | Recovered | Injected | Recovered
$M_{\mathrm{tot}}/M_{\odot}$ | $62.0$ | $61.28^{+5.32}_{-5.34}$ | $45.0$ | $46.38^{+3.93}_{-2.95}$
$q$ | $0.2$ | $0.21^{+0.06}_{-0.04}$ | $0.5$ | $0.43^{+0.17}_{-0.11}$
$\chi_{\text{eff}}$ | $0.42$ | $0.41^{+0.07}_{-0.07}$ | $-0.2$ | $-0.2^{+0.13}_{-0.13}$
$\chi_{1z}$ | $0.5$ | $0.5^{+0.07}_{-0.1}$ | $-0.6$ | $-0.23^{+0.21}_{-0.26}$
$\chi_{2z}$ | $0.0$ | $-0.05^{+0.6}_{-0.56}$ | $0.6$ | $-0.1^{+0.42}_{-0.46}$
$\iota$ | $\pi/3$ | $0.9^{+0.22}_{-0.26}$ | $\pi/3$ | $0.71^{+0.33}_{-0.37}$
$\phi_{\text{ref}}$ | $0$ | $4.97^{+1.14}_{-4.77}$ | $0.8$ | $1.63^{4.48}_{-1.47}$
$\alpha$ | $1.4$ | $1.4^{+0.02}_{-0.02}$ | $3.81$ | $3.82^{+0.04}_{-0.04}$
$\delta$ | $-0.6$ | $-0.6^{+0.09}_{-0.1}$ | $0.63$ | $0.6^{+0.04}_{-0.05}$
Table 2: NR injection recovery. For each simulated signal, corresponding to
the NR simulations SXS:BBH:0110 and SXS:BBH:0237, the injected value and the
recovered parameter with their uncertainties (90$\%$ confidence intervals) are
listed. The listed parameters correspond to the total mass $M_{\mathrm{tot}}$
in solar masses $M_{\odot}$, mass ratio $q=m_{2}/m_{1}<1$, effective spin
parameter $\chi_{\text{eff}}$, z-components of the dimensionless spin vector
of each black hole ($\chi_{1z}$ and $\chi_{2z}$), inclination $\iota$,
coalescence phase $\phi_{\text{ref}}$ at the reference frequency, right
ascension $\alpha$, and declination $\delta$.

As we can see from the results of Table 2 and the posteriors shown in Fig. 14
and Fig. 15, the model is able to reproduce quite well both the intrinsic and
the extrinsic parameters from both simulated events. The main difference in
intrinsic parameter recovery is seen for the mass ratio and individual spin
components of the SXS:BBH:0237 injection. For the mass ratio, as a cross-
check, a injection recovery with the model IMRPhenomXHM was also performed,
confirming a similar bias in the result, possibly due to a lack of mode
content in the models with respect to the full NR waveform. But in any case
the injected value is still inside the $90\%$ confidence interval (CI) of the
recovered posterior. Regarding the individual spin components, although they
are not well recovered for SXS:BBH:0237, the effective spin parameter
$\chi_{\mathrm{eff}}$, which is the dominant spin parameter in the inspiral PN
expressions, is very well recovered for both injections.

Figure 14: SXS:BBH:0110 recovered posteriors.

Regarding extrinsic parameter recovery, it can be seen from Table (2 that the
sky location, parameterised as usual by right ascension and declination, is
well recovered for both simulated signals, and the recovered inclination is
inside the $90\%$ CI for both signals, with the result for SXS:BBH:0237
slightly more biased. The orbital phase is generally not a very well recovered
quantity, but the posterior distribution indicates a good recovery trend.

Figure 15: SXS:BBH:0237 recovered posteriors.

###  V.2 Real events

After testing parameter estimation performance of the model on synthetic
injected signals, a test on real gravitational wave data is performed. We have
selected two events from the first GWTC-1 catalog Abbott  _et al._ (2019) of
detections from the first and second observing runs of the Advanced LIGO and
Advanced Virgo detectors: GW150914, which was the first gravitational wave
detection Abbott  _et al._ (2016b), and GW170729, which has previously been
studied with various higher-mode models due to its high mass Chatziioannou
_et al._ (2019). We have employed the public 32 seconds strain data from GWOSC
LIGO Scientific Collaboration, Virgo Collaboration (2019a) sampled at
$4\text{kHz}$ around each event, with the public power spectral densities
estimated for those segments of data LIGO Scientific Collaboration, Virgo
Collaboration (2019b), and the public calibration envelopes of the Hanford and
Livingston detectors (for GW150914) and the Virgo detector (for GW170729) LIGO
Scientific Collaboration, Virgo Collaboration (2019c).

Parameter estimation has been performed with the parallelBilby pipeline, with
the following settings for the Dynesty sampler: nlive=2048, walks=200, nact=10
and maxmcmc=15000. Results for GW150914 are compared in Fig. 16 with the
published results from GWTC-1 Abbott  _et al._ (2019) obtained for this event
with the phenomenological waveform model IMRPhenomPv2, which does not contain
subdominant spherical harmonic modes information in its non-precessing core
model, and allows for generic spin configurations through the twisting-up
technique of the non-precessing model to describe precessing systems. Due to
the different physical content of both models, complete agreement in the
results is not expected, but results are consistent with the published ones.
Results for GW170729 are shown in Fig. 17, compared with the results of
Chatziioannou  _et al._ (2019) where the non-precessing multimode models
IMRPhenomHM and SEOBNRv4HM were employed, and a good consistency can be seen
between the three models.

Figure 16: GW150914 recovered posteriors.

Figure 17: GW170729 recovered posteriors.

Apart from the results for these two real events, the model was employed in
the re-analysis Colleoni  _et al._ (2020) of the O3a event GW190412 Abbott
_et al._ (2019), which we analysed with the last generation of frequency
domain phenomenological waveform models, the IMRPhenomX family and the
dominant mode and multimode time domain model presented in this work.
Excellent consistency was found in the results when comparing the last
generation of calibrated multimode non-precessing waveform models for the
different model families, IMRPhenomXHM, SEOBNRv4HM, NRHybSur3dq8 and
IMRPhenomTHM, while the previous generation of phenomenological waveform
models, IMRPhenomHM, based on an analytical uncalibrated rescaling of the
dominant mode model IMRPhenomD, showed higher discrepancies.

###  V.3 Benchmarks

In a similar way as we have done with previous models Pratten  _et al._
(2020b), we estimate the speed of our new waveform model, IMRPhenomTHM,
compared to other non-precessing frequency- and time-domain waveform models
including higher order modes: IMRPhenomXHM García-Quirós  _et al._ (2020),
IMRPhenomHM London  _et al._ (2018), SEOBNRv4HM Cotesta  _et al._ (2018),
SEOBNRv4HM$\\_$ROM Cotesta  _et al._ (2020) and NRHybSur3dq8 Varma  _et al._
(2019b). Thus, we show in Fig. 3 benchmarking results for an arbitrarily
chosen non-precessing case, $(q,\chi_{1z},\chi_{2z})=(3,0.5,-0.3)$, in a
frequency range from 10 to 2048 Hz by comparing different waveform models. The
timing test is performed with GenerateSimulation (included in the
LALSimulation package), averaging over 100 repetitions. In the left panel the
dependency on the total mass is shown for a range between
$(3,\>300)\>M_{\odot}$ for the waveforms in the time-domain generated with the
ChooseTDWaveform function (included in LALSimulation) and a fixed sample rate
of 4096 Hz. In the case of frequency-domain models the curves also include the
time to perform an inverse Fourier transform. Furthermore, for the time-domain
approximants one has to take into account the fact that in order to have all
the higher order modes starting at the desired minimum frequency of the
waveform for analysis purposes, one has to start the waveform generation
earlier. This is caused by the fact that higher order modes with $m>2$ have a
higher frequency content for the same array of times. This effect is
encapsulated in the amporder parameter which changes the starting frequency of
the waveform generation in the following way,

$f_{\text{start}}=\frac{2f_{\text{min}}}{\texttt{amporder}+2}.$ (48)

Hence, in our timing test we loop also over different values of the amporder
parameter to evaluate its effect on the speed of the waveform generation. The
right plot in Fig. 3 shows the timing test performed in the Fourier domain by
calling the distinct waveform models using the SimInspiralFD function
(included in LALSimulation). In this plot the total mass is varied as in the
previous panel, but the frequency resolution is computed internally to adapt
to the length of the waveform. Moreover, in this plot we have replaced the
time-domain SEOBNRv4 model by its frequency-domain counterpart
SEOBNRv4HM$\\_$ROM.

From the left plot of Fig. 3 one can conclude that IMRPhenomTHM performs
faster than the other two time-domain waveform models for total masses above
$30M_{\odot}$, while for low total masses ($<30M_{\odot}$) the SEOBNRv4HM can
not be generated and IMRPhenomTHM becomes comparable to NRHybSur3dq8. We note
here also the important effect of the amporder parameter: as expected, the
higher its value, the lower the starting frequency, and thus the generation of
the time-domain model becomes slower. However, even with that we observe that
IMRPhenomTHM becomes comparable to the frequency domain models for total
masses $>30M_{\odot}$, or even faster in the case of IMRPhenomHM. This is due
to the fact that the symbolic phenomenological expressions constituting the
core of IMRPhenomTHM make the model computationally efficient. Similar
conclusions can be extracted from the right panel of Fig. 3, where one
observes that while for total masses below $30M_{\odot}$ the model is slower
than the frequency-domain SEOBNRv4HM$\\_$ROM and IMRPhenomXHM models, for
higher total masses it becomes comparable and even faster than
SEOBNRv4HM$\\_$ROM, but never reaching the speed of IMRPhenomXHM, which
benefits from the multibanding technique Vinciguerra  _et al._ (2017); García-
Quirós  _et al._ (2020) to become the fastest higher order mode waveform model
in the frequency domain.

We have also estimated the efficiency of our model IMRPhenomTHM compared to
other models with higher order modes by computing their mean likelihood
evaluation time in the LALInference Bayesian parameter estimation code Veitch
_et al._ (2015). This test has been performed choosing an equal-mass
configuration with 100 different total masses in two different mass ranges: a
low mass range with range $[M_{\min},M_{\max}]=[10,50]M_{\odot}$, and high
mass range with $[M_{\min},M_{\max}]=[50,100]M_{\odot}$. Dimensionless spin
magnitudes are distributed randomly between $[0,0.99]$ with a random isotropic
distribution of spin vectors, and a reference frequency of $20$ Hz is chosen.
Furthermore, two different segment lengths of $\Delta T=4$ s, 8 s are studied
as they are typical for the currently detected BBH GW signals Abbott  _et al._
(2019, 2020c). For the likelihood evaluations we employ the Advanced LIGO zero
detuned power spectral density Barsotti  _et al._ (2018) as we also do for our
match calculations in Sec. IV and NR injection studies in Sec. V.1. For each
total mass we perform 100 likelihood evaluations with randomly chosen spin
configurations. In Table 3 the average of these $10^{4}$ likelihood
evaluations for each model is shown.

The numbers in Table 3 are consistent with the results of Fig. 3 as they show
that IMRPhenomTHM is the most efficient time-domain model including higher
order modes, and for higher total masses it becomes faster or comparable to
some frequency domain models like IMRPhenomHM and SEOBNRv4HM$\\_$ROM. This
shows that IMRPhenomTHM is perfectly suited for parameter estimation analysis
of high mass GW events like the ones detected by the LIGO and Virgo
collaboration Abbott  _et al._ (2019, 2020c) and a valuable tool to complement
and disentangle systematics coming from the use of frequency domain waveform
models in parameter estimation studies of GW events.

Figure 18: Average waveform evaluation time for different multimode non-
precessing waveform models employing the GenerateSimulation interface of
LALSimulation. Left: Time-domain waveform evaluation time as a function of the
total mass of the system. Right: Frequency-domain waveform evaluation time as
a function of the total mass of the system.

$[M_{\min},M_{\max}]$ | $\Delta T$ | IMRPhenomTHM | IMRPhenomXHM | IMRPhenomHM | SEOBNRv4HM | SEOBNRv4HM$\\_$ROM | NRHybSur3dq8
---|---|---|---|---|---|---|---
[50,100] | 4 s | 9.5 | 5.0 | 67.2 | 1415.9 | 29.5 | 221.1
8 s | 9.8 | 6.6 | 121.4 | 1430.2 | 35.2 | 230.2
[10,50] | 4 s | 55.6 | 13.6 | 72.6 | 1630.4 | 38.7 | 307.1
8 s | 69.1 | 20.6 | 166.0 | 1619.8 | 65.5 | 326.8

Table 3: Mean likelihood evaluation time in milliseconds for several models
including higher order modes for equal masses. The numbers represent averages
over two different total mass ranges
$[M_{\min},M_{\max}]=\\{[10,50],[50,100]\\}M_{\odot}$ and random spin
orientations and magnitudes. The first column indicates the total mass range
in which the models are evaluated and the second one specifies the data
analysis segment length in seconds used for the calculations.

##  VI Conclusions

In this work a new phenomenological model for non-precessing binary black hole
signals has been presented, based on the construction of the IMRPhenomT model
previously presented in Estellés  _et al._ (2020). It incorporates additional
spherical harmonic modes which are a requirement of current gravitational wave
observations, and its construction is native in the time domain,
differentiating it with respect to other phenomenological models from the
current and past generations.

It has been shown that the model accurately reproduces the LVCNR Catalog of
Numerical Relativity simulations, and that there is good agreement with other
state of the art non-precessing multimode models like NRHybSur3dq8 and
IMRPhenomXHM. The model also reproduces accurate results for the final kick
velocity due to anisotropic gravitational wave emission, as compared to
existing NR fits for this quantity and to the results of other multimode
models, which is an exigent test of the correct global multimode structure of
the waveforms. The model can be used to accurately estimate parameters from
injected simulated signals and to produce consistent results for public
gravitational wave detections, being also competitively fast for waveform
generation and likelihood evaluation, despite its native time domain nature,
which requires to perform a numerical Fast Fourier Transform in order to
employ the model for parameter estimation. As a further proof of the
reliability of the model, results from a preliminary version were employed in
the re-analysis of the gravitational wave event GW190412, showing very good
consistency with current generation calibrated non-precessing multimode models
and showing likelihood evaluation times comparable to the fastest
phenomenological model IMRPhenomXHM.

The model presented in this work, IMRPhenomTHM, complements the so-called
“generation X” of phenomenological waveform models, which contains the fastest
and most accurate family of frequency domain phenomenological models
IMRPhenomXAS, IMRPhenomXHM, IMRPhenomXP and IMRPhenomXPHM. In a future work,
an extension to model precessing systems, IMRPhenomTPHM will be presented,
based on the features already introduced in Estellés  _et al._ (2020) for
extending the dominant mode model IMRPhenomT into IMRPhenomTP. Having a
complementary time domain phenomenological framework will help to disentangle
possible modelling issues of the frequency domain framework, particulary
elimitating the Stationary Phase Approximation from the precession
description, and allowing to incorporate further methods of reproducing the
precessing angular momenta dynamics. As a first step, the non-precessing model
IMRPhenomTHM fulfills the requirements for the current state-of-the-art non-
precessing waveform models.

## Acknowledgements

We thank the reviewers of the LIGO Scientific Collaboration, Maria Haney,
Jonathan Thompson, Eleanor Hamilton and Jacob Lange for reviewing the LALSuite
code implementation, and for carefully reading the manuscript and valuable
feedback. We thank Alessandro Nagar, Sebastiano Bernuzzi and Enno Harms for
giving us access to $\it{Teukode}$ Harms  _et al._ (2014), which was used to
generate our extreme-mass-ratio waveforms.

This work was supported by European Union FEDER funds, the Ministry of
Science, Innovation and Universities and the Spanish Agencia Estatal de
Investigación grants PID2019-106416GB-I00/AEI/10.13039/501100011033,
FPA2016-76821-P, RED2018-102661-T, RED2018-102573-E, FPA2017-90687-REDC,
Vicepresidència i Conselleria d’Innovació, Recerca i Turisme, Conselleria
d’Educació, i Universitats del Govern de les Illes Balears i Fons Social
Europeu, Comunitat Autonoma de les Illes Balears through the Direcció General
de Política Universitaria i Recerca with funds from the Tourist Stay Tax Law
ITS 2017-006 (PRD2018/24), Generalitat Valenciana (PROMETEO/2019/071), EU COST
Actions CA18108, CA17137, CA16214, and CA16104, and the Spanish Ministry of
Education, Culture and Sport grants FPU15/03344 and FPU15/01319. M.C.
acknowledges funding from the European Union’s Horizon 2020 research and
innovation programme, under the Marie Skłodowska-Curie grant agreement No.
751492. D.K. is supported by the Spanish Ministerio de Ciencia, Innovación y
Universidades (ref. BEAGAL 18/00148) and cofinanced by the Universitat de les
Illes Balears. The authors thankfully acknowledge the computer resources at
MareNostrum and the technical support provided by Barcelona Supercomputing
Center (BSC) through Grants No. AECT-2020-2-0015, AECT-2020-1-0025,
AECT-2019-3-0020, AECT-2019-2-0010, AECT-2019-2-0017, AECT-2019-1-0022, from
the Red Española de Supercomputación (RES).

Authors also acknowledge the computational resources at the cluster CIT
provided by LIGO Laboratory and supported by National Science Foundation
Grants PHY-0757058 and PHY-0823459. This research has made use of data
obtained from the Gravitational Wave Open Science Center LIGO Scientific
Collaboration, Virgo Collaboration (2019a), a service of LIGO Laboratory, the
LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by
the U.S. National Science Foundation. Virgo is funded by the French Centre
National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale
della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by
Polish and Hungarian institutes.

##  Appendix A Post-Newtonian amplitudes

In the following we show the PN expressions used to construct the inspiral
amplitudes as described in Sec. II.2.2. The expressions for the different
modes in Eq. (12) are based on 3PN order expressions from Blanchet  _et al._
(2008) with 2PN spin corrections from Buonanno  _et al._ (2013) and 1.5PN
contributions from Arun  _et al._ (2008). For the amplitude of the $(2,2)$
mode, we also include 3.5PN corrections from Faye  _et al._ (2012). This is
mostly equivalent to the PN amplitudes employed by the frequency domain model
IMRPhenomXHM which are listed in Appendix E of García-Quirós  _et al._ (2020),
except that the expressions listed there already contain the stationary phase
approximation contributions needed for the frequency domain representation.
The difference to García-Quirós  _et al._ (2020) is that known 3PN
contributions to the $l=3$, $m=3$ mode have not been included here since we
find that a better performance with the calibrated higher order coefficients
was achieved without including these 3PN terms. Here they are listed after
factoring out the common prefactor $2\eta\sqrt{16/5}x(t)$:

$\begin{split}\bar{H}_{22}^{PN}(x(t))=&1+\Big{(}\frac{55\eta}{42}-\frac{107}{42}\Big{)}x+\Big{[}-\frac{2\delta\chi_{1}}{3\left(\frac{1-\delta}{2}+\frac{\delta+1}{2}\right)}+\frac{2\delta\chi_{2}}{3\left(\frac{1-\delta}{2}+\frac{\delta+1}{2}\right)}+\frac{2\eta\chi_{1}}{3}+\frac{2\eta\chi_{2}}{3}-\frac{2\chi_{1}}{3}-\frac{2\chi_{2}}{3}+2\pi\Big{]}x^{3/2}\\\
&+\Big{[}\frac{\delta\chi_{1}^{2}}{2}-\frac{\delta\chi_{2}^{2}}{2}+\frac{2047\eta^{2}}{1512}-\eta\chi_{1}^{2}-\eta\chi_{2}^{2}+2\eta\chi_{1}\chi_{2}-\frac{1069\eta}{216}+\frac{\chi_{1}^{2}}{2}+\frac{\chi_{2}^{2}}{2}-\frac{2173}{1512}\Big{]}x^{2}\\\
&+\Big{(}\frac{34\pi\eta}{21}-24i\eta-\frac{107\pi}{21}\Big{)}x^{5/2}-\Big{[}-\frac{856\gamma}{105}+\frac{114635\eta^{3}}{99792}-\frac{20261\eta^{2}}{2772}+\frac{41\pi^{2}\eta}{96}-\frac{278185\eta}{33264}\\\
&-\frac{428}{105}\log(16x)+\frac{2\pi^{2}}{3}+\frac{428i\pi}{105}+\frac{27027409}{646800}\Big{]}x^{3}+\Big{(}\frac{40\pi\eta^{2}}{27}-\frac{4066i\eta^{2}}{945}-\frac{2495\pi\eta}{378}+\frac{14333i\eta}{162}\\\
&-\frac{2173\pi}{756}\Big{)}x^{7/2},\end{split}$ (49)

$\begin{split}\bar{H}_{21}^{PN}(x(t))=&\frac{1}{3}i\delta\sqrt{x}-\frac{1}{2}i(m_{1}\chi_{1}+m_{2}\chi_{2})x-i\delta\Big{(}\frac{17}{84}-\frac{5\eta}{21}\Big{)}x^{3/2}\\\
&+\Big{[}\frac{1}{6}\delta(1+\log(16))+i\Big{(}\delta\frac{\pi}{3}+\frac{11}{28}\delta\eta(\chi_{1}+\chi_{2})-\frac{1}{6}(m_{1}\chi_{1}+m_{2}\chi_{2})+\frac{205}{84}(\chi_{1}-\chi_{2})\Big{)}\Big{]}x^{2}\\\
&-i\Big{(}\frac{43\delta}{378}+\frac{509\delta\eta}{378}-\frac{79\delta\eta^{2}}{504}\Big{)}x^{5/2}-\Big{[}\frac{17\delta}{168}(1+\log(16))-\delta\eta\Big{(}\frac{353}{84}+\frac{\log(2)}{7}\Big{)}+i\pi\delta\Big{(}\frac{\eta}{14}-\frac{17}{84}\Big{)}\Big{]}x^{3},\end{split}$
(50)

$\begin{split}\bar{H}_{33}^{PN}(x(t))=&-\frac{3}{4}i\sqrt{\frac{15}{14}}\delta\sqrt{x}+\Big{(}3i\sqrt{\frac{15}{14}}\delta-\frac{3}{2}i\sqrt{\frac{15}{14}}\delta\eta\Big{)}x^{3/2}\\\
&+\Big{[}-\frac{3}{16}i\sqrt{\frac{3}{70}}\left(5\chi_{1}(\delta(5\eta-4)+19\eta-4)+5\chi_{2}(5\delta\eta-4\delta-19\eta+4)+12\delta\left(-7i+5\pi\right.\right.\\\
&\left.\left.+10i\log\left(3/2\right)\right)\right)\Big{]}x^{2}+\Big{(}-\frac{887}{88}i\sqrt{\frac{3}{70}}\delta\eta^{2}+\frac{919}{22}i\sqrt{\frac{3}{70}}\delta\eta-\frac{369}{88}i\sqrt{\frac{3}{70}}\delta\Big{)}x^{5/2},\end{split}$
(51)

$\begin{split}\bar{H}_{44}^{PN}(x(t))=&\frac{8}{9}\sqrt{\frac{5}{7}}(3\eta-1)x-\frac{8}{9}\sqrt{\frac{5}{7}}\left(-\frac{175\eta^{2}}{22}+\frac{1273\eta}{66}-\frac{593}{110}\right)x^{2}\\\
&-\frac{8}{9}\sqrt{\frac{5}{7}}\Big{[}-12\pi\eta+i\left(\eta\left(\frac{1193}{40}-24\log(2)\right)-\frac{42}{5}+8\log(2)\right)+4\pi\Big{]}x^{5/2}\\\
&-\frac{8}{9}\sqrt{\frac{5}{7}}\left(-\frac{226097\eta^{3}}{17160}+\frac{146879\eta^{2}}{2340}-\frac{1088119\eta}{28600}+\frac{1068671}{200200}\right)x^{3},\end{split}$
(52)

$\begin{split}\bar{H}_{55}^{PN}(x(t))=&-\frac{3}{2}i\sqrt{\frac{15}{14}}\delta(\eta-2)x^{3/2}-\frac{1}{88}i\sqrt{\frac{3}{70}}\delta\left(887\eta^{2}-3676\eta+369\right)x^{5/2}\\\
&+\frac{\delta\left(-3645i\pi(3\eta-8)+\eta\left(21870\log\left(\frac{3}{2}\right)-96206\right)-5832\left(10\log\left(\frac{3}{2}\right)-7\right)\right)}{216\sqrt{210}}x^{3}.\end{split}$
(53)

## References

  * Veitch  _et al._ (2015) J. Veitch _et al._ , Phys. Rev. D91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Ashton  _et al._ (2019) G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter, B. Goncharov, and et al., The Astrophysical Journal Supplement Series 241, 27 (2019).
  * Aasi  _et al._ (2015) J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, and et al., Classical and Quantum Gravity 32, 074001 (2015).
  * Acernese  _et al._ (2014) F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, and et al., Classical and Quantum Gravity 32, 024001 (2014).
  * Pretorius (2005) F. Pretorius, Physical Review Letters 95 (2005), 10.1103/physrevlett.95.121101.
  * Campanelli  _et al._ (2006) M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Physical Review Letters 96 (2006), 10.1103/physrevlett.96.111101.
  * Baker  _et al._ (2006) J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), arXiv:gr-qc/0511103 [gr-qc] .
  * Husa  _et al._ (2016) S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D93, 044006 (2016), arXiv:1508.07250 [gr-qc] .
  * Khan  _et al._ (2016) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Hannam  _et al._ (2014) M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, _et al._ , Phys.Rev.Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Bohé  _et al._ (2016) A. Bohé, M. Hannam, S. Husa, F. Ohme, M. Puerrer, and P. Schmidt,  _PhenomPv2 - Technical Notes for LAL Implementation_, Tech. Rep. LIGO-T1500602 (LIGO Project, 2016).
  * Pürrer (2016) M. Pürrer, Phys. Rev. D 93, 064041 (2016), arXiv:1512.02248 [gr-qc] .
  * Abbott  _et al._ (2016a) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241102 (2016a), arXiv:1602.03840 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Class. Quant. Grav. 34, 104002 (2017), arXiv:1611.07531 [gr-qc] .
  * Abbott  _et al._ (2019) B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. Adhikari, V. Adya, C. Affeldt, and et al., Physical Review X 9 (2019), 10.1103/physrevx.9.031040.
  * Chatziioannou  _et al._ (2019) K. Chatziioannou, R. Cotesta, S. Ghonge, J. Lange, K. K. Ng, J. Calderón Bustillo, J. Clark, C.-J. Haster, S. Khan, M. Pürrer, and et al., Physical Review D 100 (2019), 10.1103/physrevd.100.104015.
  * Abbott  _et al._ (2020a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 102, 043015 (2020a), arXiv:2004.08342 [astro-ph.HE] .
  * Abbott  _et al._ (2020b) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 896, L44 (2020b), arXiv:2006.12611 [astro-ph.HE] .
  * LIGO Scientific Collaboration (2020) V. C. LIGO Scientific Collaboration, “Gwtc-2: Compact binary coalescences observed by ligo and virgo during the first half of the third observing run,” (2020), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2018) B. Abbott _et al._ (KAGRA, LIGO Scientific, VIRGO), Living Rev. Rel. 21, 3 (2018), arXiv:1304.0670 [gr-qc] .
  * Estellés  _et al._ (2020) H. Estellés, A. Ramos-Buades, S. Husa, C. García-Quirós, M. Colleoni, L. Haegel, and R. Jaume, “Imrphenomtp: A phenomenological time domain model for dominant quadrupole gravitational wave signal of coalescing binary black holes,” (2020), arXiv:2004.08302 [gr-qc] .
  * Pratten  _et al._ (2020a) G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni, A. Ramos-Buades, H. Estelles, and R. Jaume, “Setting the cornerstone for the imrphenomx family of models for gravitational waves from compact binaries: The dominant harmonic for non-precessing quasi-circular black holes,” (2020a), arXiv:2001.11412 [gr-qc] .
  * García-Quirós  _et al._ (2020) C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, “Imrphenomxhm: A multi-mode frequency-domain model for the gravitational wave signal from non-precessing black-hole binaries,” (2020), arXiv:2001.10914 [gr-qc] .
  * Pratten  _et al._ (2020b) G. Pratten, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, M. Mateu-Lucena, R. Jaume, M. Haney, D. Keitel, J. E. Thompson, and S. Husa, “Let’s twist again: computationally efficient models for the dominant and sub-dominant harmonic modes of precessing binary black holes,” (2020b), arXiv:2004.06503 [gr-qc] .
  * Ghosh  _et al._ (2016) A. Ghosh, A. Ghosh, N. K. Johnson-McDaniel, C. K. Mishra, P. Ajith, W. Del Pozzo, D. A. Nichols, Y. Chen, A. B. Nielsen, C. P. Berry, and et al., Physical Review D 94 (2016), 10.1103/physrevd.94.021101.
  * Ghosh  _et al._ (2017) A. Ghosh, N. K. Johnson-McDaniel, A. Ghosh, C. K. Mishra, P. Ajith, W. D. Pozzo, C. P. L. Berry, A. B. Nielsen, and L. London, Classical and Quantum Gravity 35, 014002 (2017).
  * Yunes and Pretorius (2009) N. Yunes and F. Pretorius, Physical Review D 80 (2009), 10.1103/physrevd.80.122003.
  * Ramos-Buades  _et al._ (2020) A. Ramos-Buades, P. Schmidt, G. Pratten, and S. Husa, Phys. Rev. D 101, 103014 (2020), arXiv:2001.10936 [gr-qc] .
  * Marsat and Baker (2018) S. Marsat and J. G. Baker, (2018), arXiv:1806.10734 [gr-qc] .
  * Goldberg  _et al._ (1967) J. N. Goldberg, A. J. MacFarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan, J. Math. Phys. 8, 2155 (1967).
  * Wiaux  _et al._ (2007) Y. Wiaux, L. Jacques, and P. Vandergheynst, J. Comput. Phys. 226, 2359 (2007), arXiv:astro-ph/0508514 [astro-ph] .
  * Cabero  _et al._ (2017) M. Cabero, A. B. Nielsen, A. P. Lundgren, and C. D. Capano, Physical Review D 95 (2017), 10.1103/physrevd.95.064016.
  * Blanchet  _et al._ (1995) L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Phys. Rev. Lett. 74, 3515 (1995), arXiv:gr-qc/9501027 .
  * Blanchet (1998) L. Blanchet, Class. Quant. Grav. 15, 113 (1998), [Erratum: Class.Quant.Grav. 22, 3381 (2005)], arXiv:gr-qc/9710038 .
  * Blanchet  _et al._ (2004) L. Blanchet, T. Damour, G. Esposito-Farèse, and B. R. Iyer, Phys. Rev. Lett. 93, 091101 (2004).
  * Marchand  _et al._ (2016) T. Marchand, L. Blanchet, and G. Faye, Class. Quant. Grav. 33, 244003 (2016), arXiv:1607.07601 [gr-qc] .
  * Buonanno  _et al._ (2009) A. Buonanno, B. R. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, Physical Review D 80, 084043 (2009).
  * Bohé  _et al._ (2017) A. Bohé _et al._ , Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Blanchet  _et al._ (2008) L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, Classical and Quantum Gravity 25, 165003 (2008).
  * Buonanno  _et al._ (2013) A. Buonanno, G. Faye, and T. Hinderer, Physical Review D 87 (2013), 10.1103/physrevd.87.044009.
  * Arun  _et al._ (2008) K. G. Arun, A. Buonanno, G. Faye, and E. Ochsner, “Higher-order spin effects in the amplitude and phase of gravitational waveforms emitted by inspiraling compact binaries: Ready-to-use gravitational waveforms,” (2008), arXiv:0810.5336 [gr-qc] .
  * Faye  _et al._ (2012) G. Faye, S. Marsat, L. Blanchet, and B. R. Iyer, Classical and Quantum Gravity 29, 175004 (2012).
  * Baker  _et al._ (2008) J. G. Baker, W. D. Boggs, J. Centrella, B. J. Kelly, S. T. McWilliams, and J. R. van Meter, Physical Review D 78, 044046 (2008).
  * Kelly  _et al._ (2011) B. J. Kelly, J. G. Baker, W. D. Boggs, S. T. McWilliams, and J. Centrella, Physical Review D 84, 084009 (2011).
  * Damour and Nagar (2014) T. Damour and A. Nagar, Physical Review D 90, 024054 (2014).
  * Kokkotas and Schmidt (1999) K. D. Kokkotas and B. G. Schmidt, Living Reviews in Relativity 2 (1999), 10.12942/lrr-1999-2.
  * Cotesta  _et al._ (2018) R. Cotesta, A. Buonanno, A. Bohé, A. Taracchini, I. Hinder, and S. Ossokine, Physical Review D 98 (2018), 10.1103/physrevd.98.084028.
  * Berti  _et al._ (2009) E. Berti, V. Cardoso, and A. O. Starinets, Classical and Quantum Gravity 26, 163001 (2009).
  * Tichy and Marronetti (2008) W. Tichy and P. Marronetti, Physical Review D 78 (2008), 10.1103/physrevd.78.081501.
  * Rezzolla  _et al._ (2008) L. Rezzolla, P. Diener, E. N. Dorband, D. Pollney, C. Reisswig, E. Schnetter, and J. Seiler, The Astrophysical Journal 674, L29–L32 (2008).
  * Buonanno  _et al._ (2008) A. Buonanno, L. E. Kidder, and L. Lehner, Physical Review D 77 (2008), 10.1103/physrevd.77.026004.
  * Barausse  _et al._ (2012) E. Barausse, V. Morozova, and L. Rezzolla, The Astrophysical Journal 758, 63 (2012).
  * Healy  _et al._ (2014) J. Healy, C. O. Lousto, and Y. Zlochower, Physical Review D 90 (2014), 10.1103/physrevd.90.104004.
  * Hofmann  _et al._ (2016) F. Hofmann, E. Barausse, and L. Rezzolla, The Astrophysical Journal 825, L19 (2016).
  * Jiménez-Forteza  _et al._ (2017) X. Jiménez-Forteza, D. Keitel, S. Husa, M. Hannam, S. Khan, and M. Pürrer, Physical Review D 95 (2017), 10.1103/physrevd.95.064024.
  * Varma  _et al._ (2019a) V. Varma, D. Gerosa, L. C. Stein, F. Hébert, and H. Zhang, Physical Review Letters 122 (2019a), 10.1103/physrevlett.122.011101.
  * Boyle  _et al._ (2019) M. Boyle, D. Hemberger, D. A. B. Iozzo, G. Lovelace, S. Ossokine, H. P. Pfeiffer, M. A. Scheel, L. C. Stein, C. J. Woodford, A. B. Zimmerman, and et al., Classical and Quantum Gravity 36, 195006 (2019).
  * Harms  _et al._ (2014) E. Harms, S. Bernuzzi, A. Nagar, and A. Zenginoğlu, Classical and Quantum Gravity 31, 245004 (2014).
  * Reisswig and Pollney (2011) C. Reisswig and D. Pollney, Classical and Quantum Gravity 28, 195015 (2011).
  * Schmidt  _et al._ (2017) P. Schmidt, I. W. Harry, and H. P. Pfeiffer, “Numerical relativity injection infrastructure,” (2017), arXiv:1703.01076 [gr-qc] .
  * SXS Collaboration (2019a) SXS Collaboration, “SXS Gravitational Waveform Database,” https://www.black-holes.org/waveforms (2019a).
  * SXS Collaboration (2019b) SXS Collaboration, “The Spectral Einstein Code (SpEC),” https://www.black-holes.org/code/SpEC.html (2019b).
  * LIGO Scientific Collaboration (2020) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL), https://doi.org/10.7935/GT1W-FZ16 (2020).
  * Barsotti  _et al._ (2018) L. Barsotti, P. Fritschel, M. Evans, and S. Gras,  _The updated Advanced LIGO design curve_, Tech. Rep. LIGO-T1800044 (LIGO Project, 2018).
  * Healy and Lousto (2017) J. Healy and C. O. Lousto, Phys. Rev. D 95, 024037 (2017).
  * Smith  _et al._ (2020) R. J. E. Smith, G. Ashton, A. Vajpeyi, and C. Talbot, Monthly Notices of the Royal Astronomical Society 498, 4492–4502 (2020).
  * Speagle (2020) J. S. Speagle, Monthly Notices of the Royal Astronomical Society 493, 3132–3158 (2020).
  * Abbott  _et al._ (2016b) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016b), arXiv:1602.03837 [gr-qc] .
  * LIGO Scientific Collaboration, Virgo Collaboration (2019a) LIGO Scientific Collaboration, Virgo Collaboration, “Gravitational Wave Open Science Center,” https://www.gw-openscience.org (2019a).
  * LIGO Scientific Collaboration, Virgo Collaboration (2019b) LIGO Scientific Collaboration, Virgo Collaboration, “Power Spectral Densities (PSD) release for GWTC-1,”  (2019b), https://dcc.ligo.org/LIGO-P1900011/public .
  * LIGO Scientific Collaboration, Virgo Collaboration (2019c) LIGO Scientific Collaboration, Virgo Collaboration, “Calibration uncertainty envelope release for GWTC-1,”  (2019c), https://dcc.ligo.org/LIGO-P1900040/public .
  * Colleoni  _et al._ (2020) M. Colleoni, M. Mateu-Lucena, H. Estellés, C. García-Quirós, D. Keitel, G. Pratten, A. Ramos-Buades, and S. Husa, “Towards the routine use of subdominant harmonics in gravitational-wave inference: re-analysis of gw190412 with generation x waveform models,” (2020), arXiv:2010.05830 [gr-qc] .
  * London  _et al._ (2018) L. London, S. Khan, E. Fauchon-Jones, C. García, M. Hannam, S. Husa, X. Jiménez-Forteza, C. Kalaghatgi, F. Ohme, and F. Pannarale, Phys. Rev. Lett. 120, 161102 (2018), arXiv:1708.00404 [gr-qc] .
  * Cotesta  _et al._ (2020) R. Cotesta, S. Marsat, and M. Pürrer, Phys. Rev. D 101, 124040 (2020), arXiv:2003.12079 [gr-qc] .
  * Varma  _et al._ (2019b) V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E. Kidder, and H. P. Pfeiffer, Physical Review D 99 (2019b), 10.1103/physrevd.99.064045.
  * Vinciguerra  _et al._ (2017) S. Vinciguerra, J. Veitch, and I. Mandel, Class. Quant. Grav. 34, 115006 (2017), arXiv:1703.02062 [gr-qc] .
  * García-Quirós  _et al._ (2020) C. García-Quirós, S. Husa, M. Mateu-Lucena, and A. Borchers, (2020), arXiv:2001.10897 [gr-qc] .
  * Abbott  _et al._ (2020c) R. Abbott _et al._ (LIGO Scientific, Virgo), (2020c), arXiv:2010.14527 [gr-qc] .