A Details of BHBH test simulations

Relativistic Simulations of Black Hole-Neutron Star Mergers: Effects of black-hole spin


Black hole-neutron star (BHNS) binary mergers are candidate engines for generating both short-hard gamma-ray bursts (SGRBs) and detectable gravitational waves. Using our most recent conformal thin-sandwich BHNS initial data and our fully general relativistic hydrodynamics code, which is now AMR-capable, we are able to efficiently and accurately simulate these binaries from large separations through inspiral, merger, and ringdown. We evolve the metric using the BSSN formulation with the standard moving puncture gauge conditions and handle the hydrodynamics with a high-resolution shock-capturing scheme. We explore the effects of BH spin (aligned and anti-aligned with the orbital angular momentum) by evolving three sets of initial data with BH:NS mass ratio : the data sets are nearly identical, except the BH spin is varied between (anti-aligned), 0.0, and 0.75. The number of orbits before merger increases with , as expected. We also study the nonspinning BH case in more detail, varying between 1, 3, and 5. We calculate gravitational waveforms for the cases we simulate and compare them to binary black-hole waveforms. Only a small disk () forms for the anti-aligned spin case () and for the most extreme mass ratio case (). By contrast, a massive (), hot disk forms in the rapidly spinning () aligned BH case. Such a disk could drive a SGRB, possibly by, e.g., producing a copious flux of neutrino-antineutino pairs.


I Introduction

Mergers of black hole-neutron star (BHNS) binaries are expected to be among the most promising sources of gravitational waves detectable by ground-based laser interferometers like LIGO Abbott and the LIGO Scientific Collaboration (2008); Brown et al. (2004), VIRGO Acernese and the VIRGO Collaboration (2006); Beauville and the LIGO-VIRGO Working Group (2008), GEO Lück and the GEO600 collaboration (2006), TAMA Ando and the TAMA collaboration (2002); Tatsumi and the TAMA collaboration (2007), and AIGO aig (), as well as by the proposed space-based interferometers LISA Heinzel et al. (2006) and DECIGO Kawamura and the DECIGO collaboration (2006). Theoretical models indicate that a neutron star-neutron star (NSNS) Duez et al. (); Stephens et al. (); Shibata et al. (2006, 2003, 2005); Shibata and Taniguchi (2006); Shibata et al. (2007) or black hole-neutron star (BHNS) Faber et al. (2006a, b); Shibata and Uryū (2006, 2007); Shibata et al. (2007); Shibata and Taniguchi (2008) binary merger may result in a hot, massive disk around a black hole, whose temperatures and densities could be sufficient to trigger a short-hard gamma-ray burst (SGRB). Indeed, SGRBs have been repeatedly associated with galaxies with extremely low star formation rates (see Gehrels et al. (2007) and references therein for a review), indicating that the source is likely to involve an evolved population, rather than main sequence stars. The number of detectable BHNS mergers in the observable universe is still an open question. The uncertainties arise from some aspects of population synthesis calculations that are only partially understood. The estimated event rate of BHNS mergers observable by an Advanced LIGO detector typically falls in the range  Kalogera et al. (2007).

Motivated by the significance of BHNS binaries for both GRB and gravitational-wave physics, many simulations of BHNS systems have been performed over the past years. In most dynamical simulations of BHNS binaries to date, the self-gravity of the neutron star (NS) and/or the tidal gravity of the black hole (BH) are treated in a Newtonian or post-Newtonian framework (see, e.g., Lee (2001); Rosswog et al. (2004); Rosswog (2005); Kobayashi et al. (2004); Rantsiou et al. (2008); see also Löffler et al. (2006), who provide fully relativistic head-on collision simulations).

Our first BHNS merger simulations were performed in an approximate relativistic framework Faber et al. (2006a, b), which assumed conformal flatness of the spatial metric throughout the evolution (see Isenberg (2008); Wilson et al. (1996)). In particular, we evolved extreme mass ratio (1) initial data developed earlier by our group Baumgarte et al. (2004); Taniguchi et al. (2005). Though this simulation technique only allows for crude estimates, we found that mergers of irrotational BHNS binaries may lead to disks with masses up to 0.3 , with sufficient heating to emit the neutrino fluxes required to launch an SGRB.

To date the only self-consistent, general relativistic, dynamical simulations of BHNS inspiral and coalescence have been those of Shibata et al. Shibata and Uryū (2006, 2007); Shibata and Taniguchi (2008); Yamamoto et al. (2008), Etienne et al. Etienne et al. (2008) (hereafter Paper I), and Duez et al. Duez et al. (2008). In many earlier calculations, especially those with initial mass ratios , significant disks are formed after the NS is disrupted, and, for very stiff nuclear equations of state (EOSs), the core of the NS may survive the initial mass transfer episode and remain bound. These findings contrast with some semi-analytic relativistic arguments, which suggest that disks with appreciable masses are difficult to form via the merger of BHNS binaries Miller (2005).

In their first set of fully relativistic BHNS simulations, Shibata and Uryū found disk masses in the range of 0.1 – 0.3 for corotating neutron stars Shibata and Uryū (2006, 2007). Subsequently, Shibata and Taniguchi (Shibata and Taniguchi (2008), hereafter ST) found smaller disk masses of 0.04 – 0.16 for irrotational neutron stars, which are more realistic. However, our own fully general relativistic simulations of irrotational BHNSs in Paper I suggest that the disk mass is no more than 0.04 . Such a small disk mass, if true for generic cases, might disfavor BHNS mergers as possible central engines for SGRBs.

Recently, Yamamoto, Shibata, and Taniguchi presented a new code which implements an adaptive mesh refinement (AMR) algorithm Yamamoto et al. (2008). In one of their test simulations of BHNS binaries, they use the same initial data as in ST, but they find that the disk mass is less than 0.001 . This result disagrees with that in ST, but is consistent with our result in Paper I. They attribute the discrepancy to different grid structures and computational methods.

Most recently, Duez et al. developed a new fully general relativistic hydrodynamics code Duez et al. (2008) that evolves the generalized harmonic formulation of Einstein’s equations using a pseudospectral method, and the equations of hydrodynamics via shock-capturing, finite-difference techniques. They use this code to perform a simulation of a single equal-mass BHNS and find that the disk mass is less than 0.03, which agrees with our result in Paper I.

As we report in this paper, our code has been modified to use the moving-box Carpet Schnetter et al. (2004) AMR infrastructure. With AMR, our code is now capable of performing BHNS simulations with resolutions equivalent to Paper I at only about the total computational cost. This enables us to simulate strong-field regions at higher resolution while simultaneously extending our grid’s outer boundaries much farther into the wavezone for more accurate wave extraction. We also evolve initial data at larger binary separations than in Paper I. Some of the BHNS initial data evolved in this paper are also new, as the BHs are now spinning. All of our initial data are constructed in the conformal thin-sandwich (CTS) formalism. The BH spin is incorporated by imposing boundary conditions on the BH horizon as in Cook and Pfeiffer (2004); Caudill et al. (2006). The initial NS is irrotational and is modeled by an polytropic EOS.

In this paper, we present a new set of BHNS simulations that probes how the initial BH spin and the binary mass ratio affect the dynamics and outcome of the merger. In particular, we study mergers of nonspinning BHNS binaries with mass ratios , , and . For the mass ratio , we study three cases in which the BH spin parameter () is (counter-rotating), 0 and 0.75. As in Paper I, we focus on binaries in which the NS tidally disrupts before reaching the innermost stable circular orbit (ISCO) and plunging into the BH. These systems are the most likely to yield a significant disk after merger. We thus do not evolve extreme mass-ratio binaries in this paper.

Starting with approximately the same initial orbital angular velocity , we find that when the BH is spinning, the merger time is delayed when the BH spin aligns with the orbital angular momentum, and is shortened when the BH spin is anti-aligned with the orbital angular momentum. For example, evolutions of initial data with [ is the initial Arnowitt-Deser-Misner (ADM) mass of the system], require about 4.25 orbits before merger for , but 6.5 orbits before merger for . This finding is expected, since for large, aligned spins there is more angular momentum that must be radiated away to bring the binary to merger. This result is consistent with similar findings in the case of BHBHs (see e.g. Campanelli et al. (2006a)).

