Initial data for black hole–neutron star binaries: A flexible, high-accuracy spectral method

We present a new numerical scheme to solve the initial value problem for black hole–neutron star binaries. This method takes advantage of the flexibility and fast convergence of a multidomain spectral representation of the initial data to construct high-accuracy solutions at a relatively low computational cost. We provide convergence tests of the method for both isolated neutron stars and irrotational binaries. In the second case, we show that we can resolve the small inconsistencies that are part of the quasiequilibrium formulation, and that these inconsistencies are significantly smaller than observed in previous works. The possibility of generating a wide variety of initial data is also demonstrated through two new configurations inspired by results from binary black holes. First, we show that choosing a modified Kerr-Schild conformal metric instead of a flat conformal metric allows for the construction of quasiequilibrium binaries with a spinning black hole. Second, we construct binaries in low-eccentricity orbits, which are a better approximation to astrophysical binaries than quasiequilibrium systems.

Figure 1

Subdomain decomposition close to the compact objects, in the equatorial plane. The apparent horizon of the BH (right) is an inner boundary of the numerical domain, while the surface of the NS (dashed line) is the boundary between the two spherical shells on the left.

Figure 2

Error in the energy density for an isolated NS as a function of the isotropic radius $R_{\mathrm{iso}}$. The reference configuration is obtained by numerical integration of the TOV equations. The spikes in the error function are due to a change in the sign of (${\rho }_{\mathrm{num}}-{\rho }_{\mathrm{an}}$).

Figure 3

Convergence of the angular velocity, the Euler constant (which controls the mass of the star), and the mass of the BH while iterating at the lowest resolution R0 for an equal-mass binary with initial separation $d/M_0=11.507$. The values plotted are the differences from the final results at the highest resolution R3. One step is defined as a passage from point 1 to point 10 in the iterative procedure described in Sec. 3c.

Figure 4

Convergence of the spin of the BH $J_{\mathrm{BH}}$ and the total linear momentum $P_{\mathrm{ADM}}$ at our lowest resolution R0 for an equal-mass binary with initial separation $d/M_0=11.507$.

Figure 5

Convergence of the Hamiltonian and momentum constraints with resolution for an equal-mass binary at a separation $d/M_0=11.507$.

Figure 6

Convergence of the surface fitting method measured as the evolution of the error in the coefficients of the expansion of $R_{\mathrm{surf}}$ in spherical harmonics, computed here as the difference with our results at our highest resolution.

Figure 7

Convergence of the ADM energy, ADM angular momentum, difference between Komar and ADM energies, and position of the center of the BH with resolution for an equal-mass binary at a separation $d/M_0=11.507$. We plot $\Delta E_b=(E_b-E_b^{*})/M_0$, $\Delta J=(J-J^{*})/M_0^2$, $\Delta (E-K)=(\delta M-\delta M^{*})/M_0$, and $\Delta c=||{\mathbf{c}}_{\mathrm{BH}}-{\mathbf{c}}_{\mathrm{BH}}^{*})|{|}_2$, where the reference results $E_b^{*}$, $J^{*}$, $\delta M^{*}$, and ${\mathbf{c}}_{\mathrm{BH}}^{*}$ are those at resolution R3 ($N^{1/3}=51.5$). The difference $\delta M$ between $M_K$ and $E_{\mathrm{ADM}}$, an indication of how close to equilibrium the system is, reaches $2.6\times {10}^{-6}$ at the highest resolution. This is significantly larger than the estimated error in either $E_{\mathrm{ADM}}$, or $\delta M$ shown in the figure.

Figure 8

Binding energies of equal-mass binaries for initial data from our solver, from Taniguchi et al. [21], and from 3PN predictions for model 3PN-M0 (see Table ).

Figure 9

Difference between our results for the binding energy of a sequence of quasiequilibrium equal-masses binaries, and the 3PN predictions from model 3PN-B. The errors represented here come from the difference between the ADM and the Komar energies, except for the 3 closest binaries, for which the numerical error can no longer be neglected. We also represent the influence of tidal effects (from model 3PN-M0) and eccentricity (model 3PN-MA). Any binary with an eccentricity $e=0.01$, initially closer to its apocenter than to its pericenter, should have a binding energy between the results from models 3PN-M0 and 3PN-MA. Model 3PN-MP, representing an eccentric binary at its pericenter, is not plotted here, but predicts even higher energies than model 3PN-M0.

Figure 10

Angular momentum of equal-mass binaries for initial data generated by our solver, initial data from Taniguchi et al. [21], and 3PN predictions. We see that both numerical results are in very good agreement, even when they begin to diverge from the 3PN models of circular orbits.

Figure 11

Evolution of equal-mass binaries after 0, 1, and 2 steps of our iterative method reducing the eccentricity. We plot the time derivative of the coordinate separation between the BH and the NS, $2v_r=\dot{d}$.