Catalog of precessing black-hole-binary numerical-relativity simulations

(

NR has been useful beyond modeling.NR waveform injections have been used in several studies, including to assess the presence of systematic bias in waveform models [64,65], and to estimate intermediate mass black hole binary merger rates [66].NR waveforms have also been used for direct comparisons and parameter estimation of gravitational waveobservations [67,68].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 [69][70][71][72][73][74][75][76][77].These fits have a number of applications, such as tests of general relativity [78][79][80][81].
Several large catalogs of BBH NR simulations exist [82][83][84][85][86][87][88].A quasicircular BBH is described by eight intrinsic parameters; the masses of each of the black holes m 1 and m 2 and their respective spins S 1 and 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 ⃗ χ i ¼ S i =m 2 i .The majority of simulations contained within these catalogs cover the precessing parameter space up to mass ratio q ¼ m 2 =m 1 ¼ 8 and dimensionless spin magnitudes χ ≤ 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 catalog is to provide a systematic covering of that parameter space.
The primary objective of this catalog was to support the development of a new precessing phenomenological model that is calibrated to numerical relativity waveforms [20].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 [15,16].For the first catalog used to inform the first precessing Phenom model we therefore chose no more than five points in each parameter direction.We also restrict to spins on the larger black hole, motivated by the results in Refs.[58,89], which suggest that two-spin effects will in general be difficult to measure in signals with signalto-noise ratios (SNRs) lower than ∼100.This choice was found to be sufficient; the PhenomPNR model of the dominant contribution to the signal [the (2,2) multipoles in the coprecessing frame] [20], constructed from 19 alignedspin waveforms and 40 precessing-binary waveforms, is of comparable accuracy to the equivalent contributions to the NRSurrogate model [29], 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 modeling community and to provide processed datasets that are appropriate for parameter estimation studies [64,66].
Our catalog contains datasets from 80 different configurations of precessing BBH systems.These configurations cover four mass ratios q ¼ m 2 =m 1 ∈ f1; 2; 4; 8g at four different spin magnitudes χ 2 ¼ j ⃗ S 2 j=m 2 2 ∈f0.2;0.4;0.6;0.8g each at five different spin orientations such that the angle between the orbital angular momentum and spin vector of the larger black hole is one of f30°; 60°; 90°; 120°; 150°g.The configurations are specified at a reference orbital frequency.The catalog can be accessed online at https:// data.cardiffgravity.org/bam-catalogue/.
In the following section we briefly summarize the methods used by the BAM code [4,5] 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 catalog.In Sec.IV we perform a waveform accuracy analysis to validate the catalog.Finally we conclude with Sec.V where we discuss what regions of parameter space and how the catalog can be used to contribute to the continuing advance of gravitational wave data analysis.

A. Simulation method
The simulations in this catalog were produced using BAM [4,5], a moving-box-based mesh-refinement numerical-relativity code that solves the 3 þ 1 decomposed Einstein equations.Specifically, for the simulations in this catalog we evolve Bowen-York wormhole data [90][91][92] via the χ variant of the moving-puncture treatment [2,3] of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation [93,94].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.[4,5] 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 Ψ 4 [95] following the procedure outlined in Ref. [96].Following these references, Ψ 4 is decomposed via the spherical multipolar decomposition where ψ 4;lm are the individual multipole moments, −2 Y lm are the spin-weighted spherical harmonics over the unitsphere defined by Θ and Φ.The quantity Ψ 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 modeling, 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 Ψ 4 data by ω 2 , where ω is the angular Fourier frequency [97].
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, 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 ¼ fN 0 ; N 1 ; …; N 1 max g, where l 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 max differs for each puncture so that d lmax;i =m i is approximately the same for each puncture.See Ref. [4] for more details on the BAM grid structure, and Ref. [98] for typical choices for numbers of levels and relative box sizes.Section II B 2 provides further details of the choices we made for the grid configurations.
The values N i , l 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 time step is specified as dt n ¼ 0.25d n .The finest level consists of two grids, one centered 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. [4] that the Berger-Oliger time stepping becomes unstable on the coarsest nonmoving 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.[4,5], except for the use of 0.25 rather than 0.5 for the Courant factor, which is necessary to sufficiently reduce the time stepping error in long simulations.
A pseudospectral elliptic solver is used to calculate binary wormhole initial data [92], with eccentricity reduced to < 2 × 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 B 1.

Initial data construction
We wish our simulations to begin at a user specified reference orbital frequency MΩ orb with spin vector S ≡ S 2 on the larger black hole (which we designate the secondary).The orientation of S 2 can be defined by the angle θ LS ¼ arccosð LN • Ŝ2 Þ between the spin vector S 2 and Newtonian orbital angular momentum vector L N , and the angle ϕ rS ¼ arccosðr • Ŝ2⊥ Þ between the projection of the spin vector on to the orbital plane S 2⊥ and the separation vector r from the larger black hole to the smaller.For all NR configurations in this work we choose ϕ rS ¼ 0 at the reference frequency MΩ orb , and a range of misalignment angles θ LS .The positions and momenta of the black holes consistent with these constraints must then be determined at MΩ orb , approximately chosen to minimize 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 blackhole positions and momenta at MΩ orb .Two methods were used for the simulations in this catalog.
For simulations with χ 2 ∈ f0.4; 0.8g the initial data parameters were determined by adapting the method used in previous work [15,99,100].For this method the physical parameters of the system ðq; S 1 ; S 2 Þ are specified at a much larger separation D start than the NR simulations will start at.The effective-one-body (EOB) equations of motion are then evolved up to MΩ 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; S 1 ; S 2 Þ at MΩ orb .During the course of inspiral from D start to MΩ orb for the single-spin precessing systems in this catalog, the angle θ LS can be seen to vary no more than ∼1°, while ϕ rS increases continuously (see Fig. 3 and Fig. 4 in Ref. [101]).To achieve a specific choice of ðθ LS ; ϕ rS Þ at a prescribed value of MΩ orb , the method was extended with an iterative refinement of the angle ϕ rS at D start until the parameters at MΩ orb are within a suitable tolerance of our desired values.Full details of this adapted method are given in Appendix A.
The simulations χ 2 ∈ f0.2; 0.6g were performed later, and were able to make use of a more recent method to produce low-eccentricity initial parameters, as described in Ref. [102].This method provides a post-Newtonian estimate of low-eccentricity parameters at a prescribed orbital frequency, making it possible to specify ðθ LS ; ϕ 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 catalog.We instead relied upon the manual perturbation approach outlined in Refs.[15,99,100] and described below.
While the initial data parameters generated from solving the EOBequations of motion will lead to low eccentricity simulations, in general this will not be low enough to satisfy our definition of a quasicircular binary.We placed an upper limit on the eccentricity at 2 × 10 −3 , based on the observation in Ref. [98] 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 ∼1000M, estimate the eccentricity, and make iterative small perturbations to the momenta of the component black holes [98,99].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%.Once initial data parameters are found that yield a sufficiently low eccentricity then a high-resolution production simulation is performed using the same parameters.
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; 700M similar to the method described in Ref. [103].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 [103].
In our production simulations we in general find that true eccentricity differs from that calculated in our lowerresolution eccentricity-reduction simulations.For a few of the cases in this catalog the eccentricity of the lowerresolution simulation was below our 2 × 10 −3 threshold, but the eccentricity of the production simulation exceeded it, as can be seen in Tables I and II.Nonetheless, only a handful of cases have eccentricities above 3 × 10 −3 , and only one is close to 4 × 10 −3 (CF_77).

Grid configurations
The simulations performed for this catalog 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.Section 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 [15,16], using the parameters ðq; 0; χ 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.[101,104,105], 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 ≤ 1=ð20f 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 × d 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.
Although details for each configuration are given in the public data release, we provide an example here to give a sense of the typical numbers of points, numbers of levels, TABLE I. Initial data parameters and relaxed properties of the precessing BBH configurations in this catalog 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ð LN • Ŝ2 Þ, the effective spins χ eff and χ 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 ½300; 800M using the method described in [103].The orbital frequency Mω orb is derived from the dynamics at relaxed times.The number of orbits N orb is from calculated from the relaxed time that Mω orb is reported at until the peak in the (2,2) multipole moment of Ψ 4 .

III. SIMULATION PROPERTIES
In this section we discuss the properties of the simulations in our catalog.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 quasicircular 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 summarize the properties of the remnant black holes.

A. Simulation configurations
Our catalog consists of dynamics and waveform data from NR simulations of 80 binary-black-hole configurations.
For each simulation two mass parameters m 1 , m 2 were chosen such that The initial data are iteratively constructed from these parameters such that the Arnowitt-Deser-Misner (ADM) mass of each puncture equals its respective mass parameter to within 0.02% [4].At subsequent times the masses of each black hole are recorded as the apparent-horizon masses M AH;1 ; M AH;2 of each puncture, which are related to the black hole hole masses m 1 , m 2 through the Christodoulou formula [106].This approach agrees well with the ADM mass of each puncture; the level of agreement is quantified further in Sec.III B.
In addition to the masses and spins, we must also choose the initial separation of the binary.For a binary undergoing noneccentric 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Ω orb .For this catalog we prefer to choose MΩ 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 B 1, with the goal of finding parameters consistent with noneccentric 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Ω orb ¼ 0.0225.This value was chosen to produce simulations of ∼2000M in length, which we expected to be sufficiently accurate for our modeling purposes, based on the experience of producing the aligned-spin simulations in Refs.[15,16].
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 ΔT ∼ 1=η, where η ¼ m 1 m 2 =M 2 is the symmetric mass ratio.Therefore, if simulations at massratios q ¼ 1 (η ¼ 0.25) and q ¼ 8 (η ≈ 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 blackhole 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.[107,108].For example, for our q ¼ 2 configurations with χ 2 ¼ 0.8, where all simulations begin at MΩ ¼ 0.0225, we see that the θ LS ¼ 30°configuration merges in 2254M, while the θ LS ¼ 150°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 [109] implementation of PhenomD [16].This provides a utility function XLALSimIMRPhenomD-ChirpTime that calculates the time until the peak in the (2,2)-strain of a specific system configuration given a starting gravitational wavefrequency, which is approximately twice the orbital frequency.The starting frequency is optimized using a simple interval bisection procedure until the peak time is ∼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 Tables I and II.Each configuration is characterized by its mass ratio q, the dimensionless spin magnitude χ 2 ¼ S 2 =m 2 2 , the spin angle θ LS ¼ arccosð LN • Ŝ2 Þ, the initial orbital frequency MΩ orb (or alternatively the initial binary separation D=M), and the binary's eccentricity, e.
For the final values reported in Tables I and II, the eccentricity e is estimated over the region ½300; 800M using the method described in [103].We also show the effective spin parameters χ eff and χ p .The effective aligned spin χ eff is defined in terms of the individual parallel spin components χ and parametrizes the dominant spin effect on the orbital phasing, as discussed in Refs.[107,108,111].The effective precession spin χ p is defined as [101], where 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 inplane spin of χ p , based on the leading-order spin-precession effects [112,113].In the single-spin configurations in this catalog, we will always have χ eff ¼ m 2 χ 2 cos θ LS =M and χ p ¼ χ 2 sin θ LS .
In Tables I and II we provide the quantities ðχ; θ LS ; χ eff ; χ p ; D=MÞ as specified in the initial data and (in brackets) at a relaxed time, t 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 modeling.We wish t rel to be as early as possible, to maximize the length of the usable waveform.We choose a relaxed time of where t peak is the time of the peak amplitude of the junk radiation in the (2,2) multipole moment of Ψ 4 and t damp is an estimate of the exponential decay time of the junk radiation, which we estimate as t damp ¼ 76m 2 .The damping time of the (2,2,0) quasinormal 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 [114]; we find that 76m is a reasonable choice for all of the configurations in this catalog.In Table I the initial orbital frequency MΩ orb is calculated at t rel and the number of orbits N orb is calculated from t rel to the merger time, which we define to be the time at which the peak in the l ¼ 2 multipole moments of Ψ 4 occurs.The retarded merger time, which we denote as t M , is given by the difference between the merger time and the tortoise coordinate where r is the distance from the punctures at which the Ψ 4 data is extracted.The final black hole has a mass of M f , a spin of χ f and a recoil velocity v R .We discuss the calculation of these remnant quantities in more detail in Sec.III D below.

B. 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 [91].The ADM puncture mass agrees well with CATALOG OF PRECESSING BLACK-HOLE-BINARY … PHYS.REV.D 109, 044032 (2024) the mass calculated from the area of the apparent horizon in the case of nonspinning black holes [115], but becomes less accurate for high spins [99].
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. [99] suggest that the error in the mass estimate could be on the order of ∼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 ∼0.01%, while for spins of 0.8 it is ∼0.04%. 2 , the spin tilt angle arccosð LN • Ŝ2 Þ, the effective spins χ eff and χ 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 ½300; 800M using the method described in [103].The orbital frequency Mω orb is derived from the dynamics at relaxed times.The number of orbits N orb is from calculated from the relaxed time that Mω orb is reported at until the peak in the (2,2) multipole moment of Ψ 4 .(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.

C. Initial black-hole spins
The black-hole spin S 2 is specified as part of the Bowen-York extrinsic curvature.The main source of uncertainty in the dimensionless spin χ 2 ¼ 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 χ 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 Table I that there is only a small discrepancy in the spin magnitude after the relaxation time.Since the initial blackhole 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 ∼0.001.
During the inspiral θ LS is not constant; it will oscillate due to the nutation of the orbital plane, as illustrated for one configuration in Fig. 2. Ideally, we would set up our simulations so that the mean value of θ 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 θ 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 θ 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 Tables I and II we report the mean value of θ LS at the start of the simulation and at the relaxed time.We estimate the value at the relaxed time by fitting to θ LS a sinusoidal ansatz of the form, where A 0 , A 1 , f 0 , f 1 , φ, 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 Tables I and II instead of the pointwise value of the NR data for θ LS .An example of this fit can be seen in Fig. 2. The resulting value of θ LS is used, along with the relaxed-time value of the spin magnitude, to calculate the relaxed-time values of χ eff and χ p .

D. Remnant properties
The final black hole that remains after the merger is characterized by its mass, spin, and recoil.We report each of these quantities in Tables I and II.
As with the relaxed-time quantities reported in Tables I and II, the mass and spin of the final black hole, M f and χ f , are calculated from the apparent horizon [116].As a consistency check we also estimate the mass and angular momentum of the final spacetime from the gravitationalwave 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 [4] and the result extrapolated to infinity.The final mass estimated in this way typically agrees with the horizon measure to within 5 × 10 −4 .Given the mass, perturbation theory provides a relationship between the black-hole spin and the frequency of the signal multipoles during ringdown [105,117].We calculate the ringdown FIG. 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, θ LS ¼ 60.The solid black line is for a configuration with dimensionless spin magnitude χ 2 ¼ 0.4 and the dashed black line represents χ 2 ¼ 0.8.FIG. 2. Linear trend of the angle θ LS between the Newtonian orbital angular momentum and spin of the spinning black hole as a function of simulation time t.θ LS as determined from NR data is plotted as a thick black line.The initial data value of θ LS is plotted as a horizontal thick black line.The dashed black line is a sinusoidal fit using Eq. ( 7).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.frequency of the ðl ¼ 2; jmj ¼ 2Þ multipoles by taking the Fourier transform of the waveform between 10 and 100M after merger [118], where merger is here defined as the time at which the sum of the square of the l ¼ 2 multipoles is maximized.The ringdown frequency is then the frequency at which the peak in this frequency domain postmerger waveform occurs [119].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 × 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 rel until the end of the simulation.We report here only the magnitude ofthe recoil velocity v R .The linear momentum is itself calculated as described in [4].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 catalog 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 χ ¼ 0.8.However, from the values presented in Tables I and II, 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 [113], and for full merger calculations in the numerical-relativity recoil studies cited in the Introduction.
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.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 [29] for those cases within the catalog that lie within the calibration region of NRSur7dq4 (q ≤ 4).To obtain this prediction, we used the value of the black hole spins, rotated into the LAL frame [120], 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 catalog, 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 catalog we studied a subset of four of the configurations described in Tables I and II.These four configurations are CF_47, CF_59, CF_66, and CF_80, with physical parameters ðq;χ;θ LS Þ ¼ fð4;0.4;60Þ;ð4;0.8;120Þ;ð8;0.4;30Þ;ð8;0.8;150Þg.The set of simulations used in the accuracy analysis of the (4,0.4,60)case were performed with a lower starting frequency of Mω 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 (8,0.8,150)configuration.These resolutions correspond to a number of grid points N ¼ f80; 96; 120; 144g in the boxes surrounding the punctures.Typically the width of the smallest box around each black hole is on the order of ∼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 B 2. We extracted the waveform data at R ext ¼ f50; 60; 70; 80; 90gM, 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 a convergence analysis, e.g., waveform amplitude and phase.In Sec.IV C 1 we sketch out how standard numerical convergence properties translate to the waveform mismatch, and provide more detailed derivations in Appendices C and D.

A. Data quality
Our accuracy analysis of the gravitational waveforms calculated from our simulations focuses on the mismatch uncertainty as detailed in Sec.IV C.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 coprecessing waveform in Sec.IV B. 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 binaryblack-hole waveforms, even for the dominant multipoles, we do not expect a convergence analysis of subdominant 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. [121].)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 postrelaxation-time ψ 4;lm multipole with an amplitude above 10 −5 .Depending on the application, of course, one may wish to apply a more or less stringent requirement.

B. 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 FIG. 4. Comparison of coprecessing frame time domain amplitudes for the l ¼ m modes for l ∈ f2; 3; 4; 5g.The left panel shows CF_7 with initial parameters ðq; χ 2 ; θÞ ¼ ð1; 0.4; 60Þ and the right panel shows CF_79 with parameters ðq; χ 2 ; θÞ ¼ ð8; 0.8; 120Þ.
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 is maximized and resampled using a constant time step of 0.1M.We then rotated the waveform into the coprecessing frame and aligned the coprecessing phases at merger.The coprecessing frame is one that precesses along with the binary and is advantageous here as it means we can focus on the error in a single multipole, (2,2), which is dominant in this frame, rather than considering the error in each of the l ¼ 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 coprecessing frame as well as two of the Euler angles α and β required to rotate the waveform from the inertial frame into the coprecessing frame.The Euler angles were calculated using the method detailed in Refs [122,123].(We do not consider accuracy of the third Euler angle, γ, because it is calculated from the other two; the accuracy of γ is of the same order as that for α.) 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; this would again lead to a factor of four increase is the extent of the radiation-extraction level (over what is currently required at R ext ¼ 90M, again requiring 2 6 times more memory.As such, we consider a closer spacing in numerical resolution and extraction radii in our convergence analysis.
We nominally expect the error due to extraction radius R ext to fall off as 1=R 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 [4,5].For both the extraction radius and numerical resolution, we determine the appropriate convergence order by studying the convergence behavior of the phase of the ðl ¼ 2; m ¼ 2Þ multipole in the coprecessing frame.We then identify the value of n for which the quantity, is minimized, where ϕ is the phase of the (2,2) multipole in the coprecessing frame, C ¼ as in Eq. (C5), and Δ i is the variable in the error expansion, i.e., numerical resolution or the inverse extraction radius.The quantity δðtÞ was calculated over the length of the waveform up to merger and the mean value δ ¼ hδðtÞi is shown in Fig. 5.This was done for both waveforms of varying resolution and extraction radius for the q ¼ 8, χ ¼ 0.8, θ LS ¼ 150°c onfiguration.In calculating the convergence order with respect to varying resolution we used waveforms with Δ f1;2;3g ¼ f1=144; 1=120; 1=96g while when considering the convergence order with respect to extraction radius we used Δ f1;2;3g ¼ f1=90; 1=70; 1=60g.From the results shown in Fig. 5, we conclude that the results are consistent with fourth-order finite differencing, which implies that the time stepping dominates the error budget.From inspecting time-dependent δð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 ∼1000M of the waveforms, but the fourth-order-accurate time stepping dominates in the last ∼500M before merger, and dominates overall.As expected, we see that the radiation extraction errors fall off as 1=R ext .
We should note that in Eq. ( 9) we have used the absolute value of the difference between the phases, which is a weaker test of convergence than considering ϕðΔ 1 Þ − ϕðΔ 2 Þ, for which we do not see the expected convergence rate.This is consistent with previous work (for example, the convergence analysis of the nonspinning q ¼ 18 configuration in Ref. [15]), and indicates that convergence is not clean enough to reliably apply Richardson extrapolation.The Richardson extrapolated quantities that we use below allow an approximate estimate of the error.This may be an overly optimistic estimate of the uncertainty in the numerically challenging ðq ¼ 8; χ ¼ 0.8Þ configurations, but we expect it to be a conservative estimate of the uncertainty for the rest of the catalog.FIG. 5.The value δ as given by Eq. ( 9) 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, χ ¼ 0.8, θ LS ¼ 150°was used in this analysis.
For the wave-extraction error, we have data at multiple extraction radii for every configuration, and can identify configurations at all mass ratios and spins (roughly half of all simulations) that satisfy clean 1=R ex falloff in the finiteextraction-radius error.In this section we consider only the ðq ¼ 8; χ ¼ 0.8; θ LS ¼ 150°Þ configuration, but we will use the full catalog when we consider mismatches below.
Assuming these convergence orders, we then calculate the Richardson-extrapolated values of the amplitude, phase, α and β, as functions of time, using Eq.(C4) in Appendix C. We used the resolutions N ¼ f120; 144g and the extraction radii R ext ¼ f80; 90g 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.This is clearly demonstrated in Fig. 6.
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 ∼0.1 radians for the medium resolution (N ¼ 96) simulation.Similarly, from Fig. 7, the dephasing due to the finite extraction radius is ∼0.4 radians for the waveform extracted at 90M.Since the simulations comprising the bulk of the catalog 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 Table III.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 coprecessing 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 Table III, we can see that the maximum relative error in the amplitude of the coprecessing 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 modeling 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.

C. 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 using the standard inner product weighted by the power spectral density of the detector S n ðfÞ, The match is then defined as the inner product optimized over various sets of parameters Θ [107], where the individual waveforms have been normalised so that hhjhi ¼ 1.We also define the mismatch, Since these are precessing configurations, we calculate precessing matches as described in Appendix B of Ref. [101].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 powerweighted match as described in Appendix E, based on the work in Ref. [124], using PhenomPv3 [57] as the model for the low-frequency part of the waveform.We then calculate the mismatch as given by Eq. ( 12).There are multiple ways to estimate the error contribution from the NR waveform at masses where the NR waveform will constitute only part of the full waveform, and both the power-weighted mismatch and a naive calculation that begins at the starting frequency of the NR waveform (which is the most common approach) are difficult to interpret.In Appendix.E we justify our use of the power-weighted mismatch in this work.
The matches performed here employ the noise spectrum of advanced LIGO at design sensitivity [125].We optimize the match over time and phase shifts and polarization.Since in all cases we are comparing numerical simulations of identical binary configurations, we do not optimize over the in-plane spin angle as is sometimes done for precessing systems, as an improvements in the match would likely be largely artificial.We report the SNR-weighted average match values, further averaged over a range of inclination angles between 0 and π.
We first extrapolate the mismatch due to finite resolution and extraction radius separately, assuming a particular falloff 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 Appendices C and D. We use this approach to combine estimates of the extraction-radius and numerical-resolution mismatches to give an estimate of the mismatch uncertainty of our NR waveforms.

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 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. [126]), 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.
We can consider the following form of the inner product, where we have explicitly included the normalizations N i for each waveform, N i ¼ R jh i j 2 df.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., where h 2 ðfÞ ¼ AðfÞe iϕðfÞ , and AðfÞ is the real amplitude and ϕðfÞ is the phase.We assume the true waveform to be normalized, so A waveform containing some amplitude error ΔA is given by Substituting Eq. ( 16) into ( 13) we find 2 df and we have assumed that ΔA is small in order to make the approximation in the final step.A mismatch based on this expression, 1 − Iðh 1 ;h 2 Þ, would be proportional to c − b 2 at leading order, and both of these scale as the amplitude squared.In a true mismatch, following Eqs.( 10) and ( 12), we would also take into account the frequency-dependent noise curve, and optimize over some parameters.The above calculation can be generalized to include the noise power spectrum by considering instead the whitened signals, hðfÞ → hðfÞ= ffiffiffiffiffiffiffiffiffiffiffi S n ðfÞ p , and the main result is the same.Optimizations will not in general change the dominant (amplitude square) scaling of the result.We can conclude, then, that the mismatch scales with the square of the amplitude error.
Similarly for a normalized waveform that contains some phase error Δϕ, where ϕ is the "true" phase.Substituting this expression into Eq.( 13) we find 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.
The waveform quantities at finite resolution or extraction radius can be expressed as a Richardson extrapolation of the appropriate expansion parameter (see Appendix C 1).The difference in the phase and amplitude between two waveforms (labeled A and B) is therefore equal to the difference between the leading-order error term (i.e.Δq ∼ Δ n B − Δ n A , see Appendix C 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 where Δ i is the value of the expansion parameter for the ith waveform and κ 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 A similar derivation to the one discussed here is also presented in Ref. [127].From this we can see that the ratio of the mismatch M between two pairs of waveforms is given by 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 MðA∶BÞ and MðB∶CÞ in order to get the mismatch between the final pair MðA∶CÞ is given by A more rigorous proof of this result is presented in Appendix D; and we note that this result also follows trivially from the triangle inequality, as shown in Ref. [128].
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 where M resolution is the mismatch due to the finite resolution of the numerical simulation and M extraction radius is the mismatch due to the finite distance from the source at which the waveforms were extracted.

Convergence order
We performed matches between waveforms extracted at R 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 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, χ ¼ 0.8, θ LS ¼ 150°) since this is the only case for which we have the very high resolution run.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 ext ¼ 80 and R ext ¼ 90 often does not follow the trend-we expect this is because the mismatch between these waveforms is so small [ð10 −6 Þ] that it is very sensitive to any data processing performed in the course of calculating the match.The mismatch between R ext ¼ 50 and R 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.

Extrapolation
Given the expected convergence order of the code, n (in this case n ¼ 4), 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.Figure 8 supports our expectation that the ðq ¼ 8; χ ¼ 0.8; θ LS ¼ 150°Þ configuration (CF_80) has the largest numerical-resolution uncertainty in the sense that the mismatches between the N ¼ 96 and N ¼ 120 simulations are larger than for any other configuration.This includes the long CF_47 simulation, which, even though it has longer to accumulate phase error than any other simulation (the merger time was at t M ∼ 2800M, compared to t M ∼ 2000M for the main simulations in the catalog), shows better agreement between the N ¼ 96 and N ¼ 120 runs than the shorter CF_80 simulations.We focus on the CF_80 simulations in the following discussion.
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 κ res using 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 κ res we can estimate the mismatch between the high or very high resolution waveforms with an infinitely FIG. 9. Mismatch between waveforms at varying extraction radii and the waveform extracted at R ext ¼ 90 M as a function of total mass.FIG. 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, χ ¼ 0.8, θ LS ¼ 150°configuration.
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 catalog of waveforms presented in Tables I and II.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, The result of this extrapolation procedure is shown in Fig. 14.We could only perform this calculation for the case (q ¼ 8, χ ¼ 0.8, θ LS ¼ 150°) 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.Based on a maximum mismatch between the medium and very high resolution waveforms MðΔ 96 ∶ Δ 144 Þ ¼ 2.45 × 10 −4 , the maximum mismatch between a medium resolution waveform and an infinitely well-resolved one is 6.0 × 10 −4 . 1  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 ext ¼ 90 M in Fig. 12, for a system with total mass 100M ⊙ .This fit gives the value of κ ext for every value of the  This value increases slightly to 6.5 × 10 −4 if we do not use power-weighted mismatches, but this small change does not change our final estimate for the total mismatch uncertainty.total mass of the system.From this we can calculate the mismatch between the waveform at Rext ¼ 90 M and the "true" waveform from MðΔ 90 ∶ Δ ∞ Þ ¼ κ ext 90 2 .The mismatch between the waveform extracted at R ext ¼ 90 M and the "true" waveform is shown in Fig. 15 for our three illustrative examples.The configuration that gives the greatest mismatch is q ¼ 4, χ ¼ 0.4, θ LS ¼ 60.From this case we take the maximum mismatch between a waveform at R ext ¼ 90 and at R ext → ∞ to be 1.4 × 10 −3 .We note that we are able to identify clean 1=R ex falloff in the uncertainty and perform clean extrapolation of the error estimate in roughly half of the configurations, and these also give results consistent with this value.
We estimate the mismatch between our medium resolution waveform extracted at R ext ¼ 90 M and the true FIG.12. Variation of the mismatch with extraction radius for a system with a total mass of 100 M. The mismatches shown are with respect to the waveform extracted at R ext ¼ 90 M. The line shows the fit based on Eq. ( 9) assuming first-order convergence.It is consistent with all extraction radii except at R ext ¼ 50 M. FIG. 13.Variation of the mismatch with resolution for a q ¼ 8, χ ¼ 0.8, θ LS ¼ 150°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 κ res assuming fourth-order convergence calculated using Eq. ( 24) and is consistent with all resolutions except the lowest one at N ¼ 80.

FIG. 14.
Projected mismatch between a waveform extracted at a resolution of N ¼ 96 and one that is infinitely well-resolved.waveform using Eq.(23).A conservative estimate of the mismatch between a waveform extracted at a finite distance of 90 M from the source for a simulation performed with a grid spacing d ¼ 0.0104 and the theoretical 'analytical' solution is therefore 3.83 × 10 −3 .This provides an estimate of the overall mismatch error of the waveforms presented in this catalog of 0.004 (or 0.4%).

V. CATALOG COMPARISON
A number of numerical relativity groups have started building larger and more comprehensive catalogs 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 [82,83] and the Rochester Institute of Technology (RIT) [84][85][86][87] and Maya catalogs [88].
The SXS Collaboration has produced the largest catalog to date with 2,019 BBH simulations spanning 1 ≤ q ≤ 10 and 0 ≤ χ ≤ 1.The RIT catalog contains 1,881 BBH simulations covering 1 ≤ q ≤ 128 and 0 ≤ χ ≤ 0.99 and the Maya catalog contains 452 unique BBH waveforms from more than 600 BBH simulations ranging between 1 ≤ q ≤ 15 and 0 ≤ χ ≤ 0.8.Unlike the simulations presented here, the SXS, RIT and Maya catalogs all contain simulations where the individual black hole spins can be zero or perfectly aligned/antialigned with the orbital angular momentum.Considering only the precessing parameter space, the SXS Collaboration has produced 1,429 simulations spanning 1 ≤ q ≤ 6 and 0 ≤ χ ≤ 0.99, the RIT catalog contains 561 simulations covering 1 ≤ q ≤ 15 and 0 ≤ χ ≤ 0.99 and the Maya catalog contains 315 waveforms ranging between 1 ≤ q ≤ 8 and 0.1 ≤ χ ≤ 0.8.
Figure 16 compares the parameter space coverage of the simulations presented here with the existing nonzero-spin simulations included in the SXS, RIT, and Maya catalogs over the mass ratio and larger black hole spin tilt angle and spin-magnitude parameter space.We see that although the existing catalogs provide good coverage for high black hole spins χ ≥ 0.5 and near equal mass ratios, there is a dearth of existing precessing simulations for low black hole spins χ < 0.25 and unequal mass ratios q ≳ 4.
Recent gravitational wave observations [see e.g., [52,53]] 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 [20].Most astrophysical models suggest that BBH at larger mass ratios will be rare (e.g., Ref. [129]), but given that there has been one observation to date at q ∼ 10 [53] extension of this parameter-space coverage to yet higher mass ratios will be necessary in the future.

VI. DISCUSSION
We have produced a catalog containing 80 waveforms from single-spin precessing systems with mass ratios up to q ¼ 8, dimensionless spin magnitudes up to χ ¼ 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 Oð0.05%Þ.We estimate the uncertainty in the initial spin magnitude to be Oð10 −3 Þ while the uncertainty in the initial spin inclination is Oð1°Þ.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 × 10 −4 , while the final spin magnitude is accurate to within 5 × 10 −3 .
The starting frequency of the simulations was chosen such that the simulations were all around a similar length (∼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 catalog presented here is sufficient to capture a wide range of single-spin precession effects.Most notably, the systems contained within it have a nonzero opening angle of the precession cone ranging from ∼1°to ∼115°radians.The cases with the largest opening angles display the most dominant precession effects.In particular, for initial configurations with mass ratios q ∼ 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 catalogs do not extend beyond q ¼ 5 so consequently, this region of parameter space is poorly covered by NR simulations.Indeed, in this catalog, 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 catalog, 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 catalog 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 modeling requirements of future gravitational wave observations.Existing catalogs (such as the SXS, Maya and RIT catalogs) provide a comprehensive coverage of the two-spin precessing parameter space up to q ¼ 4.This catalog 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 modeling that includes effects due to the in-plane spin direction, this catalog contains an incomplete and possibly random sampling of points.
There are many directions in which this catalog 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 catalog presented here is estimated to have required around 25 million CPU hours in total) and generic modeling 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 [129], 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 [44], only one has been found to have a mass ratio clearly greater than 8 [53].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 θ LS > 90°, 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 catalog, the cases CF_79 and CF_80 (at q ¼ 8 and χ ¼ 0.8) are already sufficiently short to limit model accuracy and we anticipate this will continue for decreasing spin magnitude as we go to highermass ratios [20].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 [58,89].Taking the detection threshold to be SNR 10, then we expect to be able to identify two-spin effects in 0.1% of detections.
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).
The public dataset can be accessed online [136].
calculating jΔMω 1 j as in the same way as jΔMω 0 j in the previous step.If jΔMω 1 j ≤ Mω tol the algorithm stops otherwise proceed to the next step.
(3) Third candidate parameters at D start (n ¼ 2) Set Δϕ to be 10°if jjΔMω jΔMω n j is calculated in the same way as previous steps.If jΔMω n j ≤ Mω 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Ω 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 interprocessor 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 parallelized highperformance 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 interprocessor 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-thecloud" infrastructure [138] to create container-based cluster instances using OCI Terraform on AMD64 128-core BM.Standard.E2.64 nodes, and ran at an average 1.8M/hr.Similar runs on cosma6 ran at about 3.7M/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 where Δ is the expansion parameter ( 1 N for resolution or 1 R 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 where e n Δ n is the leading-order error contribution, n is the convergence order of the simulation and n 0 > n.
Considering two waveforms computed using different expansion parameters Δ 1 and Δ 2 , we can solve the two simultaneous equations that arise from Eq. (C2) to give where is the Richardson extrapolation [116,139] of qðΔ 1 Þ.RðΔ 1 ; Δ 2 Þ has a higher-order error due to the truncation of the expansion in Δ than qðΔ gives the truncation error of the quantity q.

Convergence
Considering now three waveforms computed with expansion parameters Δ 1 > Δ 2 > Δ 3 we can eliminate q Ã in Eq. (C3).Neglecting higher-order error terms, the ratio of the difference between two sets of numerical waveforms with expansion parameter 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 [116,139] as a power series in an expansion parameter Δ.Consider two detector responses, h 1 and h 2 , at two resolutions with expansion parameters Δ 1 and Δ 2 , respectively.As seen in Eq. (C2), for an ith-order accurate finite difference method these can be represented by their truncated Richardson expansions 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 ðΔ n 2 − Δ n 1 Þ.Utilizing linearity in the inner product, we have 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 κ is defined as where we have used Eq.(C6) to re-write h 1 in terms of h and again neglected terms of higher order in the expansion parameter.κ can be seen to be bounded below by zero from the Cauchy-Schwarz inequality.The mismatch may then be approximated as It is important to note that the leading-order expansion parameter terms in the approximation Eq. (C11) are quadratic in the expansion parameters.In addition, while it is likely that the leading-order coefficient κ cannot be calculated directly, it is independent of any resolutionspecific expressions.As such κ 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, Equation (C13) 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 normalized, jh 1 j ¼ jh 2 j ¼ jh 3 j ¼ 1 according to the inner product of Eq. (10).The match additionally includes some optimizations, as in Eq. (11).We assume that h 1 and h 3 already incorporate such optimizations, so that the inner products hh 1 jh 2 i ¼ A and hh 2 jh 3 i ¼ B are equal to the matches Mðh 1 ; h 2 Þ and Mðh 2 ; h 3 Þ.We would now like to estimate an upper bound on C ¼ hh 1 jh 3 i, and on the mismatch Mðh 1 ; h 3 Þ.
We write each waveform with reference to one of the others.Choose h 2 , since that is our "middle" waveform.We can write, Both h 2⊥ and h 0 2⊥ are orthogonal to h 2 , but are not necessarily the same waveform, and the weights are chosen to ensure that all waveforms are normalized.We can also write this as A ¼ cosðθ 12 Þ and B ¼ cosðθ 23 Þ, and we therefore have and so, If the two orthogonal contributions are the same, then the combined inner product will be cosðθ , so the two waveforms are "equally far apart", but the combined inner product 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ðθÞ ≈ 1 − θ 2 =2, and so we can make the approximation, and therefore Note that, while M 12 and M 23 satisfy our definition of mismatches, M 13 does not until an additional optimization is performed on the inner product C.However, the resulting mismatch can only be lower (or equal) to M 13 , and so M 13 is an upper bound of the combined mismatch, Mðh 1 ; h 3 Þ, up to higher-order corrections.

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 ðfÞ maximized over various sets of parameters Θ [107], as given by Eq. (11).The mismatch may then by defined by Eq. (12).For precessing waveforms, the set of parameters Θ that are maximized over are a relative time shift t 0 between the waveforms, a relative phase shift ϕ 0 , and the detector response polarization angle ψ 0 [101].The precessing matches performed in this work are calculated as described in Appendix B of Ref. [101].
If the starting frequency f NR of the NR waveform is lower than the minimum frequency f min in the innerproduct calculation, then the mismatch can be used to produce a reliable estimate of the NR uncertainty in a waveform.If, however, f NR > f min , then the NR waveform does not represent the complete waveform that would be observed in the detector.Any waveform that extends to this lower frequency must include a contribution from some other source (e.g., a PN or EOB inspiral approximant), the uncertainty of which we do not know.For a given configuration f min scales with the total mass, and as we consider lower masses, the NR contribution to the total signal power will decrease; at sufficiently low mass that NR waveform will contribute negligible power, and its uncertainty will be irrelevant.One way to obtain a more realistic estimate of the NR contribution to the overall error budget is to weight the mismatch with the relative power of the NR contribution to the total signal, following the procedure described in Ref. [124].
To perform matches over a frequency range f ∈ ½f min ; f max that extends below the minimum frequency f NR of one of the waveforms in this catalog scaled to a specified total mass, the full integral from f min to the maximum NR frequency can be approximated using a power-weighted mismatch using the method described in Ref. [124].This method takes into account the missing inspiral part of the waveform between f 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 ∈ ½f NR ; f max and the The inspiral contribution khk 2 ðf min ;f NR Þ to khk can be calculated using any appropriate inspiral waveform.For this work the precessing waveform model PhenomPv3 [57] was used as the inspiral waveform.
It is important to make clear as described in [124] that M pow will be a lower bound to the mismatch M ⩾ M pow .
Although the inspiral parts are identical below f min , the optimization in the match calculation may introduce relative time and phase shifts that cause the mismatch to be nonzero over ½f min ; f NR .This means that a naive mismatch, that simply starts at f NR , is also not an upper bound; without incorporating further information about the inspiral contribution, neither the power-weighted mismatch nor a naive mismatch can be taken as a conservative estimate of the mismatch error.We can, however, consider both estimates and determine which is more appropriate.Figure 17 shows the relative contributions of PN and NR parts to the total mismatch in Eq. (E2) for the CF_66 configuration.We see that at around 60M ⊙ both contribute roughly equally.Below around 10M ⊙ the NR waveform contributes less than 10% to the final result.When we consider mismatches with respect to extraction radius, we find that the maximum mismatch (without power weighting) is at masses below 60M ⊙ , and so we do not consider the naive mismatch (starting at f NR ) to be a reliable estimate of the error.For these cases we consider the power-weighted mismatch (which peaks around 10M ⊙ ) to be more reliable.
In the case of the numerical-resolution mismatches, the power-weighting and naive estimates peak at similar masses.This is shown in Fig. 18 for the CF_80 configuration.Given that the NR contribution to the powerweighted mismatch is still above 75% at this mass, we consider the naive mismatch to be more reliable, and we find that it makes no appreciable difference to our estimate of the total mismatch error.

FIG. 3 .
FIG.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).