After the merger, we find a remnant disk with rest mass between and (assuming the NS has a rest mass of ). When the initial BH is nonspinning, the disk mass is for , for and tiny () for . For a fixed mass ratio , the disk mass is for and tiny () for . The increase in disk mass with increasing BH spin agrees qualitatively with the semi-relativistic, smoothed particle hydrodynamics (SPH) simulations of large mass-ratio () BHNS mergers reported in Rantsiou et al. (2008).

We find that the remnant disks are hot. A rough estimate suggests that the temperature in the disks is K. The results of BH disk simulations in Setiawan et al. (2006) suggest that the disk could generate a gamma-ray energy erg from neutrino-antineutrino annihilation, which is promising for BHNS mergers as plausible central engines for SGRBs.

We compute gravitational waveforms and the effective wave strain in the frequency domain for our mergers. As in Paper I, we find measurable differences between BHNS waveforms and those produced by BHBH mergers within the Advanced LIGO band. These differences appear at frequencies at which the NS is tidally disrupted and accreted by the black hole. Extracting the masses and spins of the objects from the inspiral data, the merger waveforms could provide information about the radius of the NS, which in turn may be used to constrain the NS equation of state.

In contrast with Paper I, we find substantial disks in two of our nonspinning BH cases. The reason is that in Paper I, we imposed limits on pressure to prevent spurious heating/cooling in the low-density artificial “atmosphere”. We find that imposing these artificial limits spuriously removes angular momentum, especially during the NS tidal disruption phase, thereby suppressing disk formation (see Section V.4 for further details).

As we cautioned in Paper I, our findings are fundamentally limited by the physics we do not model in these preliminary simulations (e.g., neutrino transport, magnetic fields, etc), as well as by our choice of a simple -law EOS to model the (cold and hot) nuclear EOS. Disk masses may depend sensitively on the initial structure of the neutron star, especially the low-density outer regions, which is determined by the adopted EOS. So all BHNS merger results should be viewed in light of these caveats. Our main focus has been to establish that reliable simulations in full general relativity can now be launched to address these important issues.

This paper is organized as follows. In Secs. II and III, we briefly outline the basic equations and their specific implementation in our general relativistic, hydrodynamics scheme. Here we also provide an overview of our initial data, gauge conditions, matter evolution equations, and diagnostics. Sec. IV presents code tests designed to validate our new AMR-based scheme. In Sec. V, we review the results of our BHNS merger simulations. In Sec. VI we summarize our findings and comment on future directions.

Ii Basic Equations

The formulation and numerical scheme for our simulations are basically the same as those already reported in Duez et al. (2005a); Etienne et al. (2008), to which the reader may refer for details. Here we introduce our notation, summarize our method, and point out the latest changes to our numerical technique. Geometrized units () are adopted, except where stated explicitly. Greek indices denote all four spacetime dimensions (0, 1, 2, and 3), and Latin indices imply spatial parts only (1, 2, and 3).

We use the 3+1 formulation of general relativity and decompose the metric into the following form:


The fundamental variables for the metric evolution are the spatial three-metric and extrinsic curvature . We adopt the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism Shibata and Nakamura (1995); Baumgarte and Shapiro (1998) in which the evolution variables are the conformal exponent , the conformal 3-metric , three auxiliary functions , the trace of the extrinsic curvature , and the tracefree part of the conformal extrinsic curvature . Here, . The full spacetime metric is related to the three-metric by , where the future-directed, timelike unit vector normal to the time slice can be written in terms of the lapse and shift as . The evolution equations of these BSSN variables are given by Eqs. (9)–(13) in Paper I.

It has been suggested that evolving or instead of gives more accurate results in BHBH simulations (see e.g. Brügmann et al. (2008a); Marronetti et al. (2008); Lousto and Zlochower (2008)). We have tried all three techniques. We find that in the presence of hydrodynamic matter, the evolution of the variable yields a more stable evolution inside the apparent horizon (AH) near the puncture, leading to better rest-mass conservation (see Sec. III.2). We therefore adopt -evolution in all of our BHNS simulations.

It has been suggested that Kreiss-Oliger dissipation is sometimes useful to add in the BSSN evolution equations outside the AH to reduce high-frequency numerical noise associated with AMR refinement interfaces Baker et al. (2006). We have also tried this technique, and found that the dissipation results in slightly smaller Hamiltonian and momentum constraint violations. We typically do not use it in our BHNS simulations. However, in one case (Case E, see Sec. V), we found that the presence of hydrodynamic matter causes the conformal related metric to lose positive definiteness near the puncture, which is unphysical. This behavior eventually causes the code to crash. We find that adding Kreiss-Oliger dissipation, coupled with a high but finite pressure ceiling deep inside the AH, solves this problem.

We adopt standard puncture gauge conditions: an advective “1+log” slicing condition for the lapse and a “Gamma-freezing” condition for the shift van Meter et al. (2006). Thus we have


where . The parameter is set to at the beginning of all simulations, where is the initial ADM mass of the BHNS binary. If the value for is fixed at this value throughout the simulation, both Hamiltonian and momentum constraint violation increase to marginally unacceptable levels (, as measured by Eqs. (40) and (41) in Paper I) outside the BH during the merger phase, especially in the high BH spin cases. By simply doubling to early in the inspiral phase (specifically, at ) and fixing it at that value significantly reduces constraint violations outside the BH during merger. However, late in the merger stage of Case B, we find that Hamiltonian constraint violation increases rapidly () outside the BH, unless we enable Kreiss-Oliger dissipation as well.

The fundamental matter variables are the rest-mass density , specific internal energy , pressure , and four-velocity . We adopt a -law EOS with , which reduces to an polytropic law for the initial (cold) neutron star matter. The stress-energy tensor is given by


where is the specific enthalpy. In our numerical implementation of the hydrodynamics equations, we evolve the “conservative” variables , , and . They are defined as


The evolution equations for these variables are given by Eqs. (21)–(24) in Paper I.

Iii Numerical Methods

iii.1 Initial data

Our initial data are constructed by solving Einstein’s constraint equations in the CTS formalism, which allows us to impose an approximate helical Killing vector by setting the time derivatives of the conformally related metric to zero. We model the NS as an irrotational polytrope, and impose the black hole equilibrium boundary conditions of Cook and Pfeiffer Cook and Pfeiffer (2004) on the black hole horizon. Details of this method can be found in Taniguchi et al. (2008). Our initial data here differ from those described in Taniguchi et al. (2008) only in terms of the black hole spin, which can be specified by virtue of a free (vectorial) parameter that appears in the boundary conditions of Cook and Pfeiffer (2004). In Taniguchi et al. (2008) we adopted the method in Caudill et al. (2006) by iterating over this parameter until the quasilocal spin of the black hole vanishes, making both binary components irrotational. Here we also consider rotating black hole configurations, for which we iterate over until the black hole spin equals certain specified values. We focus on black hole spins that are aligned or anti-aligned with the orbital angular momentum, and consider (aligned), (nonspinning), and (anti-aligned).

The initial data are calculated using the Lorene spectral methods numerical libraries web (). To map these spectral configurations onto our non-spectral simulation grid, we first construct our numerical grid and record the positions of each point in physical coordinates. Then we evaluate the field and hydrodynamic quantities at these points based on their spectral coefficients ft0 (). Finally, the excised BH region is filled with constraint-violating initial data, using the “smooth junk” technique we developed and validated in Etienne et al. (2007) (see also Brown et al. (2007, 2008)). In particular, we extrapolate all initial data quantities from the BH exterior into the interior with a 7th order polynomial, using a uniform stencil spacing of .

All of the NSs considered in this paper have a compaction of , where is the ADM mass and is the (circumferential) radius of the NS in isolation. Since we model the NS with an () polytropic EOS, the rest mass of the NS, , scales with the polytropic constant as . For a NS with compaction , we find the ADM mass for the isolated NS to be , with an isotropic radius and circumferential (Schwarzschild) radius of . The maximum rest-mass density of this NS is .

