# Effect of double spin-precession and higher harmonics on spin-induced
quadrupole moment measurements

Divyajyoti<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Madras, Chennai 600036, India  Centre for Strings,
Gravitation and Cosmology, Department of Physics, Indian Institute of
Technology Madras, Chennai 600036, India  N. V. Krishnendu
<EMAIL_ADDRESS>International Centre for Theoretical Sciences
(ICTS), Survey No. 151, Shivakote, Hesaraghatta Hobli, Bengaluru, 560089,
India  Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), Callinstr. 38, 30167 Hannover, Germany  Leibniz Universitat
Hannover, D-30167 Hannover, Germany  Muhammed Saleem<EMAIL_ADDRESS>School
of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
Marta Colleoni<EMAIL_ADDRESS>Departament de Física, Universitat de
les Illes Balears, IAC3 – IEEC, Crta. Valldemossa km 7.5, E-07122 Palma, Spain
Aditya Vijaykumar<EMAIL_ADDRESS>Canadian Institute for Theoretical
Astrophysics, University of Toronto, 60 St George St, Toronto, ON M5S 3H8,
Canada  International Centre for Theoretical Sciences (ICTS), Survey No. 151,
Shivakote, Hesaraghatta Hobli, Bengaluru, 560089, India  K. G. Arun
<EMAIL_ADDRESS>Chennai Mathematical Institute, Siruseri 603103, India
Chandra Kant Mishra<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Madras, Chennai 600036, India  Centre for Strings, Gravitation and
Cosmology, Department of Physics, Indian Institute of Technology Madras,
Chennai 600036, India

###### Abstract

We investigate the prospect of performing a null test of binary black hole
(BBH) nature using spin-induced quadrupole moment (SIQM) measurements. This is
achieved by constraining a deviation parameter ($\delta{\kappa}$) related to
the parameter ($\kappa$) that quantifies the degree of deformation due to the
spin of individual binary components on leading (quadrupolar) spin-induced
moment. Throughout the paper, we refer to $\kappa$ as the SIQM parameter and
$\delta\kappa$ as the SIQM-deviation parameter. The test presented here
extends the earlier SIQM-based null tests for BBH nature by employing waveform
models that account for double spin-precession and higher modes. We find that
waveform with double spin-precession gives better constraints for
$\delta\kappa$, compared to waveform with single spin-precession. We also
revisit earlier constraints on the SIQM-deviation parameter for selected GW
events observed through the first three observing runs (O1–O3) of LIGO-Virgo
detectors. Additionally, the effects of higher-order modes on the test are
also explored for a variety of mass-ratio and spin combinations by injecting
simulated signals in zero-noise. Our analyses indicate that binaries with
mass-ratio greater than 3 and significant spin precession may require
waveforms that account for spin-precession and higher modes to perform the
parameter estimation reliably.

##  I Introduction

There have been over 90 statistically significant detections of binary
coalescence events during the first three observing runs (O1-O3) of current
interferometric gravitational wave (GW) detectors Abbott  _et al._ (2016a,
2017, 2019, 2020a, 2021a). These include the LIGO Aasi  _et al._ (2015a);
Buikema  _et al._ (2020) and Virgo Acernese  _et al._ (2015a, 2019) detectors,
as well as KAGRA Akutsu  _et al._ (2021); Aso  _et al._ (2013) which has
recently joined the network. Upon detection, an elaborate analysis is
undertaken to deduce essential characteristics of the binary system,
encompassing properties like masses, spins, orientation and location in the
sky. This necessitates the utilization of accurate waveform models in
conjunction with an efficient parameter inference algorithm to guarantee the
accurate estimation of binary parameters Abbott  _et al._ (2016b). Each of
these observations engender many significant follow-up analyses. One such
endeavour involves determining the true nature of the compact objects in the
binary system. For instance, the nature of the secondary object in the binary
coalescence event GW190814 Abbott  _et al._ (2020b) is still a topic of
discussion as the secondary’s mass is consistent with the lightest black holes
(BHs) and heaviest neutron stars (NSs), along with other, more exotic
composition stars Clesse and Garcia-Bellido (2022); Vattis  _et al._ (2020);
Huang  _et al._ (2020); Roupas (2021); Biswas  _et al._ (2021); Tews  _et al._
(2021).

Various methods exist for investigating the true nature of the compact object
in a binary system Laarakkers and Poisson (1999); Ryan (1997); Poisson (1998);
Pacilio  _et al._ (2020); Gürlebeck (2015); Le Tiec and Casals (2021); Mendes
and Yang (2017); Uchikata  _et al._ (2016); Johnson-Mcdaniel  _et al._ (2020);
Jiménez Forteza  _et al._ (2018); Abdelsalhin  _et al._ (2018); Datta and Bose
(2019); Hartle (1973); Chatziioannou  _et al._ (2016); Datta  _et al._ (2021).
An analysis based on the spin-induced multipole moments is one among them.
Spin-induced multipole moments arise due to the spins of individual compact
objects in the binary Laarakkers and Poisson (1999). From observations, one
can measure these spin-induced multipole moments and then use that information
to distinguish black hole binaries from binaries composed of other compact
objects. The method based on the leading order spin-induced multipole moment,
the spin-induced quadrupole moment (SIQM) measurement, has been explored in
detail Krishnendu  _et al._ (2017, 2019a, 2019b); Krishnendu and Yelikar
(2020); Saleem  _et al._ (2022) and applied to the observed GW events from the
first three observing runs of advanced LIGO-Virgo detectors Abbott  _et al._
(2021b, c). Moreover, the possibility of measuring spin-induced quadrupole
using future detectors, and simultaneous measurement of spin-induced
quadrupole and octupole moment parameters Krishnendu  _et al._ (2019a); Saini
and Krishnendu (2023) have also been studied. A template bank for binaries of
exotic compact object searches has recently been developed, accounting for the
spin-induced quadrupole moment and tidal effects Chia  _et al._ (2022).
References Krishnendu  _et al._ (2017); Abbott  _et al._ (2021b, c) report the
tests of binary black hole nature using spin-induced quadrupole moments using
a phenomenological waveform model, IMRPhenomPv2 Husa  _et al._ (2016); Khan
_et al._ (2016); Hannam  _et al._ (2014), containing only the dominant modes
($\ell=2,|m|=2$) in the co-precessing frame and an effective spin parameter.
Recently, Ref. LaHaye  _et al._ (2023) came up with a fully precessing
waveform implementation of SIQM test for low-mass binaries, focusing on
binaries in which at least one object is in the lower mass gap ($<3M_{\odot}$)
Lyu  _et al._ (2023).

###  I.1 Current work

The effects of spin-precession and higher-order modes have been extensively
studied in the context of testing general relativity (TGR) using binary BHs
(see, for instance, Refs. Calderón Bustillo  _et al._ (2017); Krishnendu and
Ohme (2022, 2021); Puecher  _et al._ (2022); Mehta  _et al._ (2023); Islam
(2021); Breschi  _et al._ (2019)). By injecting the most up-to-date
phenomenological waveform models with full spin-precession (IMRPhenomXP) and
higher modes (IMRPhenomXHM, IMRPhenomXPHM Pratten  _et al._ (2021, 2020);
García-Quirós  _et al._ (2020)) for binary black hole signals of varying
masses and spins, we investigate the effects of spin-precession and higher
modes on $\delta\kappa$ measurements. Specifically, our injections include
binaries with mass ratios ($q=m_{1}/m_{2}$, where $m_{1}$ and $m_{2}$ are the
detector-frame component masses and $m_{1}>m_{2}$) in the range $q\in[1,5]$
and in-and-out-of-plane spin effects for a fixed mass binary
($M=30M_{\odot}$). We observe the difference in the bounds of $\delta\kappa$
between IMRPhenomXP (doubly spin-precessing) and IMRPhenomPv2 (singly spin-
precessing waveform model) in three scenarios: varying mass ratio, varying
effective aligned spin parameter ($\chi_{\text{eff}}$), and varying effective
spin-precession parameter ($\chi_{\text{p}}$). We note that in all cases
IMRPhenomXP outperforms IMRPhenomPv2. Next, we employ the higher mode waveform
models IMRPhenomXHM and IMRPhenomXPHM to study the effect of HMs on
$\delta\kappa$. We observe that for higher mass ratios, the higher mode
waveform models perform better compared to the dominant mode model
IMRPhenomXP. Following the TGR analyses on the second and third GW transient
catalogs (GWTC-2 Abbott  _et al._ (2021b) and GWTC-3 Abbott  _et al._
(2021c)), we measure $\delta\kappa$ of the binary systems using IMRPhenomXP
and IMRPhenomXPHM. We find that together with the effect of spin-precession,
the inclusion of higher modes plays a critical role when analysing binaries
with mass-asymmetries similar to that in the event GW190412 Abbott  _et al._
(2020c) (mass ratio $q\approx 3.7$). For GW190412, the bounds on
$\delta\kappa$ obtained with IMRPhenomXPHM are constrained enough to rule out
the boson star binaries, subject to the assumptions in the current work. We
also report the revised bounds on $\delta\kappa$ from selected GW events
observed through the first three observing runs of the LIGO-Virgo detectors.
This paper is organized as follows. In Sec II, we detail the waveform model
and parameter estimation method. Our results from simulated GW events are
shown in Sec. III, and real GW observations are reported in Sec. IV. We
conclude with Sec. V.