FIG. 6 .
FIG.6.Resolution dependence of the absolute error in the time domain coprecessing phase, relative to the Richardson-extrapolated phase.The phases have been aligned at merger.

FIG. 8 .
FIG. 8. Mismatch between waveforms at varying resolution against the high resolution (N ¼ 120) waveform as a function of total mass.

FIG. 11 .
FIG. 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.

FIG. 15 .
FIG. 15.Projected mismatch between a waveform extracted at a radius of R ext ¼ 90 M and one extracted infinitely far away.

FIG. 17 .
FIG. 17.Relative contributions of the NR and inspiral (denoted by PN) to the power-weighted mismatch.
Now N 0 ¼ 144, on the radiation extraction level N 6 ¼ 360 and the finest level has N 15 ¼ 96.The resolution on the coarsest level is N 0 ¼ 71.0M, the outer boundary is at R max ¼ 6816M (given that the merger time of our simulations was ∼2000M, in retrospect we could have removed the coarsest level).The resolution on the finest level is d 15 ¼ 0.0022M, or m 1 =51.2.The wave-extraction level is now l ¼ 7, and the resolution there is d 7 ¼ d 0 =2 7 ¼ 0.555M;

TABLE II .
Initial data parameters and relaxed properties of the precessing BBH configurations in this catalog 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

