# A novel signal-consistency test for gravitational-wave searches of generic
black hole binaries

Stefano Schmidt  0000-0002-8206-8089<EMAIL_ADDRESS>Nikhef, Science Park 105,
1098 XG, Amsterdam, The Netherlands  Institute for Gravitational and Subatomic
Physics (GRASP), Utrecht University, Princetonplein 1, 3584 CC Utrecht, The
Netherlands  Sarah Caudill  0000-0002-8927-6673 Department of Physics,
University of Massachusetts, Dartmouth, MA 02747, USA  Center for Scientific
Computing and Data Science Research, University of Massachusetts, Dartmouth,
MA 02747, USA

###### Abstract

We propose a novel signal-consistency test applicable to a broad search for
gravitational waves emitted by generic binary black hole (BBH) systems. The
test generalizes the time domain $\xi^{2}$ signal-consistency test currently
utilized by the GstLAL pipeline, which is restricted to aligned-spin circular
orbits and does not account for higher-order modes (HMs). After addressing the
mathematical details of the new test, we quantify its advantages in the
context of searching for precessing BBHs and/or BBHs with HM content. Our
results reveal that for precessing signals, the new test has the potential to
diminish $\xi^{2}$ by up to two orders of magnitude when compared to the
standard test. However, in the case of signals with HM content, only a modest
enhancement is observed. Recognizing the computational burden associated with
the new test, we also derive an approximated signal-consistency test. This
approximation maintains the same computational cost as the standard test while
sacrificing only a few percent of accuracy in the low SNR regime. The
approximated test can be easily implemented in any matched filtering pipeline
aimed at searching for BBH signals with strong precession content.

##  I Introduction

Gravitational-wave searches for binary black holes (BBHs) have been a
cornerstone in the activities of the LIGO-Virgo-KAGRA (LVK) collaboration [1,
2, 3] and have made possible the discovery of nearly one hundred GW events
during the first three observing runs [4, 5, 6, 7].

BBH searches can either be performed in a model-dependent way through the
technique of matched filtering [8, 9, 10, 11, 12] or a model-independent way
through excess-power methods [13, 14, 15]. The matched-filter method computes
the correlation between GW detector data and a set of templates of modelled
BBH waveforms. This method proves more sensitive for lower-mass systems, as
the power in these signals is spread over many time-frequency bins [16, 17].
Beyond the initial matched-filter or excess-power stage, searches consist of
analysis pipelines [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29] that
employ various statistical tests to improve the separation of GW candidates
from false alarms caused by the detectors’ non-stationary and non-Gaussian
noise [30, 31, 32, 33, 34].

A class of these statistical tests, known as signal-consistency tests [35, 36,
37, 22, 38, 39, 40, 25, 41, 42, 43], have been designed to further distinguish
between false alarms and real GW signals by computing the agreement between
the data and the signal model assumed by the search. Signal-consistency tests
have been implemented by all LVK matched-filtering pipelines and have played
an integral role in enabling many high-significance GW detections.

Expanding our searches to include signals from precessing binaries [44, 45,
46, 47, 48, 49] and/or binaries with higher-order modes (HMs) [50, 51, 52, 53,
54] necessitates the generalization of signal-consistency tests to a more
versatile framework. Neglecting to update these tests can result in decreased
search sensitivity, potentially offsetting the benefits of using a more
diverse set of templates.

While the $\chi^{2}$ time-frequency signal-consistency test [35] and its
variant, the sine-Gaussian $\chi^{2}$ discriminator [38], have been
successfully applied in searches including higher-order modes [51, 52] and
precessing signals [49], little attention has been given to generalizing other
types of signal-consistency tests, or even to the development of new ones.

In this paper, we introduce a new signal-consistency test that is specifically
designed for precessing and/or HM signals. It is derived from the
autocorrelation-based least-squares test, denoted $\xi^{2}$, which is
currently utilized by the GstLAL search pipeline [22, 24, 55, 29]. In Sec. II,
we provide some general background on matched filtering and on the state of
the art signal-consistency test, while in Sec. III we introduce our new
generalized signal-consistency test $\xi_{\text{sym}}^{2}$ and its
computationally convenient approximate expression $\xi_{\text{mix}}^{2}$. In
Sec. IV we discuss the range of applicability and performance of the newly
introduced test and its approximated version. Sec. V gathers some final
remarks and future outlook.

##  II Background

According to the theory of general relativity [56], a gravitational wave only
carries two physical degrees of freedom $h_{+}$ and $h_{\times}$, also called
polarizations. For a generic GW emitted by a compact binary, the polarizations
depend on the intrinsic properties of the binary such as the two compact
objects’ masses $m_{1}$, $m_{2}$ and spins $\mathbf{s}_{1}$, $\mathbf{s}_{2}$.
The signal observed at the source may also depend on the so-called extrinsic
properties including position, usually parametrized in spherical coordinates
by a distance $D$ from the origin, a polar angle called inclination angle
$\iota$, and an azimuthal angle $\varphi$.

For the purpose of modeling, it is customary to expand the gravitational
waveform’s dependence on the extrinsic parameters in terms spin-2 spherical
harmonics $Y_{\ell m}(\iota,\varphi)$ [56]

$h_{+}+ih_{\times}=\frac{1}{D}\sum^{\infty}_{\ell=0}\sum^{\ell}_{m=-\ell}Y_{\ell
m}(\iota,\varphi)e^{im\varphi}h_{\ell m}(t)$ (1)

where the modes $h_{\ell m}$ are complex functions of the intrinsic parameters
$m_{1},m_{2},\mathbf{s}_{1},\mathbf{s}_{2}$. Each mode can be decomposed into
a time-dependent amplitude $A_{\ell m}$ and a time-dependent phase $\phi_{\ell
m}$ such that:

$h_{\ell m}=A_{\ell m}e^{i\phi_{\ell m}}.$ (2)

As a consequence the imaginary part $h^{\text{I}}_{\ell m}$ of each mode is
equivalent to the real part $h^{\text{R}}_{\ell m}$ shifted by a constant
phase of $\pi/2$. In the frequency domain, this takes the simple form of

$\tilde{h}^{\text{R}}_{\ell m}=i\tilde{h}^{\text{I}}_{\ell m}$ (3)

where $\tilde{\phantom{a}}$ indicates the Fourier transform.

If the two spins $\mathbf{s}_{1},\mathbf{s}_{2}$ are misaligned with the
binary orbital angular momentum $\mathbf{L}$, the binary plane experiences
precessional motion, where the orbital angular momentum rotates around a
(roughly) constant direction [57, 58, 59, 60, 61]. If the two spins are
aligned with $\mathbf{L}$, the orbital plane will point toward a fixed
direction, establishing an axis of symmetry for the non-precessing binary
system. Mathematically, this symmetry translates into a symmetry between the
modes:

