# First gravitational-wave search for intermediate-mass black hole mergers
with higher order harmonics

Koustav Chandra<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Bombay, Powai, Mumbai 400 076, India  Juan Calderón
Bustillo<EMAIL_ADDRESS>Instituto Galego de Física de Altas
Enerxías, Universidade de Santiago de Compostela, 15782 Santiago de
Compostela, Galicia, Spain  Department of Physics, The Chinese University of
Hong Kong, Shatin, N.T., Hong Kong  Archana Pai<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Bombay, Powai, Mumbai
400 076, India  I. W. Harry  University of Portsmouth, Portsmouth, PO1 3FX,
United Kingdom<EMAIL_ADDRESS>
###### Abstract

Current matched-filter searches for gravitational waves from binary black-hole
mergers compare the calibrated detector data to waveform templates that omit
the higher-order mode content of the signals predicted by General Relativity.
However, higher-order emission modes become important for highly inclined
asymmetric sources with masses above $\simeq 100M_{\odot}$, causing current
searches to be ill-suited at detecting them. We present a new gravitational-
wave search that implements templates including higher-order modes, adapted
signal-glitch discriminators, and trigger-ranking statistics to specifically
target signals displaying strong higher modes, corresponding to nearly edge-on
sources with total redshifted masses in the intermediate-mass black-hole range
$M_{T}(1+z)\in(100,500)M_{\odot}$ and mass-ratios $q\in(1,10)$. Our search
shows a volumetric sensitivity gain of up to 450% to these signals compared to
existing searches omitting higher-order modes. We deploy our search on public
data from the third observing run of Advanced LIGO. While we find no
statistically significant candidates beyond those already reported elsewhere,
our search sets the stage to search for higher-mode rich signals in future
observing runs. The efficient detection of such signals is crucial to
performing detailed tests of General Relativity, observing strong-field
phenomena, and maximizing the chances of observing the yet uncharted realm of
intermediate-mass black hole binaries.

intermediate-mass black hole — gravitational waves

##  I Introduction

The population of binary black hole mergers identified by the current
generation ground-based gravitational-wave detector network has helped us
uncover several interesting and unexpected features about the intrinsic
properties of black holes in the local Universe (Abbott  _et al._ , 2019,
2021a, 2021b). While we expect that most of these merger components are
stellar remnants, there are certain gravitational wave events whose origin can
also be explained by a hierarchical formation pathway (Rodriguez  _et al._ ,
2020; Zevin  _et al._ , 2021; Fragione  _et al._ , 2020; Kimball  _et al._ ,
2021; Baibhav  _et al._ , 2021; Gerosa and Fishbach, 2021; Tagawa  _et al._ ,
2021). In this formation channel, second (or higher) generation black holes
merge to form increasingly massive black hole depending on the properties of
the host environment. As a result, hierarchically assembled black holes can
populate the pair-instability mass gap and can help explain the formation of
intermediate-mass black holes in the local Universe (Fishbach  _et al._ ,
2017; Gerosa and Berti, 2017; Rodriguez  _et al._ , 2019; Kimball  _et al._ ,
2020; Doctor  _et al._ , 2021; Mapelli  _et al._ , 2021; Fragione  _et al._ ,
2022).

Figure 1: Time-frequency maps of two different non-precessing binary black
hole signals with weak higher-order modes (left) and strong higher modes
(right) injected into Advanced LIGO Livingston data. The left panel
corresponds to a signal from a face-on binary black hole $(\iota=0)$ with a
mass ratio of $q\sim 1.2$ while the right panel corresponds to an edge-on
$(\iota=\pi/2)$ binary with a mass ratio of $q\sim 10.1$. In both cases, we
have used a detector frame total mass of $M_{T}(1+z)\sim 500M_{\odot}$. The
left panel displays a (short) chirp morphology, while the right shows a more
complex (multi-chirp) structure due to its higher-mode content. The lack of
inspiral in the above panels is due to the high total mass of the system.

Binaries producing intermediate-mass black hole remnants are particularly
interesting. Firstly, the merger remnants could be seeds that grow into
supermassive black holes in the galactic nuclei, thus providing essential
feedback on galaxy evolution (Quinlan and Shapiro, 1987; Greene  _et al._ ,
2020). Secondly, because of their large total mass
($\sim\mathcal{O}(10^{2}~{}M_{\odot}$)), the frequency of the gravitational
wave produced during the merger, and the ringdown phase is near the most
sensitive band of existing gravitational wave detectors, therefore offering
the best scenarios to study the behavior of gravity in its most extreme regime
(Abbott  _et al._ , 2016a; Carullo  _et al._ , 2019; Isi  _et al._ , 2019).
Furthermore, during these stages, subdominant gravitational wave emission
modes can get strongly triggered depending on the properties of the source. If
observed, these modes can provide crucial information on the behavior of the
final object, enabling to test the nature of it through black hole
spectroscopy (Dreyer  _et al._ , 2004; Berti  _et al._ , 2016; Bhagwat  _et
al._ , 2018; Carullo  _et al._ , 2019; Isi  _et al._ , 2019; Cabero  _et al._
, 2020; Bustillo  _et al._ , 2021) or the observation of phenomena of crucial
astrophysical relevance like gravitational recoil (Calderón Bustillo  _et al._
, 2018a, 2019; Varma  _et al._ , 2020; Calderón Bustillo  _et al._ , 2022).

The full inspiral-merger-ringdown gravitational wave emission of compact
binary mergers can be computed through either semi-analytical (Santamaria  _et
al._ , 2010; Khan  _et al._ , 2016; Husa  _et al._ , 2016; Buonanno and
Damour, 2000; Damour  _et al._ , 2013; Buonanno and Damour, 1999; Cotesta  _et
al._ , 2018) or numerical techniques Szilagyi  _et al._ (2009); Jani  _et al._
(2016); Healy  _et al._ (2017). Therefore, such gravitational waves can be
extracted from the noisy detector data through the optimal method of matched-
filtering (Wainstein and Zubakov, 1962; Sathyaprakash and Dhurandhar, 1991;
Allen  _et al._ , 2012). This is the cross-correlation of the detector data
with pre-computed waveform templates. For the effectiveness of matched-
filtering, the search templates need to be faithful representations of the
incoming gravitational wave signal; otherwise, the searches might miss them.

However, while General Relativity predicts that gravitational waves are a
superposition of several emission modes $h_{\ell,m}$, current template-based
searches only implement the dominant quadrupole modes, given by $h_{2,\pm 2}$
(Usman  _et al._ , 2016; Venumadhav  _et al._ , 2019; Sachdev  _et al._ ,
2019; Aubin  _et al._ , 2021; Chu  _et al._ , 2022). Such a strategy has been
demonstrated to effectively detect signals from face-on (or face-off) systems
with nearly equal mass and total redshifted masses $\lesssim 100M_{\odot}$,
for which non-quadrupolar (or higher-order) modes contribute negligibly. Also,
such optimally-oriented, symmetric sources are intrinsically luminous and
hence easier to detect. However, asymmetric black-hole binaries emit
gravitational waves with strong higher harmonics, especially during the merger
and ringdown stages (Mills and Fairhurst, 2021) (See Figure 1 for a
qualitative example). The contribution of higher modes in the observed signal
increases as the orbital inclination of the system deviates from face-on/off,
especially impacting those of high mass, for which the frequency of the
dominant harmonic lies below the optimal sensitivity of the detector. Their
omission in searches, therefore, dramatically reduces the sensitivity to these
sources (Capano  _et al._ , 2014; Calderón Bustillo  _et al._ , 2016, 2017,
2018b; Chandra  _et al._ , 2020), potentially causing an “observational bias”
against asymmetric intermediate-mass black hole binaries with large orbital
inclinations.

Harry  _et al._ (2018) developed a prototype matched-filter-based search for
non-spinning binary black hole sources containing higher-order harmonics,
demonstrating its ability to recover synthetic signals in simulated Advanced
LIGO noise that is free of instrumental transients or “glitches”. These noisy
artifacts tend to affect search efficiency severely. Here, we turn this
prototype search into a fully working one and deploy it on data collected
during the third observing run of Advanced LIGO. First, we expand the search
to the case of aligned-spin (non-precessing) sources. Second, to make the
search effective when applied to real data, we adapt existing signal-glitch
discriminators to separate the noisy transients from astrophysical ones
better. While we do not find any new statistically significant signal
candidates beyond those already reported elsewhere (Abbott  _et al._ , 2021a,
b; Chandra  _et al._ , 2021a; Nitz  _et al._ , 2021a, b; Olsen  _et al._ ,
2022), we show that our search is up to 450% more sensitive than past searches
with overlapping parameter space, mainly when the target sources are at nearly
edge-on orientation.

The rest of this paper is organized as follows. In Sec. II we review the fully
generic search method. Sec. III presents the details of our coincidence
analysis with two detectors, and Sec. IV assesses the benefits of deploying
our search. In Sec. V we present the results of our search on O3 data, and we
conclude in Sec. VI.

##  II Background

Figure 2: Ratio of the optimal signal-to-noise ratio of waveforms including
only the dominant $(\ell,m)=(2,\pm 2)$ mode and including the modes
$(\ell,|m|)=(2,1),~{}(2,2),~{}(3,3),~{}(4,4)$ and $(5,5)$, generated with
identical parameters. The top panel shows results for varying the source’s
total mass and mass ratio for fixed edge-on inclinations. The bottom panel
shows results for varying mass ratio and inclination for a fixed redshifted
total mass of $M_{T}(1+z)=300~{}M_{\odot}$. In both panels, we have set the
black hole spins to zero, the azimuth to $\phi=0$, and the antenna response
$F_{+}=1$ and $F_{\times}=0$. We consider a detector characterized by the
Advanced LIGO Livingston sensitivity during the O3 run.

