Adiabatic gravitational waveform model for compact objects undergoing quasicircular inspirals into rotating massive black holes

We present bhpwave: a new python-based, open-source tool for generating the gravitational waveforms of stellar-mass compact objects undergoing quasicircular inspirals into rotating massive black holes. These binaries, known as extreme-mass-ratio inspirals (EMRIs), are exciting mHz gravitational wave sources for future space-based detectors such as the Laser Interferometer Space Antenna (LISA). Relativistic models of EMRI gravitational wave signals are necessary to unlock the full scientific potential of mHz detectors, yet few open-source EMRI waveform models exist. Thus we built bhpwave, which uses the adiabatic approximation from black hole perturbation theory to rapidly construct gravitational waveforms based on the leading-order inspiral dynamics of the binary. In this work, we present the theoretical and numerical foundations underpinning bhpwave. We also demonstrate how bhpwave can be used to assess the impact of EMRI modeling errors on LISA gravitational wave data analysis. In particular, we find that for retrograde orbits and slowly spinning black holes we can mismodel the gravitational wave phasing by as much as $\sim 10$ radians without significantly biasing EMRI parameter estimation.

Figure 1

A $32\times 64$ grid in $(\alpha ,\chi )$ mapped to the $(\widehat{a},\widehat{\mathrm{\Omega }})$ domain. Each dot represents a sampled point in the parameter space. These points are equally spaced in $\alpha$ and in $\chi$, but cluster around $\widehat{\mathrm{\Omega }}={\widehat{\mathrm{\Omega }}}_{\mathrm{ISCO}}$ and $\widehat{a}={\widehat{a}}_{\max }$. The shading of each point reflects the relative density of points in the $(\widehat{a},\widehat{\mathrm{\Omega }})$ domain, while the dashed lines represent the domain boundaries. Note that this is the downsampled version of the numerical grid used for interpolating flux and trajectory information within bhpwave.

Figure 2

The energy flux ${\mathcal{F}}_E$ (left panel), (rescaled) orbital phase $\check{\mathrm{\Phi }}$ (middle panel), and (rescaled) orbital time to merger $\check{t}$ as a function of (dimensionless) orbital frequency $\widehat{\mathrm{\Omega }}$. Different colors (shades) correspond with different values of $\widehat{a}$. Therefore, lower spins terminate at lower ISCO frequencies. Note that the spin has a very small effect on the flux.

Figure 3

The normalized energy flux ${\mathcal{F}}_N$ (left panel), normalized orbital phase ${\check{\mathrm{\Phi }}}_N$ (middle panel), and normalized orbital time to merger ${\check{t}}_N$ as a function of normalized distance from the ISCO frequency ${\widehat{\mathrm{\Omega }}}_{\mathrm{ISCO}}$ in frequency space. We rescale each vertical axis by the leading post-Newtonian coefficient. Thus all curves asymptote to one at low frequencies.

Figure 4

The normalized energy flux ${\mathcal{F}}_N$ (left panel), normalized orbital phase ${\check{\mathrm{\Phi }}}_N$ (middle panel), and normalized orbital time to merger ${\check{t}}_N$ as a function of the final grid parameters $\alpha$ and $\chi$. The values of $\chi$ correspond to the same values of $\widehat{a}$ used in Fig. 2.

Figure 5

The logarithm of the maximized phase difference $\log ||\delta {\mathrm{\Phi }}^I|{|}_{\widehat{a},\widehat{\mathrm{\Omega }}}$ as described in Sec. 3b for different observation periods $T_{\mathrm{obs}}$. We can see that for any observation period, all the way up to 25 years, our phase spline maintains subradian phase accuracy when compared to direct numerical integration of the orbital phase.

Figure 6

Waveform amplitudes $A_{\ell ,m}$ as a function of the radius of the orbit $r_0/M$. The top panel plots $A_{\ell =2,m=2}$ for different Kerr spin parameters $\widehat{a}\ge 0.9$. The middle and bottom panels plot $A_{\ell ,m=2}$ and $A_{\ell ,m=\ell }$, respectively, for $\widehat{a}=0.9999$ and $2\le \ell \le 15$. In the two lower panels, the mode magnitudes are ordered inversely with $\ell$. In other words, the $\ell =2$ mode is the dominant mode, the $\ell =3$ is the next most dominant mode, then $\ell =4$, $\ell =5$, and so on, with $\ell =15$ being the mode with the smallest magnitude (and darkest solid line).

Figure 7