$h_{\ell m}=(-1)^{\ell}h^{*}_{\ell-m}$ (4)

where $*$ denotes complex conjugation.

In most aligned-spins systems, it turns out that only the $\ell=|m|=2$ modes
give a significant contribution to the polarizations and thus all the other
HMs are neglected [56, 62, 63, 50]. Thus, the waveform from an non-precessing
binary without imprints from HMs has a strikingly simple expression:

$\displaystyle h_{+}$
$\displaystyle=\frac{1}{D}\frac{1+\cos^{2}\iota}{2}\;\Re\\{h_{22}e^{i2\varphi}\\},$
(5) $\displaystyle h_{\times}$
$\displaystyle=\frac{1}{D}\cos\iota\;\Im\\{h_{22}e^{i2\varphi}\\}.$ (6)

Note that as a consequence of Eq. (3), the two polarizations in Fourier space
are related by the simple relation $\tilde{h}_{+}\propto i\tilde{h}_{\times}$
and, for $\iota=0$, we have trivially $\tilde{h}_{+}=i\tilde{h}_{\times}$.

In the remainder of this section, we describe how the simplicity of Eqs. (5)
and (6) enters the the expression for the currently used $\xi^{2}$ signal-
consistency test [22] and we describe how to move away from the assumption of
aligned-spin systems without HMs, tackling the most general case.

###  II.1 Overview

The core of matched filter relies on computing the correlation between two
time-series $a(t),b(t)$, weighted by the Power Spectral Density (PSD) of the
noise $S_{n}(f)$. Mathematically, we can define a time-dependent complex
scalar product:

$\langle
a|b\rangle(t)=4\int_{0}^{\infty}\operatorname{d}\\!{f}\;\frac{\tilde{a}^{*}(f)\tilde{b}(f)}{S_{n}(f)}e^{i2\pi
ft}$ (7)

It is convenient to define the real $(\cdot|\cdot)$ and imaginary
$[\cdot|\cdot]$ part of the scalar product as:

$\langle a|b\rangle(t)=(a|b)(t)+i[a|b](t)$ (8)

Given a timeseries $a(t)$, we can use the above scalar product to computed the
normalized timeseries $\hat{a}(t)$:

$\hat{a}(t)=\frac{a(t)}{\sqrt{(a|a)(t=0)}}$ (9)

Following the notation of [22], the output of the matched-filtering procedure
is a complex timeseries $z(t)$:

$z(t)=(d|h_{R})(t)+i(d|h_{I})(t)$ (10)

where $h_{R}$ and $h_{I}$ are the real and imaginary part of a normalized
complex template. The actual form of the template depends on the assumptions
about the nature of the GW signal to detect.

We may also define the Signal-to-noise ratio (SNR) timeseries $\rho(t)$ as:

$\rho(t)=|z(t)|=\sqrt{(d|h_{R})^{2}(t)+(d|h_{I})^{2}(t)}$ (11)

This is the primary quantity of interest of a GW search.

Given a trigger at time $t=0$, the $\xi^{2}$ test relies on predicting the SNR
timeseries $z(t)$ obtained by filtering a signal $h$ with a matching
templates. The predicted timeseries $R(t)$ is then compared to the measured
timeseries $z(t)$ to compute the $\xi^{2}$ quantity:

$\xi^{2}=\frac{\int_{-\delta t}^{\delta
t}\;\operatorname{d}\\!{t}\;\;|z(t)-R(t)|^{2}}{\int_{-\delta t}^{\delta
t}\;\operatorname{d}\\!{t}\;\;\left(2-2\left|\frac{R(t)}{R(0)}\right|^{2}\right)}$
(12)

where the integral extends on a short time window $[-\delta t,\delta t]$
around the trigger time. It is convenient to express $\delta t$ in terms of
the so-called autocorrelation length (ACL) [22], defined as the number of
samples in the time-window $[-\delta t,\delta t]$ at a given sample rate
$f_{\text{sampling}}$, so that $\delta
t=(\textrm{ACL}-1)/2\,f_{\text{sampling}}$.

As shown in [22], the denominator in Eq. (12) amounts to the expected value of
$|z(t)-R(t)|^{2}$ in gaussian noise. While the expression refers to a standard
search (aligned-spin/non-HM), for simplicity we keep the same expression even
for the general case. Therefore for the general case, it must be intended as a
mere normalization and loses the physical interpretation that it has in the
standard case.

The $\xi^{2}$ defined above can be used by the GW search pipelines to veto
some loud triggers. If a trigger is caused by a noise fluctuation or non-
gaussian noise transient bursts [30], the discrepancy between the expected and
measured SNR timeseries will be large, leading to a large value of $\xi^{2}$.
This can be used to downrank certain triggers, with large improvement in
sensitivity.

As it is custom, to predicted the SNR timeseries we model the data $d$ as a
superposition of gaussian noise $n$ and a GW signal $h$: $d=n+h$. For current
ground based interferometers, the GW signal is:

$h=F_{+}h_{+}+F_{\times}h_{\times}$ (13)

where $F_{+},F_{\times}$ are the Antenna Pattern Functions [64, 65], which
define the detector’s response to a signal coming from a given direction in
the sky. For our purpose, it is convenient to express the signal model in
terms of the normalized polarizations:

$h=\mathcal{F}_{+}\hat{h}_{+}+\mathcal{F}_{\times}\hat{h}_{\times}$ (14)

where we absorbed into $\mathcal{F}_{+},\mathcal{F}_{\times}$ an overall
scaling factor (which depends on the source distance and on the sky location).
Of course, $\mathcal{F}_{+},\mathcal{F}_{\times}$ are not known at the moment
of the search but must be inferred from the value of the SNR timeseries $z(0)$
at the time of a trigger.

Figure 1:  Predicted and measured SNR timeseries, $R$ and $\rho$ respectively,
for a precessing injection at $SNR=20$ filtered with a perfectly matching
template. We measure the SNR using both the standard SNR Eq. (21) and the
“symphony” SNR Eq. (25). The predicted SNR is computed with two different
prescriptions, $R_{\text{std}}$ Eq. (22) suitable for the standard SNR and
$R_{\text{sym}}$ Eq. (26) suitable for the “symphony” SNR. On the bottom
panel, we note that $R_{\text{sym}}$ and $\rho_{\text{sym}}$ show excellent
agreement between each other while $R_{\text{std}}$ and $\rho_{\text{std}}$
show mild agreement (top panel). This is quantified by the values of
$\xi_{\text{std}}^{2}$, $\xi_{\text{mix}}^{2}$ and $\xi_{\text{sym}}^{2}$
reported on the plot. The signal is injected into gaussian noise, sampled from
the PSD LIGO Livingston PSD [66] with a rate of
$4096\text{\,}\mathrm{H}\mathrm{z}$. The waveform is characterized by masses
$m_{1},m_{2}=$28\text{\,}\mathrm{M}_{\odot}$,$3\text{\,}\mathrm{M}_{\odot}$$
and spins $\mathbf{s}_{1}=(-0.8,0.02,-0.5)$ and $\mathbf{s}_{2}=0$, observed
with an inclination $\iota=2.66$. It was generated starting from a frequency
of $10\text{\,}\mathrm{H}\mathrm{z}$ with the approximant IMRPhenomXP [67].

