# First frequency-domain phenomenological model of the multipole asymmetry in
gravitational-wave signals from binary-black-hole coalescence

Shrobana Ghosh  School of Physics and Astronomy, Cardiff University, Cardiff,
CF24 3AA, United Kingdom  Max Planck Institute for Gravitational Physics
(Albert Einstein Institute), Callinstrasse 38, D-30167 Hannover, Germany
Leibniz Universitaat Hannover, 30167 Hannover, Germany  Panagiota Kolitsidou
School of Physics and Astronomy, Cardiff University, Cardiff, CF24 3AA, United
Kingdom  School of Physics and Astronomy and Institute for Gravitational Wave
Astronomy, University of Birmingham, Edgbaston, Birmingham, B15 2TT, United
Kingdom  Mark Hannam  School of Physics and Astronomy, Cardiff University,
Cardiff, CF24 3AA, United Kingdom

###### Abstract

Gravitational-wave signals from binaries that contain spinning black holes in
general include an asymmetry between the $+m$ and $-m$ multipoles that is not
included in most signal models used in LIGO-Virgo-KAGRA (LVK) analysis to
date. This asymmetry manifests itself in out-of-plane recoil of the final
black hole and its inclusion in signal models is necessary both to measure
this recoil, but also to accurately measure the full spin information of each
black hole. We present the first model of the anti-symmetric contribution to
the dominant co-precessing-frame signal multipole throughout inspiral, merger
and ringdown. We model the anti-symmetric contribution in the frequency
domain, and take advantage of the approximations that the anti-symmetric
amplitude can be modelled as a ratio of the (already modelled) symmetric
amplitude, and analytic relationships between the symmetric and anti-symmetric
phases during the inspiral and ringdown. The model is tuned to single-spin
numerical-relativity simulations up to mass-ratio 8 and spin magnitudes of
0.8, and has been implemented in a recent phenomenological model for use in
the fourth LVK observing run. However, the procedure described here can be
easily applied to other time- or frequency-domain models.

##  I Introduction

The LIGO-Virgo-KAGRA (LVK) collaboration has published $\sim$90 gravitational-
wave (GW) observations Abbott  _et al._ (2019a, 2021a, 2021b) since the first
detection in 2015 Abbott  _et al._ (2016a). The majority of these have been
from binary black holes (BBHs), from which we are beginning to infer the
astrophysical distribution of black-hole masses and spins Abbott  _et al._
(2016b, 2019b, 2021c, 2023) and references therein. So far population
inference has had to rely on limited spin information from each binary; to
measure the magnitude and orientation of both spins we typically require
louder signals than in most of those observed so far Pürrer  _et al._ (2016);
Khan  _et al._ (2020). We also require sufficiently accurate and physically
complete theoretical waveform models. One physical effect that is necessary to
measure the full spin information is an asymmetry in the signals’ multipolar
structure that is not included in the standard full inspiral-merger-ringdown
(IMR) models of PHENOM Husa  _et al._ (2016); Khan  _et al._ (2016); Hannam
_et al._ (2014); London  _et al._ (2018); Khan  _et al._ (2019, 2020); Pratten
_et al._ (2020); García-Quirós  _et al._ (2020); Pratten  _et al._ (2021);
Thompson  _et al._ (2020); Estellés  _et al._ (2021, 2022a) or SEOBNR
Taracchini  _et al._ (2012); Pan  _et al._ (2014); Taracchini  _et al._
(2014); Bohé  _et al._ (2017); Cotesta  _et al._ (2018); Ossokine  _et al._
(2020); Matas  _et al._ (2020) family used in current LVK analyses.

A non-eccentric BBH is characterised by the black-hole masses $m_{1}$ and
$m_{2}$, and each black hole’s angular momentum, $\mathbf{S}_{i}$, which are
usually represented in geometric units as the dimensionless vectors
$\boldsymbol{\chi}_{i}=\mathbf{S}_{i}/m_{i}^{2}$. The dominant effect of the
spins on the GW signal is due to the spin components aligned with the binary’s
orbital angular momentum, which affect the rate of inspiral, and can therefore
be measured through their effect on the signal’s phase. The remaining (in-
plane) spin components have little effect on the inspiral rate. They instead
induce orbital and spin precession, which lead to modulations in the signal’s
amplitude and phase Apostolatos  _et al._ (1994). In most cases this is a
weaker contribution to the signal and more difficult to measure, in turn
making it difficult to measure the full spin information of the binary. Spin
misalignment also leads to an asymmetry in the power emitted above and below
the orbital plane, and can lead to large out-of-plane recoils of the final
black hole Bruegmann  _et al._ (2008a). Most signals to date have been too
weak to observe precession and recoil (with the notable exception of several
analyses of the signal GW200129_065458 signal – which we refer to as GW200129
in the rest of the text) Abbott  _et al._ (2021b); Hannam  _et al._ (2022);
Varma  _et al._ (2022), but more signals with measurable spin misalignment are
expected as detector sensitivities improve.

Most current generic-binary models separately consist of a model of the signal
in a non-inertial frame that tracks the precession of the orbital plane (a
“co-precessing” frame), and a model of the time- or frequency-dependent
precession angles. If the signal is decomposed into spin-weighted spherical
harmonics, the dominant contributions in the co-precessing frame are the
($\ell=2,|m|=2$) multipoles. As discussed in more detail in Sec. II, current
Phenom and EOB models assume the symmetry $h^{\rm CP}_{22}=h^{{\rm
CP}*}_{2,-2}$ in the co-precessing frame. Precessing binaries also include an
anti-symmetric contribution. There have been indications for some time that
neglecting this contribution could lead to measurement biases Kalaghatgi and
Hannam (2021); Ramos-Buades  _et al._ (2020), and more recently explicit
examples of such biases have been found Kolitsidou  _et al._ (2023). One model
that _does_ include the anti-symmetric contribution is the NR surrogate model
NRSur7dq4 Varma  _et al._ (2019), and this likely plays an important role in
being able to accurately infer the primary spin in GW200129 Hannam  _et al._
(2022); Varma  _et al._ (2022), demonstrating the need to include the anti-
symmetric contribution in Phenom and EOB models.

In this paper we present a simple method to model the anti-symmetric
contribution to the $(\ell=2,|m|=2)$ co-precessing-frame multipoles, taking
into account the phenomenology of how the anti-symmetric contribution depends
on the in-plane spin direction and relates to the symmetric contribution. Note
that all of the examples shown in this paper are constructed using either
numerical relativity waveforms or post-Newtonian estimates and that the model
introduced here is versatile in that it can be integrated into any frequency-
domain approximant.

To motivate our focus on only the anti-symmetric contribution to the dominant
multipoles, Fig. 1 shows the frequency-domain amplitude of the co-precessing-
frame multipoles for a signal with mass-ratio $q=m_{1}/m_{2}=2$, spin on the
larger black hole of $\chi=0.8$, and spin misalignment with the orbital
angular momentum of $\theta_{\rm LS}=90^{\circ}$, i.e., the spin initially
lies entirely in the orbital plane. (This is case CF_38 in Ref. Hamilton  _et
al._ (2023).) We see that the anti-symmetric $(2,2)$ amplitude is of
comparable strength to the symmetric $(3,3)$; since the $(3,3)$ extends to
higher frequencies, it will contribute more power than the anti-symmetric
$(2,2)$ at high masses, and comparable power in low-mass systems. The next-
strongest anti-symmetric contribution is to the $(3,3)$, and we see that this
is significantly weaker than the symmetric $(4,4)$. This suggests that any
model that includes symmetric contributions up to $\ell\leq 4$ need only
include the dominant $(2,2)$ anti-symmetric contribution. If we wish to
accurately model the signal to the level of the symmetric $(5,5)$
contribution, then we must also include the anti-symmetric $(3,3)$. Current
Phenom models include symmetric multipoles up to $\ell=4$, and so we limit our
attention to only the dominant anti-symmetric contribution. (Note that the
anti-symmetric (3,3) multipole is also weaker than the symmetric (2,1) and
(3,2) in this configuration.)

We find that we can model the $(2,2)$ anti-symmetric contribution using
numerical relativity (NR) simulations that cover only the reduced parameter
space of the binary’s mass ratio, the larger black hole’s spin magnitude, and
its misalignment angle; to a first approximation we _do not_ need to sample
the in-plane spin direction, which can be treated analytically. A more
complete model that includes sub-dominant in-plane-spin-direction effects, and
two-spin effects, is left to future work.

Figure 1:  Frequency-domain amplitude of co-precessing-frame symmetric and
anti-symmetric contributions for a $(q=2,\chi=0.8,\theta_{\rm LS}=90^{\circ})$
configuration. The figure shows the symmetric contributions to $(2,2)$,
$(3,3)$, $(4,4)$ and $(5,5)$, and the anti-symmetric contributions to $(2,2)$
and $(3,3)$.

In Sec. II we explain the motivation behind our modelling approach. We
describe the preparation of the NR data that we used to calibrate our model in
Sec. III. In Sec. IV we construct a model of the ratio of the anti-symmetric
and symmetric amplitudes, and in Sec. V, we present our method to construct
the anti-symmetric phase from the symmetric phase and the precession angle,
$\alpha$. We discuss the accuracy of our prescription in Sec. VI.

