The SXS Collaboration catalog of binary black hole simulations

TL;DR

This paper updates the SXS catalog with over 2000 numerical simulations of merging black holes, improving accuracy and coverage to aid gravitational wave detection and analysis.

Contribution

It provides a significantly expanded and more accurate catalog of black hole merger simulations, including new corrections and detailed error analysis.

Findings

Catalog contains 2018 configurations, 11 times larger than previous.

Median waveform length is 39 cycles, with improved numerical accuracy.

Remnant mass and spin uncertainties are an order of magnitude better than models.

Abstract

Accurate models of gravitational waves from merging black holes are necessary for detectors to observe as many events as possible while extracting the maximum science. Near the time of merger, the gravitational waves from merging black holes can be computed only using numerical relativity. In this paper, we present a major update of the Simulating eXtreme Spacetimes (SXS) Collaboration catalog of numerical simulations for merging black holes. The catalog contains 2018 distinct configurations (a factor of 11 increase compared to the 2013 SXS catalog), including 1426 spin-precessing configurations, with mass ratios between 1 and 10, and spin magnitudes up to 0.998. The median length of a waveform in the catalog is 39 cycles of the dominant $\ell=m=2$ gravitational-wave mode, with the shortest waveform containing 7.0 cycles and the longest 351.3 cycles. We discuss improvements such as correcting for moving centers of mass and extended coverage of the parameter space. We also present a thorough analysis of numerical errors, finding typical truncation errors corresponding to a waveform mismatch of $\sim 10^{-4}$. The simulations provide remnant masses and spins with uncertainties of 0.03% and 0.1% ($90^{\text{th}}$ percentile), about an order of magnitude better than analytical models for remnant properties. The full catalog is publicly available at https://www.black-holes.org/waveforms .

Tables

Table 1. Table 1: A comparison of BBH simulation catalogs. The mass of the larger black hole is m 1 subscript 𝑚 1 m_{1} , and the mass of the smaller black hole is m 2 subscript 𝑚 2 m_{2} . We use the convention that the mass ratio is m 1 / m 2 > 1 subscript 𝑚 1 subscript 𝑚 2 1 m_{1}/m_{2}>1 . The dimensionless spin magnitudes of the black holes are denoted | χ 1 , 2 | subscript 𝜒 1 2 |\chi_{1,2}| . The “SXS (2018)” row corresponds to the number of publicly available simulations at the end of 2018.

Catalog	Started	Updating?	Simulations	$m_1/m_2$ range	$\left|{\chi }_1\right|$ range	$\left|{\chi }_2\right|$ range	Precessing?	$\mathrm{Median}{N}_{\mathrm{cyc}}$	Public?
NINJA [98, 115]	2008	✗	63	1–10	0–0.95	0–0.95	✗	15	✗
NRAR [120]	2013	✗	25	1–10	0–0.8	0–0.6	✓	24	✗
Georgia Tech [122]	2016	✓	452	1–15	0–0.8	0–0.8	✓	4	✓
RIT (2017) [123]	2017	✓	126	1–6	0–0.85	0–0.85	✓	16	✓
RIT (2019) [124]	2017	✓	320	1–6	0–0.95	0–0.95	✓	19	✓
NCSA (2019) [125]	2019	✗	89	1–10	0	0	✗	20	✗
SXS (2018)	2013	✓	337	1–10	0–0.995	0–0.995	✓	23	✓
SXS (2019)	2013	✓	2018	1–10	0–0.998	0–0.998	✓	39	✓

Table 2. Table 2: Names, types, and descriptions of meanings of fields in metadata.json files.

Field name(s)	Type	Description
Identification
simulation_name	string	A SXS-assigned identifier of the simulation
alternative_names	string[ ]	Comma-separated list of alternative names, longer, more descriptive, and/or indicating the specific series of simulations this configuration belongs to. One of these alternative names is the ‘SXS:BBH:dddd’ id-number, which is guaranteed to be unique.
keywords	string[ ]	Deprecated
point_of_contact_email	string	Contact information for questions
authors_emails	string[ ]	Deprecated
References
simulation_bibtex_keys	string[ ]	References which should be cited if this simulation is used
code_bibtex_keys	string[ ]	Deprecated
initial_data_bibtex_keys	string[ ]	Deprecated
quasicircular_bibtex_keys	string[ ]	Deprecated
Input parameters for initial data
object{1,2}	string	Keyword description to identify the object type. One of {bh, ns}
initial_data_type	string	Type of initial data. One of
		BBH_CFMS – conformally flat, maximal slice
		BBH_SKS – superposed Kerr-Schild
initial_separation	double	Coordinate separation $D_0$ between centers of compact objects, as passed to the initial data solver [233, 46, 133]
initial_orbital_frequency	double	Initial orbital frequency ${\Omega }_0$ passed to the initial-data solver [46, 133]
initial_adot	double	Radial velocity parameter ${\dot{a}}_0$ passed to the initial data solver [46, 134]
Measurements of initial data
initial_ADM_energy	double	ADM energy of the initial data
initial_ADM_linear_momentum	double[3]	ADM linear momentum of the initial data
initial_ADM_angular_momentum	double[3]	ADM angular momentum of the initial data
initial_mass{1,2}	double	Christodoulou mass of the apparent horizon of each body in initial data (Eq. (5); code units)
initial_dimensionless_spin{1,2}	double[3]	Dimensionless spins of the BHs in the initial data (cf. Sec. 2.2)
initial_position{1,2}	double[3]	Initial coordinates of the center of each body
Reference quantities
relaxation_time	double	Time at which we deem junk radiation to have sufficiently decayed (code units)
reference_time	double	Time at which the quantities are extracted from the evolution (code units)
reference_mass{1,2}	double	Christodoulou masses of the black holes at reference time, cf. Eq. (5) (code units)
reference_dimensionless_spin{1,2}	double[3]	BH spins at reference time (cf. Sec. 2.2.1)
reference_position{1,2}	double[3]	Coordinates of the centers of the two bodies at reference time
reference_orbital_frequency	double[3]	Orbital angular frequency vector at reference time
reference_mean_anomaly	double	Mean anomaly at reference time
reference_eccentricity	double	Orbital eccentricity at reference time [45]
Merger/remnant quantities
number_of_orbits	double	Number of orbits until formation of a common apparent horizon
common_horizon_time	double	Evolution time at which common horizon is first detected
remnant_mass	double	Final mass of the remnant black hole after merger, cf. Sec. 2.2.3
remnant_dimensionless_spin	double[3]	Spin of the remnant black hole after merger, cf. Sec. 2.2.3
remnant_velocity	double[3]	Linear velocity of the remnant black hole after merger, cf. Sec. 2.2.3
Code information
metadata_version	int	Version of the metadata format itself. All simulations in this release of the catalog have version number 1
spec_revisions	string[ ]	(Array of) git revisions of the evolution code
spells_revision	string[ ]	(Array of) git revisions of initial data solver

Table 3. Table 3: Data files that contain GW modes. The first two entries are the preferred ones; all other files should be used only when a clearly understood need arises. The “time” and “data” labels indicate the contents of the time and data columns, where u 𝑢 u denotes the retarded time, corrected for finite-radius effects via Eq. ( 12 ), and T 𝑇 T is the raw unnormalized time coordinate of the underlying NR simulation. The last two columns indicate whether extrapolation and center-of-mass corrections were applied [cf. Secs.  2.4.1 and  2.4.2 ].

Filename	time	data	extrap?	COM?
rhOverM_Asymptotic_GeometricUnits_CoM.h5	$u/M$	$\frac{r}{M}{h}^{\ell ,m}$	✓	✓
rMPsi4_Asymptotic_GeometricUnits_CoM.h5	$u/M$	$rM{\psi }_4^{\ell ,m}$	✓	✓
rhOverM_Asymptotic_GeometricUnits.h5	$u/M$	$\frac{r}{M}{h}^{\ell ,m}$	✓	✗
rMPsi4_Asymptotic_GeometricUnits.h5	$u/M$	$rM{\psi }_4^{\ell ,m}$	✓	✗
rh_FiniteRadii_CodeUnits.h5	$T$	$r{h}^{\ell ,m}$	✗	✗
rPsi4_FiniteRadii_CodeUnits.h5	$T$	$r{\psi }_4^{\ell ,m}$	✗	✗

Table 4. Table 4: The contents of each group of the Horizons.h5 file; these datasets are provided for each of the individual apparent horizons and the common apparent horizon.

Dataset name	columns	time	data
ArealMass.dat	2	$T$	$M_{\mathrm{irr}}$	Eq. (6)
ChristodoulouMass.dat	2	$T$	$M$	Eq. (5)
CoordCenterInertial.dat	4	$T$	$\overrightarrow{x}$	Eq. (8)
DimensionfulInertialSpin.dat	4	$T$	$\overrightarrow{S}$	Eq. (4)
DimensionfulInertialSpinMag.dat	2	$T$	$S$	Eq. (1)
chiInertial.dat	4	$T$	$\overrightarrow{\chi }$	Eq. (7)
chiMagInertial.dat	2	$T$	$\left|\overrightarrow{\chi }\right|$	(Euclidean norm)

Equations

$$S_{\varphi }\equiv \frac{1}{8\pi }{\int }_{\mathcal{H}}{\varphi }^i{s}^j{K}_{ij}dA.$$

$$S\equiv \sqrt{S_{{\varphi }_{(1)}}^2+S_{{\varphi }_{(2)}}^2+S_{{\varphi }_{(3)}}^2}.$$

$${\hat{\chi }}_{\Omega m}\equiv \frac{1}{N}{\int }_{\mathcal{H}}\vec{r}\Omega dA,$$

$$\vec{S}\equiv S{\hat{\chi }}_{\Omega m}.$$

$$M^2\equiv M_{\text{irr}}^2+\frac{S^2}{4M_{\text{irr}}^2},$$

$$M_{\mathrm{irr}}^2\equiv \frac{A}{16\pi }=\frac{1}{16\pi }{\int }_{\mathcal{H}}dA.$$

$$\vec{\chi }\equiv \frac{\vec{S}}{M^2}.$$

$$\vec{x}=\frac{1}{A}{\int }_{\mathcal{H}}\vec{r}dA.$$

$${\Psi }_4=\sum _{\ell ,m}{\Psi }_4^{\ell ,m}{}_{-2}{Y}_{\ell ,m}(\theta ,\varphi ),$$

$$r{h}^{\ell ,m}=\sqrt{(\ell -1)\ell (\ell +1)(\ell +2)}({\Phi }_{\ell ,m}^{(+)}+i{\Phi }_{\ell ,m}^{(-)}).$$

$$h=(h_{+}-i{h}_{\times })=\sum _{\ell ,m}{h}^{\ell ,m}{}_{-2}{Y}_{\ell ,m}(\theta ,\varphi ).$$

$$u_{i,j}={\int }_0^{T_i}\sqrt{\frac{-1/g_j^{TT}(T)}{1-2M_{\text{ADM}}/r_j(T)}}dT-r_j(T_i)-2M_{\text{ADM}}\ln \left[\frac{r_j(T_i)}{2M_{\text{ADM}}}-1\right],$$

$${\hat{h}}^{\ell ,m}(u_k,r)\approx \sum _{j=0}^N\frac{{\hat{h}}_{(j)}^{\ell ,m}(u_k)}{r^{j+1}}.$$

$${\vec{x}}_{\text{COM}}\equiv \frac{m_1(t){\vec{x}}_1(t)+m_2(t){\vec{x}}_2(t)}{m_1(t)+m_2(t)},$$

$$⟨Q⟩\equiv \frac{1}{t_f-t_i}{\int }_{t_i}^{t_f}Q(t)dt.$$

$$\vec{\alpha }=\frac{⟨t^2⟩⟨{\vec{x}}_{\text{COM}}⟩-⟨t⟩⟨t{\vec{x}}_{\text{COM}}⟩}{⟨t^2⟩-⟨t{⟩}^2}$$

$$\vec{\beta }=\frac{⟨t{\vec{x}}_{\text{COM}}⟩-⟨t⟩⟨{\vec{x}}_{\text{COM}}⟩}{⟨t^2⟩-⟨t{⟩}^2}$$

$$e_0=\frac{B_{\Omega }}{2{\Omega }_0{\omega }_{\Omega }},$$