iii.2 Evolution of the metric and hydrodynamics

We evolve the BSSN equations with fourth-order accurate, centered finite-differencing stencils, except on shift advection terms, where we use fourth-order accurate upwind stencils. We apply Sommerfeld outgoing wave boundary conditions to all BSSN fields, as in Paper I. Our code is embedded in the Cactus parallelization framework Cac (), and our fourth-order Runge-Kutta timestepping is managed by the MoL (Method of Lines) thorn, with a Courant-Friedrichs-Lewy (CFL) factor set to 0.45 in all BHNS simulations. We use the Carpet Schnetter et al. (2004) infrastructure to implement the moving-box adaptive mesh refinement. In all AMR simulations presented here, we use second-order temporal prolongation, coupled with fifth-order spatial prolongation. The apparent horizon (AH) of the black hole is computed with the AHFinderDirect Cactus thorn Thornburg (2004).

We write the general relativistic hydrodynamics equations in conservative form. They are evolved by a high-resolution shock-capturing (HRSC) technique Duez et al. (2005a) that employs the monotonized central (MC) reconstruction scheme van Leer (1977) coupled to the Harten, Lax, and van Leer (HLL) approximate Riemann solver Harten et al. (1983). The adopted hydrodynamic scheme is second-order accurate for smooth flows, and first-order accurate when discontinuities (e.g. shocks) arise. To stabilize our hydrodynamic scheme in regions where there is no matter, we maintain a tenuous atmosphere on our grid, with a density floor set equal to times the initial maximum density on our grid. The initial atmospheric pressure is set equal to the cold polytropic value , where is the polytropic constant at . Throughout the evolution, we impose limits on the atmospheric pressure to prevent spurious heating and negative values of the internal energy . Specifically, we require , where and . Whenever exceeds or drops below , we reset to or , respectively. Applying these limits everywhere on our grid would artificially sap the angular momentum in the tidally disrupted NS, allowing matter to fall spuriously into the BH and thereby suppressing disk formation (see Sec. V). To effectively eliminate this spurious angular momentum loss, we impose these pressure limits only in regions where the rest-mass density remains very low () or deep inside the AH, where (see Sec. V.4 for more details).

At each timestep, we need to recover the “primitive variables” , , and from the “conservative” variables , , and . We perform the inversion as specified in Eqs. (57)–(62) of Duez et al. (2005a), but with a slightly modified analytic quartic solver from the GNU Scientific Library that outputs only the real roots. We use the same technique as in Paper I to ensure that the values of and yield physically valid primitive variables, except we reset to (where is the maximum value of initially) when either or is unphysical [i.e., violate one of the inequalities (34) or (35) in Paper I]. The restrictions usually apply only to the region near the puncture or in the low-density atmosphere.

iii.3 Diagnostics

During the evolution, we monitor the Hamiltonian and momentum constraints calculated by Eqs. (40)–(43) of Paper I. Gravitational waves are extracted using both the Regge-Wheeler-Zerilli-Moncrief formalism and the Newman-Penrose Weyl scalar . The extraction technique is summarized in Sec. IIIF of Paper I.

We also monitor the ADM mass and angular momentum of the system, which can be calculated by surface integrals at a large distance (Eqs. (37) and (39) of Paper I). However, with our AMR grid, the resolution is rather low in regions very far from the binary, which causes our angular momentum measurement to suffer from numerical noise. This problem can be overcome by using Gauss’s law to split the distant surface integral into a sum of integrals over an inner surface and the volume between the inner and distant surfaces. The reasons are twofold: (1) the inner surface integrals are computed more accurately, as the inner surface resides in the region where numerical resolution is higher, and (2) most of the contribution in the volume integral is from the strong-field domain, which is also in the inner (high-resolution) region. The equations we use to calculate ADM mass and angular momentum with minimal numerical noise are as follows (see e.g. Appendix A in Yo et al. (2002) for a derivation):


Here is a surface surrounding the BH, is the volume between the inner surface and a distant surface, , , , is the Ricci scalar associated with , and are Christoffel symbols associated with .

When hydrodynamic matter is evolved on a fixed uniform grid, our hydrodynamic scheme guarantees that the rest mass is conserved to machine roundoff error. This strict conservation is no longer maintained in an AMR grid, where spatial and temporal prolongation is performed at the refinement boundaries. Hence we also monitor the rest mass


during the evolution. Rest-mass conservation is also violated whenever is reset to the atmosphere value. This usually happens only in the very low-density atmosphere or deep inside the AH where accuracy is not maintained. The low-density regions do not affect rest-mass conservation significantly. However, this is not the case near the puncture, where the rest-mass density can be very high, especially after the NS is tidally disrupted and matter falls into the BH and accumulates near the puncture. The evolution near the puncture is sensitive to the choice of conformal variable used in the BSSN evolution, as discussed in Sec. II. We find that the evolution leads to better rest-mass conservation than or evolutions.

We measure the thermal energy generated by shocks via the quantity , where is the pressure associated with the cold EOS that characterizes the initial matter. The specific internal energy can be decomposed into a cold part and a thermal part: with


Hence the relationship between and is


For shock-heated gas, we always have (see Appendix B).

Iv AMR Code Tests

Our general relativistic magnetohydrodynamic (GRMHD) code has been thoroughly tested by passing a robust suite of tests. These tests include maintaining stable rotating stars in stationary equilibrium, reproducing the exact Oppenheimer-Snyder solution for collapse to a black hole, and reproducing analytic solutions for MHD shocks, nonlinear MHD wave propagation, magnetized Bondi accretion, and MHD waves induced by linear gravitational waves Duez et al. (2005a, b). Our code has also been compared with Shibata & Sekiguchi’s GRMHD code Shibata and Sekiguchi (2005) by performing simulations of identical magnetized hypermassive neutron stars Duez et al. (2006a, b) and of the magnetorotational collapse of identical stellar cores Shibata et al. (2006). We obtain good agreement between these two independent codes. Our code has also been used to simulate the collapse of very massive, magnetized, rotating stars to black holes Liu et al. (2007); merging BHBH binaries Etienne et al. (2007), BHNS binaries (Paper I), magnetized NSNS binaries Liu et al. (2008); and relativistic hydrodynamic matter in the presence of puncture black holes Faber et al. (2007). Recently, our code has been generalized to incorporate (optically thick) radiation transport and its feedback on MHD fluids Farris et al. (2008).

All of the above-mentioned tests and simulations were performed on grids with uniform spacing. In some of the simulations, we utilized the multiple-transition fisheye transformation Campanelli et al. (2006b) so that a uniform computational grid spacing corresponds to physical coordinates with spatially varying resolution. Recently, we have modified our code so that we can use the moving-box AMR infrastructure provided by Carpet Schnetter et al. (2004). To test our new code, we have performed shock-tube tests and 3+1 simulations of linear gravitational waves, single stationary and boosted puncture BHs, puncture BHBH binaries, stable, rapidly and differentially rotating relativistic stars, and relativistic Bondi accretion onto a Schwarzschild BH. All of our 3+1 AMR code tests were performed assuming equatorial symmetry (i.e., symmetry about the orbital plane), which we assume in all BHNS evolutions presented in this paper. In many of the tests we have analytic solutions with which to compare. For the BHBH tests, we compare our results with the numerical results reported in the literature. We have confirmed that our new code can evolve these systems reliably with AMR. Below we briefly report results from both our BHBH and differentially rotating relativistic star AMR code tests.

iv.1 Binary black holes

In this subsection we summarize an important testbed for calibrating the vacuum sector of our code: BHBH merger simulations. We provide details on the setup of our initial data and numerical grid, as well as on the adopted numerical techniques and simulation results, in the tables and figures of Appendix A. Here we point out some highlights.