##  II Asymmetry background

An aligned-spin binary is invariant under reflection across the orbital plane.
If we choose a coordinate system where the orbital plane is the $x$-$y$ plane
and perform a decomposition of the gravitational-wave signal into spin-
weighted spherical harmonics, then this symmetry arises in the signal
multipoles as

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

This relationship is useful when constructing a model of the multipoles of an
aligned-spin binary: we need only explicitly model the positive-$m$
multipoles, and the negative-$m$ multipoles follow from Eq. (1).

When the spins are not aligned with the orbital angular momentum, both the
spins and the orbital plane may precess Apostolatos  _et al._ (1994); Kidder
(1995). In systems with mis-aligned spins Eq. (1) no longer holds, even if
there is no orbital precession. The simplest illustration of this is the
“superkick” configuration Gonzalez  _et al._ (2007); Campanelli  _et al._
(2007); Bruegmann  _et al._ (2008a): here the black holes are of equal mass
and of equal spin, and the spins lie in the orbital plane but point in
opposite directions. The symmetry of this system implies that the orbital
plane does not precess, and although the spins will precess _in the orbital
plane_ , they both precess at the same rate and so remain oppositely directed
to each other, so that the vector sum of the two spins is zero at all times.
Although the direction of the orbital plane remains fixed, emission of linear
momentum perpendicular to the orbital plane causes the entire system to move
up and down. This linear momentum emission and “bobbing” Keppel  _et al._
(2009) manifests itself in the gravitational-wave signal as an asymmetry in
the positive- and negative-$m$ multipoles, i.e., a violation of Eq. (1)
Bruegmann  _et al._ (2008a).

The symmetry of Eq. (1) will remain broken regardless of any rotations
performed on the multipoles Boyle  _et al._ (2014). This point becomes
important when constructing signal models, where we regularly make use of a
“co-precessing frame”. In this frame the signal during the inspiral can be
approximated as that of a non-precessing binary Schmidt  _et al._ (2012) and
so many current waveform models are split into a model for aligned-spin
binaries and a model for the precession dynamics, and the precession dynamics
are then used to “twist up” a non-precessing-binary waveform to produce the
complete precessing-binary waveform Hannam  _et al._ (2014); Khan  _et al._
(2019, 2020); Pratten  _et al._ (2021); Estellés  _et al._ (2022b); Pan  _et
al._ (2014); Taracchini  _et al._ (2014); Ossokine  _et al._ (2020); Ramos-
Buades  _et al._ (2023). However, since the non-precessing-binary waveform
respects the symmetry in Eq. (1), the model cannot reproduce the asymmetry
that should be present in the true precessing-binary signal.

Several studies have considered the impact of neglecting these multipole
asymmetries. Ref. Ramos-Buades  _et al._ (2020) compares the multipoles from
precessing-binary waveforms with those from nominally equivalent non-
precessing binaries, to test a number of assumptions that go into the
construction of many commonly used waveform models, including neglecting the
multipole asymmetry. Ref. Kalaghatgi and Hannam (2021) compares NR waveforms
from configurations with different in-plane spin directions and magnitudes,
and argues that neglecting the multipole asymmetry may lead to parameter
biases even at moderate SNRs, and that including the multipole asymmetry will
be necessary to clearly measure in-plane spins and identify precessing
systems. Finally, Ref. Kolitsidou  _et al._ (2023) uses the surrogate model
NRSur7dq4 to identify the level of bias in binary measurement examples, and
confirms that neglecting the multipole asymmetry leads to biases in measuring
in-plane spins (but the masses and effective aligned spin $\chi_{\rm eff}$ are
less affected). They also confirm the importance of the multipole asymmetry in
precession measurements, showing that it had a significant impact on
measurement of the properties of the signal GW200129 Hannam  _et al._ (2022);
Varma  _et al._ (2022).

In the next section we will summarise the leading-order PN contribution to the
asymmetry, which provides some insight into the phenomenology of the multipole
asymmetry, and also motivate our modelling procedure. Although the multipole
asymmetry has been known for some time, and indeed is included in the standard
PN expressions that we use here, and is also discussed in detail in Ref. Boyle
_et al._ (2014), we are not aware of any prior treatment that discusses the
amplitude and phasing of the anti-symmetric (2,2) contribution in relation to
the symmetric contribution, or notes the simple dependence of the relative
phase between different in-plane spin directions, which is a key feature of
the asymmetry that we exploit in constructing our model.

###  II.1 Inspiral

To gain insight into the phenomenology of the multipole asymmetry during the
inspiral, we consider the leading-order post-Newtonian contributions to a
binary where only one black hole is spinning and the spin lies entirely in the
orbital plane. The binary consists of two black holes with masses $m_{1}$ and
$m_{2}$ and the dimensionless spin on the primary is $\chi=S_{1}/m_{1}^{2}$,
where $S_{1}$ is the magnitude of the black hole’s angular momentum. We use
the post-Newtonian expressions from Ref. Arun  _et al._ (2009), where in this
single-spin case the symmetric and anti-symmetric spin contributions are
$\chi_{s}=\chi_{a}=\chi/2$. The in-plane spin components incline the total
angular momentum $\mathbf{J}$ with respect to the normal to the orbital plane
(and direction of the Newtonian orbital angular momentum, $\mathbf{L}$) by an
angle $\iota$, and the azimuthal precession angle of $\mathbf{L}$ around
$\mathbf{J}$ is $\alpha$; this is also the azimuthal angle of the total in-
plane spin. As such, if we choose the instantaneous orbital plane to coincide
with the $x$-$y$ plane, then the entirely in-plane spin can be written as
$\boldsymbol{\chi}=\chi(\cos(\alpha),\sin(\alpha),0)$.

We start with the multipole $h_{22}$ as given in Eq. (B1) in Ref. Arun  _et
al._ (2009). Requiring symmetry due to exchanging black holes (see the
discussion prior to Eq. (4.15) in the same paper), leads to the relation
$h_{\ell m}(\Phi)=(-1)^{\ell+m}h_{\ell-m}^{*}(\Phi+\pi)$, where $\Phi$ is the
orbital phase. We can enter the instantaneous orbital plane by setting
$\iota=\alpha=0$, and from our choice of spin we can then substitute
$\chi_{ax}=\chi_{sx}=\chi\cos(\alpha)/2$ and
$\chi_{ay}=\chi_{sy}=\chi\sin(\alpha)/2$; see sec. VI.B of Ref. Hamilton  _et
al._ (2021) for details. We then have an approximation to the symmetric and
anti-symmetric contributions to the signal in the co-precessing frame,
$h_{22}^{\rm{CP},s}=(h_{22}^{\rm CP}+h_{2-2}^{*\rm CP})/2$ and
$h_{22}^{\rm{CP},a}=(h_{22}^{\rm CP}-h_{2-2}^{*\rm CP})/2$, up to $O(v^{4})$,

$\displaystyle h_{22}^{\rm{CP},s}$ $\displaystyle=$ $\displaystyle
A\left(1+\frac{(55\eta-107)v^{2}}{42}\right)e^{-2i\Phi},$ (2) $\displaystyle
h_{22}^{\rm{CP},a}$ $\displaystyle=$ $\displaystyle
A\frac{(1+\delta)\chi\,v^{2}}{4}e^{-i(\Phi+\alpha)},$ (3)

where the overall amplitude is $A=\sqrt{64\pi/5}M\eta v^{2}/D_{L}$, $M$ is the
total mass, $\eta=m_{1}m_{2}/M^{2}$ is the symmetric mass ratio,
$\delta=(m_{1}-m_{2})/M$, $v$ is the relative speed of the two black holes,
and $D_{L}$ is the luminosity distance to the source. Note that in Ref. Arun
_et al._ (2009) the symbol $\Phi$ denotes the orbital phase in the
instantaneous orbital plane, but here we use it to denote the total orbital
phase that enters in the waveform.

We may immediately note several important features from Eqs. (2) and (3). The
spin does not enter the amplitude of the symmetric contribution at this order.
The anti-symmetric contribution enters at $O(v^{2})$ lower than the symmetric
contribution. We see that we may also consider the anti-symmetric amplitude
$|h_{22}^{\rm{CP},a}|$ as a simple rescaling of the symmetric amplitude,
$|h_{22}^{\rm{CP},s}|$.

The in-plane spin direction, $\alpha$, does not enter into the amplitude of
the anti-symmetric contribution, but does modify the phase. The physical
interpretation is that the phase of the anti-symmetric contribution depends on
the direction of the in-plane spin relative to the separation vector of the
two black holes. This will vary with the orbital phase, $\Phi$, but also the
(slower) precession rotation of the spin, given by $\alpha$. This is also
consistent with the observation in studies of out-of-plane recoil, that the
recoil amplitude depends sinusoidally on the initial direction of the in-plane
spin Bruegmann  _et al._ (2008a).

