Survey of four precessing waveform models for binary black hole systems

Angular momentum and spin precession are expected to be generic features of a significant fraction of binary black hole systems. As such, it is essential to have waveform models that faithfully incorporate the effects of precession. Here, we assess how well the current state-of-the-art models achieve this for waveform strains constructed only from the $\ell =2$ multipoles. Specifically, we conduct a survey on the faithfulness of the waveform models seobnrv5phm, teobresums, imrphenomtphm, imrphenomxphm to the numerical relativity (NR) surrogate nrsur7dq4 and to NR waveforms from the SXS catalog. The former assessment involves systems with mass ratios up to 6 and dimensionless spins up to 0.8. The latter employs 317 short and 23 long SXS waveforms. For all cases, we use reference inclinations of zero and 90°. We find that all four models become more faithful as the mass ratio approaches unity and when the merger-ringdown portion of the waveforms are excluded. We also uncover a correlation between the coprecessing $(2,\pm 2)$ multipole mismatches and the overall strain mismatch. We additionally find that for high inclinations, precessing $(2,\pm 1)$ multipoles that are more faithful than their $(2,\pm 2)$ counterparts, and comparable in magnitude, improve waveform faithfulness. As a side note, we show that use of uniformly filled parameter spaces may lead to an overestimation of precessing model faithfulness. We conclude our survey with a parameter estimation study in which we inject two precessing SXS waveforms (at low and high masses) and recover the signal with seobnrv5phm, imrphenomtphm, and imrphenomxphm. As a bonus, we present preliminary multidimensional fits to model unfaithfulness for Bayesian model selection in parameter estimation studies.

Figure 1

Typical waveform amplitudes in the frequency domain from the sample of BBHs generated for the survey in Sec. 4a. In particular, we display in blue (orange) a trio of light (heavy) precessing binaries starting from $f_0=37.5\mathrm{Hz}$ (12 Hz) with $Q=1.\overline{1},2,4$ (dot dashed, dotted, solid) at a fiducial luminosity distance of 500 Mpc and initial inclination of ${\vartheta }_{\mathrm{LN},0}=0$ shown against both Advanced LIGO design sensitivity [196] (black solid curve) and LIGO Livingston detector strain sensitivity during O3 [197] (faint gray curve). We used the approximant imrphenomxphm to generate these waveforms. For each case, we have $\max ({\chi }_{\mathrm{p}})=0.8$. The orange and blue vertical bands mark the regions of $f_{\mathrm{peak}}$ for each triplet of waveforms. Since $f_{\mathrm{peak}}$ varies for each set of parameters, we opted to show bands here instead of six additional vertical lines marking individual $f_{\mathrm{peak}}$.

Figure 2

The coverage of the spin space for all our parameter sets shown using the parallel and perpendicular effective spin scalars computed at the initial time $t_0$. In each panel, ${\chi }_{\mathrm{eff}}(t_0)$ of Eq. (8) is plotted in the vertical axis. From left to right, ${\chi }_{\mathrm{p}}(t_0)$ [Eq. (9)], ${\chi }_{\mathrm{p}}^{\mathrm{Gen}}(t_0)$ [Eq. (10)], and ${\chi }_{\perp ,\mathrm{J}}(t_0)$ [Eq. (12)] are plotted in the horizontal axes, respectively. The blue disks correspond to the values of these quantities coming from our discrete parameter set of Sec. 4a (see Table 1). The green squares, orange diamonds, and the black crosses represent the same quantities obtained from our random-uniformly-filled set (Sec. 4c), the short and the long sxs sets (Secs. 5a and 5b), respectively. As ${\chi }_{\mathrm{p}}$ has no dependence on the azimuthal components of the spins, the parameters of the discrete and the short sxs sets yield degenerate values clustered at ${\chi }_{\mathrm{p}}\approx 0.4$, 0.72, 0.8, whereas both ${\chi }_{\mathrm{p}}^{\mathrm{Gen}}$ and ${\chi }_{\perp ,\mathrm{J}}$ cover their respective ranges better since they are ${\varphi }_i$ dependent.

Figure 3

The main result of this section: the distributions of the sky-optimized mismatch, ${\overline{\mathcal{M}}}_{\mathrm{opt}}$, given by Eq. (21) between the models listed in the legend and nrsur7dq4 which we take as a proxy for numerical relativity. The parameters for the waveforms are given in Table 1 with the $Q=6$ subset excluded from the figure and we consider only the $\ell =2$ strain for the mismatches. The left (right) panels display the results for low (high) mass, $M=37.5M_{\odot }$ ($M=150M_{\odot }$), systems, respectively. seob (v5), being the most recent of the four models, shows improvement, especially for heavier systems. The dashed vertical line in both panels marks mismatch of $1-{0.9}^{1/3}\approx 0.035$ translating to an event loss rate of 10% [199, 200].

