A Surrogate model of gravitational waveforms from numerical relativity simulations of precessing binary black hole mergers

We present the first surrogate model for gravitational waveforms from the coalescence of precessing binary black holes. We call this surrogate model NRSur4d2s. Our methodology significantly extends recently introduced reduced-order and surrogate modeling techniques, and is capable of directly modeling numerical relativity waveforms without introducing phenomenological assumptions or approximations to general relativity. Motivated by GW150914, LIGO’s first detection of gravitational waves from merging black holes, the model is built from a set of 276 numerical relativity (NR) simulations with mass ratios $q\le 2$, dimensionless spin magnitudes up to 0.8, and the restriction that the initial spin of the smaller black hole lies along the axis of orbital angular momentum. It produces waveforms which begin $\sim 30$ gravitational wave cycles before merger and continue through ringdown, and which contain the effects of precession as well as all $\ell \in \{2,3\}$ spin-weighted spherical-harmonic modes. We perform cross-validation studies to compare the model to NR waveforms not used to build the model and find a better agreement within the parameter range of the model than other, state-of-the-art precessing waveform models, with typical mismatches of ${10}^{-3}$. We also construct a frequency domain surrogate model (called NRSur4d2s_FDROM) which can be evaluated in 50 ms and is suitable for performing parameter estimation studies on gravitational wave detections similar to GW150914.

Figure 1

A schematic of the method for building a surrogate model for a function $X(t;\mathit{\lambda })$. The red dotted lines show $X(t)$ evaluated at a selected set of greedy parameters ${\mathbf{\Lambda }}_i$ used to build a linear basis, and the blue dots show the associated empirical nodes in time from which $X_S(t;\mathit{\lambda })$ can be reconstructed by interpolation with high accuracy. The blue lines indicate fits for $X(t;\mathit{\lambda })$ as a function of $\mathit{\lambda }$ at each of the empirical time nodes. The cyan dot shows a generic parameter ${\mathit{\lambda }}_0$ that is not in the set of greedy parameters. To compute $X_S(t;{\mathit{\lambda }}_0)$, each fit is evaluated at ${\mathit{\lambda }}_0$ (the yellow diamonds), and then the empirical interpolant is used to construct $X_S(t;{\mathit{\lambda }}_0)$ at arbitrary times (the dotted black line).

Figure 2

Surrogate workflow. A greedy “PN-sampler” selects the most informative parameter values $\{{\mathbf{\Lambda }}_i{\}}_{i=1}^N$ for a fixed parametric and temporal range. For each selected value ${\mathbf{\Lambda }}_i$, SpEC generates a gravitational waveform. A surrogate model building algorithm (cf. Fig. 1) is applied to a set of suitably aligned and decomposed (cf. Fig. 6) numerical relativity waveforms thereby producing a trial surrogate. A handful of validation tests are performed to assess the surrogate’s quality. If the surrogate performs poorly for some parameter values, one could produce additional numerical relativity waveforms near those values and rebuild a more accurate surrogate.

Figure 3

Diagram of the 4 spin components in the 5D parameter subspace. We attempt to obtain ${\varphi }_{\chi }=0$ at $t_0=4500M$ before peak amplitude, but in practice the NR simulations have arbitrary values of ${\varphi }_{\chi }$.

Figure 4

A “triangle plot” showing all possible two-dimensional projections and one-dimensional histograms of the greedy parameters $G$ selected by the procedure of Sec. 3b. These are the parameters used for the numerical relativity simulations. Made using the Python package corner.py [65].

Figure 5

Maximum surrogate errors found during the second greedy algorithm (see Sec. 3b) for determining ${\mathbf{\Lambda }}_i$ using trial PN surrogates. The noise is due to the random resampling, as well as the possibility of the parameter space fits becoming worse by adding a data point. The finite order of the fits leads to an error floor of ${10}^{-3}$, so we keep and perform NR simulations for only the first 300 greedy parameters.

Figure 6

Waveform decomposition schematic. A series of decompositions are applied to a set of NR waveform modes $\{h^{\ell m}\}$ yielding easier-to-approximate waveform data pieces (shown as cyan ellipses) for which we ultimately fit. Two types of objects are shown: time series data as an ellipse and operators/maps as rectangles. A red outlining border identifies an object which uses a modeling approximation which will not go away with additional NR waveforms. These decomposition errors are quantified and shown to be smaller than other sources of error in Sec. 6. An additional source of error that will not converge away with more NR waveforms results from the assumption that each data piece transform in a simplistic way with changes of ${\varphi }_{\chi }$.