Finally, we note a key observation for the model that we will produce: if we
modify the initial in-plane spin direction by $\Delta\alpha$, this will induce
a simple overall phase shift in the anti-symmetric contribution,
$h_{22}^{\rm{CP},a}$. This suggests that, given a set of single-spin
numerical-relativity (NR) waveforms that cover the parameter space of mass
ratio, aligned-spin magnitude and in-plane spin magnitude, we will have enough
information to build a model of the anti-symmetric contribution to single-spin
waveforms _without the need to also sample multiple initial in-plane spin
directions_. We have just such a set of waveforms to hand, as used to
construct the first NR-tuned full inpiral-merger-ringdown model of (the
symmetric contribution to) precessing-binary waveforms Hamilton  _et al._
(2021), and discussed in Ref. Hamilton  _et al._ (2023).

###  II.2 Merger and ringdown

Before proceeding to construct a model, we consider the phenomenology of the
anti-symmetric contribution through merger and ringdown, and inspect which
inspiral features hold for the entire waveform.

One of the main features of the anti-symmetric contribution that we see in the
leading-order inspiral single-spin expressions (2) and (3) is that an in-plane
rotation of the spin by an angle $\Delta\alpha$ results in a corresponding
shift in the anti-symmetric (2,2) phase by $\Delta\alpha$. This is evident
from Fig. 2; the two configurations considered here correspond to the
superkick configuration described earlier. $\vec{S}^{\perp}_{1}$ denotes the
initial in-plane spin vector of the primary and $\vec{r}_{12}$ is the initial
separation vector pointing from the primary to the secondary. It is clear that
the asymmetry phase for a configuration with
$\vec{S}^{\perp}_{1}\perp\vec{r}_{12}$ can be easily produced by applying a
phase shift of $\pi/2$ to the anti-symmetric waveform of a configuration with
$\vec{S}^{\perp}_{1}\parallel\vec{r}_{12}$. We note that the simple phase
relationship does not appear to hold as well through merger and ringdown, but
the deviation is small enough that this could be due to numerical error, and
requires more detailed study in future.

Figure 2: Anti-symmetric waveform constructed from NR waveforms for two
superkick configurations in the time-domain,
$\vec{S}^{\parallel}_{1}\perp\vec{r}_{12}$ in dashed black and
$\vec{S}^{\perp}_{1}\perp\vec{r}_{12}$ in solid grey, shows a constant phase
offset of $\pi/2$. Inset: Dashed blue line shows that the anti-symmetric
waveform for $\vec{S}^{\perp}_{1}\perp\vec{r}_{12}$ can be constructed by just
applying a $\pi/2$ phase shift to the
$\vec{S}^{\perp}_{1}\parallel\vec{r}_{12}$ waveform even in the strong-field
regime close to merger $(t_{\rm{merger}}=1784M)$.

A second key feature of the anti-symmetric contribution is that its frequency
is roughly half that of the symmetric contribution (plus a small correction
from the spin-precession rate $\dot{\alpha}$). Fig. 3 shows the time-domain GW
frequency of the symmetric and anti-symmetric contributions for a
configuration where only the larger black hole is spinning, with the spin
($\chi=0.7$) entirely in-plane. We see that during the inspiral the anti-
symmetric frequency is approximately half that of the symmetric, as we expect.

Figure 3: The symmetric and anti-symmetric frequencies obtained from NR data
of a $(q=2,\chi=0.7,\theta_{\rm LS}=90^{\circ})$ configuration. During
inspiral the anti-symmetric frequency, $\omega_{a}$, is about half of the
symmetric frequency, $\omega_{s}$ (or nearly equal to the orbital frequency),
and close to merger quickly catches up with the symmetric frequency. As
expected from perturbation theory, $\omega_{s}=\omega_{a}$ during ringdown.

During ringdown the two frequencies are equal. This is consistent with our
expectations from perturbation theory. In the ringdown, where perturbation
theory results are applicable, the $(\ell=2,m=\pm 2)$ multipoles will have the
same (constant) frequency, and will decay exponentially at the same rate. We
therefore expect both the symmetric and anti-symmetric combinations of
$h_{22}$ and $h_{2,-2}$ to have the same frequency, and for the ratio of the
symmetric and anti-symmetric amplitudes to be constant throughout the
ringdown.

The third property we took from the inspiral expressions (2) and (3) was that
the anti-symmetric amplitude can be considered as a rescaling of the symmetric
amplitude. Fig. 4 illustrates that this holds for the entire waveform. It
shows the frequency-domain ampiltude of the symmetric and anti-symmetric (2,2)
contributions for a configuration with $(q=1,\chi=0.4,\theta_{\rm
LS}=60^{\circ})$, case CF_7 in Ref. Hamilton  _et al._ (2023). We see in
particular that the turnover to ringdown occurs at the same frequency (the
ringdown frequency is at $fM\sim 0.09$ for this configuration). We also see
that, as discussed above, the symmetric and anti-symmetric contributions decay
at the same rate.

Figure 4: The amplitude of the $(l=2,m=2)$ symmetric and anti-symmetric
waveform components in the frequency-domain co-precessing frame for a
$(q=1,\chi_{1}=0.4,\theta_{LS}=60^{\circ})$ NR simulation.

We will now exploit these features to construct a model for the anti-symmetric
amplitude and phase.

###  II.3 Structure of the model

Based on the observations in the previous section we construct an approximate
model of the anti-symmetric contribution to the (2,2) multipole in the co-
precessing-frame as follows. We start with a model for the symmetric
contribution, which provides us with the symmetric amplitude, $A_{s}(f)$, and
symmetric phase, $\phi_{s}(f)$. In the examples we consider here the symmetric
contribution is calculated from NR waveforms. In a full model, we would start
with the symmetric contribution from an already existing model. An explicit
example is the multi-mode precessing-binary model described in Ref. Thompson
_et al._ (2023), but the anti-symmetric model described in this paper can be
applied to any existing frequency-domain precessing-binary model that
separately provides $A_{s}(f)$, and symmetric phase, $\phi_{s}(f)$, and the
precession angle $\alpha(f)$, as we describe below.

To construct the anti-symmetric amplitude, $A_{a}(f)$, we model the ratio
$\kappa(f)=A_{a}(f)/A_{s}(f)$. In the inspiral we use a post-Newtonian
estimate of the amplitude ratio, and find that we can model the amplitude
ratio accurately through to the ringdown by adding only one higher-order term,
which we fit to our NR data. In the ringdown we treat the amplitude ratio as a
constant, as motivated in the previous section. The amplitude model is
presented in Sec. IV

To construct the anti-symmetric phase, during the inspiral we combine the
symmetric phase, $\phi_{s}(f)$, and the precession angle $\alpha(f)$, as
prescribed by Eq. (3), i.e., $\phi_{a}(f)\sim\phi_{s}(f)/2+\alpha(f)$. In the
merger-ringdown the anti-symmetric phase will behave as
$\phi_{a}(f)\sim\phi_{s}(f)$. We apply an overall time and phase shift to
smoothly connect these two functional forms at some transition frequency. We
note that we find that it is possible to produce an accurate model of the
anti-symmetric phase $\phi_{a}(f)$ using the same transition function for any
binary, i.e., the parameters of the transition function do not need to be fit
across the binary parameter space. The phase model is presented in Sec. V

##  III Numerical relativity data

The model is tuned to 80 single-spin NR waveforms generated using the BAM code
Bruegmann  _et al._ (2008b); Husa  _et al._ (2008). They cover mass ratios
$q=\\{1,2,4,8\\}$, spin magnitudes on the larger black holes of
$\chi_{1}=\\{0.2,0.4,0.6,0.8\\}$, and angles of spin misalignment with the
orbital angular momentum of $\theta_{\rm
LS}=\\{30^{\circ},60^{\circ},90^{\circ},120^{\circ},150^{\circ}\\}$. In
labelling the configurations, the cases are ordered according to the mass
ratios, then the spin magnitudes, then the misalignment angles; for example,
case 57 corresponds to $(q=4,\chi_{1}=0.8,\theta_{\rm LS}=60^{\circ})$. This
is the same indexing as in Ref. Hamilton  _et al._ (2023), which provides full
details on the catalog of simulations. To motivate the model and test our
modelling assumptions we have also used families of simulations that consider
variations in the initial in-plane-spin direction, based on those in Ref.
Kalaghatgi and Hannam (2021). One notable addition to this set were
simulations of superkick configurations where the black holes were given an
additional out-of-plane momentum, to remove a secular drift in the centre-of-
mass.

Several processing steps are performed to prepare the NR data for modelling,
using the tools in Ref. London (2015); London  _et al._ (2020). We start with
the waveform data for the $\ell=2$ multipoles of the Weyl scalar,
$\psi_{4,2m}$, in an inertial frame. We apply a Hann window to remove “junk”
radiation from the beginning of the waveforms (a non-physical artefact of the
initial data), and to remove numerical noise in the post-ringdown waveform.
Furthermore, to ensure that the frequency-domain step-size is sufficiently
small, the time-domain data are padded with zeros to the right. More details
are given in Ref. Hamilton  _et al._ (2023).