###  II.2 The original $\xi^{2}$ test

As discussed above, for the case of a spin-aligned system where HM are not
considered, there is a well known symmetry between the two polarizations:
$\tilde{h}_{+}\propto i\tilde{h}_{\times}$. As a consequence [48], the complex
template is simply given by $\hat{h}_{22}$ and the detection statistics does
not depend on $\iota$ and $\varphi$:

$\displaystyle h_{R}=\hat{h}^{\text{R}}_{22}$ (15) $\displaystyle
h_{I}=\hat{h}^{\text{I}}_{22}$

We note that the template $\hat{h}^{\text{R}}_{22}$ is equivalent to the
normalized plus polarization, evaluated at zero reference phase. Moreover,
from Eq. (6) we see that $\hat{h}^{\text{I}}_{22}$ is equal to the normalized
cross polarization up to a minus sign, depending on the value of inclination
being considered: $\hat{h}^{\text{I}}_{22}=\pm\hat{h}_{\times}$. We call
“standard” the SNR timeseries obtained using such templates:

$\rho_{\text{std}}(t)=\sqrt{(h|\hat{h}^{\text{R}}_{22})^{2}+(h|\hat{h}^{\text{I}}_{22})^{2}}$
(16)

In this symmetric case, the observed signal Eq. (14) can be expressed as [48]:

$h=\mathcal{A}\;\Re\\{h_{22}e^{i\phi_{C}}\\}=\mathcal{A}\,\left(\hat{h}^{\text{R}}_{22}\cos\phi_{C}-\hat{h}^{\text{I}}_{22}\sin\phi_{C}\right)$
(17)

where $\mathcal{A}$ is an amplitude factor and $\phi_{C}$ is a constant phase
shift, both depending on $\mathcal{F}_{+},\mathcal{F}_{\times},\iota$ and
$\varphi$.111 The factorization Eq. (17) is very convenient, since all the
extrinsic parameters are absorbed into an overall amplitude and a phase, which
the matched filtering is able to analytically maximise over. The same does not
happen for the general case and the maximisation over the inclination and
reference phase is performed by brute force (with a template bank). Hence, the
predicted timeseries $R_{\text{std}}(t)$ around a trigger at $t=0$ is given
by:

$\displaystyle R_{\text{std}}(t)$
$\displaystyle=(h|\hat{h}^{\text{R}}_{22})(t)+i(h|\hat{h}^{\text{I}}_{22})(t)$
$\displaystyle=\mathcal{A}\,(\cos\phi_{C}-i\sin\phi_{C})\\{(\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{R}}_{22})(t)+i(\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{I}}_{22})(t)\\}$
$\displaystyle=z(0)\\{(\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{R}}_{22})(t)+i(\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{I}}_{22})(t)\\}$
(18)

where, to simplify the expression above, we used the fact that
$z(0)=\mathcal{A}\,e^{-i\phi_{C}}$ and the following symmetries, arising from
Eq. (3):

$\displaystyle(\hat{h}^{\text{I}}_{22}|\hat{h}^{\text{I}}_{22})(t)=(\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{R}}_{22})(t)$
(19)
$\displaystyle(\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{I}}_{22})(t)=-(\hat{h}^{\text{I}}_{22}|\hat{h}^{\text{R}}_{22})(t)=[\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{R}}_{22}](t)$
(20)

By comparing the measured timeseries $z(t)$ and the expected timeseries
$R_{\text{std}}(t)$ with Eq. (12), one can compute $\xi_{\text{std}}^{2}$ for
a non-precessing search.

The standard SNR Eq. (16) can also be expressed in terms of
$\hat{h}^{\text{R}}_{22}$ only:

$\rho_{\text{std}}(t)=\sqrt{(h|\hat{h}^{\text{R}}_{22})^{2}+[h|\hat{h}^{\text{R}}_{22}]^{2}}$
(21)

We can then also rewrite $R_{\text{std}}(t)$, using the last equality of Eq.
(20):

$R_{\text{std}}(t)=z(0)\\{(\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{R}}_{22})(t)+i[\hat{h}^{\text{R}}_{22}|\hat{h}^{\text{R}}_{22}](t)$
(22)

For standard signals, the two expressions are equivalent. However, for
precessing and/or HM they may give very different results, due to the breaking
of the symmetry Eq. (19-20).

##  III A New Generalized $\xi^{2}$ Signal-consistency Test

In the generic case, where $\tilde{h}_{+}\not\propto i\tilde{h}_{\times}$, it
is necessary to use a different set of templates to filter the interferometric
data [68, 69, 51]:

$\displaystyle h_{R}=\hat{h}_{+}$ (23) $\displaystyle
h_{I}=\hat{h}_{\perp}=\frac{1}{\sqrt{1-\hat{h}_{+\times}^{2}}}(\hat{h}_{\times}-\hat{h}_{+\times}\hat{h}_{+})$

where

$\hat{h}_{+\times}=(\hat{h}_{+}|\hat{h}_{\times})(t=0)$ (24)

The real quantity $\hat{h}_{+\times}$ is a crucial measure of the precession
and/or HM content of a template and will be of primary interest in the rest of
the paper. The non-precessing non-HM limit can be recovered by
$\hat{h}_{+\times}=0$.

Note that if $h_{+}$ and $h_{\times}$ are normalized, the template
$\hat{h}_{\perp}$ is normalized by definition:
$\langle\hat{h}_{\perp}|\hat{h}_{\perp}\rangle=1$. Moreover, $\hat{h}_{+}$ and
$\hat{h}_{\perp}$ are orthogonal vectors, i.e.
$(\hat{h}_{+}|\hat{h}_{\perp})=0$. Indeed, the vectors for
$\hat{h}_{+},\hat{h}_{\perp}$ follows the Gram-Schmidt “orthogonalization”
prescription to create a set of orthonormal basis from the set of basis
vectors $\\{\hat{h}_{+},\hat{h}_{\times}\\}$.

We call symphony SNR222 The term “symphony” comes from the title of a paper
describing the statistics [51]. the timeseries Eq. (11) produced with the
templates above:

