High-accuracy waveforms for black hole-neutron star systems with spinning black holes

The availability of accurate numerical waveforms is an important requirement for the creation and calibration of reliable waveform models for gravitational wave astrophysics. For black hole-neutron star binaries (BHNS), very few accurate waveforms are however publicly available. Most recent models are calibrated to a large number of older simulations with good parameter space coverage for low-spin nonprecessing binaries but limited accuracy, and a much smaller number of longer, more recent simulations limited to nonspinning black holes. In this paper, we present long, accurate numerical waveforms for three new systems that include rapidly spinning black holes, and one precessing configuration. We study in detail the accuracy of the simulations, and in particular perform for the first time in the context of BHNS binaries a detailed comparison of waveform extrapolation methods to the results of Cauchy characteristic extraction. The new waveforms have $<0.1\mathrm{rad}$ phase errors during inspiral, rising to $\sim (0.2–0.4)\mathrm{rad}$ errors at merger, and $\lesssim 1\%$ error in their amplitude. We compute the faithfulness of recent analytical models to these numerical results for the late inspiral and merger phases covered by the numerical simulations, and find that models specifically designed for BHNS binaries perform well (faithfulness $F>0.99$) for binaries seen face on. For edge-on observations, particularly for precessing systems, disagreements between models and simulations increase, and models that include precession and/or higher-order modes start to perform better than BHNS models that currently lack these features.

Figure 1

Simulation Q4S9 at a time when $\sim 30\%$ of the neutron star mass has been accreted by the black hole. Left: matter with density above $6\times {10}^7\mathrm{g}/{\mathrm{cm}}^3$, and apparent horizon of the black hole. Center: finite difference grid, showing both mesh refinement and the fact that the grid does not cover vacuum regions. Right: inner region of the pseudospectral grid, showing mesh refinement close to the black hole and in regions where dense matter is present. We only show points below the orbital plane and within $\sim 15M$ of the black hole’s center. The pseudospectral grid extends to much larger distances ($500M$), using spherical shells not shown on the figure. All three figures are taken from our lowest resolution simulation.

Figure 2

Distribution of the BHNS simulations in the SXS catalog projected in the $Q$-${\chi }_{\parallel }^{\mathrm{BH}}$ plane. Simulations from [16] are red dots (SXS2019), while simulations from this work are green dots (SXS2020). The region of parameter space covered by 134 short waveforms from [10] is shown in gray (LEA). For context, we also show the regions of parameter space where $(1.2–1.6)M_{\odot }$ neutron stars satisfying the equations of state constraints of [27] always disrupt, disrupt for some equations of state and/or neutron star masses only, or never disrupt, according to the fitting formula from [34].

Figure 3

GW strain for simulations Q3S9 (left), Q4S9 (center) and Q3S75p (right). Each plot shows results for an observer in the direction of the total angular momentum (top, “face on” for nonprecessing system) and from a direction that, at $t=0$, is orthogonal to the total angular momentum and within the orbital plane of the binary (bottom, “edge on” at $t=0$). We show both polarizations $h_{+}$ (solid lines) and $h_{\times }$ (dashed lines) and shift the average value of the strain for readability. The edge-on signals clearly show the impact of mass asymmetry (oscillations in the maximum of the strain, due to a significant $l=3$, $m=3$ mode) and, for Q3S75p, precession (slow oscillation in the amplitude of the signal, nonzero value of $h_{\times }$). Q3S75p goes through slightly more than half a precession cycle over the course of the simulation.

Figure 4

Estimates of the phase error for the (2,2) mode of the GW signal in each simulation. We include the impact of mass losses (${\mathrm{\Phi }}_{\text{dM}}$), extrapolation error (${\mathrm{\Phi }}_{\text{ext}}$), and grid resolution (${\mathrm{\Phi }}_{dis}$), following the methods described in [16] and Sec. 3b. The vertical dot-dashed curves indicate $t_{\text{peak}}$.

Figure 5

