Surrogate models for precessing binary black hole simulations with unequal masses

TL;DR

This paper introduces two new surrogate models for binary black hole mergers that are faster and more accurate than previous models, covering higher mass ratios and generic spins, aiding gravitational wave data analysis.

Contribution

The paper presents NRSur7dq4 and emnantModel, extending surrogate modeling to higher mass ratios and generic spins with improved accuracy and efficiency.

Findings

Models outperform existing ones by at least an order of magnitude.

Accurately predict waveform and remnant properties within training range.

Effective even when extrapolated beyond training parameters.

Abstract

Only numerical relativity simulations can capture the full complexities of binary black hole mergers. These simulations, however, are prohibitively expensive for direct data analysis applications such as parameter estimation. We present two new fast and accurate surrogate models for the outputs of these simulations: the first model, NRSur7dq4, predicts the gravitational waveform and the second model, \RemnantModel, predicts the properties of the remnant black hole. These models extend previous 7-dimensional, non-eccentric precessing models to higher mass ratios, and have been trained against 1528 simulations with mass ratios $q\leq4$ and spin magnitudes $\chi_1,\chi_2 \leq 0.8$, with generic spin directions. The waveform model, NRSur7dq4, which begins about 20 orbits before merger, includes all $\ell \leq 4$ spin-weighted spherical harmonic modes, as well as the precession frame dynamics and spin evolution of the black holes. The final black hole model, \RemnantModel, models the mass, spin, and recoil kick velocity of the remnant black hole. In their training parameter range, both models are shown to be more accurate than existing models by at least an order of magnitude, with errors comparable to the estimated errors in the numerical relativity simulations. We also show that the surrogate models work well even when extrapolated outside their training parameter space range, up to mass ratios $q=6$.

Equations

$$\text{mathpzc}h(t,\iota ,{\phi }_0)=\sum _{\ell =2}^{\infty }\sum _{m=-l}^l\text{mathpzc}{h}_{\ell m}(t)\text{leavevmode}{}_{-2}{Y}_{\ell m}(\iota ,{\phi }_0).$$

$$\text{mathpzc}{h}_{\ell m}^{\mathrm{coorb}}(t)=\text{mathpzc}{h}_{\ell m}^{\mathrm{copr}}(t)\text{leavevmode}{e}^{im\varphi (t)}.$$

$$\varphi (t)=\frac{\text{arg}[\text{mathpzc}{h}_{2,-2}^{\mathrm{copr}}(t)]-\text{arg}[\text{mathpzc}{h}_{2,2}^{\mathrm{copr}}(t)]}{4},$$

$$d{s}^2={\left(\frac{q^a-q^b}{\Delta q}\right)}^2+\sum _{i\in \{1,2\}}{\left(\frac{|{\mathbf{\chi }}_i^a-{\mathbf{\chi }}_i^b|}{\Delta \chi }\right)}^2,$$

$$A(t)=\sqrt{\sum _{\ell m}|\text{mathpzc}{h}_{\ell m}(t){|}^2}.$$

$$\text{mathpzc}{h}_{\ell m}^{\pm }=\frac{\text{mathpzc}{h}_{\ell ,m}^{\mathrm{coorb}}\pm \text{mathpzc}{h}_{\ell ,-m}^{\mathrm{coorb}}{}^{\ast }}{2}.$$

$$\mathbf{\Lambda }=[\log (q),{\chi }_{1x}^{\mathrm{coorb}},{\chi }_{1y}^{\mathrm{coorb}},{\hat{\chi }}^{\mathrm{coorb}},{\chi }_{2x}^{\mathrm{coorb}},{\chi }_{2y}^{\mathrm{coorb}},{\chi }_a^{\mathrm{coorb}}],$$

$${\hat{\chi }}^{\mathrm{coorb}}=\frac{{\chi }_{\mathrm{eff}}^{\mathrm{coorb}}-38\eta ({\chi }_{1z}^{\mathrm{coorb}}+{\chi }_{2z}^{\mathrm{coorb}})/113}{1-76\eta /113},$$

$${\chi }_{\mathrm{eff}}^{\mathrm{coorb}}=\frac{q\text{leavevmode}{\chi }_{1z}^{\mathrm{coorb}}+{\chi }_{2z}^{\mathrm{coorb}}}{1+q},$$

$$\eta =\frac{q}{(1+q{)}^2},$$

$${\chi }_a^{\mathrm{coorb}}=\frac{1}{2}({\chi }_{1z}^{\mathrm{coorb}}-{\chi }_{2z}^{\mathrm{coorb}}).$$

$$\mathcal{MM}=1-\frac{⟨\text{mathpzc}{h}_1,\text{mathpzc}{h}_2⟩}{\sqrt{⟨\text{mathpzc}{h}_1,\text{mathpzc}{h}_1⟩⟨\text{mathpzc}{h}_2,\text{mathpzc}{h}_2⟩}},$$

$$⟨\text{mathpzc}{h}_1,\text{mathpzc}{h}_2⟩=4\mathrm{Re}{\int }_{f_{\min }}^{f_{\max }}\frac{{\widetilde{\text{mathpzc}h}}_1(f){\widetilde{\text{mathpzc}h}}_2^{\ast }(f)}{S_n(f)}df,$$

Full text

I Introduction

As the LIGO Aasi et al. (2015) and Virgo Acernese et al. (2015) detectors reach their design sensitivity, gravitational wave (GW) detections Abbott et al. (2016a, 2017a, b, 2017b, 2017c, 2017d, 2018a) are becoming routine Abbott et al. (2018b, c). To maximize the science output of the data collected by the network of detectors, it is crucial to accurately model the source of the GWs. Among the most important sources for these detectors are binary black hole (BBH) systems, in which two black holes (BHs) lose energy through GWs, causing them to inspiral and eventually merge.

Numerical relativity (NR) simulations are necessary to accurately model the late inspiral and merger stages of the BBH evolution. These simulations accurately solve Einstein’s equations to predict the evolution of the BBH spacetime. The most important outputs of NR simulations are the gravitational waveform and the mass, spin, and recoil kick velocity of the remnant BH left after the merger.

For interpreting detected signals, model waveforms are used to compare with detector data and infer the properties of the source Cutler and Flanagan (1994); Abbott et al. (2016c); Veitch et al. (2015). The mass and spin of the remnant determine the black hole ringdown frequencies, which are used in testing general relativity Abbott et al. (2016d); LIG (2019); Ghosh et al. (2018). In addition, the recoil kick is astrophysically important because it can cause the remnant BH to be ejected from its host galaxy Campanelli et al. (2007a); Gonzalez et al. (2007a); Gerosa and Sesana (2015).