$\rho_{\text{sym}}(t)=\sqrt{(h|\hat{h}_{+})^{2}+(h|\hat{h}_{\perp})^{2}}$ (25)

With such templates, the measured timeseries $R_{\text{sym}}(t)$ in terms of
$\hat{h}_{+},\hat{h}_{\times}$ is given by:

$\displaystyle R_{\text{sym}}(t)$
$\displaystyle=(h|\hat{h}_{+})(t)+i(h|\hat{h}_{\perp})(t)$
$\displaystyle=\mathcal{F}_{+}(\hat{h}_{+}|\hat{h}_{+})(t)+\mathcal{F}_{\times}(\hat{h}_{\times}|\hat{h}_{+})(t)$
$\displaystyle+i\frac{1}{\sqrt{1-\hat{h}^{2}_{+\times}}}\Bigl{\\{}\mathcal{F}_{+}(\hat{h}_{+}|\hat{h}_{\times})(t)+\mathcal{F}_{\times}(\hat{h}_{\times}|\hat{h}_{\times})(t)\Bigr{.}$
$\displaystyle\Bigl{.}-\hat{h}_{+\times}\mathcal{F}_{+}(\hat{h}_{+}|\hat{h}_{+})(t)-\hat{h}_{+\times}\mathcal{F}_{\times}(\hat{h}_{\times}|\hat{h}_{+})(t)\Bigr{\\}}$
(26)

Here the values of the scaled antenna patterns
$\mathcal{F}_{+},\mathcal{F}_{\times}$ are not known and they have to be
estimated starting from the value $z(0)$ of the complex SNR timeseries at the
time of the trigger. Solving for ${z(0)=R_{\text{sym}}(t=0)}$, we find that:

$\displaystyle\mathcal{F}_{+}$
$\displaystyle=\Re\\{z(0)\\}-\Im\\{z(0)\\}\;\frac{\hat{h}_{+\times}}{\sqrt{1-\hat{h}^{2}_{+\times}}}$
(27) $\displaystyle\mathcal{F}_{\times}$
$\displaystyle=\Im\\{z(0)\\}\;\frac{1}{\sqrt{1-\hat{h}^{2}_{+\times}}}$ (28)

This expression can be plugged into Eq. (26) to compute the expected SNR
timeseries $R_{\text{sym}}(t)$, ready to be used for a novel general signal-
consistency test $\xi_{\text{sym}}^{2}$.

Of course, in the limit $\hat{h}_{+\times}\to 0$ the newly defined consistency
test must recover the original computation Eq. (18). Indeed, for
$\hat{h}_{+\times}=0$, we have

$\displaystyle\hat{h}_{+}$ $\displaystyle=\hat{h}^{\text{R}}_{22}$
$\displaystyle\hat{h}_{\times}$
$\displaystyle=\begin{cases}\hat{h}^{\text{I}}_{22}\;\;\;&\text{if}\;\;\cos\iota>0\\\
-\hat{h}^{\text{I}}_{22}\;\;\;&\text{if}\;\;\cos\iota<0\\\ \end{cases}$

and replacing $\hat{h}_{+},\hat{h}_{\times}$ with
$\hat{h}^{\text{R}}_{22},\hat{h}^{\text{I}}_{22}$ in Eq. (26) gives the
desired result. Note that the sign ambiguity in $\hat{h}_{\times}$ does not
impact the final result.

In Fig. 1, we show a visual example of how well the predicted SNR timeseries
matches the real one, both for the “standard” and the “symphony” case. We
inject a precessing BBH with $SNR=20$ into gaussian noise and filter the data
with the same signal; for both cases, we plot the measured SNR timeseries
$\rho_{\text{std}}(t)$ and $\rho_{\text{sym}}(t)$ and the expected timeseries
$R_{\text{std}}(t)$ and $R_{\text{sym}}(t)$.

We note that $R_{\text{std}}(t)$ and $\rho_{\text{std}}(t)$ show mild
agreement with each other and the recovered SNR (i.e. the magnitude of the
peak of the SNR timeseries) is $\sim 10\%$ lower than the injected value333
This would be the outcome of performing a precessing search using a standard
pipeline, optimized for non-precessing signals.. On the other hand,
$R_{\text{sym}}(t)$ is able to capture all the complicated features of
$\rho_{\text{sym}}(t)$, as shown by the lower $\xi^{2}$ value.

###  III.1 Approximating the new $\xi^{2}$ test

Although the timeseries $R_{\text{sym}}(t)$ is an accurate approximation of
the “symphony” SNR, it is not expressible as a multiplication between the
measured complex SNR $z(0)$ and a template-dependent complex timeseries and
for this reason, its adoption increases the computational cost of the
consistency test and requires substantial changes to the existing software
infrastructure.

To avoid this shortcoming, we introduce here a generalization $R_{\text{mix}}$
of Eq. (18) that is suitable for precessing and/or HM signals.
$R_{\text{mix}}$ is intended to approximate $\rho_{\text{sym}}(t)$, by keeping
the simple functional form of the standard test. This is achieved by
recognizing that the scalar $\hat{h}_{+\times}$ Eq. (24) is typically small.
As a consequence, we discard from Eq. (26) all the timeseries with a prefactor
$\propto\hat{h}_{+\times}$ and we obtain:

$\displaystyle R_{\text{sym}}(t)\simeq$ $\displaystyle\;\Re
z(0)\\{(\hat{h}_{+}|\hat{h}_{+})(t)+i(\hat{h}_{+}|\hat{h}_{\times})(t)\\}$
$\displaystyle+i\Im
z(0)\\{(\hat{h}_{\times}|\hat{h}_{\times})(t)+i(\hat{h}_{\times}|\hat{h}_{+})(t)\\}$
(29)

We can then achieve the desired expression by recognizing that
${(\hat{h}_{+}|\hat{h}_{+})(t)\simeq(\hat{h}_{\times}|\hat{h}_{\times})(t)}$
and
${(\hat{h}_{\times}|\hat{h}_{+})(t)\simeq-(\hat{h}_{+}|\hat{h}_{\times})(t)}$.
Thus, Eq. (III.1) becomes:

$\displaystyle R_{\text{mix}}(t)=$ $\displaystyle
z(0)\,\frac{1}{2}\Big{\\{}\left[(\hat{h}_{+}|\hat{h}_{+})(t)+(\hat{h}_{\times}|\hat{h}_{\times})(t)\right]$
$\displaystyle+i\left[(\hat{h}_{+}|\hat{h}_{\times})(t)-(\hat{h}_{\times}|\hat{h}_{+})(t)\right]\Big{\\}}$
(30)