$${\chi }_{\mathrm{eff}}\equiv \frac{(m_1{\vec{\chi }}_1+m_2{\vec{\chi }}_2)\cdot \hat{L}}{m_1+m_2}=\frac{m_1{\chi }_{1\parallel }+m_2{\chi }_{2\parallel }}{m_1+m_2}.$$

$${\vec{\chi }}_{A\perp }\equiv {\vec{\chi }}_A-({\vec{\chi }}_A\cdot \hat{L})\hat{L}$$

$$\cos {\theta }_{AL}\equiv {\hat{\chi }}_A\cdot \hat{L},$$

$$\cos \Delta \Phi \equiv {\hat{\chi }}_{1\perp }\cdot {\hat{\chi }}_{2\perp },$$

$$⟨h^1,h^2⟩={\int }_{-\infty }^{+\infty }\frac{{\tilde{h}}^1(f){\tilde{h}}^2{}^{\ast }(f)}{S_n(|f|)}df,$$

$$\mathcal{O}(h^1,h^2)=\text{Re}\left[\frac{⟨h^1,h^2⟩}{\sqrt{⟨h^1,h^1⟩⟨h^2,h^2⟩}}\right].$$

$$\mathcal{M}(h^1,h^2)=1-\underset{\delta \varphi ,\delta t}{\max }\mathcal{O}(h^1,h^2,\delta \varphi ,\delta t).$$

$$\mathcal{M}<\frac{D}{2{\rho }^2},$$

$$\frac{u}{M}$$

$$h^1(t,\theta ,\varphi )$$

$$h^2(t,\theta ,\varphi )$$

$${\log }_2{N}_{\text{resamp}}=1+\lceil {\log }_2\underset{\text{LevN}}{\max }{N}_{\text{samp}}\rceil .$$

$$\mathcal{O}(\delta \varphi ,\delta t)=$$

Full text

1 Introduction

Advanced LIGO [1] and Virgo [2] inaugurated the era of gravitational-wave astronomy in 2015 by observing gravitational waves passing through Earth for the first time [3, 4, 5]. This first gravitational-wave signal, named GW150914, was emitted during the merger of two black holes [6]. Subsequently, gravitational-wave signals have been observed from a merger of two neutron stars, GW170817 [7], and from nine further black-hole mergers [8, 9, 10, 11, 12].

Making the most sensitive searches for binary coalescence in noisy detector data requires accurate gravitational-wave templates. Further, inferring properties of the sources of these signals requires comparing the data against millions of accurate templates. During the late stages of a compact binary merger, when the components move at relativistic speeds and spacetime becomes nonlinearly dynamical, analytic approximations to the binary dynamics [13] break down. This strong-gravity regime reveals the behavior of curved spacetime under the most extreme conditions, such as the nonlinear dynamics of merging black holes, the formation and relaxation of dynamical horizons [14], and the nature of the remnant black hole left behind following the merger of a binary black hole (BBH) [15]. The strong-gravity regime also has the potential to place strong upper bounds on deviations from general relativity (or to reveal such deviations if they exist) [16, 17, 18, 19].

In this highly nonlinear regime, accurate solutions of Einstein’s equations of general relativity require numerical-relativity calculations: direct solution of the full dynamical field equations using high-performance computing (for summaries, see [20, 21, 22, 23] and references therein), which became possible in 2005 [24, 25, 26]. Such simulations are essential in exploring the dynamics of spacetime curvature itself. They have revealed the simplicity of the merger phase [27] and the potentially strong recoil of the remnant (e.g. [28, 29, 30]), motivating studies of the interplay between the linear momenta of the black holes and of the surrounding spacetime [31, 32]. Simulations have also been used for visualizations of curved spacetime [33, 34, 35, 36, 37, 38, 39, 40], investigations of spin quantities [41], and the relaxation of spacetime to the Kerr solution following merger [42, 43, 44]. The motion of the black hole horizons and horizon curvature quantities have been used to explore eccentric dynamics [45, 46, 47, 48], spin precession [49, 50, 51, 52], and the first law of binary black hole mechanics [53, 54, 55, 56, 57]. These in turn have been compared to analytic post-Newtonian and self-force approximations (see also [58, 59, 60]), mapping out the bounds of validity of these approximations.

A key application of BBH simulations is the accurate modeling of gravitational waves emitted by these systems during their late inspiral, merger, and final ringdown. Waveforms extracted from BBH simulations are essential for analyzing observed gravitational-wave signals from black hole binaries. Indeed, all BBH observations by LIGO and Virgo were analyzed using waveform families that rely on numerical relativity for their construction, most notably effective-one-body waveform models [61, 62, 63, 64, 65] and phenomenological waveform models [66, 67, 68]. Numerical simulations are also central in validating such waveform models [69, 70, 71, 72, 73, 74, 75, 76], and were used to validate GW searches [77, 78, 79]. Waveforms from numerical relativity are also used directly in parameter estimation [80, 81], to construct template banks [82], and to construct waveform families without intermediate analytical models, using methods such as reduced order modeling [83, 84, 85, 86]. Today’s simulations span ground-based detectors’ frequency bands only for total masses $\gtrsim 50M_{\odot}$. For smaller total masses, numerical-relativity waveforms can be “hybridized” by attaching them to the end of analytic inspiral waveforms to produce waveforms that span detectors’ sensitive bands [84, 87, 88, 89, 90].

Most applications of numerical relativity to gravitational-wave astronomy [91, 86, 64, 67, 79] require merger simulations for a large number of different BBH configurations. Several groups have answered this challenge by creating public catalogs of numerical-relativity simulations for merging black holes [92, 93, 94]. Even with this considerable progress, computing simulation catalogs that meet the needs of current and future gravitational-wave observatories remains quite challenging. Simulations must be long enough to span LIGO’s sensitive frequency band or to allow reliable hybridization. This is straightforward for black-hole binaries with sufficiently high total masses (such as GW150914, which remained in LIGO’s band for only $\sim$0.2 seconds [3]) but not yet possible for BBHs with lower masses (such as GW170608, which remained in LIGO’s band for $\sim$2 seconds [10]). Moreover, simulations must achieve sufficient accuracy [as discussed further in the text surrounding Eq. (25)], which is especially challenging when computing the waves’ higher-order spherical-harmonic modes. Additionally, catalogs must span a vast, 7-dimensional parameter space, including the different possible mass ratios and black-hole spins. While a small number of simulations with extremely high mass ratios (up to $m_{1}/m_{2}\sim 100$ [95, 96]) and nearly extremal spins [97] have been achieved, much of the high-mass-ratio, high-spin parameter space remains completely unexplored.

Early efforts to explore the extensive parameter space of BBHs with numerical relativity highlighted the importance of collaborative efforts for development of numerical-relativity codes, production of numerical-relativity waveforms, error estimation, analysis of the resulting waveforms, and in storing and making broadly available the resulting data. This in turn inspired the creation of several collaborative projects rooted in numerical relativity and gravitational wave analysis, with the goal of creating the simulation catalogs that gravitational-wave astronomy requires. The first such collaboration resulting in a catalog was the Numerical INJection Analysis (NINJA) project that began in 2008 [98]. The goal of NINJA was to test how well ground-based gravitational-wave detectors could find gravitational waves in their data, by injecting realistic gravitational-wave signals (created from numerical-relativity waveforms) into detector noise [77, 98]. NINJA was a pioneering example of collaboration between members of the numerical relativity and gravitational wave analysis communities. Its catalog included simulations generated with a variety of numerical-relativity codes, including BAM [99], CCATIE [100, 101, 102, 103, 104], Hahndol [105, 106], LAZEV [107], LEAN [108], MayaKranc [109, 110], the Princeton University code [111, 24, 27, 112], the University of Illinois code [113], and our code SpEC [114] (cf. Section 2.1). Members of the NINJA project created a catalog of 23 simulations for creating injected signals, but this initial catalog spanned a very limited subspace of the 7-dimensional BBH parameter space.

In 2012, a followup effort, “NINJA-2” yielded a catalog [115, 78] of 63 simulations spanning a broader region of the parameter space, but still not including any simulations with precessing black-hole spins. As the simulations were created by various groups and collaborations, the simulations in the NINJA-2 catalog use a variety of codes, including SpEC [114] as well as the moving puncture/BSSN [116] codes BAM [99], LAZEV [107], LEAN [108], Llama [117], and Maya/MayaKranc [109, 110]. Note that many of the moving puncture codes utilize a common Cactus infrastructure [118] and the Einstein Toolkit [119], a finite-volume discretization code. In 2013 the Numerical Relativity Analytical Relativity (NRAR) collaboration [120] computed and published 25 simulations, with the main focus being the comparison of the numerical relativity waveforms with the analytical waveform families in use by LIGO.

Building on the NINJA and NRAR efforts, a number of numerical-relativity groups have begun building larger, more comprehensive catalogs, spanning more of the parameter space (including spin precession) and including more orbits before merger (enabling the waveforms to span the detectors’ frequency bands for lower total masses). These catalogs are summarized in Table 1. In 2013 the SXS collaboration released a catalog [121] with 174 BBH simulations created using the Spectral Einstein Code (SpEC) [114]. By the end of 2018, this catalog had grown to include 337 simulations for BBHs, seven from binaries with a black hole and a neutron star, and two from a pair of neutron stars. The full SXS catalog (including the new simulations presented here) is publicly accessible at https://www.black-holes.org/waveforms [92]. In May 2016 the Georgia Tech group released a catalog [94] of 452 distinct BBH simulations (from a pool of more than 600 BBH simulations [122]). The Rochester Institute of Technology (RIT) group released a catalog [93] in 2017 that included 126 simulations [123], as well as an updated catalog in 2019 containing a total of 320 simulations [124].

In this paper, we present a major update to the SXS Collaboration’s catalog. Our catalog, created using SpEC, now includes 2018 simulations, an increase of a factor of 11 over our 2013 catalog. The median waveform length is now $N_{\rm cyc}=$ 39 cycles of $\ell=m=2$ gravitational waves, while in our initial catalog [121], only half of the simulations had more than 24 gravitational-wave cycles. Here $N_{\rm cyc}$ is approximated by doubling the number of orbits during inspiral up to merger (when a common horizon forms), as determined by the coordinate trajectories of the black holes. The increased number of cycles means that typical waveforms in our catalog now tend to span LIGO’s sensitive band over a broader range of total masses. We estimate the numerical uncertainty for most waveforms in the catalog, finding a typical mismatch of $\sim 10^{-4}$ between the highest and second-highest resolutions. Our catalog now includes spins up to 0.998 and mass ratios up to 10, with considerably better coverage in the space of mass ratios up to 2 and spins up to 0.8 (motivated by the estimated properties of GW150914). In addition we have re-run a number of earlier simulations with a more modern version of SpEC. In some cases this improves the precision of higher order modes and removes the imprint of gauge changes from the coordinate trajectories. Since our past simulations have been widely used, we retain all versions of the simulations in our catalog with different labels; further details are given in App. A. Our catalog is publicly available at https://www.black-holes.org/waveforms [92], in a format based on that of the NRAR project [126].

The rest of this paper is organized as follows. Section 2 summarizes the numerical methods that we employ in SpEC. Then, Sec. 3 summarizes the areas of the parameter space our catalog covers (and what areas it does not yet cover). Section 4 estimates the accuracy of the catalog’s waveforms. After comparing the remnant properties to analytic fits in Sec. 5, we conclude in Sec. 6. We document the formats of our publicly available data in App. A, our definitions for calculation of mismatches in App. B, and sign conventions in App. C.

2 Summary of methods

2.1 Spectral Einstein Code

We use the Spectral Einstein Code (SpEC) [114] to model merging black holes and the gravitational waves they emit. The first step in a binary black-hole simulation is constructing initial data. We construct constraint-satisfying BBH initial data using the Extended Conformal Thin Sandwich (XCTS) [127, 128] equations. In the XCTS formulation of the BBH initial value problem, the initial spatial slice has i) a spatial metric proportional to a freely chosen conformal metric and ii) a freely chosen trace of its extrinsic curvature. Typically, when solving the XCTS equations, we choose the conformal metric and trace of extrinsic curvature to be weighted superpositions of the analytic solutions for two single black holes in Kerr-Schild coordinates [129], but some (typically older) simulations in the SXS catalog instead are conformally flat and have a vanishing trace of the extrinsic curvature (i.e. maximal slicing) [130]. In the XCTS formulation, the conformal metric and trace of the extrinsic curvature have freely chosen time derivatives; we construct quasi-equilibrium initial data by setting these time derivatives to zero. We then solve the XCTS equations on a grid with two excised regions using a spectral elliptic solver [131], with boundary conditions on the excision boundaries chosen to ensure that these boundaries are apparent horizons [130, 129]. The solution yields initial data for a BBH evolution, including the initial spatial metric, the initial extrinsic curvature, and the initial lapse and shift, which determine the initial coordinate choice.