##  II Analysis Setup

In this section, we review the details of waveform models used and the basics
of Bayesian parameter estimation and hypothesis testing. Typically, the
evolution of an inspiralling compact binary can roughly be divided into three
stages: an early inspiral, late inspiral & merger, and the final ringdown.
During the early inspiral stage, the separation between the compact objects in
the binary is large, and hence, their evolution can be modeled as a
perturbation series in the velocity parameter. The post-Newtonian (PN) theory
provides an analytic expression for the inspiral phase incorporating various
physical effects such as the spin-orbit effects, self-spin effects, cubic and
higher order spin-effects, spin-precession effects, orbital eccentricity
effects, etc. (see, for instance, Blanchet (2002); Mishra  _et al._ (2016);
Henry  _et al._ (2022)). On the other hand, one needs to invoke numerical
relativity techniques to model the highly non-linear relativistic merger phase
(see Lehner (2001) for a review on numerical relativity modeling techniques).
Further, the ringdown part can be modeled perturbatively using the BH
perturbation theory techniques Sasaki and Tagoshi (2003); Pretorius (2007).

In the inspiral phase, the effect of spin-induced quadrupole moment starts to
appear at 2 PN, together with the other spin-spin terms. More precisely, the
leading order PN coefficient is of the schematic form Poisson and Will (1995)
$Q=-\kappa\chi^{2}m^{3}$, where the negative sign indicates the _oblate_
deformation due to the spinning motion. The proportionality constant,
$\kappa$, can take different values for different compact objects. For black
holes, $\kappa_{\text{BH}}$ is 1. Slowly spinning neutron stars can have
$\kappa$ values in the range $\rm{\kappa_{NS}}\sim 2-14$ Laarakkers and
Poisson (1999); Pappas and Apostolatos (2012a, b), whereas for more exotic
stars like boson stars, this range can be $\rm{\kappa_{BS}\sim}10-100$ Ryan
(1997) depending on internal composition. There also exist gravastar proposals
where the value $\rm{\kappa_{GS}}$ can match the BH value but also allows for
negative values and prolate deformations Uchikata and Yoshida (2016).
Measuring the SIQM parameter from GW observations can thus provide unique
information about the nature of the compact object.

For a binary system composed of two BHs with $\kappa_{i}$, following
Krishnendu  _et al._ (2019b), we define $\kappa_{i}=1+\delta\kappa_{i}$, where
$i=1,2$, and $\delta\kappa_{i}=0$ gives the BH limit. Since the simultaneous
measurements of both $\delta\kappa_{i}$ lead to uninformative results,111This
is with reference to the current detector sensitivities and dominant mode
waveform models. A detailed investigation can be carried out using future
detector sensitivities and/or higher mode waveform models, but it is beyond
the scope of this paper. we stick to the proposal of Krishnendu  _et al._
(2019b), where a symmetric combination of $\delta\kappa_{i}$ is measured
keeping the anti-symmetric combination ($\delta\kappa_{a}$) to zero. We call
this symmetric combination SIQM deviation parameter henceforth and use the
definition $\delta\kappa_{s}=(\delta\kappa_{1}+\delta\kappa_{2})/2$.

###  II.1 Waveform models used for the test

GW data analysis has routinely employed phenomenological waveform models with
varying properties. For instance, the SIQM analysis in the past was carried
out using the precessing dominant mode phenomenological waveform model
IMRPhenomPv2 Khan  _et al._ (2020). The two spin parameters, one for the out-
of-plane spin effects – the effective spin parameter ($\chi_{\mathrm{eff}}$)
Ajith  _et al._ (2011); Santamaria  _et al._ (2010), and another for the in-
plane spin effects – the spin-precession parameter ($\chi_{\mathrm{p}}$)
Schmidt  _et al._ (2012); Hannam  _et al._ (2014); Schmidt  _et al._ (2015),
are among the best-measured spin parameters in the parameter estimation of a
GW signal. It has been demonstrated that $\chi_{\mathrm{eff}}$ captures the
spin effects along the direction of the angular momentum axis Ajith  _et al._
(2011), and $\chi_{\mathrm{p}}$ measures the spin effects in the orbital plane
of the binary Schmidt  _et al._ (2015). The effective spin parameter for a
binary with dimensionless spin components,
$\chi_{i}=(\vec{S_{i}}\dotproduct\hat{L})/m_{i}^{2}$, can be defined as

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

Here, $\vec{S_{i}}$ is the individual spin angular momentum vector of the
compact object in the binary with mass $m_{i}$, and $\hat{L}$ represents the
unit vector along the angular momentum axis of the binary. In terms of the
perpendicular spin vectors,
$S_{i\perp}=|\hat{L}\times(\vec{S_{i}}\times\hat{L})|$, the effective spin-
precession parameter can be written as

${\chi_{\mathrm{p}}}=\frac{1}{A_{1}m_{1}^{2}}\max(A_{1}S_{1\perp},A_{2}S_{2\perp}),$
(2)

where $A_{1}=2+(3/2q)$ and $A_{2}=2+(3q/2)$ are mass parameters defined in
terms of the mass ratio $q=m_{1}/m_{2}>1$.

With recent developments in phenomenological modeling, it is now possible to
describe quasi-circular binaries involving generic spin components.
Specifically, the waveform IMRPhenomXP is the current state-of-the-art
phenomenological model where the two-spin effects are introduced Pratten  _et
al._ (2020). Further, the IMRPhenomXPHM model (which is also fully precessing)
includes higher order modes $(\ell,|m|)=(3,3),(4,4),(2,1),(3,2)$ in the co-
precessing frame, in addition to the dominant mode $(\ell,|m|)=(2,2)$ García-
Quirós  _et al._ (2020); Pratten  _et al._ (2021). We modify the inspiral
phase coefficients of these waveform models at 2 PN and 3 PN orders by
introducing explicit dependence on SIQM parameters, and study the measurement
probabilities for simulated GW signals from BBHs and the detected GW events.

###  II.2 Brief overview of Bayesian parameter estimation

Bayesian stochastic sampling algorithms are routinely employed to perform
parameter estimation of detected compact binary signals Veitch and Vecchio
(2010); Veitch  _et al._ (2015); Thrane and Talbot (2019). To review the
details of parameter estimation, we start with the definition of likelihood
function assuming Gaussian noise,

$\mathcal{L}(d|\theta)\propto\exp\quantity(-\frac{1}{2}\innerproduct{d-h\quantity(\theta)}{d-h\quantity(\theta)}),$
(3)

where the inner product $\innerproduct{a}{b}$ is defined as

$\innerproduct{a}{b}=4\int_{0}^{\infty}\frac{a^{*}(f)b(f)}{S_{n}(f)}.$ (4)