We then transform to a frame that tracks $\hat{J}(t)$; as described in Ref.
Hamilton  _et al._ (2021), this retains the approximation that the direction
of $\hat{J}(t)$ is constant throughout the waveform. Modelling deviations from
this approximation are left to future work, and are discussed further in Ref.
Hamilton  _et al._ (2021), along with the procedure to track $\hat{J}(t)$ and
perform this transformation. At the level of approximation and accuracy of the
anti-symmetric model presented here, we do not expect this approximation to
have any appreciable effect on the final model. We then transform the
$\psi_{4,2\pm 2}$ multipoles to a co-precessing frame, the “quadrupole-
aligned” (QA) Schmidt  _et al._ (2011) or “optimal emission direction”
O’Shaughnessy  _et al._ (2011); Boyle  _et al._ (2011) frame. In this frame
the multipoles are significantly simplified, with the $\left(l=2,|m|=2\right)$
multipoles having the strongest amplitude and the precession-induced
modulations minimised.

The $(2,\pm 2)$ multipoles of time-domain co-precessing-frame $\Psi_{4}$ are
now transformed to the frequency domain,

$\tilde{\psi}^{CP}_{4,2\pm 2}(f)=\int\psi^{CP}_{4,2\pm 2}(t)e^{-2\pi
ift}\,dt.$ (4)

To compute the strain, we note that $\Psi_{4}(t)=-\ddot{h}(t)$, so in the
frequency domain we may write,

$\tilde{h}^{CP}_{2,\pm 2}(f)=-\frac{\tilde{\psi}^{CP}_{4,2\pm
2}(f)}{\omega^{2}},$ (5)

where $\omega=2\pi f$. The anti-symmetric and symmetric components of the
waveform in the QA frame are computed from

$\displaystyle\tilde{h}_{s}^{NR}(f)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\tilde{h}^{CP}_{2,2}+\tilde{h}^{*CP}_{2,-2}),$ (6)
$\displaystyle\tilde{h}_{a}^{NR}(f)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\tilde{h}^{CP}_{2,2}-\tilde{h}^{*CP}_{2,-2}).$ (7)

The symmetric and anti-symmetric strains are complex quantities that can be
written as

$\displaystyle\tilde{h}_{s}^{NR}(f)$ $\displaystyle=$ $\displaystyle
A_{s}^{NR}(f)e^{i\phi^{NR}_{s}(f)},$ (8) $\displaystyle\tilde{h}_{a}^{NR}(f)$
$\displaystyle=$ $\displaystyle A_{a}^{NR}(f)e^{i\phi^{NR}_{a}(f)}$ (9)

and we can easily calculate their amplitude, $A^{NR}$, and phase, $\phi^{NR}$,
as their absolute value and argument, respectively. We denote the ratio
between the anti-symmetric and symmetric amplitudes as $\kappa_{\rm
NR}(f)=A_{a}^{NR}(f)/A_{s}^{NR}(f)$.

##  IV Model of the amplitude ratiio

We model the amplitude of the anti-symmetric contribution $A_{a}(f)$ as a
ratio of the symmetric contribution $A_{s}(f)$, i.e.,
$A_{a}(f)=\kappa(f)A_{s}(f)$.

Our first step in constructing the ratio model is to compute the ratio in the
framework of PN theory as a PN expansion in terms of $v/c$, where $v$ is the
relative velocity of the two black holes and $c$ is the speed of light, and we
choose geometric units where $c=1$. We again restrict our analysis to single-
spin binaries.

To compute the PN ratio, we follow a procedure similar to the illustrative
calculation in Sec. II. We obtain from Ref. Arun  _et al._ (2009) the complex
1.5PN expressions of the ($\ell=2,|m|=2)$ multipoles, $h^{PN}_{2\pm 2}$, for
spinning, precessing binaries with generic inclination angle $\iota$ moving on
nearly circular orbits. The strains of the $\ell=|m|=2$ multipoles can then be
transformed to a co-precessing frame that follows the instantaneous orbital
plane. To achieve this, we choose a simple co-precessing frame where we set to
zero the inclination angle of the orbital angular momentum relative to the
total angular momentum, $\iota=0$, and also set the precession angle,
$\alpha=0$. We then use an approximate reduction to a single-spin system
Hamilton  _et al._ (2021),

$\displaystyle\chi_{s/a,x}$ $\displaystyle=$
$\displaystyle-\chi\sin(\theta_{LS}-\iota)\cos(\alpha)/2$ (10)
$\displaystyle\chi_{s/a,y}$ $\displaystyle=$
$\displaystyle-\chi\sin(\theta_{LS}-\iota)\sin(\alpha)/2$ (11)
$\displaystyle\chi_{s/a,z}$ $\displaystyle=$
$\displaystyle\chi\cos(\theta_{LS}-\iota)/2$ (12)

where $\chi_{s}=(\chi_{1}+\chi_{2})/2$, $\chi_{a}=(\chi_{1}-\chi_{2})/2$ are
the symmetric and anti-symmetric spins defined in Ref. Arun  _et al._ (2009)
and $\theta_{LS}$ is the inclination of the spin from the orbital angular
momentum vector.

We then find that the symmetric and anti-symmetric amplitudes are,

$\displaystyle A_{s}^{PN}(f)$ $\displaystyle=$
$\displaystyle-\frac{4M}{21D_{L}}\sqrt{\frac{\pi}{5}}v^{2}\eta\left(42+84\pi
v^{3}+v^{2}(-107+55\eta)\right.$ (13) $\displaystyle\left.\ \ \ \ \ \ \ \ \ \
\ \ \ \ \ -28v^{3}(1+\delta-\eta)\chi\cos\theta_{LS}\right),$ $\displaystyle
A_{a}^{PN}(f)$ $\displaystyle=$
$\displaystyle-\frac{2M}{D_{L}}\sqrt{\frac{\pi}{5}}v^{4}(1+\delta)\eta\chi\sin\theta_{LS}.$
(14)

The ratio of the two amplitudes is then given by,

$\kappa_{\rm PN}(f)=\frac{21v^{2}(1+\delta)\chi\sin\theta_{LS}}{2(42+84\pi
v^{3}+v^{2}(-107+55\eta)-28v^{3}(1+\delta-\eta)\chi\cos\theta_{LS})},$ (15)

where $\eta=m_{1}m_{2}/M^{2}$ is the symmetric mass ratio and
$\delta=(m_{1}-m_{2})/M>0$ is a fractional mass difference. The expression of
the PN ratio of the anti-symmetric amplitude over the symmetric amplitude
depends on the symmetric mass ratio, $\eta$, the spin magnitude, $\chi$, and
the angle $\theta_{LS}$ of the system. Consequently, Eq. (15) can be used to
compute the PN ratio of any configuration as a function of frequency since
$v=(\pi fM)^{1/3}$.

We cannot expect the PN estimate of the amplitude ratio, Eq. (15), to be
accurate through merger and ringdown. To capture this behaviour we
investigated the addition of unknown higher-order terms and fit their
coefficients to the NR data. The simplest approach is to add only one
additional term, for example,

$\kappa(f)=\kappa_{\rm PN}(f)\,(1+bv^{n}),$ (16)

where $b$ is fit to the NR data. To do this it was necessary to choose an
appropriate frequency range over which to perform the fit. The NR data are
noisy in the early inspiral and in the post-ringdown frequencies, so we first
identified a consistent definition of a frequency range that could be used for
all 80 NR simulations, based on the ringdown frequency, $f_{RD}$ of each NR
configuration. The frequency range that we used in the final fits was given by
$Mf_{min}=\left(Mf_{RD}-0.01\right)/5$ and $Mf_{max}=Mf_{RD}-0.002$.

We investigated fits to the amplitude ratio of the form Eq. (16) with
$n=3,5,7$. Fig. 5 shows the results for one configuration, and illustrates
that $n=5$ provides the most accurate fit. In fact, we find that $n=5$
consistently provides the most accurate fit across all 80 NR simulations.

Figure 5: Amplitude ratio for the $(q=1,\chi=0.4,\theta_{LS}=60^{o})$
configuration, for the NR data, the PN ratio $\kappa_{\rm PN}$, and fits to
higher-order corrections as in Eq. (16), with $n=3,5,7$. The vertical dashed
line indicates the ringdown frequency.

It is also clear from Fig. 5 that our fit is not accurate beyond the ringdown
frequency. Beyond this point the amplitude ratio appears to plateau. This is
consistent with our expectations from perturbation theory, as discussed in
Sec. II.2, which predicts that the amplitude ratio will be constant throughout
ringdown. To include this feature in our model, we fix $\kappa(f)$ to the
value $\kappa(f_{RD})$ at frequencies $f>f_{RD}$. To avoid a sharp transition
we use a moving average algorithm such that,

$\kappa_{n}(f)=\frac{1}{(2k+1)}\sum_{i=n-k+1}^{n+k+1}\kappa(f_{i}),\;\;n\in[k,N-k].$
(17)

We use a symmetric window of an equal number of points ($k=40$) on either side
of the frequency $f$ to calculate the moving average. Here $N$ is the length
of the frequency series and the algorithm updates $\kappa(f)$ for $40\leq
n\leq N-40$.

We fit the coefficient $b$ in Eq. (16) (with $n=5$) to each of the $80$
single-spin NR waveforms from the NR catalogue in Ref. Hamilton  _et al._
(2023). Fig. 6 shows the resulting fit for the amplitude ratio for the
$(q=1,\chi=0.4,\theta_{LS}=60^{\circ})$ configuration. We see that the final
fit, including the levelling off of the amplitude ratio through ringdown,
agrees well with the NR data. The lower panel of the figure shows the
resulting estimate for the anti-symmetric amplitude. We see that this agrees
well with the NR data up to the point where the NR amplitude becomes dominated
by noise.

