# A Survey of Four Precessing Waveform Models for Binary Black Hole Systems

Jake Mac Uilliam1 Sarp Akçay1 Jonathan E. Thompson2 1University College
Dublin, Belfield, D4, Dublin, Ireland  2Theoretical Astrophysics Group,
California Institute of Technology, Pasadena, CA 91125, U.S.A.

###### Abstract

Angular momentum and spin precession are expected to be generic features of a
significant fraction of binary black hole systems. As such, it is essential to
have waveform models that faithfully incorporate the effects of precession.
Here, we assess how well the current state of the art models achieve this for
waveform strains constructed only from the $\ell=2$ multipoles. Specifically,
we conduct a survey on the faithfulness of the waveform models SEOBNRv5PHM,
TEOBResumS, IMRPhenomTPHM, IMRPhenomXPHM to the numerical relativity (NR)
surrogate NRSur7dq4 and to NR waveforms from the SXS catalog. The former
assessment involves systems with mass ratios up to six and dimensionless spins
up to 0.8. The latter employs $317$ short and $23$ long SXS waveforms. For all
cases, we use reference inclinations of zero and $90^{\circ}$. We find that
all four models become more faithful as the mass ratio approaches unity and
when the merger-ringdown portion of the waveforms are excluded. We also
uncover a correlation between the co-precessing $(2,\pm 2)$ multipole
mismatches and the overall strain mismatch. We additionally find that for high
inclinations, precessing $(2,\pm 1)$ multipoles that are more faithful than
their $(2,\pm 2)$ counterparts, and comparable in magnitude, improve waveform
faithfulness. As a side note, we show that use of uniformly-filled parameter
spaces may lead to an overestimation of precessing model faithfulness. We
conclude our survey with a parameter estimation study in which we inject two
precessing SXS waveforms (at low and high masses) and recover the signal with
SEOBNRv5PHM, IMRPhenomTPHM and IMRPhenomXPHM. As a bonus, we present
preliminary multidimensional fits to model unfaithfulness for Bayesian model
selection in parameter estimation studies.

###### pacs:

04.25.D-, 04.30.Db, 04.30.-w, 04.80.Nn, 04.25.D-, 04.25.dg 04.25.Nx, 95.30.Sf,
97.60.Jd

###### Contents

  1. I Introduction
  2. II Review of Precessing Waveform Construction
  3. III Quantifying Waveform Faithfulness
  4. IV Faithfulness Survey I: Comparisons with NRSur7dq4
    1. IV.1 Discrete Parameter Set with zero Inclination
      1. IV.1.1 Effects of the Merger-Ringdown Portion of GWs on Faithfulness
      2. IV.1.2 The Dependence of Waveform Unfaithfulness on the Mass Ratio
      3. IV.1.3 The Effect of AS/Co-precessing multipoles on Faithfulness
      4. IV.1.4 Mismatches in the Extrapolation Region of NRSur7dq4
    2. IV.2 Behavior of Model Unfaithfulness under changing Inclination
    3. IV.3 Comparisons over a uniformly filled Parameter Space
  5. V Faithfulness Survey II: Comparisons with Numerical Relativity Waveforms
    1. V.1 Short SXS Waveforms
      1. V.1.1 Effects of Merger-Ringdown on Faithfulness
      2. V.1.2 Dependence of Inspiral-only Faithfulness on Intrinsic Parameters
      3. V.1.3 Model Faithfulness under increasing Inclination
      4. V.1.4 Effects of Co-precessing Multipoles on Faithfulness
      5. V.1.5 $Q\geq 5$ Comparisons
    2. V.2 Long SXS Waveforms
  6. VI Injection/Recovery Study
    1. VI.1 Main Results
    2. VI.2 Inferring Precession from the Injections
  7. VII Summary
  8. VIII Discussion
  9. A Multidimensional Fits of Mismatches for Bayesian Model Selection
  10. B A brief Comparison of SEOBNRv5PHM with SEOBNRv4PHM
  11. C Another Comparison: IMRPhenomXO4a vs. MSA and SpinTaylor Implementations of IMRPhenomXPHM
  12. D Timing Benchmarks

##  I Introduction

Over the course of three observing periods, the terrestrial network of
gravitational wave interferometers has detected roughly 100 compact binary
inspiral-merger events Abbott  _et al._ (2016a, b, c, d, e, f, g, h, 2017a,
2017b, 2017c, 2017d, 2017e, 2019); Venumadhav  _et al._ (2019, 2020); Nitz
_et al._ (2019a, 2020); Abbott  _et al._ (2020a, b, c, 2021a); Nitz  _et al._
(2021); Abbott  _et al._ (2017f, 2020d, 2020e, 2021b); Olsen  _et al._ (2022);
Williams (2024). An overwhelming majority of these involved binary black hole
systems where each compact body is expected to have non-negligible spin
angular momentum, which has been confirmed by the analysis of the binary black
hole population thus far detected Abbott  _et al._ (2019, 2021c, 2021a).
Depending on the binary formation scenario, a significant subpopulation of
these binaries can have spins misaligned with the orbital angular momentum
Mandel and Farmer (2022); Mapelli (2021); Gerosa  _et al._ (2013); Vitale  _et
al._ (2017); Rodriguez  _et al._ (2016); Stevenson  _et al._ (2017).

The leading-order general relativistic effect of such misaligned spins is the
precession of the orbital angular momentum vector around the total angular
momentum of the system, which leaves a faint, but detectable imprint on the
gravitational waveform as amplitude and phase modulations on a timescale 1.5
post-Newtonian (PN) orders longer than the orbital timescale Apostolatos  _et
al._ (1994); Kidder (1995), and at an order 1.5PN higher in the phase than the
leading order terms Cutler and Flanagan (1994); Poisson and Will (1995). As
such, spin effects are harder to infer, or even detect, than the chirp mass
and the [symmetric] mass ratio of the binary system. Nonetheless, information
about black hole spins has been obtained successfully Abbott  _et al._ (2016i,
f, j, 2017g, 2017e, 2019); Nitz  _et al._ (2019a); Venumadhav  _et al._
(2020); Abbott  _et al._ (2020a, b, e, 2021c, d, 2021a, 2021d); Nitz  _et al._
(2021); Olsen  _et al._ (2022). One way to convey this information is to
provide distributions for the inferred spin magnitudes and tilts, the latter
with respect to some reference frame. However, thus far, such posteriors have
been mostly uninformative with a few exceptions as illustrated by Figs. 6, 10,
and 11 of the three Gravitational Wave Transient Catalogs, respectively Abbott
_et al._ (2019, 2021c, 2021a).

A more fruitful way of gleaning spin information has been via the construction
of specific projections of the spins parallel and perpendicular to the orbital
angular momentum at a reference frequency, usually 20 Hz. These projections
reduce the seven-dimensional intrinsic parameter space (mass ratio and six
spin components) to three dimensions, which was shown to essentially capture
the phenomenology of precessing waveforms Hannam  _et al._ (2014). The
parallel scalar has come to be known as the effective spin parameter
$\chi_{\text{eff}}$ Damour (2001); Racine (2008). Though the term effective
(or reduced or PN) spin has been used for similar scalars Ajith (2011);
Schmidt  _et al._ (2011), we employ the phenomenological definition given in
our Eq. (8) Hannam  _et al._ (2014); Pürrer  _et al._ (2013) which has become
standard. It was shown in Ref. Racine (2008) that $\chi_{\text{eff}}$ is a
conserved quantity up to 1.5PN, i.e., neglecting spin-spin and higher-order
interactions. The magnitude of $\chi_{\text{eff}}$ changes very little even
with the inclusion of these interactions. As such, it has become a very useful
quantity in parameter estimation, especially since it emerges in the PN series
for the waveform phase at a more dominant order than all other parameters
except for the chirp mass and the symmetric mass ratio Cutler and Flanagan
(1994); Poisson and Will (1995); Baird  _et al._ (2013). In fact, in the
strong field regime, the contribution of the $\chi_{\text{eff}}$ term is large
enough to cause a well known partial degeneracy between the symmetric mass
ratio and the parallel component of the total spin Cutler and Flanagan (1994);
Poisson and Will (1995); Baird  _et al._ (2013); Hannam  _et al._ (2013);
Pürrer  _et al._ (2013).

$\chi_{\text{eff}}$ also provides a way to gather information about the
properties of the binary black hole (BBH) population since isotropically
distributed spins would result in a normal distribution for
$\chi_{\text{eff}}$ centered on zero. Several studies have already looked at
this over the ensemble of the detected BBHs and more or less agree that the
distribution is somewhat asymmetric and peaks at values slightly above zero
Galaudage  _et al._ (2021); Golomb and Talbot (2022); Adamcewicz and Thrane
(2022); Tong  _et al._ (2022); Biscoveanu  _et al._ (2022); Callister  _et
al._ (2022); Baibhav  _et al._ (2022).

The perpendicular scalar is known as the effective precession parameter (or
spin) $\chi_{\text{p}}$ Hannam  _et al._ (2014); Schmidt  _et al._ (2015) and
is given below in Eq. (9). It was shown by Refs. Hannam  _et al._ (2014);
Schmidt  _et al._ (2015) that the dominant effects of precession on the GW
waveform can be characterized by mapping the four perpendicular components of
the binary’s spin vectors to just one parameter, $\chi_{\text{p}}$, assigned
to be the sole perpendicular component of the spin of the larger (primary)
component of the binary. Since $\chi_{\text{p}}$ provides information about
the perpendicular component of the spins, specifically, the perpendicular
component of the primary’s spin, it is taken as an indication of spin
precession.

The inference of $\chi_{\text{p}}$ is more challenging than that of
$\chi_{\text{eff}}$ because, for instance, in the limit that the orbital
angular momentum is much larger than the total spin, the perpendicular
components can be considered to be 0.5PN order higher than the parallel
components Arun  _et al._ (2009). This can also be understood in terms of the
respective contributions of the parallel and the perpendicular components to
the waveform phase. Comparing the contribution of the $\chi_{\text{eff}}$-term
to the overall phase with that of the accumulated precession phase (cf. Eq.
(45) of Ref. Apostolatos  _et al._ (1994)) for a canonical neutron star binary
inspiralling from 10 Hz, we see that the former is roughly two orders of
magnitude larger. Accordingly, at moderate signal-to-noise ratios (SNRs), it
is not possible to obtain a clear measurement of the individual spins Schmidt
_et al._ (2015). Moreover, the inclination of the orbit also matters; it is
harder to detect precession in binaries where the line of sight is aligned
with the total angular momentum as the precession-induced modulations are
minimized Schmidt  _et al._ (2012). Therefore, it is not surprising that there
have thus far been only three GW events for which there is evidence for
precession: GW190412 Abbott  _et al._ (2020b), GW190521 Abbott  _et al._
(2020c, e); Miller  _et al._ (2023), and GW200129_065458 Abbott  _et al._
(2021a); Hannam  _et al._ (2022); Payne  _et al._ (2022); Macas  _et al._
(2023).

Discerning precession in compact binary inspirals (and mergers) is important
for breaking parameter degeneracies and understanding binary formation
scenarios. To this end, accurate modeling of spin effects has become crucial
in GW astronomy. Though nonprecessing, i.e., aligned-spin, quasi-circular
waveform models such as those given by Refs. Cotesta  _et al._ (2018); Varma
_et al._ (2019a); Estellés  _et al._ (2022a); Pratten  _et al._ (2020);
García-Quirós  _et al._ (2020); Riemenschneider  _et al._ (2021); Pompili  _et
al._ (2023); Nagar  _et al._ (2023) have matured to the faithfulness level of
$\lesssim 10^{-3}$ with respect to numerical relativity waveforms111We quote
the approximate median value for a specific mismatch with respect to numerical
relativity simulations over a large sample., the precessing models are
approximately one half to an order of magnitude worse in faithfulness and much
more sensitive to modeling systematics. As the detectors’ sensitivities
improve steadily through O4, O5 and beyond, systematic errors in parameter
estimation due to mismodelling the effects of precession will dominate over
the statistical error, especially for signals with large SNR Cutler and
Flanagan (1994); Flanagan and Hughes (1998a). In fact, recent research has
shown that even the current accuracy of the NR waveforms may not be enough for
bias-free parameter estimation at the projected sensitivities of the third
generation ground-based interferometers Pürrer and Haster (2020). For this
reason, the waveform modeling community has been endeavoring toward building
ever more faithful precessing models. There are now several precessing quasi-
spherical waveform models borne out of various waveform “families” and all
have achieved $\lesssim 10^{-2}$ faithfulness. The state of the art among
these is NRSur7dq4 Varma  _et al._ (2019b), the latest precessing member of
the numerical relativity surrogate family Blackman  _et al._ (2014, 2017);
Varma  _et al._ (2019a); Williams  _et al._ (2019); Walker  _et al._ (2022);
Yoo  _et al._ (2023) with an NR faithfulness of $\sim 10^{-4}$. It was
recently employed in a re-analysis of most O3 events Islam  _et al._ (2023).

Another set of models is the NR-informed effective-one-body (EOB) family
Buonanno and Damour (1999, 2000); Damour  _et al._ (2000); Damour (2001);
Buonanno  _et al._ (2006); Damour  _et al._ (2015). This contains two similar,
but distinct subfamilies: SEOBNR Bohé  _et al._ (2017); Babak  _et al._
(2017); Cotesta  _et al._ (2018); Hinderer  _et al._ (2016); Steinhoff  _et
al._ (2016); Lackey  _et al._ (2019); Matas  _et al._ (2020); Ossokine  _et
al._ (2020); Gadre  _et al._ (2022); Thomas  _et al._ (2022); Ramos-Buades
_et al._ (2023); van de Meent  _et al._ (2023); Pompili  _et al._ (2023);
Khalil  _et al._ (2023); Mihaylov  _et al._ (2023) and TEOBResumS Damour and
Nagar (2014); Bernuzzi  _et al._ (2015); Nagar  _et al._ (2017, 2018); Akcay
_et al._ (2019, 2021); Nagar  _et al._ (2020); Chiaramello and Nagar (2020);
Gamba  _et al._ (2022); Nagar and Rettegno (2021); Albertini  _et al._ (2021);
Nagar and Albanesi (2022); Albertini  _et al._ (2022); Gonzalez  _et al._
(2022); Nagar  _et al._ (2023); Andrade  _et al._ (2023). A third major family
is the phenomenological inspiral-merger-ringdown, IMRPhenom waveforms Ajith
_et al._ (2007, 2008, 2011); Santamaria  _et al._ (2010); Hannam  _et al._
(2014); Schmidt  _et al._ (2015); Husa  _et al._ (2016); Khan  _et al._
(2016); Dietrich  _et al._ (2017); London  _et al._ (2018); Khan  _et al._
(2019a); Dietrich  _et al._ (2019); Khan  _et al._ (2019b); Pratten  _et al._
(2021, 2020); Thompson  _et al._ (2020); Estellés  _et al._ (2021); García-
Quirós  _et al._ (2020); Estellés  _et al._ (2022a, b); Hamilton  _et al._
(2021); Yu  _et al._ (2023), which have mostly been built in the frequency
domain. This has enabled faster parameter estimation runs compared with the
time domain EOB models. However, there are now also a few time domain
IMRPhenom models Estellés  _et al._ (2022a, b).

Specific waveform models from each family are usually referred to as
approximants. The state of the art approximant from each of the families
listed in the previous paragraph are, respectively, SEOBNRv5PHM Ramos-Buades
_et al._ (2023); Khalil  _et al._ (2023); Mihaylov  _et al._ (2023); SEO
(recent upgrade from SEOBNRv4PHM Ossokine  _et al._ (2020); Gadre  _et al._
(2022); Thomas  _et al._ (2022)), TEOBResumS222Nonprecessing and precessing
approximants in this family are all called by the same name.Akcay  _et al._
(2021); Gamba  _et al._ (2022), IMRPhenomTPHM Estellés  _et al._ (2022b) (time
domain) and IMRPhenomXPHM Pratten  _et al._ (2021) (frequency domain). Note,
we do not employ here the more recent IMRPhenomXODE Yu  _et al._ (2023) or
IMRPhenomXO4a Thompson  _et al._ (2023), nor the upgraded IMRPhenomXPHM
Colleoni  _et al._ (2024) though we do include a brief comparison involving
these two models in App. C. Essentially, these models are “too new” as the
bulk of our work was already completed by the time they appeared.