We iteratively construct BBH initial data, tuning the initial data to achieve a BBH with the desired properties. Our iterative scheme uses two nested loops. The inner loop solves the XCTS equations, adjusting our choices for the free data (conformal metric, trace of extrinsic curvature, and the time derivatives of each) and boundary conditions, until the resulting BBHs have the desired mass ratio and spins [132, 133]. The outer loop briefly (typically for a few orbits) evolves the initial data resulting from the inner loop, and adjusts the initial coordinate velocities to yield a BBH with small orbital eccentricity [134, 46, 47], typically $e_{0}\sim 10^{-4}$ as defined in Eq. (17). For some simulations in the SXS catalog, we intentionally omit the eccentricity-reduction loop, to obtain initial data for BBHs with non-negligible orbital eccentricity.

We evolve the initial data using a first-order version of the generalized harmonic formulation [135, 136, 111, 137] of Einstein’s equations with constraint damping [138, 111, 137]. We choose an initial gauge that approximates a time-independent solution in a co-rotating frame, and then we smoothly change to damped harmonic gauge [139, 140, 141], which we have found to work well numerically near the time of merger.

We evolve the initial data using a multidomain spectral method [142, 137, 143, 141, 144]. Timestepping is done via the method of lines using a fifth-order Dormand-Prince integrator with a proportional-integral adaptive timestepping control system that chooses an appropriate step size while achieving a desired time-stepping error [145]. The computational domain extends from pure-outflow excision boundaries conforming to the shapes of the apparent horizons (AH) [143, 141, 144, 146] to an artificial outer boundary where we impose constraint-preserving boundary conditions [137, 147, 148]. We also impose boundary conditions on the incoming characteristic fields at each internal boundary of the computational domain [149, 150]. After a common AH forms, the simulation automatically stops, interpolates onto a new grid with only a single excision boundary [143, 144], and continues evolving on this new domain through ringdown, until the ringdown gravitational waves have left the domain.

Spectral methods are exponentially convergent, i.e. spatial truncation errors in a given subdomain of fixed size and shape decrease exponentially with the number of collocation points. Our simulations use multiple subdomains, and the size, shape, and even the number of these subdomains changes dynamically as a simulation proceeds ($h$-refinement), as controlled by our spectral adaptive mesh refinement (AMR) procedure [151, 152]. Moreover, we choose the accuracy of a simulation not by picking the number of grid points directly, but instead by specifying a tolerance parameter that governs when AMR should add or subtract grid points within a given subdomain ($p$-refinement), or when it should split or join subdomains ($h$-refinement). As a result, we should not always expect strict convergence as a function of the AMR tolerance parameter. Convergence may fail in several ways. For example, two otherwise-identical simulations with different AMR tolerances may happen to have the same number of grid points in a particular subdomain at a particular time, because their local truncation errors are below or above both thresholds. Alternatively, the two simulations may have entirely different subdomain boundaries in a given region at some other time. Furthermore, AMR’s decisions exhibit hysteresis. Despite these issues, most of our simulations do exhibit convergence with AMR tolerance, as shown in Sec. 4.

2.2 Black hole masses, spins, and centers

2.2.1 Quasi-local definitions

Defining mass and spin in general relativity is nontrivial; for recent reviews, see [153, 154, 155]. The simplest definitions apply to asymptotically flat spacetime, in the limit of approaching either spacelike infinity $i^{0}$ (i.e. the ADM mass and angular momentum [156]) or future null infinity $\mathscr{I}^{+}$ (the Bondi mass and angular momentum [157]). These ADM and Bondi quantities give global values for the entire spacetime, including binding energy and energy in gravitational waves, but they do not yield the masses and spins of the individual component black holes. One possible approach for determining component masses and spins, while the black holes are far from merger, is to perform asymptotic matching in a buffer region [158]. This approach, however, has not been explored in numerical simulations, and becomes invalid as the black holes approach merger.

Instead, NR simulations rely on “quasi-local” mass and spin [153] measurements from AHs. Quasi-local masses and spins recover the Kerr mass and spin when evaluated on the AH of a Kerr black hole, and also evolve in agreement with tidal torquing and heating approximations [97]. The definition of quasi-local mass that we adopt relies on our chosen measure of quasilocal spin. Given an apparent horizon $\mathcal{H}$ within the current constant-time hypersurface $\Sigma$, and a vector field $\phi^{i}$ tangent to $\mathcal{H}$, the component of spin angular momentum inside this 2-surface generated by $\phi^{i}$ is given by [159] (following the conventions of [129])

[TABLE]

Here $s^{i}$ denotes the outgoing spacelike unit normal to $\mathcal{H}$ tangent to $\Sigma$, $K_{ij}$ is the extrinsic curvature of $\Sigma$ (with sign conventions given in Appendix C), and $dA$ is the induced proper area element on $\mathcal{H}$.

In axisymmetry, making $\phi^{i}$ the symmetry’s rotational Killing vector field yields the corresponding conserved angular momentum. But black holes in a binary merger are only approximately axially symmetric long before merger and after ringdown. With no exact rotational Killing vector available, instead we follow Refs. [160, 129] and compute approximate Killing vectors on the apparent horizon. For details on our procedure for measuring quasi-local spin, and a discussion of other approaches, see [51]; here, we briefly summarize our method.

We solve an eigenvector problem to find the three tangent vectors $\{\phi_{(1)}^{A},\phi_{(2)}^{A},\phi_{(3)}^{A}\}$ that come closest to solving the Killing equation. Here, capital Latin indices $A,B,\dots$ run over the two-dimensional tangent space of $\mathcal{H}$. The Killing equation implies that a Killing vector is divergence-free, which means the vector has vanishing expansion and shear. We choose to begin with expansion-free vectors $\phi_{(i)}^{A}=\epsilon^{AB}D_{B}z_{(i)}$, writing them as the “curl” of a potential $z_{(i)}$. Here $\epsilon^{AB}$ and $D_{A}$ are the induced Levi-Civita tensor and covariant derivative of $\mathcal{H}$, respectively. Minimizing the average squared shear of $\phi^{A}$ over the surface then yields an eigenvector problem for the three $z_{(i)}$. As in Ref. [129], we fix the normalization of each $z_{(i)}$ by requiring the variance of each $z_{(i)}$ to agree with the same variance in the Kerr metric (cf. Eq. (A22) and the surrounding discussion in [129] for details). Finally, we define the spin magnitude $S$ as the Euclidean magnitude of a vector with three components of angular momentum found by inserting the three $\phi_{(i)}^{A}$ with smallest eigenvalues into Eq. (1):

[TABLE]

We define the “spin function” $\Omega\equiv\epsilon^{AB}D_{A}\omega_{B}$, where $\omega_{A}$ is the projection of the 1-form $K_{ij}s^{j}$ into the tangent space of $\mathcal{H}$. Then, we take the three moments of the spin function [51] to compute the direction of the spin,

[TABLE]

where $\vec{r}$ is the Euclidean position vector in the coordinate system of the simulation. We choose the normalization factor $N$ so that the Euclidean norm of $\hat{\chi}_{\Omega m}$ is 1. Finally, we define the full dimensionful spin vector as

[TABLE]

With the dimensionful spin in hand, we define the mass interior to $\mathcal{H}$ from the Christodoulou formula for (uncharged) Kerr BHs [161],

[TABLE]

where the irreducible mass $M_{\rm irr}$ depends only on the apparent horizon’s area:

[TABLE]

Even though the Christodoulou relation is, strictly speaking, only justified for stationary BHs, we also use it on dynamical AHs as a quasi-local mass.

From the dimensionful spin vector $\vec{S}$ given by Eq. (4) and the mass $M$ given by Eq. (5), we define the dimensionless spin vector

[TABLE]

We compute the magnitude of this vector $\vec{\chi}$ using the Euclidean norm. Equation (5) implies $0\leq|\vec{\chi}|<1$ for all of our black holes. For a discussion of how well this relation is satisfied in simulations of merging black holes with nearly extremal spins, see Ref. [41].

Finally, we define a coordinate center $\vec{x}$ for each AH as the surface-area weighted average of the location of the AH,

[TABLE]

This means that the (irreducible) mass dipole moment of the surface vanishes in a coordinate system centered on $\vec{x}$. In practice, we require that this condition be only approximately satisfied by truncating our representation of the shape of the AH at $\ell=1$ in an expansion in spherical harmonics, as shown in Eqs. (37)–(40) in [144].

2.2.2 Definition of relaxation time and reference time

for measuring initial BH quantities

The initial data do not perfectly describe two black holes in quasi-equilibrium. At the start of each simulation, therefore, the geometry relaxes to equilibrium on the dynamical time-scale of the individual BHs, changing the masses and spins of each BH by a fractional amount of order $10^{-5}$, and emitting a spurious pulse of gravitational radiation (often referred to as “junk radiation”). Our simulations do not attempt to resolve this pulse of short-wavelength gravitational radiation. Therefore, BH quantities (like their masses) fluctuate with an amplitude of about $\sim 10^{-5}$ for a few $100M_{0}$, before the fluctuations subside. Here, $M_{0}$ denotes the sum of the two Christodoulou masses at $t=0$. Subsequently, BH quantities vary on the inspiral time-scale.

To avoid the impact of junk radiation on our output quantities, we define a reference time $t_{\rm ref}>0$ as early in the simulation as possible,111In principle, one can choose any time to define the BBH parameters, since the dynamics provide a unique map from one time to the next. In practice, for analytical understanding and comparison with post-Newtonian theory, we want to choose the earliest possible time, where the post-Newtonian approximation is most accurate. but after the initial transients have decayed. We define a relaxation time $t_{\rm relax}>0$ at which we deem that transients have decayed, and then we set $t_{\rm ref}$ to be at least $t_{\rm relax}$. We extract the “initial” BH properties we ascribe to each black hole in our simulations at $t_{\rm ref}$. For most simulations in the catalog $t_{\rm ref}=t_{\rm relax}$, but we allow the two times to be different because they represent different concepts and because for some simulations (e.g. comparisons with waveform models or other NR codes) we desire to specify parameters at particular times. Since we do not attempt to resolve the initial transients, waveforms computed at different resolutions correspond to binaries with slightly different physical parameters. Because of our high precision, these small differences complicate our convergence testing, as discussed further in Sec. 4. We also recommend to only use the gravitational waveforms for retarded time $u>t_{\rm ref}$.

In practice, we compute $t_{\rm relax}$ as follows: We begin by defining a window of size $T_{\mathrm{window}}=300M_{0}$ and a time interval $\delta t=10M_{0}$. Considering the time series $M_{\rm irr,1}(t)$ of irreducible masses for the primary black hole,222In some earlier simulations, we used the timeseries of the magnitude of the dimensionless spin. For most simulations, these two quantities give comparable relaxation times. In practice, the areal mass tends to oscillate less than the spin magnitude. we compute its standard deviation $\sigma_{n}$ for sliding time windows $t\in[n\delta t,n\delta t+T_{\rm window}]$, $n=0,\ldots 30$. We compute a running average of the $\sigma_{n}$ using sequential sets of $10$ segments. As the junk radiation propagates away, the running average of $\sigma_{n}$ decreases with time. We identify the earliest time at which the running average stops decreasing, calling the center of this time window $t_{\text{relax,1}}$ (if this condition fails, we set $t_{\rm relax,1}=600M_{0}$). We repeat this calculation for the secondary BH, and take the larger of the two values as the final relaxation time, $t_{\rm relax}=\max(t_{\rm relax,1},t_{\rm relax,2})$.

2.2.3 Remnant masses, spins, and recoil velocities

After merger and ringdown, there remains a single remnant BH with its own mass, spin, and a recoil velocity (a “kick” caused by asymmetry in momentum carried by GWs).