The fit cannot reproduce the NR amplitude ratio in all cases. In many cases
the NR amplitude ratio oscillates during the inspiral. Fig. 7 shows an extreme
example of this feature, from the $(q=1,\chi=0.8,\theta_{LS}=30^{\circ})$
configuration. It is not clear what causes these oscillations. Oscillations in
the co-precessing-frame amplitude ratio during the inspiral can be due to the
choice of co-precessing frame, as we will discuss later. However, if that were
the cause of the oscillations in the NR data, we would expect there to be some
correlation with the degree of precession in the configuration. Instead we do
not find any relationship between the amplitude of the oscillations, which in
some cases are negligible (as in Fig. 6), and in others (like Fig. 7) lead to
significant variations in the amplitude ratio. We do find, though, that our
model passes through a reasonable estimate of the mean of the oscillations.
The largest impact is on the constant value of the amplitude that the model
settles on for the ringdown regime; in the example in Fig. 7 the model’s
estimate of the anti-symmetric amplitude during the ringdown is roughly 20%
below the NR value.

Finally, in some cases we also found that the amplitude ratio in the NR data
did not level off during the ringdown. Fig. 8 is an example of this. We have
not been able to determine the reason for this. As noted previously, we expect
that since the $(2,2)$ and $(2,-2)$ multipoles decay at the same rate, that
the ratio between their amplitudes would remain constant during the ringdown.
It is possible that this effect is obscured by numerical noise. Regardless of
the cause, and in the absence of any compelling explanation for alternative
behaviour, we have chosen to impose the behaviour that we expect from
perturbation theory.

Figure 6: The amplitude ratio (top) and amplitude (bottom) for the
$(q=1,\chi=0.4,\theta_{LS}=60^{\circ})$ configuration. The vertical line
indicates the ringdown frequency for the $(\ell=2,|m|=2)$ multipoles.

Figure 7:  The amplitude ratio (top) and amplitude (bottom) for the
$(q=1,\chi=0.8,\theta_{LS}=30^{\circ})$ configuration. The vertical line
indicates the ringdown frequency for the $(\ell=2,|m|=2)$ multipoles.

Figure 8:  The amplitude ratio (top) and amplitude (bottom) for the
$(q=2,\chi=0.4,\theta_{LS}=90^{\circ})$ configuration. The vertical line
indicates the ringdown frequency for the $(\ell=2,|m|=2)$ multipoles.

The coefficient $b$ is shown as a function of symmetric mass ratio and
misalignment angle in Fig. 9; the values for all four spin magnitudes are
shown together. This plot illustrates the general trend of $b$ across the
parameter space. We find that there is no general trend with respect to the
spin magnitude, and this is illustrated more clearly in Fig. 10, which shows
$b$ as a function of symmetric mass ratio, for each spin magnitude, for the
subset of cases with $\theta_{\rm LS}=90^{\circ}$. The PN amplitude ratio
already includes a linear dependence on the spin magnitude, and given the
uncertainty in our fits, we do not attempt to identify a higher-order spin
dependence. We motivate this point further in Sec. VI (Fig. 14). We then
include all simulations into a two-dimensional parameter-space fit of the
form,

$b(\eta,\theta_{\rm
LS})=b_{0}+b_{1}\eta+b_{2}\theta_{LS}+b_{3}\eta\theta_{LS}$ (18)

where we find $b_{0}=18.0387$, $b_{1}=15.4509$, $b_{2}=55.1140$ and
$b_{3}=-203.6290$. From Eq. (18), we notice that $b$ does not go to zero when
$\theta_{LS}$ is $0^{\circ}$ or $180^{\circ}$. However, the presence of the
$\sin\theta_{LS}$ term at the numerator of the ansatz of the ratio model in
Eq. (15) ensures that the multipole asymmetry goes to zero at these points.
The amplitude as predicted by this fit is shown on each of our figures and
labelled as “final model”.

Two caveats to this approach are worth noting. One is the choice of co-
precessing frame. Previous work has shown that the symmetric (2,2)
contribution takes a simple form in the QA co-precessing frame; indeed, the
amplitude evolution can be approximated by the (monotonic) amplitude of an
equivalent aligned-spin configuration. This is not necessarily the case for
the anti-symmetric contribution, and this is one possible cause of the
oscillations that we see (although, as noted previously, it does not show a
clear correlation with the strength of precession). Conversely, we found that
the amplitude evolution of the anti-symmetric (2,2) multipole was monotonic
for PN waveforms if we choose $\iota=\alpha=0$ in their construction (which is
equivalent to choosing a co-precessing frame that tracks the Newtonian orbital
angular momentum, i.e., the normal to the orbital plane), but if we consider
the PN waveforms in the QA frame then the anti-symmetric (2,2) amplitude shows
strong modulations. This illustrates that the anti-symmetric component can
depend strongly on the choice of co-precessing frame, and although we do not
expect this to be a significant issue at the level of accuracy or
approximation of the current model, it may require better understanding in
future refinements of anti-symmetric models.

The second point is that our model is constructed based on the phenomenology
of single-spin binaries. If the model is to be used for generic binaries, one
much choose a method to treat two-spin configurations. One option, which is
employed in the PhenomX implementation Hamilton  _et al._ (2021); Thompson
_et al._ (2023), is as follows.

We can describe the spin of an equivalent single-spin system using Eqs. (16)
and (17) in Ref. Hamilton  _et al._ (2021), but diverge slightly in the
definition of the in-plane spin as follows:

$\chi^{\perp}=\begin{cases}\chi_{\rm{as}}\cos^{2}(\theta_{q})+\chi_{\rm{p}}\sin^{2}(\theta_{q}),&1\leq
q\leq 1.5\\\ \chi_{\rm{p}},&q>1.5\end{cases}$ (19)

where $\chi_{\rm{p}}$ is the effective precession spin as defined in Ref.
Schmidt  _et al._ (2015). The anti-symmetric amplitude for an equal-mass
binary with both spins equal in magnitude, entirely in-plane and pointing in
the same direction, must drop to zero. The symmetric in-plane spin combination
of Eq. (15) in Ref. Hamilton  _et al._ (2021) cannot account for this effect
in superkick configuration. Therefore, we can instead use an anti-symmetric
in-plane spin combination,

$\chi_{\rm{as}}=\frac{|\bf{S}^{\perp}_{1}-\bf{S}^{\perp}_{2}|}{m_{1}^{2}},$
(20)

to appropriately map two-spin to single-spin systems for generating the anti-
symmetric waveform.

Figure 9: Surface
$b(\eta,\theta_{LS})=b_{0}+b_{1}\eta+b_{2}\theta_{LS}+b_{3}\eta\theta_{LS}$
fit of the model’s coefficient, $b$, to the two-dimensional parameter space
${\eta,\theta_{LS}}$. The red points denote the $80$ computed $b$ coefficients
of the multipole asymmetry amplitude model.  Figure 10: The $b$ coefficient as
a function of the symmetric mass ratio, $\eta$, for a selected angle
$\theta_{LS}=90^{o}$ and all the available spin values,
$\chi=[0.2,0.4,0.6,0.8]$. The grey line shows the surface fit,
$b(\eta,90^{\circ})$ from Eq. (18).

##  V Phase model

GW phases for chirp signals in the frequency domain are typically quite
featureless and therefore not conducive to modeling. Following a suite of
models for the symmetric phase Husa  _et al._ (2016); Khan  _et al._ (2016);
Pratten  _et al._ (2020), we focus first on the anti-symmetric phase
derivative. Fig. 11 demonstrates that the anti-symmetric phase derivative
(i.e., frequency) behaviour in the time-domain as discussed in Sec. II.2 is
preserved in the frequency domain as well. Remarkably, we find that it is
possible to construct a map of the symmetric phase derivative to the anti-
symmetric phase derivative that is independent of the binary’s parameters.
Therefore, we do not need to produce any parametric fits for the map over the
intrinsic parameter space of BH binaries, which makes this model extremely
simple.

Figure 11: The anti-symmetric phase derivative (thick black line), the
symmetric phase derivative (grey dashed line) and the combination of half of
the symmetric phase derivative with the derivative of the precession angle
$\alpha$ (thin grey line). The specific configuration shown here is for a
single-spin binary with $(q=1$, $\chi=0.2$, $\theta_{LS}=30^{\circ})$.

Our model of the anti-symmetric phase derivative is defined by the piecewise
function,

$\phi_{a}^{\prime}(f)=\begin{cases}\frac{1}{2}\phi^{\prime}_{s}(f)+\alpha^{\prime}(f)+A_{0},&f\leq
f_{T}-\frac{f_{w}}{2}\\\
\phi^{\prime}_{\rm{int}}(f),&f_{T}-\frac{f_{w}}{2}<f\leq
f_{T}+\frac{f_{w}}{2}\\\ \phi^{\prime}_{s}(f),&f_{T}+\frac{f_{w}}{2}\leq
f<0.12,\end{cases}$ (21)