The fundamental assumption behind most gravitational wave data analysis is
that the detector output $d(t)$ is composed of two additive components, namely
the noise $n(t)$ and the signal strain $s(t)$ as:

$d(t)=n(t)+s(t)$ (1)

The noise $n(t)$, a stochastic process, varies randomly with time. Assuming
that the noise is wide-sense stationary and Gaussian, its statistical
properties are fully described by its one-sided power spectral density
$S_{n}(f)$. The gravitational wave signal, on the other hand, is
deterministic, and it will impart a strain $s(t)$ that is parameterized by a
vector $\bm{\lambda}$ and is given by the following linear combination:

$s(t;\bm{\lambda})=F_{+}h_{+}(t-t_{c};D_{L},\bm{\theta})+F_{\times}h_{\times}(t-t_{c};D_{L},\bm{\theta})~{}.$
(2)

Here, $F_{+/\times}$ are the sky-location and polarisation angle-dependent
antenna response patterns of the detector to the two gravitational wave
polarisations, $h_{+/\times}$. $t_{c}$ is the merger time of the signal, and
$D_{L}$ is the luminosity distance to the source. The morphology and the
evolution of these two polarisation states depend on the properties and
orientation of the source, which we denote using the vector $\bm{\theta}$.

We can measure the loudness of a signal as it appears in a detector by
calculating the optimal signal-to-noise ratio:

$\rho_{\mathrm{opt}}=\sqrt{(s|s)}$ (3)

Here, $(a|b)$ is the real part of the noise-weighted inner product:

$\langle
a|b\rangle=4\int^{f_{\mathrm{max}}}_{f_{\mathrm{min}}}df~{}\frac{\tilde{a}^{\ast}(f)\tilde{b}(f)}{S_{n}(f)},$
(4)

between two real-valued time series. The tilde here denotes the Fourier
transform of the corresponding time-domain data. $f_{\mathrm{min}}$ and
$f_{\mathrm{max}}$ represent upper and lower frequency cutoffs, and the symbol
$\ast$ denotes complex conjugation.

In principle, we can model this signal strain for any given gravitational wave
source. This way, we can construct waveform templates $h(t;\bm{\lambda})$ for
the expected signal. When templates are available, the optimal way to retrieve
the signals from the noise is via matched-filtering whose output is the
matched-filter signal-to-noise ratio (SNR) (Wainstein and Zubakov, 1962;
Sathyaprakash and Dhurandhar, 1991; Allen  _et al._ , 2012):

$\rho=\frac{(d|h(\bm{\lambda}))}{(h(\bm{\lambda})|h(\bm{\lambda}))^{1/2}}~{}.$
(5)

However, the parameters $\bm{\lambda}$ are not known a priori. Therefore
matched-filter-based searches construct and use discrete banks of template
waveforms spanning the search space to compare the data and numerically
maximize the SNR. In certain situations, we can make simplifying assumptions
about the morphology of the incoming signal and can hence analytically
maximize over some of these parameters. This greatly reduces the computational
cost of the search, for it reduces the dimensionality of the parameter space
over which we need to maximize the SNR numerically.

###  II.1 Gravitational Wave Higher Harmonics

The gravitational wave emission from a compact merger can be expressed as a
superposition of different GW emission modes $h_{\ell,m}$ weighted by spin-2
spherical harmonics ${}^{-2}Y_{\ell,m}$ (Goldberg  _et al._ , 1967; Blanchet,
2014):

$h_{+}-ih_{\times}=\frac{1}{D_{L}}\sum_{\ell=2}^{\infty}\sum_{m=-l}^{l}~{}^{-2}Y_{\ell,m}(\iota,\phi)h_{\ell,m}(t-t_{c};\bm{\Xi})$
(6)

Above, $(\iota,\phi)$ define the polar and azimuthal angles of a spherical
coordinate system centered at the center of mass of the system, and $\bm{\Xi}$
collectively denotes the intrinsic parameters of the source, namely the
individual masses $m_{1,2}$ and spins $\vec{\chi}_{1,2}$.

In order to isolate the contribution from each orientation parameter, is it
useful to decompose the harmonics $Y^{-2}_{\ell,m}$ in terms of an overall
amplitude term that depends on the source inclination and an overall phase
term that depends on the azimuth as:

${}^{-2}Y_{\ell,m}=A_{\ell,m}(\iota)e^{-im\phi}.$

This makes it obvious that the inclination angle $\iota$ determines the
amplitude of each mode while the azimuth $\phi$ determines the way the modes
combine, which can dramatically change the signal morphology observed by
different observers around the source. Similarly, we can also express the
$(\ell,m)$ emission mode of a gravitational wave signal as

$h_{\ell,m}=\mathcal{A}_{\ell,m}(t-t_{c};\bm{\Xi})e^{-i\Phi_{\ell,m}(t-t_{c};\bm{\Xi})}$

.

During most of the inspiral part of non-precessing quasi-circular mergers, the
emission is vastly dominated by the quadrupolar mode or $(\ell,|m|)=(2,2)$.
However, the impact of higher-order modes in the signal becomes significant
during the merger and ringdown stages. This effect becomes enhanced for high-
mass systems, for which the frequency of the dominant harmonic can lie below
the detector’s sensitive band. In addition, the amplitude
$\mathcal{A}_{\ell,m}$ of these modes relative to the $\mathcal{A}_{2,\pm 2}$
grows with increasing mass ratio, making the higher-order modes very relevant
for asymmetric mass sources (Berti  _et al._ , 2007; Mills and Fairhurst,
2021). Finally, while the spherical harmonics other than $Y^{-2}_{2,\pm 2}$ is
mostly zero for face-on (or face-off) orientations $(\iota=0,\pi)$, they reach
their maxima for intermediate inclinations, making the higher-order modes have
a stronger impact for large orbital inclinations, e.g. for edge-on
orientations $(\iota=\pi/2)$.

To visualize the above, we, in Figure 2, shows the ratio of the optimal
signal-to-noise ratio of the $h_{+}$ polarisations of simulated waveforms from
non-spinning binary black hole systems when these are generated using only the
dominant quadrupole and using the modes
$(\ell,|m|)=(2,1),~{}(2,2),~{}(3,3),~{}(4,4)$ and $(5,5)$. The top panel shows
how the ratio varies as a function of the redshifted total mass and mass ratio
for (fixed) edge-on inclinations. The bottom panel shows the same as a
function of mass ratio and inclination for a fixed redshifted total mass of
$300M_{\odot}$. For approximately $51\%$ of the simulated signals, we find
that the quadrupole mode contributes less than $75\%$ to the SNR, indicating
that the sub-dominant harmonics significantly contribute to the overall signal
loudness, especially for asymmetric, nearly edge-on sources. This makes it
obvious that higher modes are crucial to represent such signals accurately.

###  II.2 Generic SNR Statistic

Irrespective of the parameters of the template, we can always write the
template strain, in the Fourier domain, as:

$\displaystyle\tilde{h}(f)$
$\displaystyle=F_{+}\sqrt{(h_{+}|h_{+})}~{}\hat{h}_{+}+F_{\times}\sqrt{(h_{\times}|h_{\times})}~{}\hat{h}_{\times}$
(7) $\displaystyle=A(u\hat{h}_{+}+\hat{h}_{\times})$

by defining:

$\displaystyle u=$
$\displaystyle\frac{F_{+}\sqrt{(h_{+}|h_{+})}}{F_{\times}\sqrt{(h_{\times}|h_{\times})}}$
(8) $\displaystyle A=$ $\displaystyle
F_{\times}\sqrt{(h_{\times}|h_{\times})}$ $\displaystyle\hat{h}_{+/\times}=$
$\displaystyle\frac{h_{+/\times}}{\sqrt{(h_{+/\times}|h_{+/\times})}}$

This indicates that it is always possible to express a generic gravitational
wave transient in terms of an overall amplitude term $A$ and a weighted linear
combination of the unit-normalized gravitational wave polarisations.

If we substitute Eq. (7) in Eq. (5), we get:

$\rho=\frac{(d|u\hat{h}_{+}+\hat{h}_{\times})}{(u\hat{h}_{+}+\hat{h}_{\times}|u\hat{h}_{+}+\hat{h}_{\times})^{1/2}}$
(9)

Maximising $\rho$ over the $u-$dependence yields (Harry  _et al._ , 2018):

$\max_{u}\rho^{2}=\frac{(d|\hat{h}_{+})^{2}+(d|\hat{h}_{\times})^{2}-2(d|\hat{h}_{+})(d|\hat{h}_{\times})(\hat{h}_{+}|\hat{h}_{\times})}{1-(\hat{h}_{+}|\hat{h}_{\times})^{2}}~{},$
(10)

which by construction is effectively maximized over the luminosity distance to
the source, the sky location of the source, and the polarisation angle. Using
an inverse Fast Fourier transform routine, we can also carry out the
maximization over $t_{c}$. Therefore, we can use this SNR statistic to search
for any gravitational wave transient irrespective of their morphology,
provided the waveform of the transient is well-modeled. For that reason, we
use this statistic to perform our search with higher-order mode templates that
are completely parameterized by the vector
$\bm{\theta}=(m_{1},~{}m_{2},~{}\chi_{1z},~{}\chi_{2z},~{}\iota,~{}\phi)$.

