# A phenomenological gravitational waveform model of binary black holes
incorporating horizon fluxes

Samanwaya Mukherjee  0000-0002-9055-5784<EMAIL_ADDRESS>Inter-University
Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune - 411007,
India  Sayak Datta<EMAIL_ADDRESS>Max-Planck-Institut für
Gravitationsphysik (Albert-Einstein-Institut), D-30167 Hannover, Germany
Leibniz Universität Hannover, D-30167 Hannover, Germany  Khun Sang Phukon
School of Physics and Astronomy and Institute for Gravitational Wave
Astronomy,
University of Birmingham, Edgbaston, Birmingham, B15 9TT, United Kingdom
Sukanta Bose  Department of Physics and Astronomy, Washington State
University, 1245 Webster, Pullman, Washington 99164-2814, U.S.A.

###### Abstract

Subjected to the tidal field of its companion, each component of a coalescing
binary suffers a slow change in its mass (tidal heating) and spin (tidal
torquing) during the inspiral and merger. For black holes, these changes are
associated with their absorption of energy and angular momentum fluxes. This
effect modifies the inspiral rate of the binary, and consequently, the phase
and amplitude of its gravitational waveform. Numerical relativity waveforms
contain these effects inherently, whereas analytical approximants for the
early inspiral phase have to include them manually in the energy balance
equation. In this work, we construct a frequency-domain gravitational waveform
model that incorporates the effects of tidal heating of black holes. This is
achieved by recalibrating the inspiral phase of the waveform model IMRPhenomD
to incorporate the phase corrections for tidal heating. We also include
corrections to the amplitude, but add them directly to the inspiral amplitude
model of IMRPhenomD. We show that the new model is faithful, with less than 1%
mismatch, against a set of hybrid waveforms, except for one outlier that
barely breaches this limit. The recalibrated model shows mismatches of up to
$\sim 16\%$ with IMRPhenomD for high mass ratios and spins. Amplitude
corrections become less significant for higher mass ratios, whereas the phase
corrections leave more impact – suggesting that the former is practically
irrelevant for gravitational wave data analysis in Advanced LIGO (aLIGO, Virgo
and KAGRA. Comparing with a set of 219 numerical relativity waveforms, we find
that the median of mismatches decreases by $\sim 4\%$ in aLIGO zero-detuned
high power noise curve, and by $\sim 2\%$ with a flat noise curve. This
implies a modest but notable improvement in waveform accuracy.

##  I Introduction

The launch of gravitational wave (GW) astronomy has had a stellar start, with
the detection of over 90 compact binary coalescences (CBCs) so far [1, 2, 3].
The fourth (O4) observation run of the ground-based GW detectors LIGO [4],
Virgo [5] and KAGRA [6] are expected to detect many more CBCs, which enables
us to subject General Relativity (GR) to unprecedented tests. Such tests
demand high precision CBC waveforms. In the field of GW data analysis,
waveform models serve as templates against which the real data are compared
and estimates of the source parameters are made by the matched filtering
technique [7] and Bayesian inference [8, 9, 10]. Since a sizable subset of
CBCs are high-mass binaries, which have a good fraction of their signal power
arriving from the late-inspiral and merger parts, these tests benefit from
employing inspiral-merger-ringdown (IMR) waveforms. One needs to solve the
full Einstein equations numerically to extract gravitational waveforms
predicted by GR to meet the accuracy standard imperative for describing the
merger-ringdown parts of a CBC. However, such simulations of numerical
relativity (NR) are of extreme computational cost, impeding the usage of this
scheme for creating long waveforms that span the entire frequency range of GW
detectors for intermediate to stellar-mass binaries. On the other hand, post-
Newtonian (PN) approximation (see, e.g., Ref. [11] and the references therein)
or the effective-one-body (EOB) framework [12, 13, 14] describes the inspiral
part in analytical forms, but fail to be reliably accurate in the late
inspiral to the merger regime.

Efforts to construct accurate IMR models have been based on the EOB formalism
with calibrations to NR data (SEOBNR [13, 15, 16] and TEOBResumS [17, 18, 19]
models), phenomenological waveform models built by combining PN and NR
waveforms [20, 21, 22, 23, 24, 25], and more recently, surrogate models [26,
27, 28]. The computational cost of EOB models, however, is burdened by the
need to tackle the orbital dynamics through solving a complex system of
ordinary differential equations. Phenomenological models, in comparison, are
equipped with closed-form expressions for the phase and amplitude of the
waveforms for a given set of binary parameters, and are much faster to
evaluate. The latter models contain PN-inspired phase and amplitude behaviors
augmented with higher-order terms calibrated against a set of “hybrid
waveforms”. The hybrids are constructed by stitching either PN (e.g., in
IMRPhenomC [21]) or EOB (e.g. IMRPhenomD [25, 24], IMRPhenomXAS [23] for
aligned-spin cases) waveforms with NR ones at a suitable frequency, to allow
smooth transitions of phase and amplitude from the inspiral to the merger-
ringdown regime. Utmost care needs to be borne when producing complete IMR
waveforms by combining these parts, to limit introducing waveform systematics,
which can reduce the effectiveness of tests of GR. The increasing modeling
errors and computational cost of gravitational waveforms with increasing mass
ratio and spins [29] have inspired the advent of reduced-order models (ROMs)
and surrogate models. By fitting interpolated decomposed waveform data pieces
over the binary parameter space, surrogate models can significantly accelerate
NR (e.g. [30, 28, 31]) or EOB waveforms (e.g. [32, 33, 34]), while maintaining
high accuracy within their parameter space of validity.

In phenomenological models for binary black holes (BBHs), the inspiral
waveform is usually constructed with the “point-particle” approximation, which
also acts as a baseline for creating generic waveforms with tidal effects
arising from the finite size of the components, e.g. for binary neutron stars
(BNSs) [35, 36]. However, even for black holes (BHs) there exists a finite-
size effect that needs to be consistently incorporated within a complete BBH
waveform model, namely, the absorption of energy and angular momentum by their
horizons. Due to their causal structure, BHs in GR are perfect absorbers that
behave as dissipative systems [37, 38]. A significant feature of a BH is its
horizon, which is a null surface and a “one-way membrane” that does not allow
energy to escape outward. In a binary, a BH experiences tidal interactions
from its companion, which cause changes in its mass, angular momentum, and
horizon area. This phenomenon is called tidal heating (TH) [39, 40, 41, 42,
43]. This effect is weak but important in many respects. Recent studies [44,
45] have found that the presence of TH can be measurable, especially in the
future era of GW detectors with orders of magnitude higher signal-to-noise
ratios (SNRs) compared to the current ones. Combined with the fact that TH is
much more significant for BHs than horizonless compact objects, this makes TH
a viable discriminator for horizons to identify BHs against theoretically
possible exotic compact objects (ECOs) [46], which can mimic BHs in their GW
signals [47, 48].

In NR simulations, TH arises naturally due to the presence of the BH apparent
horizons. While evolving the binary, BH singularities are expunged from the
computational domain by excision techniques or by the “moving puncture”
method. The properties of the apparent horizons – the masses, spins, and the
horizon areas – continue to get impacted by the energy and angular momentum
fluxes throughout the binary evolution up to the merger. Changes in BH masses
and spins impact the inspiral rate of the binary, leaving its imprint in the
phase of GW signals. Complete IMR waveforms created with NR data for late
inspiral and merger, then, would have incomplete information about BBH systems
if they are joined with inspiral approximants devoid of the imprints of
horizon fluxes. Ignoring this effect may also lead to ambiguities in the tests
of GR from GW observations of CBCs, where small deviations from GR predictions
are probed under the strong gravity conditions of a CBC. Lesser systematics
will facilitate a more accurate evaluation of GR’s validity in such cases.

In the EOB formalism, progress has been made to include the effects of horizon
absorption for nonspinning binaries [49, 50, 51] and for aligned spins [52].
These works use numerical evaluation of horizon fluxes using a frequency-
domain perturbative approach along stable and unstable orbits in the test-mass
limit. This test-mass knowledge is then hybridized with lower-order analytical
information that is valid for comparable masses. These calculations were
implemented in the construction of TEOBResumS [17], an NR-informed EOB
waveform model which includes the horizon-absorption effects in a PN-inspired
resummed form. Alvi [41] has calculated, under PN approximation, leading-order
expressions of mass and spin evolution of BBH systems. Alvi’s expressions are
explicitly valid for any mass ratio. A resummed version of these expressions
was used by Damour and Nagar [52] to describe them under the EOB framework.

In Phenomenological models, horizon absorption was taken into account in the
flux calculations of IMRPhenomC [21], up to 2.5PN order. However, later models
with more accurate calibrations spanning larger regions of the parameter space
(e.g. IMRPhenomD or IMRPhenomXAS) do not account for this effect explicitly in
their inspiral parts.

In this work, we construct a phenomenological BBH waveform model with a
consistent inclusion of TH effects, and study its contribution to an
improvement in waveform systematics. A complete BBH model, which includes the
phase and amplitude modifications for TH from the early inspiral to the
merger, not only provides a comprehensive foundation but also helps to lay the
groundwork for a more general waveform model with the horizon parameters,
which can act as discriminators for BHs from ECOs or NSs [44, 45].

We first describe some preliminaries and notations of gravitational waveforms
in Sec. II. Then in Sec. III.1, we build an inspiral model containing explicit
corrections due to TH at 2.5PN, 3.5PN, and 4PN orders in their phase and
amplitude. The point-particle baseline for the inspiral waveform is based on
an aligned-spin EOB approximant. In Sec. III.2, we describe the publicly
available NR waveforms of the SXS catalog [53] used for merger and ringdown
phases. In Sec. IV, we describe the construction of new hybrid waveforms by
stitching together the TH-corrected inspiral and merger-ringdown parts.
Following IMRPhenomD (henceforth referred to as PhenomD), we model the hybrids
in the frequency domain in Sec. V and calibrate the inspiral part with
augmented PN expressions, but against the new set of hybrid waveforms. Section
VI discusses the faithfulness of the new model, and in Sec. VII, we compare it
with PhenomD and NR data within the parameter range of calibration. We
conclude in Sec. VIII with discussions and future directions of this work.

Throughout the article, we will use geometric units by considering the
gravitational constant ($G$) and the speed of light in vacuum ($c$) set to
unity, except while calculating physical quantities.

##  II Waveform conventions and notation

We consider $l=|m|=2$ spherical harmonic modes of gravitational waves from
coalescing binary black holes. In the time domain, the complex GW strain can
be written as,

$h_{2,2}(\Theta,t)=A(\Theta,t)e^{-i\phi(\Theta,t)},$ (1)

where $\Theta$ is the parameter vector containing the parameters of the
binary. The intrinsic parameters are the dimensionless spin vectors
${\bf\it\chi_{1}},{\bf\it\chi_{2}}$ and the two masses $m_{1}$ and $m_{2}$.
Since we are concerned with non-precessing binaries only, the spin vectors
have nonzero components only along (or opposite to) the direction of the
orbital angular momentum. In terms of the normalized orbital angular momentum
$\hat{L}$ and the dimensionful spin vectors ${\bf\it S_{i}}$ ($i=1,2$ for the
two objects), the components of dimensionless spins can be expressed as

$\chi_{i}=\frac{{\bf\it S_{i}}\cdot\hat{L}}{m_{i}^{2}}\,.$ (2)

We define the mass ratio as $q=m_{1}/m_{2}\geqslant 1$, total mass
$M=m_{1}+m_{2}$, and the symmetric mass ratio $\eta=m_{1}m_{2}/M^{2}$.

The complex strain in Eqn. (1) can be constructed from the plus and cross
polarization states as

$h(\Theta,t)=h_{+}(\Theta,t)-ih_{\crossproduct}(\Theta,t)\,.$ (3)

The Fourier transform of the complex strain is defined as

$\tilde{h}(f)=\int_{-\infty}^{\infty}h(t)e^{-i2\pi ft}\differential t\,,$ (4)

where $f$ denotes frequency. For the non-eccentric waveforms that we use, $f$
can be written as a monotonically increasing function of time $t$ until the
merger.

In gravitational waveform modeling, the dimensionless frequency $Mf$ plays an
important role, since the total mass $M$ acts as a scaling factor. For better
readability, we define $F\equiv Mf$.

##  III Input waveforms

###  III.1 Inspiral waveforms

####  III.1.1 Effective-one-body description of a compact binary coalescence

In the effective-one-body approach to the two-body problem in GR, originally
developed by Buonanno and Damour [12], the dynamics of two compact objects of
masses $m_{1}$ and $m_{2}$ and spins ${\bf\it\chi_{1}}$ and ${\bf\it\chi_{2}}$
is mapped onto the dynamics of an effective particle of mass
$\mu=m_{1}m_{2}/(m_{1}+m_{2})$ and spin ${\bf\it\chi_{\ast}}$ moving in the
background of a deformed Kerr geometry with mass $M=m_{1}+m_{2}$ and spin
${\bf\it\chi}_{\rm Kerr}$. The spin mapping
$\\{{\bf\it\chi_{1}},{\bf\it\chi_{2}}\\}\to{\bf\it\chi_{\ast}}$ and the
deformation of the Kerr background, parametrized by the symmetric mass ratio
$\eta$, implies that the test-particle dynamics reproduces the PN-expanded
dynamics of the original two-body system. Free parameters are introduced into
the models that represent unknown, higher-order PN terms, or additional
physical effects like corrections due to eccentricity. Such free parameters
are calibrated to NR simulations. With the EOB system specified, its
conservative dynamics can be described by an EOB Hamiltonian [16, 15], while
the non-conservative dynamics is contained in a parametrized radiation-
reaction term that is inserted in the equations of motion. This term sums over
the outgoing GW modes and is calibrated to reproduce NR simulations. The
combination of these two pieces describes the binary inspiral through to the
merger, at which point a ringdown waveform is attached to the inspiral-merger
waveform. The ringdown waveform is constructed as a linear superposition of
the dominant quasinormal modes (QNMs) of the Kerr BH formed at merger [54,
52], with the amplitude and phase of each QNM mode determined by the process
of stitching the ringdown part with the inspiral-merger parts.

For the purpose of creating hybrid waveforms in our work, we choose SEOBNRv2
[16], an aligned-spin EOB approximant, as the point-particle baseline. More
recent versions of the SEOBNR family of waveforms are available [13] at
present, but we create our model based on the former one for a direct
comparison with PhenomD. We also mention here that for PhenomD, NR
calibrations of SEOBNRv2 were removed before creating hybrid waveforms,
referred to as SEOBv2. In this work, however, we use the original calibrated
version of the EOB model since in the hybrids, stitched at a frequency much
lower than the innermost stable circular orbit (ISCO), we expect little to no
NR information from the inspiral approximant to be present.

####  III.1.2 Corrections for tidal heating of black holes

The early inspiral part of a CBC can be described by the analytical formalism
of post-Newtonian (PN) expansions [11], under the approximation that the
source is slowly moving and weakly gravitating. In this formalism, the
evolution of the orbital phase $\Psi(t)$ of a compact binary is computed as a
perturbative expansion in a small parameter, typically taken to be the
characteristic velocity $v=(\pi Mf)^{1/3}$. This analytical procedure demands
$v\ll 1$, which makes it useful in the early inspiral phase of a CBC.

An electrically neutral spinning black hole in GR – the Kerr black hole (KBH)
– is stationary when it is isolated. On the other hand, when a KBH is a member
of a binary, it feels its companion’s tidal field, which acts as a non-
axisymmetric perturbation [39]. This perturbation causes changes in the mass,
spin, and horizon area of the KBH over time [41]. Since the KBH experiences
the tidal field of its orbiting companion, it absorbs (emits) energy from
(into) the orbit. The absorption part is present in non-spinning BHs as well.
Additionally, for a KBH, a spin frequency higher than the angular frequency of
the tidal field slows down due to tidal interactions, which in turn makes the
KBH lose its rotational energy. The slowing down of a rotating BH due to the
gravitational dissipation produced by exterior mass is analogous to the
slowing down of a rotating planet by viscous dissipation due to tides raised
by an exterior moon that increases its internal thermal content - a phenomenon
known as tidal heating. Due to this similarity, the energy and angular
momentum flux in BBHs is also termed tidal heating [40].

During the early stages of the binary evolution, the change in the orbital
frequency due to emission of GWs is much smaller than the orbital frequency
itself. In this adiabatic inspiral regime, loss of binding energy $E(v)$ of
the two-body system with time equals the GW flux emitted to future null
infinity ($\mathcal{F}^{\infty}(v)$) plus the energy flux absorbed by the two
BH horizons ($\mathcal{F}^{\rm H}(v)$). So the energy balance condition
becomes

$-\derivative{E(v)}{t}=\mathcal{F}^{\infty}(v)+\mathcal{F}^{\rm H}(v).$ (5)

Evolution of the orbital phase $\phi$ and the characteristic velocity $v$,
obtained from this equation, read

$\derivative{\phi}{t}=\frac{v^{3}}{M}\,,\quad\quad\quad\derivative{v}{t}=-\frac{\mathcal{F}(v)}{E^{\prime}(v)},$
(6)

where $\mathcal{F}(v)=\mathcal{F}^{\infty}(v)+\mathcal{F}^{\rm H}(v)$. These
equations yield a solution for the phase $\Psi(f)$ of the frequency-domain
waveform $\tilde{h}(f)=\tilde{A}(f)e^{-i\Psi(f)}$ [55]:

$\Psi(f)=2(t_{c}/M)v^{3}-2\phi_{c}-\pi/4-\frac{2}{M}\int(v^{3}-\bar{v}^{3})\frac{E^{\prime}(\bar{v})}{\mathcal{F}(\bar{v})}\differential\bar{v},$
(7)

where $E^{\prime}(v)=\derivative*{E(v)}{v}$.

The frequency-domain amplitude $\tilde{A}(f)$ can be written as111Note that we
define $E(v)$ as the total binding energy. Ref [56] treats $E(v)$ as the
specific binding energy (binding energy per unit mass), so the power of $M$ in
our treatment is different from theirs in this particular equation. [56]

$\tilde{A}(f)=\mathcal{C}\frac{2\eta
M^{3/2}}{d_{L}}v\sqrt{\frac{\pi}{3}}\left[-\frac{E^{\prime}(v)}{\mathcal{F}(v)}\right]^{1/2},$
(8)

where $\mathcal{C}$ is a constant that depends on the inclination of the
orbital plane with the line of sight, and $d_{L}$ denotes the distance to the
binary.

Splitting $\mathcal{F}$ into $\mathcal{F}^{\infty}$ and $\mathcal{F}^{\rm H}$
in Eq. (7) enables one to write the frequency-domain phase in the form

$\Psi(f)=\Psi_{\rm PP}(f)+\Psi_{\rm TH}(f),$ (9)

with $\Psi_{\rm PP}(f)$ being the GW phase under point particle (PP)
approximation that corresponds to $\mathcal{F}^{\infty}$, and $\Psi_{\rm
TH}(f)$ being the phase correction due to TH. The absorbed flux can be
expressed as the sum of the rates of mass increment [41] for the two BHs:

$\displaystyle\mathcal{F}^{\rm H}(v)=$
$\displaystyle{}\frac{32}{5}\eta^{2}\frac{v^{15}}{4}\sum_{i=1}^{2}\left(\frac{m_{i}}{M}\right)^{3}\left(1+3\chi^{2}_{i}\right)$
(10)
$\displaystyle\times\left\\{-(\hat{L}\cdot\hat{S}_{i})\chi_{i}+2\left[1+\sqrt{1-\chi^{2}_{i}}\right]\frac{m_{i}}{M}v^{3}\right\\}\
,$

where the indices $i=1,2$ denote the two BHs, and $\hat{S}_{i}$ and $\hat{L}$
are the unit vectors along the directions of the $i$th object’s spin and the
orbital angular momentum, respectively.

Using this expression, one can calculate the phase contribution of tidal
heating in gravitational waveforms of BBHs to be [44, 45]

$\displaystyle\Psi_{\rm TH}=$
$\displaystyle{}\frac{3}{128\eta}\left[-\frac{10}{3}H^{(\rm BBH)}_{\rm
eff5}\,\log\left(v\right)\right.$ (11)
$\displaystyle-\frac{5}{168}v^{2}H^{(\rm BBH)}_{\rm
eff5}\left(952\eta+995\right)$
$\displaystyle\left.+\frac{5}{3}v^{3}\log\left(v\right)(-4H^{(\rm BBH)}_{\rm
eff8}+H^{(\rm BBH)}_{\rm eff5}\Psi_{\text{SO}})\right],$

where the quantities $H^{(\rm BBH)}_{\rm eff5}$ and $H^{(\rm BBH)}_{\rm eff8}$
are defined as

$\displaystyle H^{(\rm BBH)}_{\rm eff5}\equiv$
$\displaystyle{}\sum_{i=1}^{2}\left(\frac{m_{i}}{M}\right)^{3}\left(\hat{L}\cdot\hat{S}_{i}\right)\chi_{i}\left(3\chi_{i}{}^{2}+1\right)\,,$
(12a) $\displaystyle H^{(\rm BBH)}_{\rm eff8}\equiv$
$\displaystyle{}\leavevmode\nobreak\ 4\pi H^{(\rm BBH)}_{\rm
eff5}+\sum^{2}_{i=1}\left(\frac{m_{i}}{M}\right)^{4}\left(3\chi_{i}{}^{2}+1\right)$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\times\left(\sqrt{1-\chi_{i}{}^{2}}+1\right)\,.$
(12b)

$\Psi_{\text{SO}}$ is a “spin-orbit” term, given by

$\begin{split}\Psi_{\text{SO}}=&\leavevmode\nobreak\
\frac{73}{3(1+q)^{2}}\left\\{q^{2}\big{(}\hat{L}\cdot\hat{S}_{1}\big{)}\chi_{1}+\big{(}\hat{L}\cdot\hat{S}_{2}\big{)}\chi_{2}\right\\}\\\
&+\frac{15q}{(1+q)^{2}}\left\\{\big{(}\hat{L}\cdot\hat{S}_{1}\big{)}\chi_{1}+\big{(}\hat{L}\cdot\hat{S}_{2}\big{)}\chi_{2}\right\\}\,.\end{split}$
(13)

We find the PN expansion of the total amplitude by expanding
$\left[-E^{\prime}(v)/\mathcal{F}(v)\right]^{1/2}$ in Eq. (8) in powers of $v$
up to $v^{7}$, and separate the contribution due to TH:

$\tilde{A}(f)=\tilde{A}_{\rm PP}(f)+\tilde{A}_{\rm TH}(f).$ (14)

$\tilde{A}_{\rm PP}(f)$ is mentioned in Eq. (5.7) of Ref. [56], and we get the
expression for $\tilde{A}_{\rm TH}(f)$ as

$\displaystyle\tilde{A}_{\rm TH}(f)=\leavevmode\nobreak\ \mathcal{C}$
$\displaystyle\frac{M^{5/6}}{d_{L}\pi^{2/3}}\sqrt{\frac{5\eta}{24}}f^{-7/6}\left[\frac{1}{8}H^{(\rm
BBH)}_{\rm eff5}v^{5}\right.$ (15) $\displaystyle\left.\leavevmode\nobreak\
+\left(\frac{1079}{1792}+\frac{103}{192}\eta\right)H^{(\rm BBH)}_{\rm
eff5}v^{7}\right.$ $\displaystyle\left.+\frac{1}{32}\left\\{\tilde{A}_{\rm
SO}H^{(\rm BBH)}_{\rm eff5}\right.\right.$
$\displaystyle\left.\left.-8\left(H^{(\rm BBH)}_{\rm eff8}-\pi H^{(\rm
BBH)}_{\rm eff5}\right)\right\\}v^{8}\right],$

where $\tilde{A}_{\rm SO}$ is the spin-orbit term

$\displaystyle\tilde{A}_{\rm SO}=$ $\displaystyle\leavevmode\nobreak\
\frac{179}{6(1+q)^{2}}\left\\{q^{2}(\hat{L}\cdot\hat{S}_{1})\chi_{1}+(\hat{L}\cdot\hat{S}_{2})\chi_{2}\right\\}$
(16)
$\displaystyle+\frac{35q}{2(1+q)^{2}}\left\\{(\hat{L}\cdot\hat{S}_{1})\chi_{1}+(\hat{L}\cdot\hat{S}_{2})\chi_{2}\right\\}.$

###  III.2 Numerical relativity waveforms for merger-ringdown

As a BBH progresses towards the merger phase, the BHs come so close to each
other that the system becomes too compact for the weak-gravity condition to
hold. In addition, their orbital velocities culminate to values comparable to
$c$, breaching the slow-motion approximation. The analytical PN results, with
their highest order of expansion available currently, become increasingly
inaccurate to describe the system evolving under such extreme conditions.
Starting from this late inspiral phase to the merger and ringdown phase where
the two BHs merge and the final remnant BH settles to equilibrium, one needs
to solve the full Einstein equations numerically for extracting GWs to conform
to the accuracy standards of the waveforms required for searches of GW signals
or their parameter estimation (PE). This is the field of NR, which has seen
exciting breakthroughs [57, 58, 59] in the early 2000s, and great strides
towards simulating merging black holes with different masses and spins [60,
61, 62, 63].

We use NR data from the publicly available SXS catalog [53], computed using
the Spectral Einstein Code (SpEC) [64, 65]. SpEC uses the excision techniques
to remove the BH singularity for extracting gravitational waves from binary
black hole systems. Unlike the PN framework where the fluxes of energy and
angular momentum absorbed by the BHs have to be considered explicitly in the
energy balance equation, numerical simulations of full Einstein equations
capture this effect inherently. Scheel et al. [66] have demonstrated the slow
change in the BH masses and spins through the binaries’ temporal evolution in
NR simulations with SpEC. They track the apparent horizons as a function of
time, and at frequent time intervals, they measure both the surface area and
the spin of the horizons. The spin computation has been carried out using the
approximate Killing vector formalism [67]. The mass of the black hole is then
computed using Christodoulou’s formula [68],

$M^{2}=M_{\rm irr}^{2}+\frac{S^{2}}{4M_{\rm irr}^{2}},$ (17)

where $M_{\rm irr}$ is the irreducible mass of the Kerr BH,

$M_{\rm irr}^{2}=\frac{1}{16\pi}\int_{\mathcal{H}}\differential A.$ (18)

Here $S$ is the (dimensionful) spin of the KBH and the integration covers the
total horizon area $A$.

# |  Simulation label |  $q$ |  $\chi_{1}$ |  $\chi_{2}$ |  $\chi_{\rm PN}$ |  $N_{\rm orb}$ |  $F_{\rm stitch}$
---|---|---|---|---|---|---|---
1 |  SXS:BBH:0156 |  1 |  $-$0.95 |  $-$0.95 |  $-$0.79 |  13 |  0.0057
2 |  SXS:BBH:0151 |  1 |  $-$0.6 |  $-$0.6 |  $-$0.5 |  15 |  0.0054
3 |  SXS:BBH:0001 |  1 |  0 |  0 |  0 |  28 |  0.0041
4 |  SXS:BBH:0152 |  1 |  0.6 |  0.6 |  0.5 |  23 |  0.0052
5 |  SXS:BBH:0153 |  1 |  0.85 |  0.85 |  0.7 |  25 |  0.0052
6 |  SXS:BBH:0234 |  2 |  $-$0.85 |  $-$0.85 |  $-$0.72 |  28 |  0.0037
7 |  SXS:BBH:0238 |  2 |  $-$0.5 |  $-$0.5 |  $-$0.43 |  32 |  0.0037
8 |  SXS:BBH:0169 |  2 |  0 |  0 |  0 |  16 |  0.0062
9 |  SXS:BBH:0253 |  2 |  0.5 |  0.5 |  0.43 |  29 |  0.0046
10 |  SXS:BBH:2131 |  2 |  0.85 |  0.85 |  0.72 |  25 |  0.0056
11 |  SXS:BBH:1936 |  4 |  $-$0.8 |  $-$0.8 |  $-$0.71 |  17 |  0.0057
12 |  SXS:BBH:1418 |  4 |  $-$0.4 |  $-$0.5 |  $-$0.37 |  67 |  0.005
13 |  SXS:BBH:0167 |  4 |  0 |  0 |  0 |  16 |  0.0071
14 |  SXS:BBH:1417 |  4 |  0.4 |  0.5 |  0.37 |  80 |  0.0057
15 |  SXS:BBH:1907 |  4 |  0 |  0.8 |  0.12 |  21 |  0.0062
16 |  SXS:BBH:1423 |  8 |  $-$0.6 |  $-$0.75 |  $-$0.57 |  18 |  0.0066
17 |  SXS:BBH:0064 |  8 |  $-$0.5 |  0 |  $-$0.43 |  19 |  0.0066
18 |  SXS:BBH:0063 |  8 |  0 |  0 |  0 |  26 |  0.0065
19 |  SXS:BBH:0065 |  8 |  0.5 |  0 |  0.43 |  34 |  0.0063
20 |  SXS:BBH:1426 |  8 |  0.48 |  0.75 |  0.47 |  26 |  0.0078
Table 1: Hybrid waveforms used to calibrate the model. The first column lists
the simulation IDs of the corresponding NR data in the SXS catalog. The last
column reports the dimensionless frequency at the midpoint of the stitching
region.

Figure 1: Hybrid waveforms for four different configurations in the parameter
space of $q,\chi_{1},\chi_{2}.$ Inspiral waveforms are generated by adding the
phase and amplitude corrections due to TH to the SEOBNRv2 model, shown in
blue. NR waveforms are shown in orange, and the hybrid waveforms are shown in
black dashed-dotted lines. The $x$ and $y$ axes denote time (in units of total
mass) and the real part of time-domain strain, $h_{+}(t)$ (in units of
$M/D_{L}$), respectively. The blue shaded areas denote the stitching regions.

During the evolution, each SXS NR waveform is extracted at a series of times
on a set of concentric coordinate spheres surrounding the binary, decomposed
in modes of spin-weighted spherical harmonic functions. Then the waveforms are
extrapolated to future null infinity $\mathscr{I}^{+}$. The dataset of each
SXS waveform contains several gravitational waveform modes $(l,m)$, and the
orders of extrapolation $N=2,3,4$. A higher order of extrapolation is
preferred when accurate waveforms are needed in the inspiral, and lower order
extrapolation is preferred for more accuracy at the merger-ringdown phase
[53]. Since we use these waveforms for constructing hybrids, we need more
accuracy for the merger-ringdown part, due to the fact that some part of the
inspiral has to be replaced by the PN waveforms in the hybrids. We use 20
aligned-spin non-eccentric BBH waveforms from SXS with $1\leqslant q\leqslant
8$ and $-0.95\leqslant\chi_{1},\chi_{2}\leqslant 0.85$ for calibrating our
model. We choose $(2,2)$ modes, and the extrapolation order $N=2$ for our
purpose.

##  IV Hybrid construction

For constructing a hybrid, the NR waveform has to be ‘stitched’ with the
analytical inspiral waveform in a frequency region where both the waveforms
have sufficient accuracy, which also demands that the morphology of these
waveforms do not deviate significantly from one another. For PN waveforms, one
conventionally chooses ISCO of the binary to be the endmost point of validity.
The corresponding GW frequency $f_{\rm ISCO}$ is considered to be the highest
frequency for PN expansions. In our work, we ensured that the frequency at the
stitching region is below this threshold for all the hybrids - which also
ensures that the NR calibrations in SEOBNRv2 are minimally present before the
stitching starts.

While stitching at a frequency as low as possible (considering ISCO) should
enable one to incorporate the maximum number of NR cycles into the hybrid,
there are two more factors that affect the choice of the stitching region:

  * •

The initial NR data do not perfectly describe two black holes in quasi-
equilibrium. At the start of each simulation, the geometry relaxes to
equilibrium on the dynamical time-scale of the individual BHs, changing the
mass and spin of each BH by a fractional amount of order $10^{-5}$, and
emitting a spurious pulse of gravitational radiation (often referred to as
‘junk radiation’).

  * •

A robust hybrid should depend weakly on the small changes to the stitching
region. A monochromatic signal would be completely degenerate under a shift in
the coalescence time $t_{c}$ and the coalescence phase $\phi_{c}$; this
degeneracy is broken by the increase in frequency with time. The start and end
of the stitching interval should reflect enough change in the frequency to
break this degeneracy.

The first issue is addressed in our work by choosing the stitching region
beyond the specified relaxation time for each NR simulation in the SXS
catalog, the time interval (in units of $M$) by which the junk radiation dies
out. To address the second requirement, MacDonald et al. [69] have recommended
the stitching interval to satisfy $\delta\omega/\omega_{m}\geqslant 0.1$,
where $\omega_{m}$ is the GW frequency at the midpoint of the stitching
interval, and $\delta\omega$ is the change in frequency over the interval.
This choice ensures that the residual oscillations in $t_{c}$ with changing
$\omega_{m}$ are below $1M$. In our work, we perform the stitching over 4 GW
cycles, and we place the stitching region in a way to satisfy this condition.
In the literature, the construction of hybrids has been performed both in time
[20, 24, 23] and frequency domains [21, 70]. In our work, we choose the
former, primarily because the NR waveforms are available in the time domain,
and typically they have too few numbers of GW cycles for a reliable Fourier
transform.

In time domain, any two non-precessing binary waveforms $h^{A}(\Theta,t)$ and
$h^{B}(\Theta,t)$, with the same set of intrinsic parameters $\Theta$, only
differ by a relative time shift ($\Delta t$) and an overall constant phase
difference ($\phi_{0}$):

$h^{A}(\Theta,t)=e^{i\phi_{0}}h^{B}(\Theta,t+\Delta t).$

In the region where we expect the PN and NR waveforms to be morphologically
similar, we can write

$h^{\rm PN}(\Theta,t)=e^{i\phi_{0}}h^{\rm NR}(\Theta,t+\Delta t).$ (19)

However, in reality, $h^{\rm PN}(\Theta,t)$ and $h^{\rm NR}(\Theta,t)$ differ
from each other morphologically due to waveform systematics, so Eq. (19) does
not hold in general [69]. To generate a hybrid waveform, then, one needs to
find a suitable region where PN and NR waveforms are approximately the same
and minimize the square of the difference between the left-hand and right-hand
sides in Eq. (19), by varying $\Delta t$ and $\phi_{0}$. Since the GW
frequency of a CBC increases with time, one can alternatively write the
angular GW frequency as a function of time ($\omega(t)$) and minimize the
quantity

$\delta=\int_{t_{1}}^{t_{2}}\left(\omega^{(\rm NR)}(t)-\omega^{(\rm
PN)}(t+\Delta t)\right)^{2}\differential t$ (20)

with respect to $\Delta t$. Here $t_{1}$ and $t_{2}$ denote the start and end
of the stitching region.

We list the steps taken to construct the hybrids as follows:

  1. 1.

SEOBNRv2 is a time-domain waveform, and the phase and amplitude corrections
for TH are calculated in the frequency domain. So, we have to first convert
the time-domain data to the Fourier domain to apply the corrections. Since the
time-domain data is finite, we use a tapered window function to minimize Gibbs
artifacts, called the Planck-taper window [69]:

$\displaystyle w(x)=\left\\{\begin{array}[]{ll}0,&x\leqslant x_{1}\\\
\left[e^{y(x)}+1\right]^{-1},&x_{1}<x<x_{2}\\\ 1,&x_{2}\leqslant x\leqslant
x_{3}\\\ \left[e^{z(x)}+1\right]^{-1},&x_{3}<x<x_{4}\\\ 0,&x_{4}\leqslant
x,\end{array}\right.$ (21)

where $y(x)=(x_{2}-x_{1})/(x-x_{1})+(x_{2}-x_{1})/(x-x_{2})$, and
$z(x)=(x_{3}-x_{4})/(x-x_{3})+(x_{3}-x_{4})/(x-x_{4})$.

The frequency-domain waveform, after correcting for TH, is then converted back
to the time domain for hybridization.

  2. 2.

NR data from the SXS catalog sample the GW strain $h^{(\rm NR)}(t)$ in non-
uniform timesteps, with a higher sampling rate in regions of higher GW
frequency. We resample the data with the coarsest sampling rate by using cubic
interpolation.

  3. 3.

Having obtained $h^{(\rm PN)}(t)$ and $h^{(\rm NR)}(t)$, we minimize the
quantity $\delta$ defined in Eq. (20) with respect to $\Delta t$. Phase
alignments are done at the temporal midpoint of the stitching region,
$t=(t_{1}+t_{2})/2$.

  4. 4.

The hybrid waveform is then constructed as

$\displaystyle h^{(\rm Hyb)}(t)=$ $\displaystyle\leavevmode\nobreak\
\mathcal{F}(t)h^{(\rm PN)}(t)$ (22)
$\displaystyle+\quantity[1-\mathcal{F}(t)]e^{i\phi^{\prime}_{0}}h^{(\rm
NR)}(t+\Delta t^{\prime}),$

where $\phi^{\prime}_{0}$ and $\Delta t^{\prime}$ are the phase offset at the
midpoint of the stitching region and the value of $\Delta t$ that minimizes
$\delta$, respectively. $\mathcal{F}(t)$ is a blending function defined by

$\displaystyle\mathcal{F}(t)=\left\\{\begin{array}[]{ll}1,&t<t_{1}\\\
\frac{1}{2}\quantity(1+\cos\frac{\pi(t-t_{1})}{t_{2}-t_{1}}),&t_{1}\leqslant
t\leqslant t_{2}\\\ 0,&t_{2}<t\,.\end{array}\right.$ (23)

In the time domain hybrids, parts of the waveforms at $t<t_{1}$ are purely
post-Newtonain, and for $t>t_{2}$ the waveforms contain purely NR data. In
between, the blending function $\mathcal{F}(t)$ smoothly stitches the two
waveforms together after proper alignment in phase and time.

Figure 1 shows four of the total twenty hybrids made for calibrating the
waveform model, at the extreme points of the calibration parameter space. The
blue shaded regions correspond to the stitching intervals, $t_{1}\leqslant
t\leqslant t_{2}$. In Table 1, we list the NR data used to create the set of
hybrids for calibration. We report the simulation labels of aligned-spin NR
data from SXS, dimensionless spin magnitudes for the two component BHs, the
effective spin parameter used for modeling the inspiral (discussed in Sec.
V.3), the number of orbits present in the NR data, and the value of
dimensionless frequency $Mf$ at the midpoint of the stitching regions.

##  V Modeling the hybrids in the frequency domain

We perform phenomenological parametrization of the hybrid waveforms in the
frequency domain. Template waveforms in the frequency domain are of particular
preference since the parametrization is easier, and also a search using
frequency-domain templates is computationally inexpensive compared to one
using time-domain templates. Our parametrization follows the construction of
PhenomD [24, 25] for modeling the phase and amplitude of the frequency-domain
hybrids.

PhenomD used uncalibrated SEOBNRv2 as the inspiral approximant, while 19 NR
waveforms extracted by SpEC and BAM [71] codes were used for hybrid
construction and calibration. The model was constructed in the frequency
domain by dividing the entire waveform into three parts – inspiral,
intermediate, and merger-ringdown. This modular nature allows one to use a
different inspiral model while keeping the merger and ringdown intact. In that
model, inspiral is defined to be the region below $F=0.018$. Merger and
ringdown phases are separated by identifying the ringdown frequencies of
different hybrids that are used to calibrate the model. In PhenomD, the end of
inspiral ($F=0.018$) was chosen in such a way that beyond that frequency, each
of the hybrids contained purely NR data. The merger-ringdown model of PhenomD,
then, is a model of purely NR data, which we do not need to recalibrate in
this work. In fact, as discussed in Ref. [24], the modular nature of PhenomD
allows independent modeling of the inspiral, intermediate, and merger-ringdown
parts separately, without even needing any hybrid waveform. The hybrid
dataset, however, serves as a benchmark to compare the final model with and
provide long time-domain waveforms for better frequency-domain modeling.
Additionally, the hybrids are created with a smooth transition of phase and
amplitude by the stitching process, producing accurate values of the pseudo-PN
parameters capturing higher-order terms in frequency.

###  V.1 Inspiral phase model

Figure 2: Fits of $\sigma_{i}$ calculated by the NonLinearModelFit module of
Mathematica. The surfaces correspond to the 2D fits with $\eta$ and $\chi_{\rm
PN}$ as described in Eq. (30).

The inspiral approximant used in the hybrid data is SEOBNRv2 with corrections
in phase and amplitude due to tidal heating, as discussed earlier. Since the
final model is a phenomenological one, it requires ready-to-use analytical
forms for the phase and amplitude with extra coefficients calibrated to the
hybrid data. We write the frequency-domain phase of a waveform in the inspiral
as

$\displaystyle\Psi_{\rm\tiny{INS}}(F;\Theta)=$
$\displaystyle\leavevmode\nobreak\ \Psi_{\rm\tiny{TF2}}(F;\Theta)+\Psi_{\rm
TH}(F;\Theta)$ (24)
$\displaystyle+\frac{1}{\eta}\left[\sigma_{0}+\sigma_{1}F+\frac{3}{4}\sigma_{2}F^{4/3}\right.$
$\displaystyle\left.\quad\quad+\frac{3}{5}\sigma_{3}F^{5/3}+\frac{1}{2}\sigma_{4}F^{2}\right],$

where $\Psi_{\rm\tiny{TF2}}(F;\Theta)$ and $\Psi_{\rm TH}(F;\Theta)$ are 3.5PN
phase of TaylorF2 (described in Appendix A) and the TH correction (Eq. (11)),
respectively. $\sigma_{i}\>\>(i=0-4)$ are 5 phenomenological pseudo-PN
parameters. The analytical ansatz in Eq. (24) is the same as in PhenomD
inspiral model. Since we want to leverage the modular nature of PhenomD to
concatenate its merger-ringdown model with our own inspiral model, the end of
the inspiral phase in the current work needs to remain unchanged, which is
$F=0.018$. We also note here that one can, in principle, directly add the
correction $\Psi_{\rm TH}(F;\Theta)$ to the PhenomD inspiral phase (as is done
for the amplitude, discussed in Sec. V.2), without any recalibration. However,
since we intend to avoid adding the correction to the NR information present
in the last few cycles of inspiral, recalibrating the phase allows a smooth
transition from PN to NR. To find the parameter values corresponding to each
hybrid waveform, we fit the first derivative of the phase with respect to $F$:

$\displaystyle\Psi^{\prime}_{\rm\tiny INS}(F;\Theta)=$
$\displaystyle\leavevmode\nobreak\
\Psi^{\prime}_{\rm\tiny{TF2}}(F;\Theta)+\Psi^{\prime}_{\rm TH}(F;\Theta)$ (25)
$\displaystyle+\frac{1}{\eta}\left[\sigma_{1}+\sigma_{2}F^{1/3}\right.$
$\displaystyle\left.\quad\quad+\sigma_{3}F^{2/3}+\sigma_{4}F\right],$

where $\Psi^{\prime}=\partial\Psi/\partial F$. $\sigma_{0}$ is determined by
imposing $C^{(1)}$ continuity in phase at the boundary between inspiral and
merger. The fits are performed over the frequency range $F\in[0.0035,0.019]$,
ending at a slightly higher frequency, to reduce boundary effects at the
interface and find robust fits for $\sigma_{i}$. For all the 20 hybrids used
for calibration, the stitching regions are placed within this frequency range,
which warrants that there is purely NR data beyond $F=0.018$.

###  V.2 Inspiral and intermediate amplitude model

To model the amplitude within the frequency range of the inspiral, we build
the inspiral amplitude model by adding the correction due to TH, given by Eq.
(15), to the PhenomD amplitude:

$\tilde{A}_{\rm INS}(F)=\tilde{A}_{\rm D}(F)+\tilde{A}_{\rm TH}(F)\,.$ (26)

We do not recalibrate the amplitude pseudo-PN parameters of PhenomD, since the
inspiral amplitude model ends (at $F=0.014$) before the phase model does (at
$F=0.018$), and the amount of NR information is less than that of phase.
Moreover, we find that the effect of the amplitude correction is significantly
less than the phase correction, diminishing the necessity of a recalibration.

At the interface of two frequency intervals modeled separately, $C^{(1)}$
continuity on the amplitude cannot be imposed in a straightforward manner. For
the frequency-domain phase, one has the freedom to tune the coalescence time
$t_{c}$ and the coalescence phase $\phi_{c}$ to impose continuity of the phase
and its derivative across an interface. For amplitude, however, there is no
such freedom. If the amplitude is modeled separately in different frequency
regimes, a different strategy has to be followed. To ensure that the amplitude
and its derivative are continuous throughout the three phases, one needs a
frequency interval to fit a polynomial that satisfies these conditions at the
beginning and the end. This ‘intermediate’ interval is defined to be
$F\in[0.014,F_{\rm peak}]$, where $F_{\rm peak}\equiv Mf_{\rm peak}$ is the
frequency corresponding to the peak amplitude. For $F<0.014$, the amplitude is
given by Eq. (26). In the intermediate region, the amplitude is approximated
as a polynomial in $f$:

$\displaystyle\tilde{A}_{\rm int}=A_{0}$
$\displaystyle\left[\delta_{0}+\delta_{1}f+\delta_{2}f^{2}\right.$ (27)
$\displaystyle\left.+\delta_{3}f^{3}+\delta_{4}f^{4}\right].$

Here $A_{0}$ includes the leading order $f^{-7/6}$ behaviour. Evaluation of
$\delta_{i}$ follows the steps of PhenomD, described in Sec. V(C) of Ref [25].
We briefly summarize it here:

  * •

Equation (27) has five parameters, requiring five independent equations for
unique solutions. Two of them come from the $C^{(1)}$ continuity of
$\tilde{A}_{\rm int}$ with $\tilde{A}_{\rm INS}(f)$ at the beginning
($F=0.014$), and two from $C^{(1)}$ continuity with PhenomD amplitude at
$F_{\rm peak}$. The fifth equation appears from an additional condition that
the polynomial coincides with PhenomD amplitude at the mid-frequency (chosen
as a collocation point), $F=(0.014+F_{\rm peak})/2$.

  * •

Solving the aforementioned set of equations, one finds the form of
$\tilde{A}_{\rm int}$ which smoothly connects the intermediate region with
inspiral and merger.

  * •

It is worth mentioning that while choosing more collocation points would
improve the model accuracy, it would also encumber the model with a larger set
of equations to solve. Keeping in mind that amplitude errors are less
consequential than phase errors, we refrain from considering more collocation
points.

In Appendix B, we describe the evaluation of the peak frequency $f_{\rm peak}$
and the ringdown frequency $f_{\rm RD}$ as functions of the BH masses and
spins.

###  V.3 Correspondence between the phenomenological and physical parameters

Figure 3: Phase, phase derivative, and phase difference between the hybrids
(same as in Fig. 1), individual fits obtained from Eq. (25), and the final
model prediction. Top panels: Inspiral phase $\Psi(Mf)$. Middle panels:
$-\Psi^{\prime}(Mf)=-\partial\Psi/\partial(Mf)$ plotted against $Mf$. Bottom
panels: $\Delta\Psi_{\rm Fit}=(\Psi_{\rm Hybrid}-\Psi_{\rm Fit})$, and
$\Delta\Psi_{\rm Model}=(\Psi_{\rm Hybrid}-\Psi_{\rm Model})$. All the plots
are shown within the inspiral frequency range, $0.0035\leqslant Mf\leqslant
0.018$.

Figure 4: Rescaled amplitude of the full IMR waveforms as functions of $Mf$.
The blue shaded regions denote the intermediate frequency range,
$0.014\leqslant Mf\leqslant F_{\rm peak}$.

Our model has 4 phenomenological parameters corresponding to the inspiral
phase. These parameters play an intermediate role in generating gravitational
waveforms of aligned-spin compact binaries. To generate waveforms for
arbitrary values of the masses and spins, one needs to build a correspondence
between the set of physical parameters $\\{\eta,\chi_{1},\chi_{2}\\}$ and the
set of phenomenological parameters $\\{\sigma_{i}\\}$. The total mass $M$ of
the binary works as a trivial scaling factor.

In PN expansions, the leading order contribution of spin in the phase of the
waveform appears as a function of the combination

$\chi_{\rm PN}=\chi_{\rm eff}-\frac{38\eta}{113}(\chi_{1}+\chi_{2}),$ (28)

where

$\chi_{\rm eff}=\frac{m_{1}\chi_{1}+m_{2}\chi_{2}}{M}.$ (29)

In this model, $\chi_{\rm PN}$ is used as a single spin parameter to generate
the phenomenological parameters. This “effective spin approximation” works
well for aligned-spin binaries, especially in the inspiral regime [25].

To establish the correspondence between physical and phenomenological
parameters, we fit $\sigma_{i}$ obtained from different hybrids with a
polynomial of $\eta$ (up to second order) and the effective spin parameter
$\chi_{\rm PN}$ (up to third order):

$\displaystyle\sigma_{i}=\leavevmode\nobreak\
\lambda_{00}^{i}+\lambda_{10}^{i}\eta\leavevmode\nobreak$
$\displaystyle+(\chi_{\rm
PN}-1)(\lambda_{01}^{i}+\lambda_{11}^{i}\eta+\lambda_{21}^{i}\eta^{2})$ (30)
$\displaystyle+(\chi_{\rm
PN}-1)^{2}(\lambda_{02}^{i}+\lambda_{12}^{i}\eta+\lambda_{22}^{i}\eta^{2})$
$\displaystyle+(\chi_{\rm
PN}-1)^{3}(\lambda_{03}^{i}+\lambda_{13}^{i}\eta+\lambda_{23}^{i}\eta^{2}).$

For each $\sigma_{i}$, this polynomial generates a set
$\\{\lambda^{i}_{jk}\\}$ to build the correspondence between physical and
phenomenological parameter sets. Figure 2 shows the surfaces of $\sigma_{i}$
as functions of $\\{\eta,\chi_{\rm PN}\\}$ with the ansatz in Eq. (30). These
best-fit surfaces generate a constant set of $\\{\lambda^{i}_{jk}\\}$ which
define the inspiral model, along with the TaylorF2 phase and the TH
correction. In Appendix C, we list the coefficient values corresponding to
each $\sigma_{i}$. A schematic flowchart for generating the inspiral phase
model starting from the intrinsic binary parameters can be described as:

${\\{\eta,\chi_{1},\chi_{2}\\}}\xrightarrow[\text{spin}]{\text{Effective}}\\{\eta,\chi_{\rm
PN}\\}\xrightarrow{\lambda^{i}}{\\{\sigma_{i}\\}}\xrightarrow{}{\text{Inspiral
model}}.$

###  V.4 Full IMR model

For the phase, we have separately modeled the frequency region
$0.0035\leqslant F\leqslant 0.018$, and we have the phase model of PhenomD
thereafter. For the amplitude, we have a corrected model in the range
$F\leqslant 0.014$, and a polynomial function for $0.014\leqslant F\leqslant
F_{\rm peak}$. This gives us one interface for the phase and two interfaces
for the amplitude to impose $C^{(1)}$ continuity. For the amplitude, however,
$C^{(1)}$ continuity is ensured by the intermediate polynomial itself. For the
phase, we can vary $t_{c}$ and $\phi_{c}$ of the inspiral phase. After
imposing $C^{(1)}$ continuity, we can generate the complete IMR phase and
amplitude by defining a step function

$\displaystyle\theta(f-f_{0})=\left\\{\begin{array}[]{ll}-1,&f<f_{0}\\\
+1,&f\geqslant f_{0}\,.\end{array}\right.$ (31)

Using this function, one can define

$\theta^{\pm}(f;f_{0})=\frac{1}{2}\quantity[1\pm\theta(f-f_{0})],$ (32)

so that the IMR phase can be written as

$\Phi_{\rm IMR}(F)=\Phi_{\rm
INS}(F)\theta^{-}(F;F_{2})+\theta^{+}(F;F_{2})\Phi_{\rm D}(F),$ (33)

where $F_{2}=0.018$, and $\Phi_{\rm D}(F)$ corresponds to the phase model of
PhenomD. $\Phi_{\rm INS}(F)$ is given by Eq. (24).

The IMR amplitude model follows a similar treatment, which can be expressed as

$\displaystyle\tilde{A}_{\rm IMR}(F)$ $\displaystyle=\tilde{A}_{\rm
INS}(F)\theta^{-}(F;F_{3})$ (34) $\displaystyle+$
$\displaystyle\theta^{+}(F;F_{3})\tilde{A}_{\rm
int}(F)\theta^{-}(F;F_{4})+\theta^{+}(F;F_{4})\tilde{A}_{\rm D}(F),$

where $F_{3}=0.014$, and $F_{4}=F_{\rm peak}$. $\tilde{A}_{\rm INS}(F)$ is the
inspiral amplitude in Eq. (26), $\tilde{A}_{\rm int}(F)$ is the intermediate
amplitude defined in Eq. (27), and $\tilde{A}_{\rm D}(F)$ is the PhenomD
amplitude.

Finally, the plus and cross polarization states of the complete frequency-
domain waveform read,

$\displaystyle\tilde{h}_{+}(f)=$ $\displaystyle\tilde{A}_{\rm
IMR}(f)\quantity(\frac{1+\cos^{2}\iota}{2})$ (35)
$\displaystyle\times\exp{-i(\Phi_{\rm IMR}(f)-\phi_{0}-2\pi ft_{0})},$

and

$\displaystyle\tilde{h}_{\crossproduct}(f)=$ $\displaystyle-i\tilde{A}_{\rm
IMR}(f)\cos\iota$ (36) $\displaystyle\times\exp{-i(\Phi_{\rm
IMR}(f)-\phi_{0}-2\pi ft_{0})}.$

Here $\iota$ is the angle of inclination of the binary plane to the line of
sight, $\phi_{0}$ and $t_{0}$ carry the overall phase and timeshift freedom of
the complete waveform. We call the final waveform model IMRPhenomD_Horizon,
and abbreviate it to PhenomD_Horizon.

In Fig. 3, we show the inspiral phase $\Psi(Mf)$, its derivative with respect
to $Mf$, and the difference between the phase of the hybrid data and the
model, for the same configurations as in Fig. 1. We show both the individual
fits from the ansatz in Eq. (24) and the final model prediction obtained from
the best fit values of ${\lambda^{i}_{jk}}$ from Eq. (30). In Fig. 4, we show
the inspiral and intermediate amplitude, rescaled by the leading order factor
$A_{0}=\frac{1}{\pi^{2/3}}\sqrt{5\eta/24}f^{-7/6}$. The intermediate regions
are shown in a blue shaded colour. We also show the amplitude of PhenomD on
the same plots.

##  VI Faithfulness of the model

# |  Simulation label |  $q$ |  $\chi_{1}$ |  $\chi_{2}$
---|---|---|---|---
1 |  SXS:BBH:0159 |  1 |  $-$0.9 |  $-$0.9
2 |  SXS:BBH:0154 |  1 |  $-$0.8 |  $-$0.8
3 |  SXS:BBH:0148 |  1 |  $-$0.44 |  $-$0.44
4 |  SXS:BBH:0150 |  1 |  0.2 |  0.2
5 |  SXS:BBH:0170 |  1 |  0.44 |  0.44
6 |  SXS:BBH:0155 |  1 |  0.8 |  0.8
7 |  SXS:BBH:0160 |  1 |  0.9 |  0.9
8 |  SXS:BBH:0157 |  1 |  0.95 |  0.95
9 |  SXS:BBH:0014 |  1.5 |  $-$0.5 |  0
10 |  SXS:BBH:0008 |  1.5 |  0 |  0
11 |  SXS:BBH:0013 |  1.5 |  0.5 |  0
12 |  SXS:BBH:0046 |  3 |  $-$0.5 |  $-$0.5
13 |  SXS:BBH:0036 |  3 |  $-$0.5 |  0
14 |  SXS:BBH:0168 |  3 |  0 |  0
15 |  SXS:BBH:0031 |  3 |  0.5 |  0
16 |  SXS:BBH:0047 |  3 |  0.5 |  0.5
17 |  SXS:BBH:0056 |  5 |  0 |  0
18 |  SXS:BBH:0181 |  6 |  0 |  0
19 |  SXS:BBH:1424 |  6.465 |  $-$0.66 |  $-$0.8
20 |  SXS:BBH:0298 |  7 |  0 |  0
Table 2: List of the extra hybrids created for validating the model. Figure 5:
Mismatches (%) between IMRPhenomD_Horizon and the hybrid waveforms (Table 1
and Table 2) in aLIGO ZDHP noise PSD with a lower cutoff of 10 Hz, as a
function of the total mass.

A phenomenological model, calibrated with a set of “target hybrids”, should be
compared with a larger set of “test hybrids” to assess its accuracy. The test
hybrids should contain both the target hybrids and a new set of hybrids
interspersed within the parameter space of calibration. To study the model’s
behavior beyond the calibration range, test hybrids can be extended further in
the parameter space. How well the final model can predict these test hybrids
is denoted by its effectualness and faithfulness [20, 72]. A model is
effectual if it is accurate enough to predict a GW signal from the detector
noise, while to be faithful it also has to have enough accuracy to estimate
the binary parameters. A faithful model is also effectual, but the converse
may not be true. Given any two signals $h(t)$ and $g(t)$, buried in the noise
$n(t)$, one can define a noise-weighted inner product

$\innerproduct{h}{g}=2\int_{0}^{\infty}\frac{\tilde{h}(f)\tilde{g}^{\ast}(f)+\tilde{h}^{\ast}(f)\tilde{g}(f)}{S_{n}(f)}\differential
f\,,$ (37)

where $\tilde{h}(f)$ is the Fourier transform of $h(t)$ defined by

$\tilde{h}(f)=\int_{-\infty}^{\infty}h(t)e^{-2\pi ift}\differential t\,,$ (38)

and $S_{n}(f)$ is the power spectral density (PSD) of the noise.

Equipped with this definition of an inner product between two waveforms, one
can define the match between them by maximizing the inner product with respect
to the constant phase and timeshift freedom between the waveforms
$\tilde{h}(f)$ and $\tilde{g}(f)$, normalized by their individual norms:

$\mathcal{M}=\underset{t_{c},\phi_{c}}{\text{min}}\,\,\frac{\innerproduct{h(\lambda)}{g(\lambda^{\prime})}}{\sqrt{\innerproduct{h(\lambda)}{h(\lambda)}\innerproduct{g(\lambda^{\prime})}{g(\lambda^{\prime})}}}\,.$
(39)

Here $\lambda$ and $\lambda^{\prime}$ denote the intrinsic parameters. The
mismatch, defined as $1-\mathcal{M}$, quantifies how much the waveform
$\tilde{h}(f)$ morphologically differs from $\tilde{g}(f)$. The faithfulness
of a model is quantified by demanding that the mismatches between the model
and the validation dataset (test hybrids) do not exceed a certain threshold.
For PhenomD, the threshold was 1%.

Since the total mass of the binary works as a scaling factor for the
gravitational waveforms, it defines how many GW cycles are contained within
the sensitive frequency band of a detector. The lower the total mass, the more
time the binary spends in its inspiral phase, resulting in a longer waveform
given a lower cutoff for the frequency. In Advanced LIGO, complete IMR models
become important above a critical value $\sim 12M_{\odot}$ [73] of the binary
mass. To assess our model’s faithfulness, we create a set of test hybrids that
are long enough to cover the sensitive frequency band of Advanced LIGO for a
total binary mass as low as $10M_{\odot}$. At 10 Hz, the dimensionless
frequency corresponding to a $10M_{\odot}$ binary is $Mf\approx 0.0005$.

Figure 6: Mismatches (%) between PhenomD_Horizon and PhenomD. (a) Mismatches
with only the phase correction but no amplitude correction in PhenomD_Horizon.
(b) Mismatches with only the amplitude correction without any phase
correction. (c) Mismatches with corrections in both the phase and amplitude.
The plots show that the dominant contributor to these mismatches is the phase
correction.

In Table 2 we list the extra test hybrids created for assessing the model
accuracy. Figure 5 shows the mismatches between the entire set of test hybrids
(Table 1 and Table 2) against the BBH waveform model presented in this paper.
We have used the Advanced LIGO zero-detuned high power [74] (ZDHP) noise curve
with a lower frequency cutoff of 10 Hz, and upper cutoff of $Mf=0.3$. The
mismatches are plotted against the total binary mass. We find that the model
accuracy is at par with the accuracy standard of PhenomD, with almost all the
hybrids showing mismatches well below 1%, the majority of them having
mismatches around 0.1% within the mass range $10M_{\odot}\leqslant M\leqslant
100M_{\odot}$. The worst mismatch we find is for the
$\\{q=4,\chi_{1}=0,\chi_{2}=0.8\\}$ configuration, which crosses the 1% level
above $\sim 50M_{\odot}$, and reaches a maximum of 1.21% at $70M_{\odot}$.

##  VII Model comparison

Figure 6 shows the mismatches between PhenomD_Horizon and PhenomD, keeping
$\chi_{1}=\chi_{2}$. We explore the mismatches within the parameter space
$q\in[1,8]$ and $\chi_{1},\chi_{2}\in[-0.9,0.9]$, for six values of the binary
mass, $M/M_{\odot}\in[10,20,30,40,50,60]$. In Fig. 6, we show the mismatches
including only the phase correction - i.e. using the recalibrated inspiral
phase, but keeping the IMR amplitude the same as PhenomD. Since the phase
correction due to TH written in Eq. (11) is $\propto 1/\eta$ at the leading
order, higher mismatches are found for higher mass ratios, for a fixed value
of the spins. With increasing spins, mismatches rise due to the increased
significance of TH. The mismatch reaches a maximum of $\sim 16\%$ for
$30M_{\odot}$ binaries with $q\gtrsim 7$ and $\chi_{1}=\chi_{2}\gtrsim 0.8$,
and a minimum of $\sim 10^{-5}\%$ for heavier binaries. For $M\leqslant
30M_{\odot}$, we see substantial regions of the parameter space with
mismatches between 1-10% for spin magnitudes of $|\chi_{1}|=|\chi_{2}|\gtrsim
0.5$. For $M>30M_{\odot}$, the 1% level is crossed only for $q\gtrsim 4$.

Figure 6 shows mismatches when only the amplitude correction is added to the
inspiral and the intermediate region is re-evaluated accordingly, but the
phase is kept the same as PhenomD. The mismatches, in this case, are orders of
magnitude less than the phase-corrected model, and never exceed 0.001% within
this parameter space. This result also vindicates our choice of not
recalibrating the amplitude model, as discussed in Sec V.2. We also note from
Eq. (15) that at the leading order, $\tilde{A}_{\rm TH}\propto\sqrt{\eta}$,
which implies that the TH corrections are more significant for more symmetric
masses. This is in contrast to the behaviour of the phase correction, and
renders the amplitude correction ineffective for high mass-ratio binaries with
notable TH contribution in the phase. Nevertheless, $\tilde{A}_{\rm TH}$ still
rises with increasing spins.

In Fig. 6, we use the complete IMRPhenomD_Horizon model, by including both the
phase and amplitude corrections. The contours of mismatches $\geqslant 0.1\%$
are almost the same as in Fig. 6, but amplitude corrections elevate them above
0.01% throughout other regions of the parameter space. However, mismatches
$<0.1\%$ are below the accuracy level of the model, for which we considered a
1% tolerance level of modeling errors (Fig. 5). This suggests that including
the amplitude correction in the TH corrected model has negligible effect
compared to the phase correction in GW data analysis.

To demonstrate the improvement in accuracy that our model introduces over
PhenomD, we use 219 NR waveforms from the SXS catalog and compare them against
PhenomD and PhenomD_Horizon separately by computing mismatches. In Fig. 7, we
show the distributions of the mismatches, by using the Advanced LIGO ZDHP
noise curve. We plot histograms of $\log_{10}(1-\mathcal{M})$, where
$(1-\mathcal{M})$ is the average mismatch in the binary mass range
$12M_{\odot}\leqslant M\leqslant 100M_{\odot}$. The dashed lines show the
medians of the two histogram plots, where a slight overall improvement is
seen. The fractional shift in the medians (denoted by an overbar), is

$\frac{\overline{(1-\mathcal{M})}_{\rm
PhenomD}-\overline{(1-\mathcal{M})}_{\rm
PhenomD\\_Horizon}}{\overline{(1-\mathcal{M})}_{\rm PhenomD}}\approx 0.04,$
(40)

showing an improvement of $\sim 4\%$. Figure 8 shows similar histograms with a
flat noise curve. In this case, also we see a shift in the median value
towards lower mismatches, with an improvement of $\sim 2\%$.

Figure 7: Mismatches of 219 non-precessing noneccentric NR waveforms from SXS
with PhenomD_Horizon and PhenomD, in aLIGO ZDHP noise curve, averaged in the
binary mass range $12-100M_{\odot}$. The dashed lines show the median values
of the distributions. Figure 8: Same as in Fig. 7, but with a flat noise
curve.

##  VIII Discussion and future work

We have presented a phenomenological gravitational waveform model for binary
black holes that carries the signatures of tidal heating of the holes
explicitly throughout the binary evolution, starting from the early inspiral
phase and ending at their merger. In the frequency domain, it contains
corrections due to this effect both in the phase and amplitude. The entire IMR
model for the phase is divided into two parts. The inspiral phase model,
defined by the frequency range $0.0035\leqslant Mf\leqslant 0.018$, is a
recalibrated version of IMRPhenomD [24, 25] to account for the modifications
due to tidal heating. The merger-ringdown parts, defined as $Mf\geqslant
0.018$ in terms of the dimensionless frequency, is the same as PhenomD since
this region contains purely NR data, which inherently incorporates the
dephasing due to tidal heating. The amplitude corrections are added by
dividing the model into three parts. The inspiral model, with $0.0035\leqslant
Mf\leqslant 0.014$, is constructed by directly adding the correction term to
the amplitude model of PhenomD in this region. The intermediate region,
$0.014\leqslant Mf\leqslant Mf_{\rm peak}$, $f_{\rm peak}$ being the frequency
corresponding to the peak amplitude, is approximated as a polynomial in $f$
and evaluated by imposing $C^{(1)}$ continuity at the two ends and a
collocation point at their mid-frequency. The amplitude model is identical to
PhenomD for $f>f_{\rm peak}$. The final model is calibrated within the range
in the mass ratio $q\in[1,8]$ and the effective spin parameter $\chi_{\rm
PN}\in[-0.79,0.72]$.

We conducted a study of faithfulness, by calculating mismatches between the
model and a validation dataset of 40 hybrid waveforms constructed with tidal
heating included. We demonstrated that modeling errors do not lead to
mismatches higher than 1% (barring a single outlier), with most of the
mismatch values lying around 0.1%. This standard of faithfulness is comparable
to PhenomD, which also considered a maximum tolerance level of 1% mismatch.

How the new model differs from its parent model PhenomD was assessed by
calculating mismatches between them within the calibration range of the
intrinsic parameter space, for binary masses ranging from $10M_{\odot}$ to
$60M_{\odot}$. Phase corrections in the inspiral result in mismatches between
these two models as high as $\sim 16\%$, whereas amplitude corrections alone
are unable to produce mismatches above $\sim 0.001\%$. With increasing mass
asymmetry, the impact of the phase correction increases, but the amplitude
correction becomes less effective. Combined with the fact that the former has
more significance in GW searches and PE, this reciprocal nature of these two
correction terms leaves the latter mostly superfluous for practical purposes.

We compared the model with a set of 219 NR waveforms from the SXS catalog
[53], and show the mismatch distribution along with the same for PhenomD,
averaged over the binary mass range $12M_{\odot}\leqslant M\leqslant
100M_{\odot}$. With the Advanced LIGO ZDHP noise curve, we see an improvement
of $\sim 4\%$ in the medians of the two distributions. With a flat noise
curve, the improvement is relatively modest, about 2%. These overall
improvements point towards a mild yet non-negligible increase in waveform
accuracy.

Tidal heating of black holes, after all, is a weak effect in the context of
comparable-mass binaries. Nevertheless, an absence of its signatures from the
complete waveform models can bias the estimated parameters of a binary. Its
absence can also mimic a deviation from GR predictions, raising complications
in tests of GR. We leave these studies to future work. Additionally, tidal
heating has the important property of varying significantly according to the
nature of the component objects, which makes it a potential discriminator for
BHs from other objects. However, to utilize the full power of TH as a BH
identifier, one needs to model the energy and angular momentum fluxes down the
BH horizons also in the late inspiral regime described by NR, where the tidal
fields are the strongest. The feasibility and effectiveness of modeling these
quantities can be explored in future studies, where the model presented in
this work can be used directly.

The current work, in its entirety, is based on the frequency-domain waveform
model IMRPhenomD. For this reason, we perform all the comparisons of the model
presented here with the said approximant only. More recent and state-of-the-
art phenomenological waveform models are available currently, namely, the
IMRPhenomX [23, 75] family of waveforms. These waveforms are calibrated with a
much larger set of hybrid waveforms, and use superior techniques like
hierarchical modeling [76] to reduce modeling errors, and also to generate
more accurate waveforms with unequal spins. These models boast 1 to 2 orders
of magnitude better faithfulness than PhenomD, and a dramatic shift in the
mismatch distribution against NR waveforms towards lower values [23].
Including the effects of TH in these models by recalibrating the inspiral is
beyond the scope of this work, which uses only publicly available NR data from
the SXS catalog. We leave them as future possibilities.

## Acknowledgments

SM thanks the University Grants Commission (UGC), India, for providing
financial support as a senior research fellowship. The high performance
computing facility Pegasus at IUCAA was used for this work. The authors thank
Vaishak Prasad, Apratim Ganguly, Anuj Mishra, Mark Scheel, Nathan Johnson-
McDaniel, P. Ajith, and Shrobana Ghosh for useful comments and suggestions.
The Python software packages PyCBC [77], Scipy [78], Numpy [79], Matplotlib
[80] were used for computational works, analyses, and plots. Mathematica [81]
was used for data fitting purposes. SB acknowledges the hospitality of IUCAA
with gratitude and support from the NSF under Grant PHY-2309352.

##  Appendix A Post-Newtonian inspiral phase

We write the 3.5PN frequency-domain phase (without TH) as

$\Psi_{\rm TF2}(f)=2\pi ft_{c}-\phi_{c}-\frac{\pi}{4}+\frac{3}{128\eta
v^{5}}\sum_{n=0}^{7}\psi_{(n)}v^{n}\,,$ (41)

where $v=(\pi Mf)^{1/3}$.

The coefficients $\\{\psi_{(n)}\\}$ ($n$=0-7) are given by [82, 83, 84],

$\displaystyle\psi_{(0)}=\leavevmode\nobreak\ $ $\displaystyle 1\,,$ (42)
$\displaystyle\psi_{(1)}=\leavevmode\nobreak\ $ $\displaystyle 0\,,$ (43)
$\displaystyle\psi_{(2)}=\leavevmode\nobreak\ $
$\displaystyle\frac{3715}{756}+\frac{55}{9}\eta\,,$ (44)
$\displaystyle\psi_{(3)}=\leavevmode\nobreak\ $
$\displaystyle-16\pi+\frac{113}{3}\delta\chi_{a}+\left(\frac{113}{3}-\frac{76}{3}\eta\right)\chi_{s}\,,$
(45) $\displaystyle\psi_{(4)}=\leavevmode\nobreak\ $
$\displaystyle\frac{15293365}{508032}+\frac{27145}{504}\eta+\frac{3085}{72}\eta^{2}-\frac{405}{4}\delta\chi_{a}\chi_{s}+\left(-\frac{405}{8}+200\eta\right)\chi_{a}^{2}+\left(-\frac{405}{8}+\frac{5}{2}\eta\right)\chi_{s}^{2}\,,$
(46) $\displaystyle\psi_{(5)}=\leavevmode\nobreak\ $
$\displaystyle\left[\frac{38645\pi}{756}-\frac{65\pi}{9}\eta-\left(\frac{732985}{2268}-\frac{24260}{81}\eta-\frac{340}{9}\eta^{2}\right)\chi_{s}-\left(\frac{732985}{2268}+\frac{140}{9}\eta\right)\delta\chi_{a}\right]\left(1+3\ln{v}\right)\,,$
(47) $\displaystyle\psi_{(6)}=\leavevmode\nobreak\ $
$\displaystyle\frac{11583231236531}{4694215680}-\frac{640\pi^{2}}{3}-\frac{6848}{21}\gamma_{E}+\eta\left(-\frac{15737765635}{3048192}+\frac{2255\pi^{2}}{12}\right)+\frac{76055}{1728}\eta^{2}-\frac{127825}{1296}\eta^{3}-\frac{6848}{21}\ln{(4v)}$
$\displaystyle+\pi\left\\{\frac{2270}{3}\delta\chi_{a}+\left(\frac{2270}{3}-520\eta\right)\chi_{s}\right\\}+\left(\frac{75515}{144}-\frac{8225}{18}\eta\right)\delta\chi_{a}\chi_{s}+\left(\frac{75515}{288}-\frac{263245}{252}\eta-480\eta^{2}\right)\chi_{a}^{2}$
$\displaystyle+\left(\frac{75515}{288}-\frac{232415}{504}\eta+\frac{1255}{9}\eta^{2}\right)\chi_{s}^{2}\,,$
(48) $\displaystyle\psi_{(7)}=\leavevmode\nobreak\ $
$\displaystyle\frac{77096675\pi}{254016}+\frac{378515\pi}{1512}\eta-\frac{74045\pi}{756}\eta^{2}+\left\\{-\frac{25150083775}{3048192}+\frac{10566655595}{762048}\eta-\frac{1042165}{3024}\eta^{2}+\frac{5345}{36}\eta^{3}\right.$
$\displaystyle+\left.\left(\frac{14585}{8}-7270\eta+80\eta^{2}\right)\chi_{a}^{2}\right\\}\chi_{s}+\left\\{\left(\frac{-25150083775}{3048192}+\frac{26804935}{6048}\eta-\frac{1985}{48}\eta^{2}\right)\chi_{a}\right.$
$\displaystyle\left.+\left(\frac{14585}{24}-2380\eta\right)\chi_{a}^{3}+\left(\frac{14585}{8}-\frac{215}{2}\eta\right)\chi_{a}\chi_{s}^{2}\right\\}\delta+\left(\frac{14585}{24}-\frac{475}{6}\eta+\frac{100}{3}\eta^{2}\right)\chi_{s}^{3}\,.$
(49)

Here $\delta=\sqrt{1-4\eta}$, $\chi_{s}=(\chi_{1}+\chi_{2})/2$,
$\chi_{a}=(\chi_{1}-\chi_{2})/2$. $\gamma_{E}$ is the Euler’s constant,
$\gamma_{E}=0.5772156\cdots$.

##  Appendix B Ringdown and peak frequency

The quasi-normal mode (QNM) frequency of a BH can be expressed in terms of its
oscillatory (real) and damping (imaginary) parts as,

$f_{\rm QNM}=f_{\rm RD}-if_{\rm damp}\,.$ (50)

$f_{\rm RD}$ can be obtained from the fitting formula [85]

$f_{\rm RD}=\frac{1}{2\pi
M_{f}}\left[f_{1}+f_{2}(1-\chi_{f})^{f_{3}}\right]\,,$ (51)

where $M_{f}$ and $\chi_{f}$ are the mass and spin of the merger remnant BH.
Evaluation of these quantities is described in Appendix C of Ref [86]. For the
dominant harmonic mode $l=m=2$, the coefficients are given by,
$f_{1}=1.5251,\,f_{2}=-1.1568,\,f_{3}=0.1292$.

$f_{\rm damp}$ can be written as

$f_{\rm damp}=\frac{f_{\rm RD}}{2Q},$ (52)

where the quality factor $Q$ can be expressed by the fitting formula

$Q=q_{1}+q_{2}(1-\chi_{f})^{q_{3}}.$ (53)

For the $l=m=2$ mode, $q_{1}=0.7,\,q_{2}=1.4187,\,q_{3}=-0.499$.

The frequency corresponding to the peak amplitude in the Fourier domain,
$f_{\rm peak}$, can then be obtained from the phenomenological relation

$f_{\rm peak}=\left|f_{\rm RD}+\frac{f_{\rm
damp}\gamma_{3}(\sqrt{1-\gamma_{2}^{2}}-1)}{\gamma_{2}}\right|\,,$ (54)

$\gamma_{2}$ and $\gamma_{3}$ being two phenomenological parameters that can
be evaluated for given values of $\\{\eta,\chi_{1},\chi_{2}\\}$ from the
ansatz in Eq. (30). The corresponding coefficients are given in Appendix C of
Ref. [25].

##  Appendix C Coefficients for the pseudo-PN parameters

|  $\sigma_{1}$ |  $\sigma_{2}$ |  $\sigma_{3}$ |  $\sigma_{4}$
---|---|---|---|---
$\lambda_{00}$ |  1650.156 |  14056.472 |  $-$86832.697 |  146221.791
$\lambda_{10}$ |  $-$3531.389 |  $-$75378.11 |  412796.518 |  $-$667838.551
$\lambda_{01}$ |  $-$2661.071 |  137530.39 |  $-$623322.321 |  920614.957
$\lambda_{11}$ |  15385.008 |  $-$944130.978 |  4.1416$\times 10^{6}$ |  $-$5.959$\times 10^{6}$
$\lambda_{21}$ |  $-$219.668 |  1.7166$\times 10^{6}$ |  $-$7.1273$\times 10^{6}$ |  9.8223$\times 10^{6}$
$\lambda_{02}$ |  $-$6845.194 |  169601.655 |  $-$737033.591 |  1.0493$\times 10^{6}$
$\lambda_{12}$ |  38637.011 |  $-$1.2802$\times 10^{6}$ |  5.4343$\times 10^{6}$ |  $-$7.5726$\times 10^{6}$
$\lambda_{22}$ |  $-$45127.308 |  2.648$\times 10^{6}$ |  $-$1.09142$\times 10^{7}$ |  1.4893$\times 10^{7}$
$\lambda_{03}$ |  $-$2076.335 |  50611.846 |  $-$213200.291 |  288184.3
$\lambda_{13}$ |  17433.08 |  $-$468141.986 |  1.9402$\times 10^{6}$ |  $-$2.641$\times 10^{6}$
$\lambda_{23}$ |  $-$31696.712 |  1.10645$\times 10^{6}$ |  $-$4.5663$\times 10^{6}$ |  6.2489$\times 10^{6}$
Table 3: Best-fit values of $\lambda^{i}_{jk}$ parameters in Eq. (30) for the
pseudo-PN parameters $\\{\sigma_{i}\\}$ in the inspiral phase.

## References

  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, (2021), arXiv:2111.03606 [gr-qc] .
  * Nitz  _et al._ [2023a] A. H. Nitz, S. Kumar, Y.-F. Wang, S. Kastha, S. Wu, M. Schäfer, R. Dhurkunde, and C. D. Capano, 4-OGC: Catalog of Gravitational Waves from Compact Binary Mergers, Astrophys. J. 946, 59 (2023a), arXiv:2112.06878 [astro-ph.HE] .
  * Olsen  _et al._ [2022] S. Olsen, T. Venumadhav, J. Mushkin, J. Roulet, B. Zackay, and M. Zaldarriaga, New binary black hole mergers in the LIGO-Virgo O3a data, Phys. Rev. D 106, 043009 (2022), arXiv:2201.02252 [astro-ph.HE] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, _et al._ , Advanced Virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Akutsu  _et al._ [2019] T. Akutsu _et al._ (KAGRA), KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector, Nature Astron. 3, 35 (2019), arXiv:1811.08079 [gr-qc] .
  * Owen and Sathyaprakash [1999] B. J. Owen and B. S. Sathyaprakash, Matched filtering of gravitational waves from inspiraling compact binaries: Computational cost and template placement, Phys. Rev. D 60, 022002 (1999), arXiv:gr-qc/9808076 .
  * Veitch  _et al._ [2015] J. Veitch _et al._ , Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library, Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Ashton  _et al._ [2019] G. Ashton _et al._ , BILBY: A user-friendly Bayesian inference library for gravitational-wave astronomy, Astrophys. J. Suppl. 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Romero-Shaw  _et al._ [2020] I. M. Romero-Shaw _et al._ , Bayesian inference for compact binary coalescences with bilby: validation and application to the first LIGO–Virgo gravitational-wave transient catalogue, Mon. Not. Roy. Astron. Soc. 499, 3295 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * Blanchet [2014] L. Blanchet, Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries, Living Rev. Rel. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Buonanno and Damour [1999] A. Buonanno and T. Damour, Effective one-body approach to general relativistic two-body dynamics, Phys. Rev. D 59, 084006 (1999), arXiv:gr-qc/9811091 .
  * Bohé  _et al._ [2017] A. Bohé _et al._ , Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors, Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Nagar and Rettegno [2019] A. Nagar and P. Rettegno, Efficient effective one body time-domain gravitational waveforms, Phys. Rev. D 99, 021501 (2019), arXiv:1805.03891 [gr-qc] .
  * Taracchini  _et al._ [2012] A. Taracchini, Y. Pan, A. Buonanno, E. Barausse, M. Boyle, T. Chu, G. Lovelace, H. P. Pfeiffer, and M. A. Scheel, Prototype effective-one-body model for nonprecessing spinning inspiral-merger-ringdown waveforms, Phys. Rev. D 86, 024011 (2012), arXiv:1202.0790 [gr-qc] .
  * Taracchini  _et al._ [2014] A. Taracchini _et al._ , Effective-one-body model for black-hole binaries with generic mass ratios and spins, Phys. Rev. D 89, 061502 (2014), arXiv:1311.2544 [gr-qc] .
  * Nagar  _et al._ [2018] A. Nagar _et al._ , Time-domain effective-one-body gravitational waveforms for coalescing compact binaries with nonprecessing spins, tides and self-spin effects, Phys. Rev. D 98, 104052 (2018), arXiv:1806.01772 [gr-qc] .
  * Gamba  _et al._ [2022] R. Gamba, S. Akçay, S. Bernuzzi, and J. Williams, Effective-one-body waveforms for precessing coalescing compact binaries with post-Newtonian twist, Phys. Rev. D 106, 024020 (2022), arXiv:2111.03675 [gr-qc] .
  * Nagar  _et al._ [2020] A. Nagar, G. Riemenschneider, G. Pratten, P. Rettegno, and F. Messina, Multipolar effective one body waveform model for spin-aligned black hole binaries, Phys. Rev. D 102, 024077 (2020), arXiv:2001.09082 [gr-qc] .
  * Ajith  _et al._ [2008] P. Ajith _et al._ , A Template bank for gravitational waveforms from coalescing binary black holes. I. Non-spinning binaries, Phys. Rev. D 77, 104017 (2008), [Erratum: Phys.Rev.D 79, 129901 (2009)], arXiv:0710.2335 [gr-qc] .
  * Santamaria  _et al._ [2010] L. Santamaria _et al._ , Matching post-Newtonian and numerical relativity waveforms: systematic errors and a new phenomenological model for non-precessing black hole binaries, Phys. Rev. D 82, 064016 (2010), arXiv:1005.3306 [gr-qc] .
  * Hannam  _et al._ [2014] M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Pratten  _et al._ [2020] G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni, A. Ramos-Buades, H. Estelles, and R. Jaume, Setting the cornerstone for a family of models for gravitational waves from compact binaries: The dominant harmonic for nonprecessing quasicircular black holes, Phys. Rev. D 102, 064001 (2020), arXiv:2001.11412 [gr-qc] .
  * Husa  _et al._ [2016] S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal, Phys. Rev. D 93, 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é, Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Tiglio and Villanueva [2022] M. Tiglio and A. Villanueva, Reduced order and surrogate models for gravitational waves, Living Rev. Rel. 25, 2 (2022), arXiv:2101.11608 [gr-qc] .
  * Field  _et al._ [2014] S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, and M. Tiglio, Fast prediction and evaluation of gravitational waveforms using surrogate models, Phys. Rev. X 4, 031006 (2014), arXiv:1308.3565 [gr-qc] .
  * Blackman  _et al._ [2017] J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, C. D. Ott, M. Boyle, L. E. Kidder, H. P. Pfeiffer, and B. Szilágyi, Numerical relativity waveform surrogate model for generically precessing binary black hole mergers, Phys. Rev. D 96, 024058 (2017), arXiv:1705.07089 [gr-qc] .
  * Kumar  _et al._ [2016] P. Kumar, T. Chu, H. Fong, H. P. Pfeiffer, M. Boyle, D. A. Hemberger, L. E. Kidder, M. A. Scheel, and B. Szilagyi, Accuracy of binary black hole waveform models for aligned-spin binaries, Phys. Rev. D 93, 104050 (2016), arXiv:1601.05396 [gr-qc] .
  * Blackman  _et al._ [2015] J. Blackman, S. E. Field, C. R. Galley, B. Szilágyi, M. A. Scheel, M. Tiglio, and D. A. Hemberger, Fast and Accurate Prediction of Numerical Relativity Waveforms from Binary Black Hole Coalescences Using Surrogate Models, Phys. Rev. Lett. 115, 121102 (2015), arXiv:1502.07758 [gr-qc] .
  * Varma  _et al._ [2019] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E. Kidder, and H. P. Pfeiffer, Surrogate model of hybridized numerical relativity binary black hole waveforms, Phys. Rev. D 99, 064045 (2019), arXiv:1812.07865 [gr-qc] .
  * Pürrer [2016] M. Pürrer, Frequency domain reduced order model of aligned-spin effective-one-body waveforms with generic mass-ratios and spins, Phys. Rev. D 93, 064041 (2016), arXiv:1512.02248 [gr-qc] .
  * Lackey  _et al._ [2017] B. D. Lackey, S. Bernuzzi, C. R. Galley, J. Meidam, and C. Van Den Broeck, Effective-one-body waveforms for binary neutron stars using surrogate models, Phys. Rev. D 95, 104036 (2017), arXiv:1610.04742 [gr-qc] .
  * Lackey  _et al._ [2019] B. D. Lackey, M. Pürrer, A. Taracchini, and S. Marsat, Surrogate model for an aligned-spin effective one body waveform model of binary neutron star inspirals using Gaussian process regression, Phys. Rev. D 100, 024002 (2019), arXiv:1812.08643 [gr-qc] .
  * Barkett  _et al._ [2016] K. Barkett _et al._ , Gravitational waveforms for neutron star binaries from binary black hole simulations, Phys. Rev. D 93, 044064 (2016), arXiv:1509.05782 [gr-qc] .
  * Dietrich  _et al._ [2019] T. Dietrich _et al._ , Matter imprints in waveform models for neutron star binaries: Tidal and self-spin effects, Phys. Rev. D 99, 024029 (2019), arXiv:1804.02235 [gr-qc] .
  * Poisson [2009] E. Poisson, Tidal interaction of black holes and Newtonian viscous bodies, Phys. Rev. D 80, 064029 (2009), arXiv:0907.0874 [gr-qc] .
  * Cardoso and Pani [2013] V. Cardoso and P. Pani, Tidal acceleration of black holes and superradiance, Class. Quant. Grav. 30, 045011 (2013), arXiv:1205.3184 [gr-qc] .
  * Hartle [1973] J. B. Hartle, Tidal Friction in Slowly Rotating Black Holes, Phys. Rev. D8, 1010 (1973).
  * Poisson [2004] E. Poisson, Absorption of mass and angular momentum by a black hole: Time-domain formalisms for gravitational perturbations, and the small-hole / slow-motion approximation, Phys. Rev. D70, 084044 (2004), arXiv:gr-qc/0407050 [gr-qc] .
  * Alvi [2001] K. Alvi, Energy and angular momentum flow into a black hole in a binary, Phys. Rev. D64, 104020 (2001), arXiv:gr-qc/0107080 [gr-qc] .
  * Chatziioannou  _et al._ [2016] K. Chatziioannou, E. Poisson, and N. Yunes, Improved next-to-leading order tidal heating and torquing of a Kerr black hole, Phys. Rev. D 94, 084043 (2016), arXiv:1608.02899 [gr-qc] .
  * Datta [2023] S. Datta, Horizon fluxes of binary black holes in eccentric orbits, (2023), arXiv:2305.03771 [gr-qc] .
  * Datta  _et al._ [2021] S. Datta, K. S. Phukon, and S. Bose, Recognizing black holes in gravitational-wave observations: Challenges in telling apart impostors in mass-gap binaries, Phys. Rev. D 104, 084006 (2021), arXiv:2004.05974 [gr-qc] .
  * Mukherjee  _et al._ [2022] S. Mukherjee, S. Datta, S. Tiwari, K. S. Phukon, and S. Bose, Toward establishing the presence or absence of horizons in coalescing binaries of compact objects by using their gravitational wave signals, Phys. Rev. D 106, 104032 (2022), arXiv:2202.08661 [gr-qc] .
  * Cardoso  _et al._ [2016] V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, Gravitational-wave signatures of exotic compact objects and of quantum corrections at the horizon scale, Phys. Rev. D 94, 084031 (2016), arXiv:1608.08637 [gr-qc] .
  * Datta  _et al._ [2020] S. Datta, R. Brito, S. Bose, P. Pani, and S. A. Hughes, Tidal heating as a discriminator for horizons in extreme mass ratio inspirals, Phys. Rev. D 101, 044004 (2020), arXiv:1910.07841 [gr-qc] .
  * Datta and Bose [2019] S. Datta and S. Bose, Probing the nature of central objects in extreme-mass-ratio inspirals with gravitational waves, Phys. Rev. D 99, 084001 (2019), arXiv:1902.01723 [gr-qc] .
  * Nagar and Akcay [2012] A. Nagar and S. Akcay, Horizon-absorbed energy flux in circularized, nonspinning black-hole binaries and its effective-one-body representation, Phys. Rev. D 85, 044025 (2012), arXiv:1112.2840 [gr-qc] .
  * Bernuzzi  _et al._ [2012] S. Bernuzzi, A. Nagar, and A. Zenginoglu, Horizon-absorption effects in coalescing black-hole binaries: An effective-one-body study of the non-spinning case, Phys. Rev. D 86, 104038 (2012), arXiv:1207.0769 [gr-qc] .
  * Damour  _et al._ [2013] T. Damour, A. Nagar, and S. Bernuzzi, Improved effective-one-body description of coalescing nonspinning black-hole binaries and its numerical-relativity completion, Phys. Rev. D 87, 084035 (2013), arXiv:1212.4357 [gr-qc] .
  * Damour and Nagar [2014] T. Damour and A. Nagar, New effective-one-body description of coalescing nonprecessing spinning black-hole binaries, Phys. Rev. D 90, 044018 (2014), arXiv:1406.6913 [gr-qc] .
  * Boyle  _et al._ [2019] M. Boyle _et al._ , The SXS Collaboration catalog of binary black hole simulations, Class. Quant. Grav. 36, 195006 (2019), arXiv:1904.04831 [gr-qc] .
  * Buonanno  _et al._ [2009a] A. Buonanno, Y. Pan, H. P. Pfeiffer, M. A. Scheel, L. T. Buchman, and L. E. Kidder, Effective-one-body waveforms calibrated to numerical relativity simulations: Coalescence of non-spinning, equal-mass black holes, Phys. Rev. D 79, 124028 (2009a), arXiv:0902.0790 [gr-qc] .
  * Tichy  _et al._ [2000] W. Tichy, E. E. Flanagan, and E. Poisson, Can the postNewtonian gravitational wave form of an inspiraling binary be improved by solving the energy balance equation numerically?, Phys. Rev. D 61, 104015 (2000), arXiv:gr-qc/9912075 .
  * Ajith [2011] P. Ajith, Addressing the spin question in gravitational-wave searches: Waveform templates for inspiralling compact binaries with nonprecessing spins, Phys. Rev. D 84, 084037 (2011), arXiv:1107.1267 [gr-qc] .
  * Pretorius [2005] F. Pretorius, Evolution of binary black hole spacetimes, Phys. Rev. Lett. 95, 121101 (2005), arXiv:gr-qc/0507014 .
  * Campanelli  _et al._ [2006] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Accurate evolutions of orbiting black-hole binaries without excision, Phys. Rev. Lett. 96, 111101 (2006), arXiv:gr-qc/0511048 .
  * Baker  _et al._ [2006] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter, Gravitational wave extraction from an inspiraling configuration of merging black holes, Phys. Rev. Lett. 96, 111102 (2006), arXiv:gr-qc/0511103 .
  * Pfeiffer [2012] H. P. Pfeiffer, Numerical simulations of compact object binaries, Class. Quant. Grav. 29, 124004 (2012), arXiv:1203.5166 [gr-qc] .
  * Centrella  _et al._ [2010] J. Centrella, J. G. Baker, B. J. Kelly, and J. R. van Meter, Black-hole binaries, gravitational waves, and numerical relativity, Rev. Mod. Phys. 82, 3069 (2010), arXiv:1010.5260 [gr-qc] .
  * Le Tiec [2014] A. Le Tiec, The Overlap of Numerical Relativity, Perturbation Theory and Post-Newtonian Theory in the Binary Black Hole Problem, Int. J. Mod. Phys. D 23, 1430022 (2014), arXiv:1408.5505 [gr-qc] .
  * Hannam [2014] M. Hannam, Modelling gravitational waves from precessing black-hole binaries: Progress, challenges and prospects, Gen. Rel. Grav. 46, 1767 (2014), arXiv:1312.3641 [gr-qc] .
  * Scheel  _et al._ [2009] M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews, and H. P. Pfeiffer, High-accuracy waveforms for binary black hole inspiral, merger, and ringdown, Phys. Rev. D 79, 024003 (2009), arXiv:0810.1767 [gr-qc] .
  * Ossokine  _et al._ [2013] S. Ossokine, L. E. Kidder, and H. P. Pfeiffer, Precession-tracking coordinates for simulations of compact-object-binaries, Phys. Rev. D 88, 084031 (2013), arXiv:1304.3067 [gr-qc] .
  * Scheel  _et al._ [2015] M. A. Scheel, M. Giesler, D. A. Hemberger, G. Lovelace, K. Kuper, M. Boyle, B. Szilágyi, and L. E. Kidder, Improved methods for simulating nearly extremal binary black holes, Class. Quant. Grav. 32, 105009 (2015), arXiv:1412.1803 [gr-qc] .
  * Cook and Whiting [2007] G. B. Cook and B. F. Whiting, Approximate Killing Vectors on S**2, Phys. Rev. D 76, 041501 (2007), arXiv:0706.0199 [gr-qc] .
  * Christodoulou and Ruffini [1971] D. Christodoulou and R. Ruffini, Reversible transformations of a charged black hole, Phys. Rev. D 4, 3552 (1971).
  * MacDonald  _et al._ [2011] I. MacDonald, S. Nissanke, H. P. Pfeiffer, and H. P. Pfeiffer, Suitability of post-Newtonian/numerical-relativity hybrid waveforms for gravitational wave detectors, Class. Quant. Grav. 28, 134002 (2011), arXiv:1102.5128 [gr-qc] .
  * Damour  _et al._ [2011] T. Damour, A. Nagar, and M. Trias, Accuracy and effectualness of closed-form, frequency-domain waveforms for non-spinning black hole binaries, Phys. Rev. D 83, 024006 (2011), arXiv:1009.5998 [gr-qc] .
  * Bruegmann  _et al._ [2008] B. Bruegmann, J. A. Gonzalez, M. Hannam, S. Husa, U. Sperhake, and W. Tichy, Calibration of Moving Puncture Simulations, Phys. Rev. D 77, 024027 (2008), arXiv:gr-qc/0610128 .
  * Lindblom  _et al._ [2008] L. Lindblom, B. J. Owen, and D. A. Brown, Model Waveform Accuracy Standards for Gravitational Wave Data Analysis, Phys. Rev. D 78, 124020 (2008), arXiv:0809.3844 [gr-qc] .
  * Buonanno  _et al._ [2009b] A. Buonanno, B. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, Comparison of post-Newtonian templates for compact binary inspiral signals in gravitational-wave detectors, Phys. Rev. D 80, 084043 (2009b), arXiv:0907.0700 [gr-qc] .
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (KAGRA, LIGO Scientific, Virgo, VIRGO), Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Rev. Rel. 21, 3 (2018), arXiv:1304.0670 [gr-qc] .
  * Pratten  _et al._ [2021] G. Pratten _et al._ , Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes, Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Jiménez-Forteza  _et al._ [2017] X. Jiménez-Forteza, D. Keitel, S. Husa, M. Hannam, S. Khan, and M. Pürrer, Hierarchical data-driven approach to fitting numerical relativity data for nonprecessing binary black holes with an application to final spin and radiated energy, Phys. Rev. D 95, 064024 (2017), arXiv:1611.00332 [gr-qc] .
  * Nitz  _et al._ [2023b] A. Nitz, I. Harry, D. Brown, C. M. Biwer, J. Willis, T. D. Canton, C. Capano, T. Dent, L. Pekowsky, S. De, M. Cabero, A. R. Williamson, B. Machenschalk, D. Macleod, G. S. C. Davies, P. Kumar, F. Pannarale, S. Reyes, dfinstad, S. Kumar, M. Tápai, L. Singer, S. Khan, S. Fairhurst, veronica villa, K. Chandra, S. Wu, A. Nielsen, S. Singh, and T. Massinger, gwastro/pycbc: v2.0.6 release of pycbc (2023b).
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Meth. 17, 261 (2020), arXiv:1907.10121 [cs.MS] .
  * Harris  _et al._ [2020] C. R. Harris _et al._ , Array programming with NumPy, Nature 585, 357 (2020), arXiv:2006.10256 [cs.MS] .
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering 9, 90 (2007).
  * Inc. [2023] W. R. Inc., Mathematica, Version 13.3 (2023), Champaign, IL.
  * Arun  _et al._ [2009] 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, Phys. Rev. D 79, 104023 (2009), [Erratum: Phys.Rev.D 84, 049901 (2011)], arXiv:0810.5336 [gr-qc] .
  * Mishra  _et al._ [2016] C. K. Mishra, A. Kela, K. G. Arun, and G. Faye, Ready-to-use post-Newtonian gravitational waveforms for binary black holes with nonprecessing spins: An update, Phys. Rev. D 93, 084054 (2016), arXiv:1601.05588 [gr-qc] .
  * Wade  _et al._ [2013] M. Wade, J. D. E. Creighton, E. Ochsner, and A. B. Nielsen, Advanced LIGO’s ability to detect apparent violations of the cosmic censorship conjecture and the no-hair theorem through compact binary coalescence detections, Phys. Rev. D 88, 083002 (2013), arXiv:1306.3901 [gr-qc] .
  * Nakano  _et al._ [2019] H. Nakano, T. Narikawa, K.-i. Oohara, K. Sakai, H.-a. Shinkai, H. Takahashi, T. Tanaka, N. Uchikata, S. Yamamoto, and T. S. Yamamoto, Comparison of various methods to extract ringdown frequency from gravitational wave data, Phys. Rev. D 99, 124032 (2019), arXiv:1811.06443 [gr-qc] .
  * Favata  _et al._ [2022] M. Favata, C. Kim, K. G. Arun, J. Kim, and H. W. Lee, Constraining the orbital eccentricity of inspiralling compact binary systems with Advanced LIGO, Phys. Rev. D 105, 023003 (2022), arXiv:2108.05861 [gr-qc] .