Unfortunately, NR simulations are too expensive to be directly used in data analysis applications and incorporated into astrophysical models. As a result, several approximate models that are much faster to evaluate have been developed for both waveforms Khan et al. (2018); Cotesta et al. (2018); London et al. (2018); Pan et al. (2014); Bohé et al. (2017); Khan et al. (2016); Hannam et al. (2014); Taracchini et al. (2014); Pan et al. (2011); Mehta et al. (2017); Babak et al. (2017) and remnant properties Hofmann et al. (2016); Barausse et al. (2012); Jiménez-Forteza et al. (2017); Healy and Lousto (2017); Healy et al. (2014); Gonzalez et al. (2007b); Campanelli et al. (2007b); Lousto and Zlochower (2008); Lousto et al. (2012); Lousto and Zlochower (2013); Gerosa and Kesden (2016); Healy and Lousto (2018); Herrmann et al. (2007); Campanelli et al. (2007a); Gonzalez et al. (2007a); Rezzolla et al. (2008a, b); Kesden (2008); Tichy and Marronetti (2008); Barausse and Rezzolla (2009); Zlochower and Lousto (2015). These models typically assume an underlying phenomenology based on physical motivations, and calibrate any remaining free parameters to NR simulations.

Among BBHs, systems with BH spins that are misaligned with respect to the orbital angular momentum are complicated to model analytically or semi-analytically. For these systems, the spins interact with both the orbital angular momentum and each other, causing the system to precess about the direction of the total angular momentum Apostolatos et al. (1994). This precession is imprinted on the waveform as characteristic modulations in the amplitude and frequency of the GWs, and can be used to extract information about the spins of the source. One important application of the extracted spins is to distinguish between formation channels of BBHs Gerosa et al. (2013); Vitale et al. (2017); Farr et al. (2018); Gerosa et al. (2018a).

The precessing BBH problem for quasicircular orbits is parametrized by seven parameters: the mass ratio $q=m_{1}/m_{2}\geq 1$ and two spin vectors $\bm{\chi}_{1,2}$, where the index 1 (2) refers to the heavier (lighter) BH. The total mass scales out of the problem and does not constitute an additional parameter for modeling. The surrogate models of Ref. Blackman et al. (2017a) for the gravitational waveform, and Ref. Varma et al. (2019a) for the remnant properties, were the first to model the $7-$dimensional space of generically precessing BBH systems, albeit restricted to mass ratios $q\leq 2$, and dimensionless spin magnitudes $\chi_{1,2}\leq 0.8$. Trained directly against numerical simulations, these models do not need to introduce additional assumptions about the underlying phenomenology of the waveform or remnant properties that necessarily introduces some systematic error. Through cross-validation studies, it was shown that both these models achieve accuracies comparable to the numerical simulations themselves Blackman et al. (2017a); Varma et al. (2019a), and as a result, are the most accurate models currently available for precessing systems, within their parameter space of validity.

In this paper, we present extensions of the above surrogate models to larger mass ratios. Our new surrogate models are called NRSur7dq4 and NRSur7dq4Remnant, for the gravitational waveform and remnant properties, respectively. They are trained against 1528 precessing NR simulations with mass ratios $q\leq 4$, spin magnitudes $\chi_{1},\chi_{2}\leq 0.8$, and generic spin directions. Both models are made publicly available through the gwsurrogate Blackman et al. and surfinBH Varma et al. (a) Python packages; example evaluation codes are provided at Ref. SpE (2018) and Ref. Varma et al. (a), respectively, for NRSur7dq4 and NRSur7dq4Remnant.

The rest of the paper is organized as follows. Section. II covers some preliminaries to set up the modeling problem for precessing BBH systems. Section III describes the training simulations. Sec. IV describes the NRSur7dq4 waveform surrogate model. Section V describes the NRSur7dq4Remnant remnant properties surrogate model. Section VI compares these models against NR simulations to assess their accuracy. Finally, Sec. VII presents some concluding remarks. In App. A we examine how accurate these models are when extrapolated beyond mass ratio $q=4$, and in App. B we investigate some features in the error distribution of the NR simulations.

II Preliminaries and notation

It is convenient to combine the two polarizations of the waveform into a single complex, dimensionless strain $\mathpzc{h}=h_{+}-ih_{\times}$, and to represent the waveform on a sphere as a sum of spin-weighted spherical harmonic modes:

[TABLE]

Here ${}_{-2}Y_{\ell m}$ are the spin$\,=\!\!-2$ weighted spherical harmonics, and $\iota$ and $\varphi_{0}$ are the polar and azimuthal angles on the sky in the source frame.

For nonprecessing systems, the direction of orbital angular momentum ($\bm{L}$) is fixed and the $\hat{z}$ direction of the source frame is chosen to be along $\hat{\bm{L}}$ by convention. The gravitational radiation is strongest along the directions parallel and antiparallel to $\hat{\bm{L}}$. Therefore, for nonprecessing systems the quadrupole modes ($\ell=2,m=\pm 2$) dominate the sum in Eq. (1), but the nonquadrupole modes can become important at large mass ratios or $\iota$ close to $\pi/2$ Varma and Ajith (2017); Capano et al. (2014); Littenberg et al. (2013); Calderón Bustillo et al. (2017); Brown et al. (2013); Varma et al. (2014); Graff et al. (2015); Harry et al. (2018); Calderón Bustillo et al. (2016); Pekowsky et al. (2013).

By contrast, for precessing systems the direction of $\bm{L}$ varies due to precession Apostolatos et al. (1994) and so there is not a fixed axis along which the radiation is dominant. The standard practice is to choose $\hat{z}$ of the source frame along the direction of $\bm{L}$ (or the total angular momentum) at a reference time or frequency.

Heuristically, one can think of a precessing system as a nonprecessing system with time-dependent frame rotations applied to it. In this non-inertial frame the rotation causes mixing of power between modes of fixed $\ell$. For example, the power of the $(2,\pm 2)$ modes leaks into the $(2,\pm 1)$ and $(2,0)$ modes. This means that all $\ell=2$ modes can be dominant in Eq. (1). While this rotating-frame picture ignores some dynamical features such as nutation, it accounts for most of the effects of precession in the waveform.

By the same logic, one could apply a time-dependent rotation to a precessing system such that $\hat{z}$ always lies along $\hat{\bm{L}}(t)$. In this non-inertial frame, referred to as the coprecessing frame Schmidt et al. (2011); O’Shaughnessy et al. (2011); Boyle et al. (2011), the radiation is always strongest along $\hat{z}$, and the ($\ell=2,m=\pm 2$) modes are dominant. In fact, since most precessional effects are accounted for by the frame rotation, the waveform in the coprecessing frame is qualitatively similar to that of a nonprecessing system (cf. Fig. 1). This observation has been exploited in the literature Hannam et al. (2014); Pan et al. (2014); Khan et al. (2018); Blackman et al. (2017b, a) to simplify the modeling of precessing systems. Here we proceed similarly, using the coprecessing frame described in Ref. Boyle et al. (2011) and denoting the strain in this frame as $\mathpzc{h}^{\mathrm{copr}}_{\ell m}$.

The waveform can be made even simpler, and therefore easier to model, by applying an additional rotation about the $z-$axis of the coprecessing frame by an amount equal to the instantaneous orbital phase:

[TABLE]

Here we define the orbital phase,

[TABLE]