###  II.3 SNR statistic used in current searches: Non-precessing Quadrupolar
Limit

Figure 3: Magnitude of the imaginary part of the complex overlap $\langle
h_{+}|h_{\times}\rangle$ between the two polarisations for a simulated
population of quasi-circular binary black hole signals with higher-mode
content. The top panel shows this as a function of the total mass and mass
ratio, while the middle panel shows the same as a function of effective
aligned spin and mass ratio. We have numerically minimized the overlap over
the inclination angle and the azimuth for each system. The bottom panel shows
the minimum overlap as a function of the corresponding inclination angle and
effective aligned spin.

As already stated, most current matched-filter-based searches make several
assumptions about the signal model, which is not true when using waveforms
with higher harmonics (Usman  _et al._ , 2016; Venumadhav  _et al._ , 2019;
Sachdev  _et al._ , 2019; Aubin  _et al._ , 2021; Chu  _et al._ , 2022).
Specifically, they use a SNR statistic that assumes that the polarisation
states are related as $\tilde{h}_{\times}\propto i\tilde{h}_{+}$. The generic
SNR statistic in Eq. (10) for such an assumption reduces to:

$\max_{u}\rho^{2}=(d|\hat{h}_{+})^{2}+(d|\hat{h}_{\times})^{2}=|\langle
d|\hat{h}_{+}\rangle|^{2},$ (11)

as $(\hat{h}_{+}|\hat{h}_{\times})=0$. If we limit to dominant harmonics of
quasi-circular black hole binaries, the constant of proportionality is:

$\epsilon=-\frac{2\cos{\iota}}{1+\cos^{2}{\iota}}$

and its magnitude varies between $0$ (for edge-on) and $1$ (face-on/off).

Also, current searches assume that we can absorb the binary’s sky location,
orientation, luminosity distance, and corresponding polarisation angle by
applying an appropriate amplitude and phase scaling to the observed waveform,
both of which can be analytically maximized. Therefore, current matched-filter
searches only iteratively search over the source’s intrinsic parameters
$\bm{\Xi}$. Both assumptions hold if and only if we restrict to the dominant
harmonics of a quasi-circular binary black hole merger. However, they break
when we use waveforms, including higher-order modes as templates.

###  II.4 Evaluating the need for the generic SNR statistic for higher-order
mode searches

We assess the necessity of the generic SNR statistic for higher-order mode
searches by computing the magnitude of the imaginary component of the complex
overlap between the unit-normalized gravitational-wave polarisations:

$O=\mathrm{Im}\langle\hat{h}_{+}|\hat{h}_{\times}\rangle$ (12)

for O3 Advanced LIGO noise sensitivity. We perform this study with waveforms
generated using the reduced-order representation of the aligned-spin
effective-one-body model with higher-order modes, SEOBNRv4HM that includes the
spherical harmonics $(\ell,|m|)=(2,1),~{}(3,3),~{}(4,4)$ and $(5,5)$ beyond
the dominant quadrupolar mode (Cotesta  _et al._ , 2018, 2020). The simulated
waveforms imitate gravitational waves from a synthetic population of quasi-
circular black hole binaries with detector frame (redshifted) total mass
$M_{T}(1+z)\in(100,500)M_{\odot}$, mass-ratio $q=m_{1}/m_{2}\in(1,10)$ and
spins $\chi_{1z,2z}\in(-0.998,0.998)$. We distribute these sources
isotropically over the inclination angle and reference orbital phase, and then
for each of these waveforms, we compute $O$. Finally, following (Harry  _et
al._ , 2018), we numerically minimize $O$ over the inclination angle and
phase. Fig 3 shows the minimum value of $O$ as a function of total mass and
mass ratio in the top panel and as a function of mass ratio and effective
aligned spin in the middle panel.

We find that for certain configurations, the minimum value of the overlap goes
to $\sim 0.87$. Also, for $\sim 3.4\%$ of the cases, $O<0.95$, indicating that
we cannot, in general, assume $\tilde{h}_{+}(f)\propto i\tilde{h}_{\times}(f)$
for quasi-circular binary waveforms. Further, most of these low overlap
binaries have either a high total mass and/or a negative
$\chi_{\mathrm{eff}}$, indicating that their duration is short within the
detector bandwidth. Also, these low overlap binaries are oriented nearly edge-
on, as shown in the bottom panel of Fig 3. Therefore, we target binaries with
$\iota\in(75^{\circ},105^{\circ})$ and use Eq. (10) for SNR calculation.

##  III Methods for observing generic black hole binaries

There are two matched-filter analyses that explicitly search for binary black
hole mergers producing lower mass range $(100,600)~{}M_{\odot}$ intermediate-
mass black hole remnants. The first of the two, namely PyCBC-IMBH (Chandra
_et al._ , 2021a), targets quasi-circular binaries with detector frame total
masses between 100 and 600 $M_{\odot}$, with component masses greater than
$40M_{\odot}$ and mass ratio $q$ between $1$ and $10$. The search, however,
uses waveforms with just the dominant harmonics and a matched-filter SNR
statistic defined in Eq (11). Also, to reduce the number of false alarms due
to short-duration glitches (Cabero  _et al._ , 2020), this search does not use
any template with a duration less than $70$ms, measured from the fixed
starting frequency of $15$Hz.

The other intermediate-mass black hole binary-specific search is constructed
using the GstLAL software package (Sachdev  _et al._ , 2019; Cannon  _et al._
, 2021), and it also uses a template bank of quasi-circular dominant (2,2)
harmonic with $q\in(1,10)$ and assumes that $\tilde{h}_{\times}\propto
i\tilde{h}_{+}$. However, unlike PyCBC-IMBH, it targets systems with
$M_{T}(1+z)$ between 50 and 600 $M_{\odot}$, and it uses a starting frequency
of $10$ Hz for its matched-filter operation. Other than these differences, the
searches also use different signal-noise discriminators and rank coincident
triggers differently (For details, see Sec. 3.2 of Abbott  _et al._ (2022) and
the references therein).

Both of these analyses need to incorporate knowledge of the higher harmonics
into a search for gravitational waves, making them poor at detecting nearly
edge-on quasi-circular binaries, which are mass asymmetric and/or massive. In
what follows, we describe our search strategy that involves the construction
of a higher-order mode bank and adapting existing signal-glitch discriminators
to perform a search on real gravitational wave data.

Finally, we note that search algorithms beyond those relying on matched-
filtering have also been used to search for intermediate-mass black hole
binaries in LIGO-Virgo data. In particular, the template-independent search
algorithm Coherent WaveBurst (Klimenko  _et al._ , 2016; Szczepańczyk  _et
al._ , 2021) identifies gravitational wave signals by looking for coherent
power excess across different detectors using minimal assumptions on the
morphology of the expected signal. This makes this search potentially
sensitive toward a wider variety of intermediate-mass black hole binaries.
Currently, the sensitivity of the Coherent WaveBurst search (in its
intermediate-mass black hole configuration) is comparable with that of the
PyCBC-IMBH search (Chandra  _et al._ , 2021b).

###  III.1 Search space

Motivated by Figure 2, we build a template bank to target highly inclined
sources with redshifted total mass beyond $100~{}M_{\odot}$. The specifics of
our target domain are summarised in Table 1.

Parameter |
---|---
Total redshifted mass | $M_{T}(1+z)\in(100,~{}500)M_{\odot}$
Mass ratio | $q=m_{1}/m_{2}\in(1,~{}5)~{}\&~{}(5,~{}10)$
Spin z-component | $\chi_{1z,2z}\in(-0.998,~{}0.998)$
Orbital inclination | $\iota\in(75^{\circ},~{}105^{\circ})$
Azimuth | $\phi\in(0,2\pi)$
Table 1: Summary of the target parameter space covered by our template banks.
Note that we consider two banks ${\cal{B}}_{1}$ with $q\in(1,5)$ and
${\cal{B}}_{2}$ with $q\in(5,10)$.

While this search can be expanded to larger regions of the parameter space of
black hole binaries, particularly with lower orbital inclinations, we have
opted to restrict to the above constraints for the following reasons. Firstly,
templates with lower orbital inclinations will identify intrinsically louder
sources than templates with higher inclinations. Therefore, the net
sensitivity of our search would be dominated by that to face-on binaries
(which do not display higher modes), damaging our ability to evaluate how our
search improves on existing ones when targeting higher-mode rich signals. In
addition, including more templates in our bank would increase our false alarm
rate due to the increased number of templates, further hindering the
sensitivity to the sources we want to target. Secondly, as mentioned before,
sources at a high inclination produce complex, higher-mode-rich signals, while
face-on ones produce rather morphologically simple ones. Therefore, it is
reasonable to expect the corresponding two types of templates to show a
different propensity to be triggered by glitches.

Consequently, we preferred to isolate these two types of potential background
populations. Thirdly, the upper bound on the detector frame total mass is
chosen to ensure a minimum template duration, which should be larger than one
cycle for effective matched filtering. In fact, in preliminary analyses, we
found that templates for heavier systems are more susceptible to glitches,
resulting in a poorer overall search sensitivity. Similarly, we keep the mass
ratio within $q\leq 10$. Finally, due to similar arguments, we divide our
search into two separate template banks: the first – ${\cal{B}}_{1}$ –
targeting systems with $q\in(1,5)$ and the second – ${\cal{B}}_{2}$ –
targeting systems with $q\in(5,10)$. This helps protect the lower mass ratio
space from the high penalty – due to increased glitchiness – of the high mass
ratio region. Lastly, while less physically relevant, enlarging the target
parameter space would have led to a higher number of templates, equating to
higher computational costs.

