# A catalogue of precessing black-hole-binary numerical-relativity simulations

Eleanor Hamilton  School of Physics and Astronomy, Cardiff University, Queens
Buildings, Cardiff, CF24 3AA, United Kingdom  Physik-Institut, Universität
Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland  Edward Fauchon-
Jones  School of Physics and Astronomy, Cardiff University, Queens Buildings,
Cardiff, CF24 3AA, United Kingdom  Imperial College London, South Kensington
Campus, London SW7 2AZ, United Kingdom  Mark Hannam  School of Physics and
Astronomy, Cardiff University, Queens Buildings, Cardiff, CF24 3AA, United
Kingdom  Charlie Hoy  School of Physics and Astronomy, Cardiff University,
Queens Buildings, Cardiff, CF24 3AA, United Kingdom  University of Portsmouth,
Portsmouth, PO1 3FX, United Kingdom  Chinmay Kalaghatgi  School of Physics and
Astronomy, Cardiff University, Queens Buildings, Cardiff, CF24 3AA, United
Kingdom  Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands
Institute for Gravitational and Subatomic Physics (GRASP), Utrecht University,
Princetonplein 1, 3584 CC Utrecht, The Netherlands  Institute for High-Energy
Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The
Netherlands  Lionel London  School of Physics and Astronomy, Cardiff
University, Queens Buildings, Cardiff, CF24 3AA, United Kingdom  Institute for
High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands  King’s College London, Strand, London WC2R 2LS,
United Kingdom  MIT-Kavli Institute for Astrophysics and Space Research and
LIGO Laboratory, 77 Massachusetts Avenue, 37-664H, Cambridge, MA 02139, USA
Jonathan E. Thompson  School of Physics and Astronomy, Cardiff University,
Queens Buildings, Cardiff, CF24 3AA, United Kingdom  Theoretical Astrophysics
Group, California Institute of Technology, Pasadena, CA 91125, USA  Dave
Yeeles  School of Physics and Astronomy, Cardiff University, Queens Buildings,
Cardiff, CF24 3AA, United Kingdom  Institute of Physics, Academia Sinica,
Taipei 115201, Taiwan  Shrobana Ghosh  School of Physics and Astronomy,
Cardiff University, Queens Buildings, Cardiff, CF24 3AA, United Kingdom
Sebastian Khan  School of Physics and Astronomy, Cardiff University, Queens
Buildings, Cardiff, CF24 3AA, United Kingdom  Panagiota Kolitsidou  School of
Physics and Astronomy, Cardiff University, Queens Buildings, Cardiff, CF24
3AA, United Kingdom  Alex Vano-Vinuales  School of Physics and Astronomy,
Cardiff University, Queens Buildings, Cardiff, CF24 3AA, United Kingdom
CENTRA, Departamento de Fìsica, Instituto Superior Técnico IST, Universidade
de Lisboa UL, Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal

###### Abstract

We present a public catalogue of numerical-relativity binary-black-hole
simulations. The catalogue contains datasets from 80 distinct configurations
of precessing binary-black-hole systems, with mass ratios up to
$m_{2}/m_{1}=8$, dimensionless spin magnitudes on the larger black hole up to
$|\vec{S}_{2}|/m_{2}^{2}=0.8$ (the small black hole is non-spinning), and a
range of five values of spin misalignment for each mass-ratio/spin
combination. We discuss the physical properties of the configurations in our
catalogue, and assess the accuracy of the initial configuration of each
simulation and of the gravitational waveforms. We perform a careful analysis
of the errors due to the finite resolution of our simulations and the finite
distance from the source at which we extract the waveform data and provide a
conservative estimate of the mismatch accuracy. We find that the upper limit
on the mismatch uncertainty of our waveforms is $0.4\%$. In doing this we
present a consistent approach to combining mismatch uncertainties from
multiple error sources. We compare this release to previous catalogues and
discuss how these new simulations complement the existing public datasets. In
particular, this is the first catalogue to uniformly cover this parameter
space of single-spin binaries and there was previously only sparse coverage of
the precessing-binary parameter space for mass ratios $\gtrsim 5$. We discuss
applications of these new data, and the most urgent directions for future
simulation work. The public dataset can be accessed online at
https://data.cardiffgravity.org/bam-catalogue/.

##  I Introduction

After several decades of research to solve the binary-black-hole (BBH)
problem, the first numerical relativity (NR) BBH simulations through one
orbit, merger and ringdown were produced in 2005 Pretorius (2005); Campanelli
et al. (2006); Baker et al. (2006). Since then many independent numerical
relativity codes Bruegmann et al. (2008); Husa et al. (2008a); Scheel et al.
(2006); Hemberger et al. (2013); Herrmann et al. (2007a); Zlochower et al.
(2005); Sperhake (2007); Loffler et al. (2012) have been developed to simulate
BBH systems for many orbits and added support for more complex configurations
such as extremely high mass ratios and highly spinning black holes Lousto and
Zlochower (2011); Sperhake et al. (2011); Scheel et al. (2015).

The data products from NR have been crucial for the field of gravitational-
wave astronomy, including as input to develop approximate gravitational wave
models Husa et al. (2016); Khan et al. (2016); Pratten et al. (2020); García-
Quirós et al. (2020); Estellés et al. (2020); Hamilton et al. (2021); Buonanno
et al. (2009); Taracchini et al. (2014); Pan et al. (2014); Bohé et al.
(2017); Babak et al. (2017); Cotesta et al. (2018); Blackman et al. (2017a,
b); Varma et al. (2019a), to calculate remnant properties of binary mergers
Gonzalez et al. (2007a, b); Herrmann et al. (2007b); Campanelli et al. (2007a,
b); Kesden (2008); Lousto and Zlochower (2013); Healy et al. (2014); Healy and
Lousto (2017, 2018), and used directly for gravitational wave injection
studies Aylott et al. (2009); Aasi et al. (2014). These NR-dependent tools
have played a central role in the direct detection of gravitational waves and
the measurement of their source properties over the last six years Abbott et
al. (2019a, 2021a, 2021b).

Gravitational wave observations during the first three LIGO-Virgo-Kagra (LVK)
observing runs Abbott et al. (2016a, b, c, 2017a, 2017b, 2017c, 2019a, 2020a,
2020b, 2020c, 2021a, 2021c, 2021d, 2021b) have relied on theoretical models
from three families; Phenom, SEOBNR and NRSurrogate. The Phenom and SEOBNR
families use NR waveforms to inform the merger-ringdown part of the model by
calibrating a theoretically-motivated ansatz to the numerical data; the
NRSurrogate models are constructed entirely from NR input. The simplest
gravitational wave models calibrated to NR data are aligned-spin models Husa
et al. (2016); Khan et al. (2016); Taracchini et al. (2014), which capture the
most important features of the waveform but require NR waveforms that cover
only a three-dimensional parameter space. More recent aligned spin models have
benefitted from further calibration to expanded data sets Pratten et al.
(2020); García-Quirós et al. (2020); Estellés et al. (2020); Bohé et al.
(2017); Cotesta et al. (2018). Many other subsequent models have been based on
these aligned-spin models, thus indirectly benefitting from calibration to NR
data London et al. (2018); Khan et al. (2019, 2020); Thompson et al. (2020);
Pratten et al. (2021); Estellés et al. (2021); Ossokine et al. (2020); Matas
et al. (2020). Recent expansions of NR catalogues to cover the precessing
parameter space have enabled the construction of the first generic-spin models
calibrated entirely to NR data Varma et al. (2019a).

NR has been useful beyond modelling. NR waveform injections have been used in
several studies, including to assess the presence of systematic bias in
waveform models Abbott et al. (2017d); Pürrer and Haster (2019), and to
estimate intermediate mass black hole binary merger rates Abbott et al.
(2019b). NR waveforms have also been used for direct comparisons and parameter
estimation of gravitational wave observations Abbott et al. (2016d); Lange et
al. (2017). Further, NR data can be used in the construction of fits to
predict the remnant properties of a BBH merger, namely the final mass and spin
as well as the gravitational recoil Rezzolla et al. (2008a, b); Tichy and
Marronetti (2007); Barausse and Rezzolla (2009); Barausse et al. (2012);
Lousto et al. (2012); Jiménez-Forteza et al. (2017); Varma et al. (2019b);
Zappa et al. (2019). These fits have a number of applications, such as tests
of general relativity Ghosh et al. (2018); Abbott et al. (2016e, 2021e,
2021f).

Several large catalogues of BBH NR simulations exist Mroue et al. (2013);
Boyle et al. (2019); Healy et al. (2017, 2019); Healy and Lousto (2020, 2022);
Jani et al. (2016). A quasi-circular BBH is described by 8 intrinsic
parameters; the masses of each of the black holes $m_{1}$ and $m_{2}$ and
their respective spins $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$. The total mass
$M=m_{1}+m_{2}$ sets the overall frequency scale and can be factored out. We
therefore choose to set $M=1$. The dimensionless spin is defined as
$\vec{\chi}_{i}=\mathbf{S}_{i}/m_{i}^{2}$. The majority of simulations
contained within these catalogues cover the precessing parameter space up to
mass ratio $q=m_{2}/m_{1}=8$ and dimensionless spin magnitudes $\chi\leqslant
0.8$. However, the existence of simulations beyond $q=4$ of sufficient length
and accuracy to be useful in the construction of gravitational wave models is
fairly sparse. There therefore exists no broad systematic covering of the
precessing parameter space up to $q=8$ with NR simulations. One purpose of the
current catalogue is to provide a systematic covering of that parameter space.

The primary objective of this catalogue was to support the development of a
new precessing phenomenological model that is calibrated to numerical
relativity waveforms Hamilton et al. (2021). Experience with producing
previous phenomenological models suggests that we do not require an extremely
dense sampling of the parameter space to produce a reasonably accurate model
Husa et al. (2016); Khan et al. (2016). For the first catalogue used to inform
the first precessing Phenom model we therefore chose no more than five points
in each parameter direction. This choice was found to be sufficient: the
PhenomPNR model of the dominant contribution to the signal (the $(2,2)$
multipoles in the co-precessing frame) Hamilton et al. (2021), constructed
from 19 aligned-spin waveforms and 40 precessing-binary waveforms, is of
comparable accuracy to the equivalent contributions to the NRSurrogate model
Varma et al. (2019a), which was constructed from more than 1000 simulations
over a smaller volume of parameter space. (It remains to be seen how many NR
simulations are required to accurately include two-spin effects, higher
multipoles, and mode asymmetries.) Secondary objectives were to contribute
data that is useful to the waveform modelling community and to provide
processed datasets that are appropriate for parameter estimation studies
Abbott et al. (2017d, 2019b).

Our catalogue contains datasets from 80 different configurations of precessing
BBH systems. These configurations cover four mass ratios
$q=m_{2}/m_{1}\in\\{1,2,4,8\\}$ at four different spin magnitudes
$\chi_{2}=|\vec{S}_{2}|/m_{2}^{2}\in\\{0.2,0.4,0.6,0.8\\}$ each at five
different spin vectors such that the angle between the orbital angular
momentum and spin vector of the larger black hole is one of
$\\{30^{\circ},60^{\circ},90^{\circ},120^{\circ},150^{\circ}\\}$. The
configurations are specified at a reference orbital frequency. The catalogue
can be accessed online at https://data.cardiffgravity.org/bam-catalogue/.

In the following section we briefly summarise the methods used by the BAM code
Bruegmann et al. (2008); Husa et al. (2008a) to perform numerical simulations
of BBH systems and describe the workflow we use to produce low eccentricity
initial data. In Sec. III we provide a description of the properties of the
simulations contained within the catalogue. In Sec. IV we perform a waveform
accuracy analysis to validate the catalogue. Finally we conclude with Sec. V
where we discuss what regions of parameter space and how the catalogue can be
used to contribute to the continuing advance of gravitational wave data
analysis.

##  II Summary of Methods

###  II.1 Simulation method

The simulations in this catalogue were produced using BAM Bruegmann et al.
(2008); Husa et al. (2008a), a moving-box-based mesh-refinement numerical-
relativity code that solves the $3+1$ decomposed Einstein equations.
Specifically, for the simulations in this catalogue we evolve Bowen-York
wormhole data Bowen and York (1980); Brandt and Bruegmann (1997); Ansorg et
al. (2004) via the $\chi$ variant of the moving-puncture treatment Campanelli
et al. (2006); Baker et al. (2006) of the BSSN formulation Baumgarte and
Shapiro (1998); Shibata and Nakamura (1995). Spatial derivatives are
approximated algebraically through sixth-order finite differencing in the
bulk, which in turn are evolved in time through fourth-order Runge-Kutta time
stepping; see Refs. Bruegmann et al. (2008); Husa et al. (2008a) for full
details on the treatment of boundaries, buffer zones, advection derivatives,
and numerical dissipation. Finally, the gravitational wave content of the
system is extracted at some finite distance using the Newman-Penrose scalar
$\Psi_{4}$ Newman and Penrose (1962) following the procedure outlined in Ref.
Baker et al. (2002). Following these references, $\Psi_{4}$ is decomposed via
the spherical multipolar decomposition

$\displaystyle\Psi_{4}\left(\Theta,\Phi\right)={}$
$\displaystyle\sum_{\ell=2}^{\infty}\sum_{m=\ell}^{\ell}\psi_{4,\ell
m}\,^{-2}Y_{\ell m}\left(\Theta,\Phi\right),$ (1)

where $\psi_{4,\ell m}$ are the individual multipole moments, ${}^{-2}Y_{\ell
m}$ are the spin-weighted spherical harmonics over the unit-sphere defined by
$\Theta$ and $\Phi$. The individual multipole moments can be written

$\displaystyle\psi_{4,\ell m}={}$ $\displaystyle A_{\ell m}e^{-im\phi_{\ell
m}},$ (2)

where $A_{\ell m}$ and $\phi_{\ell m}$ are the amplitude and phase of the
multipoles respectively. The quantity $\Psi_{4}$ can be converted into the
gravitational wave observable known as the strain $h$ via a double time
integral. This is generally the more useful quantity to consider for
applications of the numerical data, such as waveform modelling, since it is
the quantity measured by gravitational wave detectors and so forms the
starting point for all gravitational wave astronomy. Where the strain is
required, we obtain it in the frequency domain by dividing the frequency
domain $\Psi_{4}$ data by $\omega^{2}$, where $\omega$ is the angular Fourier
frequency Reisswig and Pollney (2011).

The numerical domain consists of nested Cartesian grids of successively finer
spacing, nested in the sense that the grid at each level $n$ is encompassed by
that of level $n-1$. The grid spacing $d_{n}$ for each refinement level $n$
follows the scaling,

$d_{n}=\tfrac{d_{0}}{2^{n}}\qquad\\{n\in\mathbb{Z}^{+}\\},$ (3)

where $d_{0}$ is the spacing on the coarsest level. The coarsest levels
(largest boxes) encompass both black holes and are fixed, while for the finer
levels (smaller boxes) there is a box around each black hole, and these boxes
move with the punctures. The boxes are initially specified as cubes, where the
user provides the number of points along one side for each level, i.e., if
$N_{n}$ is the number of points in each direction on level $n$, then the user
specifies a list of $N_{i}=\\{N_{0},N_{1},...,N_{\rm l_{\rm max}}\\}$, where
$l_{\rm max}$ is the finest level. During evolution the code will dynamically
adjust the number of points, in particular near merger when individual boxes
around each puncture will be merged when they are about to overlap. In
addition, $l_{\rm max}$ differs for each puncture so that $d_{{\rm
lmax},i}/m_{i}$ is approximately the same for each puncture. See Ref.
Bruegmann et al. (2008) for more details on the BAM grid structure, and Ref.
Purrer et al. (2012) for typical choices for numbers of levels and relative
box sizes. Sec. II.2.2 provides further details of the choices we made for the
grid configurations. The values $N_{i}$, $l_{\rm max}$ for each puncture, and
the coarsest grid spacing $d_{0}$ are all provided in the public data release.

The temporal resolution is subject to a Berger-Oliger refinement scheme in
which the spacing between successive time steps halves with each successive
level. The finest level consists of two grids, one centred on each puncture,
an arrangement that is maintained as we move up the levels so long as the
grids are not so large that they would overlap. These nested grids around each
individual puncture move with the punctures as they orbit. Beyond this level
the two grids are replaced by a single grid that encompasses all of the moving
boxes and is centered on the origin. It was found in Ref. Bruegmann et al.
(2008) that the Berger-Oliger timestepping becomes unstable on the coarsest
non-moving grids, and for these we revert to a single time step specified by a
Courant factor of 0.25 applied to the finest such grid. The details are the
same as those used in Refs. Bruegmann et al. (2008); Husa et al. (2008a),
except for the use of 0.25 rather than 0.5 for the Courant factor, which is
necessary to sufficiently reduce the timestepping error in long simulations.

A pseudo-spectral elliptic solver is used to calculate binary wormhole initial
data Ansorg et al. (2004), with eccentricity reduced to $<2\times 10^{-3}$
through a series of manual iterations of the linear momenta of the punctures
in the initial parameters. This process is described in more detail in Sec.
II.2.1.

###  II.2 Simulation workflow

####  II.2.1 Initial data construction