using the coprecessing frame strain. The waveform $\mathpzc{h}^{\mathrm{coorb}}_{\ell m}(t)$ corresponds to a new frame, called the coorbital frame, in which the BHs are always on the $x-$axis, with the heavier BH on the positive $x-$axis111Here the BH positions are defined from the waveform at future null infinity and do not necessarily correspond to the (gauge-dependent) coordinate BH positions in the NR simulation.. More importantly, the waveform in the coorbital frame is nearly nonoscillatory, simplifying the modeling problem greatly. Figure 1 shows an example of a waveform in the inertial, coprecessing, and coorbital frames.

III NR simulations

Our NR simulations are performed using the Spectral Einstein Code (SpEC) SpE ; Pfeiffer et al. (2003); Lovelace et al. (2008); Lindblom et al. (2006); Szilagyi et al. (2009); Scheel et al. (2009) developed by the SXS SXS collaboration.

III.1 Parameter space coverage

We use 890 precessing NR simulations used in the construction of the surrogate models of Refs. Blackman et al. (2017a); Varma et al. (2019a), which provide coverage in the $q\leq 2$ and $\chi_{1},\chi_{2}\leq 0.8$ regions of the parameter space. We also make use of 64 aligned-spin simulations with $q\leq 4$ and $\chi_{1},\chi_{2}\leq 0.8$ used in the construction of the surrogate model presented in Ref. Varma et al. (2019b). Finally, we performed 574 new simulations with $2<q\leq 4$, $\chi_{1},\chi_{2}\leq 0.8$ and generic spin directions—these simulations are presented here for the first time. The parameters for the first 204 of these are chosen based on sparse grids as detailed in Appendix A of Ref. Blackman et al. (2017a). The remaining parameters are chosen as follows. We randomly sample 1000 points uniformly in mass ratio, spin magnitude, and spin direction on the sphere. We compute the distance between points a and b using the metric

[TABLE]

where $\Delta q=4-1=3$ and $\Delta\chi=0.8$ are the ranges of these parameters. These normalization factors are somewhat arbitrary, although any choice of order unity should provide a reasonable criteria for point selection. For each sampled parameter, we compute the minimum distance to all previously chosen parameters. We then add the sampled parameter maximizing this minimum distance to the set of chosen parameters. This is done iteratively for 370 additional parameters. The new simulations have identifiers SXS:BBH:1346-1350 and SXS:BBH:1514-2082, and are made publicly available through the SXS public catalog SXS Collaboration . The parameter space covered by the 890+64+574=1528 NR simulations used in this work is shown in Fig. 2. Note that not all of these are independent simulations: for 154 of these cases we have $q=1$, with $\bm{\chi}_{1}\neq\bm{\chi}_{2}$; for each of these cases we effectively obtain an additional simulation by exchanging the labels of the two BHs.

The start time of these simulations varies between $4693M$ and $5234M$ before the peak of the waveform amplitude, where $M=m_{1}+m_{2}$ is the total Christodoulou mass measured close to the beginning of the simulation at the “relaxation time” Boyle et al. (2019). The initial orbital parameters are chosen through an iterative procedure Buonanno et al. (2011) such that the orbits are quasicircular; the largest eccentricity for these simulations is $9.8\times 10^{-4}$, while the median value is $3.8\times 10^{-4}$.

III.2 Data extracted from simulations

We make use of the following quantities extracted from the NR simulations: the waveform modes $\mathpzc{h}_{\ell m}(t)$, the component spins $\bm{\chi}(t)$, the mass ratio $q$, and the remnant mass $m_{f}$, spin $\bm{\chi}_{f}$, and kick velocity $\bm{v}_{f}$.

The waveform is extracted at several extraction spheres at varying finite radii from the origin and then extrapolated to future null infinity Boyle et al. (2019); Boyle and Mroue (2009). Then the extrapolated waveforms are corrected to account for the initial drift of the center of mass Boyle (2016, ). The time steps during the simulations are chosen nonuniformly using an adaptive time-stepper Boyle et al. (2019). Using cubic splines, we interpolate the real and imaginary parts of the waveform modes to a uniform time step of $0.1M$; this is dense enough to capture all frequencies of interest, including near merger. The interpolated waveform at future null infinity, scaled to unit mass and unit distance, is denoted as $\mathpzc{h}_{\ell m}(t)$ in this paper.

The component spins $\bm{\chi}_{1,2}(t)$ and masses $m_{1,2}$ are evaluated on the apparent horizons Lovelace et al. (2008) of the BHs. The masses at the relaxation time Boyle et al. (2019) are used to define the mass ratio $q=m_{1}/m_{2}$. Unless otherwise specified, all masses in this paper are given in units of the total mass $M=m_{1}+m_{2}$ at relaxation. The spins are interpolated onto the same time array222The waveforms at future null infinity use a time coordinate $t$ that is different from the simulation time $\tilde{t}$ at which the spins are measured in the near zone Boyle et al. (2019). In this paper, we identify $t$ with $\tilde{t}$. While this identification is gauge-dependent, the spin directions are already gauge-dependent. We, however, note that the spin and orbital angular momentum vectors in the damped harmonic gauge used by SpEC agree quite well with the corresponding vectors in post-Newtonian (PN) theory Ossokine et al. (2015a). as used for the waveform, using cubic splines.

The remnant mass $m_{f}$ and spin $\bm{\chi}_{f}$ are determined from the common apparent horizon long after ringdown, as detailed in Ref. Boyle et al. (2019). The remnant kick velocity is derived from conservation of momentum, $\bm{v}_{f}=-\bm{P}^{\mathrm{rad}}/m_{f}$ Gerosa et al. (2018b). The radiated momentum flux $\bm{P}^{\mathrm{rad}}$ is integrated Ruiz et al. (2008) from the strain $\mathpzc{h}_{\ell m}$.

III.3 Post-processing the output of NR simulations

After extracting the strain and spins from the simulations, we apply the following post processing steps before building the surrogate models.

First, we shift the time arrays of all waveforms such that $t=0$ occurs at the peak (see Ref. Blackman et al. (2017a) for how the peak is determined) of the total waveform amplitude, defined as:

[TABLE]

Then we rotate the waveform modes such that at a reference time $t_{0}=-4300M$, the inertial frame coincides with the coorbital frame. This means that the $\hat{z}$ direction of the inertial frame is along the principal eigenvector of the angular momentum operator Boyle et al. (2011) at the reference time. In addition, the $\hat{x}$ direction of the inertial frame is along the line of separation from the lighter BH to the heavier BH (in other words, the orbital phase is zero). The spin vectors $\bm{\chi}_{1,2}(t)$ are also transformed into the same inertial frame.

We then truncate the waveform and spin time series by dropping all times $t<-4300M$ to exclude the initial transients known as “junk radiation”. After the truncation, the reference time $t=-4300M$ is also the start time of the data.

For $t>-100M$, the spin measurements from the apparent horizons start to become unreliable as the horizons become highly distorted. Following Ref. Blackman et al. (2017a), starting at $t\!=\!-100M$, we extend the spins to later times using PN spin evolution equations. This evolution is done even past the merger stage, into the ringdown. We stress that the extended spins are unphysical but are a useful parametrization to construct fits at late times.