Figure 4: The top panel shows the distribution of template waveforms in the
component mass space for our two banks. The red templates belong to the low-
mass-ratio ($q\in[1,5]$) region (${\cal{B}}_{1}$), and the blue templates
belong to the high mass-ratio ($q\in[5,10]$) region (${\cal{B}}_{2}$). The
bottom panel shows the distribution of the orbital inclination angle $\iota$
of the templates, which is limited to $\iota\in(75^{\circ},105^{\circ})$.

###  III.2 Template Bank

Figure 5: Effectualness of our template banks to a family of randomly
generated waveforms within their respective target parameter spaces. For
comparison, we also show the effectualness of the PyCBC-IMBH bank to the same
set of waveforms, which spans all of our parameter space but lacks higher-
order modes in its templates. As can be seen, the PyCBC-IMBH bank owing to its
lack of higher-order modes, is largely un-effectual (fitting factor below
0.97) within our target space.

While several such models exist, in this work, we choose the aligned-spin
effective-one-body model SEOBNRv4HM-ROM. This model includes modes
$(\ell,|m|)=(2,1),~{}(3,3),~{}(4,4)$ and $(5,5)$ beyond the quadrupolar mode
(Cotesta  _et al._ , 2018, 2020). We choose a minimum frequency cutoff for the
matched filter of $15$Hz. We note that, in Gaussian background noise, a lower
low-frequency cutoff would increase our sensitivity to binary black hole
signals simply due to increased SNR due to contributions from low frequencies.
In non-Gaussian noise, however, this also increases the chances that our
templates are triggered by glitches, yielding a larger background that would
compensate (and even overcome) the aforementioned gain in sensitivity.

We build each of our banks, namely ${\cal{B}}_{1}$ and ${\cal{B}}_{2}$, with a
stochastic template placement algorithm (Harry  _et al._ , 2008; Babak, 2008;
Harry  _et al._ , 2009) and a minimal match criterion of $0.97$. The template
placement algorithm relies on choosing a random “template ” $h_{T}$ from our
target search space, checking whether the bank, $\mathfrak{B}$ has a fitting
factor (FF):

$\mathrm{FF}=\max_{t_{c},~{}u,~{}h_{i}\in\mathfrak{B}}(h_{T}|h_{i})$ (13)

less than $0.97$ toward the template and accepting it based on it (Apostolatos
_et al._ , 1994). If otherwise, we reject the template. We repeat this process
unless we reach a sufficiently large rejection rate. The minimal match choice
of $0.97$ ensures that the maximum SNR loss due to the discreteness in the
template bank is not more than $3\%$. While a denser template bank with a
larger minimal match would lead to a smaller SNR loss – e.g., the PyCBC-IMBH
bank uses a minimum match of 0.99 – and, in principle, raise the sensitivity
of the search, this would also increase its computational cost, which is
already high due to the increase in the number of templates coming from the
addition of the orientation parameters to the bank and the need to perform two
filters per template.

Figure 4 shows the resulting template banks. The top panel shows the two banks
in terms of the component masses of the binary black hole. The bottom panel
shows the orbital inclination distribution of the templates. While
${\cal{B}}_{1}$ has 8626 templates, ${\cal{B}}_{2}$ has 40915 templates, i.e.,
almost $4.7$ times larger. Consequently, ${\cal{B}}_{2}$ will have a larger
background, making it more prone to false alarms. We discuss this impact in
terms of the search sensitivity in Sec. IV. We will call the search
implementing ${\cal{B}}_{1}$ as search-1 and the one implementing
${\cal{B}}_{2}$ as search-2.

We check the “effectualness” of our banks to their target space by calculating
fitting factor towards $\sim 100,000$ SEOBNRv4HM waveforms randomly
distributed across the target parameter space of the banks, and we summarise
our findings in Figure 5 where we plot our corresponding fitting factor
distributions. We note that for both our banks, the effectualness is higher
than the target value of $0.97$ for $\sim 98\%$ of the target signals. For
comparison, we also show the effectualness of the bank used by the PyCBC-IMBH
search towards the same set of simulated signals. The latter has a minimum
recovered fitting factor of $0.64$ and falls below $0.97$ for more than $15\%$
of the target signals showing that such a bank is ineffectual in our target
search space.

###  III.3 Dealing with instrumental transients

Figure 6: Impact of inclusion and omission of the higher-order modes on
signal-glitch consistency tests for current searches and the one presented
here. In the top panels, we plot the SNR\- $\chi^{2}_{\mathrm{r}}$
distribution for a simulated population of binary black hole obeying the
constraints of ${\cal{B}}_{1}$. In the bottom panels, we do the same for a
simulated population of binary black hole obeying the constraints of
${\cal{B}}_{2}$. The left panels show results corresponding to our banks,
including higher modes, while the right panels show the results corresponding
to the PyCBC-IMBH search. The dashed lines show contours of constant ranking
statistic $\tilde{\rho}$. The simulated signals, which contain higher modes,
have a greater mismatch with the waveforms implemented in the PyCBC-IMBH bank,
leading to significantly larger $\chi^{2}_{\mathrm{r}}$-values for the same
SNR. As a result, the search will misinterpret the simulated astrophysical
signals as noise triggers, damaging its sensitivity. The lack of triggers
beyond $\chi_{\mathrm{r}}^{2}=10$ is due to our choice of vetoing such
triggers.

####  III.3.1 Single-detector signal-glitch discriminator

In the presence of wide-sense stationary Gaussian noise, a consistently high
SNR across the detector network would have sufficed to assess the presence of
gravitational wave in the data. Advanced LIGO noise is, however, known to be
neither wide-sense stationary nor Gaussian (Cabero  _et al._ , 2020). Instead,
the detector data contains short-duration noise transients or glitches that
can produce large SNRs, mimicking transient gravitational wave, therefore
significantly damaging the search sensitivity. Consequently, it is necessary
to implement signal-glitch discrimination techniques or “vetoes” to identify
and penalize such glitches. PyCBC-based searches check whether a trigger: (a)
is an outlier in the calibrated whitened data stream (Usman  _et al._ , 2016),
(b) has a morphology that is consistent with the best-matched template Allen
(2005), (c) has any excess power beyond the maximum frequency of the best-
matched template Nitz (2018), (d) has any excess power (summed over bands) on
particular timescales Mozzon  _et al._ (2020) and (e) is consistent across
detectors Usman  _et al._ (2016); Nitz  _et al._ (2017); Davies  _et al._
(2020).

These tests’ outputs are numbers used to appropriately amend the trigger SNR
to suppress the noisy triggers. For instance, the output of the
$\chi_{\mathrm{r}}^{2}$ test (b above) is commonly combined with the SNR to
yield the “re-weighted SNR”:

$\tilde{\rho}=\begin{cases}\rho\Big{[}\frac{1}{2}\big{(}1+(\chi^{2}_{\mathrm{r}})^{3}\big{)}\Big{]}^{-1/6}&\mathrm{for\,}\chi^{2}_{\mathrm{r}}>1\\\
\rho&\mathrm{for\,}\chi^{2}_{\mathrm{r}}\leq 1\end{cases}$ (14)

This statistic down-ranks triggers whose morphology is not consistent with the
template (Babak  _et al._ , 2013), characterized by values of the
$\chi_{\mathrm{r}}^{2}>1$. Such a situation can occur if the trigger is due to
noise or it is due to an astrophysical signal that is not well-modeled by the
template, as will be the case when signals with higher-order modes are
filtered with templates that do not contain these.

In Figure 6 we have plotted $\chi^{2}_{\mathrm{r}}$ values as a function of
the SNR for a set of simulated waveforms that have been added to a
representative section of Advanced LIGO Livingston data. These simulated
waveforms are generated using the SEOBNRv4HM model, corresponding to signals
from simulated binaries within the respective target spaces of the two
template banks. The parameters of these simulated sources are uniformly
distributed in $M_{T}(1+z)$ and $m_{1}/(m_{1}+m_{2})$ space, isotropically
across the sky sphere, uniformly over the polarisation angles and uniformly in
comoving volume between bounds of $0.26~{}\mathrm{Gpc}^{3}$ and
$40~{}\mathrm{Gpc}^{3}$. The top panels correspond to waveforms within the
space spanned by ${\cal{B}}_{1}$, while the bottom ones correspond to
${\cal{B}}_{2}$. The left panels show the results of our searches, while the
right ones correspond to the case where we use the PyCBC-IMBH search. It is
noticeable that for the same SNR the latter search returns much larger values
of the $\chi_{\mathrm{r}}^{2}$ than ours as the templates used by it are the
dominant harmonics of a quasi-circular binary black hole as against the former
where higher harmonics are present in the template waveforms. This will lead
to a “false identification” of signals as glitches, damaging the search
sensitivity. For better separation between the background triggers and the
simulation, we discard triggers that are highly inconsistent with the
templates by placing a threshold of 10 on $\chi^{2}_{\mathrm{r}}$ value.
Additionally, we put a similar threshold on the short-term variation of the
noise power spectral density (Mozzon  _et al._ , 2020) following Chandra  _et
al._ (2021a) to alleviate the effects of loud broadband detector noise.

####  III.3.2 Background-Dependent Reweighing

Figure 7: Complementary cumulative distribution of single detector triggers
(solid lines) binned by detector-frame total mass $M_{T}(1+z)$, effective
inspiral spin parameter $\chi_{\mathrm{eff}}$ and symmetric mass ratio
$\eta=q/(1+q)^{2}$, fitted with an exponential polynomial distribution (dotted
lines) parameterized by $a,~{}b$ and $c$. The triggers correspond to those
identified in a representative portion of Advanced LIGO Livingston data from
the third observing run by a representational set of template groups.