We wish our simulations to begin at a user specified reference orbital
frequency $M\Omega_{\mathrm{orb}}$ with spin vector $\bm{S}\equiv\bm{S}_{2}$
on the larger black hole (which we designate the secondary). The orientation
of $\bm{S}_{2}$ can be defined by the angle
$\theta_{\text{LS}}=\arccos(\hat{\bm{L}}_{\mathrm{N}}\cdot\hat{\bm{S}}_{2})$
between the spin vector $\bm{S}_{2}$ and Newtonian orbital angular momentum
vector $\bm{L}_{\mathrm{N}}$, and the angle
$\phi_{\mathrm{rS}}=\arccos(\hat{\bm{r}}\cdot\hat{\bm{S}}_{2\perp})$ between
the projection of the spin vector on to the orbital plane $\bm{S}_{2\perp}$
and the separation vector $\bm{r}$ from the larger black hole to the smaller.
The positions and momenta of the black holes consistent with these constraints
must then be determined at $M\Omega_{\mathrm{orb}}$, approximately chosen to
minimise eccentricity. Bowen-York wormhole data can then be generated from
these parameters. The main task of initial data construction is therefore
reduced to identifying the appropriate black-hole positions and momenta at
$M\Omega_{\mathrm{orb}}$. Two methods were used for the simulations in this
catalogue.

For simulations with $\chi_{2}\in\\{0.4,0.8\\}$ the initial data parameters
were determined by adapting the method used in previous work Hannam et al.
(2010); Schmidt et al. (2012); Husa et al. (2016). For this method the
physical parameters of the system $(q,\bm{S}_{1},\bm{S}_{2})$ are specified at
a much larger separation $D_{\mathrm{start}}$ than the NR simulations will
start at. The effective-one-body (EOB) equations of motion are then evolved up
to $M\Omega_{\mathrm{orb}}$ and the parameters at this frequency are used as
input to a Bowen-York initial data solver. However for precessing systems this
method does not allow the user to specify the exact system configuration
$(q,\bm{S}_{1},\bm{S}_{2})$ at $M\Omega_{\mathrm{orb}}$. During the course of
inspiral from $D_{\text{start}}$ to $M\Omega_{\mathrm{orb}}$ for the single-
spin precessing systems in this catalogue, the angle $\theta_{\rm LS}$ can be
seen to vary no more than $\sim 1^{\circ}$, while $\phi_{\mathrm{rS}}$
increases continuously (see Fig. 3 and Fig. 4 in Schmidt et al. (2015)). To
achieve a specific choice of $(\theta_{\rm LS},\phi_{\mathrm{rS}})$ at a
prescribed value of $M\Omega_{\mathrm{orb}}$, the method was extended with an
iterative refinement of the angle $\phi_{\mathrm{rS}}$ at $D_{\mathrm{start}}$
until the parameters at $M\Omega_{\mathrm{orb}}$ are within a suitable
tolerance of our desired values. Full details of this adapted method are given
in Appendix A.

The simulations $\chi_{2}\in\\{0.2,0.6\\}$ were performed later, and were able
to make use of a more recent method to produce low-eccentricity initial
parameters, as described in Ref. Ramos-Buades et al. (2019). This method
provides a post-Newtonian estimate of low-eccentricity parameters at a
prescribed orbital frequency, making it possible to specify $(\theta_{\rm
LS},\phi_{\mathrm{rS}})$ without the need for any iterative steps. This method
also supports additional iteration steps to further reduce the eccentricity
based on NR dynamics, however this additional iteration was not used for the
simulations in this catalogue. We instead relied upon the manual perturbation
approach outlined in Sec. II.2.1 to reduce eccentricity when using either
approach to obtain the initial data.

For all the NR configurations described in this work the azimuthal angle for
the spin vector $\mathbf{S}_{2}$ placed on the larger secondary component
black hole was chosen to be $\phi=0^{\circ}$ at $M\Omega_{\mathrm{orb}}$.

While the initial data parameters generated in Sec. II.2.1 will lead to low
eccentricity simulations, in general this will not be low enough to satisfy
our definition of a quasi-circular binary. We placed an upper limit on the
eccentricity at $2\times 10^{-3}$, based on the observation in Ref. Purrer et
al. (2012) that the puncture dynamics do not give reliable eccentricity
estimates below this value, due to gauge effects. A standard iterative method
to further reduce eccentricity is to perform a low resolution simulation for
$\sim$1000$M$, estimate the eccentricity, and make iterative small
perturbations to the momenta of the component black holes Hannam et al.
(2010); Purrer et al. (2012). The method employed to estimate the eccentricity
is described below. For most of the simulations in this work a perturbation of
0.1-0.8% is applied to the magnitude of the momenta. This is normally
sufficient to reduce the eccentricity below the desired threshold. However in
cases where this is not sufficient the radial component of the momenta is also
reduced by 25-75%. The eccentricity reduction procedure is performed using low
resolution simulations in order to reduce both the computing resources and
wall time required. Once initial data parameters are found that yield a
sufficiently low eccentricity then a high resolution production simulation is
performed using the same parameters. The higher resolution simulations tend to
have higher eccentricity than the associated low resolution simulations.
Consequently, a number of the simulations presented in this paper have
eccentricities marginally above the $2\times 10^{-3}$ threshold.

There are two different ways that eccentricity is estimated for the
simulations in this work. For the shorter iterative eccentricity reduction
simulations where the merger time is not known, the puncture separation $D$ is
fit using a quadratic function with data typically in the range $[200,700]M$
similar to the method described in Ref. Husa et al. (2008b). The eccentricity
is then estimated by the maximum absolute relative difference between the fit
and the data. For production simulations the eccentricity is estimated using a
fit that also incorporates the merger time Husa et al. (2008b).

In our production simulations we in general find that true eccentricity
differs from that calculated in our lower-resolution eccentricity-reduction
simulations. For a few of the cases in this catalogue the eccentricity of the
lower resolution simulation was below our $2\times 10^{-3}$ threshold, but the
eccentricity of the production simulation exceeded it, as can be seen in Tab.
1. Nonetheless, only a handful of cases have eccentricities above $3\times
10^{-3}$, and only one is close to $4\times 10^{-3}$ (CF_8).

####  II.2.2 Grid configurations

The simulations performed for this catalogue are all computationally
expensive, requiring $O(10^{5})$ CPU hours for each production run, and we do
not have the luxury of exhaustive experiments to identify a choice of
numerical grids that provides a good balance between computational efficiency
and physical accuracy. In 3D simulations of this scale it is impractical to
perform standard convergence tests where the grid spacing $d_{0}$ is halved
between successive runs, and indeed clean convergence has rarely been observed
in binary simulations with any code, and even given promising convergence
results for one binary configuration, there exists no robust algorithm to
determine the resolution requirements to guarantee clean convergence for a
second configuration. Sec. IV presents a convergence study of several of our
configurations. In this section we discuss the heuristic requirements we place
on our grid configurations, based on past experience with BAM binary
simulations.

Our first requirement is that the width of the smallest moving box following
each component black hole should be between 1.2 and 1.5 times the maximum
effective coordinate diameter of the apparent horizon of its respective black
hole before merger. This requirement is achieved by changing the values of the
grid spacing $d_{0}$ on the coarsest level, and the finest level that exists
for the larger black hole. The number of grid points $N_{L}$ on the finest
level can also be used to adjust the size of the finest box around the black
hole, if necessary, but we find in most cases that adjusting $d_{0}$ is
sufficient.

The second requirement is to have at least ten grid points per wavelength of
the $(4,4)$ multipole moment on the level where gravitational waves are
extracted. The maximum frequency is estimated by doubling the $(2,2)$ ringdown
frequency calculated by the aligned spin gravitational wave model PhenomD Husa
et al. (2016); Khan et al. (2016), using the parameters $(q,0,\chi_{2})$. In
precessing configurations the ringdown frequency will always be lower than
this estimate, and therefore this provides a conservative estimate of the
resolution requirements. One could use a more accurate estimate of the
ringdown frequency for each precessing configuration using, for example, the
method described in Refs. Hannam et al. (2014); Schmidt et al. (2015); London
and Fauchon-Jones (2019), but for this work we found no need to do this. The
required grid spacing on level $n$ where gravitational waves are extracted is
then approximated as $d_{n}\leq 1/(20f_{\mathrm{RD}})$. This requirement is
achieved by changing the values of the finest grid spacing $d_{0}$. If this
requirement cannot be satisfied on level $n$ and level $n$ is not the last
fixed box level, then the number of grid points $N_{n+1}$ on level $n+1$ is
increased until the box size $N_{n+1}\times h_{n+1}$ is large enough to
support gravitational wave extraction at the radius required. The use of much
larger numbers of points on the wave-extraction level means that the wave-
extraction resolution requirements are also a strong determinant of the
overall computational cost, along with the resolution requirements local to
the black hole.

For most configurations both of the requirements outlined in the preceding two
paragraphs can be satisfied. However it is not always possible to satisfy both
requirements, and for such cases the smallest box sizes and extraction level
grid spacing are balanced to achieve the best possible result.

##  III Simulation properties

In this section we discuss the properties of the simulations in our catalogue.
We first motivate our coverage of the single-spin parameter space, our choice
of starting frequency for each binary, and our procedure to estimate initial
black-hole momenta and spins to achieve quasi-circular inspiral with a
prescribed spin orientation. We then discuss in detail the accuracy with which
our desired configurations are achieved, in particular the accuracy of our
specification of the black-hole masses and spins, and the spin orientations.
Finally, we summarise the properties of the remnant black holes.

###  III.1 Simulation configurations