When a new waveform model is complete, it may undergo a review where it is
extensively compared with numerical relativity simulations and with various
other approximants as well as undergoing parameter estimation tests. Such
studies are also detailed in the articles introducing the specific models
(see, e.g., the works cited above). Recently, the LIGO-Virgo-KAGRA
Collaboration (LVK) also conducted a detailed study on the faithfulness of the
waveforms generated by the models NRSur7dq4, SEOBNRv5PHM, TEOBResumS,
IMRPhenomTPHM, IMRPhenomXPHM and IMRPhenomXO4a using a set of $\gtrsim 1500$
NR simulations, but this is not publicly available. However, some of this work
has been documented in Ref. Ramos-Buades  _et al._ (2023) (also see Ref.
Thompson  _et al._ (2023) for a similar study using BAM waveforms Hamilton
_et al._ (2023); Husa  _et al._ (2008); Brügmann  _et al._ (2008). A separate
study focusing on the parameter estimation performance of the precessing
models was conducted in Ref. Puecher  _et al._ (2023). This work compared the
system parameters inferred by NRSur7dq4, IMRPhenomTPHM, IMRPhenomXPHM and
SEOBNRv4PHM for an ensemble of nearly $60$ “pure” BBH O3 events.

Here, we undertake a more systematic survey based on simulated data. First, in
Sec. IV, using the NR surrogate NRSur7dq4 (henceforth NRSur) as a proxy for
NR, we compute the unfaithfulnesses of $\\{$SEOBNRv5PHM, TEOBResumS,
IMRPhenomTPHM, IMRPhenomXPHM$\\}$ to it for a discretely spaced and a random-
uniformly filled set of intrinsic parameters with the binary mass fixed to
both a light and a heavy value, and inclination fixed to $0$ and $\pi/2$. We
focus solely on the $\ell=2$ waveform strain faithfulness throughout our work.
For the discrete parameter set, we document in detail the deterioration of
faithfulness with increasing mass asymmetry, which is well known, (see e.g.,
in Refs. Estellés  _et al._ (2022b); Hamilton  _et al._ (2021); Ramos-Buades
_et al._ (2023)), but we also uncover multimodalities in the unfaithfulness
distributions coming from strongly precessing cases. We further reveal a
correlation involving the unfaithfulness of the co-precessing $(2,\pm 2)$
multipoles. We additionally contrast the unfaithfulness results from the
discrete set with those of a uniformly-filled parameter set. Using this, we
show that the latter type of parameter set may result in the overestimation of
model faithfulness.

We then move on to direct comparisons with NR waveforms from the SXS catalog
Boyle  _et al._ (2019), choosing simulations that were not used in the
calibration of the models of interest here. We employ a set of 317 SXS
simulations with a mass ratio grouping similar to the discrete set above so
that we can re-apply the same analyses which result in our reaching very
similar conclusions. We further employ a smaller set consisting of 23 SXS
waveforms containing more than 100 GW cycles each in order to briefly assess
the models’ faithfulness for longer inspirals.

Finally, to complete our survey, we perform zero-noise injections of a
moderately and a strongly precessing numerical relativity waveform and recover
their parameters with SEOBNRv5PHM, IMRPhenomTPHM and IMRPhenomXPHM. As we
explain in detail in Sec. VI, we are unable to present PE results for
TEOBResumS for these injections.

Our faithfulness survey provides us with enough data to construct fits for
unfaithfulness over the intrinsic parameter space. These fits can be used to
generate a weighted categorical prior to inform model choice in a joint
Bayesian analysis Hoy (2022). We reduce the dimensionality of the parameter
space by employing the parallel and perpendicular effective spin projections
for our fits. We then construct three dimensional fits to the logarithm of the
unfaithfulness as functions of the mass ratio and the spin projections. These
fits can subsequently be employed to assign weights to each waveform model in
various regions of the parameter space. Our results shown in App. A are
preliminary as we intend to pursue this line of research elsewhere.

The remaining sections are organized as follows. In Sec. II, we briefly review
precession dynamics and precessing waveform construction, and introduce the
parallel and perpendicular spin projections. Sec. III introduces the various
metrics which we employ to assess waveform faithfulness. We summarize our work
in Sec. VII and conclude with a retrospective discussion in Sec. VIII.
Appendices B and C contain additional model comparisons. Finally in App. D, we
present timing benchmarks for the main models that we consider here. We
provide our data and the notebooks/codes/scripts used to obtain it in a git
repository333https://github.com/akcays2/Survey_Precessing_Models.

We work with geometrized units where $G=c=1$. $m_{1}$ and $m_{2}$ denote the
masses of the primary and secondary components of the compact binary system,
with $m_{1}\geq m_{2}$. Accordingly, we define small and large mass ratios as
$q=m_{2}/m_{1}\leq 1$ and $Q=1/q\geq 1$, the symmetric mass ratio
$\eta=q/(1+q)^{2}$, and the total mass $M=m_{1}+m_{2}$. $\mathbf{S}_{1}$ and
$\mathbf{S}_{2}$ denote the spin vectors of each binary component with the
respective dimensionless spin vectors given by
$\bm{\chi}_{i}=\mathbf{S}_{i}/m_{i}^{2}$ and $\chi_{i}:=|\bm{\chi}_{i}|$,
$i=1,2$. Unless otherwise noted, we set $M=1$. Overdots denote time
derivatives whereas hats denote unit vectors. Throughout this article, we use
the terms aligned and parallel interchangeably, as well as perpendicular/in-
plane/planar where the former direction is along the orbital angular momentum
vector and the latter is in the orbital plane.

##  II Review of Precessing Waveform Construction

For this article, we focus only on quasi-spherical binary systems as most
binaries are expected to have circularized by the time they enter the LIGO-
Virgo-KAGRA band Peters and Mathews (1963) which is supported by the
gravitational wave (GW) data so far Romero-Shaw  _et al._ (2019); Nitz  _et
al._ (2019b); Yun  _et al._ (2020); Ramos-Buades  _et al._ (2020a); Abbott
_et al._ (2021e); Lenon  _et al._ (2020); Wu  _et al._ (2020); Iglesias  _et
al._ (2022) (see Refs. Romero-Shaw  _et al._ (2020a, 2022) for a few
exceptions). However, there are binary formation scenarios in which a moderate
amount of eccentricity survives beyond the decihertz regime Breivik  _et al._
(2016); Samsing and Ramirez-Ruiz (2017); Samsing  _et al._ (2018); Rodriguez
_et al._ (2018); Gondán and Kocsis (2019); Banerjee (2018); Liu  _et al._
(2019); Zevin  _et al._ (2021); Michaely and Perets (2020); Hamers  _et al._
(2021); Spera  _et al._ (2022); Samsing  _et al._ (2022); Hellström  _et al._
(2022); Dall’Amico  _et al._ (2023). Therefore, it is important to have
template banks of sufficiently faithful eccentric and precessing waveforms in
the near future.

In general relativity, GWs have only two propagating degrees of freedom, $+$
and $\times$ polarizations, which can be obtained from the following multipole
sum in terms of spin weight $=-2$ harmonics

$h_{+}-ih_{\times}=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}{}_{-2}Y^{\ell
m}(\vartheta_{\text{LN},0},\varphi_{\text{ref}})h_{\ell m}(t),$ (1)

where
$\vartheta_{\text{LN},0}:=\cos^{-1}(\hat{\,\mathbf{L}}_{\text{N},0}\cdot\hat{\mathbf{N}})$
is the orbital inclination with $\hat{\mathbf{N}}$ being the line of sight
vector from the binary’s center of mass to the observer.
$\varphi_{\text{ref}}=\pi/2-\varphi_{c}$ is a constant reference phase Schmidt
_et al._ (2017), where $\varphi_{\text{c}}$ denotes the phase at coalescence.

The polarizations $h_{+,\times}$ couple with the detectors’ antenna patterns
on Earth. For an L-shaped interferometer, the GW strain in the time domain is
given by

$\displaystyle h(t)=$ $\displaystyle
F_{+}(\theta_{s},\phi_{s},\psi_{s})h_{+}(t,\vartheta_{\text{LN},0},\varphi_{\text{ref}})$
$\displaystyle+F_{\times}(\theta_{s},\phi_{s},\psi_{s})h_{\times}(t,\vartheta_{\text{LN},0},\varphi_{\text{ref}}),$
(2)

where $\theta_{s},\phi_{s}$ are the source sky location angles, $\psi_{s}$ is
the source polarization with respect to the detector and $F_{+,\times}$ are
the detector antenna pattern functions that can be found in, e.g., Ref.
Maggiore (2007). The GW strain at the detector (henceforth just the strain)
given by Eq. (2) is the quantity which we use in determining the faithfulness
of the EOB and the phenomenological approximants via a certain sky-maximized
mismatch defined in Sec. IV.

It was shown by Refs. Schmidt  _et al._ (2011, 2012); Boyle  _et al._ (2011)
that precessing waveform multipoles, $h_{\ell m}$, can be built to a good
approximation from the Euler rotation of aligned-spin (AS) multipoles via the
following expression

$h_{\ell
m}=\sum_{m^{\prime}=-\ell}^{\ell}D^{(\ell)\ast}_{m^{\prime},m}(-\gamma,-\beta,-\alpha)\,h^{\text{AS}}_{\ell
m^{\prime}},$ (3)

where $D^{(\ell)}_{m^{\prime},m}$ are Wigner’s D matrices Sakurai (1994);
Brown  _et al._ (2007). $\alpha,\beta$ and $\gamma$ are the Euler angles of
the frame rotation with $\alpha,\beta$ being the spherical angles of the
Newtonian orbital angular momentum vector $\mathbf{L}_{\text{N}}$ with respect
to a chosen frame, e.g.,
$\mathbf{L}_{\text{N},0}:=\mathbf{L}_{\text{N}}(f_{0})$ where $f_{0}$ is a
reference frequency. The third rotation by $\gamma$ uniquely fixes the frame
to the so-called minimal-rotation frame Boyle  _et al._ (2011).

The Euler angles $\alpha,\beta$ are obtained from the time evolution of
$\mathbf{L}_{\text{N}}=\text{L}_{\text{N}}\\!\hat{\,\mathbf{L}}_{\text{N}}$
governed by the PN precession equations. At next-to-leading order (NLO), the
orbit-averaged evolution equations can be written in the following form

$\displaystyle\dot{\mathbf{S}}_{i}$
$\displaystyle=\bm{\Omega}_{i}\times\mathbf{S}{}_{i},$ (4a)
$\displaystyle\dot{\hat{\,\mathbf{L}}}_{\text{N}}$
$\displaystyle=\bm{\Omega}_{\text{NLO}}\times\\!\hat{\,\mathbf{L}}_{\text{N}},$
(4b)

with $i=1,2$ denoting the primary and the secondary, respectively. The
precession frequencies are given by

$\displaystyle\frac{\bm{\Omega}_{i}}{v^{5}}=\eta\left(2+\frac{3}{2}q_{i}\right)\\!\hat{\,\mathbf{L}}_{\text{N}}+\frac{v}{2}\left\\{\mathbf{S}_{j}-3[(q_{i}\mathbf{S}_{i}+\mathbf{S}_{j})\cdot\\!\hat{\,\mathbf{L}}_{\text{N}}]\\!\hat{\,\mathbf{L}}_{\text{N}}\right\\},$
(5a)
$\displaystyle\bm{\Omega}_{\text{NLO}}=-\frac{v}{\eta}\left(\bm{\Omega}_{1}+\bm{\Omega}_{2}\right),$
(5b)

where $v$ is the relative speed of the binary components in the center of mass
frame, $q_{1}=1/q,\,q_{2}=q$, and $j=3-i$. Radiation reaction is incorporated
through the decay of the magnitude of the angular momentum,
$\dot{\,\text{L}}_{\text{N}}$, while the direction of the total angular
momentum vector
$\mathbf{J}=\mathbf{L}_{\text{N}}+\mathbf{S}_{1}+\mathbf{S}_{2}$ is kept fixed
at a desirable PN order. It is standard to rewrite
$\dot{\,\text{L}}_{\text{N}}$ as $\dot{v}(v)$ or $\dot{\omega}(\omega)$ using
$\text{L}_{\text{N}}=\eta/v=\eta\omega^{-1/3}$ where $\omega=v^{3}$ is the
orbital frequency. Further details can be found in, e.g., Refs. Apostolatos
_et al._ (1994); Kidder (1995); Gerosa and Kesden (2016); Akcay  _et al._
(2021); Gerosa  _et al._ (2023); Khalil  _et al._ (2023).

The coupled system of ODEs consisting of the decay of $\text{L}_{\text{N}}$
and Eqs. (4a), (4b) can be numerically solved straightforwardly. Analytic
solutions have also been found at this PN order via the multi scale approach
Kesden  _et al._ (2015); Gerosa  _et al._ (2015); Chatziioannou  _et al._
(2017); Gerosa  _et al._ (2023). PN information exists up to 3.5PN in the
radiation reaction sector Buonanno  _et al._ (2003, 2009); Chatziioannou  _et
al._ (2013) and next5LO in the precession dynamics Khalil  _et al._ (2023).
Once the solution $\mathbf{L}_{\text{N}}(t)$ is known, $\alpha(t)$ and
$\beta(t)$ can be computed immediately, and, subsequently, $\gamma$ from Boyle
_et al._ (2011)

$\dot{\gamma}=\pm\dot{\alpha}\cos\beta,$ (6)

where there is a sign freedom in the right hand side.

Note that Eq. (3) is an approximation for the true precessing multipoles,
$h_{\ell m}^{\text{prec}}$ which can be extracted from NR simulations for
example. One can then obtain the co-precessing multipoles via inverse of the
transformation in Eq. (3)

$h^{\text{coprec}}_{\ell
m}=\sum_{m^{\prime}=-\ell}^{\ell}D^{(\ell)}_{m^{\prime},m}(\alpha,\beta,\gamma)\,h_{\ell
m^{\prime}}.$ (7)

Though $h_{\ell m}^{\text{coprec}}\neq h_{\ell m}^{\text{AS}}$, the AS
multipoles are, in general, a good approximation for the co-precessing
multipoles Ramos-Buades  _et al._ (2020b) as we will show in Sec. IV.1.3.

The AS multipoles satisfy the relation
$h^{\text{AS}}_{\ell,-m}=(-1)^{\ell}h^{\text{AS}\ast}_{\ell m}$, which does
not hold for $h_{\ell m}^{\text{coprec}}$ due to the asymmetric emission of
GWs above and below the orbital plane leading to the “bobbing” of the binary
and eventually resulting in the well-known kick of the final black hole Keppel
_et al._ (2009); Bruegmann  _et al._ (2008); Ramos-Buades  _et al._ (2020b);
Kalaghatgi and Hannam (2021). Although the difference in the $\pm m$
multipoles may be small, it seems to be non-negligible for unbiased parameter
estimation of precessing binaries even at moderate SNRs Kalaghatgi and Hannam
(2021), and most definitely so at high SNRs Kolitsidou  _et al._ (2024). The
only model used in this work that does not neglect this multipole asymmetry is
NRSur7dq4 Varma  _et al._ (2019b), though we note that the newest
phenomenological approximant PhenomXO4a also contains this feature for the
dominant, quadrupolar multipole Ghosh  _et al._ (2023); Thompson  _et al._
(2023).

Another approximation made above is the use of orbit-averaged ODEs to describe
the precession dynamics which removes the nutation of the spins from the
evolution of the binary444For non-orbit-averaged versions, see Eqs. (2.8) and
(2.10) of Ref. Racine (2008).. Though wide nutation angles may leave a unique
signature on future GW events Gerosa  _et al._ (2019), the work of Ref.
Gangardt  _et al._ (2022) has found no strong evidence of nutation in any
detected event thus far.

Returning to the binary with the orbital and spin angular momenta vectors
$\\{\mathbf{L}_{\text{N}},\mathbf{S}_{1},\mathbf{S}_{2}\\}$, let
$\theta_{i}=\cos^{-1}(\\!\hat{\,\mathbf{L}}_{\text{N}}\cdot\hat{\mathbf{S}}_{i})$
denote the polar (tilt) and $\phi_{i}$ the azimuthal angle of each spin vector
with $i=1,2$. The components of the spins parallel and perpendicular to
$\\!\hat{\,\mathbf{L}}_{\text{N}}$ are given by
$\mathbf{S}_{i,||}=(\mathbf{S}_{i}\cdot\\!\hat{\,\mathbf{L}}_{\text{N}})\\!\hat{\,\mathbf{L}}_{\text{N}}=m_{i}^{2}\chi_{i}\cos\theta_{i}\\!\hat{\,\mathbf{L}}_{\text{N}}$
and $\mathbf{S}_{i,\perp}=\mathbf{S}_{i}-\mathbf{S}_{i,||}$ with magnitude
$m_{i}^{2}\chi_{i}\sin\theta_{i}$.

We have already discussed the most commonly employed parallel and
perpendicular scalars, $\chi_{\text{eff}}$ and $\chi_{\text{p}}$, in Sec. I.
$\chi_{\text{eff}}$ is given by Damour (2001)

$\chi_{\text{eff}}=\frac{1}{1+q}(\chi_{1}\cos\theta_{1}+q\chi_{2}\cos\theta_{2})$
(8)

and its perpendicular counterpart by Hannam  _et al._ (2014); Schmidt  _et
al._ (2015)