The result of these consistency tests is a set of numbers, $v_{i}$ that are
combined with the SNR to produce a single-detector statistic
$\mathfrak{r}(\rho;v_{i})$ 111For details see Sec.II C of Chandra  _et al._
(2021a) and references therein.. We expect the probability distribution of
$\mathfrak{r}$ for triggers associated with noise in each detector $I$ beyond
a threshold $\mathfrak{r}_{0}$ for each template $k$ is a falling exponential,
with template-dependent parameter $a$. The complementary cumulative
distribution of such a probability distribution is (Nitz  _et al._ , 2017)

$\bar{F}_{1}(\mathfrak{r})=e^{-a(\mathfrak{r}-\mathfrak{r}_{0})}~{}.$ (15)

The above follows from the fact that the noise probability obeys an
inhomogeneous Poisson process and that different templates identify different
types of triggers. Consequently, certain regions of the bank will produce
large values of $\mathfrak{r}$, potentially reducing the statistical
significance of triggers coming from regions of the bank less prone to loud
glitches. Therefore, we use an additional template-dependent parameter,
$\mu_{I}^{k}$, that accounts for the total noise-trigger rate in a detector.
The two are combined to obtain a model for the noise rate density for each
template (Nitz  _et al._ , 2017):

$r_{I}^{k}=-\mu_{I}^{k}~{}\partial_{\mathfrak{r}}\bar{F}_{1}(\mathfrak{r})=\mu_{I}^{k}ae^{-a(\mathfrak{r}-\mathfrak{r}_{0})}$
(16)

which is used to separate the noise triggers from signal triggers further.
Therefore template-dependent background reweighing uses the likelihood of a
template to identify a trigger associated with noise of certain loudness to
discriminate it. Generally, this fit could be performed separately for each
template, but PyCBC-based searches choose to perform the fitting over a group
of templates that behave similarly in the presence of noise. This increases
the number of triggers for the fit and hence improves the fit.

However, we find that templates with higher harmonics produce noise trigger
distribution that is not well reproduced by Eq. 15. Instead, we find
experimentally that the noise probability is better fitted to a model whose
complementary cumulative distribution is:

$\bar{F}_{2}(\mathfrak{r})=e^{-[a(\mathfrak{r}-\mathfrak{r}_{0})^{2}+b(\mathfrak{r}-\mathfrak{r}_{0})+c]}$
(17)

We show this resulting fit in Figure 7. The solid lines correspond to single-
detector triggers we obtained while performing the search in $\sim 30$ days of
O3 data from Advanced LIGO-Livingston. The triggers are grouped based on the
templates’ total mass, effective spins, and symmetric mass ratio. The dotted
lines are the estimated fits that follow the noise trigger distribution well.
We note that the triggers correspond to a representative set of templates. For
example, the blue line shows the distributions of triggers identified by
templates with $M_{T}(1+z)\in(400,500)M_{\odot}$,
$\eta=q/(1+q)^{2}\in(0.11,~{}0.12)$ and $\chi_{\mathrm{eff}}\in(-0.6,-0.4)$
while the red line shows the same for templates with
$M_{T}(1+z)\in(100,200)M_{\odot}$, $\eta=q/(1+q)^{2}\in(0.08,~{}0.09)$ and
$\chi_{\mathrm{eff}}\in(0.7,0.9)$.

We use these noise rate densities to estimate the total rate, $r_{N}^{k}$, of
coincident events in a template. Under the assumption that the noise in both
detectors is independent, this total noise rate estimate is proportional to
the product of the rate of noise triggers in each detector for a template:

$r^{k}_{N}=-\prod_{I}\mu_{I}^{k}\partial_{\mathfrak{r}}\bar{F}_{2}(\mathfrak{r})$
(18)

One can use this total noise rate estimate at the multi-detector level to
improve the statistical significance of astrophysical candidates, as described
in the next subsection.

####  III.3.3 Multi-Detector Ranking Statistics

Under the assumption that the detector noise is Gaussian and uncorrelated, it
is sufficient to use the network SNRs,
$\sqrt{\rho_{\mathrm{H}}^{2}+\rho_{\mathrm{L}}^{2}}$, to rank coincident
triggers from Hanford and Livingston detector (Pai  _et al._ , 2001).
Similarly, in the presence of non-gaussian glitches, the expression
$\sqrt{\mathfrak{r}_{\mathrm{H}}^{2}+\mathfrak{r}_{\mathrm{L}}^{2}}$ provides
a more suitable ranking of the coincident triggers (Babak  _et al._ , 2013).
We improve upon this by using a ranking statistic, $\mathcal{R}$ (Nitz  _et
al._ , 2017):

$\mathcal{R}\propto\frac{1}{\sqrt{2}}(\mathfrak{r}_{H,0}+\mathfrak{r}_{L,0}-\ln{r}_{N}^{k})$
(19)

that includes the network noise rate estimate for a given template. It is
designed to reduce to:

$\mathcal{R}\propto\frac{1}{\sqrt{2}}(\mathfrak{r}_{H}+\mathfrak{r}_{L}+\mathcal{O}(\mathfrak{r}^{2}))$
(20)

when the fit coefficient $b$ in both detectors is one. The pre-factor of
$1/\sqrt{2}$ follows from the fact that if the rescaled SNRs in both the
detectors is the same, then
$\sqrt{\mathfrak{r}_{H}^{2}+\mathfrak{r}_{L}^{2}}\sim\sqrt{2}\mathfrak{r}_{H}$.
In both detectors, improvements to the ranking statistic using the relative
probability distribution of nuisance parameters, such as $u$ and $A$, are left
for future work.

###  III.4 Estimation of statistical significance

(a) Ratio of sensitive VT of our search-1 implementing the bank
${\cal{B}}_{1}$ for and the PyCBC-IMBH search.

(b) Ratio of sensitive VT of our search-1 implementing the bank
${\cal{B}}_{2}$ and the PyCBC-IMBH search.

Figure 8: Ratio of the volume-time sensitivity (VT) at a fixed inverse-false-
alarm rate (IFAR) of the search presented in this work, which includes higher-
order modes, with respect to the existing PyCBC-IMBH search, which omits
higher modes. We show results as a function of the total redshifted mass of
the target signals. The upper panel shows results corresponding to our low
mass-ratio bank $\cal{B}_{1}$ targeting mass ratios $q\in(1,5)$ while the
bottom one corresponds to our large mass-ratio bank $\cal{B}_{2}$ targeting
mass ratios $q\in(5,10)$. In both cases, our search improves on existing ones,
especially for short-lived high-mass target signals.

We use $\mathcal{R}$ to estimate the statistical significance of coincident
triggers. Ideally, we would like to calculate it using data free of
gravitational waves. But there is no way to shield our detectors from incoming
gravitational waves. However, it is possible to destroy the signal coherence
across the detector network by shifting the data from one of the detectors
with respect to other detector(s) by a time duration greater than the light
travel time between the detectors (Abbott  _et al._ , 2016b). Because the
noise in each detector is assumed to be uncorrelated, time slides preserve the
noise properties. It also helps simulate “background triggers”, none of which
are due to a given astrophysical signal. The distribution of the ranking
statistic of these background triggers is used to compute inverse false alarm
rates (IFARs), and any “foreground” or “zero-lag” trigger that exceeds a pre-
determined inverse false alarm rate (IFAR) threshold is deemed as a detection
candidate.

Finally, since there can be chance coincidences between noise and signal
triggers, our search hierarchically removes loud triggers like other PyCBC-
based searches to minimize their impact (Usman  _et al._ , 2016) as well as
background triggers that lie within a given time window around events with
large IFAR.

##  IV Search Sensitivity

To assess the benefits of our new search against searches that do not use
higher-order mode waveforms, we estimate the sensitive volume-time product,
$\mathrm{VT}_{\mathrm{sen}}$, using a set of simulated waveforms that we have
added to $\sim 100$ days of O3 Advanced LIGO data. As described in Sec.
III.3.1, our simulation set consists of binary black-hole waveforms with
strong higher-order mode content. We compare our search performance against
the PyCBC-IMBH search, used to search for intermediate-mass black-hole mergers
in O3 data whose templates cover a parameter space including ours but ignore
higher-order modes. We choose to compare the volumetric search sensitivity
against PyCBC-IMBH (Chandra  _et al._ , 2021a) as this search performs better
than or comparable to other searches for intermediate-mass black hole binaries
(Chandra  _et al._ , 2021b).

The sensitive volume-time product of a given search measures the number of
expected signals from a population of binaries that the search can detect
beyond a pre-determined statistical significance. We estimate this sensitive
volume time using the Monte Carlo method with importance sampling. This
involves sampling from a proxy or proposal distribution rather than the
nominal or target distribution such that most of the simulated signals closely
straddle the boundary between being detected and being missed. One can find
details of the method used in the Appendix of Capano  _et al._ (2016).

Figure 8(a) & Figure 8(b) compare our searches’ sensitive volume time estimate
to the estimate from the PyCBC-IMBH search at different thresholds of IFAR.
The first plot is for a population of binary black holes with $q\lesssim 5$
while the latter is for binaries with $q\in(5,10)$, and we have used search-1
and search-2, respectively. In both cases, and independently of the IFAR
threshold, we observe that the searches implementing higher modes are more
sensitive than those ignoring them. At low IFAR a clear separation from
injections and glitches is unnecessary, and the SNR recovery dominates the
ranking statistic. Therefore, the difference in sensitivity is mainly driven
by the difference in fitting factor caused by the lack of higher modes in the
quadrupole search. In contrast, at large IFAR, better separation between
glitches and injections is needed. Therefore the detected simulated signals
need to pass all the signal-glitch discriminator checks, e.g., the
$\chi_{\mathrm{r}}^{2}$ test. For this reason, the quadrupole search further
downweighs, resulting in an increased relative sensitivity of the higher-mode
search.