Name | $q$ | $\chi$ | $\theta_{\mathrm{LS}}\left({}^{\circ}\right)$ | $\chi_{\mathrm{eff}}$ | $\chi_{\rm p}$ | $D/M$ | $e$ | $M\Omega_{\mathrm{orb}}$ | $t_{\rm M}$ | $N_{\mathrm{orb}}$ | $M_{f}$ | $\chi_{f}$ | $v_{R}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  | $\left(\times 10^{-3}\right)$ | $\left(\times 10^{-2}\right)$ |  |  |  |  | $(\mathrm{kms}^{-1})$
CF_1 | $1$ | 0.2 (0.200) | $30.0\,(29.9)$ | $0.087\,(0.087)$ |  $0.100$ $(0.100)$ |  $11.6$ $(11.4)$ | $1.327$ | $2.28$ | $1710$ | $9.61$ | $0.949$ | $0.713$ | $188$
CF_2 | $1$ |  $0.2$ $(0.200)$ | $60.0\,(59.9)$ | $0.050\,(0.050)$ |  $0.173$ $(0.173)$ |  $11.6$ $(11.4)$ | $0.931$ | $2.28$ | $1685$ | $9.43$ | $0.950$ | $0.703$ | $96$
CF_3 | $1$ | 0.2 (0.200) | $90.0\,(89.8)$ | $0.000\,(0.000)$ |  $0.200$ $(0.200)$ |  $11.6$ $(11.4)$ | $1.084$ | $2.27$ | $1667$ | $9.25$ | $0.951$ | $0.688$ | $176$
CF_4 | $1$ | 0.2 (0.200) | $120.0\,(119.8)$ | $-0.050\,(-0.050)$ |  $0.173$ $(0.173)$ |  $11.6$ $(11.4)$ | $0.759$ | $2.28$ | $1628$ | $9.01$ | $0.953$ | $0.672$ | $329$
CF_5 | $1$ | 0.2 (0.200) | $150.0\,(149.9)$ | $-0.087\,(-0.086)$ |  $0.100$ $(0.100)$ |  $11.6$ $(11.4)$ | $0.723$ | $2.27$ | $1610$ | $8.86$ | $0.953$ | $0.661$ | $147$
CF_6 | $1$ | 0.4 (0.400) | $30.0\,(29.9)$ | $0.173\,(0.174)$ |  $0.200$ $(0.200)$ |  $11.6$ $(11.3)$ | $1.227$ | $2.28$ | $1768$ | $10.0$ | $0.946$ | $0.740$ | $120$
CF_7 | $1$ | 0.4 (0.400) | $60.0\,(59.9)$ | $0.100\,(0.101)$ |  $0.346$ $(0.346)$ |  $11.6$ $(11.3)$ | $2.478$ | $2.28$ | $1723$ | $9.69$ | $0.948$ | $0.720$ | $735$
CF_8 | $1$ | 0.4 (0.400) | $90.0\,(89.9)$ | $0.000\,(0.001)$ |  $0.400$ $(0.400)$ |  $11.5$ $(11.3)$ | $3.898$ | $2.30$ | $1626$ | $9.11$ | $0.951$ | $0.691$ | $757$
CF_9 | $1$ |  $0.4$ $(0.400)$ | $120.0\,(119.9)$ | $-0.100\,(-0.099)$ |  $0.346$ $(0.347)$ |  $11.6$ $(11.3)$ | $3.213$ | $2.29$ | $1577$ | $8.72$ | $0.954$ | $0.660$ | $93$
CF_10 | $1$ | 0.4 (0.400) | $150.0\,(149.9)$ | $-0.173\,(-0.173)$ |  $0.200$ $(0.201)$ |  $11.6$ $(11.4)$ | $2.675$ | $2.29$ | $1542$ | $8.45$ | $0.956$ | $0.634$ | $324$
CF_11 | $1$ | 0.6 (0.600) | $30.0\,(29.9)$ | $0.260\,(0.260)$ |  $0.300$ $(0.299)$ |  $10.0$ $(9.79)$ | $1.269$ | $2.79$ | $1064$ | $7.34$ | $0.943$ | $0.767$ | $189$
CF_12 | $1$ | 0.6 (0.601) | $60.0\,(59.8)$ | $0.150\,(0.151)$ |  $0.520$ $(0.519)$ |  $11.6$ $(11.3)$ | $1.446$ | $2.29$ | $1737$ | $9.81$ | $0.946$ | $0.740$ | $977$
CF_13 | $1$ | 0.6 (0.600) | $90.1\,(89.8)$ | $0.000\,(0.001)$ |  $0.600$ $(0.600)$ |  $11.5$ $(11.3)$ | $1.261$ | $2.30$ | $1625$ | $9.11$ | $0.950$ | $0.697$ | $1170$
CF_14 | $1$ |  $0.6$ $(0.600)$ | $120.0\,(119.8)$ | $-0.150\,(-0.149)$ |  $0.519$ $(0.521)$ |  $11.7$ $(11.4)$ | $1.200$ | $2.27$ | $1589$ | $8.65$ | $0.955$ | $0.651$ | $98$
CF_15 | $1$ | 0.6 (0.600) | $150.0\,(149.9)$ | $-0.260\,(-0.259)$ |  $0.300$ $(0.301)$ |  $11.6$ $(11.4)$ | $1.761$ | $2.29$ | $1480$ | $8.03$ | $0.958$ | $0.609$ | $159$
CF_16 | $1$ | 0.8 (0.801) | $30.0\,(29.9)$ | $0.346\,(0.347)$ |  $0.400$ $(0.399)$ |  $11.6$ $(11.3)$ | $2.958$ | $2.28$ | $1902$ | $10.9$ | $0.939$ | $0.792$ | $874$
CF_17 | $1$ | 0.8 (0.801) | $60.0\,(59.8)$ | $0.200\,(0.202)$ |  $0.693$ $(0.692)$ |  $11.6$ $(11.4)$ | $2.691$ | $2.27$ | $1832$ | $10.3$ | $0.943$ | $0.758$ | $1690$
CF_18 | $1$ | 0.8 (0.801) | $90.1\,(89.8)$ | $0.000\,(0.002)$ |  $0.800$ $(0.801)$ |  $11.6$ $(11.3)$ | $1.027$ | $2.29$ | $1644$ | $9.17$ | $0.950$ | $0.707$ | $1220$
CF_19 | $1$ | 0.8 (0.801) | $120.1\,(119.8)$ | $-0.200\,(-0.198)$ |  $0.692$ $(0.696)$ |  $11.6$ $(11.3)$ | $1.402$ | $2.30$ | $1503$ | $8.26$ | $0.955$ | $0.641$ | $1110$
CF_20 | $1$ | 0.8 (0.801) | $150.0\,(149.9)$ | $-0.347\,(-0.346)$ |  $0.399$ $(0.403)$ |  $11.6$ $(11.3)$ | $0.552$ | $2.32$ | $1374$ | $7.48$ | $0.959$ | $0.584$ | $394$
CF_21 | $2$ | 0.2 (0.200) | $30.0\,(29.9)$ | $0.115\,(0.116)$ |  $0.100$ $(0.100)$ |  $11.6$ $(11.4)$ | $1.718$ | $2.27$ | $1934$ | $10.6$ | $0.958$ | $0.680$ | $114$
CF_22 | $2$ | 0.2 (0.200) | $60.0\,(59.9)$ | $0.067\,(0.067)$ |  $0.173$ $(0.173)$ |  $11.6$ $(11.4)$ | $1.257$ | $2.28$ | $1887$ | $10.3$ | $0.959$ | $0.659$ | $374$
CF_23 | $2$ | 0.2 (0.200) | $90.0\,(89.8)$ | $0.000\,(0.000)$ |  $0.200$ $(0.200)$ |  $11.6$ $(11.4)$ | $0.655$ | $2.28$ | $1823$ | $9.84$ | $0.961$ | $0.629$ | $207$
CF_24 | $2$ | 0.2 (0.200) | $120.0\,(119.8)$ | $-0.067\,(-0.066)$ |  $0.173$ $(0.174)$ |  $11.6$ $(11.4)$ | $1.159$ | $2.29$ | $1759$ | $9.42$ | $0.963$ | $0.596$ | $287$
CF_25 | $2$ | 0.2 (0.200) | $150.0\,(149.9)$ | $-0.116\,(-0.115)$ |  $0.100$ $(0.100)$ |  $11.6$ $(11.4)$ | $0.981$ | $2.29$ | $1713$ | $9.13$ | $0.964$ | $0.569$ | $227$
CF_26 | $2$ | 0.4 (0.400) | $30.0\,(29.9)$ | $0.231\,(0.231)$ |  $0.200$ $(0.200)$ |  $11.6$ $(11.4)$ | $1.237$ | $2.29$ | $2005$ | $11.2$ | $0.954$ | $0.737$ | $209$
CF_27 | $2$ | 0.4 (0.400) | $60.1\,(59.9)$ | $0.133\,(0.134)$ |  $0.347$ $(0.346)$ |  $11.6$ $(11.5)$ | $2.131$ | $2.26$ | $1979$ | $10.8$ | $0.956$ | $0.700$ | $713$
CF_28 | $2$ | 0.4 (0.400) | $90.1\,(89.9)$ | $-0.001\,(0.001)$ |  $0.400$ $(0.400)$ |  $11.6$ $(11.4)$ | $0.580$ | $2.28$ | $1838$ | $9.91$ | $0.961$ | $0.646$ | $169$
CF_29 | $2$ | 0.4 (0.400) | $120.1\,(119.8)$ | $-0.134\,(-0.133)$ |  $0.346$ $(0.347)$ |  $11.6$ $(11.4)$ | $2.226$ | $2.30$ | $1692$ | $9.02$ | $0.964$ | $0.577$ | $609$
CF_30 | $2$ | 0.4 (0.400) | $150.1\,(149.9)$ | $-0.231\,(-0.231)$ |  $0.200$ $(0.201)$ |  $11.6$ $(11.4)$ | $1.436$ | $2.29$ | $1626$ | $8.53$ | $0.966$ | $0.518$ | $270$
CF_31 | $2$ | 0.6 (0.601) | $30.1\,(29.9)$ | $0.346\,(0.347)$ |  $0.301$ $(0.300)$ |  $11.5$ $(11.4)$ | $1.726$ | $2.28$ | $2118$ | $12.0$ | $0.948$ | $0.795$ | $154$
CF_32 | $2$ | 0.6 (0.600) | $60.1\,(59.8)$ | $0.199\,(0.201)$ |  $0.520$ $(0.519)$ |  $11.5$ $(11.3)$ | $0.876$ | $2.30$ | $1950$ | $10.9$ | $0.952$ | $0.746$ | $1280$
CF_33 | $2$ | 0.6 (0.601) | $90.1\,(89.8)$ | $-0.001\,(0.002)$ |  $0.600$ $(0.601)$ |  $11.7$ $(11.5)$ | $0.681$ | $2.25$ | $1908$ | $10.2$ | $0.958$ | $0.669$ | $1270$
CF_34 | $2$ | 0.6 (0.601) | $120.1\,(119.8)$ | $-0.201\,(-0.199)$ |  $0.519$ $(0.522)$ |  $11.5$ $(11.3)$ | $1.531$ | $2.31$ | $1610$ | $8.57$ | $0.965$ | $0.571$ | $390$
CF_35 | $2$ | 0.6 (0.600) | $150.1\,(149.9)$ | $-0.347\,(-0.346)$ |  $0.299$ $(0.301)$ |  $11.9$ $(11.6)$ | $0.525$ | $2.24$ | $1657$ | $8.39$ | $0.968$ | $0.471$ | $395$
CF_36 | $2$ | 0.8 (0.802) | $30.1\,(29.9)$ | $0.461\,(0.463)$ |  $0.401$ $(0.399)$ |  $11.5$ $(11.2)$ | $2.096$ | $2.31$ | $2156$ | $12.4$ | $0.940$ | $0.851$ | $866$
CF_37 | $2$ | 0.8 (0.802) | $60.1\,(59.7)$ | $0.265\,(0.269)$ |  $0.694$ $(0.692)$ |  $11.5$ $(11.3)$ | $2.046$ | $2.30$ | $2014$ | $11.4$ | $0.948$ | $0.800$ | $903$
CF_38 | $2$ | 0.8 (0.802) | $90.2\,(89.7)$ | $-0.002\,(0.003)$ |  $0.800$ $(0.802)$ |  $11.6$ $(11.4)$ | $3.134$ | $2.28$ | $1844$ | $10.0$ | $0.958$ | $0.707$ | $720$
CF_39 | $2$ | 0.8 (0.802) | $120.2\,(119.7)$ | $-0.268\,(-0.264)$ |  $0.691$ $(0.696)$ |  $11.5$ $(11.3)$ | $0.869$ | $2.32$ | $1554$ | $8.29$ | $0.964$ | $0.571$ | $1180$
CF_40 | $2$ | 0.8 (0.802) | $150.1\,(149.9)$ | $-0.462\,(-0.462)$ |  $0.398$ $(0.403)$ |  $11.6$ $(11.3)$ | $0.801$ | $2.32$ | $1408$ | $7.26$ | $0.969$ | $0.431$ | $443$
Table 1: Initial data parameters and relaxed properties of the precessing BBH configurations in this catalogue with mass ratio 1 or 2. The smaller black hole has no initial spin. The associated properties of the larger black hole are identified with a subscript $2$. The spin magnitude $S_{2}/m_{2}^{2}$, the spin tilt angle $\arccos(\hat{\bm{L}}_{\mathrm{N}}\cdot\hat{\bm{S}}_{2})$, the effective spins $\chi_{\mathrm{eff}}$ and $\chi_{p}$ and the separation $D/M$ are derived from the initial conditions of the simulations and relaxed times given in brackets. The eccentricity $e$ is estimated over the region $[200,1000]M$ using the method described in Husa et al. (2008b). The orbital frequency $M\omega_{\mathrm{orb}}$ is derived from the dynamics at relaxed times. The number of orbits $N_{\mathrm{orb}}$ is from calculated from the relaxed time that $M\omega_{\mathrm{orb}}$ is reported at until the peak in the $(2,2)$ multipole moment of $\Psi_{4}$.  Name | $q$ | $\chi$ | $\theta_{\mathrm{LS}}\left({}^{\circ}\right)$ | $\chi_{\mathrm{eff}}$ | $\chi_{\rm p}$ | $D/M$ | $e$ | $M\Omega_{\mathrm{orb}}$ | $t_{\rm M}$ | $N_{\mathrm{orb}}$ | $M_{f}$ | $\chi_{f}$ | $v_{R}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  | $\left(\times 10^{-3}\right)$ | $\left(\times 10^{-2}\right)$ |  |  |  |  | $(\mathrm{kms}^{-1})$
CF_41 | $4$ | 0.2 (0.200) | $30.0\,(29.9)$ | $0.139\,(0.139)$ |  $0.100$ $(0.100)$ |  $10.9$ $(10.9)$ | $1.083$ | $2.48$ | $1998$ | $11.4$ | $0.975$ | $0.566$ | $175$
CF_42 | $4$ | 0.2 (0.200) | $60.0\,(59.8)$ | $0.080\,(0.080)$ |  $0.173$ $(0.173)$ |  $11.0$ $(10.9)$ | $2.708$ | $2.47$ | $1982$ | $11.1$ | $0.976$ | $0.536$ | $154$
CF_43 | $4$ | 0.2 (0.200) | $90.1\,(89.9)$ | $0.000\,(0.000)$ |  $0.200$ $(0.200)$ |  $10.5$ $(10.5)$ | $1.386$ | $2.61$ | $1610$ | $9.34$ | $0.977$ | $0.488$ | $170$
CF_44 | $4$ | 0.2 (0.200) | $120.0\,(119.8)$ | $-0.080\,(-0.080)$ |  $0.173$ $(0.173)$ |  $11.2$ $(11.2)$ | $1.390$ | $2.38$ | $2004$ | $10.7$ | $0.978$ | $0.433$ | $241$
CF_45 | $4$ | 0.2 (0.200) | $150.0\,(149.9)$ | $-0.139\,(-0.138)$ |  $0.100$ $(0.100)$ |  $11.3$ $(11.3)$ | $0.701$ | $2.37$ | $1980$ | $10.4$ | $0.979$ | $0.385$ | $192$
CF_46 | $4$ | 0.4 (0.400) | $30.1\,(29.9)$ | $0.277\,(0.277)$ |  $0.201$ $(0.201)$ |  $10.6$ $(10.6)$ | $1.501$ | $2.58$ | $1947$ | $11.7$ | $0.972$ | $0.664$ | $234$
CF_47 | $4$ | 0.4 (0.400) | $60.2\,(59.9)$ | $0.159\,(0.160)$ |  $0.347$ $(0.347)$ |  $10.7$ $(10.7)$ | $1.424$ | $2.53$ | $1929$ | $11.2$ | $0.974$ | $0.616$ | $164$
CF_48 | $4$ | 0.4 (0.400) | $90.2\,(89.9)$ | $-0.001\,(0.000)$ |  $0.400$ $(0.400)$ |  $11.6$ $(11.5)$ | $1.608$ | $2.28$ | $2378$ | $12.4$ | $0.977$ | $0.533$ | $448$
CF_49 | $4$ | 0.4 (0.400) | $120.2\,(119.9)$ | $-0.161\,(-0.160)$ |  $0.346$ $(0.346)$ |  $11.6$ $(11.5)$ | $0.758$ | $2.30$ | $2148$ | $11.0$ | $0.980$ | $0.424$ | $336$
CF_50 | $4$ | 0.4 (0.400) | $150.1\,(149.9)$ | $-0.277\,(-0.277)$ |  $0.199$ $(0.200)$ |  $11.5$ $(11.5)$ | $1.761$ | $2.32$ | $1954$ | $9.87$ | $0.981$ | $0.313$ | $272$
CF_51 | $4$ | 0.6 (0.600) | $30.1\,(29.9)$ | $0.415\,(0.415)$ |  $0.301$ $(0.301)$ |  $10.4$ $(10.3)$ | $1.644$ | $2.63$ | $1982$ | $12.5$ | $0.967$ | $0.762$ | $413$
CF_52 | $4$ | 0.6 (0.601) | $60.2\,(59.8)$ | $0.238\,(0.239)$ |  $0.521$ $(0.521)$ |  $10.8$ $(10.7)$ | $1.817$ | $2.52$ | $2054$ | $12.1$ | $0.970$ | $0.704$ | $782$
CF_53 | $4$ | 0.6 (0.602) | $90.3\,(89.8)$ | $-0.002\,(-0.001)$ |  $0.600$ $(0.602)$ |  $11.0$ $(10.9)$ | $2.232$ | $2.46$ | $1935$ | $10.8$ | $0.975$ | $0.600$ | $764$
CF_54 | $4$ | 0.6 (0.601) | $120.3\,(119.8)$ | $-0.242\,(-0.241)$ |  $0.518$ $(0.520)$ |  $11.4$ $(11.4)$ | $0.438$ | $2.34$ | $1949$ | $10.0$ | $0.980$ | $0.451$ | $524$
CF_55 | $4$ | 0.6 (0.600) | $150.2\,(149.8)$ | $-0.416\,(-0.416)$ |  $0.298$ $(0.300)$ |  $10.1$ $(9.88)$ | $3.600$ | $2.88$ | $945$ | $5.45$ | $0.983$ | $0.270$ | $296$
CF_56 | $4$ | 0.8 (0.801) | $30.2\,(29.9)$ | $0.553\,(0.554)$ |  $0.402$ $(0.402)$ |  $10.1$ $(10.0)$ | $1.638$ | $2.72$ | $1951$ | $13.1$ | $0.959$ | $0.859$ | $722$
CF_57 | $4$ | 0.8 (0.802) | $60.4\,(59.7)$ | $0.317\,(0.318)$ |  $0.695$ $(0.696)$ |  $10.5$ $(10.4)$ | $0.751$ | $2.60$ | $1968$ | $12.2$ | $0.965$ | $0.800$ | $1150$
CF_58 | $4$ | 0.8 (0.802) | $90.4\,(89.7)$ | $-0.004\,(-0.002)$ |  $0.800$ $(0.802)$ |  $11.5$ $(11.4)$ | $1.675$ | $2.31$ | $2342$ | $12.3$ | $0.973$ | $0.684$ | $1160$
CF_59 | $4$ | 0.8 (0.802) | $120.4\,(119.7)$ | $-0.324\,(-0.323)$ |  $0.690$ $(0.693)$ |  $11.6$ $(11.5)$ | $1.225$ | $2.30$ | $1984$ | $9.92$ | $0.979$ | $0.506$ | $710$
CF_60 | $4$ | 0.8 (0.801) | $150.2\,(149.8)$ | $-0.555\,(-0.556)$ |  $0.397$ $(0.399)$ |  $11.6$ $(11.5)$ | $0.565$ | $2.31$ | $1690$ | $8.18$ | $0.983$ | $0.271$ | $383$
CF_61 | $8$ |  $0.2$ $(0.200)$ | $30.0\,(29.9)$ | $0.154\,(0.154)$ |  $0.100$ $(0.100)$ |  $9.84$ $(9.96)$ | $1.359$ | $2.86$ | $2006$ | $12.4$ | $0.988$ | $0.437$ | $81$
CF_62 | $8$ |  $0.2$ $(0.200)$ | $60.0\,(59.8)$ | $0.089\,(0.089)$ |  $0.173$ $(0.173)$ |  $9.95$ $(10.1)$ | $0.656$ | $2.82$ | $2008$ | $12.1$ | $0.988$ | $0.402$ | $87$
CF_63 | $8$ |  $0.2$ $(0.200)$ | $90.0\,(89.7)$ | $0.000\,(0.000)$ |  $0.200$ $(0.200)$ |  $10.1$ $(10.2)$ | $1.147$ | $2.77$ | $1989$ | $11.7$ | $0.989$ | $0.345$ | $98$
CF_64 | $8$ |  $0.2$ $(0.200)$ | $120.0\,(119.8)$ | $-0.089\,(-0.089)$ |  $0.173$ $(0.174)$ |  $10.3$ $(10.4)$ | $0.739$ | $2.70$ | $2001$ | $11.3$ | $0.990$ | $0.273$ | $100$
CF_65 | $8$ |  $0.2$ $(0.200)$ | $150.0\,(149.9)$ | $-0.154\,(-0.154)$ |  $0.100$ $(0.100)$ |  $10.4$ $(10.5)$ | $1.909$ | $2.67$ | $1987$ | $11.0$ | $0.990$ | $0.199$ | $100$
CF_66 | $8$ | 0.4 (0.400) | $30.2\,(29.9)$ | $0.307\,(0.309)$ |  $0.201$ $(0.199)$ |  $9.54$ $(9.67)$ | $1.360$ | $2.96$ | $2022$ | $13.3$ | $0.986$ | $0.573$ | $118$
CF_67 | $8$ |  $0.4$ $(0.400)$ | $60.3\,(59.9)$ | $0.176\,(0.178)$ |  $0.347$ $(0.346)$ |  $9.69$ $(9.84)$ | $1.921$ | $2.90$ | $1960$ | $12.4$ | $0.987$ | $0.525$ | $83$
CF_68 | $8$ | 0.4 (0.400) | $90.3\,(89.8)$ | $-0.002\,(0.000)$ |  $0.400$ $(0.400)$ |  $10.1$ $(10.2)$ | $1.222$ | $2.76$ | $2007$ | $11.8$ | $0.988$ | $0.439$ | $193$
CF_69 | $8$ | 0.4 (0.400) | $120.3\,(119.8)$ | $-0.179\,(-0.178)$ |  $0.345$ $(0.346)$ |  $10.5$ $(10.6)$ | $1.402$ | $2.65$ | $2022$ | $11.1$ | $0.990$ | $0.318$ | $166$
CF_70 | $8$ | 0.4 (0.400) | $150.1\,(149.9)$ | $-0.308\,(-0.308)$ |  $0.199$ $(0.200)$ |  $10.8$ $(10.8)$ | $1.578$ | $2.56$ | $2040$ | $10.7$ | $0.991$ | $0.169$ | $103$
CF_71 | $8$ | 0.6 (0.601) | $30.2\,(29.9)$ | $0.461\,(0.463)$ |  $0.302$ $(0.299)$ |  $9.28$ $(9.34)$ | $2.038$ | $3.07$ | $2040$ | $14.3$ | $0.983$ | $0.711$ | $161$
CF_72 | $8$ | 0.6 (0.601) | $60.4\,(59.8)$ | $0.263\,(0.268)$ |  $0.522$ $(0.519)$ |  $9.60$ $(9.71)$ | $2.220$ | $2.92$ | $2046$ | $13.3$ | $0.985$ | $0.658$ | $330$
CF_73 | $8$ | 0.6 (0.601) | $90.5\,(89.8)$ | $-0.004\,(0.001)$ |  $0.600$ $(0.601)$ |  $10.0$ $(10.1)$ | $2.370$ | $2.78$ | $1990$ | $11.8$ | $0.988$ | $0.561$ | $320$
CF_74 | $8$ | 0.6 (0.600) | $120.4\,(119.7)$ | $-0.270\,(-0.271)$ |  $0.518$ $(0.517)$ |  $10.6$ $(10.7)$ | $1.100$ | $2.60$ | $2008$ | $10.8$ | $0.990$ | $0.420$ | $124$
CF_75 | $8$ | 0.6 (0.601) | $150.2\,(149.9)$ | $-0.463\,(-0.462)$ |  $0.298$ $(0.301)$ |  $10.6$ $(10.6)$ | $0.482$ | $2.64$ | $1637$ | $8.57$ | $0.991$ | $0.247$ | $116$
CF_76 | $8$ | 0.8 (0.802) | $30.3\,(29.8)$ | $0.614\,(0.619)$ |  $0.404$ $(0.397)$ |  $8.96$ $(8.97)$ | $1.145$ | $3.20$ | $2044$ | $15.5$ | $0.977$ | $0.848$ | $210$
CF_77 | $8$ | 0.8 (0.801) | $60.5\,(59.7)$ | $0.350\,(0.362)$ |  $0.697$ $(0.690)$ |  $9.37$ $(9.36)$ | $1.334$ | $3.04$ | $2013$ | $13.8$ | $0.982$ | $0.799$ | $193$
CF_78 | $8$ | 0.8 (0.801) | $90.6\,(89.6)$ | $-0.008\,(0.003)$ |  $0.800$ $(0.801)$ |  $10.0$ $(10.0)$ | $2.868$ | $2.80$ | $2006$ | $11.9$ | $0.987$ | $0.700$ | $363$
CF_79 | $8$ | 0.8 (0.802) | $120.5\,(119.7)$ | $-0.361\,(-0.359)$ |  $0.689$ $(0.693)$ |  $10.9$ $(10.9)$ | $2.327$ | $2.50$ | $2145$ | $11.1$ | $0.990$ | $0.548$ | $207$
CF_80 | $8$ | 0.8 (0.802) | $150.3\,(149.8)$ | $-0.618\,(-0.618)$ |  $0.397$ $(0.398)$ |  $11.5$ $(11.4)$ | $0.701$ | $2.36$ | $2169$ | $10.3$ | $0.991$ | $0.371$ | $178$
Table 2: Initial data parameters and relaxed properties of the precessing BBH
configurations in this catalogue with mass ratio 4 or 8. The smaller black
hole has no initial spin. The associated properties of the larger black hole
are identified with a subscript $2$. The spin magnitude $S_{2}/m_{2}^{2}$, the
spin tilt angle $\arccos(\hat{\bm{L}}_{\mathrm{N}}\cdot\hat{\bm{S}}_{2})$, the
effective spins $\chi_{\mathrm{eff}}$ and $\chi_{p}$ and the separation $D/M$
are derived from the initial conditions of the simulations and relaxed times
given in brackets. The eccentricity $e$ is estimated over the region
$[200,1000]M$ using the method described in Husa et al. (2008b). The orbital
frequency $M\omega_{\mathrm{orb}}$ is derived from the dynamics at relaxed
times. The number of orbits $N_{\mathrm{orb}}$ is from calculated from the
relaxed time that $M\omega_{\mathrm{orb}}$ is reported at until the peak in
the $(2,2)$ multipole moment of $\Psi_{4}$.