$\chi_{\text{p}}=\max\left(\chi_{1}\sin\theta_{1},q\frac{4q+3}{4+3q}\chi_{2}\sin\theta_{2}\right).$
(9)

The Kerr spin limit $\chi_{i}\leq 1$ imposes the condition
$0\leq\chi_{\text{p}}\leq 1$ where $\chi_{\text{p}}=0$ corresponds to a non-
precessing (aligned spin or spinless) and $\chi_{\text{p}}=1$ to a maximally
precessing binary. There is an elegant discussion in Ref. Gerosa  _et al._
(2021) on the physical meaning of this definition which leads to a generalized
version of this parameter

$\begin{split}\chi_{\text{p}}^{\text{Gen}}=&\left[(\chi_{1}\sin\theta_{1})^{2}+\left(q\frac{4q+3}{4+3q}\chi_{2}\sin\theta_{2}\right)^{2}\right.\\\
&+\left.2q\frac{4q+3}{4+3q}\chi_{1}\chi_{2}\sin\theta_{1}\sin\theta_{2}\cos(\Delta\phi)\right]^{1/2},\end{split}$
(10)

where

$\cos(\Delta\phi)=({\hat{\mathbf{S}}_{1}\times\\!\hat{\,\mathbf{L}}_{\text{N}}})\cdot({\hat{\mathbf{S}}_{2}\times\\!\hat{\,\mathbf{L}}_{\text{N}}})$
(11)

which is the span of the planar angle between $\mathbf{S}_{1}$ and
$\mathbf{S}_{2}$. $\chi_{\text{p}}^{\text{Gen}}$ is in fact equal to
$|\dot{\hat{\,\mathbf{L}}}_{\text{N}}|/\Omega_{1}$ Gerosa  _et al._ (2021) and
a comparison of Eq. (10) with Eq. (9) reveals that
$\chi_{\text{p}}^{\text{Gen}}$ can exceed 1 ($\chi_{\text{p}}^{\text{Gen}}\leq
2$) which can only be achieved by systems in which both spins are large and
mostly planar Gerosa  _et al._ (2021). This fact has already been used for a
simulated study of O4 events showing that for moderate to high SNRs, double-
spin precession can be inferred De Renzis  _et al._ (2022).

There are several other recently-proposed perpendicular scalars. For example,
Ref. Akcay  _et al._ (2021) introduced the following perpendicular scalar

$\chi_{\perp,\text{J}}=\frac{|\mathbf{S}_{1,\perp}+\mathbf{S}_{2,\perp}|}{M^{2}}.$
(12)

Ref. Hamilton  _et al._ (2021) introduced a modification to $\chi_{\text{p}}$
for improved single-spin mapping in the merger-ringdown regime for binaries
where $Q\lesssim 1.5$. Ref. Thomas  _et al._ (2021) replaced $\chi_{\text{p}}$
with a two-dimensional vector in the phenomenological waveform mapping to
improve the faithfulness of the precessing $(2,\pm 1),(3,\pm 3),(4,\pm 4)$
multipoles. Here, we concern ourselves only with
$\chi_{\text{p}},\chi_{\text{p}}^{\text{Gen}}$ and $\chi_{\perp,\text{J}}$ all
of which we collectively refer to as the $\chi_{\perp}$’s.

As can be seen from the equations defining the $\chi_{\perp}$’s, each one
depends on quantities that evolve in time (in other words, in frequency).
Traditionally, the values for these quantities are usually quoted at a
reference frequency of $20$ Hz, but this choice is somewhat arbitrary.
Alternative reference times (or frequencies) have been suggested such as
$t_{\text{p-100}}:=t_{\text{peak}}-100M$ and
$t_{\text{ISCO}}:=t(f=6^{-3/2}/(M\pi))$, where $t_{\text{peak}}$ denotes the
time at which the co-precessing $(2,2)$ multipole amplitude peaks, and
$t_{\text{ISCO}}$ is the time at which the GW frequency equals twice the
Schwarzschild innermost stable circular orbit (ISCO) frequency. Refs. Johnson-
McDaniel  _et al._ (2022); Mould and Gerosa (2022) have recently adopted yet
another alternative: the $t\to-\infty\,(f\to 0\,\text{Hz})$ limit where the
spin tilt angles $\theta_{1},\theta_{2}$ can be unambigiously obtained.
However, the angle $\Delta\phi$ can not be uniquely determined which means
that we can not compute $\chi_{\text{p}}^{\text{Gen}}$ or
$\chi_{\perp,\text{J}}$ at past infinity. One can also work with averaged
quantities instead. For example, Ref. Gerosa  _et al._ (2021) introduces a
precession cycle averaged version of $\chi_{\text{p}}$ which they further
extend to precession- and RMS-averaged versions of
$\chi_{\text{p}}^{\text{Gen}}$ in Ref. Gerosa  _et al._ (2023). Unless
otherwise noted, we quote the values of the various spin scalars at the
initial time (frequency) $t_{0}(f_{0})$.

##  III Quantifying Waveform Faithfulness

The faithfulness of a given template waveform [strain] $h^{\prime}(t)$ to a
target waveform $h_{0}(t)$ is measured in terms of the waveform match given by

${\cal M}=\max\limits_{t_{\text{c}},\varphi_{\text{c}}}\frac{\langle
h_{0}|h^{\prime}\rangle}{\sqrt{\langle h_{0}|h_{0}\rangle\langle
h^{\prime}|h^{\prime}\rangle}}.$ (13)

The match is maximized over constant time and phase shifts which we take to be
the time and phase shifts at coalescence $t_{\text{c}},\varphi_{\text{c}}$.
The angular brackets denote the noise-weighted inner product

$\hskip 56.9055pt\langle
h_{0}|h^{\prime}\rangle:=4\Re\int_{f_{i}}^{f_{f}}\frac{\tilde{h}_{0}(f)\,\tilde{h}^{\prime\ast}(f)}{S_{\text{n}}(f)}df,$
(14)

where $\tilde{h}(f)$ are the Fourier transforms of the time domain GW strains.
The inner product is weighted by the one-sided power spectral density (PSD)
$S_{\text{n}}(f)$ of the detector noise for which we use the Advanced LIGO
Aasi  _et al._ (2015) zero-detuned high-power design sensitivity aLI (a). The
integration is performed from the initial frequency $f_{i}$ to the final
frequency $f_{f}$ which we set equal to 1024 Hz. As can be seen from Fig. 1,
this choice of final frequency is large enough to encompass the merger-
ringdown portions of the waveforms for BBHs with
$M=37.5M_{\odot},150M_{\odot}$ which will be our chosen values for much of
this article. To mediate the effects caused by finite signal duration, we
choose $f_{i}$ such that there are at least three waveform cycles between it
and $f_{0}$, the initial co-precessing $(2,2)$-multipole frequency. We set the
luminosity distance of the sources to the fiducial value of
$d_{\text{L}}=500\,$Mpc.

Figure 1: Typical waveform amplitudes in the frequency domain from the sample
of BBHs generated for the survey in Sec. IV.1. In particular, we display in
blue (orange) a trio of light (heavy) precessing binaries starting from
$f_{0}=37.5\,$Hz (12 Hz) with $Q=1.\bar{1},2,4$ (dot dashed, dotted, solid) at
a fiducial luminosity distance of 500 Mpc and initial inclination of
$\vartheta_{\text{LN},0}=0$ shown against both Advanced LIGO design
sensitivity aLI (a) (black solid curve) and LIGO Livingston detector strain
sensitivity during O3 aLI (b) (faint gray curve). We used the approximant
IMRPhenomXPHM to generate these waveforms. For each case, we have
$\max(\chi_{\text{p}})=0.8$. The orange and blue vertical bands mark the
regions of $f_{\text{peak}}$ for each triplet of waveforms. Since
$f_{\text{peak}}$ varies for each set of parameters, we opted to show bands
here instead of six additional vertical lines marking individual
$f_{\text{peak}}$.

If the system of interest is precessing then the detector antenna patterns of
Eq. (2) become time-dependent. Consequently, the waveform strain at the
detector becomes dependent on constant and time-varying extrinsic parameters
(and of course the intrinsic parameters) which we can write as

$\begin{split}h(t)&=F_{+}(\theta_{s},\phi_{s},\psi_{s}(t))\;h_{+}(\vartheta_{\text{LN},0},\varphi_{\text{ref}},t_{\text{ref}})\\\
&+F_{\times}(\theta_{s},\phi_{s},\psi_{s}(t))\;h_{\times}(\vartheta_{\text{LN},0},\varphi_{\text{ref}},t_{\text{ref}}).\end{split}$
(15)

By introducing an effective amplitude $\mathcal{A}$

$\mathcal{A}(\theta_{s},\phi_{s})=\sqrt{F^{2}_{+}(\theta_{s},\phi_{s},\psi_{s})+F^{2}_{\times}(\theta_{s},\phi_{s},\psi_{s})}$
(16)

and an effective polarizability $\kappa$ via

$e^{i\kappa(\theta_{s},\phi_{s},\psi_{s})}=\left[F_{+}(\theta_{s},\phi_{s},\psi_{s})+iF_{\times}(\theta_{s},\phi_{s},\psi_{s})\right]/\mathcal{A}(\theta_{s},\phi_{s})\,,$
(17)

we can rewrite Eq. (15) as

$\begin{split}h(t)=\mathcal{A}(\theta_{s},\phi_{s})\\{&\cos\left[\kappa(\theta_{s},\phi_{s},\psi_{s})\right]h_{+}(\vartheta_{\text{LN},0},\varphi_{\text{ref}},t_{\text{ref}})\\\
&+\sin\left[\kappa(\theta_{s},\phi_{s},\psi_{s})\right]h_{\times}(\vartheta_{\text{LN},0},\varphi_{\text{ref}},t_{\text{ref}})\\}\,,\end{split}$
(18)

where we suppressed the $t$ dependence in the right-hand side.

We can now define the sky-maximized (optimized) faithfulness (match) between
the target strain and the waveform template as

${\cal
M}_{\text{opt}}=\max_{t^{\prime}_{\text{ref}},\varphi^{\prime}_{\text{ref}},\kappa^{\prime},\phi^{\prime}}\frac{\langle
h_{0}|h^{\prime}\rangle}{\sqrt{\langle h_{0}|h_{0}\rangle\langle
h^{\prime}|h^{\prime}\rangle}},$ (19)

where $t^{\prime}_{\text{ref}},\varphi^{\prime}_{\text{ref}},\kappa^{\prime}$
are template parameters to be optimized over. The details of the optimization
can be found in Refs. Harry  _et al._ (2016); Gamba  _et al._ (2022). In
particular, the $\kappa^{\prime}$ optimization is performed analytically,
while $t^{\prime}_{\text{ref}}$ is maximized via the inverse fast Fourier
transform. The maximization over $\varphi^{\prime}_{\text{ref}}$ is performed
numerically using a dual annealing algorithm Pratten  _et al._ (2021); Gamba
_et al._ (2022). $\phi^{\prime}$ represents the final degree of freedom to be
maximized over, i.e., the freedom to shift $\phi_{1},\phi_{2}$ by a constant
amount which leaves $\Delta\phi$ unchanged. We perform this maximization using
another dual annealing algorithm Gamba  _et al._ (2022). It is evident from
Eq. (18) that the match is also a function of the inclination angle
$\vartheta_{\text{LN},0}$, but since we fix the inclinations in our
comparisons at a given reference frequency, we do not need to optimize over
them.

Note that ${\cal M}_{\text{opt}}$ is a function of $\varphi_{\text{ref}}$ and
$\kappa$ so, as a final step, we compute its average over an evenly spaced
grid for $\\{\varphi_{\text{ref}},\kappa\\}\in[0,2\pi)\times[0,\pi/2)$ with
$6\times 7=42$ elements to obtain

${\cal M}_{\text{opt,av}}:=\frac{1}{42}\sum_{i=1}^{42}{\cal
M}_{\text{opt}}(\varphi_{\text{ref},i},\kappa_{i}).$ (20)

This is done to marginalize over any dependence of the match on the sky
position and obtain values which depend exclusively on the intrinsic
parameters of the source. This quantity is similar to the sky-and-polarization
averaged faithfulness given by Eq. (35) of Ref. Ramos-Buades  _et al._ (2023),
but we do not average over the inclination.

For the remainder of this article, we employ the sky-averaged, optimized
waveform mismatch

$\bar{{\cal M}}_{\text{opt}}:=1-{\cal M}_{\text{opt,av}}$ (21)

as our faithfulness gauge. The optimized mismatches that we quote henceforth
will always be this average which we may occasionally refer to as the “full
mismatch”. As a check, we also store the minimum and the maximum values of
${\bar{\cal M}}_{\text{opt}}$ over the grid as well as the standard deviation
all which we present in the data provided in our git repository.

To further disentangle possible causes of waveform mismodelling, we
additionally compute mismatches where we truncate the mismatch integral before
the transition to plunge, specifically at $f_{f}={0.6}f_{\text{peak}}$ with
$f_{\text{peak}}:=f\left(t=t_{\text{peak}}\right)$ denoting the GW frequency
at the maximum (peak) amplitude of the co-precessing $\ell=2$ strain of
NRSur7dq4, i.e., the maximum of

$\mathcal{A}_{\text{coprec}}(t):=\sqrt{\sum_{m=-2}^{m=2}h_{2m}^{\text{coprec}}(t)}$
(22)

with $t_{0}<t_{\text{peak}}$ Schmidt  _et al._ (2017). We have checked that
${0.6}f_{\text{peak}}\geq f_{\text{MECO}}$, the minimum energy circular orbit
(MECO) frequency Pratten  _et al._ (2020). Accordingly, we present results for
not only ${\bar{\cal M}}_{\text{opt}}$, but also for the merger-ringdown
excluded (inspiral-only) mismatch

${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}:={\bar{\cal
M}}_{\text{opt}}(f_{f}=0.6f_{\text{peak}}).$ (23)

##  IV Faithfulness Survey I: Comparisons with NRSur7dq4

We begin our survey with Sec. IV.1 where we assess the faithfulness of SEOB,
TEOB, TPHM and XPHM to NRSur7dq4 using a discrete grid for the intrinsic
parameters at zero inclination, by which we mean
$\vartheta_{\text{LN},0}:=\vartheta_{\text{LN}}(f_{0})=0$. Since the direction
of $\\!\hat{\,\mathbf{L}}_{\text{N}}$ evolves in time, so does
$\vartheta_{\text{LN}}$. Thus, the best we can do is specify its value at some
reference frequency. In principle, we could have employed the [approximately]
fixed inclination with respect to the total angular momentum,
$\theta_{\text{JN}}$. However, the $\phi^{\prime}$ optimization of Eq. (19)
becomes non-trivial in this case.

We also include in Sec. IV.1 a comparison in the extrapolation region of
NRSur7dq4, where the model is stated to be robust Varma  _et al._ (2019b). In
Sec. IV.2, we consider the same discrete set at
$\vartheta_{\text{LN},0}=\pi/2$ thus changing the multipole content of the
strain in Eq. (1). We perform a similar analysis in Sec. IV.3 over a random-
uniformly filled intrinsic parameter space and compare our findings with those
of Secs. IV.1 and IV.2. We should add that the most recent version of SEOB
(v5) was released near the completion of this work so we had initially
conducted our survey on SEOBNRv4PHM. As a useful sidenote, we present a brief
comparison of v5 and v4 in App. B.

Total mass | $M\in\\{37.5\;M_{\odot},150\;M_{\odot}\\}$
---|---
Mass ratio | $Q\in\\{1.\bar{1},2,4,6\\}$
Initial frequency | $f_{0}=\begin{cases}\\{30,31,35,38.5\\}\text{Hz}\ \,\text{for}\ M=37.5M_{\odot},\\\ \\{7.5,8,9,10\\}\text{Hz}\qquad\text{for}\ M=150M_{\odot}\\\ \end{cases}$
Spin magnitudes |  $\chi_{1}=\chi_{2}=0.8$
Tilt angles | $\cos\theta_{1,2}=\\{-\tfrac{\sqrt{3}}{2}+k\tfrac{\sqrt{3}}{4}\\}\cup\\{-1\\},k=0,\ldots,4$
| if $\theta_{1}=\pi,\theta_{2}\neq\pi$
Azimuthal angles |  | $\phi_{1}=0,\phi_{2}=k\pi/4,k=1,\ldots,8$
---
Table 1: Relevant parameters for the binary systems which we use in our
faithfulness survey in Sec. IV.1. $f_{0}$ denotes the initial frequency from
which we start the binary evolution. We fix $\phi_{1}=0$ at $f_{0}$ without
loss of generality.  Figure 2: The coverage of the spin space for all our
parameter sets shown using the parallel and perpendicular effective spin
scalars computed at the initial time $t_{0}$. In each panel,
$\chi_{\mathrm{eff}}(t_{0})$ of Eq. (8) is plotted in the vertical axis. From
left to right, $\chi_{\text{p}}(t_{0})$ [Eq.(9)],
$\chi_{\text{p}}^{\text{Gen}}(t_{0})$ [Eq.(10)] and
$\chi_{\perp,\text{J}}(t_{0})$ [Eq.(12)] are plotted in the horizontal axes,
respectively. The blue disks correspond to the values of these quantities
coming from our discrete parameter set of Sec. IV.1 (see Table 1). The green
squares, orange diamonds and the black crosses represent the same quantities
obtained from our random-uniformly filled set (Sec. IV.3), the short and the
long SXS sets (Secs.V.1 and V.2), respectively. As $\chi_{\text{p}}$ has no
dependence on the azimuthal components of the spins, the parameters of the
discrete and the short SXS sets yield degenerate values, clustered at
$\chi_{\text{p}}\approx 0.4,0.72,0.8$, whereas both
$\chi_{\text{p}}^{\text{Gen}}$ and $\chi_{\perp,\text{J}}$ cover their
respective ranges better since they are $\phi_{i}$-dependent.