Finally we apply a smoothing filter (see Eq. (6) of Ref. Blackman et al. (2017a)) on the spin time series to remove fast oscillations taking place on the orbital timescale. This smoothing helps improve the numerical stability of the ordinary differential equation (ODE) integrations described in Sec. IV.2. Note that we use the filtered spins for the waveform surrogate (Sec. IV) but not for the remnant surrogate (Sec. V), for which we just use the unfiltered spins since there are no ODE integrations involved.

IV Waveform surrogate

To construct the waveform surrogate, we closely follow the model of Ref. Blackman et al. (2017a), with some modifications to adapt it to higher mass ratios. We refer to the new waveform model as NRSur7dq4.

IV.1 Coorbital frame surrogate

We find that the surrogate accuracy improves when working with slowly varying functions, rather than oscillatory ones. Therefore, we first decompose the strain into several “data pieces”, each of which is a slowly varying function of time, and build a surrogate for each of them. At evaluation time, we combine the various data pieces to reconstruct the inertial frame strain. To reduce the cost of these transformations, we first downsample the inertial frame strain onto a set of 2000 time values $t^{\mathrm{coorb}}_{i}$ that are approximately uniformly spaced in the orbital phase (using the method described in App. B of Ref. Blackman et al. (2017a)).

As described in Sec. II, the waveform is simpler in the coorbital frame. A further simplification is possible by considering combinations of $m>0$ and $m<0$ counterparts of a fixed $\ell$ mode:

[TABLE]

Figure 3 shows an example of the simplification obtained with this combination. For all $m\neq 0$ modes we model the real and imaginary parts of $\mathpzc{h}_{\ell m}^{\pm}$. For $m=0$ modes, we directly model the real and imaginary parts of the coorbital frame strain $\mathpzc{h}^{\mathrm{coorb}}_{\ell,m}$. We construct an independent surrogate model for each of these data pieces and refer to the combination of these models as the coorbital frame surrogate.

As described in Ref. Blackman et al. (2017a), for each waveform data piece, we construct a linear basis using singular value decomposition with an RMS tolerance of $3\times 10^{-4}$. We then construct an empirical time interpolant with the same number of empirical time nodes as basis functions for that data piece Field et al. (2014); Maday et al. (2009); Chaturantabut and Sorensen (2010). The empirical time nodes are chosen as a subset of the 2000 coorbital time values ($t^{\mathrm{coorb}}_{i}$) and are specific to each data piece. Finally, for each empirical time node, we construct a parametric fit for the waveform data piece. The fits are parametrized as functions of the mass ratio and the spins in the coorbital frame at that time. We describe our fitting procedure in Sec. IV.3. At evaluation time, the coorbital frame spins at any time are obtained using the dynamics surrogate described in Sec. IV.2.

IV.2 Dynamics surrogate

The surrogate described in Sec. IV.1 only models the strain in the coorbital frame. We also need to model the following quantities:

The orbital phase in the coprecessing frame, which is required to transform the strain from the coorbital frame to the coprecessing frame [cf. Eq. (2)]; 2. 2.

The quaternions describing the coprecessing frame, which are required to transform the strain from the coprecessing frame to the inertial frame; 3. 3.

The spins as a function of time, which are used in the evaluation of the parametric fits described in Sec. IV.3.

We refer to the model for these quantities as the dynamics surrogate. Using the fitting method of Sec. IV.3, we first construct parametric fits for $\omega(t)$, $\Omega^{\mathrm{coorb}}_{x,y}(t)$, and $\dot{\bm{\chi}}^{\mathrm{coorb}}_{1,2}(t)$ at selected time nodes referred to as the dynamical time nodes $t^{\mathrm{dyn}}_{i}$. Here $\dot{\bm{\chi}}^{\mathrm{coorb}}_{1,2}(t)$ are the time derivatives of the coprecessing frame spins transformed to the coorbital frame, $\omega(t)$ is ${d\phi}/{dt}$ (cf. Eq. (3)), and $\Omega^{\mathrm{coorb}}_{x,y}(t)$ is the angular velocity of the coprecessing frame, transformed to the coorbital frame. These quantities are described in more detail in Sec. III of Ref. Blackman et al. (2017a). Note that $\Omega^{\mathrm{coorb}}_{z}(t)\sim 0$. For the dynamical time nodes $t^{\mathrm{dyn}}_{i}$ we chose 238 time values such that there are approximately $10$ nodes per orbit (see App. B of Ref. Blackman et al. (2017a) for details).

We use a fourth-order Adams-Bashforth scheme to integrate $\omega(t^{\mathrm{dyn}}_{i})$, $\Omega^{\mathrm{coorb}}_{x,y}(t^{\mathrm{dyn}}_{i})$, and $\dot{\bm{\chi}}^{\mathrm{coorb}}_{1,2}(t^{\mathrm{dyn}}_{i})$ over the set of dynamical time nodes $t^{\mathrm{dyn}}_{i}$ providing the time evolution of the orbital phase $\phi(t^{\mathrm{dyn}}_{i})$, the coprecessing frame quaternions $\hat{Q}(t^{\mathrm{dyn}}_{i})$, and the component spins in the coorbital frame $\bm{\chi}^{\mathrm{coorb}}_{1,2}(t^{\mathrm{dyn}}_{i})$. This involves solving a coupled ODE as described in Sec. V of Ref. Blackman et al. (2017a). At each step of the ODE integration, the coorbital frame spins at the current node $t^{\mathrm{dyn}}_{i}$ are first obtained. These are then used to evaluate the parametric fits for the derivative quantities mentioned above. Note that the spins used in the dynamics surrogate are the filtered spins mentioned in Sec. III.3; this improves the accuracy of the ODE integration by making the spin time derivatives easier to model.

IV.3 Parametric fits

For the coorbital frame surrogate of Sec. IV.1, we need to construct parametric fits at various empirical time nodes for the different data pieces. Similarly, for the dynamics surrogate of Sec. IV.2, we need to construct fits for various time derivatives at the dynamical time nodes $t^{\mathrm{dyn}}_{i}$. We use the same procedure for each of these fits. Let us refer to the data to be fitted as $y(\bm{\Lambda})$, where $\bm{\Lambda}$ is a seven-dimensional set of parameters.

For each of these fits, the seven parameters $\bm{\Lambda}$ must contain information on mass ratio $q$ and coorbital frame spins $\bm{\chi}^{\mathrm{coorb}}_{1,2}(t_{i})$ at the time corresponding to the fit. Following Ref. Varma et al. (2019a), we parametrize the fits using

[TABLE]

where $\hat{\chi}^{\mathrm{coorb}}$ is the spin parameter entering the GW phase at leading order Khan et al. (2016); Ajith (2011); Cutler and Flanagan (1994); Poisson and Will (1995) in the PN expansion

[TABLE]

and $\chi^{\mathrm{coorb}}_{a}$ is the “anti-symmetric spin”,

[TABLE]