By construction, $R_{\text{mix}}\simeq R_{\text{sym}}$, up to terms of order
$\hat{h}_{+\times}$. Finally, we introduce the additional “mixed” signal-
consistency test $\xi_{\text{mix}}^{2}$, which is obtained by comparing
$R_{\text{mix}}$ with the “symphony” SNR Eq. (25). The goodness of the
approximation given by $R_{\text{mix}}$ will be assessed in the next section.

In closing, we note that $R_{\text{mix}}$ is equivalent, up to terms of order
$\hat{h}_{+\times}$, to

$z(0)\\{(\hat{h}_{+}|\hat{h}_{+})(t)+i(\hat{h}_{+}|\hat{h}_{\times})(t)\\}$
(31)

The latter expression is not symmetric for an exchange between $\hat{h}_{+}$,
$\hat{h}_{\times}$ and for this reason we prefer to use Eq. (III.1). However,
Eq. (31) has a cleaner physical interpretation, since it amounts to Eq. (18)
with the natural replacement $\hat{h}^{\text{R}}_{22}\to\hat{h}_{+}$ and
$\hat{h}^{\text{I}}_{22}\to\hat{h}_{\times}$.

Figure 2:  Values of $\xi^{2}$ Eq. (12) as a function of the absolute value of
$(\hat{h}_{+}|\hat{h}_{\times})$, which quantifies the precession and/or HM
content of a signal. Each value is computed on random precessing BBHs,
injected into gaussian noise at a constant SNR; different panels show
different SNR. $\xi^{2}$ is computed using three different prescriptions.
$\xi_{\text{std}}^{2}$ is obtained from $R_{\text{std}}(t)$ and
$z_{\text{std}}(t)$ (label “standard”). $\xi_{\text{sym}}^{2}$ is computed
using $R_{\text{sym}}(t)$ and $z_{\text{sym}}(t)$ (label “symphony”), while
$\xi_{\text{mix}}^{2}$ uses $R_{\text{mix}}(t)$ and $z_{\text{sym}}(t)$. For
this study we set ACL $=701$.

##  IV Assessing the validity of the new statistics

In this section, we quantify the benefits of introducing the new statistics
$\xi_{\text{sym}}^{2}$ over the standard statistics $\xi_{\text{std}}^{2}$ and
compare it with the approximated test $\xi_{\text{mix}}^{2}$. We compute
$\xi_{\text{std}}^{2}$, $\xi_{\text{mix}}^{2}$ and $\xi_{\text{sym}}^{2}$ for
$15000$ randomly sampled BBH signals, injected into gaussian noise at
different values of SNR. We uniformly sample the total mass
${M\in[10,50]$\text{\,}\mathrm{M}_{\odot}$}$ and mass ratio
${q={m_{1}}/{m_{2}}\in[1,15]}$, while reference phase and inclination, as well
as the sky location, are drawn from a uniform distribution on the sphere.s We
also sample the starting frequency $f_{\text{min}}$ in the range
$[5,20]$\text{\,}\mathrm{H}\mathrm{z}$$.

In our study, we explore two scenarios. In one case, we focus on precessing
systems with both spins sampled isotropically inside the unit sphere. In the
other case, we consider aligned-spin systems but include higher modes (HMs) in
the waveform. For the latter experiment, we sample both z-components of the
spins uniformly between $[-0.99,0.99]$ and we consider the HM with
$(\ell,|m|)=(2,2),(2,1),(3,3),(3,2),(4,4)$. We utilize the frequency domain
approximants IMRPhenomXP [67] and IMRPhenomXHM [70] for the two scenarios
respectively. We employ the PSD computed over the first three month of the
third observing run at the LIGO Livingston detector [66] and sample
$100\text{\,}\mathrm{s}$ of gaussian noise at a sample rate of
$4096\text{\,}\mathrm{H}\mathrm{z}$ for each signal under study. We compute
the $\xi^{2}$ using a window (ACL) of $701$ points centered around the
injection time.

Results pertaining to the precessing signals are presented in Fig. 2, while
those concerning HMs are shown in Fig. 3. In both figures for varying SNR, we
plot the $\xi^{2}$ values as a function of
$\hat{h}_{+\times}=(\hat{h}_{+}|\hat{h}_{\times})$.

First of all, we note that in the “zero noise” case $\xi_{\text{std}}^{2}$ and
$\xi_{\text{mix}}^{2}$ are both non-zero. This means that $R_{\text{std}}(t)$
and $R_{\text{mix}}(t)$ are not able is not able to predict exactly the
behaviour of the SNR timeseries. This is expected, since in the precessing/HM
regime, they are both approximations to the true SNR. On the other hand,
$\xi_{\text{sym}}^{2}$ is always zero (up to numerical noise), showing the
appropriateness of the newly introduced $\xi_{\text{sym}}^{2}$.

In the presence of noise, the three tests trivially agree when
$|\hat{h}_{+\times}|$ goes to zero. Indeed, as noted several times, for
${|\hat{h}_{+\times}|\to 0}$, it holds
${\rho_{\text{sym}}\to\rho_{\text{std}}}$, hence the three tests yield
identical results.

By looking at the injected precessing signals in Fig. 2, we see that
$\xi_{\text{sym}}^{2}$ is always superior to the standard test, with an
improvement as large as two orders of magnitude in the SNR $=100$ case. As
long as $\xi_{\text{mix}}^{2}$ is considered, we observe that
$\xi_{\text{mix}}^{2}$ and $\xi_{\text{sym}}^{2}$ show similar dispersions in
the low SNR case. Therefore, in the presence of a substantial amount of noise,
the accuracy improvement provided by $\xi_{\text{sym}}^{2}$ is negligible over
the approximation given by $\xi_{\text{mix}}^{2}$. The discrepancy between the
two $\xi^{2}$ tests increases for $SNR=100$; in that case, the noise level is
lower and it must be of the same order of magnitude of the terms neglected to
obtain $\xi_{\text{mix}}^{2}$.

The fact that $\xi_{\text{mix}}^{2}$ degrades its performance at high SNR
should not be of concern, as the signal-consistency test is less crucial for
the high SNR region. Indeed, due to the rarity of very loud signals, it is
feasible to perform targeted follow up and ad hoc studies, hence assessing the
significance of a trigger with other strategies. For this reason, we conclude
that $\xi_{\text{mix}}^{2}$ is likely to perform close to optimality for the
vast majority of the practical applications and we recommend its
implementation in any pipeline aiming to search for precessing signals.

Figure 3:  Values of $\xi^{2}$ Eq. (12) as a function of the absolute value of
$(\hat{h}_{+}|\hat{h}_{\times})$, which quantifies the precession and/or HM
content of a signal. Each value is computed on random aligned spin BBHs with
HM content, injected into gaussian noise at a constant SNR; different panels
show different SNR. $\xi^{2}$ is computed using three different prescriptions.
$\xi_{\text{std}}^{2}$ is obtained from $R_{\text{std}}(t)$ and
$z_{\text{std}}(t)$ (label “standard”). $\xi_{\text{sym}}^{2}$ is computed
using $R_{\text{sym}}(t)$ and $z_{\text{sym}}(t)$ (label “symphony”), while
$\xi_{\text{mix}}^{2}$ uses $R_{\text{mix}}(t)$ and $z_{\text{sym}}(t)$. For
this study we set ACL $=701$.