As described in Sec. 2.2.1, our simulations extract at regular intervals the BH mass $M(t)$, spin $\vec{\chi}(t)$ and center $\vec{x}(t)$ from local quantities on each apparent horizon. This is also true after merger, when we extract the mass, spin, and center for the AH of the remnant from local quantities on the common apparent horizon. These quantities are highly dynamical immediately after merger because the remnant horizon is strongly, dynamically curved when it forms. After the common AH forms, it relaxes, quickly at first, and then rings down over a characteristic timescale determined by the remnant’s spin [162, 163]. After merger, we compute the remnant AH mass and spin on a dense (though not uniformly sampled) set of times. To compute the final AH mass and spin, we simply split these time samples into thirds, and choose the mass and spin averaged over the final third of the time samples. This approach seeks to mitigate small, residual time-dependent variations in the remnant mass and spin caused by numerical noise.

We also employ a simple procedure to estimate the coordinate recoil velocity (this turns out to be very close to the well-defined recoil velocity arising from gravitational-wave momentum flux integrals [164, 165]). We compute the coordinate center of the remnant AH as we do for the individual AHs during inspiral, using Eq. (8). Taking the last third of time samples of $\vec{x}(t)$, we model each of its components with a least-squares fit to a linear function of time. We then interpret the slopes of these fits as the coordinate velocities of the remnant BH.

2.3 Gravitational wave extraction

We extract the emitted gravitational waves from our simulations through two independent methods. The first computes the Newman-Penrose scalar $\Psi_{4}$ on a set of coordinate spheres centered at the initial (coordinate) center of mass of the two black holes. This is done by computing the Weyl tensor, projecting with a flat-space orthonormal null tetrad to form $\Psi_{4}$, and expanding in terms of spherical harmonics of spin weight $-2$. See Refs. [134, 166] for details. We do not use a properly orthonormalized null tetrad in computing $\Psi_{4}$ from the Weyl tensor, nor do we use anything other than coordinate-sphere extraction surfaces. Therefore, our computation of $\Psi_{4}$ at a finite radius differs from the standard definition by a multiplicative factor $1+\mathcal{O}(1/r)$. We eliminate these differences by extrapolating the waveforms to future null infinity (discussed in Sec. 2.4.1) to remove these and other near-zone effects. We also remove some artifacts of our choice of coordinates in the initial data via center-of-mass correction (discussed in Sec. 2.4.2).

Applications using the waveforms from the SXS catalog should use the extrapolated waveforms with center-of-mass corrections; however, for diagnostic purposes, we also make available the raw, finite-radius spherical-harmonic modes $r\Psi_{4}^{\ell,m}$ (available in the SXS catalog as files named rPsi4_FiniteRadii_CodeUnits.h5). In each case, the value of $\Psi_{4}$ can be evaluated at a point $(\theta,\phi)$ using

[TABLE]

where ${}_{-2}Y_{\ell,m}$ are the spin-weight $s=-2$ spherical harmonics, using the conventions given in Ref. [126].

The second gravitational-wave extraction method is independent of the first one. In this method, we compute the metric perturbation directly using Sarbach and Tiglio’s [167] formulation of the Regge-Wheeler and Zerilli equations [168, 169]. First, we compute the metric perturbation $\delta g_{ab}=g_{ab}-\eta_{ab}$ about a Minkowski background $\eta_{ab}$. From the variables of our first-order formulation of Einstein’s equations, we also read off the first time and spatial derivative of $\delta g_{ab}$. We evaluate $\delta g_{ab}$ and its derivatives on a set of coordinate spheres centered at the initial (coordinate) center of mass of the binary, and on each of these coordinate spheres we expand these quantities in terms of spin-weighted spherical harmonics. We then compute the spin-weighted spherical-harmonic modes of the Regge-Wheeler quantity $\Phi^{(-)}$ and the Zerilli quantity $\Phi^{(+)}$, which are combinations of the metric perturbations and their derivatives given in Eqs. (16–18), (22–29), and (A12–A21) of Ref. [170].333We use the opposite sign convention for $\Phi^{(+)}$ as does Ref. [170]; in particular, we replace the minus sign in front of Eq. (29) of Ref. [170] with a plus sign. Our sign convention agrees with that of Ref. [164], and ensures that a linearized wave in TT gauge satisfies Eqs. (52) and (53) of Appendix C, assuming Eqs. (10) and (11). Finally, we compute the modes of the strain using (cf. Eq. (83) of [171] and Eq. (4.34) of [164])

[TABLE]

Note that $\Phi^{(\pm)}$ are not gauge invariant in the general sense, but are sometimes referred to as such in the context of perturbation theory [172, 173, 167, 174]. Specifically, the definitions of $\Phi^{(\pm)}$ involve quantities that are invariant at first order under infinitesimal gauge transformations about a fixed background, assuming small metric perturbations from that background. However, gauge changes that cannot be treated as infinitesimal affect the waveforms at a significant level. In Sec. 2.4.2, we describe and remove some such gauge effects that are present in the waveforms even after extrapolation to future null infinity.

For diagnostic purposes, we make the modes $rh^{\ell,m}$ of the raw finite-radius quantities available in the SXS catalog as files rh_FiniteRadii_CodeUnits.h5. However, as with the $\Psi_{4}$ waveforms described above, the raw finite-radius quantities contain near-zone and gauge effects. For applications using the waveforms in the SXS catalog, one should instead use the version of these waveforms that are extrapolated to future null infinity and center-of-mass corrected (see Sec. 2.4 below). The strain at a point $(\theta,\phi)$ can be evaluated as

[TABLE]

Finally, note that each file in the SXS catalog describing raw, extrapolated, or center-of-mass corrected $h$ contains the real and imaginary parts of $h^{\ell,m}$, which differs by a minus sign in the imaginary component from the format described in Ref. [126]. More significantly, all waveforms in our catalog prior to this release have had an overall sign change in the definition of the strain. The strain as given in Eq. (11) reflects the current definition of the waveforms in our catalog, for more information see Appendix C. For details on checking the sign convention of the waveform files, see Appendix A.3.1.

2.4 Waveform post-processing

Our catalog contains waveforms that are extrapolated to future null infinity and corrected for center-of-mass motion. Here, we detail the extrapolation and center-of-mass corrections applied to each $h$ and $\Psi_{4}$ waveform in this catalog.

2.4.1 Extrapolation

During the evolution, we extract each waveform at a series of times on a set of concentric coordinate spheres surrounding the binary, decomposed in modes of spin-weighted spherical harmonic functions. We then extrapolate the waveforms to future null infinity, $\mathscr{I}^{+}$. Our method is similar to the one described previously in Refs. [175, 176], but we modify it to permit accurate extrapolation of precessing systems as follows. We transform the waveform modes into a corotating frame [177], in which the rotation is factored out, so that the corotating waveform modes vary slowly in time, even for precessing systems. We then simply extrapolate the real and imaginary parts of the corotating waveform modes. Ref. [175] shows that such slow time-varying behavior is crucial to convergence of the extrapolation process. We then transform the extrapolated result back into an inertial frame. Our previous extrapolation method [175, 176] decomposed the complex waveform modes into phase and amplitude, which results in slow temporal variations only for non-precessing binaries.

For each simulation, we compute both $h$ and $\Psi_{4}$ independently (Sec. 2.3). We extrapolate both quantities to $\mathscr{I}^{+}$ by the same method, though for simplicity we will only describe extrapolation of $h$ here. We extract the strain waveform as $h^{\ell,m}(T_{i},R_{j})$ on a grid of coordinate times $T_{i}$ (where $T=0$ is the start of the simulation) and on coordinate spheres of fixed radii $R_{j}$. The extraction radii are chosen between $R=R_{\rm min}$ and the outer boundary of the simulation domain, typically with about 24 extraction radii spaced uniformly in $1/R$. For most simulations in the catalog, we choose the innermost extraction radius $R_{\rm min}$ to be $100M_{0}$. For newer simulations, we chose $R_{\rm min}$ to be at least $\pi/\Omega_{0}$, where $\Omega_{0}$ is the initial orbital frequency; thus $R_{\rm min}$ is at least one gravitational-wave wavelength from the origin, and the waveform there is not dominated by near-field effects. We explicitly compute the areal radius of each sphere $r_{j}$, which depends on the time, by integration using the evolved metric. We also extract the average value of the metric component $g^{TT}$ over this sphere.444Each quantity used in the extrapolation is made available for download in files named rh_FiniteRadii_CodeUnits.h5 and rPsi4_FiniteRadii_CodeUnits.h5, documented in App. A.3. These quantities allow us to define the retarded time

[TABLE]

where $M_{\textrm{ADM}}$ is the ADM mass measured in the initial data. The second and third terms are essentially the familiar tortoise coordinate of the Schwarzschild solution, while the first term is a correction to the time coordinate. We choose this retarded time so that the one-form $du$ is approximately null with respect to the evolved metric [175].

Using these quantities, we can also express the waveform modes $h^{\ell,m}(u_{i,j},r_{i,j})$. Then, we interpolate the data to a common set of retarded times, $u_{k}$. We construct this set to be the largest subset of $u_{0,j}$ such that the waveform at each radius has known values at each time $u_{k}$. We also interpolate the areal radius of each sphere to the set of common times, so that the waveform can now be expressed as $h^{\ell,m}(u_{k},r_{k,j})$.

The next step is to rotate the waveform at each radius into a corotating frame [177]. To avoid the complication of extrapolating the transformation to the corotating frame itself, we simply choose the outermost extraction radius to define the corotating frame. We cannot expect the waveform at any other radius to be in precisely its own corotating frame, but in practice we still achieve our objective of ensuring that the waveform at each radius is slowly varying. We denote the waveforms in this frame as $\hat{h}^{\ell,m}(u_{k},r_{k,j})$.

Now, with slowly varying data tabulated on a common set of retarded times and a series of radii, we can extrapolate the waveform to infinite radius by approximating each mode with a polynomial of order $N$:

[TABLE]

At each time step, we choose the coefficients $\hat{h}^{\ell,m}_{(j)}$ to minimize the sum of the squared differences between the numerical data at that time and the polynomial value—real and imaginary parts being treated separately. The asymptotic waveform in the rotating frame is simply $\hat{h}^{\ell,m}_{(0)}(u)$. We then obtain the final asymptotic waveform by inverting the rotation that was applied above.

Note that the $m=0$ modes in SpEC waveforms appear to be generally unreliable, in the sense that they do not appear to converge with increasing extrapolation order or varying extraction radii, and they do not agree with CCE results [176]. This is true in the inertial frame for non-precessing systems, and typically true in co-precessing [66, 110, 178] or co-rotating [177] frames for precessing systems. This means that other modes in the inertial frame may be polluted by these inaccurate co-rotating frame $m=0$ modes for precessing systems.

The code that we use to perform this entire extrapolation procedure is available in the open-source python module GWFrames [179]. We provide all finite-radius waveforms, along with the extrapolated waveforms for $N=2$, $3$, and $4$, for both $h$ and $\Psi_{4}$. The extrapolated data files also contain a waveform from the outermost extraction radius, which is also given as a function of the corrected retarded time $u$ and scaled by the areal radius; we expect this will remove some (though not all) of the gauge artifacts present in raw waveforms extracted at finite radius. The waveform measured at the outermost extraction radius is different from the waveform computed using $N=0$ in Eq. (13), because the latter essentially averages the contributions from data extracted at smaller radii—which can be worse than doing no extrapolation at all. For this reason, we do not provide a waveform extrapolated using $N=0$, nor do we provide one for $N=1$ [175].

The choice of extrapolation order $N$ for a particular purpose must be informed by the behavior of near-field effects. If $\lambda$ is a typical wavelength present in a given mode, we expect the higher-order terms in the polynomial to scale not just as $1/r^{j}$ relative to the lowest-order term, but as $(\lambda/r)^{j}$ [175, 180]. Thus, as the binary spirals in toward merger and the length scales $\lambda$ become smaller, the polynomial will converge much more quickly with $N$. Thus, for example, if extrapolation with $N=4$ is required for accurate results early in the waveform, then $N=2$ may be sufficient closer to merger. On the other hand, using a large value of $N$ when the polynomial converges quickly can lead to overfitting of features that are poorly modeled as functions of retarded time and polynomials in $1/r$—so that $N=2$ extrapolation may actually be better than higher-order extrapolation during the merger and ringdown. The general rule of thumb, then, is to use higher-order extrapolation for analyses that require more accuracy during the inspiral, and lower-order extrapolation for analyses that require more accuracy during the merger and ringdown; there is no single choice that is best for all applications. In all cases, it is best to test the dependence of results on the extrapolation order by running an analysis multiple times using each of the various extrapolation orders. For the application of building waveform surrogate models [84, 181, and references therein], we use $N=2$ for the entire waveform.