We empirically found this parameterization to perform more accurately than the more intuitive choice $\bm{\Lambda}_{\mathrm{ref56}}=[q,\chi^{\mathrm{coorb}}_{1x},\chi^{\mathrm{coorb}}_{1y},\chi^{\mathrm{coorb}}_{1z},\chi^{\mathrm{coorb}}_{2x},\chi^{\mathrm{coorb}}_{2y},\chi^{\mathrm{coorb}}_{2z}]$ used in Ref. Blackman et al. (2017a).

Fits are constructed using the forward-stepwise greedy fitting method described in App. A of Ref. Blackman et al. (2017b). We choose the basis functions to be a tensor product of 1D monomials in the components of $\bm{\Lambda}$. The components of $\bm{\Lambda}$ are first affine mapped to the interval $[-1,1]$ before constructing the tensor product. We consider up to cubic powers in $\log(q)$ and up to quadratic powers in the spin parameters. We find that going to higher powers does not significantly improve the fit accuracy within the training region, but the mass ratio extrapolation errors estimated in App. A become much larger.

It is always possible to improve the accuracy of a fit by adding more basis functions. However, this can lead to over-fitting when the data contain some noise. Our source of noise is mostly due to NR truncation error, but also systematic errors such as waveform extrapolation and residual eccentricity. In order to safeguard against over-fitting, we perform 10 trial fits, leaving a random $10\%$ of the dataset out as validation points in each trial, to determine the set of basis functions used in constructing the final fit. We allow a maximum of 100 basis functions for each fit. See App. A of Ref. Blackman et al. (2017b) for more details.

IV.4 Surrogate evaluation

To evaluate the surrogate, we begin with a user-specified mass ratio $q$ and spins $\bm{\chi}^{\mathrm{coorb}}_{1,2}$ at the initial time $t=-4300M$. Note that at this time, the inertial frame coincides with the coorbital frame. These values are used to initialize the dynamics surrogate described in Sec. IV.2, which predicts the coprecessing frame quaternions $\hat{Q}(t^{\mathrm{dyn}}_{i})$, the orbital phase $\phi(t^{\mathrm{dyn}}_{i})$ in the coprecessing frame, and the coorbital frame spins $\bm{\chi}^{\mathrm{coorb}}_{1,2}(t^{\mathrm{dyn}}_{i})$ at the dynamic time nodes $t^{\mathrm{dyn}}_{i}$. We then use cubic splines to interpolate these quantities on to the time array for the coorbital frame surrogate $t^{\mathrm{coorb}}_{i}$, giving us $\hat{Q}(t^{\mathrm{coorb}}_{i})$, $\phi(t^{\mathrm{coorb}}_{i})$, and $\bm{\chi}^{\mathrm{coorb}}_{1,2}(t^{\mathrm{coorb}}_{i})$.

The coorbital frame surrogate described in Sec. IV.1 is used to predict the strain in the coorbital frame. This involves evaluating the fits at the empirical time nodes for this surrogate using $\bm{\chi}^{\mathrm{coorb}}_{1,2}(t^{\mathrm{coorb}}_{i})$ and $q$. Then, the orbital phase $\phi(t^{\mathrm{coorb}}_{i})$ is used to transform the strain from the coorbital frame to the coprecessing frame (cf. Eq. 2). Finally, the coprecessing frame quaternions $\hat{Q}(t^{\mathrm{coorb}}_{i})$ are used to transform the strain from the coprecessing frame to the inertial frame (this involves Wigner matrices, see App. A of Ref. Boyle et al. (2011)). This gives us $\mathpzc{h}_{\ell m}(t^{\mathrm{coorb}}_{i})$, which is interpolated onto any required time array $t$ using cubic splines to get $\mathpzc{h}_{\ell m}(t)$.

V Remnant surrogate

To construct the remnant properties surrogate, we closely follow the model of Ref. Varma et al. (2019a). We refer to the new model presented here as NRSur7dq4Remnant.

We model the remnant mass $m_{f}$, spin $\bm{\chi}_{f}$, and kick velocity $\bm{v}_{f}$. Before constructing the fits, $\bm{\chi}_{f}$ and $\bm{v}_{f}$ are transformed into the coorbital frame at $t\!=\!-100M$. We model each component of the vectors independently. The fits are parametrized by the same $\bm{\Lambda}$ of Eq. (7), but using the component spins at $t\!=\!-100M$. Unlike the waveform surrogate case, we do not filter out orbital-timescale oscillations. The filtered spins were found to be necessary for the accuracy of the time integration in Sec. IV.2, which is not necessary here because the remnant properties can evaluated from the BBH parameters at a single time $t\!=\!-100M$.

All fits are performed using Gaussian Process Regression (GPR), as described in the supplementary materials of Ref. Varma et al. (2019a). We find that GPR fitting is, in most cases, more accurate but also significantly more expensive than the polynomial fitting method described in Sec. IV.3. GPR becomes impractical to use for the waveform surrogate as there are hundreds of fits that need to be evaluated to generate the waveform. For the remnant fits, however, the additional cost of GPR is acceptable because one is only fitting 7 quantities ($m,\bm{\chi}_{f},\bm{v}_{f}$). In addition, GPR naturally provides error estimates which can be useful in data analysis applications. The efficacy of the GPR error estimate in reproducing the underlying error of the surrogate models was investigated thoroughly in the supplementary materials of Ref. Varma et al. (2019a).

Although NRSur7dq4Remnant is parameterized internally by input spins specified in the coorbital frame at $t\!=\!-100M$, we allow the user to specify input spins at earlier times, and in the inertial frame; this case is handled by two additional levels of spin evolution. Given the inertial-frame input spins at an initial orbital frequency $f_{0}$, we first evolve the spins using a post-Newtonian (PN) approximant — 3.5PN SpinTaylorT4 Buonanno et al. (2003); Boyle et al. (2007); Ossokine et al. (2015a) — until we reach the domain of validity of the more accurate NRSur7dq4 ($t=-4300M$ from the peak). We then use the dynamics surrogate of NRSur7dq4 to evolve the spins until $t\!=\!-100M$. These spins are then transformed to the coorbital frame and used to evaluate the remnant fits. Thus, spins can be specified at any given orbital frequency and are evolved consistently before estimating the final BH properties. Note that NRSur7dq4 uses the filtered spins, while NRSur7dq4Remnant expects unfiltered spins at $t\!=\!-100M$, but we find that the errors introduced by this discrepancy are negligible compared to the errors due to PN spin evolution.

VI Results

We evaluate the accuracy of our new surrogate models by comparing against the waveform and remnant properties from the NR simulations used in this work. For this, we perform a 20-fold cross-validation study to compute “out-of-sample” errors as follows. We first randomly divide the 1528 training simulations into 20 groups of $\mathchar 21016\relax\,76$ simulations each. For each group, we build a trial surrogate using the $\mathchar 21016\relax\,1452$ remaining training simulations and test against these $\mathchar 21016\relax\,76$ validation ones, which may include points on the boundary of the training set.

VI.1 Waveform surrogate errors

To estimate the difference between two waveforms, $\mathpzc{h}_{1}$ and $\mathpzc{h}_{2}$, we use the mismatch

[TABLE]