Figure 4

Merger-ringdown truncated mismatches versus full-length waveform mismatches. The dark-shaded half violins are the distributions of ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$ [Eq. (23)] between the numerical relativity surrogate nrsur7dq4 and the four approximants listed along the bottom horizontal axis of each panel. Similarly, the light-shaded half violins are the distributions of ${\overline{\mathcal{M}}}_{\mathrm{opt}}$ [Eq. (21)], which are shown in the legend as ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{X=\mathrm{MR}}$. These are the same quantities plotted in Fig. 3 as histograms. The left (right) panel contains the data from the light (heavy), $M=37.5M_{\odot }(150M_{\odot })$, binaries. We used the default Gaussian kernel density estimator of the seaborn library [201] to smooth the histograms for the violin plots. The horizontal dashed line marks the mismatch value of $1-{0.9}^{1/3}\approx 0.035$.

Figure 5

Mass ratio separated, merger-ringdown truncated mismatches, ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$, between nrsur7dq4 and the approximants {seob, teob, tphm, xphm} labeled along the bottom axes and plotted in blue, orange, green, red, respectively. For each approximant, we present four half violins representing, from left to right, the ${\vartheta }_{\mathrm{LN},0}=0$ mismatches for the subsets with mass ratios $Q=1.\overline{1},2$, 4, 6 which are shown along the top axes. The $Q=6$ subset is discussed separately in Sec. 4a4 as it involves comparisons with nrsur7dq4 in its extrapolation region. The horizontal dashed line marks the mismatch value of $1-{0.9}^{1/3}\approx 0.035$.

Figure 6

The association of certain peaks in the mismatch distributions with specific regions in the spin parameter space. In each of the four panels, we scatter plot the $\{{\chi }_{\perp ,\mathrm{J}},{\chi }_{\mathrm{eff}}\}$ values of the $Q=4$ subset of the discrete parameter space. The larger disks correspond to points which yield values of ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$ [Eq. (23)] in the rightmost peaks of the histograms plotted in the insets where the thick, vertical dashed lines separate the multiple modes (peaks) of each histogram. The smaller squares mark the complementary $Q=4$ cases. The seob, tphm histograms in the first and third insets correspond to the $Q=4$ half violins from the left panel of Fig. 5, while the rest from the right panel.

Figure 7

The correlation between ${\overline{\mathcal{M}}}_{22,\mathrm{AS}}$ and $\overline{\mathcal{M}}({\vartheta }_{\mathrm{LN},0}=0)$ for the $M=37.5M_{\odot }$ cases, where ${\overline{\mathcal{M}}}_{22,\mathrm{AS}}$ is the mismatch [via Eq. (13)] between nrsur’s coprecessing (2,2) multipole and each AS (2,2) multipole of the four models labeled in the legends where we also show the square of the Pearson correlation coefficient. The upper panels show the mismatches from the combined $Q=1.\overline{1}$ and $Q=2$ subsets, while the lower panels show the results from the $Q=4$ subsets. The cluster of red dots in the upper left corner of the $Q=4$ xphm panel are the cases for which the aforementioned MSA prescription of xphm’s precession dynamics breaks down.

Figure 8

Comparison of the mismatch distributions between face-on (${\vartheta }_{\mathrm{LN},0}=0$) and edge-on (${\vartheta }_{\mathrm{LN},0}=\pi /2$) configurations. In the top row, we display violin plots of ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}({\vartheta }_{\mathrm{LN},0}=0)$ vs ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}({\vartheta }_{\mathrm{LN},0}=\pi /2)$ for the $M=37.5M_{\odot }$ ($M=150M_{\odot }$) sets in the left (right) panels. Likewise, the bottom row shows ${\overline{\mathcal{M}}}_{\mathrm{opt}}({\vartheta }_{\mathrm{LN},0}=0)$ vs ${\overline{\mathcal{M}}}_{\mathrm{opt}}({\vartheta }_{\mathrm{LN},0}=\pi /2)$ for the same sets. Note that we exclude the $Q=6$ subset here, but we include it in Fig. 9.

Figure 9

Comparison of the mismatch distributions between “face-on” (${\vartheta }_{\mathrm{LN},0}=0$) and “edge-on” (${\vartheta }_{\mathrm{LN},0}=\pi /2$) configurations separated by mass ratio. In particular, we present results here only for the MR-truncated, sky-angle optimized mismatch ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$. The top row displays results for the light ($M=37.5M_{\odot }$) BBH set, and the bottom row the heavy ($M=150M_{\odot }$) BBHs.

Figure 10