Our catalogue consists of dynamics and waveform data from NR simulations of 80
binary-black-hole configurations. We chose configurations with four mass
ratios $q=m_{2}/m_{1}\in\\{1,2,4,8\\}$, four values of the dimensionless spin
on the larger black hole,
$\chi_{2}=|\bm{S}_{2}|/m_{2}^{2}\in\\{0.2,0.4,0.6,0.8\\}$ (the smaller black
hole has zero spin), and five values of the misalignment of the black-hole
spin with direction of the Newtonian orbital angular momentum,
$\theta_{\text{LS}}=\arccos(\hat{\bm{L}}_{\mathrm{N}}\cdot\hat{\bm{S}}_{2})\in\\{30^{\circ},60^{\circ},90^{\circ},120^{\circ},150^{\circ}\\}$.
The configurations with spins $\chi_{2}=\\{0.4,0.8\\}$ were simulated first
and used to produce the first PhenomPNR model, and the later simulations at
$\chi_{2}=\\{0.2,0.6\\}$ were used as independent verification waveforms
Hamilton et al. (2021).

For each simulation two mass parameters $m_{1},m_{2}$ were chosen such that
$M=m_{1}+m_{2}=1$. The initial data are iteratively constructed from these
parameters such that the ADM mass of each puncture equals its respective mass
parameter to within 0.02% Bruegmann et al. (2008). At subsequent times the
masses of each black hole are recorded as the apparent-horizon masses
$M_{\mathrm{AH},1},M_{\mathrm{AH},2}$ of each puncture, which are related to
the black hole hole masses $m_{1},m_{2}$ through the Christodoulou formula
Christodoulou (1970). This approach agrees well with the ADM mass of each
puncture; the level of agreement is quantified further in Sec. III.2.

In addition to the masses and spins, we must also choose the initial
separation of the binary. For a binary undergoing non-eccentric inspiral there
is a one-to-one correspondence between the black-hole separation and the
orbital frequency, so we may alternatively specify the initial orbital
frequency, $M\Omega_{\rm orb}$. For this catalogue we prefer to choose
$M\Omega_{\rm orb}$, because our primary purpose is to construct a frequency-
domain waveform model, and it would be convenient if we were able to start to
NR tuning at the same frequency for each configuration. This also motivates
the iterative procedure described in Sec. II.2.1, with the goal of finding
parameters consistent with non-eccentric inspiral for a configuration defined
at a specified starting frequency. For the first simulations we performed, at
mass ratios $q=1$ and $q=2$, we chose $M\Omega_{\rm orb}=0.0225$. This value
was chosen to produce simulations of $\sim$2000$M$ in length, which we
expected to be sufficiently accurate for our modelling purposes, based on the
experience of producing the aligned-spin simulations in Refs. Husa et al.
(2016); Khan et al. (2016).

The duration of the simulations varies with binary mass ratio and the
magnitude of the component of the spin aligned with the orbital angular
momentum. At leading post-Newtonian order the merger time from a given
starting frequency scales with $\Delta T\sim 1/\eta$, where
$\eta=m_{1}m_{2}/M^{2}$ is the symmetric mass ratio. Therefore, if simulations
at mass-ratios $q=1$ ($\eta=0.25$) and $q=8$ ($\eta\approx 0.1$) start at the
same orbital frequency, the $q=8$ simulation will take roughly 2.5 times as
long to merge. (This is a first-order approximation, and we see in the final
results that the variation is not quite so extreme.) In addition, if the
black-hole spin is aligned with the orbital angular momentum, the binary will
inspiral more slowly, and this will also increase the time to merger.
Conversely, a spin in the opposite direction to the orbital angular momentum
will decrease the time to merger. This effect of spin is most easily seen in
PN calculations, e.g., Refs. Cutler and Flanagan (1994); Poisson and Will
(1995). For example, for our $q=2$ configurations with $\chi_{2}=0.8$, where
all simulations begin at $M\Omega=0.0225$, we see that the $\theta_{\rm
LS}=30^{\circ}$ configuration merges in $2254M$, while the $\theta_{\rm
LS}=150^{\circ}$ configuration merges in only $1505M$. Since we do not wish to
perform expensive tests on the resolution requirements to achieve similar
levels of accuracy for much longer simulations, for mass ratios $q=4$ and
$q=8$, we adjust the starting frequency to limit the time to merger to
approximately $2000M$.

To meet the soft requirement of simulation merger by $2000M$ we estimate the
merger time using the LALSimulation LIGO Scientific Collaboration (2018)
implementation of PhenomD Khan et al. (2016). This provides a utility function
XLALSimIMRPhenomDChirpTime that calculates the time until the peak in the
(2,2)-strain of a specific system configuration given a starting gravitational
wave frequency, which is approximately twice the orbital frequency. The
starting frequency is optimized using a simple interval bisection procedure
until the peak time is $\sim 2000M$. A lower bound on the orbital frequency is
set at 0.0225. The average retarded merger time (calculated as detailed below)
for the simulations that required a higher starting orbital frequency was
$1983M$, with a minimum of $1610$, and a maximum of $2169M$. While
XLALSimIMRPhenomDChirpTime performed sufficiently well, overall it slightly
under estimated the merger time. One simulation (CF_55) was mistakenly
performed using an increased starting frequency, resulting in a much shorter
simulation with a retarded merger time of just $945M$.

The properties of each simulation are presented in Tabs. 1 and 2. Each
configuration is characterised by its mass ratio $q$, the dimensionless spin
magnitude $\chi_{2}=S_{2}/m_{2}^{2}$, the spin angle
$\theta_{\text{LS}}=\arccos(\hat{\bm{L}_{\mathrm{N}}}\cdot\hat{\bm{S}}_{2})$,
the initial orbital frequency $M\Omega_{\rm orb}$ (or alternatively the
initial binary separation $D/M$), and the binary’s eccentricity, $e$. For the
final values reported in Tab. 1 and 2, the eccentricity $e$ is estimated over
the region $[300,800]M$ using the method described in Husa et al. (2008b). We
also show the effective spin parameters $\chi_{\rm eff}$ and $\chi_{\rm p}$.
The effective aligned spin $\chi_{\rm eff}$ is defined in terms of the
individual parallel spin components
$\chi_{i}^{\parallel}=\bm{\chi}_{i}\cdot\hat{\mathbf{L}}_{\rm N}$ as Ajith et
al. (2011),

$\chi_{\rm
eff}=\frac{m_{1}\chi_{1}^{\parallel}+m_{2}\chi_{2}^{\parallel}}{M},$ (4)

and parameterises the dominant spin effect on the orbital phasing, as
discussed in Refs. Cutler and Flanagan (1994); Poisson and Will (1995); Baird
et al. (2013). The effective precession spin $\chi_{\rm p}$ is defined as
Schmidt et al. (2015),

$\chi_{\rm p}=\frac{S_{\rm p}}{m_{1}^{2}},$ (5)

where $A_{1}S_{\rm
p}=\text{max}\left(A_{1}S_{1}^{\perp},A_{2}S_{2}^{\perp}\right)$,
$A_{1}=2+3m_{2}/(2m_{1})$, and $A_{2}=2+3m_{1}/(2m_{2})$. In a generic two-
spin system the dominant precession effect can be approximated by a single-
spin system where the larger black hole has an in-plane spin of $\chi_{\rm
p}$, based on the leading-order spin precession effects Apostolatos et al.
(1994); Kidder (1995). In the single-spin configurations in this catalogue, we
will always have $\chi_{\rm eff}=m_{2}\chi_{2}\cos\theta_{\rm LS}/M$ and
$\chi_{\rm p}=\chi_{2}\sin\theta_{\rm LS}$.

In Tab. 1 we provide the quantities $(\chi,\theta_{\rm
LS},\chi_{\mathrm{eff}},\chi_{\rm p},D/M)$ as specified in the initial data
and (in brackets) at a _relaxed time_ , $t_{\mathrm{rel}}$. This is the time
at which we estimate that the unphysical junk radiation in the initial data
have radiated away, and the GW data can be used for analysis and modelling. We
wish $t_{\mathrm{rel}}$ to be as early as possible, to maximise the length of
the usable waveform. We choose a relaxed time of,

$t_{\mathrm{rel}}=t_{\mathrm{peak}}+2t_{\mathrm{damp}},$ (6)

where $t_{\mathrm{peak}}$ is the time of the peak amplitude of the junk
radiation in the $(2,2)$ multipole moment of $\Psi_{4}$ and
$t_{\mathrm{damp}}$ is an estimate of the exponential decay time of the junk
radiation, which we estimate as $t_{\mathrm{damp}}=76m_{2}$. The damping time
of the $(2,2,0)$ quasi-normal mode for a nonspinning black hole of mass $m$ is
approximately $71m$, and approximately $83m$ for a black hole with
dimensionless spin magnitude 0.8 Berti et al. (2006); we find that $76m$ is a
reasonable choice for all of the configurations in this catalogue. In Tab. 1
the initial orbital frequency $M\Omega_{\rm orb}$ is calculated at
$t_{\mathrm{rel}}$ and the number of orbits $N_{\rm orb}$ is calculated from
$t_{\mathrm{rel}}$ to the merger time, which we define to be the time at which
the peak in the $\ell=2$ multipole moments of $\Psi_{4}$ occurs. The retarded
merger time, which we denote as $t_{\rm M}$, is given by the difference
between the merger time and the tortoise co-ordinate

$\displaystyle r^{*}={}$ $\displaystyle r+2\ln\left|\frac{r}{2}-1\right|,$ (7)

where $r$ is the distance from the punctures at which the $\Psi_{4}$ data is
extracted.

The final black hole has a mass of $M_{f}$, a spin of $\chi_{f}$ and a recoil
velocity $v_{R}$. We discuss the calculation of these remnant quantities in
more detail in Sec. III.4 below.

###  III.2 Initial black-hole masses

We estimate the black-hole masses using the ADM mass calculated at the
puncture. The black holes are represented in the initial data as wormholes,
and the ADM mass calculated at the second asymptotically flat end of each
wormhole provides a good estimate of that black hole’s mass. In the puncture
framework, this mass estimate is easy to calculate at each black hole’s
puncture Brandt and Bruegmann (1997). The ADM puncture mass agrees well with
the mass calculated from the area of the apparent horizon in the case of
nonspinning black holes Tichy and Bruegmann (2004), but becomes less accurate
for high spins Hannam et al. (2010).

Figure 1 demonstrates the effect spin magnitude has on the initial data ADM
puncture mass as a function of time from the start of the simulation. The
results in Appendix A of Ref. Hannam et al. (2010) suggest that the error in
the mass estimate could be on the order of $\sim$0.5% for black holes with
spin 0.8. However, the estimates in that paper were made on the initial data.
Our results suggest that after the junk radiation has left the system (most
radiates to infinity, but some falls back into the black hole), the apparent-
horizon estimate of the mass is closer to the original ADM-puncture-mass
estimate. We see that for black holes with spins of 0.4 the error due to using
the ADM puncture mass is on the order of $\sim$0.01%, while for spins of 0.8
it is $\sim$0.04%. (The oscillations in Fig. 1 are due to uncertainties in the
apparent-horizon estimate.) From this we conclude that the errors in the mass
estimates are negligible.

Figure 1: Comparison of the relative percentage error between the apparent-
horizon mass and the initial data ADM puncture mass for the larger black hole,
as a function of simulation time $t$. Both simulations have initial parameters
$q=2$, $\theta_{\text{LS}}=60$. The solid black line is for a configuration
with dimensionless spin magnitude $\chi_{2}=0.4$ and the dashed black line
represents $\chi_{2}=0.8$.

###  III.3 Initial black-hole spins

The black-hole spin $\mathbf{S}_{2}$ is specified as part of the Bowen-York
extrinsic curvature. The main source of uncertainty in the dimensionless spin
$\chi_{2}=\mathbf{S}_{2}/m_{2}^{2}$ is the accuracy of the mass: as some of
the junk radiation falls into the black hole, the mass increases, and so
$\chi_{2}$ decreases. However, as we saw previously, the final value of the
mass as estimated from the area of the apparent horizon agrees well with the
nominal value for each configuration. We also see in Tab. 1 that there is only
a small discrepancy in the spin magnitude after the relaxation time. Since the
initial black-hole spins are prescribed analytically in the Bowen-York initial
data, we can reliably estimate the uncertainty in the apparent-horizon
measurement of the spin magnitude to be within $\sim$0.001.

During the inspiral $\theta_{\text{LS}}$ is not constant; it will oscillate,
as illustrated for one configuration in Fig. 2. Ideally, we would set up our
simulations so that the mean value of $\theta_{\rm LS}$ was equal to our
prescribed value at the start frequency. We see in Fig. 2 the two ways in
which our data deviate from this ideal. (1) There is an inaccuracy in the
initial value of $\theta_{\rm LS}$, which is within the tolerance set in our
initial-data construction procedure, and (2) this value is at an extremum of
the oscillations in $\theta_{\rm LS}$, and so the mean will be offset from the
target value in the initial-data construction. We also see that the mean value
slowly varies over the course of the simulation, although typically by only a
fraction of a degree over the entire inspiral.

In Tabs. 1 and 2 we report the mean value of $\theta_{\rm LS}$ at the start of
the simulation and at the relaxed time. We estimate the value at the relaxed
time by fitting to $\theta_{\rm LS}$ a sinusoidal ansatz of the form,

$\Big{(}A_{0}+A_{1}t\Big{)}\sin\Big{(}2\pi\big{(}f_{0}t+f_{1}t^{2}\big{)}+\varphi\Big{)}+C_{0}+C_{1}t,$
(8)

where $A_{0}$, $A_{1}$, $f_{0}$, $f_{1}$, $\varphi$, $C_{0}$ and $C_{1}$ are
all free parameters, from the relaxed time up to three orbits after the
relaxed time. The value of the linear part of the fit at the relaxed time is
reported in Tabs. 1 and 2 instead of the pointwise value of the NR data for
$\theta_{\rm LS}$. An example of this fit can be seen in Fig. 2. The resulting
value of $\theta_{\rm LS}$ is used, along with the relaxed-time value of the
spin magnitude, to calculate the relaxed-time values of $\chi_{\rm eff}$ and
$\chi_{\rm p}$.