where $\tilde{\mathpzc{h}}(f)$ indicates the Fourier transform of the complex strain $\mathpzc{h}(t)$, ∗ indicates a complex conjugation, $\mathrm{Re}$ indicates the real part, and $S_{n}(f)$ is the one-sided power spectral density of a GW detector. We taper the time domain waveform using a Planck window McKechan et al. (2010), and then zero-pad to the nearest power of two. We further zero-pad the waveform to increase the length by a factor of eight before performing the Fourier transform. The tapering at the start of the waveform is done over $1.5$ cycles of the $(2,2)$ mode. The tapering at the end is done over the last $30M$. Note that our model contains times up to $100M$ after the peak of the waveform amplitude, and the signal has essentially died down by the last $30M$. We take $f_{\mathrm{min}}$ to be twice the waveform angular velocity (as defined by Ref. Boyle (2013)) at the end of the initial tapering window, and $f_{\mathrm{max}}$ is chosen to be 4 times the waveform angular velocity at $t=0$; the extra factor of 4 is chosen to resolve up to $m=4$ spherical-harmonic modes, with an extra margin of a factor of 2. We compute mismatches with a flat noise curve ($S_{n}=1$) as well as with the Advanced-LIGO design sensitivity noise curve LIGO Scientific Collaboration (2018). Mismatches are computed following the procedure described in Appendix D of Ref. Blackman et al. (2017b). In particular, we optimize over shifts in time, polarization angle, and initial orbital phase. Both plus and cross polarizations are treated on an equal footing by using a two-detector setup where one detector sees only the plus and the other only the cross polarization. We compute the mismatches at 37 points uniformly distributed on the sky in the source frame, and we use all available modes of a given waveform model.

Figure 5 summarizes the out-of-sample mismatches for NRSur7dq4 against the NR waveforms. In Fig. 4(a) we show mismatches computed using a flat noise curve. We compare this with the truncation error in the NR waveforms themselves, estimated by computing the mismatch between the two highest available resolutions of each NR simulation. The errors in the surrogate model are well within the estimated truncation errors of the NR simulations. In addition, we also show the errors for the waveform model SEOBNRv3 Pan et al. (2014); Babak et al. (2017), which also includes spin precession effects 333Note that SEOBNRv3 spins are specified at a reference frequency, rather than a time before merger. We choose the reference frequency such that the waveform begins at $t=-4300M$ before the waveform amplitude peak (as defined in Eq. 5).. The surrogate errors are at least an order of magnitude lower than those of SEOBNRv3.

Apart from SEOBNRv3, another model commonly used in data analysis applications is IMRPhemomPv2 Hannam et al. (2014). IMRPhemomPv2 was shown to be comparable in accuracy to SEOBNRv3 in Ref. Blackman et al. (2017a), at least in order of magnitude. Therefore, for simplicity, we do not show comparisons of IMRPhemomPv2 to NR here. Note that updated versions of both SEOBNRv3 (based on Ref. Cotesta et al. (2018)) and IMRPhemomPv2 (see Ref. Khan et al. (2018)) are under development, but are not currently available publicly. We note that these models are calibrated only against aligned-spin NR simulations, using a much smaller set of simulations than our model. Both these factors contribute to the accuracy of these models. On the other hand, these models are expected to be valid for larger mass ratios and spin magnitudes than our model, although their accuracy in that region is unknown due to lack of sufficient number of simulations.

We note that the NR truncation mismatch distribution in Fig. 4(a) has a tail extending to $\mathcal{MM}\sim 0.1$. We find that these cases occur when the spins of the two highest resolutions of the simulation are inconsistent with each other because of unresolved effects during junk-radiation emission, meaning that the two resolutions represent different physical systems. This means that comparing the resolutions for these cases gives us an error estimate that is too conservative and does not reflect the actual truncation error of the simulations. We expect the actual truncation error to be closer to the errors reproduced by the surrogate model (which is trained on the high resolution data set) in Fig. 4(a). Evidence for these claims is provided in App. B.

Fig. 4(b) shows mismatches computed using the Advanced LIGO design sensitivity noise curve LIGO Scientific Collaboration (2018). In this case, results depend on the total mass $M$ of the system. Consequently, we show the median and 95th percentile values at different $M$, rather than full histograms. Once again, the surrogate errors are comparable to those of the NR simulations, and are at least an order of magnitude lower than that of SEOBNRv3. Over the mass range $50-200M_{\odot}$, mismatches for NRSur7dq4 are always $\lesssim 8\times 10^{-3}$ at the $95$ percentile level.

Fig. 4(b) shows a comparison of waveforms computed via NRSur7dq4, SEOBNRv3, and NR for the cases that lead to the largest error for NRSur7dq4 and SEOBNRv3 in Fig. 4(a). The surrogate shows reasonable agreement with NR, even for its worst case, while SEOBNRv3 shows a noticeably larger deviation in both cases.

In Figs. 5 and 4(b) we use all available modes for NRSur7dq4 and SEOBNRv3. NRSur7dq4 models all modes $\ell\leq 4$, while SEOBNRv3 models only the $\ell=2$ modes. For the NR waveforms in Figs. 5 and 4(b), we include all modes $\ell\leq 5$ to account for the error due to neglecting $\ell>4$ modes in NRSur7dq4. To better understand what fraction of the SEOBNRv3 error comes from neglecting modes with $\ell>2$, we repeat the calculations leading to the SEOBNRv3 histogram in Fig. 4(a) in Fig. 6, while restricting all waveforms to $\ell=2$. While there is a noticeable move towards lower mismatches when restricted to $\ell=2$, the median and $95$th percentile values change only marginally, suggesting that the main error source for SEOBNRv3 are the $\ell=2$ modes themselves.

VI.2 Remnant surrogate errors

We evaluate the accuracy of the remnant surrogate NRSur7dq4Remnant by comparing against the NR simulations through a cross-validation study as in Sec. VI.1. Out-of-sample errors for the remnant properties predicted by NRSur7dq4Remnant are shown in Fig. 7. $95$th percentile errors are $\mathchar 21016\relax\,5\!\times\!10^{-4}M$ for mass, $\mathchar 21016\relax\,2\!\times\!10^{-3}$ for spin magnitude, $\mathchar 21016\relax\,4\!\times\!10^{-3}$ radians for spin direction, $\mathchar 21016\relax\,4\!\times\!10^{-4}\,c$ for kick magnitude, and $\mathchar 21016\relax\,0.2$ radians for kick direction. Our errors are at the same level as the NR resolution error, estimated by comparing the two highest NR resolutions. The largest errors in the kick direction can be of order $\mathchar 21016\relax\,1$ radian. The bottom-right panel of Fig. 7 shows the joint distribution of kick magnitude and kick direction error for NRSur7dq4Remnant, showing that direction errors are larger at low kick magnitudes. Our error in kick direction is below $\mathchar 21016\relax\,0.2$ radians whenever $v_{f}\gtrsim 2\times 10^{-3}c$.