###  IV.1 Discrete Parameter Set with zero Inclination

As both LIGO and Virgo noise curves are frequency, and therefore mass,
dependent, we consider a set of light and heavy binaries for our survey with
$M=37.5M_{\odot}$ and $150M_{\odot}$. We divide each set into four equal-size
subsets separated by mass ratio values of $Q=1.\bar{1},2,4,6$. We have chosen
$Q=1.\bar{1}\,(q=0.9)$ to break the $Q=1$ symmetry where the nonprecessing
$(2,1)$ multipole equals zero, i.e., $h_{21}^{\text{AS}}=0$. The subset with
$Q=6$ is in the so-called extrapolation region of NRSur7dq4 so we delegate the
comparison pertaining to it to Sec. IV.1.4.

A lighter binary mostly accumulates SNR during its inspiral whereas a heavy
enough binary could have as much SNR accumulated during the merger-ringdown
stage as the inspiral assuming signals enter the detector band at $\gtrsim
20\,$Hz. We illustrate this in Fig. 1 with three light and three heavy BBH
waveforms with $Q=1.\bar{1},2,4$ and the same spin vectors. The merger-
ringdown SNR for the light binaries ranges from a quarter to $\lesssim$ half
of its inspiral counterpart (computed from 38 Hz, see $f_{0}$’s in Table 1)
while the same SNR for the heavy systems ranges from 70% to 90% of its
inspiral counterpart (from 12 Hz). Therefore, by considering light and heavy
systems, we are effectively dividing our survey sample into two halves whereby
merger-ringdown modelling is much more important for waveform faithfulness in
one half (heavy BBHs) than the other (light BBHs).

Any comparison with NRSur7dq4 is ultimately limited by the fact the surrogate
waveforms have a maximum inspiral time length of $4300M$. For this reason, we
have chosen $f_{0}=\\{30,31,35\\}\text{Hz}$ for the $Q=1.\bar{1},2,4$ subsets
of the $M=37.5M_{\odot}$ set and $f_{0}=\\{7.5,8,9\\}\text{Hz}$ for the same
subsets of the $M=150M_{\odot}$ set. These yield between 27 and 47 GW cycles
for light BBHs and between 25 and 45 GW cycles for the heavy BBHs.

Each $Q$ subset consists of a grid of spin angles
$\\{\theta_{1},\theta_{2},\Delta\phi:=\phi_{2}-\phi_{1}\\}$ where the tilt
angles are evenly spaced in $\cos\theta_{i}$ from $-\sqrt{3}/2$ to
$\sqrt{3}/2$ in steps of $\sqrt{3}/4$. We also included the grid points with
$\cos\theta_{i}=-1$ ($i=1$ _or_ 2) with the intention to create a small subset
of cases with near-transitional precession Apostolatos  _et al._ (1994) to
test the models’ robustness. Specifically, the eight cases with
$Q=4,\theta_{1}=\pi,\theta_{2}=\pi/6$ yield
$\\!\hat{\,\mathbf{L}}_{\text{N}}(t_{0})\cdot\mathbf{J}_{\text{N},0}\approx-0.02$
and $|\mathbf{J}_{\text{N},0}|\approx 0.025\approx
0.033|\mathbf{L}_{\text{N},0}|$ (in units of $M^{2}$). The addition of these
extra grid points creates a somewhat skewed coverage of the parameter space as
can be seen by the blue dots in Fig. 2.

$\Delta\phi$ runs from $\pi/4$ to $2\pi$ in steps of $\pi/4$. We fix the
dimensionless Kerr spin parameter $\chi$ to 0.8 for all cases resulting in
$\chi_{\text{p}}\leq
0.8,\chi_{\perp,\text{J}}<0.55,\chi_{\text{p}}^{\text{Gen}}<1.45$ and
$|\chi_{\text{eff}}|<0.7$. We summarize these parameter choices in Table 1,
and show the coverage of the dimensionless intrinsic parameter space in terms
of the parallel and perpendicular spin projections in Fig. 2 as the blue
disks. The figure also shows the same quantities plotted for the uniformly
filled parameter set of Sec. IV.3 and the SXS sets of Sec. V.

Figure 3: The main result of this section: the distributions of the sky-
optimized mismatch, ${\bar{\cal M}}_{\text{opt}}$, given by Eq. (21) between
the models listed in the legend and NRSur7dq4 which we take as a proxy for
numerical relativity. The parameters for the waveforms are given in Table 1
with the $Q=6$ subset excluded from the figure and we consider only the
$\ell=2$ strain for the mismatches. The left (right) panels display the
results for low (high) mass, $M=37.5M_{\odot}$ ($M=150M_{\odot}$), systems,
respectively. SEOB (v5), being the most recent of the four models, shows
improvement, especially for heavier systems. The dashed vertical line in both
panels marks mismatch of $1-0.9^{1/3}\approx 0.035$ translating to an event
loss rate of 10% Owen (1996); Flanagan and Hughes (1998b).

The three-angle $\\{\theta_{1},\theta_{2},\Delta\phi\\}$ grid has in total
$(6^{2}-1)\times 8=280$ elements. This then yields $280\times 3=840$ separate
binaries to consider per total mass separated by $Q=1.\bar{1},2,4$. Overall,
we have $840\times 2=1680$ systems for each of which we compare the waveforms
generated by $\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$
to NRSur via the averaged, optimized mismatch ${\bar{\cal M}}_{\text{opt}}$
given by Eqs. (19, 21) for which the last quantity to consider is the lower
limit of the match integral, i.e., $f_{i}$. For the $M=37.5M_{\odot}$ set, we
pick $f_{i}=f_{0}+3\,\text{Hz}$ which captures approximately 25 to 42 GW
cycles depending on the intrinsic parameters. This allows for plenty of SNR
accumulation in the inspiral stage before transitioning to plunge. On the
other hand, for the $M=150M_{\odot}$ set, we deliberately fix
$f_{i}=11\,\text{Hz}$ leading to fewer inspiral cycles captured in the match
integration, thus increasing the relative contribution of the merger-ringdown
stage to the overall SNR.

In Fig. 3, we show the main result of this section where we plot the
distributions of ${\bar{\cal M}}_{\text{opt}}$ between NRSur7dq4 and SEOB,
TEOB, TPHM, XPHM. In the left (right) panel, we present the mismatches for the
light (heavy), i.e., $M=37.5M_{\odot}\;(150M_{\odot})$, systems for all
$Q=1.\bar{1},2,4$. An apparent feature in the figure is the superior
faithfulness of SEOB exhibiting no cases of ${\bar{\cal
M}}_{\text{opt}}>0.035$ and having $\approx{22}\%$ of the heavy-binary
mismatches below $10^{-3}$. We also note comparable performances among TEOB,
TPHM and XPHM. For the heavy systems, the distributions of ${\bar{\cal
M}}_{\text{opt}}$ for these approximants become wider with more cases of
${\bar{\cal M}}_{\text{opt}}>0.035$.

Perhaps the most curious feature in Fig. 3 is the little “bump” in the heavy
XPHM histogram at ${\bar{\cal M}}_{\text{opt}}\sim 0.3$ with a similar smaller
“island” for the light-binary XPHM mismatches at $0.1\leq{\bar{\cal
M}}_{\text{opt}}\leq 0.3$. For both light and heavy systems, we find that
these are the same 40 cases which have $Q=4,\theta_{1}=\pi$ with the worst
mismatches coming from the $\theta_{2}>\pi/2$ subset. For these
configurations, which have the most negative $z$ projections of spins, i.e.,
$\chi_{\text{eff}}\lesssim-0.5$, the default MSA prescription for the
precession dynamics of XPHM breaks down Chatziioannou  _et al._ (2017) and the
model defaults to the single-spin, next-to-next-to-leading order (NNLO)
prescription for the spin dynamics with a 3PN approximation to $|\mathbf{L}|$.
This is a known shortcoming of XPHM and is expected to be fixed in the next
version of the model Colleoni  _et al._ (2024).

Figure 4:  Merger-ringdown truncated mismatches versus full-length waveform
mismatches. The dark-shaded half violins are the distributions of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ [Eq. (23)] between the numerical relavity
surrogate NRSur7dq4 and the four approximants listed along the bottom
horizontal axis of each panel. Similarly, the light-shaded half violins are
the distributions of ${\bar{\cal M}}_{\text{opt}}$ [Eq. (21)]. which are shown
in the legend as ${\bar{\cal M}}_{\text{opt}}^{X=\text{MR}}$. These are the
same quantities plotted in Fig. 3 as histograms. The left (right) panel
contains the data from the light (heavy), $M=37.5M_{\odot}(150M_{\odot})$,
binaries. We used the default Gaussian kernel density estimator of the seaborn
library Waskom (2021) to smooth the histograms for the violin plots. The
horizontal dashed line marks the mismatch value of $1-0.9^{1/3}\approx 0.035$.

####  IV.1.1 Effects of the Merger-Ringdown Portion of GWs on Faithfulness

The increased amplitude/phase disagreement between NR and waveform models in
the plunge-merger-ringdown stages of the binary evolution is a well-known
shortcoming, usually illustrated in terms of time-domain plots of the
waveforms. For this reason, we investigated how the mismatches plotted in Fig.
3 change when neglecting the merger ringdown portions of the waveforms. We
quantify this in terms of the merger-ringdown truncated (i.e., inspiral only)
mismatch given in Eq. (23). There are two reasons for which we a priori expect
the MR-truncated waveforms to agree better with NRSur, i.e., we expect the
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ distributions to occupy lower
values than their ${\bar{\cal M}}_{\text{opt}}$ counterparts. The first reason
pertains to signal morphology: the inspiral is much smoother than the MR and
therefore easier to optimize over in the match computation. The second reason
has to do with the modelling: the MR stages are more difficult to model
analytically, especially without input from NR simulations. As such, the
waveform models employ NR-informed fits for the MR regime.

Our second expectation is a larger shift between the ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ and ${\bar{\cal M}}_{\text{opt}}$
distributions for the $M=150M_{\odot}$ sample since the MR part of the signal
is much more important for these heavier binaries.