Mismatches between nrsur7dq4 and seob, teob, tphm, xphm for the uniformly filled parameter space of 1000 cases. Upper panels compare ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$ with ${\overline{\mathcal{M}}}_{\mathrm{opt}}$ at inclinations of ${\vartheta }_{\mathrm{LN},0}=0$ (left) and $\pi /2$ (right) and the lower panels compare ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}({\vartheta }_{\mathrm{LN},0}=0)$ with ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}({\vartheta }_{\mathrm{LN},0}=\pi /2)$ (left) and ${\overline{\mathcal{M}}}_{\mathrm{opt}}({\vartheta }_{\mathrm{LN},0}=0)$ with ${\overline{\mathcal{M}}}_{\mathrm{opt}}({\vartheta }_{\mathrm{LN},0}=\pi /2)$ (right).

Figure 11

Comparison of model mismatches to nrsur7dq4 between the discrete parameter set and the random-uniformly-filled one. The left half violins are the mismatch distributions from the entire $M=\{37.5M_{\odot },150M_{\odot }\},Q=\{1.\overline{1},2,4\}$ discrete set while the right ones are the mismatch distributions from the 1000-case random-uniformly-filled set. In the top (bottom) row, we plot ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$ (${\overline{\mathcal{M}}}_{\mathrm{opt}}$) at inclinations of 0 (left panels) and $\pi /2$ (right panels).

Figure 12

Mismatches defined by Eq. (21) between the four approximants and the 317 short sxs waveforms shown at an inclination of ${\vartheta }_{\mathrm{LN},0}=0$ for a direct comparison with Fig. 3. The left and right columns of panels correspond to mismatches from the light ($M=37.5M_{\odot }$) and heavy ($M=150M_{\odot }$) binaries.

Figure 13

Merger-ringdown truncated and full mismatches between the four approximants and the 317 short sxs waveforms shown for inclinations of 0 and $\pi /2$. The left (right) $2\times 2$ blocks of figures correspond to light, $M=37.5M_{\odot }$ (heavy, $M=150M_{\odot }$) binaries. The top panels compare ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$ with ${\overline{\mathcal{M}}}_{\mathrm{opt}}$ for ${\vartheta }_{\mathrm{LN},0}=0$ and $\pi /2$ and the bottom panels compare ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}({\vartheta }_{\mathrm{LN},0}=0)$ with ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}({\vartheta }_{\mathrm{LN},0}=\pi /2)$ and ${\overline{\mathcal{M}}}_{\mathrm{opt}}({\vartheta }_{\mathrm{LN},0}=0)$ with ${\overline{\mathcal{M}}}_{\mathrm{opt}}({\vartheta }_{\mathrm{LN},0}=\pi /2)$.

Figure 14

Inspiral-only mismatches at ${\vartheta }_{\mathrm{LN},0}=0$, i.e., ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}({\vartheta }_{\mathrm{LN},0}=0)$ between our entire $Q\lesssim 4$ short sxs set of waveforms and seob, teob, tphm, xphm plotted as a function of ${\chi }_{\perp ,\mathrm{J}}$ [Eq. (12)]. Since ${\chi }_{\perp ,\mathrm{J}}$ is maximized for ${\theta }_{1,2}\to \pi /2$, we observe increasing mismatches as the spins become more planar, i.e., for stronger precession.

Figure 15

Sky-optimized mismatches [Eq. (21)] between {seob, teob, tphm, xphm} and the 317 short sxs waveforms shown for inclinations of ${\vartheta }_{\mathrm{LN},0}=0$ and $\pi /2$, separated into three subsets in terms of the mass ratio $Q$. The top (bottom) panels show the mismatch distributions for binaries with total mass $M=37.5M_{\odot }(150M_{\odot })$.

Figure 16

The correlation between ${\overline{\mathcal{M}}}_{22,\mathrm{AS}}$ and $\overline{\mathcal{M}}({\vartheta }_{\mathrm{LN},0}=0)$ for the $M=37.5M_{\odot },Q\lesssim 2$ cases of our short sxs set. ${\overline{\mathcal{M}}}_{22,\mathrm{AS}}$ is the mismatch [Eq. (13)] between the sxs coprecessing (2,2) multipole and each AS (2,2) multipole of the four models. Each legend also gives the $R^2$ value for the log-log datasets. This figure is akin to the top panel of Fig. 7.

Figure 17

Trace plots of the mismatches between {seob, teob, tphm, xphm} and numerical relativity waveforms from 23 long sxs simulations as a function of the total binary mass. We computed the mismatches at the values of $M=\{37.5,75,112.5,150,187.5,225\}{M}_{\odot }$ (marked by the faint vertical lines) and linearly connected the data points. The top (bottom) two rows show ${\overline{\mathcal{M}}}_{\mathrm{opt}}^{\text{noMR}}$ (${\overline{\mathcal{M}}}_{\mathrm{opt}}$) with odd (even) rows corresponding to ${\vartheta }_{\mathrm{LN},0}=0$ (${\vartheta }_{\mathrm{LN},0}=\pi /2$) results. The four columns represent, from left to right, the results for seob, teob, tphm, xphm.