The observed strength of the higher-order modes in the simulated signals grows
as the total mass increases. Also, most of the in-band signal is dominated by
the merger and ringdown stages. For this reason, we observe that the higher-
mode search yields a larger sensitivity gain for increasing total mass. In
particular, the sensitivity gain of search-1 (search-2) goes from a factor of
$3.81\pm 0.44$ $(1.18\pm 0.17)$ for $M_{T}(1+z)\in(100,~{}200)M_{\odot}$ to
$4.48\pm 0.35$ $(2.25\pm 0.20)$ for $M_{T}(1+z)\in[300,500]M_{\odot}$. We note
that the reduced sensitivity gain of the search-2 is because its bank,
${\cal{B}}_{2}$, has a significantly larger amount of templates and,
therefore, suffers from an increased background.

Event Name | GPS Time [s] | Search-1 | Search-2
---|---|---|---
$\mathrm{IFAR}~{}[\mathrm{yr}]$ | $\rho_{\mathrm{H}}$ | $\rho_{\mathrm{L}}$ | $\mathrm{IFAR}~{}[\mathrm{yr}]$ | $\rho_{\mathrm{H}}$ | $\rho_{\mathrm{L}}$
$\mathrm{GW}190408\\_181802$ | $1238782700.3$ | 3.44 | 7.74 | 7.80 | 60.60 | 9.04 | 8.02
$\mathrm{GW}190412$ | $1239082262.2$ | - | - | - | 10.50 | 6.81 | 10.63
$\mathrm{GW}190503\\_185404$ | $1240944862.3$ | 3.52 | 8.81 | 7.41 | - | - | -
$\mathrm{GW}190513\\_205428$ | $1241816086.8$ | 1.49 | 7.31 | 6.95 | - | - | -
$\mathrm{GW}190517\\_055101$ | $1242107479.8$ | 1.88 | 5.97 | 7.82 | - | - | -
$\mathrm{GW}190519\\_153544$ | $1242315362.4$ | 31.35 | 8.52 | 9.56 | 18.82 | 7.07 | 8.96
$\mathrm{GW}190521$ | $1242442967.4$ | 1.13 | 7.56 | 11.85 | - | - | -
$\mathrm{GW}190521\\_074359$ | $1242459857.4$ | 43069.81 | 10.52 | 19.18 | 598.19 | 8.97 | 17.05
$\mathrm{GW}190602\\_175927$ | $1243533585.1$ | 2.07 | 6.49 | 10.16 | - | - | -
$\mathrm{GW}190706\\_222641$ | $1246487219.3$ | 108.26 | 8.80 | 8.24 | - | - | -
$\mathrm{GW}190727\\_060333$ | $1248242632.0$ | 27.65 | 7.50 | 6.95 | - | - | -
$\mathrm{GW}190828\\_063405$ | $1251009263.8$ | 12390.23 | 8.50 | 9.64 | 210.00 | 9.11 | 10.23
$\mathrm{GW}190915\\_235702$ | $1252627040.7$ | 16.83 | 7.49 | 6.88 | 9.80 | 8.87 | 9.66
$\mathrm{GW}191109\\_010717$ | $1257296855.2$ | 5.43 | 8.52 | 12.25 | - | - | -
$\mathrm{GW}200128\\_022011$ | $1264213229.9$ | 2.44 | 6.41 | 6.32 | - | - | -
$\mathrm{GW}200224\\_222234$ | $1266618172.4$ | 1769.27 | 10.59 | 10.77 | 106.05 | 10.90 | 10.56
$\mathrm{GW}200225\\_060421$ | $1266645879.4$ | 6.08 | 7.63 | 6.09 | 59.95 | 7.85 | 6.38
$\mathrm{GW}200311\\_115853$ | $1267963151.3$ | 2501.38 | 9.70 | 8.94 | 533.38 | 9.85 | 9.25
Table 2: gravitational wave events recovered by our search with inverse false
alarm rate $>1$ year, calculated using O3 data from Hanford and Livingston
detector. These events have all been reported in Abbott  _et al._ (2021a, b);
Nitz  _et al._ (2021a, b); Chandra  _et al._ (2021a); Olsen  _et al._ (2022)
with a comparatively higher statistical significance. We attribute their
comparatively lower significance to their lack of detectable higher-order
modes and/or being outside the bounds of our target search space.

##  V  Results from the third LIGO observing run

Given the sensitivity improvement to massive binary black hole signals with
higher harmonics, we deployed our new search to analyze the publicly available
O3 data from Advanced LIGO detectors (Abbott  _et al._ , 2021c). To do this,
we divided the data into nine independent blocks of $\sim 30$-day duration and
analyzed it using our template banks. We report all the O3 candidates reported
by our searches with an inverse false alarm rate $>1$ year in Table 2, finding
no new candidates beyond those reported in Abbott  _et al._ (2021a, b); Nitz
_et al._ (2021a, b); Olsen  _et al._ (2022). We note that the (expected)
reduced statistical significance of the reported events with respect to those
obtained by existing searches highlights that these candidates are outside our
target space.

We note that the lack of new detections is not unexpected. Firstly, signals
from binaries with the inclinations targeted by our banks are significantly
weaker than those from face-on binaries typically detected by existing
searches. Secondly, in contrast to existing searches, we have performed a two-
detector search using only Advanced LIGO Hanford and Advanced LIGO Livingston
data, reducing the loudness of the signals across the detector network. We
have checked that these two effects lead to an optimally observable volume
$\simeq 7$ times smaller than that of search targeting rather than face-on
signals. The fact that existing searches have only reported three events to
date within our targeted mass range (GW190521, GW200220, and GW190426) Abbott
_et al._ (2021b) is perfectly consistent with our results. In addition, the
inclusion of a third – Virgo – detector in our search will also improve our
ability to discriminate glitches from true signals, further increasing our
sensitivity. We leave the inclusion of further detectors in our search for
further improvement. Such work will primarily involve extending the two-
detector ranking statistic to a multi-detector one.

##  VI Conclusions

We have presented the first matched-filter gravitational-wave search for
compact mergers, including the impact of higher gravitational-wave modes. On
the one hand, detecting merger-ringdown higher modes is crucial to enabling
tests of General Relativity in the strong-field regime and is key to observing
several strong-field phenomena arising during the merger-ringdown stages. On
the other hand, due to the high mass of these systems, the emission from
intermediate-mass black-hole mergers contains strong higher-mode contributions
for orientations other than face-on/off. This makes this type of search
crucial for both fundamental physics and astrophysics.

While one can extend our search to signals from any binary merger, we have
specifically restricted to those with strong higher-mode content: asymmetric,
non-precessing, massive black-hole binaries at large orbital inclinations,
leaving the extension to more generic inclinations as future work. While we
find no new candidates with IFAR $>1$ year beyond those already reported
elsewhere, we have demonstrated that our search is up to 450% more sensitive
than existing matched-filter searches with overlapping parameters space but
omitting the higher-mode content. We stress that the lack of new, higher-mode
rich detections is somewhat expected, as the highly inclined systems we target
approximately span an observable volume $\simeq 7$ times smaller than the
rather face-on ones targeted by existing searches, which have detected three
events in our target mass-range to date.

Also, the expected merger rate density of intermediate-mass black hole
binaries is $<1/\mathrm{Gpc}^{3}/\mathrm{yr}$ as compared to $\sim
23/\mathrm{Gpc}^{3}/\mathrm{yr}$ for stellar-mass black hole binaries (Abbott
_et al._ , 2022, 2021d). This means that we expected to observe $\sim 100$
stellar-mass black holes as compared to $\sim 3$ intermediate-mass black hole
binaries during O3 given that the achieved sensitivity allowed us to observe
$30M_{\odot}+30M_{\odot}$ and $100M_{\odot}+100M_{\odot}$ to a distance of
$\sim 1$Gpc and $\sim 2$Gpc respectively. As the detector’s distance reach
improves, we are likely to observe more intermediate-mass black hole binaries,
and thus, we are more likely to observe signals from our target population in
the future (Abbott  _et al._ , 2018).

Finally, we note that our existing search needs to undergo several
improvements before being equivalently mature with respect to existing ones.
First, our search is restricted to data from only two detectors. Extending
this to an arbitrary number of detectors would greatly help better
discriminate true signals from glitches. Second, there is room to improve how
we exploit the information contained in higher harmonics to remove more
background triggers. Third, our search at the moment is meant to complement
existing searches with overlapping target parameter space. A near-term
development goal is to devise a strategy that combines the search output of
the contributing searches so that we quote a single quantitative estimate of
an event’s statistical significance. Last, as it is also the case for other
existing searches, an obvious extension would be expanding the target search
space to include orbital precession or eccentricity. The first is
straightforward to implement, given the genericness of our search method,
while the others require a combination of stronger vetoes (Dhurandhar  _et
al._ , 2017; Jadhav  _et al._ , 2021; McIsaac and Harry, 2022; Davis  _et al._
, 2022), glitch-robust search statistics and glitch subtraction Cornish  _et
al._ (2021), that can exploit the complex waveform morphology due to addition
of new physics.

## Acknowledgements