Figure 7

Waveform modes in the inertial frame for SXS:BBH:0338 with $q=2$, $|{\overrightarrow{\chi }}_1|=0.8$, ${\theta }_{\chi }=1.505$, ${\varphi }_{\chi }=-1.041$ and ${\chi }_2^z=0.8$. For precessing systems, all $\ell =2$ modes contain significant power in the inertial frame. The NR waveform is aligned to have the canonical orientation at $t=t_0$.

Figure 8

Waveform modes in the coprecessing frame for SXS:BBH:0338. The mode power hierarchy is now the same as for a nonprecessing waveform, with the $(2,\pm 2)$ modes dominating, but small effects of precession are still present in the mode amplitudes and phases. The amplitudes of the $(2,\pm 2)$ modes have small nearly opposite oscillations.

Figure 9

Top: Quaternion $q$ representing the time-dependent rotation from the coprecessing frame to the inertial frame (solid lines) and the filtered quaternion $q_{\min -\mathrm{filt}}$ (dashed lines) for SXS:BBH:0338. Bottom: Differences between the filtered and unfiltered quaternions. This difference results in an error when reconstructing the waveform in the inertial frame, contributing to a “decomposition” error in the surrogate model.

Figure 10

Phases ${\phi }_p$ (top) and ${\phi }_d$ (middle) for SXS:BBH:0338. These phases represent the total amount of precession and the instantaneous direction of precession respectively. Shown are phases computed from the unfiltered coprecessing quaternion (thick orange lines) and the filtered quaternion (thin black lines). The orbital time scale oscillation in ${\phi }_d$ is suppressed after filtering. Bottom: Differences between the filtered and unfiltered phases.

Figure 11

For the real-valued oscillatory components $X$ such as $\mathrm{Im}[{\tilde{h}}^{2,0}]$, we perform a Hilbert transform to obtain a complex signal $H(X)$ and extract an amplitude and phase. The dashed green line shows the imaginary part of $H(X)$.

Figure 12

Top: waveform differences $\delta h(t)$ investigating the removal of the ${\varphi }_{\chi }$ dependence on the waveform. Each colored band includes waveforms compared to SXS:BBH:0346 and SXS:BBH:0346 for several different values of ${\varphi }_{\chi }$. Before making any adjustment, the errors (${\varphi }_{\chi }$ differences) are large. After adjusting, the errors (${\varphi }_{\chi }$ adjusted) are comparable to resolution errors during the inspiral but grow large at merger. The decomposition errors are negligible. Bottom: differences in ${\phi }_{-}^{2,2}$. Our analytic approximation to remove the effect of ${\varphi }_{\chi }$ on the waveform does not affect ${\phi }_{-}^{2,2}$, but here we see that the orbital phase at merger can vary by nearly a radian for different values of ${\varphi }_{\chi }$, which is the most significant contribution to the ${\varphi }_{\chi }$ adjusted errors in the top figure.

Figure 13

Histograms of time domain waveform errors $\mathcal{E}$ relevant to the surrogate. Equal areas under the curves correspond to equal numbers of cases, and the curves are normalized such that the total area under each curve when integrated over ${\log }_{10}(\mathcal{E})$ is 1. Solid black: the resolution error comparing the highest and second highest resolution NR waveforms. Dotted brown: the error intrinsic to the surrogate’s waveform decomposition. Filtering out nutation in the quaternions and neglecting the small but nonzero $\delta q_z$ due to discrete time sampling leads to errors in the reconstructed waveforms. These errors are nearly zero for nonprecessing cases, and even for precessing cases they are smaller than the resolution errors. Thin solid blue: the errors when the full surrogate attempts to reproduce the set of waveforms from which it was built. Dashed purple: the errors when trial surrogates attempt to reproduce NR waveforms that were omitted during the surrogate construction.

Figure 14

Parameter dependence of the error $\mathcal{E}[h,h_S]$ when reproducing the set of NR waveforms with the surrogate. Diagonal: For each parameter plot, the black dots label the (parameter value, $\mathcal{E}[h,h_S]$) pairs. Off diagonal: For each pair of parameters, we show the 2D projection of parameters as in Fig. 4 while varying the color and size of the point based on the error $\mathcal{E}[h,h_S]$. Points are placed in order of increasing error, to ensure the small yellow points with large errors are visible. Larger spin magnitudes, especially for precessing spin configurations, correlate with larger errors.

Figure 15

Errors ${\mathcal{E}}_X$ showing the error contribution of a single surrogate component $X$.