Here, $S_{n}$ is the noise power spectral density of the detector. The
posterior distribution on model parameters ($\theta$), given the data ($d$)
and likelihood (${\mathcal{L}}$), is

$p(\theta|d)=\frac{\mathcal{L}(d|\theta)\mathcal{\pi}(\theta)}{p(d)}\,.$ (5)

Here $p(d)$ is the Bayesian evidence or marginalized likelihood and is a
normalization constant obtained by integrating the likelihood over the entire
prior region for the set of parameters, $\theta$. A generic binary is
characterised by a set of intrinsic and extrinsic parameters. Adding to this
parameter set, in the current analysis, we include $\delta\kappa_{s}$ as a
free parameter to be constrained from the data. We can extract the
$\delta\kappa_{s}$ posterior by marginalizing over all other parameters
$\theta_{\rm BBH}$ from the multi-dimensional posterior samples as

$p(\delta\kappa_{s}|d)=\int p(\theta|d)\text{d}\theta_{\text{BBH}}.$ (6)

We use LALSimulation LIGO Scientific Collaboration (2018) for generating all
waveforms and the _nested_ sampling algorithm Skilling (2004) implemented
through dynesty Speagle (2020) sampler in bilby Ashton  _et al._ (2019) and
bilby_pipe Romero-Shaw  _et al._ (2020) for parameter estimation throughout
the analysis.

##  III Results from the simulations

Figure 1: The violin plots show the posterior distributions for the SIQM-
deviation parameter for various injection studies. Top: Four different mass
ratio cases are chosen, $q=(1,1.5,3,5)$; the spin magnitudes and related
angles are fixed to the values included in Table 1. Bottom Left: Three
different values for the effective spin parameter are chosen,
$\chi_{\text{eff}}$ = $0.3,0.5,0.7$ (see Table 2 for complete information);
mass ratio is taken to be $q=1.5$ and $\rm{\chi_{\mathrm{p}}=0.3}$. Bottom
right: Three different values for spin-precession parameter are chosen,
$\chi_{\mathrm{p}}$ = 0.3, 0.5, 0.7 (see Table 3 for information on other spin
parameters); mass ratio is taken to be $q=3$ and $\chi_{\text{eff}}=0.5$. The
total mass is fixed to $M=30~{}M_{\odot}$, and the network SNR is 40 for all
cases. The injections are performed using the fully spin-precessing dominant
mode waveform (IMRPhenomXP) and recovered with the same (orange) as well as
with single spin-precessing dominant mode waveform IMRPhenomPv2 (grey). The
horizontal black-dashed lines denote the injected value, and the coloured
lines inside the violins indicate the 90% credible intervals for the
respective posterior distributions. The legend follows the pattern “injected
waveform – recovery waveform”.

We demonstrate the importance of using a waveform model with double spin-
precession and higher modes for analysing binary black hole signals with
varying properties. Specifically, we look into different binaries of varying
mass asymmetries and in-plane & out-of-plane spin parameters, while fixing all
other parameters, to look into the effect of mass and spin variations in the
$\delta\kappa_{s}$ measurements. We consider four waveform models,
IMRPhenomXPHM (waveform with two-spin effects and higher modes), IMRPhenomXP
(waveform with two-spin effects and dominant mode),222 The solutions employed
in IMRPhenomX family are precession-averaged. IMRPhenomXHM (waveform with no
spin-precession effects but higher modes), and IMRPhenomPv2 (waveform model
with single-spin precession approximation and dominant mode).

For injections, we fix the total mass of the binary to be $30~{}M_{\odot}$ and
vary mass ratios and spins. All the binaries are placed in such a way that the
network signal-to-noise ratio (SNR) is 40 in a three-detector network
consisting of two advanced LIGO Abbott  _et al._ (2018); Harry and LIGO
Scientific Collaboration (2010); Aasi  _et al._ (2015b) and one advanced Virgo
Acernese  _et al._ (2015b); Acernese and et.al. (2006) detectors with advanced
sensitivity LIGO Scientific Collaboration (2022). All the injections
considered in this section are zero-noise injections and represent BBH mergers
(i.e., $\delta\kappa_{s}=0$). The GW signal from a binary system with spin-
precession is characterized by component masses ($m_{1},m_{2}$), dimensionless
spin magnitudes ($a_{1},a_{2}$), spin angles
($\Phi_{1},\Phi_{2},\Phi_{12},\Phi_{\mathrm{j}l}$), luminosity distance
($d_{L}$), angles measuring the location of the source in the sky and the
orientation of the source with respect to the line of sight
($\iota,\theta,\phi,\psi$), and the SIQM deviation parameter
($\delta\kappa_{s}$). We assume a uniform prior ranging [-500, 500] on
$\delta\kappa_{s}$. More details about the definition of these parameters and
corresponding prior ranges are given in Appendix A and Table 4.

###  III.1 Comparison of ${\delta\kappa_{s}}$ estimates: IMRPhenomXP and
IMRPhenomPv2

The effect of $\rm{\chi_{eff}}$ on the posteriors of $\delta\kappa_{s}$ is
well established and was explored in Ref. Krishnendu  _et al._ (2019b), albeit
using IMRPhenomPv2. In this section, we wish to compare the bounds on
$\delta\kappa_{s}$ obtained from IMRPhenomXP (doubly spin-precessing waveform)
and IMRPhenomPv2 (single spin-precession waveform). Furthermore, while the
study in Krishnendu  _et al._ (2019b) was performed for aligned-spin systems,
here, we choose systems with precessing spins and hence also explore the
effect of varying $\chi_{\mathrm{p}}$ on $\delta\kappa_{s}$. We consider three
cases:

  * •

Fixed spins with varying mass ratio to investigate the effect of mass ratio on
$\delta\kappa_{s}$.

  * •

Fixed masses and spin-precession parameter ($\chi_{\mathrm{p}}$), varying
effective spin parameter ($\rm{\chi_{\text{eff}}}$).

  * •

Different spin-precession values, keeping the masses and aligned-spin
components fixed.

For all cases, we inject binary black hole ($\delta\kappa_{s}=0$) signals with
precessing spins using IMRPhenomXP model and recover these with IMRPhenomXP
and IMRPhenomPv2. We present the results in the form of violin plots in Fig.
1.

####  III.1.1 Effect of mass ratio variation on $\delta\kappa_{s}$

To see the effect of mass ratio on the $\delta\kappa_{s}$ measurements, we
inject binaries with fixed spin magnitudes and angles, and vary the mass ratio
as $q=1,1.5,3,5$. The spin parameter values are listed in Table 1. The
injection signals are created using IMRPhenomXP and are analysed with both
IMRPhenomPv2 and IMRPhenomXP.

For all the cases, we observe that IMRPhenomXP outperforms IMRPhenomPv2,
especially for binaries with higher mass asymmetry. Moreover, we observe from
the top panel of Fig. 1 that an increase in mass ratio results in better
constraints on $\delta\kappa_{s}$. This is consistent with the findings of
Ref. Krishnendu  _et al._ (2019a) where the dependence of mass ratio on the
errors of SIQM parameter ($\Delta\kappa_{s}$) are discussed in detail [see
Fig. 2 and the discussion around Eqs. (4.2)–(4.4) there]. Moreover, the values
of $\chi_{\text{eff}}$ also increase gradually as we go from equal mass case
to unequal mass while keeping spin magnitudes and angles fixed. Both of these
effects result in improvement of $\delta\kappa_{s}$ bounds with increasing
mass ratio. We also show the corresponding corner plots on various parameters
in Fig. 5 in the Appendix.

Figure 2: Posterior distributions of the SIQM deviation parameter for various
mass ratios $q=1,1.5,3,5$. The total mass is fixed to $M=30~{}M_{\odot}$, and
network SNR is 40 for all cases. The spin magnitudes and angles are fixed (see
Table 1). We use fully spin-precessing, higher mode waveform model
IMRPhenomXPHM for injections and recover using the same (orange), aligned spin
higher mode waveform IMRPhenomXHM (grey), and fully spin-precessing dominant
mode waveform model IMRPhenomXP (blue) to study the effect of higher modes on
$\delta\kappa_{s}$ measurements. The vertical black-dashed lines denote the
injected value. The legend follows the pattern “injected waveform – recovery
waveform”.