We first simulate an equal-mass BHBH system. We use puncture initial data Beig and Ó Murchadha (1994, 1996); Brandt and Brügmann (1997) with initial binary separation , where is the ADM mass of the binary. The initial “bare” mass and momentum of each puncture are chosen according to Tichy and Brügmann (2004) to make the orbit quasicircular (see also Baumgarte (2000)). This initial configuration is close to the R4 run presented in Baker et al. (2006, 2007). We use 8 refinement levels centered at each BH and place the outer boundary at . The resolution in the innermost refinement level is , so that there are about 15 grid points across the diameter of each black hole’s AH initially. This number rapidly increases during the early inspiral to a fixed value of about 24 grid points due to our gauge choice, and increases again to about 40 grid points, post-merger. We find that the binary inspiral orbit, gravitational waveform, energy and angular momentum emitted by gravitational waves, and the final black hole’s spin parameter are very similar to those reported by others in Baker et al. (2006, 2007).

To further validate our numerical results, we compute the quantities


where is the ADM angular momentum of the initial binary, and are the ADM mass and angular momentum of the final (merged) system, and and are the energy and angular momentum carried off by gravitational waves. We find that both and are very close to and , the mass and angular momentum of the final BH. We compute and from the irreducible mass and the ratio of proper polar to equatorial circumference of the final black hole’s horizon. Specifically, we first use Eq. (5.3) of Alcubierre et al. (2005) to solve for the dimensionless spin :


where is the complete elliptic integral of the second kind. We next compute by and by . These formulae are derived for a Kerr spacetime. They are applicable to the merged BH as the spacetime approaches Kerr at late times. Gravitational waves are extracted via the Newman-Penrose scalar . The quantities and are computed by integrating Eqs. (50) and (52) of Paper I. Conservation of energy and angular momentum demands that . Hence these parameters are good indicators of the accuracy of the simulations, whenever a significant amount of energy and angular momentum are radiated during the evolution. We found and in our equal-mass BHBH simulation.

We next perform an unequal-mass BHBH simulation. The initial data are generated by the TwoPunctures code Ansorg et al. (2004) with initial binary (coordinate) separation . The approximate initial momentum of each BH for a quasicircular configuration is calculated by the 3PN expression using Eq. (65) of Brügmann et al. (2008a). The “bare” masses of the punctures are adjusted so that the irreducible mass ratio of the BHs equals 3. The initial configuration approximates one of the cases studied in Berti et al. (2007). We use 9 and 10 refinement levels centered at the larger and smaller BH, respectively, with an outer boundary at . The resolution around the large and small BH is and , respectively, so that there are initially about 36 and 22 grid points across the diameter of the larger and smaller black holes’ AH, respectively. We find that our computed values of and agree with Table I of Berti et al. (2007) to within 4%, and that the kick velocity agrees with Fig. 2 of González et al. (2007). We note that exact agreement is not expected as our initial configurations are slightly different. Finally, we find excellent energy and angular momentum conservation: and .

These tests demonstrate our code’s ability to evolve dynamical vacuum spacetimes containing moving BHs. We next study our code’s ability to handle hydrodynamics in a strong-field spacetime with moving AMR boxes.

iv.2 Rotating relativistic star

Figure 1: The moving AMR grid in the rotating relativistic star test. Shown here are the snapshots of the two innermost refinement boxes (the nested squares) and part of the third refinement level in the equatorial plane. The circle denotes the boundary of the star.

In our BHNS simulations, the NS is placed inside a high resolution refinement box initially. However, when it is tidally disrupted by the black hole, the NS matter spreads into other refinement boxes, with different resolutions. The purpose of this test is to demonstrate the ability of our AMR code to maintain the equilibrium of a stable, differentially rotating relativistic star crossing refinement boundaries during the evolution, and verify second-order convergence as the resolution is improved.

The star considered in this test is the same hypermassive neutron star studied in Baumgarte et al. (2000), which is modeled by an polytropic EOS. The mass of the star is 1.7 times greater than the maximum mass of a nonrotating neutron star with the same EOS, but has been shown to be dynamically stable in previous simulations Baumgarte et al. (2000); Duez et al. (2003, 2005a). The star is highly flattened because of its rapid rotation. It has an equatorial (coordinate) radius , and a polar radius . The rapid rotation also causes a torus-like density distribution in which the maximum density occurs off-center. The central rotation period of the star is .

Our AMR grid contains four refinement levels that move in time. The side-lengths of the refinement boxes in the equatorial plane are , , , and . Since we impose equatorial symmetry, the lengths in the -direction are half of the sizes stated above. The coarsest grid is fixed and the outer boundary is placed at , and . We initially center the star at the origin, and move the center of the four refinement boxes according to the following prescription:


where , and are constants. We set the amplitudes , and the angular frequency of the moving center (). With these parameters, the bulk of the star crosses the innermost refinement boxes during the evolution (see Fig. 1). The grid spacing doubles at each successive refinement level. We evolve the star with three different resolutions, which we label low, medium and high resolution runs. The grid spacing in the innermost refinement level is , , and for the low, medium and high resolution runs, respectively.

Figure 2: Evolution of rest mass , ADM mass , and angular momentum of a rapidly rotating star using a moving AMR grid with low (black solid lines), medium (red dotted lines) and high (blue dashed lines) resolution.

Figure 2 shows the evolution of the rest mass , ADM mass , and angular momentum of the star, calculated by Eqs. (10)–(13). Since there is no black hole in the spacetime, we calculate only the volume integrals of and over the entire spatial grid. As mentioned in Sec. III.3, strict rest-mass conservation is not expected when using an AMR grid. However, the deviations in , and converge to zero at slightly higher than second order with increasing resolution. In addition, the amplitudes of the gravitational radiation decrease to zero at second order with increasing resolution. Gravitational radiation does not contribute significantly to the nonconservation of and observed in Fig. 2; and for all three runs.

Figure 3: Density profiles in the equatorial plane along the -axis (upper graph) and -axis (lower graph) for the high resolution run. The rest-mass density is normalized to the initial maximum density, . Solid (black), dotted (red) and dashed (blue) lines show the profiles at times , and , respectively.

Figure 3 demonstrates that the equilibrium profile of the rest-mass density is well-maintained in the high resolution simulation, apart from a slight drift of the star’s center.

V Results

v.1 Overview

A 0.00 3.0 8.81 0.702 0.0333
A-MSep 0.00 3.0 7.17 0.668 0.0441
A-SSep 0.00 3.0 5.41 0.638 0.0623
B 0.75 3.0 8.81 1.096 0.0328
B-SSep 0.75 3.0 5.53 1.011 0.0594
C 0.50 3.0 8.81 0.443 0.0338
D 0.00 5.0 8.79 0.518 0.0333
E 0.00 1.0 8.61 0.938 0.0347
Table 1: Initial BHNS configurations. Here, the mass ratio is defined as , where and are the ADM masses of the black hole and neutron star at infinite separation, is the total ADM mass of the binary system, the initial binary coordinate separation, is the spin parameter of the black hole (always aligned or anti-aligned with the orbital angular momentum), and is the initial orbital angular velocity. All neutron stars have the same nondimensional rest (baryon) mass , which is the maximum rest mass of a nonrotating NS with the same polytropic EOS.
Case Grid Hierarchy (in units of ) Max. resolution
A 0.0333 (196.7, 98.35, 49.18, 24.59, 12.29, 6.147, 3.073, 1.414 [1.660]) 41 85
A-MSep 0.0441 (197.0, 98.49, 49.24, 24.62, 12.31, 6.156, 3.078, 1.416 [1.662]) 40 82
A-SSep 0.0623 (197.3, 98.67, 49.33, 24.67, 12.33, 6.167, 3.083, 1.418 [1.665]) 40 76
B 0.0328 (210.2, 92.49, 46.24, 23.12, 11.56, 5.780, 2.890, 1.445 [1.642], 0.7554 [N/A]) 56 80
B-SSep-HRes 0.0594 (197.6, 98.79, 49.40, 24.70, 12.35, 6.174, 3.087, 1.544, 0.7718 [N/A]) 59 78
B-SSep-MRes 0.0594 (197.6, 98.79, 49.40, 24.70, 12.35, 6.174, 3.087, 1.544, 0.7718 [N/A]) 43 57
B-SSep-LRes 0.0594 (197.6, 98.79, 49.40, 24.70, 12.35, 6.174, 3.087, 1.544, 0.7718 [N/A]) 37 50
C 0.0338 (196.6, 98.31, 49.16, 24.58, 12.29, 6.145, 3.072, 1.413 [1.659]) 37 85
D 0.0333 (196.3, 98.12, 49.06, 24.53, 12.26, 4.415, 2.208, 1.104) 69 84
E 0.0347 (217.6, 108.8, 54.40, 27.20, 13.60, 6.801, 3.400, 1.496 [N/A], 0.7349 [N/A]) 47 78