The picture outline above changes when aligned-spin HM signals are considered.
In this case, the performance of $\xi_{\text{std}}^{2}$,
$\xi_{\text{mix}}^{2}$ and $\xi_{\text{sym}}^{2}$ are very comparable in the
low SNR case. In the high SNR case, $\xi_{\text{sym}}^{2}$ retains a slightly
better performance but the use of $\xi_{\text{mix}}^{2}$ does not bring any
additional improvement over the standard test $\xi_{\text{std}}^{2}$.
Therefore, if only HMs are considered, our results indicate that
$\xi_{\text{std}}^{2}$ will deliver close to optimal results.

In Fig. 4, we study how the three $\xi^{2}$ tests considered depend on the
choice of the integration window ACL. For brevity, we only consider the
$\xi^{2}$ values for precessing injections with SNR $=100$ and we choose three
different ACL $=351,701,1401$. Our results show that the $\xi^{2}$ values are
broadly consistent in the three cases, indicating the robustness of the test
against different choices of ACL.

In closing, we note that in the $\xi_{\text{std}}^{2}$ computation, we could
have used Eq. (16) instead of Eq. (21) to filter the data and, similarly, Eq.
(18) rather than Eq. (22) to compute the autocorrelation. Even though in the
standard case, the two expressions agree, they do not agree when precessing
and/or HM template are considered. In that case, the natural replacement
$\hat{h}^{\text{R}}_{22}\to\hat{h}_{+}$ and
$\hat{h}^{\text{I}}_{22}\to\hat{h}_{\times}$ leads Eqs. (16) and (18) to also
take into account the cross polarization. In this scenario, probably a
standard pipeline implementing Eq. (16) delivers smaller values of
$\xi_{\text{std}}^{2}$. However, as the GstLAL pipeline implements matched
filtering using Eq. (21), we made the choice to describe the current
situation.

Figure 4:  Values of $\xi^{2}$ for different window lengths (ACL) used for the
integral in Eq. (12). For each ACL, we report the values of
$\xi_{\text{std}}^{2}$ (blue), $\xi_{\text{mix}}^{2}$ (orange) and
$\xi_{\text{sym}}^{2}$ (green), following the color code introduced in Figs. 2
and 3. The data refers to $15000$ precessing signals injected at
$\text{SNR}=100$.

##  V Final remarks and future prospects

We introduced a novel signal-consistency test tailored for the search of GWs
emitted by precessing and/or HM binary systems. The new $\xi_{\text{sym}}^{2}$
test builds upon the “symphony” search statistics [51], enabling us to
formulate an expression for predicting the SNR time series obtained through
filtering interferometer data with a template. We demonstrate its
applicability in extending the $\xi^{2}$ statistical test currently employed
in non-precessing searches. To alleviate the computational cost of carrying
out the signal consistency test, we also proposed an approximation
$\xi_{\text{mix}}^{2}$ to our new $\xi_{\text{sym}}^{2}$, which requires
minimal changes to existing code platforms and is obtained without extra
computation with respect to the standard test. We measure the performance of
the two newly introduced tests on a set of precessing signals injected into
gaussian noise.

Based on our investigations, we conclude that the newly introduced test
$\xi_{\text{sym}}^{2}$ is the optimal signal-consistency test for filtering
precessing and/or HM signals with the search statistics introduced in [51].
However, the benefit introduced by the new test is less pronounced in the
aligned-spin HM case. The approximation $\xi_{\text{mix}}^{2}$ to our newly
introduced signal consistency test show a substantial improvement over the
traditional test in the precessing case, while no enhancement is observed for
the aligned-spin HM case.

Future work should repeat the injection study in a realistic noise realization
(i.e. with non gaussian detector’s noise). This will provide useful
information about the robustness of the new test in a real search. Moreover,
it will quantify the benefit brought by the approximation
$\xi_{\text{mix}}^{2}$ to our new test.

The approximated test $\xi_{\text{mix}}^{2}$ can be straightforwardly
implemented within the GstLAL pipeline, or in any other matched filtering
pipeline, thus marking a decisive step towards the implementation of a
precessing and/or HM search for BBH systems.

###### Acknowledgements.

We thank Khun Sang Phukon for the useful comments. S.S. and S.C. are supported
by the research program of the Netherlands Organization for Scientific
Research (NWO). S.C. is supported by the National Science Foundation under
Grant No. PHY-2309332. The authors are grateful for computational resources
provided by the LIGO Laboratory and supported by the National Science
Foundation Grants No. PHY-0757058 and No. PHY-0823459. This material is based
upon work supported by NSF’s LIGO Laboratory which is a major facility fully
funded by the National Science Foundation.