2.4.2 Center-of-mass correction

The catalog provides extrapolated waveforms in two versions: one without, and one with, a correction for displacement and drift of the center of mass (COM). These COM-corrected waveforms have filenames ending in “CoM”. We recommend using COM-corrected waveforms for all applications of our catalog, but provide the original data for completeness and to allow for comparisons between COM-corrected and uncorrected waveforms.

These corrections are necessary because the waveform modes depend on the origin of coordinates used to define the spherical harmonics: if we move the origin, the modes will be mixed. Naively, we expect the origin to be centered on the binary, as that is the natural choice and is used in post-Newtonian models, for example. That choice leads to desirable features like relatively slowly varying mode amplitudes, and frequencies that are roughly proportional to the orbital frequency times the azimuthal number $m$ of the mode.555Spin-orbit coupling introduces effects for precessing systems where small additional components are present, and oscillate at frequencies roughly proportional to $(m\pm 1)$ times the orbital frequency. However, if the origin coincides roughly with the center of mass, these additional components can also be modeled to high accuracy without accounting for spurious mode mixing [110, 182]. However, all binary-black-hole systems simulated with SpEC contain essentially random offsets and drifts of the origin of coordinates relative to the COM, causing mode-mixing that manifests as irregular behavior in the waveform modes. To a good approximation, these irregularities can be thought of as direction-dependent time translations that appear uncorrelated between different physical systems or different numerical resolutions of the same system. These irregularities appear as essentially random contributions to waveform modes that are discontinuous with respect to changing physical parameters. These random and discontinuous effects would have to be modeled by surrogate, EOB, and phenomenological waveform models [83, 84, 85, 86, 183, 79, 61, 64, 65, 67, 66, 69, 70, 71, 72, 73, 75], or be optimized away in direct searches of detector data [80, 82, 76, 81]. By removing these effects, we simplify such analyses.

Figure 1 shows the translations and boosts of the COM for two example simulations (the non-precessing SXS:BBH:0314 and the precessing SXS:BBH:0627) from our catalog, each shown at multiple resolutions. We define the COM at each instant in time using the Newtonian definition,

[TABLE]

where $m_{1}\mbox{ and }m_{2}$ are the Christodoulou masses (as defined in Sec. 2.2.1) of the primary and secondary black holes, and $\vec{x}_{1}\mbox{ and }\vec{x}_{2}$ are the coordinate positions of the AH centers as defined in Eq. (8). In our convention, black hole “1” is the more massive of the two.

Two causes contribute to the COM motion visible in Fig. 1. First, while the initial data achieves $P_{\rm ADM}=0$ [133], the initial transients during relaxation to equilibrium may cause asymmetric GW emission, and thus impart a net linear momentum onto the binary. Second, since we do not attempt to resolve this junk radiation, we do not expect that the COM motion in our simulations will be convergent; indeed, we observe essentially randomly varying coordinate velocities of the COM for evolutions at different resolutions of the same initial data set.

As can be seen in the right column plots of Fig. 1, the total COM displacement is generally only a fraction of the total mass $M_{0}$ of the binary, and much smaller than the radii of the gravitational wave extraction spheres ($R_{j}\gtrsim 100M_{0}$, cf. Sec. 2.3). Nevertheless, these small COM displacements do have a noticeable impact on the higher-order modes of the computed gravitational radiation, as can be seen in the left column plots in Fig. 1. During inspiral and close to merger, the uncorrected higher-order waveform amplitudes oscillate—most notably the $(l,m)$ modes $(2,1)$, $(3,1)$, and $(3,3)$. This modulation is not expected on physical grounds; it is a gauge effect caused by mode mixing that follows from the COM displacement. The strongest effect of mode mixing is the leaking of power from the dominant $(2,\pm 2)$ modes into the subdominant modes, because we decompose our waveform using spin-weighted spherical harmonics centered on an offset, moving origin, which is unnatural.

We apply COM corrections to remedy these effects as a post-processing step. We compute the parameters necessary for the correction from the simulation and, as the corrections are BMS transformations [157, 184, 185, 186], they do not alter the physically meaningful aspects of the waveform. Deciding how to correct waveforms for center-of-mass motion is complicated and is described in a separate paper [187]. The final procedure itself is relatively straightforward, and we summarize it here.

We implement the measurement of and corrections to the COM using the open-source python module scri [188, 186]. While there is significant COM motion, as seen in Fig. 1, our COM correction deals only with the offset and drift—that is, the linear motion. Reference [187] presents investigations into potential physical contributions to the total COM motion, and the epicycles seen in Fig. 1.

We remove these gauge effects using translations and boosts. To re-center the simulations, we first measure the offset and drift of the COM and then retroactively apply the opposite motion to the waveform, to cancel out that motion. We define a time average of any quantity $Q(t)$,

[TABLE]

We choose $t_{i}$ to be the relaxation time for the simulation (defined in Sec. 2.2.2) and set $t_{f}=0.9t_{\rm CAH}$, where $t_{\rm CAH}$ indicates the time of common apparent horizon formation, so that merger and ringdown are not included in the average.

Now we would like to find a translation $\vec{\alpha}$ and boost $\vec{\beta}$ that give the best linear approximation to the motion of the measured $\vec{x}_{{\textrm{\tiny{COM}}}}$. Within the interval $[t_{i},t_{f}]$ we perform a linear least-squares fit to $\vec{x}_{\rm COM}(t)$ resulting in a best-fit motion $\vec{\alpha}+\vec{\beta}t$. As described in Appendix E of [186], the fit can be performed analytically, giving

[TABLE]

where the second equality of each line comes from $\langle t\rangle=\frac{1}{2}(t_{i}+t_{f})$, $\langle t^{2}\rangle=\frac{1}{3}(t_{i}^{2}+t_{i}t_{f}+t_{f}^{2})$.

We then apply a displacement to negate the linear motion of the COM given by $\vec{\alpha}+\vec{\beta}t$, computing this displacement separately for each resolution of each simulation. Reference [186] first showed that this method of COM corrections does indeed remove a large fraction of mode mixing and remedy the COM offset and drift. Reference [187] further confirms that applying the COM correction does improve all waveforms in the SXS simulation catalog, and introduces a robust and quantifiable method for this purpose.

Figure 2 shows translation, boost, and total displacement values for spin-aligned (top row) and precessing (bottom row) simulations in the catalog. More recent simulations in the catalog use an improved initial-data method [133] that achieves $P_{\rm ADM}=0$ in the initial data even for precessing systems and that reduces the overall displacement of the COM, especially for precessing cases. For most systems $\vec{\alpha}$ and $\vec{\beta}t_{\rm CAH}$ are comparable. Further details on the COM correction method and analysis can be found in [187].

3 Parameter space coverage

Expanding the catalog from the original 174 configurations to 2018 configurations has substantially improved our coverage of the BBH parameter space. Figure 3 shows the binary mass ratio $q=m_{1}/m_{2}\geq 1$ and the dimensionless spin magnitudes $|\chi_{1}|$ and $|\chi_{2}|$ for the simulations in our catalog. Each point in the scatter plots in Fig. 3 represents a simulation, while the histograms show the relative number of simulations with the given range of mass ratio and dimensionless spin magnitudes. The masses and spins plotted here are measured at the reference time, as discussed in Sec. 2.2.2. In the scatter plots, we see a substantial number of precessing simulations with mass ratios up to $q=4$ and $|\chi_{A}|\leq 0.8$, which were produced in order to construct the surrogate models of [83, 189, 181]. The subscript $A$ corresponds to the larger ($A=1$) and smaller ($A=2$) black holes. In addition, we show improved coverage of the nonprecessing subspace with mass ratios up to $q=8$ and $|\chi_{A}|\leq 0.8$. New simulations in this part of the parameter space were produced in order to construct the surrogate model of [84].

In contrast, there remain large regions that are unexplored in all BH merger catalogs, including ours. The projections in $q-|\chi_{A}|$ space in Fig. 3 show that while we have fairly dense coverage at low mass ratios, mass ratios larger than $q=4$ remain sparsely explored or completely unexplored. Similarly, aside from a few equal-mass, equal-aligned spin cases, the region of spin magnitudes above 0.8 remains almost completely unexplored. Few simulations exist with both high spins and high mass ratio. These are especially challenging, as they require high resolution and delicate control of the computational domain, including the shapes, sizes, and positions of the excised regions inside the black-hole horizons (see, e.g. [97]). For example, most simulations in the catalog with $q>4$ are non-precessing.

Figure 4 shows a histogram summarizing the number of orbits before merger in our simulations. Most simulations have between 10 and 30 orbits. This length is sufficient for many gravitational-wave applications, particularly when the systems have higher total masses, and thus remain in LIGO’s sensitive frequency band for fewer orbits. Spanning LIGO’s sensitive frequency band for binaries with lower total mass is more difficult. This can be done either by producing longer simulations, or through hybridizing numerical simulations by attaching the final orbits to an approximate post-Newtonian waveform to cover earlier times (see, e.g. [190] and the references therein). In both cases, achieving sufficient accuracy for applications to gravitational wave science remains challenging.

Figure 5 shows a histogram of the estimated initial orbital eccentricity for each of our simulations. For most simulations, we tune our initial data via an iterative procedure [134, 46, 47] to produce nearly quasicircular orbits; the orbital eccentricities of these simulations are almost all below $6\times 10^{-4}$. For some simulations we intentionally wish to study eccentricity so we omit the eccentricity-reduction step; these can be seen as the tail in Fig. 5. As described in Ref. [45, 46], we estimate orbital eccentricity by a least squares fit of an analytic function to the time derivative of the orbital frequency $d\Omega/dt$. This fit is performed during the $\sim 2$ orbits following $t_{\rm relax}$. It captures the monotonic inspiral-driven long-term trend in $d\Omega/dt$, and overlaid oscillatory variations caused by orbital eccentricity [see Eqs. (70) and (76) of [46]]. Motivated by Eq. (68) of [46], we report an eccentricity (effectively averaged over the first two to three orbits by our fit) computed with

[TABLE]

where $B_{\Omega}$ and $\omega_{\Omega}$ represent the amplitude and frequency of oscillations in $d\Omega/dt$, and $\Omega_{0}$ is the orbital frequency at time $t=0$. We fit to $d\Omega/dt$ rather than to $\Omega$ because the time derivative magnifies eccentricity-induced oscillations, making them easier to fit when eccentricity is small. At large separation, Eq. (17) reproduces the Newtonian definition of eccentricity to linear order in $e_{0}$. Since we neglect higher-order corrections in $e_{0}$, large values of $e_{0}$ reported for our simulations are only rough estimates of the actual orbital eccentricity. We have not made an effort to precisely recover any post-Newtonian eccentricities [191], which we expect to differ from Eq. (17) by fractional corrections of order $(v/c)^{2}$.

3.1 Coverage in spin space

Of the seven dimensions of parameter space for quasi-circular mergers, six are spin components. This high dimensionality makes the parameter space difficult to sample densely and uniformly. Previous catalogs [94, 93] have discussed coverage in spin space but without a thorough exploration of the degree of coverage.

Visualizing the coverage is difficult due to the high dimensionality, but we can focus on certain physically relevant combinations of the spin parameters. A commonly used spin combination that strongly affects total waveform phase is the effective spin [192, 193, 66],

[TABLE]

Here $\hat{L}$ is the direction of the instantaneous Newtonian orbital angular momentum, and we carry out the projection of the spins onto $\hat{L}$ using the Euclidean metric. When the effective spin is positive, the black holes merge more slowly, causing the gravitational-wave frequency to increase more slowly; conversely, when the effective spin is negative, the black holes merge more quickly, causing the gravitational-wave frequency to increase more quickly. Another benefit to considering the effective spin is that while the spin directions can precess in a complicated manner, $\chi_{\rm eff}$ is conserved up to at least the 2nd post-Newtonian order [194]. By contrast, the in-plane components

[TABLE]

are more relevant for recoil kicks [29, 28, 30, 195] and precession dynamics [196, 197].