There are two sets of nested refinement boxes: one centered on the NS and one on the BH. This column specifies the half side length of the refinement boxes centered on both the BH and NS. When the side length around the NS is different, we specify the NS half side length in square brackets. If there is no corresponding NS refinement box (as is the case when the NS is significantly larger than the BH), we write [N/A] for that box.

Table 2: Grid configurations used in our simulations. Here, denotes the number of grid points covering the diameter of the (spherical) apparent horizon initially, and denotes the number of grid points covering the smallest diameter of the neutron star initially.
Figure 4: AMR grid setup during a typical BHNS simulation, as observed in the equatorial plane. Only the two innermost refinement boxes are shown for clarity. The black region denotes the apparent horizon of the black hole. Neutron star density contours are drawn according to (=0, 1,… 12). The maximum rest-mass density of the initial NS is , or g cm.
Case (km/s)
A 3.9% 0.559 0.93% 17.4% 33 4.5 2.2%
A-MSep 3.8% 0.560 0.81% 13.2% 33 2.5 1.9%
A-SSep 2.9% 0.557 0.45% 11.4% 33 1.75
B 15% 0.881 0.92% 13.3% 22 6.5 2.8%
B-SSep 14.3% 0.882 0.60% 6.4% 79 2.0 2.6%
C 0.8% 0.331 0.98% 24.4% 49 3.25 1.3%
D 0.8% 0.418 0.98% 19.2% 73 6.25 1.9%
E 2.3% 0.851 0.35% 7.2% 17 2.25 0.9%
Table 3: BHNS simulation results. Here is the rest mass of the material outside the AH at the end of the simulation, is the value of computed by solving Eq. (18) at late times. The total energy and angular momentum carried off by the gravitational radiation are given by and , respectively. is the kick velocity due to recoil. specifies the number of orbits in the inspiral phase before merger. and measure the violation of energy and angular momentum conservation, as defined in Eqs. (16) and (17), respectively.

We perform a number of BHNS simulations with varying initial configurations and numerical resolutions. Table 1 provides an overview of our chosen initial configurations, and Table 2 specifies the AMR grid structure used in each case. The first three cases in Table 1 (Cases A, A-MSep and A-SSep) correspond to the same BHNS binary system tracked from different initial separations along a quasiequilibrium sequence. The objective of these simulations is to demonstrate that the fate of the binaries, and their merger waveforms, are insensitive to the initial binary separation at which we begin the simulation, provided it is sufficiently large. For Case B-SSep, we evolve the system with three different resolutions (see Table 2) to demonstrate convergence. Fixing the mass ratio at and initial orbital angular frequency at , we then study the effects of black hole spin by choosing the spin parameter between and 0.75 (Cases C, A and B). In all cases, the black hole spin is either aligned or anti-aligned with the orbital angular momentum. Finally, we study the dependence on the mass ratio by varying between 1 and 5 (Cases E, A and D).

In all simulations, AMR refinement levels are initially centered on the BH and NS. After , the levels centered on the BH track the AH centroid, and the levels centered on the NS track the matter centroid defined by :


where is the volume outside the innermost refinement box surrounding the BH. Figure 4 shows snapshots of the inner refinement boxes during the Case A simulation.

In all of our BHNS simulations, we find that the integrated Hamiltonian constraint violation [as measured by Eq. (40) in Paper I] is 0.1%–0.6% outside the BH during the inspiral phase. During merger, it jumps to a peak value of 1%–2% in Cases A, C, D, and E and to in Case B. Post-merger, the violation stays roughly constant (Case E) or gradually decreases (all other cases) until the simulation is stopped. The integrated momentum constraint violation [as measured by Eq. (41) in Paper I] outside the BH generally oscillates around 1%–2%, except for a few peaks of 3%–6% during inspiral. We measure energy and momentum conservation by the quantities and defined in Eqs. (16) and (17). We find that and 1–2% (see Table 3) in all cases.

We also monitor the rest-mass conservation by computing the quantity


where is the NS rest mass computed at . Rest-mass conservation requires at all times. Most of the violation occurs inside the BH apparent horizon when the NS matter falls into the BH. We find that 0.06–0.08% during the inspiral phase (for all cases except E, which has ). During and after the merger phase, we find large “glitches” in inside the AH due to inaccuracies in that region, which causes large (3%), but these are easily accounted for and neglected. We find that is 0.1%–0.15% after the inspiral phase. This result is important because in some cases we find a low-mass () disk surrounding the BH. These results can only be trusted if the rest-mass violation outside the apparent horizon due to numerical error is much smaller than the disk mass, which is true in all of our simulations.

Figure 5: The BH’s irreducible mass and versus time, for Case A.

Figure 5 shows the BH’s irreducible mass and as a function of time for Case A, where and are the proper polar and equatorial circumferences, respectively. Notice that these two quantities remain almost constant at early times until the NS matter falls into the BH. They approach constant values again at late times as the system settles into a quasistationary state. We find the same qualitative behavior for all of our BHNS simulations. Given , we define using Eq. (18). When the BH is well-approximated by the Kerr metric (e.g., the early inspiral phase, at large separations), equals the BH’s spin parameter . In its final stage, the system consists of a BH surrounded by a disk. The spacetime near the BH therefore deviates from Kerr, and does not need to coincide with the value of the BH spin parameter defined in, e.g. the isolated horizon formulation Ashtekar and Krishnan (2004). Nevertheless, is still a reasonable measure of the BH’s spin when the disk is not very massive compared to the BH’s mass. Table 3 lists the final value of in each of our BHNS simulations.

v.2 Convergence studies

Figure 6: Numerical convergence of gravitational wave signals for three B-SSep simulations. The highest grid resolutions are , , and , so that initially the AH diameter is covered by 38, 44, and 59 gridpoints, and the NS is covered by 58, 67, and 91 gridpoints, respectively. In the top panel, we show the real part of the component of for the three waveforms. In the bottom panel, pairwise differences between the waveforms are plotted and rescaled to demonstrate second-order convergence.

To demonstrate that the resolution used in our simulations is sufficiently high, we perform a convergence test. It is known that as the BH spin is increased, higher resolution is required for accurate BHBH evolutions (see e.g. Brügmann et al. (2008b)). Hence we perform a convergence test for the highest spinning case with . We use three different grid resolutions with close initial binary separation (Cases B-SSep-HRes, B-SSep-MRes and B-SSep-LRes in Table 2), which is adequate for the purpose of a convergence test. The top panel in Fig. 6 shows the real part of , the mode of , computed with three different resolutions, showing good agreement. The bottom panel in Fig. 6 compares the pairwise differences between the waveforms and we see that the waveforms converge at second-order accuracy.

v.3 Effects of different initial binary separations

In addition to the convergence test, we need to verify that our initial binary separation is sufficiently large. That is, starting the simulation from any larger initial separation should give the same outcome after a simple time (and, for gravitational waves, a possible phase) shift.

Figure 7: Rest mass of NS matter outside the BH versus time for different initial separations for the , nonspinning BH case (Cases A, A-MSep, and A-SSep). Here, the time coordinate is shifted by , where is the time at which of the NS rest mass has fallen into the apparent horizon.
Figure 8: Comparison of gravitational waveforms for Case A (black solid line) and A-MSep (blue dashed line).