####  III.1.2 Effect of $\chi_{\mathrm{eff}}$ variation on $\delta\kappa_{s}$

Here we choose a nearly equal-mass binary, with mass ratio $q=1.5$ and
$\rm{\chi_{\mathrm{p}}=0.3}$. Since we are not exploring the effect of HMs in
this injection set, a nearly equal-mass system serves the purpose well. We
vary the $\chi_{\mathrm{eff}}$ parameter as
$\rm{\chi_{\mathrm{eff}}=0.3,0.5,0.7}$. This is done by fixing the $\rm{x-}$
and $\rm{y-}$ components of the two spin vectors and varying the $\rm{z-}$
components $\chi_{\text{1z}}$ and $\chi_{\text{2z}}$ to obtain three values of
$\chi_{\text{eff}}$ as $0.3$, $0.5$, and $0.7$ (see Table 2).

As observed in Ref. Krishnendu  _et al._ (2019b), the estimates on
$\delta\kappa_{s}$ improve as we choose large positive $\chi_{\mathrm{eff}}$
values. Also, for all values of $\chi_{\mathrm{eff}}$, the bounds obtained
using IMRPhenomXP are better than IMRPhenomPv2. These improvements can be
explained by looking at correlations between $\chi_{\text{eff}}$ and
$\delta\kappa_{s}$ shown in Fig. 6 in the Appendix.

####  III.1.3 Effect of $\chi_{\mathrm{p}}$ variation on $\delta\kappa_{s}$

Here we choose a mass ratio of $q=3$ and a moderate value of
$\chi_{\text{eff}}=0.5$. A slightly larger mass ratio is chosen here compared
to Sec III.1.2 to avoid the uninformative inference on the analyses due to
unconstrained spin-precession effects for near-equal mass binaries. Keeping
the $\rm{z-}$ component of spin vectors the same, we vary the $\rm{x-}$ and
$\rm{y-}$ components to obtain three distinct values of $\chi_{\mathrm{p}}$ as
0.3, 0.5, and 0.7 (see Table 3).

As the values of $\chi_{\mathrm{p}}$ increase, the bounds on
$\delta\kappa_{s}$ with IMRPhenomXP become tighter, enhancing the differences
between IMRPhenomXP and IMRPhenomPv2 waveforms. The IMRPhenomPv2 bounds shift
away from the injected value as we move from low to high $\chi_{\mathrm{p}}$
values, excluding 0 from the 90% credible interval for
$\chi_{\mathrm{p}}=0.7$. Additionally, they become increasingly worse (the
posteriors become broader) compared to IMRPhenomXP as we go to higher values
of $\chi_{\mathrm{p}}$. We suspect that the doubly spin-precessing model
IMRPhenomXP is helping to break certain degeneracies between the SIQM
parameter and the spins leading to a more symmetric estimate of
$\delta\kappa_{s}$ for all the $\chi_{\mathrm{p}}$ values compared to the
IMRPhenomPv2 waveform model.

$q=\frac{m_{1}}{m_{2}}$ | $\rm{\chi_{1x}}$ | $\rm{\chi_{1y}}$ | $\rm{\chi_{1z}}$ | $\rm{\chi_{2x}}$ | $\rm{\chi_{2y}}$ | $\rm{\chi_{2z}}$ | $\chi_{\text{eff}}$ | $\chi_{\mathrm{p}}$
---|---|---|---|---|---|---|---|---
1 | 0.0992 | 0.1008 | 0.6 | 0.3343 | 0.3397 | 0.35 | 0.48 | 0.48
1.5 | 0.1013 | 0.0987 | 0.6 | 0.3414 | 0.3326 | 0.35 | 0.5 | 0.3
3 | 0.1015 | 0.0984 | 0.6 | 0.3422 | 0.3318 | 0.35 | 0.54 | 0.14
5 | 0.0997 | 0.1003 | 0.6 | 0.3359 | 0.3381 | 0.35 | 0.56 | 0.14
Table 1: Values of dimensionless spin components ($\rm{\chi_{1x}...\chi_{2z}}$), effective spin parameter ($\chi_{\text{eff}}$), and precessing-spin parameter ($\chi_{\mathrm{p}}$) for different mass ratios when the dimensionless spin magnitude and spin angles have been fixed as follows: $a_{1}=0.6164,~{}a_{2}=0.5913,~{}\theta_{\text{jn}}=0.4606,~{}\Phi_{\text{j}l}=3.7926,~{}\Phi_{1}=0.2315,~{}\Phi_{2}=0.9374,~{}\Phi_{12}=0.0$. This has been used for studying the effect of higher modes (Sec. III.2) and the effect of mass ratio variation (Sec. III.1.1). $\chi_{\mathrm{eff}}$ | $\rm{\chi_{1x}}$ | $\rm{\chi_{2x}}$ | $\rm{\chi_{1y}}$ | $\rm{\chi_{2y}}$ | $\rm{\chi_{1z}}$ | $\rm{\chi_{2z}}$
---|---|---|---|---|---|---
0.3 | 0.1013 | 0.3414 | 0.0987 | 0.3326 | 0.4 | 0.15
0.5 | 0.6 | 0.35
0.7 | 0.8 | 0.55
Table 2: Here we fix the mass ratio $q=1.5$, inclination angle $\iota=0.5162$ rad, and $\rm{x-}$ and $\rm{y-}$ components of dimensionless spin vectors while varying the $\rm{z-}$ components in order to obtain different values of $\chi_{\text{eff}}$. The value of $\chi_{\mathrm{p}}$ remains constant at 0.3. This has been used to study the effect of $\chi_{\text{eff}}$ on $\delta\kappa_{s}$ (Sec. III.1.2). $\chi_{\mathrm{p}}$ | $\rm{\chi_{1x}}$ | $\rm{\chi_{2x}}$ | $\rm{\chi_{1y}}$ | $\rm{\chi_{2y}}$ | $\rm{\chi_{1z}}$ | $\rm{\chi_{2z}}$
---|---|---|---|---|---|---
0.3 | 0.2792 | 0.1 | 0.11 | 0.1 | 0.56 | 0.32
0.5 | 0.4 | 0.2 | 0.3 | 0.2
0.7 | 0.5524 | 0.3 | 0.43 | 0.3
Table 3: Here we fix the mass ratio $q=3$, inclination angle $\iota=0.5149$
rad, and $\rm{z-}$ components of the dimensionless spin vectors while varying
the $\rm{x-}$ and $\rm{y-}$ components in order to obtain different values of
$\chi_{\mathrm{p}}$. The value of $\chi_{\text{eff}}$ remains constant at 0.5.
This has been used to study the effect of precession on $\delta\kappa_{s}$
(Sec. III.1.3).

###  III.2 Effect of higher modes and possible systematic biases

In Fig. 2, the $\delta\kappa_{s}$ posteriors are shown for simulated binary
signals with a total mass of $30~{}M_{\odot}$ and mass ratios $q=1,1.5,3,5$.
We fix the spin magnitudes ($a_{1}=0.6164,~{}a_{2}=0.5913$), spin angles
($\Phi_{\text{j}l}=3.7926,~{}\Phi_{1}=0.2315,~{}\Phi_{2}=0.9374,~{}\Phi_{12}=0.0$),
and inclination angle ($\theta_{\text{jn}}=0.4606$), taking a different mass
ratio in each case, leading to different values for the dimensionless spin
components
($\rm{\chi_{1x},\chi_{1y},\chi_{1z},\chi_{2x},\chi_{2y},\chi_{2z}}$) and hence
different values of $\chi_{\text{eff}}$ and $\chi_{\mathrm{p}}$ as listed in
Table 1. The histograms in each plot correspond to the difference in the
waveform model used in analysing these signals. For instance, we generate
simulations assuming IMRPhenomXPHM model and analyse them using the same
(orange), IMRPhenomXHM (grey), and IMRPhenomXP (blue) as shown in Fig. 2. The
aim is to demonstrate the importance of using a waveform model with higher
modes in measuring $\delta\kappa_{s}$ and to examine the possible biases that
could arise by not including them.