## References

  * [1] J. Aasi et al., “Advanced LIGO,” Class. Quant. Grav., vol. 32, p. 074001, 2015.
  * [2] F. Acernese et al., “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav., vol. 32, no. 2, p. 024001, 2015.
  * [3] T. Akutsu et al., “Overview of KAGRA: Detector design and construction history,” Progress of Theoretical and Experimental Physics, vol. 2021, p. 05A101, 08 2020.
  * [4] B. P. Abbott et al., “GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs,” Phys. Rev. X, vol. 9, no. 3, p. 031040, 2019.
  * [5] R. Abbott et al., “GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run,” Phys. Rev. X, vol. 11, p. 021053, 2021.
  * [6] R. Abbott et al., “GWTC-2.1: Deep extended catalog of compact binary coalescences observed by LIGO and Virgo during the first half of the third observing run,” Phys. Rev. D, vol. 109, no. 2, p. 022001, 2024.
  * [7] R. Abbott et al., “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run,” Phys. Rev. X, vol. 13, no. 4, p. 041039, 2023.
  * [8] B. S. Sathyaprakash and S. V. Dhurandhar, “Choice of filters for the detection of gravitational waves from coalescing binaries,” Phys. Rev. D, vol. 44, pp. 3819–3834, 1991.
  * [9] S. V. Dhurandhar and B. S. Sathyaprakash, “Choice of filters for the detection of gravitational waves from coalescing binaries. 2. Detection in colored noise,” Phys. Rev. D, vol. 49, pp. 1707–1722, 1994.
  * [10] B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E. Creighton, “FINDCHIRP: An Algorithm for detection of gravitational waves from inspiraling compact binaries,” Phys. Rev. D, vol. 85, p. 122006, 2012\.
  * [11] K. Cannon et al., “Toward Early-Warning Detection of Gravitational Waves from Compact Binary Coalescence,” Astrophys. J., vol. 748, p. 136, 2012.
  * [12] S. Babak et al., “Searching for gravitational waves from binary coalescence,” Phys. Rev. D, vol. 87, no. 2, p. 024033, 2013.
  * [13] S. Klimenko, I. Yakushin, A. Mercer, and G. Mitselmakher, “Coherent method for detection of gravitational wave bursts,” Class. Quant. Grav., vol. 25, p. 114029, 2008.
  * [14] V. Necula, S. Klimenko, and G. Mitselmakher, “Transient analysis with fast Wilson-Daubechies time-frequency transform,” J. Phys. Conf. Ser., vol. 363, p. 012032, 2012.
  * [15] M. Drago et al., “Coherent WaveBurst, a pipeline for unmodeled gravitational-wave data analysis,” 6 2020.
  * [16] B. P. Abbott et al., “GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence,” Phys. Rev. Lett., vol. 116, no. 24, p. 241103, 2016.
  * [17] B. P. Abbott et al., “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett., vol. 119, no. 16, p. 161101, 2017.
  * [18] S. Privitera, S. R. P. Mohapatra, P. Ajith, K. Cannon, N. Fotopoulos, M. A. Frei, C. Hanna, A. J. Weinstein, and J. T. Whelan, “Improving the sensitivity of a search for coalescing binary black holes with nonprecessing spins in gravitational wave data,” Phys. Rev. D, vol. 89, no. 2, p. 024003, 2014.
  * [19] T. Adams, D. Buskulic, V. Germain, G. M. Guidi, F. Marion, M. Montani, B. Mours, F. Piergiovanni, and G. Wang, “Low-latency analysis pipeline for compact binary coalescences in the advanced gravitational wave detector era,” Class. Quant. Grav., vol. 33, no. 17, p. 175012, 2016.
  * [20] S. A. Usman et al., “The PyCBC search for gravitational waves from compact binary coalescence,” Class. Quant. Grav., vol. 33, no. 21, p. 215004, 2016.
  * [21] C. Capano, I. Harry, S. Privitera, and A. Buonanno, “Implementing a search for gravitational waves from binary black holes with nonprecessing spin,” Phys. Rev. D, vol. 93, no. 12, p. 124007, 2016.
  * [22] C. Messick et al., “Analysis Framework for the Prompt Discovery of Compact Binary Mergers in Gravitational-wave Data,” Phys. Rev. D, vol. 95, no. 4, p. 042001, 2017.
  * [23] A. H. Nitz, T. Dent, T. Dal Canton, S. Fairhurst, and D. A. Brown, “Detecting binary compact-object mergers with gravitational waves: Understanding and Improving the sensitivity of the PyCBC search,” Astrophys. J., vol. 849, no. 2, p. 118, 2017.
  * [24] S. Sachdev, S. Caudill, H. Fong, R. K. L. Lo, C. Messick, D. Mukherjee, R. Magee, L. Tsukada, K. Blackburn, P. Brady, P. Brockill, K. Cannon, S. J. Chamberlin, D. Chatterjee, J. D. E. Creighton, P. Godwin, A. Gupta, C. Hanna, S. Kapadia, R. N. Lang, T. G. F. Li, D. Meacher, A. Pace, S. Privitera, L. Sadeghian, L. Wade, M. Wade, A. Weinstein, and S. L. Xiao, “The gstlal search analysis methods for compact binary mergers in advanced ligo’s second and advanced virgo’s first observing runs,” 2019.
  * [25] C. Hanna et al., “Fast evaluation of multidetector consistency for real-time gravitational wave searches,” Phys. Rev. D, vol. 101, no. 2, p. 022003, 2020.
  * [26] F. Aubin et al., “The MBTA pipeline for detecting compact binary coalescences in the third LIGO–Virgo observing run,” Class. Quant. Grav., vol. 38, no. 9, p. 095004, 2021.
  * [27] G. S. Davies, T. Dent, M. Tápai, I. Harry, C. McIsaac, and A. H. Nitz, “Extending the PyCBC search for gravitational waves from compact binary mergers to a global network,” Phys. Rev. D, vol. 102, no. 2, p. 022004, 2020.
  * [28] Q. Chu et al., “SPIIR online coherent pipeline to search for gravitational waves from compact binary coalescences,” Phys. Rev. D, vol. 105, no. 2, p. 024023, 2022.
  * [29] B. Ewing et al., “Performance of the low-latency GstLAL inspiral search towards LIGO, Virgo, and KAGRA’s fourth observing run,” 5 2023.
  * [30] L. Blackburn et al., “The LSC Glitch Group: Monitoring Noise Transients during the fifth LIGO Science Run,” Class. Quant. Grav., vol. 25, p. 184004, 2008.
  * [31] B. P. Abbott et al., “Characterization of transient noise in Advanced LIGO relevant to gravitational wave signal GW150914,” Class. Quant. Grav., vol. 33, no. 13, p. 134001, 2016.
  * [32] M. Cabero et al., “Blip glitches in Advanced LIGO data,” Class. Quant. Grav., vol. 36, no. 15, p. 15, 2019.
  * [33] S. Soni et al., “Reducing scattered light in LIGO’s third observing run,” Class. Quant. Grav., vol. 38, no. 2, p. 025016, 2020.
  * [34] D. Davis et al., “LIGO detector characterization in the second and third observing runs,” Class. Quant. Grav., vol. 38, no. 13, p. 135014, 2021.
  * [35] B. Allen, “${\chi}^{2}$ time-frequency discriminator for gravitational wave detection,” Phys. Rev. D, vol. 71, p. 062001, 2005.
  * [36] P. Shawhan and E. Ochsner, “A New waveform consistency test for gravitational wave inspiral searches,” Class. Quant. Grav., vol. 21, pp. S1757–S1766, 2004.
  * [37] K. Cannon, C. Hanna, and J. Peoples, “Likelihood-Ratio Ranking Statistic for Compact Binary Coalescence Candidates with Rate Estimation,” 4 2015.
  * [38] A. H. Nitz, “Distinguishing short duration noise transients in LIGO data to improve the PyCBC search for gravitational waves from high mass binary black hole mergers,” Class. Quant. Grav., vol. 35, no. 3, p. 035016, 2018.
  * [39] S. Dhurandhar, A. Gupta, B. Gadre, and S. Bose, “A unified approach to $\chi^{2}$ discriminators for searches of gravitational waves from compact binary coalescences,” Phys. Rev. D, vol. 96, no. 10, p. 103018, 2017.
  * [40] V. Gayathri, P. Bacon, A. Pai, E. Chassande-Mottin, F. Salemi, and G. Vedovato, “Astrophysical signal consistency test adapted for gravitational-wave transient searches,” Phys. Rev. D, vol. 100, no. 12, p. 124022, 2019.
  * [41] P. Godwin et al., “Incorporation of Statistical Data Quality Information into the GstLAL Search Analysis,” 10 2020.
  * [42] C. McIsaac and I. Harry, “Using machine learning to autotune chi-squared tests for gravitational wave searches,” Phys. Rev. D, vol. 105, no. 10, p. 104056, 2022.
  * [43] L. Tsukada et al., “Improved ranking statistics of the GstLAL inspiral search for compact binary coalescences,” Phys. Rev. D, vol. 108, no. 4, p. 043004, 2023.
  * [44] T. Dal Canton, A. P. Lundgren, and A. B. Nielsen, “Impact of precession on aligned-spin searches for neutron-star–black-hole binaries,” Phys. Rev. D, vol. 91, no. 6, p. 062010, 2015.
  * [45] I. W. Harry, A. H. Nitz, D. A. Brown, A. P. Lundgren, E. Ochsner, and D. Keppel, “Investigating the effect of precession on searches for neutron-star–black-hole binaries with advanced ligo,” Phys. Rev. D, vol. 89, p. 024010, Jan 2014.
  * [46] S. Fairhurst, R. Green, M. Hannam, and C. Hoy, “When will we observe binary black holes precessing?,” Phys. Rev. D, vol. 102, p. 041302, Aug 2020.
  * [47] N. Indik, K. Haris, T. Dal Canton, H. Fehrmann, B. Krishnan, A. Lundgren, A. B. Nielsen, and A. Pai, “Stochastic template bank for gravitational wave searches for precessing neutron-star–black-hole coalescence events,” Phys. Rev. D, vol. 95, no. 6, p. 064056, 2017.
  * [48] I. Harry, S. Privitera, A. Bohé, and A. Buonanno, “Searching for Gravitational Waves from Compact Binaries with Precessing Spins,” Phys. Rev. D, vol. 94, no. 2, p. 024012, 2016.
  * [49] C. McIsaac, C. Hoy, and I. Harry, “Search technique to observe precessing compact binary mergers in the advanced detector era,” Phys. Rev. D, vol. 108, no. 12, p. 123016, 2023.
  * [50] J. Calderón Bustillo, S. Husa, A. M. Sintes, and M. Pürrer, “Impact of gravitational radiation higher order modes on single aligned-spin gravitational wave searches for binary black holes,” Phys. Rev. D, vol. 93, no. 8, p. 084019, 2016.
  * [51] I. Harry, J. Calderón Bustillo, and A. Nitz, “Searching for the full symphony of black hole binary mergers,” Phys. Rev. D, vol. 97, no. 2, p. 023004, 2018.
  * [52] K. Chandra, J. Calderón Bustillo, A. Pai, and I. Harry, “First gravitational-wave search for intermediate-mass black hole mergers with higher order harmonics,” 7 2022.
  * [53] C. Mills and S. Fairhurst, “Measuring gravitational-wave higher-order multipoles,” Phys. Rev. D, vol. 103, p. 024042, Jan. 2021.
  * [54] D. Wadekar, T. Venumadhav, A. K. Mehta, J. Roulet, S. Olsen, J. Mushkin, B. Zackay, and M. Zaldarriaga, “A new approach to template banks of gravitational waves with higher harmonics: reducing matched-filtering cost by over an order of magnitude,” 10 2023.
  * [55] K. Cannon, S. Caudill, C. Chan, B. Cousins, J. D. E. 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, “Gstlal: A software framework for gravitational wave discovery,” 2020.
  * [56] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford Master Series in Physics, Oxford University Press, 2007.
  * [57] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, “Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries,” Phys. Rev. D, vol. 49, pp. 6274–6297, 1994.
  * [58] L. E. Kidder, C. M. Will, and A. G. Wiseman, “Spin effects in the inspiral of coalescing compact binaries,” Phys. Rev. D, vol. 47, no. 10, pp. R4183–R4187, 1993.
  * [59] L. E. Kidder, “Coalescing binary systems of compact objects to postNewtonian 5/2 order. 5. Spin effects,” Phys. Rev. D, vol. 52, pp. 821–847, 1995\.
  * [60] A. Buonanno, Y.-b. Chen, and M. Vallisneri, “Detecting gravitational waves from precessing binaries of spinning compact objects: Adiabatic limit,” Phys. Rev. D, vol. 67, p. 104025, 2003.  [Erratum: Phys.Rev.D 74, 029904 (2006)].
  * [61] M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krishnan, and D. Merritt, “Spin Flips and Precession in Black-Hole-Binary Mergers,” Phys. Rev. D, vol. 75, p. 064030, 2007.
  * [62] L. Pekowsky, J. Healy, D. Shoemaker, and P. Laguna, “Impact of higher-order modes on the detection of binary black hole coalescences,” Phys. Rev. D, vol. 87, no. 8, p. 084008, 2013.
  * [63] J. Healy, P. Laguna, L. Pekowsky, and D. Shoemaker, “Template Mode Hierarchies for Binary Black Hole Mergers,” Phys. Rev. D, vol. 88, no. 2, p. 024034, 2013.
  * [64] L. S. Finn and D. F. Chernoff, “Observing binary inspiral in gravitational radiation: One interferometer,” Phys. Rev. D, vol. 47, pp. 2198–2219, 1993.
  * [65] P. Jaranowski, A. Krolak, and B. F. Schutz, “Data analysis of gravitational - wave signals from spinning neutron stars. 1. The Signal and its detection,” Phys. Rev. D, vol. 58, p. 063001, 1998.
  * [66] Abbott, R. and others, “Noise curves used for Simulations in the update of the Observing Scenarios Paper.” https://dcc.ligo.org/LIGO-T2000012/public, 2022.  [Online; accessed 29-January-2024].
  * [67] G. Pratten et al., “Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes,” Phys. Rev. D, vol. 103, no. 10, p. 104056, 2021.
  * [68] C. Capano, Y. Pan, and A. Buonanno, “Impact of higher harmonics in searching for gravitational waves from nonspinning binary black holes,” Phys. Rev. D, vol. 89, no. 10, p. 102003, 2014.
  * [69] P. Schmidt, F. Ohme, and M. Hannam, “Towards models of gravitational waveforms from generic binaries II: Modelling precession effects with a single effective precession parameter,” Phys. Rev. D, vol. 91, no. 2, p. 024043, 2015.
  * [70] C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, “Multimode frequency-domain model for the gravitational wave signal from nonprecessing black-hole binaries,” Phys. Rev. D, vol. 102, no. 6, p. 064002, 2020.