Comparisons between waveforms produced by the FEW and bhpwave models for an EMRI with intrinsic parameters $(M,\mu ,a,p_0,e_0,x)=({10}^6{M}_{\odot },10M_{\odot },0,12.05,0,1)$. The top panels show the evolution of the strain $h_{+}$ with time $t$. In the top plots, the solid (blue) lines correspond to the waveform computed by FEW, while the dashed (orange) lines are produced by bhpwave. The three different panels focus on three small time windows within the full four year signal. The bottom panel plots the difference in the gravitational wave phases between the waveform models. The solid (cyan) line shows the phase difference when the waveforms are evaluated at the same set of parameters. The dashed (black) line is the phase difference after rescaling the FEW parameters to account for different definitions of $G{M}_{\odot }$.

Figure 8

Comparison between the discrete Fourier transform of the time-domain gravitational wave strain $h_{+}$ produced by bhpwave (solid blue line) and the frequency-domain waveform for $h_{+}$ produced by bhpwave using the stationary phase approximation (dashed orange line). For reference, the sensitivity curve $S_n$ is also plotted as the dot-dashed (black) curve. Note that we downsample the number of the frequency samples included in this plot by a factor of 10000 in order to more clearly see how well the two waveforms agree, even in the highly oscillatory region between ${10}^{-3}–{10}^{-2}\mathrm{Hz}$.

Figure 9

Comparisons of the statistical uncertainties $\mathrm{\Delta }{\theta }_{\text{stat}}^i$ (solid lines) and systematic biases $\mathrm{\Delta }{\theta }_{\text{bias}}^i$ (dashed lines) in the intrinsic parameters $(M,\mu ,\widehat{a},p_0)$ as a function of spin $\widehat{a}$ for a range of binaries with different masses. Different mass combinations $(M,\mu )$ corresponds to different colors and markers, as given by the legend in the third panel. The systematic biases are introduced by downsampling the trajectory data to introduce larger interpolation errors into our waveform model.

Figure 10

The impact of flux errors on parameter estimation as a function of spin $\widehat{a}$, as discussed in Sec. 4. The circular markers (purple line) plot the mismatch between $h_E(\theta )$ and $h_B(\theta -{\theta }_{\text{bias}})$, where $h_B$ is our accurate model, $h_E$ is our biased model, $\theta$ are the parameters of the GW source predicted by $h_B$ and ${\theta }_{\text{bias}}$ are the biases in those parameters from the true values $\theta -{\theta }_{\text{bias}}$. The diamond markers (blue line) plot the maximum ratio between these biases and the intrinsic uncertainty in the extrinsic parameters, while the triangle markers (green line) plot the same ratio for the intrinsic parameters. The square markers (yellow line) report the overall dephasing between the two models. The dashed line marks a magnitude of 1.

Figure 11

The fractional interpolation error ${\mathrm{\Delta }}_{(7,6)}^{\mathcal{F}}$ of the flux spline ${\mathcal{F}}_{N(7,6)}^I$ as a function of the grid parameters $\alpha$ and $\chi$.

Figure 12

Convergence of the maximum fractional interpolation error $\Vert {\mathrm{\Delta }}_{(n+1,n)}^{\mathcal{F}}{\Vert }_{\infty }$ as we increase $n$. Increasing $n$ by one amounts to doubling the sampling density in each dimension of our grid. Thus a power law decay of $2^{-4n}$ (dashed blue line) indicates fourth-order convergence of our bicubic spline.

Figure 13

Fractional error between our interpolated flux data ${\mathcal{F}}_E^I$ and the flux data reported within the Circular Orbit Self-force Data repository in the BHPToolkit. The left panel displays errors across the entire domain covered by our spline. The middle plot focuses on a region where we find strong disagreement between our model and the repository. The right plot compares our data to flux calculations from a high-precision ($>100$ digits) Mathematica code in this problem region, demonstrating that our results are consistent with high-precision calculations.

Figure 14

The same comparison as the left plot of Fig. 13, only this time we construct our flux interpolant by imposing not-a-knot boundary conditions (top panel) or natural boundary conditions (bottom panel), leading to worse agreement with the pretabulated flux data.

Figure 15

Convergence of the maximum absolute interpolation error for the frequency-weighted time and phase splines, $\Vert {\mathrm{\Delta }}_{(n,n)}^{\mathrm{\Omega }t}{\Vert }_{\infty }$ (orange squares) and $\Vert {\mathrm{\Delta }}_{(n,n)}^{\mathrm{\Phi }}{\Vert }_{\infty }$ (green diamonds), along with the fractional interpolation error for the unweighted time-spline $\Vert {\mathrm{\Delta }}_{(n,n)}^t{\Vert }_{\infty }$ (blue diamonds). The dashed lines correspond to a power-law convergence of $2^{4n}$.

Figure 16

Convergence of the maximum fractional interpolation error for the frequency spline, $\Vert {\mathrm{\Delta }}_{(n,n)}^{\mathrm{\Omega }}{\Vert }_{\infty }$.

Figure 17

Posterior distribution based on an parallel-tempered MCMC sampling (dark blue lines) and a Fisher matrix analysis (lighter orange lines).