Figure 7 shows the rest mass of NS matter outside the AH as a function of time for Cases A, A-MSep, and A-SSep. These cases correspond to initial binaries along the same quasiequilibrium sequence but with different initial separations. The time origin is shifted by , where denotes the time at which of the NS rest mass has fallen into the AH. We see that the results in Cases A and A-MSep agree very well, but the result for A-SSep deviates from the other two. Figure 8 shows the gravitational waveforms Re() emitted in Cases A and A-MSep. We again see that the waveforms agree very well after a simple time and phase shift. These results indicate that initial separations corresponding to those in Cases A-MSep and A are sufficiently large. For almost all other simulations, the initial value of is similar to or smaller than that for Case A (), suggesting that the initial binary separation should be sufficiently large in all of these cases. The only exception is Case E, for which initially. This value should be compared with that of Case A-SSep, which had initially. Therefore, our results for Case E may be inaccurate, since the initial separation may be too small. From an astrophysical perspective, our Case E may also be the least relevant, since the formation of a BH with a NS of equal mass seems very unlikely.

v.4 Effects of pressure ceiling

Figure 7 shows that a disk of rest mass remains outside the AH at the end of simulation, where is the rest mass of the initial NS. This finding disagrees with the result reported in Paper I. The reason for this discrepancy is due to artificially imposed limits on pressure that were applied everywhere in the fluid in Paper I. Here, we apply these limits only in regions where the density is tiny (i.e. inside the artificial “atmosphere”) or where the conformal factor is large (i.e., deep inside the AH), as discussed in Sec. III.2. Setting the pressure ceiling everywhere has little effect on the inspiral dynamics. However, when the NS is tidally disrupted and forms a relatively low-density “tail”, shock heating becomes a significant effect, causing to exceed (see Appendix B). In our original formulation, we reset to . Artificially lowering the pressure in this way forces the low-density material to lose a significant amount of torque and angular momentum and fall into the BH. Even though the density of the tail is low, it is still many orders of magnitude larger than the artificial atmospheric density, and hence has substantial influence on the dynamics of the merging NS.

To more accurately model the physics of the merger phase and maintain a stable evolution inside the BH, we now impose the limits on only in very low-density regions (where the rest-mass density ) or when . We set in most simulations, except Case E where we need to set to stabilize the evolution. We note that even for this relatively small , holds true only deep inside the AH, so this readjustment should not affect the dynamics outside the BH. By relaxing this artificial pressure ceiling, we now observe a small remnant disk in many cases and achieve significantly better angular momentum conservation. For example, for Case A we previously had but now after relaxing the pressure ceiling.

v.5 Effects of black-hole spin

Figure 9: Trajectories of the BH and NS coordinate centroids for Cases A (left) and B (right). The BH coordinate centroid corresponds to the centroid of the AH, and the NS coordinate centroid is computed via: , where is the total simulation volume. Note that this equation contains a different than Eq. (20).
Figure 10: Snapshots of density and velocity profiles at selected times for Case A. The contours represent the density in the orbital plane, plotted according to , (=0, 1, … 12), with darker greyscaling for higher density. The maximum initial NS density is , or . Arrows represent the velocity field in the orbital plane. The black hole AH interior is marked by a filled black circle. In cgs units, the initial ADM mass for this case is skm.
Figure 11: Snapshots of density and velocity profiles at selected times for Case B. The contours represent the density in the orbital plane, plotted according to (=0, 1, … 12), with darker greyscaling for higher density. The maximum initial NS density is , or . Arrows represent the velocity field in the orbital plane. We specify the black hole AH interior in each snapshot with a filled black circle. In cgs units, the initial ADM mass for this case is skm.
Figure 12: Snapshots of density and velocity profiles at (top row) and (bottom row) for Cases C, A, and B, which have initial BH spins of , , and , respectively. The contours represent the density in the orbital plane, plotted according to (=0, 1, … 12), with darker greyscaling for higher density. The maximum initial NS density is for all cases, or . Arrows represent the velocity field in the orbital plane. The black hole AH interior is marked by a filled black circle. In cgs units, the initial ADM mass for these cases is skm.
Figure 13: Rest-mass fraction outside the BH for different initial BH spins (Cases C, A, and B). Here, the time coordinate is shifted by , the time at which of the NS rest mass has fallen into the apparent horizon.

In Cases A–C, we vary the initial spin parameter of the black hole (parallel to the orbital angular momentum) from = to 0.75, fixing and the initial orbital angular velocity . Notice in Table 1 that as the initial BH spin parameter increases, the total initial angular momentum increases, requiring more gravitational-wave cycles to emit angular momentum and bring the BH and NS close enough to merge. Thus we expect the binary to undergo more orbits before merger as increases. This is precisely what we find: for spins , 0.0, and 0.75, the binary inspiral phase lasts for 3.25, 4.5, and 6.5 orbits, respectively. Figure 9 shows the coordinate trajectories of the AH and NS centroids for Cases A and B.

Figure 10 shows snapshots of rest-mass density contours for Case A. The upper panels show the binary at the start of the simulation, two orbits into the simulation, and 3.5 orbits into the simulation, corresponding to the time at which the first few density contours have crossed into the AH. Notice that the equilibrium shape of the NS is for the most part undisturbed after 2 orbits. The lower panels in Fig. 10 are snapshots taken at , , and . They demonstrate how the NS tail deforms into a quasistationary disk, as the bulk of the matter is accreted onto the hole.

Similarly, the upper plots of Fig. 11 demonstrate that for Case B, the initial NS (upper left) retains its shape after 4 orbits (upper middle), and begins shedding its outer layers due to tidal disruption at about 5 orbits (upper right). Notice however, in this case, the tail at (lower left) is quite massive, and at (lower center) it is much larger than the tail at in Case A (lower left plot in Fig. 10). The lower right plot is a snapshot taken near the end of the simulation, when a quasistationary disk resides outside the BH, at .

To compare the effects of BH spin for all cases in this study, we plot the density contours at and for Cases A–C in Fig. 12. In the case (Case C), the NS is basically “swallowed whole” by the BH during merger, leaving of the NS rest mass as a disk. As the spin increases, NS tidal disruption becomes more pronounced, resulting in long tidal tails that eventually form disks with rest mass and 15% the rest mass of the NS in Cases A and B, respectively (see Fig. 13). This result is not surprising. The ISCO decreases as the BH spin increases, and hence the tidal disruption of the NS occurs farther from the ISCO as the BH spin increases. Thus BHs with higher spin would likely lead to a more massive disk.

Figure 13 shows the rest mass outside the AH as a function of time. We see that there are two phases of matter falling into the BH: a plunge phase and an accretion phase. The plunge phase occurs early in the merger as part of the NS matter streams onto the BH and the rest deforms into a tidal tail. The plunge phase ends when 70%–90% of the NS matter falls into the BH, resulting in a sudden drop in the slope of the exterior vs. time plot in Fig. 13. The matter in the tail, having larger specific angular momentum, spreads out and forms a disk. Material with lower specific angular momentum in the disk accretes onto the BH. Since there is neither viscosity nor magnetic fields in our simulation, the accretion should eventually cease as the evolution continues. However, in realistic astrophysical environments, magnetic fields do exist and could have substantial influence on the dynamics of the disk on secular timescales. In addition, some material in the disk is shock heated to high temperature (see Sec. V.7). Copious amounts of neutrinos and antineutrinos may be generated as a result. Both magnetic fields and neutrino-antineutrino pairs could drive energetic jets, resulting in an SGRB (see Sec. V.7 for further discussion).

v.6 Effects of varying the mass ratio

Figure 14 demonstrates how different mass ratios alter the NS material falling into the BH. For the and cases (Cases E and A), we clearly identify the plunge and accretion phases mentioned in the previous subsection. However, for the case (Case D), the NS basically plunges into the BH, leaving of its rest mass outside the BH at the end of the simulation. This result is not surprising since, at the ISCO, the tidal effect of the BH is smaller for larger , resulting in tidal disruption occurring closer to the ISCO. Moreover, for a fixed NS compaction, the horizon size of the BH is larger for larger . Hence more NS matter is expected to fall into the BH.