where the phase derivative in the intermediate region is given by,

$\begin{split}\phi^{\prime}_{\rm
int}(f)=\frac{1}{2}\bigg{[}1-\frac{3}{2f_{w}}((f-f_{T})-\frac{(f-f_{T})^{3}}{3f_{w}^{2}})\bigg{]}\big{(}\frac{\phi^{\prime}_{s}}{2}+\alpha^{\prime}+A_{0}\big{)}\\\
+\frac{1}{2}\bigg{[}1+\frac{3}{2f_{w}}((f-f_{T})-\frac{(f-f_{T})^{3}}{3f_{w}^{2}})\bigg{]}(\phi^{\prime}_{s}+B_{0})\end{split}$
(22)

As is evident from Fig. 11, the functional form of region I transitions to the
functional form of region II in the intermediate region. To ensure smooth
transition in the intermediate region an obvious choice would be a $tanh$
windowing function. Noting that a $tanh$ function can be computationally
inefficient during model evaluation, we instead use the Taylor expansion of
the $tanh$ function up to second order with appropriate normalization to
construct the phase derivative functional form of Eq. (22).

The phase derivative ansatz was calibrated to NR simulations by treating the
transition frequency, $f_{T}$, the width of the transition window, $f_{w}$,
and the phase coefficients ($A_{0}$ and $B_{0}$) as free parameters. The
window parameters were not particularly sensitive to the tuning and the minor
variations could be attributed to the noise in the anti-symmetric phase
derivatives; $f_{T}=0.85f_{RD}$ and $f_{w}=0.005$ was found to be an optimal
choice across the entire single-spin parameter space.

Once the parameters of the window function were fixed, we investigated the
impact of fixing the phase coefficients. A best fit to data for $B_{0}=0$ was
found to be consistent across the parameter space; $A_{0}$ on the other hand
showed some variation, but no specific trend. Furthermore, choosing the
algebraic mean of $A_{0}$ for the set of $80$ simulations did not
significantly impact the quality of the fit. This shows (1) the fitting
algorithm tries to find the best $A_{0}$ for continuity of the phase
derivative at $f_{T}$, and (2) the variation in optimal $A_{0}$ across the
parameter space is more likely due to noisy data and not a function of
intrinsic parameters of the binary.

Applying a shift to the phase derivative is equivalent to an overall time
shift of the waveform. We exploit this freedom by fixing the symmetric phase
derivative minima to be $0$ at $f_{RD}$. This imposes,

$A_{0}=\frac{1}{2}\phi^{\prime}_{s}(f_{T})-\alpha^{\prime}(f_{T}).$ (23)

Figure 12 shows that $A_{0}$ obtained from the fitting algorithm and from NR
data using Eq. (23) reasonably well for a majority of the cases (cf. y-axis on
Fig. 11). The anti-symmetric model will be used with a phenomenological
symmetric waveform model, so an $A_{0}$ derived from the symmetric waveform
model makes the phase construction self-consistent and robust. Furthermore,
Fig. 12 highlights that the gain in accuracy by making a model to capture the
near-stochastic behaviour of $A_{0}$ may be outweighed by errors introduced by
over-modelling. As such, our model of the anti-symmetric phase is NR-informed
but, noting that the data for the anti-symmetric waveform is often close to
numerical noise, we prioritised our understanding of the physics rather than
model optimization.

Figure 12: Comparison of phase coefficient $A_{0}$, obtained from fitting the
phase ansatz to NR data (light purple triangles) across the parameter space
and from NR data using Eq. (23) (blue circles), for cases with $\chi=0.4$ and
0.8. The algebraic mean of the set of coefficients are shown by horizontal
lines in corresponding colors. The case index is as described in Sec. III.

The phase of the anti-symmetric waveform is obtained by integrating the two
pieces,

$\phi_{a}(f)=\begin{cases}\frac{1}{2}\phi_{s}(f)+\alpha(f)+A_{0}f+\phi_{A0},&f\leq
f_{T}\\\ \phi_{s}(f)+\phi_{B0},&f_{T}\leq f<0.12,\end{cases}$ (24)

where the integration constant $\phi_{B0}$ is determined by continuity at
$f_{T}$

$\phi_{B0}=\alpha(f_{T})-\phi_{s}(f_{T})+A_{0}f_{T}.$ (25)

Finally, the phase of the asymmetry is modulated by the in-plane spin
direction $\alpha$; therefore, the initial phase, $\phi_{A0}$, is the value of
$\alpha$ at a reference frequency i.e., $\phi_{A0}=\alpha(f_{\rm{ref}})$.

##  VI Model accuracy

A standard measure of waveform accuracy used extensively in the literature is
the match of the waveform model with NR waveforms, defined as,

$M(h_{NR},h_{M})=4\,\rm{Re}\int\limits_{0}^{\infty}\frac{h_{NR}(f)h^{*}_{M}(f)}{S_{n}(f)}df,$
(26)

where $h(f)=h_{+}(f)-ih_{\times}(f)$ is a complex frequency sequence
constructed from the two polarizations of the waveform. As such, calculating
matches of just the anti-symmetric waveform is not physically meaningful.
Furthermore, a true measure of performance of precessing waveforms in data
analysis can only be obtained by calculating matches of the full waveform in
the inertial frame, making considerations for precession as well as extrinsic
parameters. Therefore, matches of just the anti-symmetric waveform in the co-
precessing frame provide some indication of the accuracy of one ingredient in
the full waveform, but do not indicate the overall accuracy of the
corresponding precessing waveform.

Figure 13: Mismatches of the anti-symmetric waveform model in the co-
precessing frame with NR data. The black triangles show the mismatches for the
combined amplitude and phase model while the magenta squares (cyan circles)
show the mismatches for just the amplitude (phase) model with the phase
(amplitude) constructed from NR data. The dashed magenta and cyan lines show
the average mismatch for the amplitude and phase model, respectively; on
average they are of comparable accuracy. The black solid line shows the
average mismatch for the overall model. The x-axis denotes the case index of
the NR simulation as usual i.e., five different $\theta_{\rm{LS}}$ for each
spin magnitude shown in figure, for $q=1$,$2$,$4$ and $8$.

However, an inner product like that in Eq. (26) is a useful measure of
agreement between two complex frequency series. Since a match-like calculation
for the anti-symmetric contribution in the co-precessing frame cannot be
interpreted in terms of either signal detection efficiency or parameter
measurement accuracy, there is no reason to include the detector sensitivity,
and so we will use a simpler inner product of the form,

$\langle
h_{NR}|h_{M}\rangle=\rm{Re}\int\limits_{f_{1}}^{f_{2}}h_{NR}(f)h^{*}_{M}(f)df,$
(27)

where $f_{1}$ is the starting frequency of the NR waveform in geometric units,
$Mf_{1}=0.02$, and $Mf_{2}=0.15$, after which point the amplitude of the NR
waveform is below the noise floor of the data. We consider normalised
waveforms, $\hat{h}=h/\sqrt{\langle h|h\rangle}$, so that the maximum value of
the inner product is one. We used the standard implementation of this inner
product in pycbc Nitz  _et al._ (2021), a python software package for GW data
analysis, for our match computations. We then consider the “mismatch” between
the two waveforms,

$\mathcal{M}=1-\frac{\langle h_{NR}|h_{M}\rangle}{\sqrt{\langle
h_{NR}|h_{NR}\rangle\langle h_{M}|h_{M}\rangle}}.$ (28)

In Fig. 13 we show the mismatches of the anti-symmetric waveform constructed
from our model with the 80 NR waveforms that were used to calibrate the model.
To determine the accuracy of the individual components, we also computed
matches of the amplitude (phase) model complemented by phase (amplitude)
constructed from NR data, with the full anti-symmetric waveform constructed
from NR data. The overall accuracy of both the models are comparable. We note
that although the model was verified using the same waveforms as used for
modelling, since the NR tuning was relatively simple — i.e., a single co-
efficient in Eq. (16) fit to the four-parameter ansatz Eq. (18) across a two-
dimensional parameter space — using a much smaller subset of waveforms would
have produced a model of similar accuracy, and the simplicity of this model
and the single-spin parameter space obviates any concerns about over-fitting
or unexplored corners of parameter space.

To investigate the quality of the surface fit in Eq. (18), we also computed
mismatches for the amplitude model using fit coefficients $b$ of Eq. (16). As
is evident from Fig. 14, for most cases the performance is unchanged and for
the handful of cases where the mismatch changes, the difference is not very
significant. This further illustrates that capturing the non-linear dependence
on spin magnitude is unlikely to make significant improvement to the amplitude
model. In addition to the argument made for using Eq. (23) to calculate
$A_{0}$ from the symmetric waveform and precession angle, $\alpha$, we
calculated mismatches for the different choices of $A_{0}$ in the phase model
– i.e., $A^{fit}_{0}$and $A^{equation}_{0}$ in Fig. 12 – to confirm that there
was no reduction in model accuracy.

Figure 14: Mismatches (same as Fig. 13) showing comparison of the amplitude
model constructed using the spin magnitude-independent surface fit of Eq. (16)
(blue squares) with the amplitude model constructed from the true fit
coefficients (green circles).