Estimates of the relative error in the amplitude of the (2,2) mode of the GW signal for each simulation. From top to bottom, we show results for Q3S9, Q4S9 and Q3S75p. Dashed curves show comparisons between the amplitude obtained with different numerical resolutions ($L0$, $L1$, $L2$), while the solid curves show comparisons between the amplitude for different extrapolation orders ($N2$, $N3$, $N4$). The vertical dot-dashed curves indicate $t_{\text{peak}}$. We clearly see that extrapolation is the main source of error when estimating the GW amplitude. Amplitude errors are $\lesssim 1\%$, except for simulation Q4S9 at the time of merger.

Figure 6

Difference in the phase of the (2,2) mode for various estimates of ${\mathrm{\Psi }}_4$ in system Q3S9, after application of a time and phase shift minimizing the phase error in the window $t\in [500,2300]$. We show phase differences between our numerical resolutions ($L0$,$L1$,$L2$), between different extrapolation orders ($N2$,$N3$,$N4$, for the highest resolution simulation), and between the highest resolution simulation with $N=3$ and the waveforms obtained using CCE from data extracted at $R=200M$ (Cce2), $R=300M$ (Cce3) and $R=400M$ (Cce4). We see very good agreement between the various CCE waveforms, and between $N=2$ extrapolation (solid blue curve) and CCE.

Figure 7

Same as Fig. 6, but for the (2,2) mode of ${\mathrm{\Psi }}_4$ in the Q4S9 configuration. We find very similar results, including very good agreement between all CCE results and the $N=2$ extrapolation results.

Figure 8

Same as Fig. 6, but for the (2,2) mode of ${\mathrm{\Psi }}_4$ in the Q3S75p configuration.

Figure 9

Imaginary part of the (2,2) mode of the GW strain ($h_{\times }$ polarization) using CCE extraction from $R=400M$ and extrapolation with $N=3$. The dotted line show the amplitude of the signal ($\sqrt{|h{|}^2}$). In the CCE results, the average value of the strain is shifted away from 0 during the first $\mathrm{\Delta }t\sim 1000M$.

Figure 10

Same as Fig. 6, but for the (2,2) mode of the GW strain. For all CCE waveforms, we subtract a constant value from the strain (see text). Early time oscillation in the CCE results are due to the initial drift in the average value of the GW strain.

Figure 11

Same as Fig. 10, but after applying a highpass filter to all signals to remove the drift in the average value of the strain (tenth-order Butterworth filter with critical frequency of $0.002M^{-1}\approx 72\mathrm{Hz}$). We note that this filter leads to changes in the phase of the signal that are larger than any of the errors displayed here; this plot indicates that most of the error observer in Fig. 10 is indeed due to a slow drift in the average value of the strain, but the filtering does not provide us with a better waveform template.

Figure 12

Same as Fig. 10, but for the (3,3) mode of the GW strain. The CCE method provides much better result for this mode, with great self-consistency between CCE results using different extraction radii, and small differences between CCE and $N=2$ extrapolation.

Figure 13

Amplitude of the (2,2) mode of the GW strain at a retarded time $t=3000M$ for simulation Q3S9-L2. We show the values estimated at finite radii (normalized to $R_{\min }=100M$), and the extrapolation polynomials of order $N=2$, 3, 4, 5. From left to right, we show extrapolation functions fitted to data at all 20 radii used by our simulation, the 16 smallest radii of that set, and the 12 smallest radii of that set (only the $N=3$ results is distinguishable from the others on the middle panel). While low-order extrapolation provides consistent and visually reasonable results, higher order extrapolation is less reliable. This can be contrasted with results in vacuum simulations, where results converge to a well-defined answer as $N$ increases.

Figure 14

Top: amplitude of the complex strain $h=h_{+}+i{h}_{\times }$ for system Q3S9 seen face on. We show results for our low- and high-resolution simulations, for the two BHNS models, as well as for one tidal model that does not include neutron star disruption and one BBH model that includes precession and higher-order modes. Bottom: phase difference between the models and the high-resolution numerical result. All signals are aligned by minimizing the phase difference in the time interval $t\in [500,2000]M$.

Figure 15

Same as Fig. 14, but for system Q3S75p.