Figure 14: Rest-mass fraction outside the BH for different BHNS mass ratios. Here, the time coordinate is shifted by , the time at which of the NS rest mass has fallen into the apparent horizon.
Figure 15: Snapshots of density and velocity profiles at (top row) and (bottom row) for Cases E, A, and D, which have initial mass ratios of , 3, and 5, respectively. The contours represent the density in the orbital plane, plotted according to (=0, 1, … 12), with darker greyscaling for higher density. The maximum initial NS density is for all cases, or . Arrows represent the velocity field in the orbital plane. The black hole AH is marked by a filled black circle. In cgs units, the total ADM mass for these cases is skm.
Figure 16: Snapshots of rest-mass density and contours for Case E ( mass ratio). The first and third rows show snapshots of density contours and velocity profiles in the orbital plane. The second and fourth rows show snapshots of . The density is plotted according to (=0, 1, … 12), with darker greyscaling for higher density. The maximum density of the initial NS is , or . Arrows represent the velocity field in the orbital plane. is plotted according to (=0, 1, … 12), with darker greyscaling for higher . These figures neglect in regions where the density is less than . The black hole AH interior is marked by a filled black circle. In cgs units, s=km.

Figure 15 shows snapshots of density and velocity profiles of the three cases at and . We again see different structures of the NS tail among the three cases. The merger of Case E ( case) is particularly interesting. Figure 16 shows snapshots of rest-mass density and contours. Since the BH is less massive in Case E, tidal disruption occurs at a farther binary separation (upper middle plots). The disrupted NS matter curls around the BH, forming a hot, low-density spiral (upper right plots) that winds around the AH and smashes into the tidal tail, generating a large amount of shock heating. Some of the heated NS matter transfers angular momentum to the other part of the tail and falls into the BH. The remaining matter in the tail deforms into an inhomogeneous disk (lower center plots) before settling into a quasistationary state, in which a high density, relatively low temperature torus of matter surrounds the BH (lower right plots). Analysis of the accretion versus time (Fig. 14) demonstrates that the plunge phase is relatively long compared to Case A. Because of the complicated interaction of the disrupted NS matter in the merging process, the disk mass is less than that of Case A.

v.7 Disk profile

Figure 17: Snapshots of rest-mass density , , and velocity profiles at the end of the simulation, for Cases E (left column), A (middle column), and B (right column) on the orbital (, first and third rows) and meridional ( and , second and fourth rows) planes. In the top six plots, density contours are shown according to (=0, 1, … 13) for Cases A and B, (=0, 1, … 12) for Case E. The maximum density of the initial NS is for all cases, or . The bottom six plots show contours of , plotted according to (=0, 1, … 12). The figures neglect regions where the density is less than . All contour plots show darker greyscaling for higher and . Arrows represent the velocity field in the given plane. The AH interior in the orbital plane is marked by a filled black circle. Length scales are specified in km, assuming the NS has a rest mass of 1.4; it can be converted to units of via the formula km.

In Cases A, B, and E, a substantial disk forms after the merger. Figure 17 shows the density and thermal energy profiles of the disk in both the equatorial and meridional plane at the end of the simulation for Cases A, B, and E. For ease of comparison, we specify the length scale in km, assuming the NS has a rest mass of 1.4; it can be converted to units of via the formula km. The most distinguishable features of the disks are their characteristic radius and final mass, . The characteristic disk radii lie between 20km and 100km, as seen in the top row of plots. The mass ratio (Case E) produces the smallest (km), densest (maximum density ), and least massive disk (). Meanwhile, the resulting disk in the canonical case (Case A) is about 2.5 times larger than Case E (km), but possesses a lower maximum density by an order of magnitude (). Due in part to its larger volume, the mass of the disk in Case A is about 70% larger than in Case E (). The maximum density of Case B’s disk is , which is similar to Case A. However, the disk is about twice as large in size, with a characteristic radius of km, yielding a total disk mass of .

Despite their differences, the disks in Fig. 17 share many common features. For example, the top row of plots show that the disk forms a torus whose density plummets at the BH ISCO (roughly at coordinate radius 7km, 25km, and 12km for Cases E, A, and B, respectively, in the equatorial plane ft2 ()). Also, the characteristic height of the disks is about 15%–20% of the characteristic radius in each case (top two rows of plots). Finally, the higher density regions tend to be colder than lower density regions, but everywhere in the disks in Cases A and B, and everywhere in Case E’s disk. The characteristic in the disks in Cases A, B, and E are roughly 200, 200, and 10, respectively.

As in Paper I, the temperature can be estimated roughly by modeling the temperature dependence of the specific thermal internal energy density as


(compare Popham et al. (1999)), where is the mass of a nucleon, is the Boltzmann constant, and is the radiation constant. The first term represents the approximate thermal energy of the nucleons, and the second term accounts for the thermal energy due to radiation and (thermal) relativistic particles. The factor reflects the number of species of ultrarelativistic particles that contribute to thermal radiation. When K, where is the mass of electron, thermal radiation is dominated by photons and . When , electrons and positrons become ultrarelativistic and also contribute to radiation, and . At sufficiently high temperatures and densities (K, ), thermal neutrinos are copiously generated and become trapped, so, taking into account three flavors of neutrinos and anti-neutrinos, .

The characteristic density and in the disk for Cases A and B are (or ) and . We find (or 4MeV). For Case E, we have (or ). We find K (or 1MeV) for Case E’s disk. According to numerical models of rotating BHs with ambient disks in Setiawan et al. (2006), the remnant disks could produce a total gamma-ray energy erg from neutrino-antineutrino annihilation. This result is promising for generating a SGRB.

v.8 Gravitational wave emission

Figure 18: Gravitation wave signal from Case A. Shown here are the (2,2) and (3,3) modes of and extracted at radii (black solid lines), (red dotted lines) and (blue dashed lines). Note that the three lines almost overlap in each case.
Figure 19: Gravitational wave signal from (top to bottom) Cases C, A and B. Black solid (blue dash) line denotes the (2,2) mode of () extracted at .
Figure 20: Same as Fig. 19 but for cases (top to bottom) E, A and D.
Figure 21: Gravitational wave spectrum for the Cases C, A and B, computed via the method summarized in Sec. IIIF of Paper I. The solid curve shows the combined waveform found by attaching the restricted 2.5 order PN waveform to our numerical signal (including only the dominant and modes), while the dotted curve shows the contribution from the latter only. The dashed curve is the analytic fit derived by Ajith et al. (2008) from analysis of multi-orbit nonspinning BHBH binaries with the same mass ratios as the BHNS, which maintain significantly more power at higher frequencies. The heavy solid curve is the effective strain of the Advanced LIGO detector, defined such that . To set physical units, we assume a NS rest mass of and a source distance of =100Mpc.
Figure 22: Same as Fig. 21 but for Cases E, A and D.

Following the literature, we decompose the gravitational waveform and into spin-weighted spherical harmonics as follows


where and are real functions. Each mode is a function of radius and time only. We extract gravitational waves using both the Regge-Wheeler-Zerilli-Moncrief formalism and the Newman-Penrose Weyl scalar . We find good agreement between waveforms computed by both methods, as demonstrated in Paper I. Here we only present the waveforms computed from .

Gravitational waves are extracted at several different radii. We find that the extracted waveforms overlap extremely well provided the extraction radius , as shown in Fig. 18.

Figures 19 and 20 demonstrate that the maximum amplitude of the waves is similar in Cases A, B, C, D, and E. With approximately the same initial , we observe more wave cycles as the BH spin is increased (Fig. 19) and as the mass ratio is decreased (Fig. 20), consistent with the discussion in Secs. V.5 and V.6. Table 3 specifies the radiated energy and angular momentum in each case, as well as the kick velocity due to recoil.