Note that the anti-symmetric waveform model is downstream from the symmetric
waveform model as well as the precession angle models. Therefore, enhancement
in performance of the overall model due to the addition of an asymmetry model
must always be discussed in the context of the underlying symmetric,
precessing waveform model. This is beyond the scope of present work and will
be discussed in the context of the current generation frequency-domain
precessing-binary PhenomXPNR model Thompson  _et al._ (2023).

##  VII Conclusion

We have presented a method to model the anti-symmetric contribution to the
$(\ell=2,|m|=2)$ multipoles in the co-precessing frame. This is a general
approach that can be applied to any frequency-domain model that provides the
symmetric contribution to the (2,2) amplitude, phase and precession angle,
$\alpha$. We expect that the main insights of this model can be used to easily
generalise the procedure to the time-domain and be useful in including
multipole asymmetry in current generation time-domain models Estellés  _et
al._ (2022b); Ramos-Buades  _et al._ (2023). We summarise the key insights
here.

For the amplitude, we observe that the anti-symmetric amplitude can be easily
modeled as a rescaling of the symmetric amplitude. As shown in Fig. 4, both
amplitudes have the same basic structure, in particular the same ringdown
frequency and decay rate. In the inspiral our model of the amplitude ratio is
based on a PN expression for single-spin systems, and in the late inspiral and
merger a higher-order correction to the PN expression is calibrated to 80 NR
simulations of single-spin binaries Hamilton  _et al._ (2023). In the ringdown
we make use of the prediction from perturbation theory that the ratio of the
symmetric and anti-symmetric amplitudes will be constant. The amplitude model
ratio is presented in Eqs. (15), (16) and the fit to the higher-order
contribution given in Eq. (18).

For the phase, we use the facts that in the inspiral the frequency of the
anti-symmetric contribution equals the orbital frequency plus the spin
precession frequency, and that during the ringdown the symmetric and anti-
symmetric frequencies are the same. We are able to construct a mapping from
the symmetric phase and precession angle, $\phi_{s}$ and $\alpha$, to the
anti-symmetric phase, $\phi_{a}$, motivated by the 80 single-spin simulations,
but without the need of any explicit tuning. The model of the phase is given
in Eqs. (23), (24) and (25). An additional crucial observation is that a
rotation of the initial in-plane spin direction of $\Delta\alpha$ introduces a
corresponding shift of $\Delta\alpha$ into the anti-symmetric phase. It is
this observation that allows us to model the dependence on in-plane-spin
direction, even though we do not have a set of NR simulations that span this
subspace.

The procedure presented here is not in itself a signal model. As already
noted, one must also provide the symmetric amplitude, phase and precession
angle. In addition, having constructed a model for the anti-symmetric
contribution to the (2,2) multipoles in the co-precessing frame, one must then
“twist them up” to produce the multipoles in the inertial frame. This has been
done for the recent extension of the multi-mode frequency-domain precessing-
binary IMRPhenomXPNR model Thompson  _et al._ (2023).

This is the first phenomenological full inspiral-merger-ringdown model of the
anti-symmetric multipole contributions, and there are many directions for
improvement and issues to be resolved. The first limitation of the model is
that it is based only on single-spin binaries. We expect that generic two-spin
systems can in most cases be modelled to a good approximation by single-spin
systems, but given that the anti-symmetric contribution is itself a weak
effect, it is unclear how well this approximation can be used. It would be
useful to study the applicability of the single-spin approximation for the
anti-symmetric contribution, and, indeed, to extend the model to two-spin
systems. The anti-symmetric model is also limited to the $(\ell=2,|m|=2)$
multipoles. We argue in Fig. 1 that the (2,2) contribution will be sufficient
for most applications, since the anti-symmetric amplitude is far weaker than
the symmetric, but for high signal-to-noise-ratio (SNR) systems, higher-order
anti-symmetric multipoles will ultimately be required.

When decomposing the signal into a co-precessing frame and corresponding time-
or frequency-dependent angles, one will find (depending on the co-precessing
frame chosen) oscillations in the amplitude and/or phase of the anti-symmetric
contribution; recall that the co-precessing-frame maps exactly to a non-
precessing aligned-spin waveform only for the leading-order quadrupole terms
Schmidt  _et al._ (2012); Hamilton  _et al._ (2021); Thompson  _et al._
(2023). In this work we have removed oscillations in the inspiral amplitude by
simply setting the inclination angle $\iota$ and precession angle $\alpha$ to
zero in the PN expressions we have used for the signal multipoles Arun  _et
al._ (2009); oscillations remain in the co-precessing-frame signals for the NR
waveforms, but our model consists of only a monotonic fit through these
oscillations. A better understanding of these oscillations, or at least a
model that captures them, is a necessary next step in modelling the multipole
asymmetries.

We have assumed that the only affect of the initial in-plane spin direction is
to introduce an overall offset in the anti-symmetric phase. This approximation
will not be exact during the merger-ringdown. In particular, in the ringdown
we have made the approximation that the amplitude ratio will be independent of
the initial in-plane spin direction. However, the amplitude ratio will be
determined by the relative phase of the symmetric and anti-symmetric
contributions when ringdown begins, and we have not attempted to model this
effect.

These limitations aside, we find that our method to construct the anti-
symmetric contribution has an average mismatch error better than 0.03. Given
that the anti-symmetric contribution is in general less than 10% of the total
SNR, we expect that its contribution to the total mismatch uncertainty of a
model will be less than $3\times 10^{-3}$, which is below the average mismatch
error of current Phenom and SEOBNR precessing-binary models Ramos-Buades  _et
al._ (2023); Thompson  _et al._ (2023). With the addition of this first anti-
symmetric model, we expect that current models will be able to make more
precise measurements of black-hole spins and gravitational recoil. The
accuracy of these measurements will depend significantly on the underlying
symmetric model, as well as where a signal lies in the binary parameter space,
and we leave such accuracy studies with individual models to future work.

##  VIII Acknowledgements

We thank Jonathan Thompson, Eleanor Hamilton and Lionel London for many
thought-provoking discussions on generic precessing binaries. We are indebted
to Lionel London for sharing his numerical relativity data utility package
with us, this was used to extract the anti-symmetric waveform data in the co-
precessing frame. We thank Frank Ohme and Jannik Mielke for insightful
discussions on multipole asymmetry, and Chinmay Kalaghatgi for discussions on
in-plane spin rotation in superkick configurations.

S.G thanks Sebastian Khan for stimulating discussions on frequency-domain
waveform modeling, Duncan Macleod for excellent inputs on computing practices
and Edward Fauchon-Jones for support in using supercomputing facilities.

The authors were supported in part by Science and Technology Facilities
Council (STFC) grant ST/V00154X/1 and European Research Council (ERC)
Consolidator Grant 647839. S.G acknowledges support from the Max Planck
Society’s Independent Research Group program. P.K was also supported by the GW
consolidated grant: STFC grant ST/V005677/1.

Simulations used in this work were performed on the DiRAC@Durham facility,
managed by the Institute for Computational Cosmology on behalf of the STFC
DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital
funding via STFC capital grants ST/P002293/1 and ST/R002371/1, Durham
University and STFC operations grant ST/R000832/1. In addition, several of the
simulations used in this work were performed as part of an allocation
graciously provided by Oracle to explore the use of our code on the Oracle
Cloud Infrastructure.

This research also used the supercomputing facilities at Cardiff University
operated by Advanced Research Computing at Cardiff (ARCCA) on behalf of the
Cardiff Supercomputing Facility and the HPC Wales and Supercomputing Wales
(SCW) projects. We acknowledge the support of the latter, which is part-funded
by the European Regional Development Fund (ERDF) via the Welsh Government. In
part the computational resources at Cardiff University were also supported by
STFC grant ST/I006285/1.

Various plots and analyses in this paper were made using Python software
packages LALSuite LIGO Scientific Collaboration (2018), PyCBC Nitz  _et al._
(2021), gwsurrogate Field  _et al._ (2014), Matplotlib Hunter (2007), Numpy
Harris  _et al._ (2020), lmfit and Scipy Virtanen  _et al._ (2020).