We see the significance of higher modes when we go to higher mass ratio
binaries. More precisely, the posteriors are tightly constrained to the true
value ($\delta\kappa_{s}=0$) as the binary becomes asymmetric, and the $90\%$
bound on $\delta\kappa_{s}$ improves from 1.3 to 0.18 (nearly seven times)
when moving from near-equal mass binary ($q=1.5$) to the most asymmetric
binary ($q=5$) when using the higher mode waveform IMRPhenomXHM. For
comparison, the improvement is only three times when we use the dominant mode
waveform IMRPhenomXP. For both $q=3$ and $q=5$ cases, IMRPhenomXHM constraints
are better than the other models. This is because the relative contribution
from the higher modes increases for higher mass ratios, and hence IMRPhenomXHM
performs better than IMRPhenomXP. While IMRPhenomXPHM is also a higher mode
waveform, the six spin parameters of the binary system increase the
dimensionality of the parameter space. This is not the case with IMRPhenomXHM
as it is an aligned-spin waveform model and has only two spin parameters.
Since the system is only mildly precessing ($\rm{\chi_{\mathrm{p}}\sim
0.14}$),333See Table 1 for values of other spin parameters. the additional
spin parameters do not contribute much toward improving the bounds on
$\delta\kappa_{s}$, and the four extra parameters in IMRPhenomXPHM add a
disadvantage compared to recovery with IMRPhenomXHM. We have already discussed
the effects of spin-precession in detail in Sec. III.1.3.

While all three waveforms are unable to recover the injected
$\delta\kappa_{s}=0$ value for $q=1$ case, we believe that this is because of
the prior railing effect, which arises when the injection is exactly at $q=1$.
Since nearly the entire posterior volume lies outside the injected value of
mass ratio (see Fig. 7), this causes a bias in the recovery of the chirp mass
parameter. Given that both $\chi_{\text{eff}}$ and $\delta\kappa_{s}$ are
correlated with chirp mass, it causes a bias in both of these parameters, and
the injected value of $\delta\kappa_{s}$ is not recovered. These correlations
can be seen in the corner plot (Fig. 7) where we see biases in the values of
$\mathcal{M}$ and $\chi_{\text{eff}}$, resulting in a biased
$\delta\kappa_{s}$ posterior.444In fact, by fixing q=1 for recovery waveforms,
we see that the IMRPhenomXPHM injections are recovered with both IMRPhenomXP
and IMRPhenomXPHM templates. Note also that the $\delta\kappa_{s}$ posteriors
are well constrained on the positive side compared to the poorly constrained
(with a long tail) negative side. This is likely due to our choice of aligned
spin systems (i.e, with $\chi_{\text{eff}}>0$) in our investigations, and
reverse trends are expected for anti-aligned systems ($\chi_{\text{eff}}<0$)
as was observed in Ref. Krishnendu  _et al._ (2019b) (see Sec. III B 1 there
for a detailed discussion).

Figure 3: Posterior distributions on SIQM-deviation parameter for observed GW
events. The curves labeled with IMRPhenomPv2 (blue) correspond to the previous
results Abbott  _et al._ (2021b, c) using IMRPhenomPv2 waveform model. These
are being compared with the new models IMRPhenomXP (grey) and IMRPhenomXPHM
(orange). The vertical dotted line indicates the BBH limit
$\delta\kappa_{s}=0$, and the numbers written inside the plots denote the 50%
quantiles with error bounds at 5% and 95% quantiles for different waveform
model recoveries in respective colours. Figure 4: Corner plot for GW190412
showing posteriors on the SIQM deviation parameter ($\delta\kappa_{s}$),
effective spin ($\chi_{\text{eff}}$), and mass ratio ($q=m_{1}/m_{2}$). We
show on the same plot results from single spin-precessing dominant mode
waveform model IMRPhenomPv2 (blue), fully spin-precessing dominant mode
waveform model IMRPhenomXP (grey), and fully spin-precessing HM waveform model
IMRPhenomXPHM (orange). The histograms shown on the diagonal of the plot are
1D marginalized posteriors for the respective parameters with vertical dashed
lines denoting $90\%$ credible intervals. The contours in the 2D plots are
also drawn for 90% credible interval. The titles on the 1D marginalized
posteriors for respective parameters and recoveries indicate 50% quantiles
with error bounds at 5% and 95% quantiles.

##  IV Effect of spin-induced orbital precession and higher modes on real
events

We compare the performance of different waveform models in constraining
$\delta\kappa_{s}$ from the GW transient catalogs in Fig. 3. We choose the
events with best bounds on $\delta\kappa_{s}$ from GWTC-1, 2, and 3 to
demonstrate the effect of waveform models on constraining the
$\delta\kappa_{s}$ parameter. The IMRPhenomPv2 results Abbott  _et al._
(2021b, c) are shown in blue dot-dashed lines along with the IMRPhemomXPHM
(orange) and IMRPhemomXP (grey) results.

For events with nearly equal mass ($q\approx 1.4$ \- 1.7) and slow spins
($\chi_{\text{eff}}\approx 0.05$ \- 0.2), GW151226 and GW170608, even if we
use the more informed models IMRPhemomXPHM and IMRPhemomXP, the bounds do not
alter considerably compared to IMRPhenomPv2. On the other hand, the posteriors
show a considerable difference from the IMRPhenomPv2 counterpart for GW190412.
Note that GW190412 is the first asymmetric BBH event ($q\approx 3.75$) with an
indication for moderate spins Abbott  _et al._ (2020c) and higher modes.
Hence, we expect the most noticeable effect on the bounds of
$\delta\kappa_{s}$ from this event. The requirement of using waveform models
with higher modes is also evident from the IMRPhemomXPHM and IMRPhemomXP
comparison for GW190412 as shown in Fig. 3. This can also be seen in the
corner plot shown in Fig. 4, where the bounds on mass ratio and
$\chi_{\text{eff}}$ are considerably different for IMRPhenomXPHM compared to
the dominant mode waveform models. As discussed in the previous sections,
bounds on $\delta\kappa_{s}$ strongly depend on these parameters, and hence,
for an event with non-negligible higher mode content, we observe that
IMRPhenomXPHM performs much better. In fact, the bounds for GW190412 obtained
using IMRPhenomXPHM exclude boson star binaries as the source, subject to the
assumptions made in our study (such as neglect of the tidal corrections,
assumption that $\delta\kappa_{a}=0$, and that spin-induced effects are
accounted for only in the inspiral part of the waveform).

##  V Summary

Spin-induced multipole moment-based tests were routinely employed to determine
the nature of compact binary signal during the first three observing runs of
the advanced LIGO and advanced Virgo detectors Abbott  _et al._ (2021b, c). In
this study, we extend the applicability of the test to binaries with spin-
precession effects not considered in previous versions of the test and discuss
the possible improvements in the measurement of the SIQM deviation parameter
using a more informed waveform model containing two spin-precession effects
and higher modes. Starting with a simulation study, we demonstrate the
applicability of the SIQM test on binaries with large spin-precession and
moderate mass asymmetries. We find that there are considerable differences in
the bounds of $\delta\kappa_{s}$ obtained using IMRPhenomPv2 compared to
IMRPhenomXP for highly spinning systems. We also report on the improvements
and biases observed in the $\delta\kappa_{s}$ bounds with the choice of
different mass ratios and compare them between IMRPhenomXP and IMRPhenomPv2
waveform models. Further, by injecting higher mode spin-precessing signals, we
find that higher mode waveforms are essential when analysing GW signals with
high mass asymmetries. Finally, we re-analyse selected events from GWTC-1, 2,
and 3 with the most up-to-date waveform models including two-spin precession
effects and higher modes. Our findings show that IMRPhemomXPHM may be
preferred for analysing GW events such as GW190412 Abbott  _et al._ (2020c),
where there is evidence for mass asymmetry and non-negligible spin effects.
While the current paper studies waveform systematics on the SIQM tests for
various spin and mass ratio configurations by injecting BBH waveforms
consistent with GR, a detailed follow-up study with non-GR injections may be
carried out in a future work.