Figure 2: Linear trend of the angle $\theta_{\text{LS}}$ between the Newtonian
orbital angular momentum and spin of the spinning black hole as a function of
simulation time $t$. $\theta_{\text{LS}}$ as determined from NR data is
plotted as a thick black line. The initial data value of $\theta_{\text{LS}}$
is plotted as a horizontal thick black line. The dashed black line is a
sinusoidal fit using Eq. 8. The vertical black lines mark the fit bounds. The
linear part of the sinusoidal fit is plotted as a solid black line. The dotted
black line marks the value of the linear trend line at the lower bound. Lower
opacity lines of the NR data and extrapolated fit are plotted outside of the
fit region.

###  III.4 Remnant properties

The final black hole that remains after the merger is characterised by its
mass, spin, and recoil. We report each of these quantities in Tabs. 1 and 2.

As with the relaxed-time quantities reported in Tab. 1 and 2, the mass and
spin of the final black hole, $M_{f}$ and $\chi_{f}$, are calculated from the
apparent horizonAlcubierre (2008). As a consistency check we also estimate the
mass and angular momentum of the final spacetime from the gravitational-wave
signal. The mass can be calculated by subtracting the radiated energy from the
initial ADM mass of the spacetime. The radiated energy is in turn calculated
from the gravitational wave signal measured at a series of extraction radii
Bruegmann et al. (2008) and the result extrapolated to infinity. The final
mass estimated in this way typically agrees with the horizon measure to within
$5\times 10^{-4}$. Given the mass, perturbation theory provides a relationship
between the black-hole spin and the frequency of the signal multipoles during
ringdown London and Fauchon-Jones (2019); Leaver (1985). We calculate the
ringdown frequency of the $(\ell=2,|m|=2)$ multipoles by taking the Fourier
transform of the waveform between $10$ and $100M$ after merger London (2015),
where merger is here defined as the time at which the sum of the square of the
$\ell=2$ multipoles is maximised. The ringdown frequency is then the frequency
at which the peak in this frequency domain post-merger waveform occurs London
et al. (2014). This then allows us to make an independent estimate of the
final spin. We find that this estimate of the final spin typically agrees with
the horizon measure to within $5\times 10^{-3}$.

We calculated the recoil (or kick) velocity of the final black hole by
integrating the radiated linear momentum from the relaxed time
$t_{\text{rel}}$ until the end of the simulation. We report here only the
magnitude of the the recoil velocity $v_{R}$ . The linear momentum is itself
calculated as described in Bruegmann et al. (2008). We used the value of the
linear momentum extracted at a distance $90M$ from the source. Since the
recoil velocity is very sensitive to the in-plane spin directions, this
catalogue does not comprehensively explore the range of recoil velocities that
can be seen for systems with mass ratios up to $q=8$ and dimensionless spin
magnitudes up to $\chi=0.8$. However, from the values presented in Tabs. 1 and
2, we can see that the largest magnitude kick velocities tend to be seen for
systems where the two initial black holes are equal in mass, a general trend
that can be seen even in inspiral post-Newtonian estimates Kidder (1995), and
for full merger calculations in the numerical-relativity recoil studies cited
in the Introduction.

Figure 3: A comparison of the calculation of the recoil velocity from the
radiated linear momentum calculated by the BAM code (black crosses) with the
prediction by the NRSur7dq4 model (blue dots).

We investigated the effect of the finite resolution of the simulations and the
radius at which the rate of change of the linear momentum was extracted on the
calculation of the final recoil velocity and found it to be negligible. The
effect of the extraction radius was found to be less than $10\%$ of the final
recoil velocity, while the effect of the resolution was even smaller.

We compared the results of our calculation of the recoil velocity with the
prediction given by the NRSurrogate model NRSur7dq4 Varma et al. (2019a) for
those cases within the catalogue that lie within the calibration region of
NRSur7dq4 ($q\leq 4$). To obtain this prediction, we used the value of the
black hole spins, rotated into the LAL frame Schmidt et al. (2017), and the
orbital frequency $100M$ prior to merger. This comparison is shown in Fig. 3.
As can be seen from these results, for most of the cases contained within the
catalogue, the calculation from the radiated linear angular momentum agrees
well with the prediction by NRSur7dq4. However, in a small number of cases
(most notably CF_37 and CF_38) the two values differ by around $50\%$ of the
value calculated from the radiated linear momentum. However, these values
remain within the bounds predicted by NRSur7dq4 for an equivalent
configuration but with a different value for the in-plane spin angle. We
therefore do not find these discrepancies too concerning, and we leave
determining their exact cause to a future investigation.

##  IV Waveform accuracy

In order to assess the accuracy of the data that comprise this catalogue we
studied a subset of four of the configurations described in Tabs. 1 and 2.
These four configurations are CF_47, CF_59, CF_66, and CF_80, with physical
parameters
$\left(q,\chi,\theta_{\text{LS}}\right)=\left\\{\left(4,0.4,60\right),\left(4,0.8,120\right),\left(8,0.4,30\right),\left(8,0.8,150\right)\right\\}$.
The set of simulations used in the accuracy analysis of the
$\left(4,0.4,60\right)$ case were performed with a lower starting frequency of
$M\omega_{\text{orb}}=0.023$ to provide an assessment of the accuracy of a
longer simulation.

The two main sources of error in our waveforms are the finite resolution of
the simulation and the finite radius at which the data are extracted. In order
to assess the effect of the finite resolution, we performed a set of three
simulations with low, medium and high resolution for each of the four
configurations listed above. We also performed an additional simulation with
very high resolution for the $\left(8,0.8,150\right)$ configuration. These
resolutions correspond to a number of grid points $N=\\{80,96,120,144\\}$ in
the boxes surrounding the punctures. Typically the width of the smallest box
around each black hole is on the order of $\sim 2m/N$, where $m$ is the mass
of that black hole; the details of how the grid is determined for each
configuration are given in Sec. II.2.2. We extracted the waveform data at
$R_{\text{ext}}=\left\\{50,60,70,80,90\right\\}M$, which were all on the same
refinement level.

In quantifying the error in the waveforms due to these two sources we focus on
estimating the mismatch between the medium resolution waveforms extracted at a
distance of 90M from the source and the “true” waveform at infinitely good
resolution and infinitely far from the source. We calculated the convergence
order of the BAM code with respect to the resolution and extraction radius
then used this to extrapolate the mismatch. We also used Richardson
extrapolation to estimate the truncation error due to resolution and
extraction radius.

Mismatches are calculated from a noise-weighted inner product between
waveforms, and extrapolate differently to the quantities that are usually
considered in an convergence analysis, e.g., waveform amplitude and phase. In
Sec. IV.3.1 we sketch out how standard numerical convergence properties
translate to the waveform mismatch, and provide more detailed derivations in
Appendix C.

Figure 4: Comparison of coprecessing frame time domain amplitudes for the
$l=m$ modes for $l\in\\{2,3,4,5\\}$. The left panel shows CF_7 with initial
parameters $(q,\chi_{2},\theta)=(1,0.4,60)$ and the right panel shows CF_79
with parameters $(q,\chi_{2},\theta)=(8,0.8,120)$.

###  IV.1 Data quality

Our accuracy analysis of the gravitational waveforms calculated from our
simulations focuses on the mismatch uncertainty as detailed in section IV.3.
This is because it is the overall mismatch uncertainty that is most relevant
to most gravitational-wave astronomy applications. We also consider the error
in the phase and amplitude of the dominant multipole of the co-precessing
waveform in section IV.2. However, we often also wish to know the accuracy of
the individual signal multipoles, for example when using them to construct
waveform models, or when using the NR waveforms as proxy signals to test
gravitational-wave data-analysis pipelines. In this paper we do not perform a
separate convergence analysis of the individual multipoles; given that clean
convergence is rare in any binary-black-hole waveforms, even for the dominant
multipoles, we do not expect a convergence analysis of sub-dominant multipoles
to be informative.

Here we simply note that the phasing accuracy of the waveforms is dominated by
the phase accuracy of the inspiral dynamics, and this can be assessed through
an accuracy analysis of the dominant multipole. (An important exception is the
signal near merger, as discussed in Ref. Calderón Bustillo et al. (2015).) For
the signal amplitude we assess the accuracy by the presence of noise in the
data. For example, Fig. 4 compares the relative strength of the gravitational-
wave multipoles for two simulations. We see that it is not possible to
conclude that a particular set of multipoles will always be reliable. In the
CF_7 simulation the $(3,3)$ and $(5,5)$ multipoles cannot be trusted before
merger; we would not expect these to be useful, for example, to calibrate a
model of the signal amplitude. On the other hand, in the CF_79 simulation we
see that, despite a low level of noise at early times, all of the multipoles
in the figure could well be used to model the amplitude. Rather than choose a
set of “trustworthy” multipoles, we instead suggest that for most applications
one should use only the parts of a post-relaxation-time $\psi_{4,\ell m}$
multipole with an amplitude above $10^{-5}$. Depending on the application, of
course, one may wish to apply a more or less stringent requirement.

###  IV.2 Amplitude and Phase accuracy

In order to estimate the numerical error in the waveform quantities due to the
finite resolution of the simulation and the finite radius at which the data
were extracted, we performed Richardson extrapolation; see Appendix C. This
requires first estimating the convergence order of the code with respect to
these quantities. We first processed the data, removing the junk radiation
from the waveform in the inertial frame in which the simulation was performed.
We aligned the waveforms at merger, where merger is defined to be the time at
which the quantity

$\displaystyle\mathcal{A}^{2}={}$
$\displaystyle\sum^{2}_{m=-2}\left|A_{2m}\left(t\right)\right|^{2},$ (9)