We also compare the performance of our fits against several existing fitting formulae for remnant mass, spin, and kick which we denote as follows: HBMR (Hofmann et al. (2016); Barausse et al. (2012) with $n_{M}\!=\!n_{J}\!=\!3$), UIB Jiménez-Forteza et al. (2017), HL Healy and Lousto (2017), HLZ Healy et al. (2014), and CLZM (Gonzalez et al. (2007b); Campanelli et al. (2007b); Lousto and Zlochower (2008); Lousto et al. (2012); Lousto and Zlochower (2013) as summarized in Gerosa and Kesden (2016)). To partially account for spin precession, these fits are corrected as described in Ref. Johnson-McDaniel et al. (2016) and used in current LIGO/Virgo analyses Abbott et al. (2016e, 2017b): spins are evolved using PN from relaxation to the Schwarzschild innermost stable circular orbit, and final UIB and HL spins are post-processed by adding the sum of the in-plane spins in quadrature. Figure 7 shows that our procedure to predict remnant mass, spin magnitude, and kick magnitude for precessing systems is more accurate than these existing fits by at least an order of magnitude.

Our fits appear to outperform the NR simulations when estimating the spin direction. Once again, this is due to the post-junk-radiation initial spins of the two highest resolutions being inconsistent with each other for some of our simulations, so that different resolutions represent different physical systems (cf. App. B). Therefore, the errors estimated by comparing the two highest resolutions is a poor estimate of the actual truncation error for these cases. The actual truncation error is likely to be close to the errors reproduced by the surrogate.

The NRSur7dq4Remnant fits in Fig. 7 are evaluated using the NR spins at $t\!=\!-100M$ as inputs. In typical applications, one may have access to the spins only at the start of the waveform, rather than at $t\!=\!-100M$. For this case, as described in Sec. V, we use a combination of PN and NRSur7dq4 to evolve the spins from any given starting frequency to $t\!=\!-100M$. These spins are then used to evaluate the NRSur7dq4Remnant fits. Thus, spins can be specified at any given orbital frequency and are evolved consistently before estimating the final BH properties. This is a crucial improvement (introduced by Ref. Varma et al. (2019a)) over previous results, which, being calibrated solely to non-precessing systems, suffer from ambiguities regarding the time/frequency at which spins are defined.

Figure 8 shows the errors in NRSur7dq4Remnant when the spins are specified at an orbital frequency $f_{0}\!=\!10\,{\rm Hz}$. These errors are computed by comparing against 23 long NR ($3\times 10^{4}M$ to $10^{5}M$ in length) simulations Boyle et al. (2019) with mass ratios $q\leq 4$ and generically oriented spins with magnitudes $\chi_{1},\chi_{2}\sim 0.5$. None of these simulations were used to train the fits. Longer PN evolutions are needed at lower total masses, and the errors are therefore larger. These errors will decrease with an improved spin evolution procedure. Note, however, that our predictions are still more accurate than those of existing fitting formulae (cf. Fig. 7).

VII Conclusion

We present new NR surrogate models for precessing BBH systems with generic spins and unequal masses. In particular, we model the two most-used outputs of NR simulations: the gravitational waveform and the properties (mass, spin, and recoil kick) of the final BH formed after the merger. Trained against 1528 NR simulations with mass ratios $q\leq 4$, spin magnitudes $\chi_{1,2}\leq 0.8$, and generic spin directions, both these models are shown to reproduce the NR simulations with accuracies comparable to those of the simulations themselves.

The waveform model, NRSur7dq4, includes all spin-weighted spherical harmonic modes up to $\ell=4$. The precession frame dynamics and spin evolution of the BHs are also modeled as byproducts. Through a cross-validation study, we show that the mismatches for NRSur7dq4 against NR computed with the Advanced LIGO design sensitivity noise curve are always $\lesssim 8\times 10^{-3}$ at the $95$ percentile level over the mass range $50-200M_{\odot}$. This is at least an order of magnitude improvement over existing waveform models. NRSur7dq4 is made publicly available through the gwsurrogate Blackman et al. Python package, with example evaluation code at Ref. SpE (2018).

For the final BH model, NRSur7dq4Remnant, the $95$th percentile errors are $\mathchar 21016\relax\,5\!\times\!10^{-4}M$ for mass, $\mathchar 21016\relax\,2\!\times\!10^{-3}$ for spin magnitude, $\mathchar 21016\relax\,4\!\times\!10^{-4}\,c$ for kick magnitude. Once again, these are lower than that of existing models by at least an order of magnitude. In addition, we also model the spin and kick directions. Moreover, the GPR methods employed here naturally provide error estimates along with the fitted values. These uncertainty estimates can be incorporated into data analysis applications to marginalize over systematic uncertainties. NRSur7dq4Remnant is made publicly available through the surfinBH Varma et al. (a) Python package, which includes an example evaluation code.

Futher, we provide a Python package, binaryBHexp Varma et al. (b), to visualize the complex precessing dynamics as predicted by these surrogate models Varma et al. (2019c).

VII.1 Future work

In App. A we test the performance of these surrogate models when extrapolated outside their training range to $q=6$. We find that our models become worse at these mass ratios, but are still comparable or better than existing models. Unfortunately, suitable precessing simulations are currently not available for testing at intermediate mass ratios $4<q<6$. In general, we advice caution with extrapolation. A natural improvement of both NRSur7dq4 and NRSur7dq4Remnant is to extend their range of validity with new training simulations at higher mass ratios and spin magnitudes. We note, however, that both these regimes are increasingly expensive to model in NR.

Another important limitation of these models is that they are restricted to the same length as the NR simulations (starting time of $\sim 4300M$ before the peak or about 20 orbits). For LIGO, assuming a starting GW frequency of 20 Hz, the (2, 2) mode of the surrogate is valid for total masses $M\gtrsim 66M_{\odot}$. This number, however, depends on the mass ratio. Fig. 9 shows the mass range of validity of NRSur7dq4 as a function of mass ratio. We compare this with the parameters of the 10 BBH detections seen by LIGO and Virgo in the first two observing runs Abbott et al. (2018a). NRSur7dq4 sufficiently covers the posterior spread of most but not all of these detections, the main limitation being the number of orbits covered by the model. However, see Ref. Kumar et al. (2018) for an example of NR surrogates used in data analysis with GW signals.

A promising avenue to extend the length of the waveforms is to “hybridize” the simulations using PN waveforms in the early inspiral. This approach already was found to be successful for the case of aligned-spin BBH Varma et al. (2019b), but still needs to be generalized to precessing spins. Furthermore, it is not clear if the current length of the NR simulations is sufficient to guarantee good attachment of the PN and NR waveforms for precessing BBH.

Despite these limitations, in their regime of validity, the models presented in the paper are the most accurate models currently available for precessing BBHs. As shown in this paper, our models rival the accuracy of the NR simulations, while being very cheap to evaluate. As more and more BBHs are detected at higher signal-to-noise ratios, fast yet accurate models such as these will contribute to turning GW astronomy into high precision science.

Acknowledgements.