Figure 6 is one view of the distribution of black-hole spins, measured at the reference time. Shown are the effective spin $\chi_{\rm eff}$ and the magnitudes of the in-plane vectors, $\chi_{A\perp}\equiv|\vec{\chi}_{A\perp}|$. The catalog contains a large number of non-spinning and aligned-spin simulations, leading to a spike at low in-plane spins. There is a population of simulations with $|\vec{\chi}_{A\perp}|\approx 0.8$, which were used to build the surrogate models of [85, 83].

A different view of the parameter space is relevant for understanding precession dynamics. Namely, a convenient combination of parameters is given by the two “tilt angles” $\theta_{A\,L}$ [197],

[TABLE]

and the in-plane angle $\Delta\Phi$ between the two $\vec{\chi}_{A\perp}$ vectors,

[TABLE]

with the sign determined according to $\operatorname{sgn}\Delta\Phi\equiv\operatorname{sgn}\{\hat{L}\cdot[(\chi_{1}\times\hat{L})\cdot(\chi_{2}\times\hat{L})]\}$ [197]. Figure 7 plots these three parameters. Because of the large number of aligned-spin simulations, there is pileup at values of $\cos\theta_{A\,L}=\pm 1$. There are also larger number of simulations with purely in-plane spins, leading to another pileup at values $\cos\theta_{A\,L}=0$. The distribution in $\Delta\Phi$ is relatively flat.

3.2 Parameter space coverage and LIGO measurements

Figure 8 compares the parameter space coverage of the SXS catalog to selected astrophysical measurements of coalescing black holes by the LIGO detectors, obtained from [198, 199]. The left panel illustrates the magnitudes of the dimensionless spin vectors of the two black holes and their corresponding spin tilt angles $\theta_{AL}$ for each simulation in the catalog as measured at the reference time (dots, color coded by the simulation’s mass ratio). These are plotted over the marginalized posterior distributions for these quantities for GW151226 [8], one of only two observed black hole binaries, to date, with evidence for non-zero spin [12] (denoted by greyscale pixels). The right panel shows the effective spin and the mass ratio of each system in the catalog and, for GW151226 and three additional gravitational wave detections, the 90% credible contours of the marginalized 2-dimensional posterior distributions.

In both cases, the SXS catalog covers a large part of the relevant parameter space that is consistent with the LIGO measurements. In particular, LIGO observations point towards black-hole binaries with small effective spins, a region probed well with the majority of the SXS simulations. However, the comparison also suggests that there are regions of the parameter space where the coverage is sparse, especially for unequal mass systems. Future simulations will help fill in this region of the parameter space and serve to improve the accuracy of waveform models.

4 Waveform quality

As discussed in Sec 2.1, during a simulation SpEC employs dynamical adaptive $p$- and $h$-refinement, adding or subtracting grid points, or splitting or joining subdomains, according to local measures of truncation error and a globally specified AMR tolerance. The vast majority of the simulations in the catalog presented here have been run at multiple levels of this AMR tolerance (which we will henceforth call “multiple resolutions” for brevity). As discussed in Sec 2.1, because of the adaptivity, a small number of our simulations do not show convergence with resolution; these are noted in Figs. 9 and 10.

It is not always straightforward how to compare two waveforms. Some comparison criteria, such as phase or amplitude differences for a particular spherical-harmonic waveform mode, might be important for certain applications and not for others. Following [200, 85], here we take an approach motivated by the practical application of the waveforms to gravitational wave science. We define an overlap $\mathcal{O}$ between waveforms and examine overlaps between waveforms from simulations with the same initial parameters but differing resolutions. This provides one measure for the accuracy of our waveforms. A second accuracy measure involves comparing waveforms extracted using the extrapolated metric perturbation $h$ to those extracted at the same resolution with the extrapolated Weyl scalar $\Psi_{4}$.

To define the overlap $\mathcal{O}$, we employ a standard inner product $\langle h^{1},h^{2}\rangle$ between complex waveforms $h^{1}$ and $h^{2}$ given by

[TABLE]

where the tilde denotes a frequency domain signal. This product is real if $h^{1}(t)$ and $h^{2}(t)$ are real. The quantity $S_{n}(|f|)$ is the noise power spectral density, which quantifies the spectrum of the colored Gaussian noise, and is also used to whiten the signals in a gravitational wave detector. Here, we assume a flat noise spectrum, and set $S_{n}(|f|)=1$. This choice has the benefit that the results described here are independent of the total mass of the binaries. These overlaps are in qualitative agreement with those obtained using the Advanced LIGO noise curve, when the total mass of the system is in the range $50-100M_{\odot}$, as shown in Fig. 4 of Ref. [181]. Details of our implementation of the inner product are given in App. B.

Given this inner product, we define an overlap that accounts for the information in both polarizations $h_{+}$ and $h_{\times}$ in $h(t,\theta,\phi)=h_{+}-ih_{\times}$, as if measured by two ideally oriented detectors located at a given $(\theta,\phi)$ position in the source frame,

[TABLE]

For waveforms that are identical up to a re-scaling, $\mathcal{O}=1$, otherwise $\mathcal{O}<1$. We are interested in quantifying the difference in numerical simulations of the identical physical system at different numerical resolutions, so $h^{1}$ and $h^{2}$ correspond to evolutions of the same initial data set, but with different numerical resolution. As described in detail in Appendix B, when comparing two such waveforms $h(t,\theta,\phi)$ we allow for an overall rotation $\delta\phi$ of one relative to the other, and a time offset $\delta t$. For each pair of waveforms and direction $(\theta,\phi)$ from the binary to the detector, this results in a mismatch

[TABLE]

For every configuration we include in Figures 9 and 10, we evaluate the mismatch at 20 distinct source frame directions $(\theta,\phi)$.

Figure 9 shows a histogram of the resulting mismatches between the two highest resolution simulations for the 1872 simulations in the catalog with more than one resolution. The majority of the simulations, 1777, have more than two resolutions and in these cases we can assess the convergence of the waveforms by comparing the mismatch between the two highest resolutions with that between the second and third highest resolutions. We expect the former mismatch to be smaller than the latter mismatch in waveforms that converge with increasing resolution settings, and when this fails to occur we label the waveform “nonconvergent.”

The top panel of Fig. 9 depicts the mismatches for the extrapolated metric perturbation $h$ computed using the Regge-Wheeler-Zerilli extraction technique described in Sec. 2.3. The bottom panel depicts the mismatches between the Weyl scalar $\Psi_{4}$, weighted by $f^{-2}$ in order to give the same frequency weighting as $h$. We see that overall the mismatches are small, with the mismatches broadly lying between $\sim 10^{-6}$ and $10^{-3}$, appropriate for many current applications to gravitational wave science. This is true for both the mismatches computed using $h$ and $\Psi_{4}f^{-2}$. There are 37440 mismatches plotted in each panel of Fig. 9, corresponding to 20 detector directions $(\theta,\phi)$ for 1872 simulations; only 3558 of these mismatches are nonconvergent as determined by $h$ and 4851 as determined by $\Psi_{4}f^{-2}$.

The error estimates in Fig. 9 have a tail that extends to rather high mismatches. This feature was investigated in App. B of Ref. [181], and shown to be a consequence of unresolved initial transients. As mentioned in Sec. 2.2.2, these initial transients cause simulations at different resolutions to correspond to binaries with slightly different physical parameters. In particular, we find that the in-plane spins of different resolutions can be inconsistent with each other. Ref. [181] used surrogate-modeling tools to show that the high-$\mathcal{M}$ tail in Fig. 9 is dominated by the small differences in system parameters. By training a surrogate model on high-resolution simulations, and evaluating it with spins of medium-resolution simulations, Ref. [181] found mismatches (between the surrogate and medium-resolution NR waveforms) to always be below $\lesssim 10^{-2}$. Thus our error estimate is overly conservative and does not reflect the actual truncation error of the simulations. We expect that the actual truncation error tail in Fig. 9 should only extend to $\lesssim 10^{-2}$ rather than $\sim 10^{-1}$.

Furthermore, the mismatches in Fig. 9 are somewhat pessimistic measures of the accuracy of the waveforms, since they actually measure the error in the second highest resolution, not the highest. Given the refinement scheme in SpEC, we cannot use extrapolation of the convergent waveforms to provide an error measurement on the highest resolution simulation. The mismatches also combine all possible sources of error together, so that we cannot distinguish truncation error from other sources of error, such as inaccuracies in our prescription for extrapolating the waveforms to infinity.

Figure 10 provides another perspective on the accuracy of the waveforms in the catalog. The top panel provides a histogram of the mismatches between $\ddot{h}f^{-2}$ and $\Psi_{4}f^{-2}$ from for the highest resolution for each waveform. This mismatch is between two quantities that are equal in theory, and so it provides an independent assessment of the numerical accuracy of the waveform, and especially of the two waveform extraction techniques. The factor $f^{-2}$ ensures the same frequency weighting for the signals used in the mismatches of both Figs. 9 and 10. The overall level of the mismatches is lower here than in Fig. 9, with fewer cases extending beyond a mismatch of $10^{-3}$. Instead of comparing $\ddot{h}$ vs. $\Psi_{4}$, one could perform mismatches of $h$ against $\iint\Psi_{4}$, but by using the time derivative $\ddot{h}$ we avoid the difficulty of needing to fix two integration constants, which is not straightforward [201].

The bottom panel of Fig. 10 is a scatter plot of the mismatch between $\ddot{h}f^{-2}$ and $\Psi_{4}f^{-2}$ at the highest resolution and the same mismatch at the second-highest resolution. The dashed line has unit slope and helps to quickly assess which mismatch is larger. While we see a few more points above the line, the plot indicates that the mismatch between $\ddot{h}f^{-2}$ and $\Psi_{4}f^{-2}$ is roughly independent of resolution, but the scatter is wide with many outliers. Furthermore, the histogram in the top panel of Fig. 10 is visually unchanged when computed with the second-highest resolution rather than the highest resolution (as is shown). Therefore Fig. 10 shows an additional source of error, perhaps caused by differences in waveform extraction and extrapolation, which is independent of numerical resolution and is smaller on average than numerical truncation error. This also suggests numerical truncation error affects the strong-field evolution more than it does the wave propagation in the far zone; otherwise the different wave extraction methods would show larger differences with numerical resolution.

All waveforms in Figs. 9 and 10 have been extrapolated to infinity as described in Secs 2.4.1. To quantify errors in our extrapolation procedure, Fig. 11 shows mismatches between waveforms that are identical except for details of extrapolation. In particular, we compare the highest-resolution waveform from each of our simulations, extrapolated to infinity using extrapolation order $N=2$ (the standard choice we make, e.g., in surrogate models [181]), versus the same waveform extrapolated to infinity using extrapolation order $N=3$. The mismatches peak below $\mathcal{M}\sim 10^{-6}$ with a tail that extends to $10^{-3}$, demonstrating that by this measure, errors in our extrapolation procedure are on average unimportant compared to numerical truncation error.

A sufficient condition for two waveform models to be indistinguishable is [202, 203, 204, 205]

[TABLE]

where $\mathcal{M}$ is the mismatch [cf. Eq. (24)] and $\rho$ is the signal-to-noise ratio (SNR) of the observation the models are describing. Here $D$ is the number of relevant model parameters, with $D=8$ for spin-precessing systems. For $\rho=24$ (the case for GW150914 [3]), this corresponds to a mismatch of $7\times 10^{-3}$. Comparing to Fig. 9, we find that for most of our simulations, different numerical resolutions are indistinguishable in this sense, using a flat noise curve. This gives us confidence that our numerical waveforms are well suited for interpreting gravitational-wave observations as loud as GW150914, provided that the total mass is sufficiently high that the waveform is long enough to span the observed signal. However, note that using Eq. (25) to determine whether two waveforms are indistinguishable in a particular detector would require using that detector’s noise curve in the mismatch calculations and choosing a total mass for each mismatch.

5 Remnant properties

The remnant properties of the final black hole are of great interest both astrophysically and for constructing semi-analytical waveform models. On the astrophysical side, the final mass, spin and kick velocity give important information about the possible progenitors of the system and may help distinguish BBH formation channels [206]. In waveform modelling, the final mass and spins are important ingredients in constructing the full waveform, determining the quasi-normal modes whose superposition creates the ringdown signal. This requires one to connect the parameters of the black holes far from merger to those of the remnant. Thus the task of inferring remnant properties has received a lot of attention [207, 28, 29, 208, 209, 210, 211, 212, 91, 213, 214, 215, 216, 217, 195, 218, 219, 220, 197, 221, 222, 223, 189].