TABLE III .
Relative error in the waveform quantities compared with the Richardson extrapolated quantities for the q ¼ 8, χ ¼ 0.8, θ LS ¼ 150°configuration.
1 j − Mω tol j > Mω tol =2 otherwise set Δϕ to be 5°.ϕ 2 is chosen to be the target azimuthal spin angle ϕ c þ Δϕ. ϕ ω;2 and jΔMω 2 j are calculated in the same way as previous steps.If jΔMω 2 j ≤ Mω tol the algorithm stops otherwise proceed to the next step.(4)Further candidate parameters at D start (n > 2)Set Δϕ to be 10°if jjΔMω n−1 j − Mω tol j > Mω tol =2 otherwise set Δϕ to be 5°.If ϕ ω;2 > ϕ ω;1 this indicates that the azimuthal spin angle is being rotated in the wrong direction.As such, if ϕ ω;2 12 − θ 23 Þ, which is the best match we can have; ifA ¼ B then C ¼ 1, i.e., h 1 ¼ h 3 .Alternatively, if hh 2⊥ jh 0 2⊥ i ¼ 0, then C ¼ AB.In general we are interested in cases where all inner products are close to unity, and so if A ¼ 1 − M 12 and B ¼ 1 − M 23 , where M 12 and M 23 are the mismatches, then we will have C ≈ 1 − M 12 − M 23 ; the mismatches add linearly to produce C.The lowest value of C occurs when hh 2⊥ jh 0 2⊥ i ¼ −1 and we have C ¼ cosðθ 12 þ θ 23 Þ.This allows the extreme case whereh 1 ¼ ðh 2 þ h 2⊥ Þ=