where $A_{2m}$ are the amplitudes of the $\ell=2$ multipoles, is maximised and
resampled using a constant time step of 0.1M. We then rotated the waveform
into the co-precessing frame and aligned the co-precessing phases at merger.
The co-precessing frame is one which precesses along with the binary and is
advantageous here as it means we can focus on the error in a single multipole
(the $\left(2,2\right)$ multipole, which is dominant in this frame, rather
than considering the error in each of the $\ell=2$ multipoles (which all have
appreciable power in the inertial frame) independently. The quantities for
which we are interested in quantifying the numerical error are the amplitude
and phase of the (2,2) multipole in the co-precessing frame as well as two of
the Euler angles $\alpha$ and $\beta$ required to rotate the waveform from the
inertial frame into the co-precessing frame. The Euler angles were calculated
using the method detailed in Ref. O’Shaughnessy et al. (2011); Boyle et al.
(2011).

The standard way to perform a convergence test with respect to the resolution
is to perform a set of three simulations where the resolution improves by a
factor of two between each of the simulations. This is computationally
prohibitive — the high resolution run would use $2^{6}$ times as much memory
as the low resolution run. Similarly, with BAM’s box-based mesh refinement we
cannot extract a waveform at three different radii on the same level a
reasonable distance from the source if each of the extraction radii is twice
as far away from the source as the previous one.

We nominally expect the error due to extraction radius $R_{\rm ext}$ to fall
off as $1/R_{\rm ext}$, although we will confirm that in our analysis. The
numerical-resolution convergence order is less clear. The spatial finite-
differencing in the bulk is sixth-order, but the time-evolution is fourth-
order; either may dominate the error budget, depending on the resolution
choices and length of the simulation Bruegmann et al. (2008); Husa et al.
(2008a). For both the extraction-radius and numerical-resolution, we determine
the appropriate convergence order by studying the convergence behaviour of the
phase of the $(\ell=2,m=2)$ multipole in the co-precessing frame. We then
identify the value of $n$ for which the quantity,

$\displaystyle\delta={}$
$\displaystyle\left(\phi\left(\Delta_{1}\right)-\phi\left(\Delta_{2}\right)\right)-\mathcal{C}\left(\phi\left(\Delta_{2}\right)-\phi\left(\Delta_{3}\right)\right),$
(10)

is minimised, where where $\phi$ is the phase of the $(2,2)$ multipole in the
co-precessing frame,
$\mathcal{C}=\frac{\Delta_{1}^{n}-\Delta_{2}^{n}}{\Delta_{2}^{n}-\Delta_{3}^{n}}$
as in Eq. (32), and $\Delta_{i}$ is the variable in the error expansion, i.e.,
numerical resolution or the inverse extraction radius. The quantity
$\delta(t)$ was calculated over the length of the waveform up to merger and
the mean value $\delta=\langle\delta(t)\rangle$ is shown in Fig. 5. This was
done for both waveforms of varying resolution and extraction radius for the
$q=8$, $\chi=0.8$, $\theta_{\text{LS}}=150^{\circ}$ configuration. In
calculating the convergence order with respect to varying resolution we used
waveforms with $\Delta_{\\{1,2,3\\}}=\\{1/144,1/120,1/96\\}$ while when
considering the convergence order with respect to extraction radius we used
$\Delta_{\\{1,2,3\\}}=\\{1/90,1/70,1/60\\}$. From the results shown in Fig. 5,
we make the conservative conclusion that the code is consistent with fourth
order finite-differencing, which implies that the time-stepping dominates the
error budget. From inspecting time-dependent $\delta(t)$ calculated over the
length of the waveform, we also find that the spatial differencing (with
sixth-order accuracy) dominates the error over the first $\sim 1000M$ of the
waveforms, but the fourth-order-accurate time stepping dominates in the last
$\sim 500M$ before merger, and dominates overall. As expected, we see that the
radiation extraction errors fall off as $1/R_{\text{ext}}$.

Figure 5: The value $\delta$ as given by equation 10 as a function of convergence order. The circle markers represent waveforms of differing resolution. The square markers represent waveforms at differing extraction radii. The configuration $q=8$, $\chi=0.8$, $\theta_{\text{LS}}=150^{\circ}$ was used in this analysis. | % Error
---|---
| Resolution | Extraction Radius | Total
$\phi$ | 0.08 | 0.4 | 0.4
$A$ | 2.5 | 0.9 | 2.7
$\alpha$ | 0.4 | 0.05 | 0.4
$\beta$ | 0.2 | 0.06 | 0.2
Table 3: Relative error in the waveform quantities compared with the
Richardson extrapolated quantities for the $q=8$, $\chi=0.8$,
$\theta_{\text{LS}}=150^{\circ}$ configuration. Figure 6: Resolution
dependence of the absolute error in the time domain co-precessing phase,
relative to the Richardson-extrapolated phase. The phases have been aligned at
merger. Figure 7: Extraction radius dependence of the absolute error in the
time domain co-precessing phase, relative to the Richardson-extrapolated
phase. The phases have been aligned at merger.

Assuming these convergence orders, we then calculate the Richardson-
extrapolated values of the amplitude, phase, $\alpha$ and $\beta$, as
functions of time, using Eq. (31) in App. C. We used the resolutions
$N=\\{120,144\\}$ and the extraction radii $R_{\text{ext}}=\\{80,90\\}$ in
calculating these Richardson-extrapolated values, which were used to estimate
the error in each of these quantities. The error in the waveform quantities
does not monotonically increase with decreasing resolution since not all of
the resolutions lie in the convergence regime. The error in the low and medium
resolution waveforms is therefore estimated from the combination of the error
between these waveforms and the very high resolution waveform and between the
very high resolution waveform and the “true” waveform.

Since a time shift was performed to align the waveforms at merger, where the
phases were then aligned, rather than aligning the frequencies at merger, the
phase difference does not show a quadratic fall off to zero but rather tends
to a constant value and then falls rapidly at merger. As can be seen from Fig.
6, the dephasing of the waveform due to the finite resolution is $\sim 0.1$
radians for the medium resolution ($N=96$) simulation. Similarly, from Fig. 7,
the dephasing due to the finite extraction radius is $\sim 0.4$ radians for
the waveform extracted at $90M$. Since the simulations comprising the bulk of
the catalogue were performed with medium resolution and we recommend using the
waveform extracted at $90M$, these are the key values to focus on. The total
phase error in the waveform is therefore estimated to be about 0.4 radians by
combining the errors in quadrature.

The relative error for each of the quantities we are interested in are given
in Tab. 3. The quantities presented in this table are calculated as follows;
the relative error is taken to be the maximum error, found from Richardson
extrapolation, as described above, divided by the maximum value of the
quantity, over the length of the waveform. Since we aligned the phases at
merger, both the error in the phase and the phase itself are maximum at the
start of the waveform. In contrast, both the amplitude and the error in the
amplitude peak at merger. We therefore report here the relative error in the
peak of the (2,2) multipole in the co-precessing frame. This gives an error
around an order of magnitude larger than during the inspiral, where we see a
total relative error of 0.1%. The order of the relative error in the
precession angles is fairly consistent over the length of the waveform.

From the values given in Tab. 3, we can see that the maximum relative error in
the amplitude of the co-precessing waveform is of the order of a few percent,
while the relative error in the phase and in the precession angles is around
half a percent. This is relevant for the production of a tuned precessing
model using data from these simulations since it implies that the model for
the precession angles cannot be accurate to more than 0.5%. Similar results
were seen for the other simulations for which we have multiple resolutions.

The errors in the amplitude and the precession angles are affected by the
dephasing in the waveform. Therefore, although these results are a good
diagnostic for the reliability of the code and a good way to compare accuracy
between different simulations performed with the same code, they are difficult
to translate into meaningful measures of the accuracy for waveform modelling
or other gravitational wave applications. In order to get a more meaningful
estimate of the accuracy of the waveform we performed the mismatch analysis
presented in the following section.

###  IV.3 Matches

The waveform quantities examined in the previous section are the standard
quantities used when estimating the convergence order and accuracy of a NR
code. While useful when comparing the accuracy between simulations and codes,
these accuracy measures are difficult to interpret in gravitational-wave
astronomy applications — the sensitivity of a search, or the accuracy of a
measurement of the properties of a binary system. When assessing the accuracy
of a waveform it is usually more useful to consider an estimate of the
mismatch error.

The match between two waveforms is defined to be the standard inner product
weighted by the power spectral density of the detector $S_{n}\left(f\right)$
optimised over various sets of parameters $\Theta$ Cutler and Flanagan (1994),

$\displaystyle M\left(h_{1},h_{2}\right)={}$
$\displaystyle\max_{\Theta}\left[4\text{Re}\int^{f_{\text{max}}}_{f_{\text{min}}}\frac{\tilde{h}_{1}\left(f\right)\tilde{h}^{*}_{2}\left(f\right)}{S_{n}\left(f\right)}\text{d}f\right],$
(11)

where the individual waveforms have been normalised so that $M(h,h)=1$. We
also define the mismatch:

$\displaystyle\mathcal{M}={}$ $\displaystyle 1-M\left(h_{1},h_{2}\right).$
(12)

Since these are precessing configurations, we calculate precessing matches as
described in Appendix B of Ref. Schmidt et al. (2015). In order to see how the
match varies over a range of total masses that might be observed by current
ground-based detectors, we further calculate the power-weighted match as
described in Appendix E, based on the work in Ref. Ohme et al. (2011), using
PhenomPv3 Khan et al. (2019) as the model for the low frequency part of the
waveform. We then calculate the mismatch as given by Eq. (12).

We first extrapolate the mismatch due to finite resolution and extraction
radius separately, assuming a particular fall-off in the respective errors. In
order to then find the overall mismatch due to _both_ finite resolution and
extraction radius– i.e. the mismatch with the infinitely far away, infinitely
well resolved “true” waveform– we need to correctly combine these
calculations. The motivation for the correct way to combine such errors is
sketched out in the following section and given in more detail in Appendix D.
We then perform an independent calculation to confirm the mismatches follow
the same convergence relation as the waveform quantities discussed in the
previous section and finally calculate the extrapolated mismatch.

####  IV.3.1 Dependence on expansion parameter and addition of mismatch

In the following we look at how the mismatch behaves with respect to an
expansion parameter, e.g., the numerical resolution or the radius at which the
gravitational-wave signal is extracted. We then consider the addition of
mismatch errors. The calculation below, where the mismatch is expanded in
terms of either the amplitude or phase, is a standard calculation, but we
discuss it in detail here to help motivate the final result, which is somewhat
surprising: although contributions to the error in the amplitude or phase of
the signal combine in quadrature, as one might expect (see, for example, Sec.
II.A of Ref. Lindblom et al. (2008)), separate mismatches should added
according to Eq. (22) below.

To find how the ratio of two matches between waveforms of differing expansion
parameter depends on the expansion parameter, we can examine how the match
depends on the amplitude and phase of the waveform. From Eq. (11), the match
goes as

$\displaystyle M\sim{}$ $\displaystyle\text{Re}\left[\frac{1}{N_{1}N_{2}}\int
h_{1}\left(f\right)h_{2}^{*}\left(f\right)\text{d}f\right],$ (13)

where $N_{i}$ are the normalisation of each of the waveforms respectively. We
take $h_{1}$ to be the waveform containing either the amplitude or phase error
and $h_{2}$ to be the “true” waveform– i.e.,

$\displaystyle h_{1}\left(f\right)={}$ $\displaystyle
h_{2}\left(f\right)+\Delta h\left(f\right),$ (14)

where $h_{2}\left(f\right)=A\left(f\right)e^{i\phi\left(f\right)}$, and $A(f)$
is the real amplitude and $\phi(f)$ is the phase. We assume the true waveform
to be normalised, so

$\displaystyle N_{2}={}$
$\displaystyle\left[\int\left|h_{2}\left(f\right)\right|^{2}\text{d}f\right]^{\frac{1}{2}}=\left[\int
A^{2}\left(f\right)\text{d}f\right]^{\frac{1}{2}}=1.$ (15)

A waveform containing some amplitude error $\Delta A$ is given by

$\displaystyle h_{1}\left(f\right)={}$
$\displaystyle\left(A\left(f\right)+\Delta
A\left(f\right)\right)e^{i\phi\left(f\right)}.$ (16)

Substituting Eq. (16) into (13) we find

$\displaystyle M\sim{}$ $\displaystyle\left[\int\left(A+\Delta
A\right)^{2}\text{d}f\right]^{-\frac{1}{2}}\text{Re}\int A\left(A+\Delta
A\right)\text{d}f$ $\displaystyle={}$
$\displaystyle\left(1+2b+c\right)^{-\frac{1}{2}}\left(1+b\right)$
$\displaystyle\simeq{}$ $\displaystyle 1+\frac{1}{2}\left(b^{2}-c\right),$
(17)

where $b=\int A\Delta A\text{d}f$, $c=\int\left(\Delta A\right)^{2}\text{d}f$
and we have assumed that $\Delta A$ is small in order to make the
approximation in the final step. The mismatch, as defined in Eq. (12),
therefore goes as $b^{2}+c$ and so is dominated by the square of the amplitude
error.

Similarly for a normalised waveform that contains some phase error
$\Delta\phi$,

$\displaystyle h_{1}\left(f\right)={}$ $\displaystyle
A\left(f\right)e^{i\left(\phi\left(f\right)+\Delta\phi\left(f\right)\right)},$
(18)

where $\phi$ is the “true” phase. Substituting this expression into Eq. (13)
we find

$\displaystyle M\sim{}$ $\displaystyle\text{Re}\left[\int
A^{2}e^{i\Delta\phi}\text{d}f\right]$ $\displaystyle\simeq{}$
$\displaystyle\text{Re}\int
A^{2}\left(1+i\Delta\phi-\frac{1}{2}\left(\Delta\phi\right)^{2}\right)\text{d}f$
$\displaystyle={}$ $\displaystyle\int
A^{2}\left(1-\frac{1}{2}\left(\Delta\phi\right)^{2}\right)\text{d}f$
$\displaystyle={}$ $\displaystyle 1-\frac{1}{2}\int
A^{2}\left(\Delta\phi\right)^{2}\text{d}f,$ (19)

where again we have assumed that the error in the phase is small in order to
perform the expansion in the middle step. The mismatch is therefore dominated
by the square of the phase error.

Figure 8: Mismatch between waveforms at varying resolution against the high
resolution ($N=120$) waveform as a function of total mass.  Figure 9: Mismatch
between waveforms at varying extraction radii and the waveform extracted at
$R_{\text{ext}}=90$M as a function of total mass.

The waveform quantities at finite resolution or extraction radius can be
expressed as a Richardson extrapolation of the appropriate expansion parameter
(see Appendix 1). The difference in the phase and amplitude between two
waveforms (labelled $A$ and $B$) is therefore equal to the difference between
the leading order error term (i.e. $\Delta
q\sim\Delta_{B}^{n}-\Delta_{A}^{n}$, see Appendix 2 for more detail). Since
the mismatch is proportional to the _square_ of the error in these waveform
quantities, we find that the convergence relation for the mismatch takes the
form

$\displaystyle\mathcal{M}\left(\Delta_{A}:\Delta_{B}\right)={}$
$\displaystyle\kappa\left(\Delta_{A}^{n}-\Delta_{B}^{n}\right)^{2},$ (20)

where $\Delta_{i}$ is the value of the expansion parameter for the $i$th
waveform and $\kappa$ is a co-efficient to be found. If one of the waveforms
being considered is the “true” waveform and thus contains no numerical error,
then the mismatch between any reference waveform and this true waveform will
be given by
$\mathcal{M}\left(\Delta_{\text{ref}}:\Delta_{\text{true}}=0\right)=\kappa\Delta_{\text{ref}}^{2n}$.
A similar derivation to the one discussed here is also presented in Ref.
Ferguson et al. (2021). From this we can see that the ratio of the mismatch
$\mathcal{M}$ between two pairs of waveforms is given by

$\displaystyle\frac{\mathcal{M}\left(A:B\right)}{\mathcal{M}\left(B:C\right)}={}$
$\displaystyle\frac{\left(\Delta_{A}^{n}-\Delta_{B}^{n}\right)^{2}}{\left(\Delta_{B}^{n}-\Delta_{C}^{n}\right)^{2}}.$
(21)

This result will be used in the following sections to study the convergence
properties of our numerical-relativity waveforms via their mismatch error.

We can also see that the correct way to combine the mismatches between two
sets of waveforms $\mathcal{M}\left(A:B\right)$ and
$\mathcal{M}\left(B:C\right)$ in order to get the mismatch between the final
pair $\mathcal{M}\left(A:C\right)$ is given by

$\displaystyle\mathcal{M}\left(A:C\right)={}$
$\displaystyle\left(\sqrt{\mathcal{M}\left(A:B\right)}+\sqrt{\mathcal{M}\left(B:C\right)}\right)^{2}.$
(22)

A more rigorous proof of this result is presented in Appendix D.

As stated above, our main use for this result is to combine the mismatch due
to different sources of error in our numerical waveforms. We estimate the
mismatch between our waveforms at finite extraction radius and finite
resolution and the true waveform using

$\displaystyle\mathcal{M}={}$
$\displaystyle\left(\sqrt{\mathcal{M}_{\text{resolution}}}+\sqrt{\mathcal{M}_{\text{extraction
radius}}}\right)^{2},$ (23)

where $\mathcal{M}_{\text{resolution}}$ is the mismatch due to the finite
resolution of the numerical simulation and $\mathcal{M}_{\text{extraction
radius}}$ is the mismatch due to the finite distance from the source at which
the waveforms were extracted.

####  IV.3.2 Convergence order

We performed matches between waveforms extracted at $R_{\text{ext}}=90M$ for
the high resolution simulations against all other resolutions available for a
given configuration. These results are shown in Fig. 8. We also performed
matches between waveforms extracted at $R_{\text{ext}}=90$M and all other
available extraction radii for the medium resolution simulations for each
configuration. These results are shown in Fig. 9. In both of these
comparisons, we have calculated the match against a single resolution or
extraction radius. We therefore expect that the matches will improve for cases
where the values of the resolution or extraction radius are closer to each
other. From Fig. 9 we can see that the matches generally follow this trend,
implying that it is reasonable to assume the waveform is being extracted
sufficiently far from the source that we may be in the convergence regime.
This is not true for the mismatches with respect to resolution shown in Fig.
8. The matches between (i) the low and high resolutions and (ii) the medium
and high resolutions clearly do not follow any trend for most of the
configurations. From this we can see that it is not reasonable to treat the
low and medium resolutions as if they lie within the convergence regime.

That the medium resolution does not lie quite within the convergence regime is
demonstrated clearly in Fig. 10, where we show the mismatch between the medium
and high and the high and very high resolutions using Eq. (21) for varying
convergence order. From this analysis it is clear that the mismatch is closest
to being fourth-order convergent. This analysis could only be done for the
case ($q=8$, $\chi=0.8$, $\theta_{\text{LS}}=150^{\circ}$) since this is the
only case for which we have the very high resolution run.

Figure 10: Mismatches demonstrating fourth order convergence of the BAM code
with respect to resolution. The solid lines show the calculated mismatch while
the dashed lines show the predicted mismatch for varying convergence order.
The mismatch was calculated with respect to the $N=120$ resolution run. This
is for the $q=8$, $\chi=0.8$, $\theta_{\text{LS}}=150^{\circ}$ configuration.

Conversely, since it seems reasonable to assume the waveforms extracted at
varying extraction radii mostly lie within the convergence regime, we
calculated the ratio of the mismatch between each of the pairs of waveforms
from different extraction radii using Eq. (21) for varying convergence order.
For each of the four configurations we investigated it was found that the
results were most consistent with first order convergence. This is
demonstrated in Fig. 11, where the solid lines show the calculated mismatch
between two waveforms of different extraction radii and the dotted red line
shows the expected value of the match for first order convergence.

Not all the waveforms from the different extraction radii show perfect
convergence for every configuration. The mismatch between $R_{\text{ext}}=80$
and $R_{\text{ext}}=90$ often does not follow the trend– we expect this is
because the mismatch between these waveforms is so small
$\left(\mathcal{O}\left(10^{-6}\right)\right)$ that it is very sensitive to
any data processing performed in the course of calculating the match. The
mismatch between $R_{\text{ext}}=50$ and $R_{\text{ext}}=90$ also often does
not follow the trend and we do not expect it to hold for small extraction
radii.

The convergence order calculated using this method agrees with the estimate
calculated in the previous section; the code is approximately fourth order
convergent with respect to resolution and first order convergent with respect
to the extraction radius.

####  IV.3.3 Extrapolation

Figure 11: Mismatches demonstrating first order convergence of the BAM code
with respect to extraction radius. The black lines indicated in the legend
show the calculated mismatch, while the dotted red line shows the predicted
mismatch for the pair of waveforms indicated by the solid line based on the
mismatch indicated by the dashed line, assuming first order convergence.
Figure 12: Variation of the mismatch with extraction radius for a system with
a total mass of 100M. The mismatches shown are with respect to the waveform
extracted at $R_{\text{ext}}=90$M. The line shows the fit based on Eq. (10)
assuming first order convergence. It is consistent with all extraction radii
except at $R_{\text{ext}}=50$M. Figure 13: Variation of the mismatch with
resolution for a $q=8$, $\chi=0.8$, $\theta_{\text{LS}}=150^{\circ}$ system
with a total mass of 100M. The mismatches shown are with respect to the medium
resolution run and so the mismatch is zero at $N=120$. The line shows the
relation in Eq. (20) with $\kappa_{\text{res}}$ assuming fourth order
convergence calculated using Eq. (24) and is consistent with all resolutions
except the lowest one at $N=80$.  Figure 14: Projected mismatch between a
waveform extracted at a resolution of $N=96$ and one that is infinitely well
resolved.

Having verified the convergence order of the code $n$, we can calculate the
convergence relation of the mismatches shown in Figs. 8 and 9. We first look
at the mismatch due to the finite resolution of the simulation. We assume that
the two highest resolution simulations we performed ($N=120$ and $N=144$) lie
in the convergence regime but we know the two lower resolution simulations do
not. Assuming fourth order convergence, we use Eq. (20) to calculate
$\kappa_{\text{res}}$ using

$\displaystyle\kappa_{\text{res}}={}$
$\displaystyle\frac{\mathcal{M}\left(\Delta_{144}:\Delta_{120}\right)}{\left(144^{4}-120^{4}\right)^{2}}.$
(24)

From Fig. 13, which shows the convergence relation for the mismatches
calculated for a system with total mass 100M⊙, we can see that this appears to
be a reasonable assumption. From $\kappa_{\text{res}}$ we can estimate the
mismatch between the high or very high resolution waveforms with an infinitely
well resolved waveform. However, we want to know the mismatch for the medium
resolution runs since this is the resolution that was used to perform the
simulations for the catalogue of waveforms presented in Tabs. 1 and 2. Since
this resolution does not lie in the convergence regime (and the phase error
does not improve monotonically from the medium resolution to the high and very
high resolution waveforms) we cannot simply use the calculated convergence
relation in order to estimate the mismatch between a waveform at this
resolution and the “true” waveform. Instead we use the formula given in Eq.
(22) to add the mismatch between the medium resolution and the very high
resolution waveforms to the mismatch between the very high resolution waveform
and the “true” waveform:

$\displaystyle\mathcal{M}\left(\Delta_{96}:\Delta_{\infty}\right)={}$
$\displaystyle\left(\sqrt{\mathcal{M}\left(\Delta_{96}:\Delta_{144}\right)}+\right.$
$\displaystyle\left.\qquad\sqrt{\mathcal{M}\left(\Delta_{144}:\Delta_{\infty}\right)}\right)^{2}$
$\displaystyle={}$
$\displaystyle\left(\sqrt{\mathcal{M}\left(\Delta_{96}:\Delta_{144}\right)}+\sqrt{\frac{\kappa_{\text{res}}}{144^{4}}}\right)^{2}.$
(25)

The result of this extrapolation procedure is shown in Fig. 14. We could only
perform this calculation for the case ($q=8$, $\chi=0.8$,
$\theta_{\text{LS}}=150^{\circ}$) since this is the only case for which we
have the very high resolution run. However, from Fig. 8 we can see that the
mismatch between the medium and high resolution runs is the worst for this
case, so this estimate should give an upper bound for the mismatch between the
medium resolution run and the “true” waveform. Fig. 14 shows the projected
mismatch between a medium resolution waveform and one that is infinitely well
resolved for a range of total masses. The maximum mismatch between a medium
resolution waveform and an infinitely well resolved one is $6.0\times
10^{-4}$.

Figure 15: Projected mismatch between a waveform extracted at a radius of
$R_{\text{ext}}=90$M and one extracted infinitely far away.

We next examine the mismatch due to the finite distance from the source at
which the waveform is extracted. To calculate the first order convergence
relation with respect to the extraction radius, we performed a fit through
each of the mismatches which were found to follow the convergence relation.
This is demonstrated for mismatches between waveforms of different extraction
radii and the waveform at $R_{\text{ext}}=90$M in Fig. 12, for a system with
total mass 100M⊙. This fit gives the value of $\kappa_{\text{ext}}$ for every
value of the total mass of the system. From this we can calculate the mismatch
between the waveform at $R\text{ext}=90$M and the “true” waveform from
$\mathcal{M}\left(\Delta_{90}:\Delta_{\infty}\right)=\frac{\kappa_{\text{ext}}}{90^{2}}$.

The mismatch between the waveform extracted at $R_{\text{ext}}=90$M and the
“true” waveform is shown in Fig. 15. The configuration that gives the greatest
mismatch is $q=4$, $\chi=0.4$, $\theta_{\text{LS}}=60$ since, as noted above,
this simulation was much longer than the others and so has greater opportunity
to accumulate phase error between the two waveforms. The maximum mismatch
between a waveform at $R_{\text{ext}}=90$ and at
$R_{\text{ext}}\rightarrow\infty$ is taken to be $1.4\times 10^{-3}$.

Figure 16: Comparison between the parameters of the new BBH simulations
presented here (CF) and the existing BBH simulations in the the SXS, RIT and
Maya catalogues. The _Top left_ spin disk shows simulations with spin on the
larger black hole $0<\chi<0.25$, _Top right_ $0.25\leq\chi<0.5$, _Bottom left_
$0.5\leq\chi<0.75$ and _Bottom right_ $\chi\geq 0.75$. The radius of each disk
shows the mass ratio of the binary and the orientation shows the spin tilt
angle of the larger black hole. Spin tilt angles of $90^{\circ}$ means that
the spin vector lies in the binary’s orbital plane.

We estimate the mismatch between our medium resolution waveform extracted at
$R_{\text{ext}}=90$M and the true waveform using Eq. (23). A conservative
estimate of the mismatch between a waveform extracted at a finite distance of
90M from the source for a simulation performed with a grid spacing $d=0.0104$
and the theoretical ‘analytical’ solution is therefore $3.83\times 10^{-3}$.
This provides a limit on the mismatch error of the waveforms presented in this
catalogue of 0.004, or 0.4%.

##  V Catalogue comparison

A number of numerical relativity groups have started building larger and more
comprehensive catalogues that span a growing region of the parameter space. At
time of publication, there are a total of 4,352 publicly available BBH
simulations in a combination of the Simulating eXtreme Spacetimes (SXS)
Collaboration Mroue et al. (2013); Boyle et al. (2019) and the Rochester
Institute of Technology (RIT) Healy et al. (2017, 2019); Healy and Lousto
(2020, 2022) and Maya catalogues Jani et al. (2016).

The SXS collaboration has produced the largest catalogue to date with 2,019
BBH simulations spanning $1\leq q\leq 10$ and $0\leq\chi\leq 1$. The RIT
catalogue contains 1,881 BBH simulations covering $1\leq q\leq 128$ and
$0\leq\chi\leq 0.99$ and the Maya catalogue contains 452 unique BBH waveforms
from more than 600 BBH simulations ranging between $1\leq q\leq 15$ and
$0\leq\chi\leq 0.8$. Unlike the simulations presented here, the SXS, RIT and
Maya catalogues all contain simulations where the individual black hole spins
can be zero or perfectly aligned/anti-aligned with the orbital angular
momentum. Considering only the precessing parameter space, the SXS
collaboration has produced 1,429 simulations spanning $1\leq q\leq 6$ and
$0\leq\chi\leq 0.99$, the RIT catalogue contains $561$ simulations covering
$1\leq q\leq 15$ and $0\leq\chi\leq 0.99$ and the Maya catalogue contains 315
waveforms ranging between $1\leq q\leq 8$ and $0.1\leq\chi\leq 0.8$.

Fig. 16 compares the parameter space coverage of the simulations presented
here with the existing non-zero-spin simulations included in the SXS, RIT and
Maya catalogues over the mass ratio and larger black hole spin tilt angle and
spin magnitude parameter space. We see that although the existing catalogues
provide good coverage for high black hole spins $\chi\geq 0.5$ and near equal
mass ratios, there is a dearth of existing precessing simulations for low
black hole spins $\chi<0.25$ and unequal mass ratios $q\gtrsim 4$.

Recent gravitational wave observations (see e.g. Abbott et al., 2020b, c) have
shown a need for BBH simulations in this low black hole spin and unequal mass
ratio region of the parameter space in order to build reliable and accurate
waveform models for use in Bayesian inference. The uniform coverage of the
single-spin space up the $q=8$ has made it possible to construct an accurate
generic precessing-binary model for future observations Hamilton et al.
(2021). Most astrophysical models suggest that BBH at larger mass ratios will
be rare (e.g., Ref. Abbott et al. (2021g)), but given that there has been one
observation to date at $q\sim 10$ Abbott et al. (2020c) extension of this
parameter-space coverage to yet higher mass ratios will be necessary in the
future.

##  VI Discussion

We have produced a catalogue containing 80 waveforms from single-spin
precessing systems with mass ratios up to $q=8$, dimensionless spin magnitudes
up to $\chi=0.8$ and a range of spin inclination angles. In all cases the spin
was placed on the larger black hole. We estimate our uncertainty in the masses
of the initial black holes to be $\mathcal{O}(0.05\%)$. We estimate the
uncertainty in the initial spin magnitude to be $\mathcal{O}(10^{-3})$ while
the uncertainty in the initial spin inclination is $\mathcal{O}(1^{\circ})$.
Similarly, we obtain estimates of the uncertainty of the remnant properties
reported in this paper. We find the final mass has an uncertainty of $5\times
10^{-4}$, while the final spin magnitude is accurate to within $5\times
10^{-3}$.

The starting frequency of the simulations was chosen such that the simulations
were all around a similar length ($\sim 2000M$) in order to limit the
dephasing in the waveform and thus ensure sufficient accuracy throughout the
evolution of the binary. We performed a careful analysis of the errors due to
the finite resolution of the simulations and due to the finite distance from
the source at which the waveforms were extracted. From this we were able to
provide a conservative estimate of the mismatch uncertainty of our waveforms
of 0.4%.

The catalogue presented here is sufficient to capture a wide range of single-
spin precession effects. Most notably, the systems contained within it have a
non-zero opening angle of the precession cone ranging from $\sim 1^{\circ}$ to
$\sim 115^{\circ}$ radians. The cases with the largest opening angles display
the most dominant precession effects. In particular, for initial
configurations with mass ratios $q\sim 8$, high spins and large spin
inclination angles, the final spin will be in the opposite direction to the
binary’s angular momentum prior to merger, thus producing a “negative” final
spin. The majority of precessing simulations in other catalogues do not extend
beyond $q=5$ so consequently, this region of parameter space is poorly covered
by NR simulations. Indeed, in this catalogue, despite having 20 precessing
simulations at $q=8$ we see only two cases where we have a negative final spin
(CF_75 and CF_80). The phenomenology of this region has therefore not yet been
thoroughly explored and a more detailed study is planned for future work.
Finally, we also see a wide range of recoil velocities for the configurations
included in this catalogue, with the highest values seen for equal mass
systems. We also see the greatest range of values for equal mass systems,
depending on the in-plane spin angle. For $q=8$ systems we see much lower
values in general across all cases.

While this catalogue was sufficient to produce the first inspiral-merger-
ringdown (IMR) model of precessing systems tuned to NR, PhenomPNR, it will
need to be greatly expanded in order to meet modelling requirements of future
gravitational wave observations. Existing catalogues (such as the SXS, Maya
and RIT catalogues) provide a comprehensive coverage of the two-spin
precessing parameter space up to $q=4$. This catalogue provides a systematic
coverage of the single-spin precessing parameter space up to $q=8$. However,
while it uses a consistent in-plane spin direction at the starting frequency
(the initial configurations all have the in-plane spin component along the
binary’s separation vector), these will translate into quite different spin
directions at merger. Consequently, for any modelling that includes effects
due to the in-plane spin direction, this catalogue contains an incomplete and
possibly random sampling of points.

There are many directions in which this catalogue can be expanded to: include
higher mass ratios, comprehensively cover rotations of the in-plane spin
component, include two-spin systems, produce longer and more accurate
waveforms and include binaries on eccentric orbits. Since the production of
these simulations are expensive (the catalogue presented here is estimated to
have required around 25 million CPU hours in total) and generic modelling at
higher mass ratios and for longer waveforms is not a completely solved
problem, it is an open question as to which direction in parameter space is
most urgent.

Assuming the mass ratio distribution reported in Abbott et al. (2021g), we
estimate that only 1.3% of observed binaries will have $q>8$. This is
supported by gravitational wave detections so far since, out of the 90
binaries reported by the LVK collaborations Abbott et al. (2021b), only one
has been found to have a mass ratio clearly greater than 8 Abbott et al.
(2020c). Similarly, from the production of PhenomPNR, we know that we will
require longer waveforms for binaries with higher mass ratios and spins,
particularly those with a spin inclination angle of
$\theta_{\text{LS}}>90^{\circ}$, since inaccuracies in post-Newtonian (PN)
estimates of precession effects become more appreciable towards lower
frequencies in this region of the parameter space. For the existing catalogue,
the cases CF_79 and CF_80 (at $q=8$ and $\chi=0.8$) are already sufficiently
short to limit model accuracy and we anticipate this will continue for
decreasing spin magnitude as we go to higher mass ratios Hamilton et al.
(2021). This is therefore a smaller fraction of binaries than the simple
requirement to extend to high mass ratios. Similarly, we expect to see two-
spin effects in signals with SNRs greater than 100 Pürrer et al. (2016); Khan
et al. (2020). Taking the detection threshold to be SNR 10, then we expect to
be able to identify two-spin effects in 0.1% of detections.

We have not considered eccentricity here, but eccentric modelling and NR
simulations have been discussed in other works Boyle et al. (2019); Ramos-
Buades et al. (2019, 2020); Huerta et al. (2018); Hinder et al. (2018); Nagar
et al. (2021); Islam et al. (2021); Liu et al. (2022).

From this, we conclude that the most urgent extension is required to systems
with higher mass ratios. This is closely followed by a systematic sampling
that will explore the most dominant physical effects of two-spin systems (such
as those that will impact the opening angle of the precession cone at merger)
or the in-plane direction (such as the recoil velocity).

##  VII Acknowledgements

We would like to thank Steve Fairhurst, Frank Ohme, Vivien Raymond for many
useful discussions; Kieran Philips for extensive optimisation of the BAM code
that resulted in a greater than 30% increase in speed and a much lower memory
footprint, and also for his work in testing the code on the Oracle Cloud
Infrastructure (OCI); Paul Hopkins for wide-ranging cluster support and also
assistance in setting up OCI runs; Phil Bates at Oracle for initiating the OCI
work and to Phil and his team for extensive technical support and guidance.

The authors were supported in part by Science and Technology Facilities
Council (STFC) grant ST/V00154X/1 and European Research Council (ERC)
Consolidator Grant 647839. E. Hamilton was supported in part by Swiss National
Science Foundation (SNSF) grant IZCOZ0-189876. L. London was supported at
Massachusetts Institute of Technology (MIT) by National Science Foundation
Grant No. PHY-1707549 as well as support from MIT’s School of Science and
Department of Physics. A. Vano-Vinuales thanks FCT for financial support
through Project No. UIDB/00099/2020.

This work used the DiRAC@Durham facility managed by the Institute for
Computational Cosmology on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University
and STFC operations grant ST/R000832/1. DiRAC is part of the National
e-Infrastructure. Additionally, this research was undertaken using the
supercomputing facilities at Cardiff University operated by Advanced Research
Computing at Cardiff (ARCCA) on behalf of the Cardiff Supercomputing Facility
and the HPC Wales and Supercomputing Wales (SCW) projects. We acknowledge the
support of the latter, which is part-funded by the European Regional
Development Fund (ERDF) via the Welsh Government. This work was also supported
in part by Oracle Cloud credits and related resources provided by the Oracle
for Research program.

##  Appendix A Iterative initial data construction for single spin precessing
systems

The initial data construction method of Husa et al. (2016) for aligned-spin
systems is extended to precessing systems using the follow iterative brute-
force algorithm. A sequence of input parameter sets

$\bm{\theta}_{i}\equiv\left(q_{i}^{D},\bm{S}_{1,i}^{D},\bm{S}_{2,i}^{D}\right)=\left(q,\bm{0},S_{2}\hat{\bm{t}}_{i}\right),$
(26)

will be iteratively refined, defined at a user-specified separation of
$D_{\text{start}}$ and used as initial conditions for the EOB solver. The
single-spin orienation unit vector $\hat{\bm{t}}$ can be expressed in terms of
angles $\theta$ and $\phi$

$\hat{\bm{t}_{i}}=\big{(}\cos(\phi_{i})\cos(\theta_{i}),\sin(\phi_{i})\cos(\theta_{i}),\cos(\theta_{i})\big{)}$
(27)

The EOB simulations are started with both component black holes placed on the
x-axis and the orbital angular momentum parallel to the $z$-axis. For this
work the black holes are placed $\sim 40M$ apart. Each successive set of
parameters are chosen based on the EOB evolution of previous parameter sets in
the sequence. This is repeated until the EOB evolution results in the required
parameters $(q,\bm{S}_{1},\bm{S}_{2})$ at $M\Omega_{\mathrm{orb}}$ using the
following algorithm

  1. (1)

Initial candidate parameters at $D_{\text{start}}$ ($n=0$)
$\phi_{0}$ is chosen to be the target azimuthal spin angle $\phi$. The EOB
solver is then run until $M\Omega_{\mathrm{orb}}$ is reached and the azimuthal
spin angle at that time is recorded as $\phi_{\omega,0}$. The EOB spin
dynamics are explored in the region around $M\Omega_{\mathrm{orb}}$ to find
the closest time when the spin angle is equal to $\phi$. The difference in
frequency $\Delta M\omega_{0}$ between this time and $M\Omega_{\mathrm{orb}}$
is recorded. If $|\Delta M\omega_{0}|\leqslant M\omega_{\mathrm{tol}}$ where
$M\omega_{\mathrm{tol}}$ is a user specified tolerance, then the algorithm
stops. For the initial data generated in this work the percentage error
tolerance of the orbital frequency is specified to be
$M\omega_{\mathrm{tol}}=1\%$. If $|\Delta M\omega_{0}|>M\omega_{\mathrm{tol}}$
then proceed to the next step.

  2. (2)

Second candidate parameters at $D_{\text{start}}$ ($n=1$)
$\phi_{1}$ is chosen to be the difference between the target azimuthal spin
angle and the azimuthal spin angle of the EOB spin dynamics at
$M\omega_{\mathrm{orb}}$, given by $\phi_{c}\equiv\phi-\phi_{\omega,0}$. The
EOB spin dynamics are explored again, recording $\phi_{\omega,1}$ at
$M\omega_{\mathrm{orb}}$, and calculating $|\Delta M\omega_{1}|$ as in the
same way as $|\Delta M\omega_{0}|$ in the previous step. If $|\Delta
M\omega_{1}|\leqslant M\omega_{\mathrm{tol}}$ the algorithm stops otherwise
proceed to the next step.

  3. (3)

Third candidate parameters at $D_{\text{start}}$ ($n=2$)
Set $\Delta\phi$ to be $10^{\circ}$ if $||\Delta
M\omega_{1}|-M\omega_{\mathrm{tol}}|>M\omega_{\mathrm{tol}}/2$ otherwise set
$\Delta\phi$ to be $5^{\circ}$. $\phi_{2}$ is chosen to be the target
azimuthal spin angle $\phi_{c}+\Delta\phi$. $\phi_{\omega,2}$ and $|\Delta
M\omega_{2}|$ are calculated in the same way as previous steps. If $|\Delta
M\omega_{2}|\leqslant M\omega_{\mathrm{tol}}$ the algorithm stops otherwise
proceed to the next step.

  4. (4)

Further candidate parameters at $D_{\text{start}}$ ($n>2$)
Set $\Delta\phi$ to be $10^{\circ}$ if $||\Delta
M\omega_{n-1}|-M\omega_{\mathrm{tol}}|>M\omega_{\mathrm{tol}}/2$ otherwise set
$\Delta\phi$ to be $5^{\circ}$. If $\phi_{\omega,2}>\phi_{\omega,1}$ this
indicates that the azimuthal spin angle is being rotated in the wrong
direction. As such, if $\phi_{\omega,2}>\phi_{\omega,1}$ then
$\phi_{n}=\phi_{c}+(2-n)\Delta\phi$, otherwise set
$\phi_{n}=\phi_{c}+(n-1)\Delta\phi$. $|\Delta M\omega_{n}|$ is calculated in
the same way as previous steps. If $|\Delta M\omega_{n}|\leqslant
M\omega_{\mathrm{tol}}$ the algorithm stops otherwise repeat this step until
this inequality is satisfied.

Once the required tolerance is met and the algorithm stops, the position,
linear momentum and spin of each black hole are taken from the EOB dynamics at
$M\Omega_{\mathrm{orb}}$ and used as input for the Bowen-York initial data
solver.

##  Appendix B NR simulations with cloud computing

For this work we have run individual BAM simulations on up to 512 processors,
and these require fast inter-processor connections to ensure that inter-
processor communication is not the dominant limitation on the calculation
speed. This is typical for NR codes, and as such these are usually run on
clusters that have been constructed primarily for highly parallelised high-
performance computing (HPC), such as the DiRAC Cosma clusters that were used
for the majority of the runs presented here. An increasingly popular source of
computing resources are cloud services. These have historically been set up
with large numbers of independent (high-throughput computing) applications in
mind. However, recently some services have improved the speed of inter-
processor communication, with the goal of making cloud computing services
competitive for HPC.

As part of the NR work presented here, we investigated the performance of BAM
on the Oracle Cloud Infrastructure (OCI). We performed a series of experiments
to determine the optimal performance we could achieve with the hardware
available at the time (2018). These tests used a bespoke “bare metal” setup
and ran at about 60% of the speed on the DiRAC cosma5 cluster. (Note that
since these tests were done, both the cosma clusters and the cores used at OCI
have been superseded by newer hardware.)

We also completed a set of production simulations: these were the five NR
simulations at mass-ratio $q=2$ and black-hole spin $S_{2}/m_{2}=0.6$
(CF_31–CF_35). Each run was performed on 128 cores and required approximately
140,000 CPU hours.

The production simulations used the “cluster-in-the-cloud” infrastructure
Woods and Williams  to create container-based cluster instances using OCI
Terraform on AMD64 128-core BM.Standard.E2.64 nodes, and ran at an average 1.8
$M$/hr. Similar runs on cosma6 ran at about 3.7 $M$/hr on double the number of
cores. These suggest that in terms of computational cost and efficiency,
cloud-based resources have the potential to be competitive to standard
clusters.

##  Appendix C Convergence estimates and Richardson extrapolation

###  1 Richardson Extrapolation

A quantity $q$ calculated at finite resolution or extraction radius can be
given by

$\displaystyle q^{*}={}$ $\displaystyle q\left(\Delta\right)+e_{i}\Delta^{i},$
(28)

where $\Delta$ is the expansion parameter ($\frac{1}{N}$ for resolution or
$\frac{1}{R_{\text{ext}}}$ for extraction radius), $e$ is the finite order
error and $i$ is the order at which the error contributes. In this paper the
quantity $q$ we are considering is the waveform extracted from the numerical
simulation. We therefore have that

$\displaystyle q^{*}={}$ $\displaystyle
q\left(\Delta\right)+e_{n}\Delta^{n}+\mathcal{O}\left(\Delta^{n^{\prime}}\right),$
(29)

where $e_{n}\Delta^{n}$ is the leading order error contribution, $n$ is the
convergence order of the simulation and $n^{\prime}>n$.

Considering two waveforms computed using different expansion parameters
$\Delta_{1}$ and $\Delta_{2}$, we can solve the two simultaneous equations
that arise from Eq. (29) to give

$\displaystyle q^{*}={}$ $\displaystyle
R\left(\Delta_{1},\Delta_{2}\right)+\mathcal{O}\left(\left(\frac{\Delta_{1}}{\Delta_{2}}\right)^{n^{\prime}}\right),$
(30)

where

$\displaystyle R\left(\Delta_{1},\Delta_{2}\right)={}$
$\displaystyle\frac{\left(\frac{\Delta_{1}}{\Delta_{2}}\right)^{n}q\left(\Delta_{2}\right)-q\left(\Delta_{1}\right)}{\left(\frac{\Delta_{1}}{\Delta_{2}}\right)^{n}-1},$
(31)

is the _Richardson extrapolation_ Richardson and Glazebrook (1911); Alcubierre
(2008) of $q\left(\Delta_{1}\right)$. $R\left(\Delta_{1},\Delta_{2}\right)$
has a higher order error due to the truncation of the expansion in $\Delta$
than $q\left(\Delta_{1}\right)$.
$R\left(\Delta_{1},\Delta_{2}\right)-q\left(\Delta_{1}\right)$ gives the
truncation error of the quantity $q$.

###  2 Convergence

Considering now three waveforms computed with expansion parameters
$\Delta_{1}>\Delta_{2}>\Delta_{3}$ we can eliminate $q^{*}$ in Eq. (30).
Neglecting higher order error terms, the ratio of the difference between two
sets of numerical waveforms with expansion parameter
$\Delta_{1}>\Delta_{2}>\Delta_{3}$ is then given by

$\displaystyle\mathcal{C}={}$
$\displaystyle\frac{q\left(\Delta_{1}\right)-q\left(\Delta_{2}\right)}{q\left(\Delta_{2}\right)-q\left(\Delta_{3}\right)}=\frac{\Delta_{1}^{n}-\Delta_{2}^{n}}{\Delta_{2}^{n}-\Delta_{3}^{n}}.$
(32)

This relation holds for features of a waveform, such as its amplitude and
phase, but not for derived quantities such as the match.

To understand how the match between a set of waveforms in a convergence series
varies, consider a detector response derived from gravitational wave strain
solutions of a finite difference approximation numerical relativity code. This
can be represented by a _Richardson expansion_ Richardson and Glazebrook
(1911); Alcubierre (2008) as a power series in an expansion parameter
$\Delta$. Consider two detector responses, $h_{1}$ and $h_{2}$, at two
resolutions with expansion parameters $\Delta_{1}$ and $\Delta_{2}$
respectively. As seen in Eq. (29), for an $i^{\mathrm{th}}$ order accurate
finite difference method these can be represented by their truncated
Richardson expansions

$\displaystyle h_{1}(f)$ $\displaystyle=h(f)+e_{n}(f)\Delta_{1}^{n},$ (33)
$\displaystyle h_{2}(f)$ $\displaystyle=h(f)+e_{n}(f)\Delta_{2}^{n},$ (34)
$\displaystyle=h_{1}(f)+e_{n}(f)(\Delta_{2}^{n}-\Delta_{1}^{n}),$ (35)

where $h(f)$ is the detector response of the exact solution, and $e_{i}(f)$
are the leading order error functions.

The match between these two detector responses can be expanded in the
expansion parameter $(\Delta_{2}^{n}-\Delta_{1}^{n})$. Utilising linearity in
the inner product, we have

$\displaystyle M(h_{1},h_{2})$
$\displaystyle=\max_{\Theta}\left[\frac{\left\langle
h_{1}\middle|h_{2}\right\rangle}{\|h_{1}\|\|h_{2}\|}\right],$
$\displaystyle=\max_{\Theta}\left[\frac{\|h_{1}\|^{2}+(\Delta_{2}^{n}-\Delta_{1}^{n})\left\langle
h_{1}\middle|e_{n}\right\rangle}{\|h_{1}\|\left(\|h_{1}\|^{2}+2(\Delta_{2}^{n}-\Delta_{1}^{n})\left\langle
h_{1}\middle|e_{n}\right\rangle+(\Delta_{2}^{n}-\Delta_{1}^{n})^{2}\|e_{n}\|^{2}\right)^{1/2}}\right],$
$\displaystyle\approx\max_{\Theta}\left[1-\frac{1}{2}\left(\frac{\|e_{n}\|^{2}}{\|h_{1}\|^{2}}-\left(\frac{\left\langle
h_{1}\middle|e_{n}\right\rangle}{\|h_{1}\|^{2}}\right)^{2}\right)(\Delta_{2}^{n}-\Delta_{1}^{n})^{2}\right],$
$\displaystyle\approx 1-\kappa(\Delta_{2}^{n}-\Delta_{1}^{n})^{2},$ (36)

where between the second and third lines we have performed a binomial
expansion of the denominator and terms of higher order in the expansion
parameter are dropped between steps. The constant coefficient $\kappa$ is
defined as

$\kappa\equiv\min_{\Theta}\left[\frac{1}{2}\left(\frac{\|e_{n}\|^{2}}{\|h\|^{2}}-\frac{\left\langle
h\middle|e_{n}\right\rangle^{2}}{\|h\|^{4}}\right)\right],$ (37)

where we have used Eq. (33) to re-write $h_{1}$ in terms of $h$ and again
neglected terms of higher order in the expansion parameter. $\kappa$ can be
seen to be bounded below by zero from the Cauchy-Schwarz inequality. The
mismatch may then be approximated as

$\mathcal{M}(h_{1},h_{2})=\kappa(\Delta_{2}^{n}-\Delta_{1}^{n})^{2}.$ (38)

It is important to note that the leading order expansion parameter terms in
the approximation Eq. (38) are quadratic in the expansion parameters. In
addition, while it is likely that the leading order coefficient $\kappa$
cannot be calculated directly, it is independent of any resolution-specific
expressions. As such $\kappa$ is constant within any convergence series. This
leads to the following two results, describing ratios of mismatches in
convergence series and the combination of mismatches in convergence series,

$\displaystyle\frac{\mathcal{M}(h_{1},h_{2})}{\mathcal{M}(h_{2},h_{3})}$
$\displaystyle=\frac{(\Delta_{1}^{n}-\Delta_{2}^{n})^{2}}{(\Delta_{2}^{n}-\Delta_{3}^{n})^{2}},$
(39) $\displaystyle\mathcal{M}(h_{1},h_{3})$
$\displaystyle=\left(\sqrt{\mathcal{M}(h_{1},h_{2})}+\sqrt{\mathcal{M}(h_{2},h_{3})}\right)^{2}.$
(40)

Eq. (40) holds generally, not just for the case of mismatches between
waveforms in a convergence series. This is shown in Appendix D.

##  Appendix D Addition of mismatches

Consider three waveforms, $h_{1}$, $h_{2}$, $h_{3}$, which are all normalised,
$|h_{1}|=|h_{2}|=|h_{3}|=1$. If $\langle h_{1}|h_{2}\rangle=A$, and $\langle
h_{2}|h_{3}\rangle=B$, we would like to estimate an upper bound on $C=\langle
h_{1}|h_{3}\rangle$.

We write each waveform with reference to one of the others. Choose $h_{2}$,
since that is our “middle” waveform. We can write,

$\displaystyle h_{1}$ $\displaystyle=$ $\displaystyle
Ah_{2}+\sqrt{1-A^{2}}h_{2\perp},$ (41) $\displaystyle h_{3}$ $\displaystyle=$
$\displaystyle Bh_{2}+\sqrt{1-B^{2}}h^{\prime}_{2\perp}.$ (42)

Both $h_{2\perp}$ and $h^{\prime}_{2\perp}$ are orthogonal to $h_{2}$, but are
not necessarily the same waveform, and the weights are chosen to ensure that
all waveforms are normalised. We can also write this as $A=\cos(\theta_{12})$
and $B=\cos(\theta_{23})$, and we therefore have,

$\displaystyle h_{1}$ $\displaystyle=$
$\displaystyle\cos(\theta_{12})h_{2}+\sin(\theta_{12})h_{2\perp},$ (43)
$\displaystyle h_{3}$ $\displaystyle=$
$\displaystyle\cos(\theta_{23})h_{2}+\sin(\theta_{23})h^{\prime}_{2\perp},$
(44)

and so,

$\langle
h_{1}|h_{3}\rangle=\cos(\theta_{12})\cos(\theta_{23})+\sin(\theta_{12})\sin(\theta_{23})\langle
h_{2\perp}|h^{\prime}_{2\perp}\rangle.$ (45)

If the two orthogonal contributions are the same, then the combined match will
be $\cos(\theta_{12}-\theta_{23})$, which is the best match we can have; if
$A=B$ then $C=1$, i.e., $h_{1}=h_{3}$. Alternatively, if $\langle
h_{2\perp}|h^{\prime}_{2\perp}\rangle=0$, then $C=AB$. In general we are
interested in cases where all matches are close to unity, and so if
$A=1-\mathcal{M}_{12}$ and $B=1-\mathcal{M}_{23}$, where $\mathcal{M}_{12}$
and $\mathcal{M}_{23}$ are the mismatches, then we will have $C\approx
1-\mathcal{M}_{12}-\mathcal{M}_{23}$, i.e., the mismatches add linearly.

The worst combined match, and therefore the upper bound on the combined
mismatch, occurs when $\langle h_{2\perp}|h^{\prime}_{2\perp}\rangle=-1$,
which can be thought of geometrically as $h_{1}$ and $h_{3}$ differing from
$h_{2}$ in opposite directions, and so we must add their differences to
calculate the combined difference. In this case, we have
$C=\cos(\theta_{12}+\theta_{23})$. This allows the extreme case where
$h_{1}=(h_{2}+h_{2\perp})/\sqrt{2}$ and $h_{3}=(h_{2}-h_{2\perp})/\sqrt{2}$,
and in this case we have $A=B=1/\sqrt{2}$, so the two waveforms are “equally
far apart”, but the combined match is $C=0$, and so $h_{1}$ and $h_{3}$ are
orthogonal to each other.

For the situations we are interested in, where the mismatches are small, we
recall that $\cos(\theta)\approx 1-\theta^{2}/2$, and so the mismatches can be
approximated as

$\displaystyle\mathcal{M}_{12}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\theta_{12}^{2},$ (46)
$\displaystyle\mathcal{M}_{23}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\theta_{23}^{2},$ (47)
$\displaystyle\mathcal{M}_{13}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\theta_{12}+\theta_{23})^{2},$ (48)