###### Acknowledgements.

We thank Michalis Agathos for his useful comments on the manuscript. DJ thanks
Sayantani Datta, Pankaj Saini, and Sajad A. Bhat for useful discussions. We
thank Anuradha Gupta, Archisman Ghosh, and Ish Gupta for the review of the
waveform code, and Pankaj Saini and Sijil Jose for help with the review
readiness. N.V.K. acknowledges the support from the Science and Engineering
Research Board (SERB), Government of India, through the National Post Doctoral
Fellowship Grant (Reg. No. PDF/2022/000379). M.S. acknowledges the support
from the National Science Foundation with Grant Nos. PHY-1806630, PHY-2010970,
and PHY-2110238. A.V. is supported by the Department of Atomic Energy,
Government of India, under Project No. RTI4001. A.V. also acknowledges the
support of the Natural Sciences and Engineering Research Council of Canada
(NSERC) (funding reference number 568580). K.G.A. acknowledges support from
the Department of Science and Technology and the Science and Engineering
Research Board (SERB) of India via Swarnajayanti Fellowship Grant No.
DST/SJF/PSA-01/2017-18 and Core Research Grant No. CRG/2021/004565. K.G.A.
also acknowledges support from the Infosys Foundation. C.K.M. acknowledges the
support of SERB’s Core Research Grant No. CRG/2022/007959. M.C. acknowledges
funding from the Spanish Agencia Estatal de Investigación, Grant No.
IJC2019-041385. Computations were performed on the CIT cluster provided by the
LIGO Laboratory. N.V.K. acknowledges Max Planck Computing and Data Facility
computing cluster Cobra for early stages of the work. We acknowledge National
Science Foundation Grants PHY-0757058 and 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. We used the following software
packages: LALSuite LIGO Scientific Collaboration (2018), bilby Ashton  _et
al._ (2019), bilby_pipe Romero-Shaw  _et al._ (2020), PyCBC Nitz  _et al._
(2020), NumPy Harris  _et al._ (2020), PESummary Hoy and Raymond (2021),
Matplotlib Hunter (2007), Seaborn Waskom (2021), jupyter Kluyver  _et al._
(2016), dynesty Speagle (2020), and corner Foreman-Mackey (2016). This
document has LIGO preprint No. LIGO-P2300337.

##  Appendix A Additional information for injection runs

The parameter space for precessing spin recoveries includes the following
parameters: inverse mass ratio ($q_{\text{inv}}=m_{2}/m_{1}$),555We use the
word inverse here to indicate that this is inverse of the mass ratio we have
used throughout the paper. Please note that bilby uses the term mass_ratio for
this. chirp mass ($\mathcal{M}$), luminosity distance ($d_{L}$), cos of
inclination angle ($\cos{\theta_{\text{jn}}}$),666Please note that this is the
inclination angle as observed in the detector frame, i.e. the angle between
the line-of-sight $\vec{\text{n}}$ and total orbital momentum vector
$\vec{\text{j}}$ of the binary. For precessing systems, value of
$\theta_{\text{jn}}$ is different from $\iota$ which is the angle between the
orbital angular momentum vector $\vec{l}$ and normal vector $\vec{\text{n}}.$
geocentric time ($t_{c}$), phase angle ($\phi_{c}$), dimensionless spin
magnitudes ($a_{1}$, $a_{2}$), spin angles ($\Phi_{1}$, $\Phi_{2}$,
$\Phi_{12}$, $\Phi_{\text{j}l}$), right ascension ($\alpha$), declination
($\delta$), polarization angle ($\psi$), and symmetric combination of SIQM
deviation parameter ($\delta\kappa_{s}$). The priors for chirp mass have been
modified according to the injection value such that the lower limit is $\sim
0.5$ times the injected value, whereas the upper limit is $\sim 3$ times the
injected value. Accordingly, we have put constraints on component masses. The
priors for the rest of the parameters are given in Table 4. For aligned spin
recovery, instead of the six spin parameters mentioned above, we use the
parameters $\chi_{1}$ and $\chi_{2}$ using the AlignedSpin prior mentioned in
bilby, which puts a uniform range of [0,0.99] on the spin magnitudes. For real
event analyses, we have used the same priors as were used for the respective
analyses of the events in GWTC-2 Abbott  _et al._ (2021b) and GWTC-3 Abbott
_et al._ (2021c) TGR papers. All the injection runs mentioned in this paper
are for a binary black hole system
($\delta\kappa_{s}=0=\rm{\delta\kappa_{a}}$) with total mass of 30
$M_{\odot}$. The distance has been scaled in each case such that we get a
fixed network SNR of 40. We have arbitrarily chosen the angles $\alpha$,
$\delta$, $\psi$, and $\phi_{c}$ as 0, and $t_{\text{gps}}=1126259462$ s.

Parameter | Prior | Range
---|---|---
$q_{\text{inv}}$ | Uniform | $0.125\text{ - }1$
$d_{L}$ |  | Uniform
---
Source Frame
$100\text{ - }1000$ Mpc
$\cos{\theta_{\text{jn}}}$ | Uniform | $-1\text{ - }1$
$t_{c}$ | Uniform |  $t_{\text{gps}}+(-2\text{ - }2)$ s
$\phi_{c}$ | Uniform | $0\text{ - }2\pi$
$a_{1}$, $a_{2}$ | Uniform | $0\text{ - }0.99$
$\Phi_{1}$, $\Phi_{2}$A | Uniform sine | $0\text{ - }\pi$
$\Phi_{12}$ | Uniform | $0\text{ - }2\pi$
$\Phi_{\text{j}l}$ | Uniform | $0\text{ - }2\pi$
$\alpha$ | Uniform | $0\text{ - }2\pi$
$\delta$ | Uniform cos | $-\pi/2\text{ - }\pi/2$
$\psi$ | Uniform | $0\text{ - }\pi$
$\delta\kappa_{s}$ | Uniform | $-500\text{ - }500$
Table 4: Priors for parameters used in precessing spin recoveries.
66footnotetext: Denoted by the spin angles tilt1 and tilt2 in bilby Ashton
_et al._ (2019).

Figure 5: Corner plot for fixed spin magnitudes and angles given in Table 1
when the mass ratio is varied as 1 (top-left), 1.5 (top-right), 3 (bottom-
left), and 5 (bottom-right). The plots show 1D and 2D posteriors for the SIQM
deviation parameter ($\delta\kappa_{s}$), chirp mass ($\mathcal{M}$),
effective spin ($\chi_{\text{eff}}$), and mass ratio ($q$). The injections are
performed using the fully spin-precessing dominant mode waveform (IMRPhenomXP)
and recovered with the same (orange) as well as with single spin-precessing
dominant mode waveform IMRPhenomPv2 (grey). The histograms shown on the
diagonal of the plots are 1D marginalized posteriors for the respective
parameters with vertical dashed lines denoting $90\%$ credible intervals and
black lines indicating the injected value of the parameters. The contours in
the 2D plots are also drawn for 90% credible interval. The titles on the 1D
marginalized posteriors for respective parameters and recoveries indicate 50%
quantiles with error bounds at 5% and 95% quantiles.

Figure 6: Same as Fig. 5 but for fixed $q=1.5$, $\rm{\chi_{\mathrm{p}}=0.3}$
when $\chi_{\text{eff}}$ is varied as 0.3 (top-left), 0.5 (top-right), and 0.7
(bottom). Figure 7: Corner plot for the equal-mass case ($q=1$) showing
posterior on the SIQM deviation parameter ($\delta\kappa_{s}$), chirp mass
($\mathcal{M}$), mass ratio ($q$), the effective spin parameter
($\chi_{\text{eff}}$), and the spin-precession parameter
($\chi_{\mathrm{p}}$). We have used a fully spin-precessing waveform model
including HMs (IMRPhenomXPHM) for both injection and recovery. The histograms
shown on the diagonal of the plot are 1D marginalized posteriors for the
respective parameters with vertical dashed lines denoting $90\%$ credible
intervals. The contours in the 2D plots are also drawn for 90% credible
interval. The black lines denote the injected value of the parameters, and the
titles on the 1D marginalized posteriors for respective parameters indicate
50% quantiles with error bounds at 5% and 95% quantiles.