In this section, we compare fits for the final mass and spin magnitude from the literature to the simulations in the catalog, restricting our attention to cases where the measured eccentricity is $\leq 2\times 10^{-3}$. We make use of the publicly available implementation of the fits in LAL [224, 225]. To estimate the error in NR data, we use the difference between the highest and second highest resolution, where more than one resolution is available. We define the errors in the fitted mass as $\Delta m=m_{{\textrm{\tiny{NR}}}}-m_{{\textrm{\tiny{fit}}}}$ and similarly for the magnitudes of the final spin.

We begin by considering the fits for the remnant mass from Healy and Lousto [222] (HL2016) and Jiménez-Forteza et al. [221] (UIB2016). As can be seen in Fig. 12, we find good agreement between fits and the NR simulations. The errors in the fits are below $|\Delta m|\lesssim 0.004M$ for 90% of the cases, more than an order of magnitude larger than the NR errors. There is also a tail that extends to larger negative errors which shows that the fits systematically overestimate the final mass.

We next consider the final dimensionless spin magnitudes, using the models from Hofmann et al. [220] (HBR2016) in addition to HL2016 and UIB2016, as shown in Fig. 13. For precessing cases we follow current LIGO/Virgo analyses and evolve the spins (using the 3.5 post-Newtonian spin evolution equations) from the relaxed time to the Schwarzschild innermost stable circular orbit (ISCO). Projections of these spins along the Newtonian orbital angular momentum direction at ISCO are used as inputs for the remnant mass and spin fits. The HL2016 and UIB2016 spin fits are augmented by the sum of the in-plane spin components at ISCO (for details see [226]). For all models, we find errors of order $\Delta|\chi|\lesssim 0.01$ for 90% of the cases, more than an order of magnitude larger than the NR errors. While there is little difference in the magnitude of the errors between different models, we find that the HBR2016 model shows the least skew around 0, most likely because of special correction factors included in the model to handle precessing cases.

More accurate fits [189, 181] for the remnant mass, spin and recoil kick velocity have recently been developed by training directly against some of these simulations. The errors in these fits are comparable to the NR errors, but the fits have been trained only against simulations with $q\leq 4$, $\chi_{1},\chi_{2}\leq 0.8$. However, they are shown to extrapolate reasonably to higher mass ratios and spins in [189, 181]. Note that when applying aligned-spin fits to precessing systems, there is an ambiguity as to what time or frequency the precessing spins are to be evaluated when using the aligned-spin model. The fits in [189, 181] resolve this ambiguity by training directly against precessing simulations. Ref. [189] also suggests that fits built only from aligned-spin NR simulations can become inadequate at SNRs $\sim 5$ times that of GW150914, so there is a need to calibrate directly to precessing simulations.

6 Conclusion

In this paper, we have presented a substantial expansion of the SXS catalog of numerical-relativity simulations of black-hole binaries, which is publicly available at https://www.black-holes.org/waveforms [92]. Our catalog now includes 2018 simulations (including 1426 that are precessing) with a median length of 39 cycles of $\ell=m=2$ gravitational waves. We have considerably expanded our coverage of the parameter space, especially for mass ratios up to 4 and spin magnitudes up to 0.8. While our catalog does include simulations with mass ratios up to 10 and simulations spin magnitudes up to 0.998, the parameter space of high mass ratios, high spins, or both remains highly challenging and largely unexplored both by our catalog and by other catalogs. The remnant masses and spins agree well with existing fits in the literature, although differences between the fits and our simulations are larger than differences between our different resolutions; as a result, improved fits have recently been constructed directly from NR simulations [189].

We also have assessed the quality of our numerical waveforms. We find that mismatches between waveforms at different numerical resolutions are smaller than $10^{-3}$ for the vast majority of simulations in our catalog, although a few simulations have larger mismatches between different resolutions. As discussed in Sec. 4, current simulations are adequate for parameter estimation with Advanced LIGO and Virgo signals detected at the current level of sensitivity. Significantly louder observations, possible with future ground- and space-based detectors, might require numerical waveforms with significantly higher accuracy.

In the future, we will continue to expand our catalog, with the ultimate goal of fully covering the parameter space of binary black holes. This will likely require novel approaches to enable simulations with both high mass ratios and nearly extremal, precessing spins. We will also work toward longer simulations, which span detectors’ frequency bands down to smaller total masses, and more thoroughly explore cases with higher eccentricity. We will develop improved initial data with less spurious “junk” gravitational radiation (e.g. [227]). Finally, we are working towards computing gravitational waves using Cauchy-characteristic extraction [228, 229, 176, 230, 231] rather than extrapolation.

We are pleased to thank Davide Gerosa and Josh Smith for helpful discussions. This work was supported in part by the Sherman Fairchild Foundation, NSF grants PHY-1708212 and PHY-1708213 at Caltech, NSF grants PHY-1606654 and DGE-1650441 at Cornell, NSF grants PHY-1606522, PHY-1654359, and AST-1559694 and by Dan Black and Goodhue-McWilliams at Cal State Fullerton. S. E. Field is partially supported by NSF grant PHY-1806665. H. Fong and C. J. Woodford acknowledge support from NSERC of Canada grant PGSD3-504366-2017, and H. Fong further acknowledges support from the University of Toronto and the Japan Society for the Promotion of Science. P. Schmidt acknowledges support from the NWO Veni grant no. 680-47-460. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. Computations were performed on the GPC supercomputer at the SciNet HPC Consortium [232]; SciNet is funded by: the Canada Foundation for Innovation (CFI) under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund (ORF) – Research Excellence; and the University of Toronto. Computations were performed on the supercomputer Briareé from the Université de Montréal, managed by Calcul Québec and Compute Canada. The operation of these supercomputers is funded by the Canada Foundation for Innovation (CFI), NanoQuébec, RMGA and the Fonds de recherche du Québec-Nature et Technologie (FRQ-NT). Computations were performed on the Wheeler cluster at Caltech, which is supported by the Sherman Fairchild Foundation and by Caltech; on NSF/NCSA Blue Waters under allocation NSF PRAC–1713694; and on XSEDE resources Bridges at the Pittsburgh Supercomputing Center, Comet at the San Diego Supercomputer Center, and Stampede and Stampede2 at the Texas Advanced Computing Center, through allocation TG-PHY990007N. Computations were performed on the Minerva high-performance computer cluster at the Max Planck Institute for Gravitational Physics in Potsdam.

Appendix A Contents of the SXS catalog

The SXS catalog contains of a set of binary black hole simulations, each labeled by an identification string of the form “SXS:BBH:dddd”, where “dddd” is a four-digit number. A given catalog entry labeled by “SXS:BBH:dddd” usually contains data from several resolutions, i.e., several simulations with identical code and identical parameters except for the AMR tolerance (cf. Sec. 2.1). These resolutions are labeled “LevN”, where “N” is an integer that increases with finer resolution. The resolution labels for different catalog entries are not necessarily comparable. Each resolution starts from identical initial data, sampled onto the appropriate grid for that resolution.

The entries “SXS:BBH:dddd” are labeled roughly (but not exactly) chronologically, so that larger numbers usually correspond to later simulations with more recent versions of SpEC. Sometimes two different catalog entries have identical physical parameters (black hole masses and spins), but are separate entries because they follow a different number of orbits, have different orbital eccentricities or outer-boundary radius, or because they use different prescriptions for initial data or gauge conditions. In addition, a number of simulations have been repeated using identical or nearly identical parameters with an improved version of SpEC, including the aligned-spin simulations of [200]. Once an entry is in the catalog, it is never replaced by a newer simulation with the same “SXS:BBH:dddd” label; instead, the newer simulation is assigned a new label. An existing entry is modified only rarely, in the case of simple corrections that can be done by postprocessing: for example, we found that some of our waveforms in the catalog had the opposite overall sign convention as intended, so the offending waveforms were changed appropriately in the catalog and the version numbers of the modified files were updated.

A.1 Available data for each simulation

Each simulation in the catalog contains a metadata.json file that lists physical and code parameters, derived quantities such as remnant properties, and informational fields. This file is documented in Table LABEL:metadata-fields of Sec. A.2.

Each simulation also includes files that contain the spherical-harmonic modes of gravitational waveforms, as described in Sec. A.3 and listed in Table 3. Finally, each simulation includes files Horizons.h5 containing masses, spins, and other properties of the apparent horizons.

All dimensionful quantities are given in arbitrary units. All vector and tensor quantities are given in an asymptotically Cartesian coordinate system in which the black holes at the initial time are along the $x$-axis (BH 1 at positive $x$), and the initial Newtonian orbital angular momentum is along the $z$-axis.

A.2 Metadata format

A.3 Format of the HDF5 data files

A.3.1 Waveform files

We provide several files in HDF5 format that contain extracted gravitational waves. Users should generally use rhOverM_Asymptotic_GeometricUnits_CoM.h5 or rMPsi4_Asymptotic_GeometricUnits_CoM.h5, which contain our best-effort gravitational wave modes, extrapolated to future null infinity (cf. Sec. 2.4.1) and corrected for center-of-mass effects (cf. Sec. 2.4.2). These files contain several groups (i.e., the HDF5-equivalent of folders), each one of them containing a separate set of GW waveform modes, corresponding to different extrapolation order. The groups are named Extrapolated_N<int>.dir, where the integer <int> indicates the polynomial order of extrapolation. See the discussion in Sec. 2.4.1 regarding how to choose extrapolation order. In addition, there is a group OutermostExtraction.dir which contains the GW modes at the largest available extraction radius, without extrapolation, but scaled and corrected as described in Sec. 2.4.1 for consistency with the extrapolated waveforms.

Each of these HDF5 groups contains one HDF5 dataset for each $(\ell,m)$ mode; this dataset is named Y_l<int1>_m<int2>.dat. The first integer <int1> indicates the value of $\ell$ for this particular mode, and the second integer <int2> indicates the value for $m$. <int1> is always $\geq 2$, whereas <int2> takes negative values for negative $m$. Each Y_l<int1>_m<int2>.dat dataset contains three columns—either

[TABLE]

Here $u/M$ is the retarded time defined in Eq. (12), made dimensionless by division with the sum of the two Christodoulou masses, $M=m_{1}+m_{2}$, where each mass is measured at the reference time via Eq. (5). The time spacing is non-uniform, with more points in regions of higher GW frequency. Note that the $h$ files contain the real and imaginary parts of $h^{\ell,m}$ as opposed to the polarizations $h_{+}$ and $h_{\times}$; see Eq. (11).

In addition to the primary files rhOverM_Asymptotic_GeometricUnits_CoM.h5 and rMPsi4_Asymptotic_GeometricUnits_CoM.h5, some files with intermediate GW data are provided. These files contain the GW modes without COM-correction and/or without extrapolation as listed in Table 3. These files may sometimes be useful for debugging purposes, but they should not generally be used. Some of these extra waveform-files contain extrapolated GW modes (extrap=yes) in the same structure and format as just described for rhOverM_Asymptotic_GeometricUnits_CoM.h5 and rMPsi4_Asymptotic_GeometricUnits_CoM.h5. The files with non-extrapolated waveforms (extrap=no in Table 3) contain groups for different extraction radii, named Rxxxx.dir, where the four-digit integer xxxx indicates the radius of the extraction sphere. SpEC chooses extraction radii always at integer values, so there is no rounding in this number. However, note that the radius is given in dimensionful code units without division by $M$.

Some of the HDF5 files will have an HDF5 dataset called VersionHist.ver in the root HDF5 group, which stores the entire version history of the file. If a file does not have this dataset, then it is on version 0. VersionHist.ver is an array of pairs where the first element in the pair is the git commit id for the parent of the commit responsible for the change, and the second element is a description of the change from the previous version.

Only version 1 of rh*.h5 files follows the sign convention for the strain $h$ in Eqs. (52,53). For version 0, there is an overall minus sign in Eqs. (52,53), and hence an overall sign difference between the waveforms contained in version 0 and version 1 of the rh*.h5 files. Notice that this also implies that relation between $\ddot{h}$ and $\Psi_{4}$ in Eq. (55) is off by a sign for rh*.h5 files on version 0. We make our best effort to ensure each type of data file is on the same version across all runs, but we still recommended checking the version when working with files across different runs.