and therefore

$\mathcal{M}_{13}\approx\left(\sqrt{\mathcal{M}_{12}}+\sqrt{\mathcal{M}_{23}}\right)^{2}.$
(49)

This is not strictly an upper bound on the match that we calculate, because we
also optimise over time and phase shifts, and the optimisation is not captured
in our calculation, and will have a different effect on each individual match.
However, we have found in toy examples that Eq. (49) provides an excellent
estimate of the combined mismatch.

##  Appendix E Power-weighted precessing mismatch

The _match_ between two real valued detector response waveforms $h_{1}(t)$ and
$h_{2}(t)$ is defined to be the standard inner product weighted by the power
spectral density of the detector $S_{n}\left(f\right)$ maximised over various
sets of parameters $\Theta$ Cutler and Flanagan (1994), as given by Eq. (11).
The _mismatch_ may then by defined by Eq. (12). For precessing waveforms, the
set of parameters $\Theta$ that are maximised over are a relative time shift
$t_{0}$ between the waveforms, a relative phase shift $\phi_{0}$, and the
detector response polarisation angle $\psi_{0}$ Schmidt et al. (2015). The
precessing matches performed in this work are calculated as described in
Appendix B of Ref. Schmidt et al. (2015).

In order to see how the match varied over a range of total masses that might
be observed by current ground based detectors, we further calculated the
power-weighted match as is described in Ref. Ohme et al. (2011)