## References

  * Abbott  _et al._ (2016a) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific and Virgo Collaborations), Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2020a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020a), arXiv:2001.01761 [astro-ph.HE] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific and Virgo Collaborations), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Aasi  _et al._ (2015a) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015a), arXiv:1411.4547 [gr-qc] .
  * Buikema  _et al._ (2020) A. Buikema _et al._ , Phys. Rev. D 102, 062003 (2020), arXiv:2008.01301 [astro-ph.IM] .
  * Acernese  _et al._ (2015a) F. Acernese _et al._ (VIRGO), Class. Quant. Grav. 32, 024001 (2015a), arXiv:1408.3978 [gr-qc] .
  * Acernese  _et al._ (2019) F. Acernese _et al._ (Virgo Collaboration), Phys. Rev. Lett. 123, 231108 (2019).
  * Akutsu  _et al._ (2021) T. Akutsu _et al._ (KAGRA), PTEP 2021, 05A101 (2021), arXiv:2005.05574 [physics.ins-det] .
  * Aso  _et al._ (2013) Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto (KAGRA), Phys. Rev. D 88, 043007 (2013), arXiv:1306.6747 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241102 (2016b), arXiv:1602.03840 [gr-qc] .
  * Abbott  _et al._ (2020b) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 896, L44 (2020b), arXiv:2006.12611 [astro-ph.HE] .
  * Clesse and Garcia-Bellido (2022) S. Clesse and J. Garcia-Bellido, Phys. Dark Univ. 38, 101111 (2022), arXiv:2007.06481 [astro-ph.CO] .
  * Vattis  _et al._ (2020) K. Vattis, I. S. Goldstein, and S. M. Koushiappas, Phys. Rev. D 102, 061301 (2020), arXiv:2006.15675 [astro-ph.HE] .
  * Huang  _et al._ (2020) K. Huang, J. Hu, Y. Zhang, and H. Shen, Astrophys. J. 904, 39 (2020), arXiv:2008.04491 [nucl-th] .
  * Roupas (2021) Z. Roupas, Astrophys. Space Sci. 366, 9 (2021), arXiv:2007.10679 [gr-qc] .
  * Biswas  _et al._ (2021) B. Biswas, R. Nandi, P. Char, S. Bose, and N. Stergioulas, Mon. Not. Roy. Astron. Soc. 505, 1600 (2021), arXiv:2010.02090 [astro-ph.HE] .
  * Tews  _et al._ (2021) I. Tews, P. T. H. Pang, T. Dietrich, M. W. Coughlin, S. Antier, M. Bulla, J. Heinzel, and L. Issa, Astrophys. J. Lett. 908, L1 (2021), arXiv:2007.06057 [astro-ph.HE] .
  * Laarakkers and Poisson (1999) W. G. Laarakkers and E. Poisson, Astrophys. J. 512, 282 (1999), arXiv:gr-qc/9709033 .
  * Ryan (1997) F. D. Ryan, Phys. Rev. D 55, 6081 (1997).
  * Poisson (1998) E. Poisson, Phys. Rev. D 57, 5287 (1998).
  * Pacilio  _et al._ (2020) C. Pacilio, M. Vaglio, A. Maselli, and P. Pani, Phys. Rev. D 102, 083002 (2020), arXiv:2007.05264 [gr-qc] .
  * Gürlebeck (2015) N. Gürlebeck, Phys. Rev. Lett. 114, 151102 (2015), arXiv:1503.03240 [gr-qc] .
  * Le Tiec and Casals (2021) A. Le Tiec and M. Casals, Phys. Rev. Lett. 126, 131102 (2021), arXiv:2007.00214 [gr-qc] .
  * Mendes and Yang (2017) R. F. P. Mendes and H. Yang, Classical Quantum Gravity 34, 185001 (2017), arXiv:1606.03035 [astro-ph.CO] .
  * Uchikata  _et al._ (2016) N. Uchikata, S. Yoshida, and P. Pani, Phys. Rev. D 94, 064015 (2016), arXiv:1607.03593 [gr-qc] .
  * Johnson-Mcdaniel  _et al._ (2020) N. K. Johnson-Mcdaniel, A. Mukherjee, R. Kashyap, P. Ajith, W. Del Pozzo, and S. Vitale, Phys. Rev. D 102, 123010 (2020), arXiv:1804.08026 [gr-qc] .
  * Jiménez Forteza  _et al._ (2018) X. Jiménez Forteza, T. Abdelsalhin, P. Pani, and L. Gualtieri, Phys. Rev. D 98, 124014 (2018), arXiv:1807.08016 [gr-qc] .
  * Abdelsalhin  _et al._ (2018) T. Abdelsalhin, L. Gualtieri, and P. Pani, Phys. Rev. D98, 104046 (2018), arXiv:1805.01487 [gr-qc] .
  * Datta and Bose (2019) S. Datta and S. Bose, Phys. Rev. D 99, 084001 (2019), arXiv:1902.01723 [gr-qc] .
  * Hartle (1973) J. B. Hartle, Phys. Rev. D 8, 1010 (1973).
  * Chatziioannou  _et al._ (2016) K. Chatziioannou, E. Poisson, and N. Yunes, Phys. Rev. D 94, 084043 (2016), arXiv:1608.02899 [gr-qc] .
  * Datta  _et al._ (2021) S. Datta, K. S. Phukon, and S. Bose, Phys. Rev. D 104, 084006 (2021), arXiv:2004.05974 [gr-qc] .
  * Krishnendu  _et al._ (2017) N. V. Krishnendu, K. G. Arun, and C. K. Mishra, Phys. Rev. Lett. 119, 091101 (2017), arXiv:1701.06318 [gr-qc] .
  * Krishnendu  _et al._ (2019a) N. V. Krishnendu, C. K. Mishra, and K. G. Arun, Phys. Rev. D 99, 064008 (2019a), arXiv:1811.00317 [gr-qc] .
  * Krishnendu  _et al._ (2019b) N. V. Krishnendu, M. Saleem, A. Samajdar, K. G. Arun, W. Del Pozzo, and C. K. Mishra, Phys. Rev. D 100, 104019 (2019b), arXiv:1908.02247 [gr-qc] .
  * Krishnendu and Yelikar (2020) N. V. Krishnendu and A. B. Yelikar, Classical Quantum Gravity 37, 205019 (2020), arXiv:1904.12712 [gr-qc] .
  * Saleem  _et al._ (2022) M. Saleem, N. V. Krishnendu, A. Ghosh, A. Gupta, W. Del Pozzo, A. Ghosh, and K. G. Arun, Phys. Rev. D 105, 104066 (2022), arXiv:2111.04135 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific and Virgo Collaborations), Phys. Rev. D 103, 122002 (2021b), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, VIRGO, and KAGRA Collaborations), (2021c), arXiv:2112.06861 [gr-qc] .
  * Saini and Krishnendu (2023) P. Saini and N. V. Krishnendu, (2023), arXiv:2308.01309 [gr-qc] .
  * Chia  _et al._ (2022) H. S. Chia, T. D. P. Edwards, R. N. George, A. Zimmerman, A. Coogan, K. Freese, C. Messick, and C. N. Setzer, (2022), arXiv:2211.00039 [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] .
  * Khan  _et al._ (2016) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Hannam  _et al._ (2014) M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * LaHaye  _et al._ (2023) M. LaHaye, H. Yang, B. Bonga, and Z. Lyu, Phys. Rev. D 108, 043018 (2023), arXiv:2212.04657 [gr-qc] .
  * Lyu  _et al._ (2023) Z. Lyu, M. LaHaye, H. Yang, and B. Bonga, (2023), arXiv:2308.09032 [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] .
  * Krishnendu and Ohme (2022) N. V. Krishnendu and F. Ohme, Phys. Rev. D 105, 064012 (2022), arXiv:2110.00766 [gr-qc] .
  * Krishnendu and Ohme (2021) N. V. Krishnendu and F. Ohme, Universe 7, 497 (2021), arXiv:2201.05418 [gr-qc] .
  * Puecher  _et al._ (2022) A. Puecher, C. Kalaghatgi, S. Roy, Y. Setyawati, I. Gupta, B. S. Sathyaprakash, and C. Van Den Broeck, Phys. Rev. D 106, 082003 (2022), arXiv:2205.09062 [gr-qc] .
  * Mehta  _et al._ (2023) A. K. Mehta, A. Buonanno, R. Cotesta, A. Ghosh, N. Sennett, and J. Steinhoff, Phys. Rev. D 107, 044020 (2023), arXiv:2203.13937 [gr-qc] .
  * Islam (2021) T. Islam, (2021), arXiv:2111.00111 [gr-qc] .
  * Breschi  _et al._ (2019) M. Breschi, R. O’Shaughnessy, J. Lange, and O. Birnholtz, Class. Quant. Grav. 36, 245019 (2019), arXiv:1903.05982 [gr-qc] .
  * Pratten  _et al._ (2021) G. Pratten _et al._ , Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Pratten  _et al._ (2020) G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni, A. Ramos-Buades, H. Estelles, and R. Jaume, Phys. Rev. D 102, 064001 (2020), arXiv:2001.11412 [gr-qc] .
  * García-Quirós  _et al._ (2020) C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, Phys. Rev. D 102, 064002 (2020), arXiv:2001.10914 [gr-qc] .
  * Abbott  _et al._ (2020c) R. Abbott _et al._ (LIGO Scientific and Virgo Collaborations), Phys. Rev. D 102, 043015 (2020c), arXiv:2004.08342 [astro-ph.HE] .
  * Blanchet (2002) L. Blanchet, Living Rev. Rel. 5, 3 (2002), arXiv:gr-qc/0202016 .
  * Mishra  _et al._ (2016) C. K. Mishra, A. Kela, K. G. Arun, and G. Faye, Phys. Rev. D 93, 084054 (2016), arXiv:1601.05588 [gr-qc] .
  * Henry  _et al._ (2022) Q. Henry, S. Marsat, and M. Khalil, (2022), arXiv:2209.00374 [gr-qc] .
  * Lehner (2001) L. Lehner, Class. Quant. Grav. 18, R25 (2001), arXiv:gr-qc/0106072 .
  * Sasaki and Tagoshi (2003) M. Sasaki and H. Tagoshi, Living Rev. Rel. 6, 6 (2003), arXiv:gr-qc/0306120 .
  * Pretorius (2007) F. Pretorius, (2007), arXiv:0710.1338 [gr-qc] .
  * Poisson and Will (1995) E. Poisson and C. M. Will, Phys. Rev. D 52, 848 (1995), arXiv:gr-qc/9502040 .
  * Pappas and Apostolatos (2012a) G. Pappas and T. A. Apostolatos, Phys. Rev. Lett. 108, 231104 (2012a), arXiv:1201.6067 [gr-qc] .
  * Pappas and Apostolatos (2012b) G. Pappas and T. A. Apostolatos, (2012b), arXiv:1211.6299 [gr-qc] .
  * Uchikata and Yoshida (2016) N. Uchikata and S. Yoshida, Class. Quant. Grav. 33, 025005 (2016), arXiv:1506.06485 [gr-qc] .
  * Khan  _et al._ (2020) S. Khan, F. Ohme, K. Chatziioannou, and M. Hannam, Phys. Rev. D 101, 024056 (2020), arXiv:1911.06050 [gr-qc] .
  * Ajith  _et al._ (2011) P. Ajith _et al._ , Phys. Rev. Lett. 106, 241101 (2011), arXiv:0909.2867 [gr-qc] .
  * Santamaria  _et al._ (2010) L. Santamaria _et al._ , Phys. Rev. D 82, 064016 (2010), arXiv:1005.3306 [gr-qc] .
  * Schmidt  _et al._ (2012) P. Schmidt, M. Hannam, and S. Husa, Phys. Rev. D 86, 104063 (2012), arXiv:1207.3088 [gr-qc] .
  * Schmidt  _et al._ (2015) P. Schmidt, F. Ohme, and M. Hannam, Phys. Rev. D 91, 024043 (2015), arXiv:1408.1810 [gr-qc] .
  * Veitch and Vecchio (2010) J. Veitch and A. Vecchio, Phys. Rev. D 81, 062003 (2010), arXiv:0911.3820 [astro-ph.CO] .
  * Veitch  _et al._ (2015) J. Veitch _et al._ , Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Thrane and Talbot (2019) E. Thrane and C. Talbot, Publications of the Astronomical Society of Australia 36, e010 (2019).
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL) (2018).
  * Skilling (2004) J. Skilling, in  _Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering_, American Institute of Physics Conference Series, Vol. 735, edited by R. Fischer, R. Preuss, and U. V. Toussaint (2004) pp. 395–405.
  * Speagle (2020) J. S. Speagle, Monthly Notices of the Royal Astronomical Society 493, 3132 (2020).
  * Ashton  _et al._ (2019) G. Ashton _et al._ , Astrophys. J. Suppl. 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Romero-Shaw  _et al._ (2020) I. M. Romero-Shaw _et al._ , Mon. Not. Roy. Astron. Soc. 499, 3295 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * 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] .
  * Harry and LIGO Scientific Collaboration (2010) G. M. Harry and LIGO Scientific Collaboration, Classical and Quantum Gravity 27, 084006 (2010).
  * Aasi  _et al._ (2015b) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015b), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015b) F. Acernese _et al._ (VIRGO Collaboration), Classical Quantum Gravity 32, 024001 (2015b), arXiv:1408.3978 [gr-qc] .
  * Acernese and et.al. (2006) F. Acernese and P. A. et.al., Classical and Quantum Gravity 23, S635 (2006).
  * LIGO Scientific Collaboration (2022) LIGO Scientific Collaboration, https://dcc.ligo.org/LIGO-P1200087-v42/public (2022).
  * Nitz  _et al._ (2020) A. Nitz, I. Harry, D. Brown, C. M. Biwer, J. Willis, T. D. Canton, C. Capano, L. Pekowsky, T. Dent, A. R. Williamson, G. S. Davies, S. De, M. Cabero, B. Machenschalk, P. Kumar, S. Reyes, D. Macleod, dfinstad, F. Pannarale, T. Massinger, S. Kumar, M. Tápai, L. Singer, S. Khan, S. Fairhurst, A. Nielsen, S. Singh, and shasvath, “gwastro/pycbc: Pycbc release v1.16.11,”  (2020).
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature (London) 585, 357 (2020), arXiv:2006.10256 [cs.MS] .
  * Hoy and Raymond (2021) C. Hoy and V. Raymond, SoftwareX 15, 100765 (2021), arXiv:2006.06639 [astro-ph.IM] .
  * Hunter (2007) J. D. Hunter, Computing in Science and Engineering 9, 90 (2007).
  * Waskom (2021) M. L. Waskom, Journal of Open Source Software 6, 3021 (2021).
  * Kluyver  _et al._ (2016) T. Kluyver, B. Ragan-Kelley, F. Pérez, B. Granger, M. Bussonnier, J. Frederic, K. Kelley, J. Hamrick, J. Grout, S. Corlay, P. Ivanov, D. Avila, S. Abdalla, C. Willing, and J. development team, in  _Positioning and Power in Academic Publishing: Players, Agents and Agendas_ , edited by F. Loizides and B. Scmidt (IOS Press, 2016) pp. 87–90.
  * Foreman-Mackey (2016) D. Foreman-Mackey, The Journal of Open Source Software 1, 24 (2016).