We thank Tito Dal Canton, Thomas Dent, and Kritti Sharma for their detailed
comments and valuable suggestions. We would also like to thank Verónica Villa-
Ortega for running an earlier version of the presented algorithm. This
material is based upon work supported by NSF’s LIGO Laboratory, a major
facility fully funded by the National Science Foundation. This research has
used data or software from the Gravitational Wave Open Science Center (gw-
openscience.org), a service of LIGO Laboratory, the LIGO Scientific
Collaboration, the Virgo Collaboration, and KAGRA. LIGO Laboratory and
Advanced LIGO are funded by the United States NSF as well as the STFC of the
United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. The Australian Research
Council provided additional support for Advanced LIGO. The authors are
grateful for the computational resources and data provided by the LIGO
Laboratory and supported by National Science Foundation Grants No. PHY-0757058
and No. PHY-0823459. The authors also acknowledge using the IUCAA LDG cluster,
Sarathi, for computational/numerical work. KC acknowledges the MHRD, the
Government of India, for the fellowship support. JCB is supported by a
fellowship from the “la Caixa” Foundation (ID100010434) and from the European
Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie grant agreement No 847648. The fellowship code is
LCF/BQ/PI20/11760016. JCB is also supported by the research grant
PID2020-118635GB-I00 from the Spain-Ministerio de Ciencia e Innovación. AP’s
research is supported by SERB-Power fellowship grant SPF/2021/000036, DST,
India. IWH acknowledge the STFC for funding through grant ST/T000333/1 and
ST/V005715/1. This document has LIGO DCC No LIGO-P2200182.