To perform matches over a frequency range
$f\in[f_{\mathrm{min}},f_{\mathrm{max}}]$ that extends below the minimum
frequency $f_{\mathrm{NR}}$ of one of the waveforms in this catalogue scaled
to a specified total mass, the full integral from $f_{\mathrm{min}}$ to the
maximum NR frequency can be approximated using a _power-weighted mismatch_
using the method described in Ref. Ohme et al. (2011). This method takes into
account the missing inspiral part of the waveform between $f_{\mathrm{min}}$
and the start of the NR waveform.

To perform a power-weighted mismatch the constituent waveforms are first split
up into contributions from NR defined over the frequency range
$f\in[f_{\mathrm{NR}},f_{\mathrm{max}}]$ and the contributions from the
inspiral below the lowest NR frequency defined over the frequency range
$f\in[f_{\mathrm{min}},f_{\mathrm{NR}})$,

$h(f)=\begin{cases}h_{\mathrm{ins}},&f\in[f_{\mathrm{min}},f_{\mathrm{NR}}),\\\
h_{\mathrm{NR}},&f\in[f_{\mathrm{NR}},f_{\mathrm{max}}].\end{cases}$ (50)

The power-weighted mismatch is then the mismatch in each region weighted by
the fraction of power in each region,

$\displaystyle\mathcal{M}_{\mathrm{pow}}$
$\displaystyle\equiv\frac{\|h\|^{2}_{(f_{\mathrm{min}},f_{\mathrm{NR}})}}{\|h\|^{2}}\mathcal{M}_{\mathrm{ins}}+\frac{\|h\|^{2}_{(f_{\mathrm{NR}},f_{\mathrm{max}})}}{\|h\|^{2}}\mathcal{M}_{\text{NR}},$
(51) $\displaystyle\mathcal{M}_{\mathrm{ins}}$
$\displaystyle\equiv\mathcal{M}_{(f_{\mathrm{min}},f_{\mathrm{NR}})}(h_{1,\mathrm{ins}},h_{2,\mathrm{ins}}),$
(52) $\displaystyle\mathcal{M}_{\mathrm{NR}}$
$\displaystyle\equiv\mathcal{M}_{(f_{\mathrm{NR}},f_{\mathrm{max}})}(h_{1,\mathrm{NR}},h_{2,\mathrm{NR}}),$
(53)