## References

  * Abbott  _et al._ (2019a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019a), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 818, L22 (2016b), arXiv:1602.03846 [astro-ph.HE] .
  * Abbott  _et al._ (2019b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 882, L24 (2019b), arXiv:1811.12940 [astro-ph.HE] .
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 913, L7 (2021c), arXiv:2010.14533 [astro-ph.HE] .
  * Abbott  _et al._ (2023) R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13, 011048 (2023), arXiv:2111.03634 [astro-ph.HE] .
  * Pürrer  _et al._ (2016) M. Pürrer, M. Hannam, and F. Ohme, Phys. Rev. D 93, 084042 (2016), arXiv:1512.04955 [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] .
  * 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] .
  * London  _et al._ (2018) L. London, S. Khan, E. Fauchon-Jones, C. García, M. Hannam, S. Husa, X. Jiménez-Forteza, C. Kalaghatgi, F. Ohme, and F. Pannarale, Phys. Rev. Lett. 120, 161102 (2018), arXiv:1708.00404 [gr-qc] .
  * Khan  _et al._ (2019) S. Khan, K. Chatziioannou, M. Hannam, and F. Ohme, Phys. Rev. D100, 024059 (2019), arXiv:1809.10113 [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] .
  * 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] .
  * Pratten  _et al._ (2021) G. Pratten _et al._ , Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Thompson  _et al._ (2020) J. E. Thompson, E. Fauchon-Jones, S. Khan, E. Nitoglia, F. Pannarale, T. Dietrich, and M. Hannam, Phys. Rev. D 101, 124059 (2020), arXiv:2002.08383 [gr-qc] .
  * Estellés  _et al._ (2021) H. Estellés, A. Ramos-Buades, S. Husa, C. García-Quirós, M. Colleoni, L. Haegel, and R. Jaume, Phys. Rev. D 103, 124060 (2021), arXiv:2004.08302 [gr-qc] .
  * Estellés  _et al._ (2022a) H. Estellés, S. Husa, M. Colleoni, D. Keitel, M. Mateu-Lucena, C. García-Quirós, A. Ramos-Buades, and A. Borchers, Phys. Rev. D 105, 084039 (2022a), arXiv:2012.11923 [gr-qc] .
  * Taracchini  _et al._ (2012) A. Taracchini, Y. Pan, A. Buonanno, E. Barausse, M. Boyle, T. Chu, G. Lovelace, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 86, 024011 (2012), arXiv:1202.0790 [gr-qc] .
  * Pan  _et al._ (2014) Y. Pan, A. Buonanno, A. Taracchini, L. E. Kidder, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Phys. Rev. D89, 084006 (2014), arXiv:1307.6232 [gr-qc] .
  * Taracchini  _et al._ (2014) A. Taracchini _et al._ , Phys. Rev. D89, 061502 (2014), arXiv:1311.2544 [gr-qc] .
  * Bohé  _et al._ (2017) A. Bohé _et al._ , Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Cotesta  _et al._ (2018) R. Cotesta, A. Buonanno, A. Bohé, A. Taracchini, I. Hinder, and S. Ossokine, Phys. Rev. D 98, 084028 (2018), arXiv:1803.10701 [gr-qc] .
  * Ossokine  _et al._ (2020) S. Ossokine _et al._ , Phys. Rev. D 102, 044055 (2020), arXiv:2004.09442 [gr-qc] .
  * Matas  _et al._ (2020) A. Matas _et al._ , Phys. Rev. D 102, 043023 (2020), arXiv:2004.10001 [gr-qc] .
  * Apostolatos  _et al._ (1994) T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D 49, 6274 (1994).
  * Bruegmann  _et al._ (2008a) B. Bruegmann, J. A. Gonzalez, M. Hannam, S. Husa, and U. Sperhake, Phys. Rev. D 77, 124047 (2008a), arXiv:0707.0135 [gr-qc] .
  * Hannam  _et al._ (2022) M. Hannam _et al._ , Nature 610, 652 (2022), arXiv:2112.11300 [gr-qc] .
  * Varma  _et al._ (2022) V. Varma, S. Biscoveanu, T. Islam, F. H. Shaik, C.-J. Haster, M. Isi, W. M. Farr, S. E. Field, and S. Vitale, Phys. Rev. Lett. 128, 191102 (2022), arXiv:2201.01302 [astro-ph.HE] .
  * Kalaghatgi and Hannam (2021) C. Kalaghatgi and M. Hannam, Phys. Rev. D 103, 024024 (2021), arXiv:2008.09957 [gr-qc] .
  * Ramos-Buades  _et al._ (2020) A. Ramos-Buades, P. Schmidt, G. Pratten, and S. Husa, Phys. Rev. D 101, 103014 (2020), arXiv:2001.10936 [gr-qc] .
  * Kolitsidou  _et al._ (2023) P. Kolitsidou, J. E. Thompson, and M. Hannam, (2023), in preparation.
  * Varma  _et al._ (2019) V. Varma, S. E. Field, M. A. Scheel, J. Blackman, D. Gerosa, L. C. Stein, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. Research. 1, 033015 (2019), arXiv:1905.09300 [gr-qc] .
  * Hamilton  _et al._ (2023) E. Hamilton _et al._ , (2023), arXiv:2303.05419 [gr-qc] .
  * Kidder (1995) L. E. Kidder, Phys. Rev. D 52, 821 (1995), arXiv:gr-qc/9506022 .
  * Gonzalez  _et al._ (2007) J. A. Gonzalez, M. D. Hannam, U. Sperhake, B. Bruegmann, and S. Husa, Phys. Rev. Lett. 98, 231101 (2007), arXiv:gr-qc/0702052 .
  * Campanelli  _et al._ (2007) M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Merritt, Phys. Rev. Lett. 98, 231102 (2007), arXiv:gr-qc/0702133 .
  * Keppel  _et al._ (2009) D. Keppel, D. A. Nichols, Y. Chen, and K. S. Thorne, Phys. Rev. D 80, 124015 (2009), arXiv:0902.4077 [gr-qc] .
  * Boyle  _et al._ (2014) M. Boyle, L. E. Kidder, S. Ossokine, and H. P. Pfeiffer, (2014), arXiv:1409.4431 [gr-qc] .
  * Schmidt  _et al._ (2012) P. Schmidt, M. Hannam, and S. Husa, Phys. Rev. D 86, 104063 (2012), arXiv:1207.3088 [gr-qc] .
  * Estellés  _et al._ (2022b) H. Estellés, M. Colleoni, C. García-Quirós, S. Husa, D. Keitel, M. Mateu-Lucena, M. d. L. Planas, and A. Ramos-Buades, Phys. Rev. D 105, 084040 (2022b), arXiv:2105.05872 [gr-qc] .
  * Ramos-Buades  _et al._ (2023) A. Ramos-Buades, A. Buonanno, H. Estellés, M. Khalil, D. P. Mihaylov, S. Ossokine, L. Pompili, and M. Shiferaw, (2023), arXiv:2303.18046 [gr-qc] .
  * Arun  _et al._ (2009) K. G. Arun, A. Buonanno, G. Faye, and E. Ochsner, Phys. Rev. D 79, 104023 (2009), [Erratum: Phys.Rev.D 84, 049901 (2011)], arXiv:0810.5336 [gr-qc] .
  * Hamilton  _et al._ (2021) E. Hamilton, L. London, J. E. Thompson, E. Fauchon-Jones, M. Hannam, C. Kalaghatgi, S. Khan, F. Pannarale, and A. Vano-Vinuales, Phys. Rev. D 104, 124027 (2021), arXiv:2107.08876 [gr-qc] .
  * Thompson  _et al._ (2023) J. E. Thompson, E. Hamilton, L. London, S. Ghosh, P. Kolitsidou, C. Hoy, and M. Hannam, (2023), in preparation.
  * Bruegmann  _et al._ (2008b) B. Bruegmann, J. A. Gonzalez, M. Hannam, S. Husa, U. Sperhake, and W. Tichy, Phys. Rev. D 77, 024027 (2008b), arXiv:gr-qc/0610128 .
  * Husa  _et al._ (2008) S. Husa, J. A. Gonzalez, M. Hannam, B. Bruegmann, and U. Sperhake, Class. Quant. Grav. 25, 105006 (2008), arXiv:0706.0740 [gr-qc] .
  * London (2015) L. London, (2015).
  * London  _et al._ (2020) L. London, E. Fauchon-Jones, and E. Hamilton, “llondon6/positive: map,”  (2020).
  * Schmidt  _et al._ (2011) P. Schmidt, M. Hannam, S. Husa, and P. Ajith, Phys. Rev. D84, 024046 (2011), arXiv:1012.2879 [gr-qc] .
  * O’Shaughnessy  _et al._ (2011) R. O’Shaughnessy, B. Vaishnav, J. Healy, Z. Meeks, and D. Shoemaker, Phys. Rev. D 84, 124002 (2011), arXiv:1109.5224 [gr-qc] .
  * Boyle  _et al._ (2011) M. Boyle, R. Owen, and H. P. Pfeiffer, Phys. Rev. D 84, 124011 (2011), arXiv:1110.2965 [gr-qc] .
  * Schmidt  _et al._ (2015) P. Schmidt, F. Ohme, and M. Hannam, Phys. Rev. D 91, 024043 (2015), arXiv:1408.1810 [gr-qc] .
  * Nitz  _et al._ (2021) A. Nitz, I. Harry, D. Brown, C. M. Biwer, J. Willis, T. D. Canton, C. Capano, T. Dent, L. Pekowsky, A. R. Williamson, G. S. C. Davies, S. De, M. Cabero, B. Machenschalk, P. Kumar, D. Macleod, S. Reyes, dfinstad, F. Pannarale, T. Massinger, S. Kumar, M. Tápai, L. Singer, S. Khan, S. Fairhurst, A. Nielsen, S. Singh, shasvath, B. U. V. Gadre, and I. Dorrington, “gwastro/pycbc: Pycbc release 1.18.1,”  (2021).
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL) (2018).
  * Field  _et al._ (2014) S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, and M. Tiglio, Phys. Rev. X 4, 031006 (2014), arXiv:1308.3565 [gr-qc] .
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * 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 585, 357 (2020).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).