Surrogate models for precessing binary black hole simulations with unequal masses

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; Abbott et al. 2017a; Abbott et al. 2016b; Abbott et al. 2017b; Abbott et al. 2017c; Abbott et al. 2017d; Abbott et al. 2018a are becoming routine Abbott et al. 2018b; Abbott et al. 2018c. 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; Rezzolla et al. 2008b; 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:

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.

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:

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.

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

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

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:

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.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)];

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

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

IV.3 Parametric fits

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 $$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

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\,76$ simulations each. For each group, we build a trial surrogate using the $\mathchar 21016\,1452$ remaining training simulations and test against these $\mathchar 21016\,76$ validation ones, which may include points on the boundary of the training set.

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

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 Note 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. 5 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.

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\,5\!\times\!10^{-4}M$ for mass, $\mathchar 21016\,2\!\times\!10^{-3}$ for spin magnitude, $\mathchar 21016\,4\!\times\!10^{-3}$ radians for spin direction, $\mathchar 21016\,4\!\times\!10^{-4}\,c$ for kick magnitude, and $\mathchar 21016\,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\,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\,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; Abbott et al. 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.

For the final BH model, NRSur7dq4Remnant, the $95$th percentile errors are $\mathchar 21016\,5\!\times\!10^{-4}M$ for mass, $\mathchar 21016\,2\!\times\!10^{-3}$ for spin magnitude, $\mathchar 21016\,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.

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.

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

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.

References

References