Figure 18

Recovery performances of seob, tphm, and xphm for a synthetic signal generated by a “zero-noise” injection of the sxs simulation 0050 into the O4a LIGO Hanford-Livingston network. In each panel, we present the results from both the light- and the heavy-mass runs as the darker and lighter halves of each violin, respectively. Shown are the recovered posteriors for $\{{\mathcal{M}}_{\mathrm{c}}^{\det },M^{\det },q,{\chi }_{\mathrm{eff}},d_{\mathrm{L}},{\chi }_1,{\theta }_1,{\chi }_2,{\theta }_2,{\vartheta }_{\mathrm{LN},0}\}$ with the dashed horizontal lines marking the injected values. For $\{{\mathcal{M}}_{\mathrm{c}}^{\det },M^{\det },d_{\mathrm{L}}\}$, we present posteriors for the relative difference between the injected and the recovered values, i.e., ${\mathrm{\Delta }}_{\mathrm{rel}}X≔(X_{\mathrm{r}}-X_{\mathrm{inj}})/X_{\mathrm{r}}$. The first panel, i.e., ${\mathrm{\Delta }}_{\mathrm{rel}}{\mathcal{M}}_{\mathrm{c}}^{\det }$, is somewhat hard to read because the models produced very narrow posteriors for the light-mass run, but much wider ones for the heavy one.

Figure 19

Same as Fig. 18, but for the sxs simulation 0628.

Figure 20

Posterior distributions recovered by tphm and xphm for $\{{\chi }_{\mathrm{eff}},{\chi }_{\mathrm{p}},{\chi }_{\mathrm{p}}^{\mathrm{Gen}},{\chi }_{\perp ,\mathrm{J}}\}$ for the injection of the sxs simulation 0050 with $M^{\det }=56.5M_{\odot }$. The black dashed distribution is the prior for each quantity. The orange line marks the injected value. The color coding of each model is consistent with the entire article.

Figure 21

Same as Fig. 20, but for the heavy 0050 PE run with an injected value of $M^{\det }=150M_{\odot }$. Note that the vertical scale is different from Fig. 20, which is why the histograms for the priors may give the illusion of looking different, when they are in fact identical.

Figure 22

Same as Fig. 20, but for the light 0628 run with an injected value of $M^{\det }=67M_{\odot }$.

Figure 23

Same as Fig. 20, but for the heavy 0628 run with an injected value of $M^{\det }=150M_{\odot }$.

Figure 24

Contour plots of the relative difference between the validation data (red points) and the third teob fit of Table 5, which is obtained from Eq. (a1) with $\{n_x,n_y,n_z\}=\{7,3,2\}$. This particular fit is a function of ${\chi }_{\perp ,\mathrm{J}},{\chi }_{\mathrm{eff}}$ and $\eta$ and represents the $M=37.5M_{\odot }$, ${\log }_{10}[{\overline{\mathcal{M}}}_{\mathrm{opt}}({\vartheta }_{\mathrm{LN},0}=0)]$ data, where ${\overline{\mathcal{M}}}_{\mathrm{opt}}$ is the sky-optimized mismatch between nrsur7dq4 and teob. The three panels from left to right show contour plots of the relative differences separated by mass ratios of $Q=1.\overline{1},2$, 4. The black dots are the data points used to construct the fit and are not covered by the contour plot shown here which only applies to the red dots.

Figure 25

Comparison of the faithfulness of seobnrv5phm with that of seobnrv4phm for the 1000-case random-uniformly-filled parameter space described in Sec. 4c. The baseline model used to measure faithfulness is nrsur7dq4.

Figure 26

Comparison of the faithfulness of the MSA implementation of imrphenomxphm with that of its more recent spintaylor implementation (ST in the figure) and the new model imrphenomxo4a. The parameters used to generate the waveforms are those of the random-uniformly-filled set of Sec. 4c. The baseline model used to measure faithfulness is nrsur7dq4.

Figure 27

Same comparison as Fig. 26, but for the discrete parameter set of Sec. 4a. Note the disappearance of a secondary mode containing high-mismatch values for the spintaylor version of xphm and imrphenomxo4a.

Figure 28

Timing benchmarks for seob, teob, tphm, and xphm. Each panel gives results from runs with key parameters displayed above with $Q=2$, ${\vartheta }_{\mathrm{LN},0}=0$, $f_0=20\mathrm{Hz}$, $d_{\mathrm{L}}=500\mathrm{Mpc}$. Each individual case was run 1000 times per model in a virtual environment on an AMD EPYC 7453 server.