These expectations are confirmed in Fig. 4 where we present the distributions
of ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ and ${\bar{\cal
M}}_{\text{opt}}$ as [half] violin plots for the four approximants for both
the light and the heavy BBH samples. We can clearly see that the distributions
of the MR-excluded mismatches (${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$),
i.e., the darker-shaded half violins, are all shifted to lower values than the
lighter-shaded half violins representing the distributions of ${\bar{\cal
M}}_{\text{opt}}$. This downward shift is more pronounced for the heavier
cases ($M=150M_{\odot}$) with the median of the ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ distributions roughly half an order of
magnitude lower than the ${\bar{\cal M}}_{\text{opt}}$ distributions (less so
for TEOB). The shift is less prominent for the light BBHs ($M=37.5M_{\odot}$)
where the MR portion of the signal is much less important than the inspiral.

One can also discern a smaller secondary peak in most of the ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ distributions which comes from the mismatches
of the $Q=4$ subset. This peak is somewhat obscured in Fig. 4 where we plot
the entire emsemble of $Q=1.\bar{1},2,4$ subsets. We investigate this in more
detail in the next section where we separate the mismatch data by mass ratio.

One can further ask whether or not there are regions in the spin space where
there is a larger gap between ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ and
${\bar{\cal M}}_{\text{opt}}$ than other regions. To this end, we define the
following ratio $\cal{R}:={\bar{\cal M}}_{\text{opt}}/{\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ and look at how it is distributed in the
$\\{\chi_{\perp}\text{'s},\chi_{\mathrm{eff}},Q\\}$ space, where
$\chi_{\perp}\text{'s}$ represents one element of the set
$\\{\chi_{\text{p}},\chi_{\text{p}}^{\text{Gen}},\chi_{\perp,\text{J}}\\}$ at
a time.

For the light BBHs, we observe $\mathcal{R}\lesssim\\{2,3.5,6\\}$ for
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM}\\}$ regardless of $Q$. For
XPHM, $\mathcal{R}\lesssim\\{4.6,6,9\\}$ for $Q=\\{1.\bar{1},2,4\\}$.
Moreover, $\mathcal{R}$ seems to increase monotonically for TEOB and TPHM as
the spin vectors become more planar, i.e., as $\theta_{1,2}\to\pi/2$, whereas
for XPHM, the large ${\mathcal{R}}$ region is mostly situated in the
$\theta_{1,2}\to 5\pi/6$ region. As for SEOB, since
$1\leq{\mathcal{R}}\lesssim 2$, no specific region stands out.

For the heavy BBHs, ${\mathcal{R}}$ increases significantly for all models.
SEOB still yields the smallest values with ${\mathcal{R}}\lesssim 6$
regardless of the mass ratio with a small region of ${\mathcal{R}}\approx 8$
in the $Q=4$ subset. The values of ${\mathcal{R}}$ output by the other three
models show more variation with mass ratio as
${\mathcal{R}}\lesssim\\{12,18,14\\}$ for TEOB,
${\mathcal{R}}\lesssim\\{6,10,12\\}$ for TPHM, and
${\mathcal{R}}\lesssim\\{6,10,17\\}$ for XPHM corresponding to
$Q=\\{1.\bar{1},2,4\\}$. Unlike the light set, we now observe
larger-${\mathcal{R}}$ regions in the $\chi_{\mathrm{eff}}>0$ half of the
parameter space for SEOB, TEOB and TPHM with the opposite relation applying to
XPHM.

Figure 5: Mass ratio separated, merger-ringdown truncated mismatches,
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$, between NRSur7dq4 and the
approximants $\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$
labeled along the bottom axes and plotted in blue, orange, green, red,
respectively. For each approximant, we present four half violins,
representing, from left to right, the $\vartheta_{\text{LN},0}=0$ mismatches
for the subsets with mass ratios $Q=1.\bar{1},2,4,6$ which are shown along the
top axes. The $Q=6$ subset is discussed separately in Sec. IV.1.4 as it
involves comparisons with NRSur7dq4 in its extrapolation region. The
horizontal dashed line marks the mismatch value of $1-0.9^{1/3}\approx 0.035$.

In order to better determine how much of the increased mismatches are
genuinely due to mismodelling of the merger-ringdown regimes, one could
construct hybrid waveforms where the inspiral-only part generated by the four
models is attached to the plunge-merger-ringdown part of the corresponding
NRSur7dq4 waveform. One could then compute ${\bar{\cal M}}_{\text{opt}}$
between the hybrid waveforms and the full NRSur7dq4 waveforms, and compare
these values to the ones obtained here. Though this strategy is, in principle,
straightforward, its implementation is highly non-trivial with precession. We
leave this for future work.

####  IV.1.2 The Dependence of Waveform Unfaithfulness on the Mass Ratio

Another known degradation of waveform faithfulness occurs for systems with
more component mass asymmetry, i.e., large (small) values of $Q\ (q)$. In
terms of our sample, this should translate to increasing mismatches with
increasing $Q$ (or decreasing $q$). In order to better reduce the potential
“contamination” of the mismatches due to merger-ringdown mismodelling, we
present only ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$, separated by $Q$
values here. Accordingly, we should expect to see ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(Q=1.\bar{1})\preceq{\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(Q=2)\preceq{\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(Q=4)\preceq{\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(Q=6)$, where we are appropriating the $\preceq$
symbol to mean that the distributions have higher occupancy per bin at lower
mismatch values. The results of this breakdown by mass ratio are shown in Fig.
5, where we plot the distributions of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ as separate half violins for the
$Q=1.\bar{1},2,4,6$ subsets, but delegate our discussion of the last subset
($Q=6$) to Sec. IV.1.4.

Let us recall that increasing $Q$ increases the number of cycles in the range
$[f_{i},f_{f}]$, so part of the increased mismatch may be simply due to having
longer waveforms. However, this is not the only factor. Given a value of $Q$,
the spin configurations that yield the longest waveforms are those with
$\theta_{i}=\pi/6$ and the shortest waveforms are the ones with
$\theta_{1}=5\pi/6\;(150^{\circ}),\theta_{2}=\pi$. It is true that the cases
with $\\{Q,\theta_{1},\theta_{2}\\}=\\{6,\pi/6,\pi/6\\}$ have higher
mismatches than, e.g., the cases with $\\{1.\bar{1},\pi/6,\pi/6\\}$, where the
former cases have nearly twice as many cycles in the $[f_{i},f_{f}]$ interval.
On the other hand, the mismatches for
$\\{Q,\theta_{1},\theta_{2}\\}=\\{6,\pi/6,\pi/6\\}$ are lower than the cases
with $\\{6,\pi/2\lesssim\theta_{1}\lesssim
5\pi/6,\pi/2\lesssim\theta_{2}\lesssim 5\pi/6\\}$, where the latter cases have
10 to 20 less cycles in band. Therefore, waveform length alone can not explain
the increased mismatches.

Fig. 5 essentially corroborates our expectations. We find that for both light
and heavy systems, the best agreement with NRSur for each approximant is
always by the $Q=1.\bar{1}$ subset with the $Q=2$ subset yielding slightly
worse mismatches. What is additionally apparent from the figure is the
relative upward shift of the mismatch distributions for the $Q=4,6$ subsets,
i.e., a clear deterioration of waveform faithfulness for more mass asymmetric
systems. We also observe an emergence of multimodalities in the distributions
of the $Q=4,6$ mismatches.

In more detail, we observe that the worst mismatches from the $Q=4$ subset
predominantly come from the cases with $\chi_{\text{p}}\gtrsim 0.7$
($0.6\lesssim\chi_{\text{p}}^{\text{Gen}}\leq 1\text{ and
}\chi_{\perp,\text{J}}\gtrsim 0.4$), where we quote the values at $f_{0}$.
This is somewhat expected as these cases precess more strongly. However, for
the heavy-mass subset, some of the worst mismatches show up by
$\chi_{\text{p}}\gtrsim 0.4$ indicating additional mismatch dependence on the
parallel projection of the spins. This might be connected with some of the
observed multimodalities. For example, we discovered that the two distinct
peaks in TEOB’s $Q=4$ mismatch distributions can be mapped to two separate
regions in the parallel-perpendicular spin projection space. Specifically, we
find that for $M=150M_{\odot}$, the mismatches of the upper peak (${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}>10^{-2}$) solely come from cases with
$\chi_{\text{eff}}<0$. Light SEOB exhibits a less prominent bimodality
separated at ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}\approx 0.004$ with the
mismatches coming from the upper peak corresponding to cases with
$\chi_{p}\gtrsim 0.7$ and mostly $\chi_{\text{eff}}<0$. Light TPHM also
exhibits a bimodality separated at ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}\approx 0.004$ corresponding to $\chi_{p}\gtrsim
0.7$ regardless of the value of $\chi_{\mathrm{eff}}$. Heavy XPHM distribution
separates at ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}\approx 0.02$ with
mismatches of the upper peak _all_ coming from $\chi_{\text{eff}}<0$ cases.

Figure 6: The association of certain peaks in the mismatch distributions with
specific regions in the spin parameter space. In each of the four panels, we
scatter-plot the $\\{\chi_{\perp,\text{J}},\chi_{\mathrm{eff}}\\}$ values of
the $Q=4$ subset of the discrete parameter space. The larger disks correspond
to points which yield values of ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$
[Eq. (23)] in the right-most peaks of the histograms plotted in the insets
where the thick, vertical dashed lines separate the multiple modes (peaks) of
each histogram. The smaller squares mark the complementary $Q=4$ cases. The
SEOB, TPHM histograms in the first and third insets correspond to the $Q=4$
half violins from the left panel of Fig. 5, while the rest from the right
panel.

We show how the multimodalities of these four cases map to
$\\{{\chi_{\perp,\text{J}}},\chi_{\text{eff}}\\}$ space in Fig. 6 for the
$Q=4$ subset. For this figure, we opted to use $\chi_{\perp,\text{J}}$ instead
of $\chi_{\text{p}}$ (or $\chi_{\text{p}}^{\text{Gen}}$) for visual clarity.
From the figure, we observe that for $\chi_{\perp,\text{J}}\gtrsim 0.4$
($\chi_{\text{p}}\gtrsim 0.5$) cases, the $\chi_{\text{eff}}<0$ corner of the
parameter space seems to be more challenging than the $\chi_{\text{eff}}>0$
corner for all approximants except for TPHM. The aforementioned breakdown of
XPHM’s MSA prescription is also exhibited in the
$\chi_{\text{eff}}\lesssim-0.5,\chi_{\perp,\text{J}}\lesssim
0.05\;(\chi_{\text{p}}\lesssim 0.2)$ corner of the parameter space.

Finally, let us finish this section with a brief breakdown of the full
mismatches (${\bar{\cal M}}_{\text{opt}}$) in terms of mass ratio so as to
complement Fig. 5. For SEOB, we find no cases of ${\bar{\cal
M}}_{\text{opt}}>0.035$ for $Q\leq 4$ as is clear from Figs. 3 and 4.
$\\{\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ are also robust for $Q\leq
2$ with very few cases of ${\bar{\cal M}}_{\text{opt}}>0.035$ for either light
or heavy BBHs, worst being 3% of the $M=150M_{\odot}\ (37.5M_{\odot})$ cases
for TEOB (TPHM), but their faithfulness degrades for higher values of $Q$. For
example, at $Q=4$, $\\{14\%,40\%,50\%\\}$ of the light-mass cases yield
${\bar{\cal M}}_{\text{opt}}>0.035$ which becomes $\\{40\%,36\%,49\%\\}$ for
the heavy mass set for $\\{\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$,
respectively. We can also contrast these percentages with their ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}>0.035$ counterparts which are
$\\{0\%,2\%,4\%\\}$ for the light and $\\{0\%,1\%,2\%\\}$ for the heavy sets.
In short, there is a noticeable degradation of model performance in going from
$Q\lesssim 2$ to $Q\gtrsim 4$ for all models including SEOB though overall it
is always better than 0.965 faithful for $Q\leq 4$.

####  IV.1.3 The Effect of AS/Co-precessing multipoles on Faithfulness

Since each precessing multipole of the four models is constructed using Eq.
(3) or some variant of it (e.g., using co-precessing multipoles instead of AS
multipoles), there are two main modelling systematics at interplay here:  (i)
systematics coming from the modelling of precession dynamics that are manifest
in the Euler angles used in the frame rotation,  (ii) systematics in the AS
(or co-precessing) multipoles that are being Euler rotated.  We focus on the
latter in this section.

Even if the Euler angles of
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ equal those
of NRSur (not actually the case), strain mismatches can still be large if the
Euler-rotated multipoles poorly match NRSur’s co-precessing multipoles.
Moreover, in the case of unequal Euler angles, the AS/co-precessing multipole
mismatch may still end up being the dominant systematic. We investigate this
here by computing AS/co-precessing multipole mismatches between the four
models and NRSur’s co-precessing multipoles.

Figure 7: The correlation between ${\bar{\cal M}}_{22,\text{AS}}$ and
${\bar{\cal M}}(\vartheta_{\text{LN},0}=0)$ for the $M=37.5M_{\odot}$ cases,
where ${\bar{\cal M}}_{22,\text{AS}}$ is the mismatch [via Eq. (13)] between
NRSur’s co-precessing $(2,2)$ multipole and each AS $(2,2)$ multipole of the
four models labelled in the legends where we also show the square of the
Pearson correlation coefficient. The upper panels show the mismatches from the
combined $Q=1.\bar{1}$ and $Q=2$ subsets, while the lower panels show the
results from the $Q=4$ subsets. The cluster of red dots in the upper left
corner of the $Q=4$ XPHM panel are the cases for which the aforementioned MSA
prescription of XPHM’s precession dynamics breaks down.

First, let us add a few details. To our knowledge, among the four models, TEOB
is the only one that rotates actual AS multipoles. SEOB twists AS multipoles
with the constant spin parameters $\chi_{i}$ replaced by the time varying
$\bm{\chi}_{i}(t)\cdot\\!\hat{\,\mathbf{L}}_{\text{N}}(t)$ evolved via SEOB-
specific dynamics. TPHM and XPHM rotate AS multipoles with modified remnant
properties Estellés  _et al._ (2022b); Pratten  _et al._ (2021). Though these
multipoles may be referred to as co-precessing, they are all approximations
and thus, not exactly equal to $h^{\text{coprec}}_{\ell m}$ of Eq. (7).
Therefore, let us denote these approximate co-precessing multipoles by
${h}^{\text{AS}}_{\ell m}$ as well. These still satisfy the
$m\leftrightarrow-m$ multipole symmetry. Though exact for non-precessing
systems, this is an approximation for the true co-precessing multipoles which
are known to violate this symmetry (see, e.g., Fig. 2 of Ref. Ramos-Buades
_et al._ (2020b)).

We can now ask: how faithful are AS multipoles to the true co-precessing
multipoles? And what is the penalty in using AS multipoles to construct
precessing waveforms? Ref. Ramos-Buades  _et al._ (2020b) provides the first
detailed answer to this question using comparisons to 72 NR simulations though
only six cases have $Q\geq 4$. They find that the AS $(2,\pm 2)$ multipoles
are faithful representations of their co-precessing counterparts with only one
(five) out of the 72 simulations resulting in mismatches larger than 0.03
(0.01) (see their Fig. 2 and Table III), but that the AS $(2,\pm 1)$
multipoles may not be considered to be so (ibid.). Especially relevant here is
their specific comparison with a short BAM Husa  _et al._ (2008); Brügmann
_et al._ (2008) simulation (ID 28) with $Q=3,\chi_{\text{p}}=0.8$, where they
identify the $m$-multipole asymmetry as the cause of the high mismatch.

In the following, we conduct a study similar to Ref. Ramos-Buades  _et al._
(2020b)’s. Specifically, we compute mismatches between the co-precessing
$(2,\pm 2)$ multipoles of NRSur7dq4 and those of
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$. We do this
by modifying the source code for the gwsurrogate package. We focus on
$h_{22}^{\text{AS}}$ which is dominant in Eq. (3), which we also treat as a
proxy for results pertaining to $h_{2,-2}^{\text{AS}}$. We proceed by first
computing a new quantity: ${\bar{\cal M}}_{22,\text{AS}}$, the standard
$\\{t_{\text{c}},\varphi_{\text{c}}\\}$-optimized mismatch between NRSur’s co-
precessing $(2,2)$ multipole and each AS $(2,2)$ multipole of the four models.
We then look for correlations between these mismatches and the strain
mismatches ${\bar{\cal M}}(\vartheta_{\text{LN},0}=0),{\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$. Indeed, as Fig. 7 highlights,
there is a very prominent correlation between $\log{\bar{\cal
M}}_{22,\text{AS}}$ and $\log{\bar{\cal M}}$ for TPHM and TEOB, with XPHM
showing a weaker trend and SEOB hardly showing any. The linear correlation of
the logs implies a powerlaw relation between these mismatches. Additionally,
we observe that

  1. (i)

the correlations are stronger between ${\bar{\cal M}}_{22,\text{AS}}$ and
${\bar{\cal M}}$ than ${\bar{\cal M}}_{\text{opt}}$. This means that the
optimization of the mismatch over extrinsic parameters partially “smears out”
this relation;

  2. (ii)

the correlations are much more prominent for the light BBHs, where more
inspiral cycles are included in the match integrals;

  3. (iii)

the correlations persist in all mass ratio subsets, but are stronger for
$Q=1.\bar{1}$ and 2.

For TPHM and TEOB, the $Q=\\{1.\bar{1},2\\}$ subset correlation coefficients
are larger than 0.9. Moreover, we observe that for the cases that yield
${\bar{\cal M}}_{22,\text{AS}}\gtrsim 0.05$, both TEOB’s and TPHM’s
${\bar{\cal M}}$ can be fit by lines parametrized by ${\bar{\cal
M}}_{22,\text{AS}}$ with slopes $\gtrsim 1$, not just their logarithms. This
is an indication that the AS multipole unfaithfulness, when large enough,
becomes the dominant systematic in waveform faithfulness for these two models.
This linearity is less obvious for SEOB and XPHM. The former model produces
the most faithful AS multipoles: ${\bar{\cal M}}_{22,\text{AS}}<0.01$ for
almost every case, whereas the latter model has the most outliers indicating
other systematics contaminating this relation such as the already mentioned
MSA-related breakdown. Despite these exceptions, when ${\bar{\cal
M}}_{22,\text{AS}}$ is large enough, the observations above also hold for
these two models.

This systematic becomes more severe with increasing mass ratio as we
determined by comparing the slopes of the linear fits to the mismatch data for
TEOB and TPHM between the $Q=1.\bar{1},2$, $Q=4$ and the $Q=6$ subsets. In
each comparison, the slope of the larger $Q$ subset was greater than the slope
of the smaller $Q$ subset, though we observed this to be more severe for TEOB
than TPHM.

When we looked at the intrinsic parameters of the cases that yield the highest
values for ${\bar{\cal M}}_{22,\text{AS}}$ for each approximant, we found that
the worst mismatches come from the $65^{\circ}\lesssim\theta_{1,2}\lesssim
115^{\circ}$ region, resulting in the largest planar spin projections. We have
already identified this region as the most challenging one consistent with the
literature. The fact that the AS multipole mismatches are also the highest in
this region tells us that the true co-precessing multipoles differ much more
than the AS multipoles for cases with mostly planar spins. We believe this to
be mostly due to the imposed $m\leftrightarrow-m$ symmetry on the AS
multipoles, which we confirmed to be more violated for the strongly precessing
cases. We did this by comparing ${\bar{\cal M}}_{22,\text{AS}}$ with
${\bar{\cal M}}^{\text{sym}}_{22,\text{AS}}$, i.e., the AS multipole
mismatches computed with respect to the $m$-symmetrized co-precesing NRSur7dq4
multipoles, and found that the quantity ${\bar{\cal
M}}_{22,\text{AS}}-{\bar{\cal M}}^{\text{sym}}_{22,\text{AS}}$ peaks for
$65^{\circ}\lesssim\theta_{1,2}\lesssim 115^{\circ}$, i.e., for the cases with
mostly in-plane spins.

When we extended the above analysis of the AS $(2,2)$ multipole to the AS
$(2,1)$ multipole, we encountered some unexpectedly high values of ${\bar{\cal
M}}_{21,\text{AS}}$ for $Q=1.\bar{1}$ for all four models. This contradicts
the thus-far observed trend of increasing mismatches with increasing $Q$. Upon
closer examination, we discovered that the problem was due to NRSur7dq4: since
for $Q=1$, $h_{21}^{\text{AS}}=h_{21}^{\text{coprec}}=0$ exactly, the model
generates a co-precessing (2,1) multipole that is barely above numerical error
for $Q=1.\bar{1}$. For this reason, we choose to disregard all $Q=1.\bar{1}$
values of ${\bar{\cal M}}_{21,\text{AS}}$.

Focusing on the $Q\geq 2$ part of the discrete set, we observe that the
previously observed tight correlation for the AS (2,2) multipole becomes much
less strong for the AS (2,1) mode. In fact, the linear trend between
$\log{\bar{\cal M}}_{21,\text{AS}}$ and $\log{\bar{\cal M}}$ is only
discernible for TPHM, with SEOB and TEOB showing some positive correlation at
the high mismatch end as well. As was the case with the AS (2,2) multipole,
the trends are more visible for the $M=37.5M_{\odot}$ BBH set. As for the
cases that yield the worst AS (2,1) mismatches, we find more of a spread in
the $\\{\theta_{1},\theta_{2}\\}$ space, but with the worst common mismatches
coming from the $\theta_{1}\leq\pi/2,\theta_{2}\geq\pi/2$ corner, while the
cases with maximum or near maximum $|\chi_{\mathrm{eff}}|$ values consistently
yield the lowest mismatches. We should keep in mind that even for the worst
cases where this correlation is tight, it need not necessarily affect the
overall unfaithfulness ${\bar{\cal M}}_{\text{opt}}$ significantly as the AS
(2,1) multipole amplitude is roughly an order of magnitude smaller than its
(2,2) counterpart. And as we discuss further below, the AS $(2,\pm 1)$
multipoles have a rather weak effect on the overall strain faithfulness.

The AS $(2,\pm 1)$ multipoles’ importance for faithful precessing waveform
construction was investigated by Ref. Ramos-Buades  _et al._ (2020b) where
they showed that while the inclusion of the AS $(2,\pm 1)$ multipoles degrades
the precessing $(2,\pm 1)$ multipole faithfulness (see their Table III and
Fig. 4), it marginally improves the precessing $(2,\pm 2)$ multipoles. How
these changes affect the faithfulness of the waveform strain is illustrated
for a $Q=5$ case in their work, where it is shown that the inclusion of the AS
$(2,\pm 1)$ multipoles lowers the $\ell=2$ strain unfaithfulness by 0.01 in
the highest mismatch regions as can be gathered by comparing the middle left
with the middle right panel in their Fig. 5. In the words of Ref. Ramos-Buades
_et al._ (2020b), this indicates that the improvement in the precessing
$(2,\pm 2)$ multipoles, due to the inclusion of the $(2,\pm 1)$ AS multipoles,
compensates for the degradation of the precessing $(2,\pm 1)$ multipoles.

We extend their work here by turning “off” $h_{2\pm 1}^{\text{coprec}}$ for
NRSur7dq4 and looking at the unfaithfulness of the resulting strain with
respect to the full-multipole content NRSur7dq4 strain. Specifically, we
compute ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ and ${\bar{\cal
M}}_{\text{opt}}$ for the $\ell=2$ strain at inclinations of zero and $\pi/2$.
Though we use the same multipole content as Ref. Ramos-Buades  _et al._
(2020b), they considered all inclinations in the one example that they present
whereas we suffice with two inclinations here, but repeat the computation for
all 1120 cases of the discrete set.

For $\vartheta_{\text{LN},0}=0$, we find that ${\bar{\cal M}}_{\text{opt}}$
ranges from $\sim~{}10^{-8}$ to $10^{-2}$ with mismatches increasing with $Q$
and the highest values coming from $\theta_{1,2}=\pi/2$ cases of the $Q=4,6$
subsets. For $\vartheta_{\text{LN},0}=\pi/2$, we observe
$10^{-5}\lesssim{\bar{\cal M}}_{\text{opt}}<0.035$, with the dependence on $Q$
becoming less important. In other words, when the co-precessing $(2,\pm 1)$
multipoles are omitted, we observe increased strain unfaithfulness at higher
inclinations. Accordingly, the co-precessing $(2,\pm 1)$ multipoles improve
the strain faithfulness more for higher inclinations, consistent with the
findings of Ref. Ramos-Buades  _et al._ (2020b). Interestingly, the largest
values of ${\bar{\cal M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$ come
from cases with $\theta_{1,2}\geq 5\pi/6$, not $\pi/2$ as was the case for
$\vartheta_{\text{LN},0}=0$.

####  IV.1.4 Mismatches in the Extrapolation Region of NRSur7dq4

Since the approximants exhibit a degradation in faithfulness as $Q$ increases,
we extended our investigation into the extrapolation regime of NRSur7dq4 by
setting $Q=6$ for the same discrete grid as before, resulting in two
additional subsets of light and heavy BBHs with 280 cases each. We have opted
to present this comparison separately from the $Q\leq 4$ cases of the previous
subsections since, as Ref. Varma  _et al._ (2019b) advises, we must exercise
caution when using NRSur7dq4 in its extrapolation region, i.e.,
$4<Q<6,0.8<\chi_{i}<1$. Nonetheless, comparisons with 100 SXS simulations with
$Q=6$ and $\chi_{i}\leq 0.8$ have yielded mismatches ranging from $\sim
10^{-4}$ to $\sim 10^{-2}$ with the median at $\sim 10^{-3}$. Overall, this is
roughly an order of magnitude worse than NRSur7dq4’s performance in its
training region, $Q\leq 4,\chi_{i}\leq 0.8$, but good enough for our purposes
here. A more direct approach would be to compare the four models with the
aforementioned 100 SXS simulations, but we were unable to retrieve these from
the SXS database. We present comparisons using different sets of SXS
simulations in Sec. V.

Returning again to Fig. 5, where the $Q=6$ ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ distributions are plotted as the rightmost
half violins for each color, we see that the worsening trends observed for the
($Q$ = 4)-case mismatches persist for the $Q=6$ cases. As before, the mismatch
distributions exhibit bi or trimodalities with the larger-mismatch peaks of
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{XPHM}\\}$’s distributions coming
from cases with either large planar spins, i.e., $\chi_{\text{p}}\gtrsim 0.7$
or with moderate $\chi_{\text{p}}$ values combined with negative
$\chi_{\text{eff}}$ values. TPHM once again exhibits mismatches that are more
symmetric with respect to $\chi_{\mathrm{eff}}$. Overall, the $Q=6$
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ values for the four models are
slightly worse than the $Q=4$ values. Interestingly, the
$\chi_{\text{eff}}\lesssim-0.5,\chi_{\text{p}}\lesssim 0.2$ corner of the
parameter space that previously yielded high mismatches for XPHM has shrunk to
a single point that has the most negative $\chi_{\mathrm{eff}}$ value of the
entire $Q=6$ set. This minor improvement in XPHM’s performance is partly due
to the fact that we no longer have the near-transitional precession cases
encountered for $Q=4$ where
$|\mathbf{J}_{\text{N},0}|\ll|\mathbf{L}_{\text{N},0}|$. which were causing
the MSA prescription to breakdown forcing the model to default to the NNLO
prescription for precession dynamics.

Taking ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ as a gauge of inspiral-only
model faithfulness, we find that
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ are
reasonably faithful at this mass ratio with the percentage of cases of
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}>0.035$ being
$\\{0\%,0\%,5\%,28\%\\}$ for the light and $\\{0\%,14\%,2\%,11\%\\}$ for the
heavy-mass sets. As for the full mismatch ${\bar{\cal M}}_{\text{opt}}$, the
percentage of cases greater than $0.035$ is, respectively,
$\\{0\%,34\%,37\%,75\%\\}$ for the light and $\\{16\%,41\%,46\%,77\%\\}$ for
the heavy-mass set.

###  IV.2 Behavior of Model Unfaithfulness under changing Inclination

Our aim in this section is to investigate how the inclusion of the precesssing
$(2,1)$ multipole, i.e., $h_{21}$, affects the faithfulness. Eq. (1) tells us
that the contribution of $h_{21}$ is maximal for $\vartheta_{\text{LN}}=\pi/2$
yielding $|{}_{-2}Y^{21}|/|{}_{-2}Y^{22}|=2$. Though
$|h_{21}^{\text{AS}}|<|h_{22}^{\text{AS}}|$ always holds true, for precessing
multipoles, one can observe that $|h_{21}|\sim|h_{22}|$ (see, e.g., Fig. 1 of
Ref. Varma  _et al._ (2019b)). This is due to the power from the AS $(2,\pm
2)$ multipoles rotating into $h_{21}$ which is time (frequency) dependent, as
is $\vartheta_{\text{LN}}$. Accordingly, the ${}_{-2}Y^{21}h_{21}$ term in the
strain (1) can be dominant. Therefore, comparisons of strain mismatches for
cases with large inclinations offer an indirect way to assess the faithfulness
of “subdominant” multipoles such as $h_{21}$ which, as we just explained, may
not be subdominant at all for certain configurations. The analysis almost
equally applies to $h_{2,-1}$ due to the precessing multipoles approximately
inheriting the $m\leftrightarrow-m$ symmetry of the AS multipoles. For a more
detailed investigation of the mismodelling of $h_{21}$ and higher multipoles
$(\ell>2)$, see Ref. Ramos-Buades  _et al._ (2020b).

This analysis also indirectly gauges the performance of the Euler angles used
in the frame rotations. Since, as shown, the faithfulness of the co-precessing
$(2,\pm 2)$ multipoles matters much more than that of the co-precessing
$(2,\pm 1)$ multipoles, the cases for which the former are very faithful, but
the resulting ${\bar{\cal M}}_{\text{opt}}$ values are high may tell us
something about the mismodeling of the Euler angles.

Figure 8: Comparison of the mismatch distributions between face-on
($\vartheta_{\text{LN},0}=0$) and edge-on ($\vartheta_{\text{LN},0}=\pi/2$)
configurations. In the top row, we display violin plots of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ vs. ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ for the
$M=37.5M_{\odot}$ ($M=150M_{\odot}$) sets in the left (right) panels.
Likewise, the bottom row shows ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$ vs. ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$ for the same sets. Note that
we exclude the $Q=6$ subset here, but we include it in Fig. 9. Figure 9:
Comparison of the mismatch distributions between “face-on”
($\vartheta_{\text{LN},0}=0$) and “edge-on” ($\vartheta_{\text{LN},0}=\pi/2$)
configurations separated by mass ratio. In particular, we present results here
only for the MR-truncated, sky-angle optimized mismatch ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$. The top row displays results for the light
($M=37.5M_{\odot}$) BBH set, and the bottom row the heavy ($M=150M_{\odot}$)
BBHs.

We proceed by re-employing the discrete parameter set of Sec. IV.1 with
$\vartheta_{\text{LN},0}=\pi/2$, and compute the mismatches ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}},{\bar{\cal M}}_{\text{opt}}$ as before. Our
results are presented in Fig. 8 where we show the distributions of the
mismatches between the $\vartheta_{\text{LN},0}=0$ and
$\vartheta_{\text{LN},0}=\pi/2$ sets for the $Q\leq 4$ portion for both light
and heavy binaries. Several features stand out in the figure which we
highlight next.

First, as in the $\vartheta_{\text{LN},0}=0$ case, the distributions of
${\bar{\cal M}}_{\text{opt}}$ are in general broader than ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ distributions and shifted to higher values.
This shift is larger for the heavy masses. Second, the inspiral-only
faithfulness of all approximants remains mostly unchanged between
$\vartheta_{\text{LN},0}=0$ and $\vartheta_{\text{LN},0}=\pi/2$, with TEOB
even showing a slight improvement for $\vartheta_{\text{LN},0}=\pi/2$. There
is also an overall trend of the high-mismatch tails being reduced.

In Fig. 9, we show a decomposition of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ and ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ in terms of
mass ratio, where the upper (lower) panels correspond to the $M=37.5M_{\odot}\
(M=150M_{\odot})$ set. We note several observations. First, the
$\vartheta_{\text{LN},0}=\pi/2$ distributions exhibit less or no modality
compared with the $\vartheta_{\text{LN},0}=0$ distributions. Second, $Q\geq 4$
$\vartheta_{\text{LN},0}=\pi/2$ mismatches have narrower distributions and
mostly occupy lower values than their $\vartheta_{\text{LN},0}=0$
counterparts. The converse holds true for the $Q\leq 2$ mismatches with the
exception of TEOB. Third, the medians of the $\vartheta_{\text{LN},0}=\pi/2$
distributions shift to lower values than the medians of the
$\vartheta_{\text{LN},0}=0$ distributions as $Q$ increases, especially for the
EOB models. The especially remarkable feature is the improved inspiral-only
faithfulness for all models for $Q=4$ at $\vartheta_{\text{LN},0}=\pi/2$ as
compared with $\vartheta_{\text{LN},0}=0$.

For $Q=6$, there is mostly a reduction of the bi/trimodalities which we
consider to be is an improvement. However, the situation is somewhat more
complicated for $Q=6$ as the values of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ of the left-most
peaks are mostly lower than the smallest values of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ whereas the
values of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ of the right-most
peaks are mostly higher than the largest values of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$. So, the
question of whether or not the models show improved inspiral-only faithfulness
for $Q=6$ may have a subjective, rather than an objective answer. Therefore,
let us focus on the $Q=4$ cases next.

Given the above observations, an obvious question is: what is the cause of the
improved faithfulness for $Q=4$ in going from $\vartheta_{\text{LN},0}=0$ to
$\vartheta_{\text{LN},0}=\pi/2$? We hypothesize that this may be due to the
$h_{2,\pm 1}$ multipoles improving the strain when they are more faithful than
$h_{2,\pm 2}$. Given that at $\vartheta_{\text{LN},0}=\pi/2$, we have
$|{}_{-2}Y^{2,\pm 1}|/|{}_{-2}Y^{2,\pm 2}|=2$, this improvement more likely
matters only for cases where $|h_{2,\pm 1}|\gtrsim|h_{2,\pm 2}|$. This may
then result in ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ being less than
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$.
Henceforth, we will simply mention the (2,1) and (2,2) multipoles with the
same results applying to the $(2,-1),(2,-2)$ multipoles.

We test our hypothesis on the TEOB mismatches first since it shows the most
$Q=4$ improvement for $\vartheta_{\text{LN},0}=\pi/2$ as well as showing
improvements for $Q=1.\bar{1},2$. We begin by computing two new quantities:
${\bar{\cal M}}_{22}^{\text{noMR}},{\bar{\cal M}}_{21}^{\text{noMR}}$, i.e.,
the inspiral-only $\\{t_{c},\varphi_{c}\\}$-maximized mismatches [see Eq.
(13)] of the precessing $(2,2)$ and $(2,1)$ multipoles between NRSur7dq4 and
TEOB. We also compute the ratio

$R:=\frac{2}{N}\sum_{i=0}^{N}\frac{|h_{21,i}|}{|h_{22,i}|},$ (24)

where we have $N+1$ samples in the time domain labelled by $i$. $R$ is the
Riemann sum of the following integral

$\int\text{d}t\frac{\left|{}_{-2}Y^{21}(\pi/2,0)\,h_{21}(t)\right|}{\left|{}_{-2}Y^{22}(\pi/2,0)\,h_{22}(t)\right|}$
(25)

$R>1$ tells us that the $(2,\pm 1)$ terms have larger magnitudes than their
$(2,\pm 2)$ counterparts in the waveform strain. Thus when $R>1$ and
${\bar{\cal M}}_{21}^{\text{noMR}}<{\bar{\cal M}}_{22}^{\text{noMR}}$, we
expect the more faithful $h_{21}$ to result in improved faithfulness for
$\vartheta_{\text{LN},0}=\pi/2$ as compared to $\vartheta_{\text{LN},0}=0$,
i.e., ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)\preceq{\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ in the respective
distributions.

We proceed by counting the cases satisfying the condition:
$[R>1]\land[{\bar{\cal M}}_{21}^{\text{noMR}}<{\bar{\cal
M}}_{22}^{\text{noMR}}]$ and comparing the resulting distributions for
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$
with those of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$. The expectation is
that we should see the distributions of the affirmative cases satisfying
${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)\preceq{\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ and that there are
enough such cases to affect the overall distributions as seen in Fig. 9. We
find our hypothesis validated by the $M=37.5M_{\odot},Q=4$ subsets of SEOB and
TEOB as well as the $M=150M_{\odot},Q=4$ subset of TPHM. We chose these
subsets because it is evident from Fig. 9 that this improvement is most
prominent for them. We also confirmed our hypothesis for the entire
$M=37.5M_{\odot},Q\leq 4$ TEOB set. We did not use XPHM because it is non-
trivial to obtain its precessing multipoles in the time domain as it is a
frequency domain approximant.

The MR-included version of Fig. 9 (which we do not include) shows less
modality in the ${\bar{\cal M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$
distributions. However, there are fewer $Q$ subsets where
$\vartheta_{\text{LN},0}=\pi/2$ mismatches are lower than the
$\vartheta_{\text{LN},0}=0$ mismatches. The exceptions to this relation are
SEOB’s light and heavy $Q=4,6$ mismatches. For these, we verified that the
above arguments held.

The fixed angle $\vartheta_{\text{JN}}$ is also at play here: precession-
induced modulations to the waveform decrease as $\vartheta_{\text{JN}}\to 0$
Apostolatos  _et al._ (1994). Using the LALSimulation function
SimInspiralTransformPrecessingWvf2PE, we computed $\vartheta_{\text{JN}}$ and
observe that as the $Q=4$ subset’s spin tilt angles go from pointing
“northward” to southward,
$\left.\vartheta_{\text{JN}}\right|_{\vartheta_{\text{LN},0}=\pi/2}$ drops
below $\left.\vartheta_{\text{JN}}\right|_{\vartheta_{\text{LN},0}=0}$. In
other words, the higher inclination cases yield smaller precession-induced
modulations, which ties in with the discussion above, but does not explain why
the precessing $(2,\pm 1)$ multipoles are sometimes more faithful than $(2,\pm
2)$ when both have the same co-precessing (AS) multipole content dominated by
$h_{2,\pm 2}^{\text{coprec}}$ ($h_{2,\pm 2}^{\text{AS}}$).

To answer this question, one would have to dismantle the individual co-
precessing (AS) multipole contribution to the twist expression (3) and compare
the interplay between various multipoles. We do not pursue this matter any
further here. Ref. Ramos-Buades  _et al._ (2020b) present additional analysis
in this direction, where they compare waveform strain and individual
precessing multipole mismatches for cases in which all [useful] AS multipole
content has been included and for cases where some AS multipoles, e.g.,
$(2,\pm 1)$ have been excluded.

The final observation relevant to this section is that just as in the
$\vartheta_{\text{LN},0}=0$ case, the distributions for ${\bar{\cal
M}}_{\text{opt}}(\pi/2)$ are shifted upward, i.e., to worse values, as
compared to ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ distributions,
especially for $Q=4,6$. The fact that ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}\leq{\bar{\cal M}}_{\text{opt}}$ is not
surprising since the merger-ringdown part of the signal is more complex as we
already mentioned. We can nonetheless ask whether or not there might be
regions in the spin space where there is a larger gap between ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ and ${\bar{\cal M}}_{\text{opt}}$ than other
regions. To this end, we again looked at how the ratio $\cal{R}={\bar{\cal
M}}_{\text{opt}}/{\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ is distributed in
the spin space. The approximants show a lot variability with respect to each
other and across $Q$ values. However, one nearly common feature is this ratio
reaching $\mathcal{O}(10)$ for $\theta_{1,2}=\pi/2$ with XPHM as the
exception. Additionally, SEOB has the smallest maxima for the ratio while XPHM
has the largest. XPHM has more occurences for large ratios in the
$\chi_{\mathrm{eff}}<0$ half whereas the other three models yield more large
ratios in the $\chi_{\mathrm{eff}}>0$ half. XPHM also behaves differently for
the $Q=4$ subset with the ratio showing very little variation in the spin
space and tending to the smallest values at $\theta_{1,2}=\pi/2$.

###  IV.3 Comparisons over a uniformly filled Parameter Space

We now consider a uniformly populated parameter space consisting of 1000 BBHs
with $M\in[{35},{225}]M_{\odot},Q\in[1,4],\chi_{i}\in[0.1,0.8]$. We choose to
make the spin tilt angles uniform in their cosine, i.e.,
$\cos\theta_{i}\in[-1,1]$ with $\phi_{1}=0,\phi_{2}\in[0,2\pi)$. To ensure
once again that the resulting inspirals are not longer than 4300M, we set
$f_{0}$ equal to a linearly decreasing function of $M$ from 35 Hz to $5.75$ Hz
as $M$ increases from 35$M_{\odot}$ to ${225}\,M_{\odot}$. We fix
$f_{i}=f_{0}+{3}$ Hz ($11$ Hz) for $M\lesssim 210M_{\odot}\,(\gtrsim
210M_{\odot})$ in the integral (14). In the following, we present the
mismatches similarly to what we have done for the discrete set of Sec. IV.1 so
that we can illustrate the differences in the approximant performance between
a more traditional parameter space (uniformly filled) and one for which
certain parameters are tailored to potentially better expose precession
related modelling systematics.

We summarize our main results in Fig. 10, where in the top row we compare
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ with ${\bar{\cal M}}_{\text{opt}}$
at inclinations of $\vartheta_{\text{LN},0}=0$ ($\pi/2$) in the left (right)
panels. As expected, the relation ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}\preceq{\bar{\cal M}}_{\text{opt}}$ holds likely
due to the combination of MR modelling systematics and simpler signal
morphology for the inspiral. The shifts between the peaks of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ and ${\bar{\cal M}}_{\text{opt}}$
distributions are larger for $\vartheta_{\text{LN},0}=\pi/2$ than
$\vartheta_{\text{LN},0}=0$ as can be discerned from the upper left and right
panels of the figure. This was also the case for the discrete set though this
is somewhat hard to extract from Fig. 8. This is a reaffirmation of the fact
that the modelling of the plunge-merger-ringdown stages of the precessing
$(2,\pm 1)$ multipoles remains challenging.

In the bottom panels of Fig. 10, we show how the mismatches shift upward more
for SEOB and XPHM than for TEOB and TPHM as the inclination goes from 0 to
$\pi/2$, more consistent with the $M=150M_{\odot}$(right) panels of Fig. 8
than $M=37.5M_{\odot}$ (left) ones. This somewhat makes sense as
$150M_{\odot}$ roughly equals the median of $M\in[35,225]M_{\odot}$, whereas
$37.5M_{\odot}$ is at the lower end.

Figure 10: Mismatches between NRSur7dq4 and SEOB, TEOB, TPHM, XPHM for the
uniformly filled parameter space of 1000 cases. Upper panels compare
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ with ${\bar{\cal M}}_{\text{opt}}$
at inclinations of $\vartheta_{\text{LN},0}=0$ (left) and $\pi/2$ (right) and
the lower panels compare ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ with ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ (left) and
${\bar{\cal M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$ with ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$ (right).

Fig. 11 better highlights the differences in model faithfulness resulting from
using a uniformly filled set vs. our purpose-built discrete set of Secs. IV.1,
IV.2. In this figure, the left half violins (darker shaded) are the mismatch
distributions from the entire
$M=\\{37.5M_{\odot},150M_{\odot}\\},Q=\\{1.\bar{1},2,4\\}$ discrete set, while
the lighter-shaded right violins correspond to the mismatch distributions from
the 1000-case random-uniformly filled set. One immediately notices that the
distributions for the latter (uniform) set have smaller minima, means and
maxima than the former (discrete) set. Though the shifts between the
respective minima are not so large, the offset between the means and the
maxima are quite prominent. The distributions for the uniform set are also
smoother, narrower and unimodal owing to the nature of the uniformly-filled
intrinsic parameter space.

Figure 11: Comparison of model mismatches to NRSur7dq4 between the discrete
parameter set and the random-uniformly filled one. The left half violins are
the mismatch distributions from the entire
$M=\\{37.5M_{\odot},150M_{\odot}\\},Q=\\{1.\bar{1},2,4\\}$ discrete set while
the right ones are the mismatch distributions from the 1000-case random-
uniformly filled set. In the top (bottom) row, we plot ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ (${\bar{\cal M}}_{\text{opt}}$) at
inclinations of 0 (left panels) and $\pi/2$ (right panels).

The main conclusion of Fig. 11 is that the model faithfulness is better for
the uniform set than our purpose-built discrete set. This is, of course, not
surprising since intrinsic parameters of the discrete set are such that it
contains many more cases with strong precession than the uniform set, where
$\chi_{\text{p}}(t_{0})\geq 0.6$ for 60% of the cases for the former vs. 13%
for the latter. We already documented this contrasting coverage of the spin
space in Fig. 2. Additionally, one third of the discrete set is made up of
$Q=4$ cases whereas the uniform set has $\sim 17\%$ of its cases with $3.5\leq
Q<4$. Recall that these high-$Q$ cases appear to be the most challenging for
the models that we consider here, as can be gathered from Figs. 5 and 9.

Perhaps the difference between the two sets is best illustrated by considering
the “extremal corner” of the five-dimensional
$\\{Q,\chi_{1},\chi_{2},\theta_{1},\theta_{2}\\}$ parameter space. Let us as,
an example, assign to this region the following bounds: $Q\in[3.5,4],$
$\chi_{i}\in[0.7,0.8],$ and $\cos(\theta_{i})\in[-0.1,0.1]$. A straightforward
computation shows that the uniform set needs to consist of 20379 points in
order to have a $50\%$ chance that one case will be inside this extremal
region, whereas the discrete parameter set has 8 out of 840 cases in this
extremal corner (see Table 1).

Let us elaborate why the differences shown in Fig. 11 matter. Our concern is
for strongly precessing events like GW200129_065458 for which the PE routines
should mostly draw samples from the large $\chi_{\text{p}}$ region of the spin
space. This in itself is not a problem. However, if one employs a model whose
faithfulness, to a target model such as NRSur7dq4, has been documented over a
uniformly filled parameter space, then one may overestimate the faitfulness of
that model. For example, a figure akin to our Fig. 10 may lead to the
conclusion that model “X” is better than 0.035 faithful in the entire
parameter space, but it may turn out that this model X actually has a large
mismatch tail extending to $\sim 0.1$ in the relevant corner of the parameter
space as illustrated in our Fig. 11. This would be quite problematic for
reliable parameter estimation from strongly precessing events. In short,
uniformly sampled intrinsic parameter spaces may not contain the best set of
parameters to properly “stress-test” precessing approximants.

Let us conclude this section by providing additional metrics. As done in the
previous section, we once again look at the ratio $\cal{R}={\bar{\cal
M}}_{\text{opt}}/{\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$. For
$\vartheta_{\text{LN},0}=0$, we observe that $\mathcal{R}\lesssim 5$
everywhere in the parameter space for all four models. For
$\vartheta_{\text{LN},0}=\pi/2$, we record that $\mathcal{R}\lesssim 10$. For
both inclinations, these upper bounds are roughly half of the corresponding
bounds of $\mathcal{R}$ for the discrete parameter set. As before XPHM has the
largest values of $\cal{R}$ and SEOB the smallest with TPHM yielding very
similar magnitudes for $\cal{R}$ in the parameter space. One major difference
with respect to the discrete set is that as we uniformly fill the mass ratio
space from $Q=1$ to 4, the clearcut regions of large $\cal{R}$ observed for
the discrete set at, e.g., $Q=4$, are no longer there. Overall, the values of
$\cal{R}$ are more uniformly spread in the spin space and show less variation.

Finally, we count the percentage of cases with ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}>0.035$ and ${\bar{\cal M}}_{\text{opt}}>0.035$
for both inclinations. As is clear from Fig. 10, this is 0% for all SEOB
mismatches and 0% for all other ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$
values regardless of inclination. As for ${\bar{\cal M}}_{\text{opt}}$
percentages, we record $\\{2\%,0.4\%,0.8\%\\}$ for ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)>0.035$, and
$\\{0.1\%,0\%,0.2\%\\}$ for ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)>0.035$, respectively for
$\\{\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$. These can be compared with
the ones from the discrete set provided at the end of Sec. IV.1.2.

##  V Faithfulness Survey II: Comparisons with Numerical Relativity Waveforms

In the second part of our survey, we assess the faithfulness of
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ to numerical
relativity waveforms from the SXS sxs  catalog. Specifically, we compare
against 317 short SXS waveforms and 23 long SXS waveforms of the LVCNR catalog
Schmidt  _et al._ (2017), where short (long) means that the total number of GW
cycles in the simulations is $\lesssim 70\,(\gtrsim 125)$. As such, the short
SXS waveforms are more comparable in length to our discrete and uniform set
NRSur7dq4 waveforms. For consistency with the previous sections, we set
$M=37.5M_{\odot},150M_{\odot}$ once again for the short waveform comparisons.
For the long waveforms, we consider several other values for $M$ as there are
only 23 cases to compare.

###  V.1 Short SXS Waveforms

Figure 12: Mismatches defined by Eq. (21) between the four approximants and
the 317 short SXS waveforms shown at an inclination of
$\vartheta_{\text{LN},0}=0$ for a direct comparison with Fig. 3. The left and
right columns of panels correspond to mismatches from the light
($M=37.5M_{\odot}$) and heavy ($M=150M_{\odot}$) binaries.

Figure 13: Merger-ringdown truncated and full mismatches between the four
approximants and the 317 short SXS waveforms shown for inclinations of 0 and
$\pi/2$. The left (right) $2\times 2$ blocks of figures correspond to light,
$M=37.5M_{\odot}$, (heavy, $M=150M_{\odot}$) binaries. The top panels compare
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ with ${\bar{\cal M}}_{\text{opt}}$
for $\vartheta_{\text{LN},0}=0$ and $\pi/2$ and the bottom panels compare
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ with
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ and
${\bar{\cal M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$ with ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$.

We employ a specific data set containing 317 short SXS waveforms, with
$\chi_{\text{p}}\lesssim{0.89}$, $\chi_{\mathrm{eff}}\in[{-0.6,0.66}]$ and
$Q\lesssim 6$. These have between 13 and 70 GW cycles. Though there are
roughly 1500 SXS BBH simulations from which we can extract waveforms, we
purposefully formed a specific subset similar to that of the discrete set of
Secs. IV.1, IV.2. Specifically, our short SXS set has 98 cases with $Q<1.3$,
114 cases with $1.9<Q<2.1$ and 101 cases with $3.99<Q<4.01$, thus having
nearly the same mass ratio split as our discrete set. As there are only three
$Q\approx 5$ and one $Q\approx 6$ simulations in our set, we discuss them
separately than the 313 $Q\lesssim 4$ cases. The relevant spin space
parameters of these simulations are plotted in Fig. 2, where it can be seen
that there are many cases with $\chi_{\text{p}}\approx 0.8$ and
$|\chi_{\mathrm{eff}}|\ll 1$. We obtained the parameters for each simulation
using LALSimulation LIGO Scientific Collaboration  _et al._ (2018)555We employ
the tools at https://git.ligo.org/waveforms/lvcnrpy/-/tree/master with public
SXS data which requires the recasting of SXS data into the LVCNR catalog
format for which we employed the tools of the sxs package Boyle and Scheel
(2023). We provide the list of the 317 SXS simulations that we use in the
linked git repository.

For the match computation, we retain the same values of
$d_{\text{L}}=500\,$Mpc, $f_{f}=1024\,$Hz and the same grid over the
$\\{\kappa,\phi_{\text{ref}}\\}$ space. As SXS simulations start from
differing frequencies, we do not employ a fixed value for $f_{0}$, but use
instead the f_lower_at_1MSUN attribute Schmidt  _et al._ (2017) of each
simulation’s data to obtain the frequency in Hertz for $1M_{\odot}$ then
rescale for $M=37.5M_{\odot},150M_{\odot}$. As before, we set
$f_{i}=f_{0}+3\,$Hz ($11\,$Hz) in the match integral for light (heavy)
systems.

The main results for zero inclination ($\vartheta_{\text{LN},0}=0$) are shown
in Fig. 12 where we plot the distributions of ${\bar{\cal M}}_{\text{opt}}$.
What stands out the most is the superior performance of SEOB as compared with
the other three models. It yields only a handful of cases of ${\bar{\cal
M}}_{\text{opt}}>0.035$ whereas TEOB, TPHM and XPHM have many mismatches above
0.035, especially for the light mass BBHs, with a few cases of even
${\bar{\cal M}}_{\text{opt}}>0.1$, which have
$\\{Q,\chi_{\text{p}}\\}=\\{1.\bar{1},0.89\\},\\{2,0.85\\},\\{4,0.8\\},\\{6,0.77\\}$.
A quick count yields that $\\{0.3,30,35,27\\}\%$ $(\\{1,31,17,9\\}\%)$ of the
317 cases have ${\bar{\cal M}}_{\text{opt}}>0.035$ for
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ for the
light (heavy) binaries.

####  V.1.1 Effects of Merger-Ringdown on Faithfulness

As we did in Sec. IV.1.1, we compare ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ with ${\bar{\cal M}}_{\text{opt}}$ for all
four models in order to investigate how the model faithfulness is affected by
the inclusion of the MR regime. When considering the inspiral-only mismatch
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$, the percentages above drop
significantly, to $\\{0.3,4,4,5\\}\%$ $(\\{0.6,1,0.6,1\\}\%)$. This drop is
shown across the top row of Fig. 13 where we display the distributions for
both ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ and ${\bar{\cal
M}}_{\text{opt}}$ at $\vartheta_{\text{LN},0}=0$ and $\pi/2$ for
$M=37.5M_{\odot},150M_{\odot}$. As was the case with our discrete set, we
observe an increased shift between the ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ and ${\bar{\cal M}}_{\text{opt}}$
distributions as the binary mass increases from $37.5M_{\odot}$ to
$150M_{\odot}$. The upward shift also increases more when the inclination goes
from 0 to $\pi/2$ though less clearly so for the heavy binaries. We also
observe that all four models behave extremely similarly for _inspiral_ -only
heavy BBHs, for which we show the distributions of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ in the upper left
panel of Fig. 14. This is not the case for the longer-inspiral light-mass
binaries except for SEOB for which we observe very similar ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}$ distributions at both masses. We discuss our
results for ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ further below.

Looking, as before, at the spread of the ratio $\cal{R}={\bar{\cal
M}}_{\text{opt}}/{\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ over spin space,
we see that $\cal{R}\lesssim\text{4}$ for light SEOB, TPHM and XPHM, with
TEOB’s range up to 7. For heavy BBHs, SEOB, XPHM yield
$\cal{R}\lesssim\text{7}$, with $\cal{R}\lesssim\text{14}$ for TPHM and
$\cal{R}\lesssim\text{22}$ for TEOB. The near-maximal and maximal values
cluster around the $\theta_{1,2}\approx\pi/2$ region of the parameter space.
Moreover, the $\theta_{1}>\pi/2$ half of the parameter space yields lower
values for $\cal{R}$ than the $\theta_{1}<\pi/2$ half though with plenty of
variability in the $\\{\theta_{1},\theta_{2}\\}$ space. The range of $\cal{R}$
shows a weaker dependence on the mass ratio here than was observed in Sec.
IV.1.1. This is due to the fact that we had previously used the same spin
parameters at all three mass ratios, whereas in this section, the spin
parameters and the mass ratios are dictated by the SXS simulations.

The mismatch distributions also exhibit some multimodality regardless of
histogram binning, which is more prominent for TEOB, TPHM and XPHM for
$M=37.5M_{\odot}$ as evident in Fig. 12. This was not the case for the $Q\leq
4$ discrete set as can be gathered from Fig. 3 though multimodalities emerged
for its $Q=4$ subset shown in Figs. 5, 9. More multimodalities are seen in
Fig. 13 in the various distributions such as TEOB and TPHM’s ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$ for $M=150M_{\odot}$ (lower
right panel). Some of the modes seem to be due entirely to the MR part of the
mismatch. For example, neither the $M=37.5M_{\odot}$, nor the $M=150M_{\odot}$
distributions for SEOB’s ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ show any of the
secondary modes present in the ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$ distributions as seen in Fig. 12,
which stand out more clearly using different histogram binning. We observe
similar, MR-related modes appearing in all the distributions for ${\bar{\cal
M}}_{\text{opt}}$ for all models. Therefore, we narrow our focus to
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ in the next few paragraphs.

Figure 14: Inspiral-only mismatches at $\vartheta_{\text{LN},0}=0$, i.e.,
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ between
our entire $Q\lesssim 4$ short SXS set of waveforms and SEOB, TEOB, TPHM, XPHM
plotted as a function of $\chi_{\perp,\text{J}}$ [Eq. (12)]. Since
$\chi_{\perp,\text{J}}$ is maximized for $\theta_{1,2}\to\pi/2$, we observe
increasing mismatches as the spins become more planar, i.e., for stronger
precession.

####  V.1.2 Dependence of Inspiral-only Faithfulness on Intrinsic Parameters

Figure 15: Sky-optimized mismatches [Eq. (21)] between
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ and the 317
short SXS waveforms shown for inclinations of $\vartheta_{\text{LN},0}=0$ and
$\pi/2$, separated into three subsets in terms of the mass ratio $Q$. The top
(bottom) panels show the mismatch distributions for binaries with total mass
$M=37.5M_{\odot}\,(150M_{\odot})$.

As done in Sec. IV.1, we investigate whether or not specific modes in the
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$ distributions can be attributed
specific regions in the $\\{Q,\theta_{1},\theta_{2}\\}$ space. As briefly
mentioned, we observe unimodal, positively skewed distributions for the heavy-
mass ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$
for all models as well as for the light-mass ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ SEOB distribution.
However, the same distributions for light-mass
$\\{\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ exhibit a clear secondary
mode peaking between ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}=0.02$ and
0.03, which emerges at $Q\approx 2$, and becomes dominant at $Q\approx 4$.
This behavior is consistent for
$\\{\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ indicating very similar
inspiral performance for these models.

The tilt angles corresponding to the $Q\approx 2$ secondary mode all come from
cases with $\theta_{1,2}\approx\pi/2,\chi_{1,2}\approx 0.8$ similar to our
results in Sec. IV.1.2. For $Q\approx 4$, it is mostly the cases with
$\theta_{1}\approx\pi/2$ and $\chi_{1}\approx 0.8$ that yield the mismatches
in the now-dominant mode of the ${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}$
distributions for $\\{\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$. There is
no strong $\theta_{2}$ dependence as $|\mathbf{S}_{2}|\ll|\mathbf{S}_{1}|$ for
$Q\geq 4$. There is however a second smaller region of high mismatches
clustered around $\theta_{1}\approx 3\pi/4$. The common theme in the above
analysis is that strongly precessing cases remain challenging for the models.
This is even true for comparable mass ratios, where we found that the highest
${\bar{\cal M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ values
come from cases with largest planar spin projections. We demonstrate this in
Fig. 14, where we plot ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=0)$ vs.
$\chi_{\perp,\text{J}}$ for the entire $M=37.5M_{\odot},Q\lesssim 4$ short SXS
set. We opt for $\chi_{\perp,\text{J}}$ as $\chi_{\text{p}}$ yields many
degenerate values and $\chi_{\text{p}}^{\text{Gen}}$ correlates less clearly.
Recall that $\chi_{\perp,\text{J}}$ becomes maximum for
$\theta_{1,2}\to\pi/2$, thus the figure tells us that the highest mismatches
come from cases with $\theta_{1,2}\approx\pi/2$.

####  V.1.3 Model Faithfulness under increasing Inclination

Returning to Fig. 13, we first focus on the bottom row where we compare
$\vartheta_{\text{LN},0}=0$ mismatches with the
$\vartheta_{\text{LN},0}=\pi/2$ ones. For light BBHs, we observe once again
that the distributions for $\vartheta_{\text{LN},0}=\pi/2$ are narrower,
consistent with the left panels of Fig. 8. Overall the inspiral-only
performance of the models is comparable modulo SEOB’s superior performance for
light BBHs. TEOB and TPHM yield very similar ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ distributions.
For light BBHs, XPHM’s distribution is also very similar to these two. For
heavy BBHs, there is less similarity in the distributions. Nonetheless, the
ranges of the mismatches for all models are very comparable modulo a secondary
mode appearing for XPHM seen in the upper right panel of Fig. 8. A quick count
yields $\\{0.3,0.6,0.6,1\\}\%$ of the light mass, and $\\{0,0,0,0\\}\%$ of the
heavy mass cases resulting in ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)>0.035$. These
percentages are lower than their $\vartheta_{\text{LN},0}=0$ counterparts. A
similar count returns $\\{0.3,28,23,20\\}\%\;(\\{0,22,4,8\\}\%)$ for the
percentage of light (heavy) cases yielding ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)>0.035$, which are also lower
than the corresponding $\vartheta_{\text{LN},0}=0$ percentages.

The large offsets between the ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2)$ and ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$ distributions for both light
and heavy masses shown in Fig. 13 is consistent with our previous results of
Sec. IV.2. As before, we also observe that changing inclination has a
subdominant effect on model faithfulness compared to the MR portion of the
waveforms. As for the ratio $\cal{R}$, we observe that
$\cal{R}\lesssim\\{\text{3, 8, 5, 5}\\}$ for the light and
$\lesssim\\{\text{5,10,9,14}\\}$ for the heavy
$\\{{\texttt{SEOB}},\texttt{TEOB},\texttt{TPHM},\texttt{XPHM}\\}$ mismatches.
Comparing these with the $\vartheta_{\text{LN},0}=0$ ratios of Sec. V.1.1
tells us that $\cal{R}$’s range has barely changed for the light binaries. In
the case of the heavy binaries, $\cal{R}$’s maximal values have decreased for
SEOB, TEOB and TPHM, while it has increased for XPHM. The large changes in the
maximal values of $\cal{R}$ are due to a few cases, each of which might be
challenging to a particular model, while not to others, and vice versa.

The $\theta_{1,2}\approx\pi/2$ region still results in larger values of
$\cal{R}$. And the $\theta_{1}>\pi/2$ half of the parameter space still yields
lower values for $\cal{R}$ than the $\theta_{1}<\pi/2$ half though this is
less distinct than it was for $\vartheta_{\text{LN},0}=0$, and with XPHM
violating this trend with the region $\theta_{1,2}\gtrsim 2\pi/3$ also
returning large values of $\cal{R}$. Consistently with our previous findings,
the maximal values of $\cal{R}$ do not vary much with changing mass ratio.

A striking feature in Fig. 13 is the prominent bimodalities in the
$\vartheta_{\text{LN},0}=\pi/2$ mismatch distributions for $M=150M_{\odot}$,
which are not present for $M=37.5M_{\odot}$. This bimodality is most prominent
for the heavy-mass TEOB and TPHM mismatches, and is seen, albeit less
strongly, for the heavy SEOB and the light TEOB, TPHM mismatches (referring
specifically to ${\bar{\cal M}}_{\text{opt}}$). Decomposing these mismatches
further into their mass ratio subsets as shown in Fig. 15 clarifies that the
bimodalities are due to the mismatches clustering into two separate regions
for $Q\lesssim 2$ vs. $Q\approx 4$. We had already observed a similar behavior
for the $\vartheta_{\text{LN},0}=0$ mismatches as $Q$ increased from 1 to 4,
which we replot in Fig. 15, mirroring our $\vartheta_{\text{LN},0}=\pi/2$
results. The cases causing these large mismatches have
$\theta_{1}\in[1.45,1.6]$ and $\chi_{\text{p}}\approx 0.8$, which, for
$Q\approx 4$, result in values of ${\bar{\cal
M}}_{\text{opt}}^{\text{noMR}}(\vartheta_{\text{LN},0}=\pi/2$) much larger
than the ones coming from the $Q\lesssim 2$ cases. These inspiral-only
mismatches further degrade by up to an order of magnitude when we consider the
full mismatches ${\bar{\cal M}}_{\text{opt}}$. Note that the values of
$\theta_{1},\theta_{2}$ do not alone explain the increased mismatches for
$Q\approx 4$ as the $Q\approx 2$ subset has even more cases with
$\theta_{1,2}\approx\pi/2$. The number of cycles do not provide a satisfactory
explanation either as our chosen $Q\approx 4$ simulations are comparable in
length to the $Q\approx 2$ subset.

Another interesting feature of Fig. 15 is the fact that for $Q\approx 4$ the
$\vartheta_{\text{LN},0}=\pi/2$ mismatches are mostly lower than
$\vartheta_{\text{LN},0}=0$ ones. Specifically, for all light BBHs and the
heavy SEOB and TPHM, we clearly have ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)\preceq{\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$. Focusing on the $Q\approx 4$
SEOB subsets alone, we confirmed that the improved faithfulness at
$\vartheta_{\text{LN},0}=\pi/2$ is due to the arguments made in Sec. IV.2,
i.e., for cases with $|h_{21}|\gtrsim|h_{22}|$ and ${\bar{\cal
M}}_{21}\lesssim{\bar{\cal M}}_{22}$, we see ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)\preceq{\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$ for the distributions. In short,
as we had observed in Sec. IV.2, when the $h_{2,\pm 1}$ multipoles are more NR
faithful than $h_{2,\pm 2}$ and their amplitudes are comparable to $|h_{2,\pm
2}|$, the resulting high-inclination strain matches are mostly better than
zero-inclination matches.

####  V.1.4 Effects of Co-precessing Multipoles on Faithfulness

Figure 16: The correlation between ${\bar{\cal M}}_{22,\text{AS}}$ and
${\bar{\cal M}}(\vartheta_{\text{LN},0}=0)$ for the $M=37.5M_{\odot},Q\lesssim
2$ cases of our short SXS set. ${\bar{\cal M}}_{22,\text{AS}}$ is the mismatch
[Eq. (13)] between the SXS co-precessing $(2,2)$ multipole and each AS $(2,2)$
multipole of the four models. Each legend also gives the $R^{2}$ value for the
$\log$-$\log$ data sets. This figure is akin to the top panel of Fig. 7.

As we had done in Sec. IV.1.3, we investigate once again whether or not the
unfaithfulness of the AS (co-precessing) multipoles have any effect on model
performance. As before, we find a correlation between $\log({\bar{\cal
M}}_{22,\text{AS}})$ and $\log({\bar{\cal M}}(\vartheta_{\text{LN},0}=0))$,
where ${\bar{\cal M}}_{22,\text{AS}}$ is the mismatch between the SXS co-
precessing $(2,2)$ multipole666 We employ the Scri package Boyle (2013); Boyle
_et al._ (2014); Boyle (2016); Boyle  _et al._ (2020) to extract the SXS co-
precessing multipoles. and each model’s AS $(2,2)$ multipole. The correlation
for our short SXS set is stronger than the discrete set of Sec. IV.1 and
persists whether or not we consider $\log{\bar{\cal M}}$ or $\log{\bar{\cal
M}}_{\text{opt}}$ or inspiral-only vs. full mismatches. It is also seen for
the $M=150M_{\odot}$ waveforms, which for the SXS set are in general longer
than the discrete set, which were limited by NRSur7dq4’s length restrictions.
As before, the correlation is strongest for the $M=37.5M_{\odot}$ TEOB and
TPHM mismatches, but also persists for SEOB and XPHM as can be seen in Fig.
16, where we plot only the data from the $Q\lesssim 2$ subset. Thus, this
figure can be directly compared with the top row of Fig. 7. In Fig. 16, we
also provide the $R^{2}$ values for the $\log$-$\log$ scatter plots, which for
TEOB and TPHM are near unity, and equals almost 0.9 for XPHM.

Following Sec. IV.1.3, we additionally compute the slopes of the linear fits
to ${\bar{\cal M}}(\vartheta_{\text{LN},0}=0)$ as a function of ${\bar{\cal
M}}_{22,\text{AS}}$ for TEOB, TPHM and XPHM. We find that the entire
$Q\lesssim 4$ set of ${\bar{\cal M}}(\vartheta_{\text{LN},0}=0)$ (313 out 317
cases) for TEOB and TPHM can be fit by lines with slopes of approximately 0.9,
and 0.5 for XPHM. This is consistent with our findings in Sec. IV.1.3 that
once the AS $(2,\pm 2)$ multipole unfaithfulness exceeds a certain value, it
becomes the dominant source of unfaithfulness for the strain mismatch.

Since the $(2,\pm 2)$ AS/co-precessing multipoles also make up the dominant
contribution to the precessing $(2,\pm 1)$ multipoles, we expect the
correlation to persist for ${\bar{\cal M}}(\vartheta_{\text{LN},0}=\pi/2)$,
albeit less strongly. Indeed, the $R^{2}$ values provided in Fig. 16 drop to
$\\{0.2,0.77,0.85,0.6\\}$, respectively for four models. And the slopes of the
linear fits are comparable. Finally, as we have already seen in Sec. IV.1.3
the AS $(2,\pm 1)$ multipoles do not yield a similar correlation: their
magnitudes are too small. As done in Sec. IV.1.3, we also checked how much
these modes matter for strain faithfulness. Specifically, we re-constructed
the SXS strains without the $h_{2,\pm 1}^{\text{coprec}}$ content and computed
the mismatch to the unaltered SXS strains. As before, we found that the
mismatch is at most 0.015 (0.07) at
$\vartheta_{\text{LN},0}=0\;(\vartheta_{\text{LN},0}=\pi/2)$. At both
inclinations, it is only few cases that yield such mismatches with the
majority of the values below 0.001 (0.01).

| ${\bar{\cal M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$ |  ${\bar{\cal M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)$
---|---|---
SXS ID | Model | $37.5M_{\odot}$ | $150M_{\odot}$ | $37.5M_{\odot}$ |  $150M_{\odot}$
0057 | SEOB | 0.027 | 0.032 | 0.011 | 0.0076
TEOB | 0.024 | 0.038 | 0.015 | 0.010
TPHM | 0.029 | 0.032 | 0.012 | 0.0081
XPHM | 0.036 | 0.038 | 0.015 | 0.0098
0058 | SEOB | 0.0069 | 0.0041 | 0.0055 | 0.0050
TEOB | 0.022 | 0.014 | 0.018 | 0.019
TPHM | 0.026 | 0.0087 | 0.013 | 0.0086
XPHM | 0.012 | 0.010 | 0.013 | 0.018
0062 | SEOB | 0.014 | 0.019 | 0.0039 | 0.0026
TEOB | 0.020 | 0.038 | 0.0079 | 0.0076
TPHM | 0.016 | 0.023 | 0.0051 | 0.0045
XPHM | 0.025 | 0.047 | 0.0062 | 0.0077
0165 | SEOB | 0.061 | 0.055 | 0.017 | 0.016
TEOB | 0.088 | 0.076 | 0.028 | 0.023
TPHM | 0.063 | 0.11 | 0.028 | 0.027
XPHM | 0.13 | 0.093 | 0.080 | 0.057
Table 2: The sky-optimized mismatches between the four $Q>4$ SXS waveforms and
the four models that we consider in this work.

####  V.1.5 $Q\geq 5$ Comparisons

As already noted, we only have four cases with $Q\geq 5$ in our short SXS set.
Three of these, namely SXS:0057,0058,0062, have $Q\approx
5,-0.18\lesssim\chi_{\mathrm{eff}}\lesssim 0$ and
$0.45\lesssim\chi_{\text{p}}\lesssim 0.5$, and one, SXS:0165, has $Q\approx
6,\chi_{\mathrm{eff}}\approx-0.45,\chi_{\text{p}}\approx 0.77$. We should
caution that we detected “V”-shaped kinks in the phases of both the precessing
and co-precessing $(2,2),(2,1)$ multipoles. The first kinks appear at $\approx
100M$ before the peak of the multipole amplitude. Additionally, SXS:0057 is
very short, i.e., its inspiral lasts $\approx 800M$. For these reasons, we
present only the full wavelength mismatch, ${\bar{\cal M}}_{\text{opt}}$ for
these cases.

We list the values of ${\bar{\cal M}}_{\text{opt}}$ for
$M=37.5M_{\odot},150M_{\odot}$ and $\vartheta_{\text{LN},0}=0,\pi/2$ for each
model in Table 2. One interesting finding is that ${\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=\pi/2)<{\bar{\cal
M}}_{\text{opt}}(\vartheta_{\text{LN},0}=0)$ for all models except for the
$M=150M_{\odot}$ mismatches for SXS:0058. Most of the improved faithfulness at