We thank Dan Hemberger, Kevin Barkett, Marissa Walker, Matt Giesler, Nils Deppe, Francois Hebert, Maria Okounkova, and Geoffrey Lovelace for helping carry out the new SpEC simulations used in this work. V.V. and M.S. are supported by the Sherman Fairchild Foundation, and NSF grants PHY–170212 and PHY–1708213 at Caltech. L.E.K. acknowledges support from the Sherman Fairchild Foundation and NSF grant PHY-1606654 at Cornell. S.E.F is partially supported by NSF grant PHY-1806665. 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. Simulations were performed on NSF/NCSA Blue Waters under allocation NSF PRAC–1713694; on the Wheeler cluster at Caltech, which is supported by the Sherman Fairchild Foundation and by Caltech; 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 for building the model were performed on Wheeler.

Appendix A Evaluating surrogates at larger mass ratios

In this Appendix we assess the performance of the NRSur7dq4 and NRSur7dq4Remnant models when evaluated at mass ratio $q=6$. Doing so is effectively an extrapolation because $q=6$ is outside the training range of the surrogates ($q\leq 4$). The surrogate models are compared against 100 NR simulations with $q=6$ and generically precessing spins with magnitudes $\chi_{1},\chi_{2}\leq 0.8$. These simulations have been assigned the identifiers SXS:BBH:2164 - SXS:BBH:2263, and are made publicly available through the SXS public catalog SXS Collaboration .

Figure 10 shows the $q=6$ extrapolation mismatches for NRSur7dq4. Also shown are the mismatches for SEOBNRv3 when compared against the same simulations. The mismatches are computed in the same manner as in Fig. 4(a), which we reproduce here for comparison. The surrogate errors become noticeably worse when extrapolating to $q=6$, but are still much smaller than the corresponding errors for SEOBNRv3.

Fig. 11 shows the performance of NRSur7dq4Remnant when extrapolating to $q=6$. We show the errors when the fits are evaluated using the NR spins at $t\!=\!-100M$ as well as when the spins are specified at the start of the NR simulations. In the latter case, we use the extrapolated dynamics surrogate of NRSur7dq4 to evolve the spins to $t\!=\!-100M$ and then evaluate the fits. We reproduce the training range errors from Fig. 7 for comparison. Also shown are the errors for the existing fitting formulae described in Sec. VI.2 when compared against the same simulations. We find that NRSur7dq4Remnant performs noticeably worse when extrapolated to $q=6$ but is still slightly better than the existing fitting formulae, except for the final spin where the existing fitting formulae perform slightly better.

In general, we find that the NRSur7dq4 and NRSur7dq4Remnant models become worse with extrapolation to $q=6$ but are still better or comparable to existing models. Unfortunately, we do not have enough suitable precessing simulations with $4<q<6$ with which to test at what mass ratio the degradation of these surrogate models becomes significant. We leave these tests, as well as extending the models to larger mass ratios by adding NR simulations, to future work.

Appendix B On the high mismatch tail in NR errors

The histogram of NR errors in Fig. 4(a) shows a significant tail to the right, i.e. at large mismatches. In Sec. VI.1, this tail was attributed to different resolutions of the same NR simulation having different physical parameters, namely the “initial” spins, which are measured at the relaxation time Boyle et al. (2019) after the poorly-resolved junk-radiation transients have settled. In this Appendix we provide some evidence for this claim. Figure 12 shows the maximum mismatch (with a flat noise curve) over points in the sky versus the difference in the relaxation-time dimensionless spins between the two highest resolutions. We refer to the two highest resolutions as HiRes and MedRes, and their corresponding relaxation-time dimensionless spins are denoted by ($\bm{\chi}^{\mathrm{HiRes}}_{1},\bm{\chi}^{\mathrm{HiRes}}_{2}$) and ($\bm{\chi}^{\mathrm{MedRes}}_{1},\bm{\chi}^{\mathrm{MedRes}}_{2}$), respectively. We note that the largest mismatch occurs when the spin difference is largest between the two resolutions. For a significant fraction of the simulations the spins can be different by about $0.1$; for these cases the two resolutions essentially represent two different physical systems, so the difference in waveforms between the two resolutions fails to be a good estimate of the truncation error in the simulations.

Figure 12 suggests that the high NR mismatch tail of Fig. 4(a) is artificially large, and if the two resolutions were to correspond to the same physical system, the tail would be shorter. We test this in Fig. 13, where we compare the surrogate against the MedRes simulations, but use the spins of the MedRes simulation ($\bm{\chi}^{\mathrm{MedRes}}_{1},\bm{\chi}^{\mathrm{MedRes}}_{2}$) to evaluate the surrogate. The surrogate mismatches against the HiRes simulations as well as the NR resolution mismatches (HiRes vs MedRes) are reproduced from Fig. 4(a) for comparison. We note that the surrogate mismatches when compared against the MedRes simulations always lie below $\sim 10^{-2}$ and do not have the high mismatch tail seen for the NR resolution mismatches. In this test, we are treating the surrogate, which is trained on the HiRes simulations, as a proxy for the HiRes dataset. Evaluating the surrogate with the parameters of a MedRes simulation is treated as a proxy for performing the HiRes simulation with the same parameters. Therefore, the green histogram in Fig. 13 can be treated as the “true” resolution error when the parameters of the resolutions are the same. As expected for this case, this estimate of the resolution error agrees with the errors for the surrogate model (red histogram).

Together, Figs. 12 and 13 show that the high NR mismatch tail in Fig. 4(a) is due to the difference in the parameters of the different NR resolutions. We believe this difference arises from spurious initial transients known as “junk radiation”. These transients result from initial data that do not precisely represent a snapshot of a binary that has evolved from $t=-\infty$. The transients quickly leave the simulation domain after about one or two binary orbits. It is computationally expensive to resolve the high spatial and temporal frequencies of the transients, so we typically choose not to resolve these transients at all, and instead we simply discard the initial part of the waveform. Because some of the transients carry energy and angular momentum down the BHs, the masses and spins are modified, so we measure “initial” masses and spins at a relaxation time Boyle et al. (2019) deemed sufficiently late that the transients have decayed away. Because we do not fully resolve the transients, their effect on the masses and spins are not always convergent with resolution.

This issue should ideally be resolved with improved, junk-free initial data (see Ref. Varma et al. (2018) for steps in this direction). In the meantime, we propose a change in how SpEC performs different resolutions for the same simulation. Currently, initial data are constructed by solving the Einstein constraint equations Lovelace et al. (2008); Ossokine et al. (2015b). The same constraint-satisfying initial data are then interpolated onto several grids of different resolution, and Einstein’s equations are evolved on each grid independently. Our proposal is to first evolve the initial data using the high resolution grid until the transients leave the simulation domain, and then interpolate the data at that time onto grids of lower resolution, and evolve Einstein’s equations on these lower-resolution grids independently. This way all resolutions start with the same initial data at a time after transients have decayed away instead of at the start of the simulation, and the masses and spins of the black holes should be convergent.

This proposal is tested in Fig. 14 for the case leading to the largest NR mismatch in Fig. 4(a). We perform the resolution branching at $t\sim 1000M$ after the start of the high resolution simulation. The outer boundary is at $\sim 600M$ and this is sufficient time for junk radiation to leave the simulation domain. We find that the mismatches decrease significantly when the resolution branching is done post-junk, as the resolutions now correspond to the same physical system.