## References

  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2111.03606 [gr-qc] .
  * Rodriguez  _et al._ (2020) C. L. Rodriguez _et al._ , Astrophys. J. Lett. 896, L10 (2020), arXiv:2005.04239 [astro-ph.HE] .
  * Zevin  _et al._ (2021) M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez, F. Antonini, D. E. Holz, and C. Pankow, Astrophys. J. 910, 152 (2021), arXiv:2011.10057 [astro-ph.HE] .
  * Fragione  _et al._ (2020) G. Fragione, A. Loeb, and F. A. Rasio, Astrophys. J. Lett. 902, L26 (2020), arXiv:2009.05065 [astro-ph.GA] .
  * Kimball  _et al._ (2021) C. Kimball _et al._ , Astrophys. J. Lett. 915, L35 (2021), arXiv:2011.05332 [astro-ph.HE] .
  * Baibhav  _et al._ (2021) V. Baibhav, E. Berti, D. Gerosa, M. Mould, and K. W. K. Wong, Phys. Rev. D 104, 084002 (2021), arXiv:2105.12140 [gr-qc] .
  * Gerosa and Fishbach (2021) D. Gerosa and M. Fishbach, Nature Astron. 5, 749 (2021), arXiv:2105.03439 [astro-ph.HE] .
  * Tagawa  _et al._ (2021) H. Tagawa, Z. Haiman, I. Bartos, B. Kocsis, and K. Omukai, Mon. Not. Roy. Astron. Soc. 507, 3362 (2021), arXiv:2104.09510 [astro-ph.HE] .
  * Fishbach  _et al._ (2017) M. Fishbach, D. E. Holz, and B. Farr, Astrophys. J. Lett. 840, L24 (2017), arXiv:1703.06869 [astro-ph.HE] .
  * Gerosa and Berti (2017) D. Gerosa and E. Berti, Phys. Rev. D 95, 124046 (2017), arXiv:1703.06223 [gr-qc] .
  * Rodriguez  _et al._ (2019) C. L. Rodriguez, M. Zevin, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio, and C. S. Ye, Phys. Rev. D 100, 043027 (2019), arXiv:1906.10260 [astro-ph.HE] .
  * Kimball  _et al._ (2020) C. Kimball, C. Talbot, C. P. L. Berry, M. Carney, M. Zevin, E. Thrane, and V. Kalogera, Astrophys. J. 900, 177 (2020), arXiv:2005.00023 [astro-ph.HE] .
  * Doctor  _et al._ (2021) Z. Doctor, B. Farr, and D. E. Holz, Astrophys. J. Lett. 914, L18 (2021), arXiv:2103.04001 [astro-ph.HE] .
  * Mapelli  _et al._ (2021) M. Mapelli _et al._ , Mon. Not. Roy. Astron. Soc. 505, 339 (2021), arXiv:2103.05016 [astro-ph.HE] .
  * Fragione  _et al._ (2022) G. Fragione, B. Kocsis, F. A. Rasio, and J. Silk, Astrophys. J. 927, 231 (2022), arXiv:2107.04639 [astro-ph.GA] .
  * Quinlan and Shapiro (1987) G. D. Quinlan and S. L. Shapiro (1987).
  * Greene  _et al._ (2020) J. E. Greene, J. Strader, and L. C. Ho, Annual Review of Astronomy and Astrophysics 58, 257 (2020), arXiv:1911.09678 [astro-ph.GA] .
  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 221101 (2016a), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Carullo  _et al._ (2019) G. Carullo, W. Del Pozzo, and J. Veitch, Phys. Rev. D 99, 123029 (2019), [Erratum: Phys.Rev.D 100, 089903 (2019)], arXiv:1902.07527 [gr-qc] .
  * Isi  _et al._ (2019) M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. Lett. 123, 111102 (2019), arXiv:1905.00869 [gr-qc] .
  * Dreyer  _et al._ (2004) O. Dreyer, B. J. Kelly, B. Krishnan, L. S. Finn, D. Garrison, and R. Lopez-Aleman, Class. Quant. Grav. 21, 787 (2004), arXiv:gr-qc/0309007 .
  * Berti  _et al._ (2016) E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, Phys. Rev. Lett. 117, 101102 (2016), arXiv:1605.09286 [gr-qc] .
  * Bhagwat  _et al._ (2018) S. Bhagwat, M. Okounkova, S. W. Ballmer, D. A. Brown, M. Giesler, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 97, 104065 (2018), arXiv:1711.00926 [gr-qc] .
  * Cabero  _et al._ (2020) M. Cabero, J. Westerweck, C. D. Capano, S. Kumar, A. B. Nielsen, and B. Krishnan, Phys. Rev. D 101, 064044 (2020), arXiv:1911.01361 [gr-qc] .
  * Bustillo  _et al._ (2021) J. C. Bustillo, P. D. Lasky, and E. Thrane, Phys. Rev. D 103, 024041 (2021), arXiv:2010.01857 [gr-qc] .
  * Calderón Bustillo  _et al._ (2018a) J. Calderón Bustillo, J. A. Clark, P. Laguna, and D. Shoemaker, Phys. Rev. Lett. 121, 191102 (2018a), arXiv:1806.11160 [gr-qc] .
  * Calderón Bustillo  _et al._ (2019) J. Calderón Bustillo, C. Evans, J. A. Clark, G. Kim, P. Laguna, and D. Shoemaker, Communications Physics 3 (2019), 10.1038/s42005-020-00446-7, arXiv:1906.01153 [gr-qc] .
  * Varma  _et al._ (2020) V. Varma, M. Isi, and S. Biscoveanu, Phys. Rev. Lett. 124, 101104 (2020), arXiv:2002.00296 [gr-qc] .
  * Calderón Bustillo  _et al._ (2022) J. Calderón Bustillo, S. H. W. Leong, and K. Chandra, (2022), arXiv:2211.03465 [gr-qc] .
  * Santamaria  _et al._ (2010) L. Santamaria _et al._ , Phys. Rev. D82, 064016 (2010), arXiv:1005.3306 [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] .
  * Husa  _et al._ (2016) S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D 93, 044006 (2016), arXiv:1508.07250 [gr-qc] .
  * Buonanno and Damour (2000) A. Buonanno and T. Damour, Phys.Rev. D62, 064015 (2000), arXiv:gr-qc/0001013 [gr-qc] .
  * Damour  _et al._ (2013) T. Damour, A. Nagar, and S. Bernuzzi, Phys.Rev. D87, 084035 (2013), arXiv:1212.4357 [gr-qc] .
  * Buonanno and Damour (1999) A. Buonanno and T. Damour, Phys. Rev. D59, 084006 (1999), arXiv:gr-qc/9811091 [gr-qc] .
  * Cotesta  _et al._ (2018) R. Cotesta, A. Buonanno, A. Bohé, A. Taracchini, I. Hinder, and S. Ossokine, Phys. Rev. D 98, 084028 (2018), arXiv:1803.10701 [gr-qc] .
  * Szilagyi  _et al._ (2009) B. Szilagyi, L. Lindblom, and M. A. Scheel, Phys.Rev. D80, 124010 (2009), arXiv:0909.3557 [gr-qc] .
  * Jani  _et al._ (2016) K. Jani, J. Healy, J. A. Clark, L. London, P. Laguna, and D. Shoemaker, Class. Quant. Grav. 33, 204001 (2016), arXiv:1605.03204 [gr-qc] .
  * Healy  _et al._ (2017) J. Healy, C. O. Lousto, Y. Zlochower, and M. Campanelli, Classical and Quantum Gravity 34, 224001 (2017), arXiv:1703.03423 [gr-qc] .
  * Wainstein and Zubakov (1962) L. A. Wainstein and V. D. Zubakov,  _Extraction of Signals from Noise_ , Dover books on physics and mathematical physics (Prentice-Hall, Englewood Cliffs, NJ, 1962).
  * Sathyaprakash and Dhurandhar (1991) B. S. Sathyaprakash and S. V. Dhurandhar, Phys. Rev. D 44, 3819 (1991).
  * Allen  _et al._ (2012) B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E. Creighton, Phys. Rev. D 85, 122006 (2012), arXiv:gr-qc/0509116 .
  * Usman  _et al._ (2016) S. A. Usman _et al._ , Class. Quant. Grav. 33, 215004 (2016), arXiv:1508.02357 [gr-qc] .
  * Venumadhav  _et al._ (2019) T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, Phys. Rev. D 100, 023011 (2019), arXiv:1902.10341 [astro-ph.IM] .
  * Sachdev  _et al._ (2019) S. Sachdev _et al._ , (2019), arXiv:1901.08580 [gr-qc] .
  * Aubin  _et al._ (2021) F. Aubin _et al._ , Class. Quant. Grav. 38, 095004 (2021), arXiv:2012.11512 [gr-qc] .
  * Chu  _et al._ (2022) Q. Chu _et al._ , Phys. Rev. D 105, 024023 (2022), arXiv:2011.06787 [gr-qc] .
  * Mills and Fairhurst (2021) C. Mills and S. Fairhurst, Phys. Rev. D 103, 024042 (2021), arXiv:2007.04313 [gr-qc] .
  * Capano  _et al._ (2014) C. Capano, Y. Pan, and A. Buonanno, Phys. Rev. D 89, 102003 (2014), arXiv:1311.1286 [gr-qc] .
  * Calderón Bustillo  _et al._ (2016) J. Calderón Bustillo, S. Husa, A. M. Sintes, and M. Pürrer, Phys. Rev. D 93, 084019 (2016), arXiv:1511.02060 [gr-qc] .
  * Calderón Bustillo  _et al._ (2017) J. Calderón Bustillo, P. Laguna, and D. Shoemaker, Phys. Rev. D 95, 104038 (2017), arXiv:1612.02340 [gr-qc] .
  * Calderón Bustillo  _et al._ (2018b) J. Calderón Bustillo, F. Salemi, T. Dal Canton, and K. P. Jani, Phys. Rev. D 97, 024016 (2018b), arXiv:1711.02009 [gr-qc] .
  * Chandra  _et al._ (2020) K. Chandra, V. Gayathri, J. C. Bustillo, and A. Pai, Phys. Rev. D 102, 044035 (2020), arXiv:2002.10666 [astro-ph.CO] .
  * Harry  _et al._ (2018) I. Harry, J. Calderón Bustillo, and A. Nitz, Phys. Rev. D 97, 023004 (2018), arXiv:1709.09181 [gr-qc] .
  * Chandra  _et al._ (2021a) K. Chandra, V. Villa-Ortega, T. Dent, C. McIsaac, A. Pai, I. W. Harry, G. S. C. Davies, and K. Soni, Phys. Rev. D 104, 042004 (2021a), arXiv:2106.00193 [gr-qc] .
  * Nitz  _et al._ (2021a) A. H. Nitz, C. D. Capano, S. Kumar, Y.-F. Wang, S. Kastha, M. Schäfer, R. Dhurkunde, and M. Cabero, Astrophys. J. 922, 76 (2021a), arXiv:2105.09151 [astro-ph.HE] .
  * Nitz  _et al._ (2021b) A. H. Nitz, S. Kumar, Y.-F. Wang, S. Kastha, S. Wu, M. Schäfer, R. Dhurkunde, and C. D. Capano, (2021b), arXiv:2112.06878 [astro-ph.HE] .
  * Olsen  _et al._ (2022) S. Olsen, T. Venumadhav, J. Mushkin, J. Roulet, B. Zackay, and M. Zaldarriaga, Phys. Rev. D 106, 043009 (2022), arXiv:2201.02252 [astro-ph.HE] .
  * 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).
  * Blanchet (2014) L. Blanchet, Living Rev. Rel. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Berti  _et al._ (2007) E. Berti, V. Cardoso, J. A. Gonzalez, U. Sperhake, M. Hannam, S. Husa, and B. Bruegmann, Phys. Rev. D 76, 064034 (2007), arXiv:gr-qc/0703053 .
  * Cotesta  _et al._ (2020) R. Cotesta, S. Marsat, and M. Pürrer, Phys. Rev. D 101, 124040 (2020), arXiv:2003.12079 [gr-qc] .
  * Cannon  _et al._ (2021) K. Cannon, S. Caudill, C. Chan, B. Cousins, J. D. Creighton, B. Ewing, H. Fong, P. Godwin, C. Hanna, S. Hooper, R. Huxford, R. Magee, D. Meacher, C. Messick, S. Morisaki, D. Mukherjee, H. Ohta, A. Pace, S. Privitera, I. de Ruiter, S. Sachdev, L. Singer, D. Singh, R. Tapia, L. Tsukada, D. Tsuna, T. Tsutsui, K. Ueno, A. Viets, L. Wade, and M. Wade, SoftwareX 14, 100680 (2021).
  * Abbott  _et al._ (2022) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Astron. Astrophys. 659, A84 (2022), arXiv:2105.15120 [astro-ph.HE] .
  * Klimenko  _et al._ (2016) S. Klimenko _et al._ , Phys. Rev. D 93, 042004 (2016), arXiv:1511.05999 [gr-qc] .
  * Szczepańczyk  _et al._ (2021) M. Szczepańczyk _et al._ , Phys. Rev. D 103, 082002 (2021), arXiv:2009.11336 [astro-ph.HE] .
  * Chandra  _et al._ (2021b) K. Chandra, A. Pai, V. Villa-Ortega, T. Dent, C. McIsaac, I. W. Harry, G. S. C. Davies, and K. Soni, in  _16th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics and Relativistic Field Theories_ (2021) arXiv:2110.01879 [gr-qc] .
  * Harry  _et al._ (2008) I. W. Harry, S. Fairhurst, and B. S. Sathyaprakash, Class. Quant. Grav. 25, 184027 (2008), arXiv:0804.3274 [gr-qc] .
  * Babak (2008) S. Babak, Class. Quant. Grav. 25, 195011 (2008), arXiv:0801.4070 [gr-qc] .
  * Harry  _et al._ (2009) I. W. Harry, B. Allen, and B. S. Sathyaprakash, Phys. Rev. D 80, 104014 (2009), arXiv:0908.2090 [gr-qc] .
  * Apostolatos  _et al._ (1994) T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D 49, 6274 (1994).
  * Allen (2005) B. Allen, Phys. Rev. D 71, 062001 (2005), arXiv:gr-qc/0405045 .
  * Nitz (2018) A. H. Nitz, Class. Quant. Grav. 35, 035016 (2018), arXiv:1709.08974 [gr-qc] .
  * Mozzon  _et al._ (2020) S. Mozzon, L. K. Nuttall, A. Lundgren, T. Dent, S. Kumar, and A. H. Nitz, Class. Quant. Grav. 37, 215014 (2020), arXiv:2002.09407 [astro-ph.IM] .
  * Nitz  _et al._ (2017) A. H. Nitz, T. Dent, T. Dal Canton, S. Fairhurst, and D. A. Brown, Astrophys. J. 849, 118 (2017), arXiv:1705.01513 [gr-qc] .
  * Davies  _et al._ (2020) G. S. Davies, T. Dent, M. Tápai, I. Harry, C. McIsaac, and A. H. Nitz, Phys. Rev. D 102, 022004 (2020), arXiv:2002.08291 [astro-ph.HE] .
  * Babak  _et al._ (2013) S. Babak _et al._ , Phys. Rev. D 87, 024033 (2013), arXiv:1208.3491 [gr-qc] .
  * Pai  _et al._ (2001) A. Pai, S. Dhurandhar, and S. Bose, Phys. Rev. D 64, 042004 (2001), arXiv:gr-qc/0009078 .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016b), arXiv:1602.03837 [gr-qc] .
  * Capano  _et al._ (2016) C. Capano, I. Harry, S. Privitera, and A. Buonanno, Phys. Rev. D 93, 124007 (2016), arXiv:1602.03509 [gr-qc] .
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, Virgo), SoftwareX 13, 100658 (2021c), arXiv:1912.11716 [gr-qc] .
  * Abbott  _et al._ (2021d) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 913, L7 (2021d), arXiv:2010.14533 [astro-ph.HE] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (KAGRA, LIGO Scientific, Virgo, VIRGO), Living Rev. Rel. 21, 3 (2018), arXiv:1304.0670 [gr-qc] .
  * Dhurandhar  _et al._ (2017) S. Dhurandhar, A. Gupta, B. Gadre, and S. Bose, Phys. Rev. D 96, 103018 (2017), arXiv:1708.03605 [gr-qc] .
  * Jadhav  _et al._ (2021) S. Jadhav, N. Mukund, B. Gadre, S. Mitra, and S. Abraham, Phys. Rev. D 104, 064051 (2021), arXiv:2010.08584 [gr-qc] .
  * McIsaac and Harry (2022) C. McIsaac and I. Harry, Phys. Rev. D 105, 104056 (2022), arXiv:2203.03449 [gr-qc] .
  * Davis  _et al._ (2022) D. Davis, M. Trevor, S. Mozzon, and L. K. Nuttall, (2022), arXiv:2204.03091 [gr-qc] .
  * Cornish  _et al._ (2021) N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse, Phys. Rev. D 103, 044006 (2021), arXiv:2011.09494 [gr-qc] .

  *[O3]: third observing run
  *[SNR]: signal-to-noise ratio
  *[FF]: fitting factor
  *[SNRs]: signal-to-noise ratio
  *[IFARs]: inverse false alarm rate
  *[IFAR]: inverse false alarm rate