Figure 16

Error contributions $\delta h(t)$ of those waveform data pieces $X$ that have the largest error $\mathcal{E}[h,h_X]$ for a selected simulation: ID 79. To compute the error, the NR waveform is decomposed into the surrogate components, and component $X$ is replaced with its surrogate evaluation. The waveform is then reconstructed, and $\delta h(t)$ is computed from Eq. (16). The solid black curve is given by Eq. (24). The dashed curve is the error in ${\phi }_p$, which is the dominant error in modeling the precession, and the dominant error source during the inspiral. The dotted curve is the error in a quantity similar to twice the orbital phase and becomes the dominant error source during the merger and ringdown. The contribution from errors in the other waveform data pieces is smaller, as shown in Table 3.

Figure 17

Mismatches, computed using a flat noise curve, versus the highest resolution NR waveforms. Histograms are normalized to show the error fraction per log-mismatch, such that the area under each curve is the same. A sufficient but not necessary condition for a mismatch to have a negligible effect is that the signal-to-noise ratio (SNR) lies below the limiting SNR ${\rho }_{*}=1/\sqrt{2\text{Mismatch}}$ given on the top axis [72]. Top: All modes available to each waveform model are included, and the NR waveforms use all $\ell \le 5$ modes. Middle: All coprecessing-frame modes other than $(2,\pm 2)$ are set to zero in all waveforms. Bottom: All coprecessing-frame modes other than $(2,\pm 1)$ and $(2,\pm 2)$ are set to zero in all waveforms. These restricted mode studies are done to compare more directly with IMRPhenomPv2 and SEOBNRv3, which retain the coprecessing-frame modes of the middle and bottom panels respectively.

Figure 18

Median (lower curves, circles) and 95th percentile (upper curves, triangles) mismatches for various total masses $M$ using the advanced LIGO design sensitivity. The median NR resolution mismatches are all below $2\times {10}^{-5}$. The “surrogate” mismatches shown here are “validation” errors described in Sec. 6a.

Figure 19

Max, mean and medians of the distributions of $\mathcal{E}$ when building a surrogate using the first $N$ waveforms and a validation set consisting the remaining $200-N$ waveforms.

Figure 20

Comparing surrogate evaluations to three NR waveforms (top, middle and bottom plots) with spins outside the 5D parameter subspace. For each case and each of the four parameter mappings, the surrogate model waveform error is shown. In all cases, the “PN” mapping performs well at the very start of the waveform and the “drop” mapping performs poorly, but there is no clear overall best mapping. Surrogate modeling errors contribute to the difficulty in assessing the quality of the mappings.

Figure 21

Cumulative error distributions of the frequency domain NRSur4d2s_FDROM surrogate waveforms compared to the time domain NRSur4d2s surrogate waveforms transformed to the frequency domain, evaluated for randomly chosen uniformly distributed parameters. The curves indicate the fraction of errors at least as large as the indicated error. The NRSur4d2s_FDROM output converges to the FFT of the NRSur4d2s output as the grid size is increased.

Figure 22

Fit residuals for the second empirical node of ${\phi }_{-}^{2,2}$ at $t=-806.5M$. Blue dashed: The maximum fit residual using all data. Thin grey lines: Maximum validation residual for individual trials. Thick black line: The root mean square (rms) of the validation residuals for $K=50$ trials. It takes its minimum value at $m=30$, which determines the number of fit coefficients to use for this node in the model. Red: The rms of the training residuals for $K=50$ trials.

Figure 23

Plots of $X(t;\lambda )$ for our toy problem evaluated at the smallest and largest values of $\lambda$ in the training set.

Figure 24

Maximum projection errors of all three reduced bases (see text for a description) versus the size of the basis.

Figure 25

Projection errors, measured in the $L_2$ norm for the three reduced bases described in the text, computed for test data generated from 1000 randomly selected parameters $\lambda$ in [1, 20]. The corresponding colored lines indicate the smallest projection errors on the training sets shown in Fig. 24. The errors for “Greedy, $L_{\infty }$” and “SVD” lie nearly on top of each other. However, the maximum projection error implied by SVD (purple line) underestimates the true errors (dots) by an order of magnitude.

Figure 26

The tenth basis element as a function of $t$ from the three reduced bases elements described in the text. The training data used is given by the parametrized function in (b2) and exhibits relatively large amplitude fluctuations. Whereas the top two plots show significant noise in the basis element, the SVD method smooths away, almost completely, the uncorrelated stochastic features to generate a basis element that is smooth in $t$.