For Case E, NS disruption occurs at relatively low orbital frequency. In fact, the total radiated angular momentum reaches 90% of its final value at (the time at which only 10% of the NS has been accreted), and 99% of its final value at . Hence there is a large amount of angular momentum still remaining in the NS matter before it is accreted, and during the accretion phase about 98% of the NS rest mass is accreted. Thus the bulk of the orbital angular momentum is used to spin up the BH. The final disk mass is only % of the binary’s ADM mass, so we can reliably estimate the final value of the BH spin parameter using the Kerr formula (Eq. (18)). We obtain , which is significantly larger than all other cases that result in a similarly small disk. Compared to the final spin parameter from equal-mass BHBH mergers, (e.g. Scheel et al. (2009) and Table 7 in Appendix A), we find that matter accretion in BHNS binary mergers is much more effective at spinning up black holes. This is not unexpected, since NS tidal disruption spreads the NS matter around the BH, diminishing the multipole moments of the system, and the associated gravitational wave emission.

Figures 21 and 22 show the effective gravitational wave strains in frequency space for all cases, computed via the method summarized in Sec. III F of Paper I, and including only the dominant and modes. For comparison against the Advanced LIGO sensitivity curve (Gustafson et al. (1999)), we plot BHNS wave strains for the cases we study, assuming the NS has a rest mass of and binary distance of 100Mpc. This is the distance required to reach one merger per year, assuming an overall rate of 10 mergers per Myr per Milky Way-equivalent galaxy (and a density of Belczynski et al. (2002). This distance is roughly that of the Coma cluster, and approximately five times the distance to the Virgo cluster. The gravitational wave spectra of nonspinning BHBH mergers (computed from Eqs. (4.12)–(4.19) of Ajith et al. (2008)) with the same mass ratio as in Cases E, A and D are plotted in the same graphs for comparison. We see that the wave signal drops significantly at and above the frequency corresponding to the onset of NS tidal disruption. The difference in wave signals between BHBH and BHNS mergers is marginally observable in the Advanced LIGO frequency band in most cases. Distinguishing BHNS from BHBH inspirals and mergers may require narrow-band detection techniques with advanced detectors. The observation of an accompanying SGRB would also serve to distinguish the events.

Vi Summary and Discussion

In this paper we perform a new set of fully general relativistic BHNS simulations using our code upgraded with AMR capability. We vary the initial binary separations and confirm that the simulation outcomes do not change when the initial BHNS separation is large enough. In most cases, a toroidal disk forms around the black hole remnant after the merger. This result revises our finding of Paper I, where we artificially imposed a uniform pressure ceiling on all the matter present in the simulation. This ceiling removed a substantial amount of angular momentum in the NS matter during the merger phase, causing the matter to fall into the BH prematurely, and suppressing disk formation.

Fixing the mass ratio , we study the effects of BH spin by evolving initial data with BH’s spin parameter (count-rotating), 0 and 0.75 with nearly the same initial orbital angular velocity . Not surprisingly, we find that the BHNS inspiral phase lasts longer as increases, and the final disk mass increases from () to () and () of the NS’s rest mass.

Next, we study the effect of varying the mass ratio by evolving BHNS initial data with and 5. For the case almost all the NS matter plunges into the BH, leaving of the NS’s rest mass to form a disk at the end of the simulation. For the case, a low-density, hot spiral region of disrupted NS matter winds around the BH and smashes into the nascent NS tidal tail, causing strong shock heating before it rapidly falls into the BH. Most of the remaining matter in the tail eventually falls into the BH, leaving about 2% of the NS’s rest mass behind to form a disk.

The disks formed after the merger of Cases A and B have a characteristic density of and a temperature of K, assuming the NS’s rest mass is . Applying the results of the simulations of BH disks in Setiawan et al. (2006), the disk may produce a gamma-ray energy of erg due to neutrino-antineutrino annihilation. Hence the merger of a BHNS binary is a promising central engine for a SGRB. But self-consistent simulations, taking into account the detailed microphysics and including the correct EOS, are necessary to fully explore the viability of and the mechanism for generating SGRBs from BHNS mergers.

We compute the BHNS gravitational waveforms and power spectra, and find that the signals are attenuated at frequencies roughly equal to double the orbital frequency at which tidal disruption begins. The resulting deviation of BHNS and BHBH waveforms is visible in the Advanced LIGO high frequency band out to distances . Should the chirp mass determination, combined with higher order PN waveform phase effects, allow for an independent determination of the component masses and spins of the binary companions, the measurement of the BHNS tidal disruption frequency should give a good estimate of the NS radius and, hence, insight into the nuclear EOS.

Carpet-based AMR has provided us with a great improvement in efficiency and accuracy, and future work will continue in this direction. To improve code performance (which is usually memory-bound), we also plan to employ new techniques (e.g., loop-tiling) for minimizing cache misses. To improve accuracy, we will implement the piecewise parabolic method (PPM) or weighted essentially non-oscillatory (WENO) reconstruction schemes in hydrodynamics/MHD and experiment with tapered grids to replace second order temporal prolongation (see Lehner et al. (2006)). We plan to study techniques for minimizing refinement boundary crossings, by adding more refinement boxes around the disrupted NS. This might help reduce our angular momentum conservation violation to below 1% for all cases.

Rich with the challenges of modeling extreme matter in the presence of intense gravitational fields, BHNS binary mergers will likely remain a topic at the forefront of theoretical astrophysics and numerical relativity for years to come. The most exciting possibility remains the simultaneous detection of a gravitational wave signal and a GRB. Our theoretical work on BHNSs is partly in anticipation and preparation for that discovery.

We would like to thank K. Taniguchi for generating the BHNS initial data. We would also like to thank M. Ansorg for providing the TwoPunctures code for generating BHBH initial data, and E. Schnetter for useful discussions about Carpet. This paper was supported in part by NSF Grants PHY02-05155 and PHY06-50377 as well as NASA Grants NNG04GK54G and NNX07AG96G to the University of Illinois at Urbana-Champaign, and NSF Grant PHY07-56514 to Bowdoin College. Simulations were performed under a TeraGrid Grant TG-MCA99S008 and on the Illinois Numerical Relativity Beowulf Cluster.

Appendix A Details of BHBH test simulations

Simulating the late inspiral, merger and ringdown of a binary black hole was for many years the “holy grail” of numerical relativity. Dozens of researchers have spent many years attempting to formulate a stable algorithm capable of solving this important problem. In this appendix we summarize two BHBH simulations in detail. Performing a simulation of binary black hole coalescence provides a highly nontrivial and reasonably comprehensive laboratory for testing 3+1 codes that solve Einstein’s vacuum field equations in the strong-field regime of general relativity. It tests the caliber of not only the basic evolution scheme, but also the all-important diagnostic routines. These routines measure globally conserved quantities like the mass and angular momentum of the system, the location and measure of all black hole horizons, the asymptotic gravitational waveforms, the recoil motion of the black hole remnant, etc. Appreciable effort must be expended to implement reliable diagnostic routines in order to extract useful physical information from the numerical output.

We consider two cases of merging, nonspinning black holes to test our code: an equal-mass () and an unequal-mass () system. Both cases are constructed via the binary puncture technique discussed in Brandt and Brügmann (1997), with parameters as listed in Table 4. Our code employs AMR with equatorial symmetry using the moving-box Carpet infrastructure, so that the boxes track the AH centroids throughout the simulations. The grid setup is described in Table 5. The basic evolution scheme and numerical techniques are described in Secs. II and III. A review of the special algorithms and parameters chosen in these simulations is summarized in Table 6, and the results are summarized in Table 7. The coordinate trajectories of the two black holes are plotted for the two cases in Fig. 23, and the gravitational waveforms are plotted in Fig. 24.

The results found here are in good agreement with those reported by previous investigators Baker et al. (2006, 2007); Berti et al. (2007); González et al. (2007). Of special significance is the degree to which total energy and angular momentum are conserved when properly accounting for losses due to gravitational wave emission. The fractional errors are seen in Table 7 to be a few times for the energy and for the angular momentum. Also, different measures of the mass and spin of the final (stationary) remnant Kerr black hole, such as the irreducible mass and the ADM mass, agree closely, given the adopted computational resources.

Equal mass Unequal mass
Location of punctures