A.3.2 Apparent horizon files

For each simulation, we provide one file, Horizons.h5, containing data computed from the apparent horizons. Each file contains 3 groups named AhA.dir, AhB.dir, and AhC.dir; these groups correspond to the two individual horizons (labeled “A” and “B”) and the common apparent horizon (labeled “C”). Typically, horizon A corresponds to the black hole with the larger initial Christodoulou mass.

Table 4 lists the seven HDF5 datasets given for each horizon (i.e., for each HDF5 group). The first column of each data set contain the time $T$ from the simulation, and the next column (for scalar quantities) or next three columns (for spatial vector quantities) contain the data. Quantities with dimension are each given in the same arbitrary units, and all vector quantities are given in the (asymptotically) inertial frame of the simulation. Note that the times $T$ are not spaced uniformly.

Appendix B Computation of mismatches

One way to assess the quality of numerical waveforms is to compute mismatches between waveforms that are supposed to be equal but are computed using different numerical resolution parameters, boundary conditions, extraction procedures, or methods of extrapolating to infinity. Our method of computing mismatches is similar but not identical to the procedure described in Appendix D of [85].

We begin with two waveforms in the frame of our simulation

[TABLE]

where the complex waveforms include both gravitational-wave polarizations, and the angular dependence of the waveforms is usually written as a sum of spin-weighted spherical-harmonic modes. We define $t=0$ as the time of maximum power in the waveform, and we truncate the first $500M$ of the waveform to eliminate effects of initial transients sometimes known as “junk radiation.” This $500M$ is uniform across all waveforms, distinct from the reference time discussed in Sec. 2.2.2, which can vary across resolutions.

If two waveforms $h^{1}$ and $h^{2}$ differ only by an overall coordinate rotation or an overall time shift, we would like the two waveforms to compare as equal, and therefore to have an overlap of unity. We accomplish this through two steps. First, before we compute the mismatch, we rotate both waveforms so that the orbital angular momentum lies along the $+z$ axis at some fiducial time $t_{0}$. In other words, at $t=t_{0}$ the coordinate frame is momentarily aligned with the minimally rotating coprecessing frame of Ref.[178]. We choose $t_{0}=t_{\rm begin}+1000M$, where $t_{\rm begin}$ is the earliest time that is covered by both waveforms, after the above-described truncation of the first $500M$ of each waveform (if we are comparing three waveforms, as is the case in Fig. 9 where we compute overlaps between three numerical resolutions, we choose $t_{\rm begin}$ to be the earliest time that is covered by all three waveforms). Second, after the rotation, we allow one of the waveforms (choose it to be $h^{2}$) to have an arbitrary azimuthal angle shift and an arbitrary time shift: $h^{2}(t,\theta,\phi)\to h^{2}(t+\delta t,\theta,\phi+\delta\phi)$. The shift $\delta\phi$ corresponds to a redefinition of the position of the two black holes in the orbit at $t=t_{0}$. Later we will minimize the mismatch over $\delta t$ and $\delta\phi$.

When evaluating the accuracy of our numerical waveforms, we do not wish to ignore the polarization information contained in these waveforms. Furthermore, while we are interested in the angular dependence of the waveforms, we do not want to concern ourselves with antenna patterns of detectors. Therefore, we compute overlaps assuming the most optimistic detector scenario: an ideal network of two detectors located at $(\theta,\phi)$ relative to our source frame and oriented normal to the direction of wave propagation, one detector measuring $h_{+}(t)$ and the other measuring $h_{\times}(t)$. This motivates the two-detector overlap defined in Eq. (23) (see also Appendix D of [85]).

To compute the Fourier transforms in Eq. (22) we use an FFT after tapering the ends of the time-domain waveforms. For the window function we use a Planck-taper window (Eq. (7) of [234]). This function depends on four parameters $t_{1},t_{2},t_{3},t_{4}$: it rises smoothly from zero at $t=t_{1}$ to unity at $t=t_{2}$, and falls smoothly from unity at $t=t_{3}$ to zero at $t=t_{4}$. We choose $t_{1}=t_{\rm begin}$ and $t_{4}=t_{\rm end}$, where $t_{\rm end}$ is the latest time that is covered by both waveforms (or all three, in the case of Fig. 9). We choose $t_{2}$ to be the time of the 10th zero-crossing of the real part of the $(2,2)$ mode after $t=t_{1}$, and we choose $t_{3}$ to be $50M$ after the peak of the waveform. Before we compute the transforms, we pad with zeros and interpolate each resolution’s time-domain waveform onto the same evenly spaced set of time samples, where the number of samples is chosen to be

[TABLE]

We then truncate the Fourier transforms at a low-frequency cutoff $f_{\rm min}$ chosen to be twice the waveform angular velocity (as defined by Ref. [177]) at $t=t_{2}$, and a high-frequency cutoff $f_{\rm max}$ chosen to be 16 times the waveform angular velocity at the time of peak waveform power; the extra factor of 8 is chosen to resolve up to $m=8$ spherical-harmonic modes, with an extra margin of a factor of 2.

The optimization over $\delta t$ can be simplified by noting that the Fourier transform of $h(t+\delta t)$ is $\tilde{h}(f)e^{2\pi i\delta t}$, so the overlap takes the form

[TABLE]

To compute $\max_{\delta t}\mathcal{O}(\delta\phi,\delta t)$ for a fixed $\delta\phi$, we evaluate the integral in Eq. (31) efficiently for many values of $\delta t$ simultaneously using an inverse FFT, and we take the maximum value. We then use standard numerical maximization techniques to maximize over $\delta\phi$, resulting in the mismatch $\mathcal{M}$ defined in Equation (24).

In order to include the effect of higher-order spherical-harmonic modes, we evaluate the mismatch at 20 points on the unit sphere, evenly spaced in $\cos\theta$ and $\phi$, that describe the direction of the detector with respect to the source. The mismatch computed using each of those 20 directions is plotted separately in Figs. 9 and 10.

Note that some mismatch computations also explicitly minimize over a polarization-angle shift $\psi$, which rotates the polarization tensor that we use to decompose the waveform into the two polarizations $h_{+}$ and $h_{\times}$. On the $z$ axis, optimization over $\psi$ is precisely degenerate with optimization over $\delta\phi$, even when all modes are included [235]; off the $z$ axis this degeneracy is broken. Here we consider $h^{1}$ and $h^{2}$ to be different even if they differ only by a polarization-angle shift, since we are considering the case of a detector network that measures both polarizations, and since our numerical waveforms contain polarization information. Hence we do not minimize over a polarization-angle shift when computing overlaps and mismatches.

Appendix C Sign conventions

With so many sign conventions in the literature, we explicitly provide an outline of sign conventions used in SpEC. Here, Greek indices represent four-dimensional spacetime coordinate indices, and Latin indices represent three-dimensional coordinate indices for a space-like hypersurface. For a spacetime metric $\psi_{\mu\nu}$ with signature $(-,+,+,+)$, we foliate the spacetime into space-like slices orthogonal to a timelike unit one-form $t_{\mu}$,

[TABLE]

where $t$ is a scalar function representing global time, and $N$ is the lapse. With a shift vector $N^{i}$, we define the spatial metric and the extrinsic curvature, respectively,

[TABLE]

where $K_{ij}$ represents the spatial components of the extrinsic curvature, and the subscript [math] indicates the time component. While the sign convention for $K_{ij}$ is rather varied in the literature, the one chosen here follows the Misner-Thorne-Wheeler convention and is found in many prominent texts [236, 237, 20, 238, 239]. There are also several texts that follow the opposite sign convention [240, 241, 242].

We define the 4-volume form $\epsilon_{\alpha\beta\gamma\delta}$ and the 3-volume form $\epsilon_{ijk}$ on the spatial slices as follows,

[TABLE]

and all others related by complete antisymmetry, $\epsilon_{abcd}=\epsilon_{[abcd]}$. Note that some texts define the 3-volume form as $\varepsilon_{ijk}=\epsilon_{ijk\mu}t^{\mu}$, which incurs a minus sign relative to the definition in Eq. (38), so that $\varepsilon_{ijk}=-\epsilon_{ijk}$.

We define the Christoffel symbols and the Riemann, Ricci, and Weyl tensors, respectively, following the Misner-Thorne-Wheeler convention,

[TABLE]

The Ricci tensor of the spatial slices is commonly given by $R_{ab}$, so we denote the spacetime Ricci tensor by ${}^{(4)}R_{\alpha\beta}$.

This allows us to define the Weyl scalars as follows,

[TABLE]

for a complex null tetrad given by,

[TABLE]

where $r^{\mu}$ is an outward pointing space-like unit vector. The orientation of the tetrad is chosen so that in Minkowski spacetime we have $\theta^{i}=\hat{x}^{i}$ and $\phi^{i}=\hat{y}^{i}$ on the $z$ axis and everywhere else $\theta^{i}$ and $\phi^{i}$ are defined in the usual way on the sphere, for more details about the tetrad see Section D in [166], but note that they use the opposite sign defining the Weyl scalars.

The sign convention of the strain polarizations $h_{+}$ and $h_{\times}$ are chosen as follows. We define $h_{\mu\nu}=\psi_{\mu\nu}-\eta_{\mu\nu}$ where $\eta_{\mu\nu}$ is the Minkowski metric. For a gravitational wave propagating in the outward radial direction, we define the strain as

[TABLE]

where $\hat{\theta}$ and $\hat{\phi}$ correspond to the usual coordinate vectors on the sphere. From Eqs. (52)–(54) we have the following relation between $\Psi_{4}$ and the second time derivative of the strain,

[TABLE]

The sign of the strain defined in Eqs. (52)–(54) is the current definition of the strain for the waveforms in our catalog. This differs by an overall sign for any strain waveform previously acquired from our catalog. Please see the end of Appendix A.3.1 for details on working with previous and current versions of the waveform files. In defining the Weyl scalars—and thus the relationship between the Weyl tensor and the strain—there is a significant representation in the literature both agreeing with our sign convention [243, 244, 35] and having the opposite sign convention [166, 20, 245, 246]. Newman and Penrose originally defined the Weyl scalars opposite to ours, but they also used a $(+,-,-,-)$ metric signature [247].

For the Regge-Wheeler and Zerilli scalars $\Phi^{(\pm)}$ we choose our sign convention so that Eqs. (10) and (11) give the same polarizations for a linearized transverse-traceless gauge wave as Eqs. (52) and (53). This means we have the same sign of $\Phi^{(-)}$ but the opposite sign of $\Phi^{(+)}$ as Ref. [170], and the same signs of $\Phi^{(\pm)}$ as Ref. [164].

Our waveform quantities are decomposed in terms of spin-weighted spherical harmonics, cf. Eqs. (9) and (11). We use the sign conventions for spin-weighted spherical harmonics as given in Ref. [126, 235]. In terms of $\theta$ and $\phi$ we give the spin-weight $-2$ spherical harmonics for $l=2$ as an example,

[TABLE]

The following is a brief selection of how our sign conventions for the strain $h$ and RWZ scalars $\Phi^{(\pm)}$ compare with a few other sources in the literature:

• •

The signs of $h$ and $\Phi^{(+)}$ differ by an overall sign compared to our previous catalog release [121].

• •

Our conventions agree with [126, 164], except that it appears that Eq. (II.5) of [126] has a sign error.

• •

Ref. [170] defines $\Phi^{(+)}$ in Eq. (29) with the opposite sign as we use here, but their definition of $\Phi^{(-)}$ in Eq. (18) agrees with ours.

• •

The same sign differences between this paper and [170] also appear in [248]. In [248], Equation (15) should have an overall sign change to match our convention and the second equation (unnumbered) in Sec. 3.3 should have the opposite sign on the $\Phi^{(+)}$ term.

• •

Ref. [171] defines vector and tensor spherical harmonics with the opposite sign to ours, which would indicate that their odd-parity RWZ function has the opposite sign of ours. However our sign for $\Phi^{(+)}$ is opposite their definition, so our overall definition of the strain $h\sim\Phi^{(+)}+i\Phi^{(-)}$ agrees with theirs up to a sign. This conclusion assumes that both papers use the same sign convention for the surface volume form $\epsilon_{AB}$, which is not made clear in their paper. There also appears to be a factor of 2 difference in the definitions of the rank-2 tensor spherical harmonics, but this might be due to unclear notation.