Gravitational waves from spinning black hole-neutron star binaries: dependence on black hole spins and on neutron star equations of state

# Gravitational waves from spinning black hole-neutron star binaries: dependence on black hole spins and on neutron star equations of state

## Abstract

We study the merger of black hole-neutron star binaries with a variety of black hole spins aligned or antialigned with the orbital angular momentum, and with the mass ratio in the range –5, where and are the mass of the black hole and neutron star, respectively. We model neutron-star matter by systematically parametrized piecewise polytropic equations of state. The initial condition is computed in the puncture framework adopting an isolated horizon framework to estimate the black hole spin and assuming an irrotational velocity field for the fluid inside the neutron star. Dynamical simulations are performed in full general relativity by an adaptive-mesh refinement code, SACRA. The treatment of hydrodynamic equations and estimation of the disk mass are improved. We find that the neutron star is tidally disrupted irrespective of the mass ratio when the black hole has a moderately large prograde spin, whereas only binaries with low mass ratios, , or small compactnesses of the neutron stars bring the tidal disruption when the black hole spin is zero or retrograde. The mass of the remnant disk is accordingly large as , which is required by central engines of short gamma-ray bursts, if the black hole spin is prograde. Information of the tidal disruption is reflected in a clear relation between the compactness of the neutron star and an appropriately defined “cutoff frequency” in the gravitational-wave spectrum, above which the spectrum damps exponentially. We find that the tidal disruption of the neutron star and excitation of the quasinormal mode of the remnant black hole occur in a compatible manner in high mass-ratio binaries with the prograde black hole spin. The correlation between the compactness and the cutoff frequency still holds for such cases. It is also suggested by extrapolation that the merger of an extremely spinning black hole and an irrotational neutron star binary does not lead to the formation of an overspinning black hole.

###### pacs:
04.25.D-, 04.30.-w, 04.40.Dg

## I introduction

Coalescing binaries composed of a black hole (BH) and/or a neutron star (NS) are among the most promising sources of gravitational waves for ground-based laser-interferometric gravitational-wave detectors such as LIGO ligo2009 () and VIRGO virgo2011 (). Detections of gravitational waves will be accomplished in a decade to come by planned next-generation detectors such as advanced LIGO, advanced VIRGO, and LCGT lcgt2010 (). Because gravitational waves are much more transparent to the absorption and scattering by the material than electromagnetic waves and even neutrinos are, gravitational-wave astronomy is expected to become a powerful and unique way to observe strongly gravitating phenomena in our Universe. Among such phenomena, the merger of a BH-NS binary plays an important role to investigate properties of the NS such as the radius and the equation of state (EOS) of a high-density nuclear matter lindblom1992 (); vallisneri2000 (); rmsucf2009 (); fgp2010 (); ptr2011 (); pror2011 (). An important constraint on the EOS is obtained from detection of a NS, which is the most massive NS currently known dprrh2010 (), by recent pulsar-timing observation. However, we still do not know the realistic EOS of the NS because there is no robust measurement of the NS radius slb2010 (). To constrain the NS radius and EOS by observing gravitational waves from the BH-NS binary, we have to prepare accurate theoretical templates of gravitational waveforms employing a wide variety of the NS EOSs and other physical parameters. For this purpose, numerical relativity is the unique approach.

The merger of a BH-NS binary is also an important target for the astrophysical study because it is a potential candidate for the progenitor of short-hard gamma-ray bursts (GRBs) in the so-called merger scenario (see Refs. nakar (); leeramirezruiz () and references therein for reviews). If a NS is tidally disrupted during the merger of a BH-NS binary, a system composed of a spinning BH and a hot, massive accretion disk of may be formed. Such an outcome could be a central engine of the GRB, because it could radiate a large amount of energy erg via neutrino emission or the Blandford-Znajek mechanism blandfordznajek1977 () in a short time scale of s and hence could launch a GRB jet. The merger scenario of the GRB is attractive when a short-hard GRB is associated with a galaxy of low star-formation rate bergeretal2005 (); fbctlfgcf2010 () because the collapsar model of GRBs woosley1993 () is not preferable. Only numerical relativity can answer quantitatively the question whether the formation of a massive accretion disk is possible in the BH-NS binary merger.

Fully general relativistic study of BH-NS binaries has achieved progress in recent years, both in computations of quasiequilibrium states grandclement2006 (); ?; tbfs2007 (); tbfs2008 (); fkpt2008 (); kst2009 () and in dynamical simulations of the merger shibatauryu2006 (); shibatauryu2007 (); shibatataniguchi2008 (); eflstb2008 (); dfkpst2008 (); elsb2009 (); skyt2009 (); dfkot2010 (); kst2010 (); cabllmn2010 (). Although these include several simulations of spinning BH-NS binaries with a qualitative ideal-gas EOS, to date, only limited number of simulations have been performed taking into account both the nuclear-theory based EOS and the BH spin dfkot2010 (). In particular, we still do not understand the dependence of the merger process and resulting gravitational waveforms on the BH spin and the EOS of the NS in detail. One goal of current numerical relativity is to clarify the effect of the BH spin on the BH-NS binary merger and subsequent NS tidal disruption adopting a wide variety of the NS EOSs.

In this paper, we report our latest results obtained by numerical-relativity simulations with a variety of the NS EOSs and the BH spins. We employ five piecewise polytropic EOSs (see Sec. II.2), all of which do not conflict with the current observation of the NS. We systematically choose physical parameters such as BH mass, NS mass, and BH spin in an astrophysically realistic range. While we only consider relatively low mass BHs in the previous work kst2010 (), we adopt a wider range of the BH mass because a large BH spin enhances the NS tidal disruption for high mass-ratio binaries. We clarify the dependence on the BH spin and on the NS EOS of the properties of the merger remnants and characteristics of gravitational waves—in particular, the gravitational-wave spectrum.

This paper is organized as follows. In Sec. II, we describe methods for a solution of initial conditions, piecewise polytropic EOSs, and models of BH-NS binaries employed in this paper. In Sec. III, the formulation and methods of numerical simulations are summarized. Section IV presents the numerical results and clarifies the effect of the BH spin and NS EOS on the tidal disruption, merger remnants, and gravitational waves. Section V is devoted to a summary. Throughout this paper, we adopt the geometrical units in which , where and are the gravitational constant and the speed of light, respectively. Our convention of notation for physically important quantities is summarized in Table 1. The nondimensional spin parameter of the BH, total mass of the system at infinite separation, mass ratio, and compactness of the NS are defined as , , , and , respectively. Latin and Greek indices denote spatial and spacetime components, respectively.

## Ii Initial condition

As in our previous works skyt2009 (); kst2010 (), we employ BH-NS binaries in quasiequilibrium states for initial conditions of our numerical simulations. In this section, we summarize the formulation and methods for the computation of a quasiequilibrium state, specifically our method of estimating the spin angular momentum of the BH in a binary and of determining the position of the rotation axis. The details of the formulation and numerical methods, except for the issues on the BH spin, are described in Ref. kst2009 (), to which the reader may refer. Computations of the quasiequilibrium states are performed using the spectral-method library LORENE LORENE ().

### ii.1 Formulation and methods

We compute a quasiequilibrium state of the BH-NS binary as a solution of the initial value problem of general relativity cook (). As far as the orbital separation is large enough, the time scale of the orbital contraction due to the gravitational radiation reaction is much longer than the orbital period, and, therefore, we may safely neglect the gravitational radiation reaction in the calculation of the quasiequilibrium state. In numerical simulations of the binary coalescences, we have to track orbits in order to calculate accurate gravitational waveforms during the late inspiral and merger phases, and hence, the orbital separation of the initial condition has to be large enough. For such initial conditions, we can neglect the gravitational radiation reaction. Thus, we give a BH-NS binary in a quasicircular orbit, i.e., the binary in an approximate equilibrium state in the corotating frame. To satisfy the quasiequilibrium requirements described above, we assume the existence of a helical Killing vector with the orbital angular velocity ,

 ξμ=(∂t)μ+Ω(∂φ)μ. (1)

We also assume that the NS is in the hydrostatic equilibrium in the corotating frame and has an irrotational velocity field, which is believed to be a reliable approximation to an astrophysically realistic configuration bildstencutler1992 (); kochanek1992 ().

We compute the three-metric , the extrinsic curvature , the lapse function , and the shift vector by a mixture of the extended conformal thin-sandwich approach york1999 (); pfeifferyork2003 () and the conformal transverse-traceless decomposition cook () in the puncture framework brandtbrugmann1997 (); clmz2006 (); bcckm2006 (). In this formalism, we assume the conformal flatness of the three-metric , the stationarity of the conformal three-metric , the maximal slicing condition , and its preservation in time , where is the flat spatial metric. Assuming that the puncture is located at , we set the conformal factor and a weighted lapse function as

 ψ=1+MP2rBH+ϕ,Φ=1−MΦ2rBH+η, (2)

where and are positive constants of mass dimension, and is a coordinate distance from the puncture. We adjust to obtain a desired mass of the BH, , and determine so as to satisfy the virial relation, i.e., the equality of the Arnowitt-Deser-Misner (ADM) mass and the Komar mass, which holds when the spacetime is stationary and asymptotically flat beig1978 (); ashtekarashtekar1979 (). , , and are determined by solving elliptic equations derived from the Hamiltonian constraint, momentum constraint, and quasiequilibrium conditions, and . We note that these quasiequilibrium conditions can be replaced by and in a conformal flatness approximation. In the puncture framework, we also decompose a conformally weighted extrinsic curvature as

 ^Aij=^∇iWj+^∇jWi−23fij^∇kWk+KPij, (3)

where is an auxiliary three-vector field and is the covariant derivative associated with . is a singular part of the extrinsic curvature, which is associated with the linear and spin angular momenta of the BH bowenyork1980 (),

 KPij = 32r2BH[liPBHj+ljPBHi−(fij−lilj)lkPBHk] (4) + 3r3BH[ϵkilSlPlklj+ϵkjlSlPlkli],

where is a unit radial vector, , and is the Levi-Civita tensor associated with the flat metric . and are parameters associated with the linear and spin angular momenta of the BH, respectively. Here we determine by the condition in which the total linear momentum of the binary vanishes, and we adjust to obtain a desired spin angular momentum of the BH, . The elliptic equation to determine is obtained by taking the derivative of Eq. (3) and using the momentum constraint.

The spin angular momentum of the BH, , is evaluated on the apparent horizon (hereafter AH), , according to the isolated horizon framework (see Refs. ashtekarkrishnan (); gourgoulhonjaramillo () for reviews). Because we do not know the position of the AH in advance in the puncture framework, we have to determine the location of the AH numerically linnovak2007 (). On the numerically determined AH, an approximate rotational Killing vector may be defined using the method developed by Cook and Whiting cookwhiting2007 () with the normalization condition proposed by Lovelace and his collaborators lopc2008 (). We focus only on the case in which the BH spin is aligned or antialigned with the angular momentum of the binary in this work, and hence, the axis of the BH spin is uniquely determined. Hereafter we set the rotational axis of the binary, equivalently the axis of the BH spin, to be the z axis and consider only the approximate Killing vector associated with the rotation in this direction. Using , we obtain the spin angular momentum of the BH via the surface integral at the AH,

 S(ϕ)BH=18π∫SKijϕidSj. (5)

We adjust to obtain a desired value of (hereafter ). We note that and do not agree exactly in the BH-NS binary spacetime due to the contribution to the extrinsic curvature from the NS, associated with .

Because we adopt a conformal flatness approximation for the induced metric, the Christodoulou mass of the BH evaluated on the AH, , and the gravitational mass evaluated at spatial infinity, , do not agree even for a single BH system due to the presence of so-called junk waves. This difference leads to an ambiguity in defining the nondimensional spin parameter of the BH. Here, we define the nondimensional spin parameter of the BH with respect to the mass evaluated at spatial infinity, i.e.,

 a≡SBHM2BH. (6)

The reason for this is that the mass and nondimensional spin parameter of the BH evaluated at the AH quickly (in our simulations, within ms) relax to and , defined at spatial infinity, respectively, as the BH absorbs the junk radiation in the vicinity of the BH dlz2008 (); lopc2008 (). We note that these values show the damping oscillation before the relaxation in the same manner as the “scalar-curvature spin” of Ref. lopc2008 () shows, because our method of evaluating these values in the simulation is basically the same as the method to define the scalar-curvature spin in Ref. lopc2008 () (see Sec. III.2).

To compute the equations of hydrostatic equilibrium for the NS matter, we assume an ideal fluid, for which the energy-momentum tensor is described by

 Tμν=ρhuμuν+Pgμν, (7)

where is the rest-mass density, is the pressure, is the specific enthalpy, is the specific internal energy, and is the four-velocity of the fluid. Basic equations for the hydrostatics are derived from the condition of irrotation, or the vanishing of the vorticity two-form bgm1997 (); asada1998 (); shibata1998 (); teukolsky1998 (),

 ωμν=∇μ(huν)−∇ν(huμ)=0, (8)

and the conservation of the specific momentum of the fluid along the helical Killing vector field . The specific enthalpy is determined from the first integral of the relativistic Euler equation (relativistic Bernoulli integral), and other thermodynamical quantities are subsequently obtained by using the EOS, which is described in Sec. II.2. The four-velocity of the fluid is calculated using the velocity potential , as in the assumption of the irrotational velocity field. The elliptic equation for is derived from the equation of continuity, .

We have no definite condition to determine the location of the center of mass of the binary in the puncture framework and we use this ambiguity to reduce an unphysical initial orbital eccentricity. We found in our previous work skyt2009 (); kst2009 () that orbits with a small eccentricity could be obtained using the “3PN-J method,” i.e., a phenomenological method to determine the location of the rotational axis in which the total angular momentum of the binary for a given value of agrees with that calculated from the third post-Newtonian (3PN) approximation. In the present case in which the BH has a finite amount of spin angular momentum, we extend the previous 3PN-J method to include the contribution from the BH spin up to the 2.5PN order fbb2006 (); bbf2006 (). Specifically, the location of the rotational axis is chosen from the condition that the orbital angular momentum of the binary agrees with a sum of nonspin terms given by Eq. (4) of Ref. blanchet2002 () and of spin terms given by Eq. (7.10) of Ref. bbf2006 () for a given value of . We show in Sec. IV.1 that this extended 3PN-J method again leads to the initial condition with a small eccentricity.

### ii.2 Piecewise polytropic equations of state

The matter inside the NS in the late inspiral phase is believed to be well-approximated by a zero-temperature nuclear matter because the cooling time scale of the NS in typical BH-NS binaries is shorter than the time scale of the gravitational radiation reaction lattimerprakash (). Hence, we employ a cold EOS, for which the rest-mass density, , determines all other thermodynamical quantities, for calculating the quasiequilibrium state of the BH-NS binary. To model nuclear-theory-based EOSs at high density with a small number of parameters, we employ a piecewise polytropic EOS. It is a phenomenologically parametrized EOS of the form

 P(ρ)=κiρΓi  for  ρi−1≤ρ<ρi  (1≤i≤n), (9)

where is the number of the pieces used to parametrize an EOS, is the rest-mass density at the boundary of two neighboring th and th pieces, is the polytropic constant for the th piece, and is the adiabatic index for the th piece. Here, , , and other parameters are freely chosen. Requiring the continuity of the pressure at each , free parameters—say —determine the EOS completely. The specific internal energy, , and hence the specific enthalpy, , are determined by the first law of thermodynamics and continuity of each variable at boundary densities, .

It was shown that piecewise polytropic EOSs with four pieces approximately reproduce most properties of the nuclear-theory-based EOSs at high density rlof2009 (). If we focus on low mass NSs with relatively low central density, the EOS at high density plays a minor role. Thus, we adopt a simplified piecewise polytropic EOS composed of two pieces, one of which models the crust EOS and the other of which the core EOS. This simplification is based on the fact that NSs in the observed binary NSs often have fairly small masses stairs2004 () and the maximum rest-mass density in such NSs may not be so high that the EOS at high density plays only a minor role in determining their structure. Furthermore, the maximum rest-mass density inside the NS should only decrease during the evolution of the BH-NS binary due to tidal elongation of the NS by the companion BH.

Table 2 lists the EOSs which we employ in this study. Following Refs. rmsucf2009 (); kst2010 (), we always fix the EOS for the crust region by the parameters below:

 Γ1 = 1.35692395, (10) κ1/c2 = 3.99873692×10−8(g/cm3)1−Γ1. (11)

The EOS of the core region is determined by two parameters. One is the adiabatic index of the core EOS, . Whereas we find that properties of gravitational waves and the merger remnants depend on in our previous study kst2010 (), we always fix in this work to focus on clarifying the effect of the BH spin aside from the difference in the adiabatic index. The other parameter is chosen to be the pressure at a fiducial density because is closely related to the radius and deformability of the NS lattimerprakash2001 (). We vary the value of systematically to investigate the effect of the stiffness of the EOS 1. With the given values of and , and are determined as

 κ2 = pρ−Γ2fidu, (12) ρ1 = (κ1/κ2)1/(Γ2−Γ1). (13)

### ii.3 Models

Numerical simulations are performed for a wide range of nondimensional BH spin parameter, , as well as for a variety of the mass ratios, . For nonspinning BH-NS binaries, we already found that the low mass ratio of is required for tidal disruption of NSs to occur sufficiently outside the innermost stable circular orbit (ISCO) of the BH unless the EOS is extremely stiff kst2010 (). If the tidal disruption occurs inside or at an orbit very close to the ISCO, we do not see strong effects of the tidal disruption. In such cases, gravitational waveforms are similar to those of a BH-BH binary even in the merger phase, and the mass of the remnant disk is negligible lksbf (). However, the allowed range of the mass ratio for the tidal disruption is modified drastically for a BH-NS binary with the prograde BH spin shibata1996 (); wigginslai2000 () because the ISCO radius 2 of the BH with a prograde spin becomes smaller by a factor of 1 – 6 bpt1972 () than that of the nonspinning BH with the same mass. Strong spin effects for the tidal disruption are also found in the numerical-relativity simulation of the spinning BH-NS binary merger with a simplified, -law EOS elsb2009 (). In this paper we perform a more systematic study of the tidal disruption for different EOSs, masses of each component, and BH spins.

Table 3 summarizes several key quantities for the initial conditions in our numerical simulations. The label for the model denotes the EOS name, the mass ratio, the NS mass, and the nondimensional spin parameter of the BH. Specifically, “a75,” “a5,” and “a-5” correspond to the spin parameters , , and , respectively. For example, HB-Q3M135a5 means that the EOS is HB and . Although we vary the NS mass systematically, the results of the merger remnant are reported only for binaries with in this paper because the difference in the NS mass complicates the properties of the remnant, such as the mass of the disk. Results for are analyzed only for gravitational waves.

For the same value of the mass ratio, we basically prepare the initial conditions with the same value of the initial angular velocity normalized by the total mass of the binary, . For 2H EOS, in which the NS radius is the largest, we exceptionally adopt a smaller value of than for other EOSs to guarantee orbits before tidal disruption occurs. The reason for this is that the tidal disruption occurs for a large orbital separation in 2H EOS. When the BH has a prograde spin, the number of orbits to the merger for a given value of increases due to spin-orbit repulsive interaction blanchet (), compared to the nonspinning BH case. On the other hand, when the BH has a retrograde spin, the number of orbits decreases due to spin-orbit attractive interaction. For , the number of orbits is typically by orbit smaller than for . For this reason, we also prepare the initial condition with a smaller value of for H EOS and .

## Iii Methods of simulations

Numerical simulations are performed using an adaptive-mesh refinement (AMR) code SACRA yst2008 (). The formulation, the gauge conditions, the numerical scheme, and the methods of diagnostics are basically the same as those described in Ref. kst2010 (), except for the correction in the treatment of hydrodynamic equations in a far region. Thus, we here only briefly review them and describe the present setup of the computational domain for the AMR algorithm and grid resolution.

### iii.1 Formulation and numerical methods

SACRA solves the Einstein evolution equations in the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism shibatanakamura1995 (); baumgarteshapiro1998 () with the moving-puncture gauge brandtbrugmann1997 (); clmz2006 (); bcckm2006 (). It evolves a conformal factor , the conformal metric , the trace of the extrinsic curvature , a conformally weighted trace-free part of the extrinsic curvature , and an auxiliary variable . Introducing an auxiliary variable and a parameter , which we typically set to be in units of , we employ a moving-puncture gauge in the form bghhst2008 ()

 (∂t−βj∂j)α = −2αK, (14) (∂t−βj∂j)βi = (3/4)Bi, (15) (∂t−βj∂j)Bi = (∂t−βj∂j)~Γi−ηsBi. (16)

We evaluate the spatial derivative by a fourth-order central finite difference, except for the advection terms, which are evaluated by a fourth-order noncentered, upwind finite difference, and employ a fourth-order Runge-Kutta method for the time evolution.

To solve the hydrodynamic equations, we evolve , , and . The advection terms are handled with a high-resolution central scheme by Kurganov and Tadmor kurganovtadmor2000 () with a third-order piecewise parabolic interpolation for the cell reconstruction. For the EOS, we decompose the pressure and specific internal energy into cold and thermal parts as

 P=Pcold+Pth,ε=εcold+εth. (17)

We calculate the cold parts of both variables using the piecewise polytropic EOS from the primitive variable , and then the thermal part of the specific internal energy is defined from as . Because vanishes in the absence of shock heating, is regarded as the finite-temperature part. In this paper, we adopt a -law ideal-gas EOS for the thermal part,

 Pth=(Γth−1)ρεth, (18)

to determine the thermal part of the pressure, and choose equal to the adiabatic index in the crust region, , for simplicity.

Because the vacuum is not allowed in any conservative hydrodynamic scheme, we put an artificial atmosphere of a small density outside the NS in the same way as done in our previous work kst2010 (). The total rest mass of the atmosphere is always less than , and hence, we can safely neglect spurious effects by accretion of the atmosphere onto the remnant disk as long as the disk mass is much larger than .

### iii.2 Diagnostics

We extract gravitational waves by calculating the outgoing part of the Weyl scalar at finite coordinate radii and by integrating twice in time as

 h+(t)−ih×(t)=−∫tdt′∫t′dt′′Ψ4(t′′). (19)

In our previous works skyt2009 (); kst2010 (), we directly perform this integration of and then subtract a quadratic function of the form to eliminate unphysical drift components in the waveform, using the least-square fitting to obtain constants , , and . In this work, we adopt a “fixed-frequency integration” method proposed by Reisswig and Pollney reisswigpollney2010 () to obtain gravitational waveforms with less unphysical components. In this method, we first perform a Fourier transformation of as

 ~Ψ4(ω)=∫dtΨ4(t)eiωt. (20)

Using this, Eq. (19) is rewritten as

 h+(t)−ih×(t)=12π∫~Ψ4(ω)ω2e−iωtdω. (21)

We then replace of the integrand with for , where is a positive free parameter in this method. By appropriately choosing , this procedure suppresses unphysical, low-frequency components of gravitational waves. As proposed in Ref. reisswigpollney2010 (), we choose to be for mode gravitational waves, where is the azimuthal quantum number. For the mode gravitational waves, we adopt and confirm that our results depend only very weakly on this choice. We also adopt this method to calculate the energy and angular momentum radiated by gravitational waves. Exceptionally, we adopt the previous method of direct time integration to estimate the orbital eccentricity in the inspiral phase because the fixed-frequency integration method may change the modulation in the gravitational waveform.

For comparisons between numerically calculated gravitational waveforms and those calculated in the PN approximations, we use the Taylor-T4 formula for two-point masses in circular orbits bbkmpsct2007 () with an additional contribution from the BH spin angular momentum santamariaetal2010 (). In this formula, the time evolution of the orbital angular velocity and orbital phase are computed using a nondimensional angular velocity by

 dXdt = 64νX55m0[1−(743336+114ν)X+(4π−11312χ+196νa)X3/2 (22) +(3410318144+5χ2+136612016ν+5918ν2)X2 −{(4159672+1898ν)π+(315711008−116524ν)χ+34χ3−(218631008ν−796ν2)a}X5/2 +{16447322263139708800−1712105γE+163π2−(56198689217728−45148π2)ν+541896ν2−56052592ν3 −856105ln(16X)−80π3χ+(641531008−45736ν)χ2+(203π−113536χ)νa}X3 −{(44154032−3586756048ν−914951512ν2)π+(252940727216−8458276048ν+41551864ν2)χ−12πχ2 +(150524+ν8)χ3−(158023954432−4515976048ν2+2045432ν3+1076νχ2)a}X7/2],
 dΘdt=X3/2m0, (23)

where , , and is the Euler constant. After and are obtained, we calculate the complex gravitational-wave amplitude of the mode and the spectrum up to the 3PN order using the formula shown in Refs. kidder2008 (); santamariaetal2010 (). Here, is

 h22 = −8√π5νm0De−2iΘX[1−(10742−5542ν)X+(2π−43χ+23νa)X3/2 (24) −(21731512+1069216ν−20471512ν2)X2−{(10721−3421ν)π+24iν}X5/2 +{27027409646800−856105γE+23π2−428105ln(16X)−(27818533264−4196π2)ν−202612772ν2 +11463599792ν3+428105iπ}X3],

where is the distance between the center of mass of the binary and an observer. Hereafter, we simply refer to this formula as the Taylor-T4 formula irrespective of the presence of the BH spin. Another way for deriving an approximate waveform is to employ an effective one-body approach (see Ref. damournagar () and references therein for reviews). In accompanied papers lksbf (), comparisons between numerical waveforms and those of the effective one body approach are extensively performed.

To estimate the mass of the remnant disk, we calculate the total rest mass outside the AH

 Mr>rAH≡∫r>rAHρ∗d3x, (25)

where is the radius of the AH as a function of the angular coordinates. We note that we systematically underestimated disk masses in our previous works performed with an old version of SACRA skyt2009 (); kst2010 (), because we evolved hydrodynamic variables and estimated disk masses only in the finer domains (described in Sec. III.3) of the size . Such a domain size is insufficient for the estimation of the disk mass if tidal disruption occurs at a distant orbit, especially for the case in which the NS radius is large ( km). In this study, we correct the treatment of hydrodynamics and the estimation of disk masses: We follow the hydrodynamics for a wide computational domain of the size . We still possibly underestimated disk masses because some of the material escapes from our computational domains and we cannot follow their return which would occur if they are bounded.

We determine key quantities of the remnant BH, i.e., the mass and nondimensional spin parameter , from the circumferential radius of the AH, assuming that the deviation from the Kerr spacetime is negligible in the vicinity of a BH horizon. We estimate the remnant BH mass, , from the circumferential radius of the AH along the equatorial plane divided by , i.e., , which gives the BH mass in the stationary vacuum BH spacetime. Similarly, the nondimensional spin parameter of the remnant BH, , is estimated from the ratio of the circumferential radius of the AH along the meridional plane to using the relation

 CpCe=√2^r+πE(a2f2^r+). (26)

This also holds for the stationary vacuum BH with the nondimensional spin parameter . Here, is a normalized radius of the horizon, and is an elliptic integral

 E(z)=∫π/20√1−zsin2θdθ. (27)

For comparison, the nondimensional spin parameter of the remnant BH is also estimated from and the irreducible mass of the remnant BH using the relation

 Mirr,f=Ce4√2π√1+√1−af2, (28)

which holds for the stationary vacuum BH. The spin parameter obtained using this relation is referred to as according to Ref. skyt2009 (). Finally, we also estimate from the values of the remnant BH computed using approximate conservation laws

 MBH,c ≡ M0−Mr>rAH−ΔE, (29) JBH,c ≡ J0−Jr>rAH−ΔJ, (30)

where the total angular momentum of the material located outside the AH, , is approximately defined by

 Jr>rAH≡∫r>rAHρ∗huφd3x. (31)

Here, we assume that the orbital angular momentum of the BH is negligible. The nondimensional spin parameter of the remnant BH is defined by , again according to Ref. skyt2009 ().

### iii.3 Setup of AMR grids

In SACRA, an AMR algorithm is implemented so that both the radii of compact objects in the near zone and the characteristic gravitational wavelengths in the wave zone can be covered with sufficient grid resolutions simultaneously. Our AMR grids consist of a number of computational domains, each of which has the uniform, vertex-centered Cartesian grids with grid points for with the equatorial plane symmetry at . We always choose for the best resolved runs in this work. We also perform simulations with and 42 for several arbitrary chosen models to check the convergence of the results and find approximately the same level of convergence as that found in the previous work (see the Appendix of Ref. kst2010 ()). In the Appendix A of this paper, we show the convergence of gravitational waveforms and the masses of the remnant disks. The AMR grids are classified into two categories: one is a coarser domain, which covers a wide region, including both the BH and NS, with its origin fixed at the approximate center of mass throughout the simulation. The other is a finer domain, two sets of which comove with compact objects and cover the region in the vicinity of these objects. We denote the edge length of the largest domain, the number of the coarser domains, and the number of the finer domains by , , and , respectively. Namely, the total number of the domains is . The grid spacing for each domain is , where is the depth of each domain.

Table 4 summarizes the parameters of the grid structure for our simulations. The structure of the AMR grids depends primarily on the mass ratio of the binary because the distances between two objects and the center of mass depend strongly on the mass ratio for our initial models. Specifically, we choose for all binaries with and , and 4. We choose for binaries with . For binaries with , we choose because we do not evaluate disk masses for them. In all the simulations, is chosen to be larger than or comparable to the gravitational wavelengths at an initial instant . One of the two finest regions covers the semimajor axis of the NS by –45 grid points. The other covers the coordinate radius of the AH typically by grid points, depending on the mass ratio and the BH spin. For the runs, the total memory required is about 11 G bytes. We perform numerical simulations with personal computers of 12 G bytes memory and of core-i7X processors with clock speeds of 3.2 or 3.33 GHz. We use 2–6 processors to perform one job with an OPEN-MP library. The typical computational time required to perform one simulation (for ms in physical time of coalescence for the case) is 4 weeks for the 6 processor case.

## Iv Numerical results

We present numerical results of our simulations, focusing, in particular, on their dependence on the BH spin and NS EOS. First, we review general merger processes in Sec. IV.1. Sections IV.2, IV.3, and IV.4 are devoted to the analysis of properties of the remnant disk and BH formed after the merger. Gravitational waveforms are shown in Sec. IV.5, their spectra in Sec. IV.6, and the energy and angular momentum radiated by gravitational waves in Sec. IV.7.

### iv.1 Overview of the merger process

Figure 1 plots the evolution of the coordinate separation defined by for models HB-Q2M135a5, HB-Q2M135, and HB-Q2M135a-5, for which takes the same values. Here, is the position of the maximum rest-mass density and is the location of the puncture, . Figure 1 shows that the numbers of orbits increases as the BH spin increases from retrograde to prograde elsb2009 (). Specifically, the number of orbit are , 5.5, and 4 for , 0, and , respectively. This difference comes primarily from the spin-orbit interaction between these two angular momenta kidder1995 (); in the PN approximation, a force proportional to the inner product of the orbital and spin angular momenta of two objects appears at 1.5PN order. Here, we do not have to consider the NS spin angular momentum in the assumption of the irrotational velocity field and, therefore, we only consider the interaction between the orbital and BH spin angular momenta throughout this paper. When these two angular momenta are parallel and the inner product is positive (), an additional repulsive force works between the BH and NS. This repulsive force reduces the orbital angular velocity because the centrifugal force associated with the orbital motion can be reduced, and hence, the luminosity of gravitational radiation, which is proportional to , is also reduced. This strong dependence of the luminosity on makes the approaching velocity smaller in the late inspiral phase, and, therefore, the number of orbits increases. Conversely, when these two angular momenta are antiparallel (), an additional attractive force increases the angular velocity and gravitational-wave luminosity in the late inspiral phase. In this case, the orbital separation decreases faster due to a larger approaching velocity, and the number of orbits becomes smaller as the retrograde BH spin increases. All these results agree qualitatively with those of Ref. elsb2009 ().

The fate of BH-NS binaries is classified into two categories. One is the case in which the NS is disrupted by the BH tidal field before the BH swallows the NS, and the other is the case in which the BH swallows the NS without tidal disruption. In this paper, we focus mainly on the former case. We plot snapshots of the rest-mass density profiles and the location of the AH on the equatorial plane at selected time slices for models HB-Q3M135a75, HB-Q3M135a5, and HB-Q3M135a-5 in Figs. 24, respectively. The NS is disrupted outside the ISCO in the cases (Figs. 2 and 3) and forms a one-armed spiral structure with a large angular momentum. The material in the inner part of the spiral arm gradually falls onto the BH due to angular momentum transport via hydrodynamic torque in the spiral arm. The material with a sufficiently large specific angular momentum escapes the capture by the BH and forms an accretion disk, which survives for a time much longer than the dynamical time scale a few ms. We note that the prompt infall of the one-armed spiral structure onto the BH occurs from a relatively narrower region for than for . The reason is that the inner edge of the spiral arm contacts the AH well before the arm becomes nearly axisymmetric due to a large radius of the AH and ISCO for . The infall of the disrupted material from a narrow region of the BH frequently occurs when the NS is tidally disrupted in a binary with a high mass ratio, whereas this is rare in a binary with a nonspinning BH because the NS is not disrupted in a high mass-ratio binary. This difference in the merger process is well-reflected in gravitational waveforms (see Sec. IV.5). By contrast, the NS is swallowed by the BH without tidal disruption, and essentially no material is left outside the ISCO for model HB-Q3M135a-5 (Fig. 4).

Note that the feature of the NS tidal disruption appears very weakly not only for model HB-Q3M135a-5 but also for model HB-Q3M135 () because the mass ratio is so high that the tidal effect is less important for the nonspinning BH with the typical NS radius –12 km. The enhancement of the tidal effect by a prograde BH spin results primarily from the decrease of the BH ISCO radius bpt1972 (). In the Boyer-Lindquist coordinates, a Kerr BH has an ISCO with a smaller radius than a Schwarzschild BH by a factor of , depending on for a prograde orbit: The ISCO radius approximately halves when the BH spin increases from to 0.75. On the other hand, the orbital separation at the onset of mass shedding depends only weakly on the BH spin in the Boyer-Lindquist coordinates fishbone1973 (); marck1983 (); ism2005 (). This decrease of the ISCO radius enhances the possibility for the disrupted material to escape capture by the BH and to form a more massive remnant disk than in the nonspinning BH case. The retrograde BH spin plays an opposite role; the ISCO radius of the Kerr BH increases by a factor of 1–1.5 for a retrograde orbit, and hence, the tidal effect is less important in the merger process.

Before closing this subsection, we estimate the degree of (undesired) orbital eccentricity in our simulations to assess the circularity of the orbital motion. For this purpose, we compute the evolution of the gauge-invariant orbital angular velocity , which is defined from the mode of by

 Ω(t)=12|Ψ4(l=m=2)||∫Ψ4(l=m=2)dt|. (32)

The evolution of the orbital angular velocity in our simulation agrees with that derived from the Taylor-T4 formula in the inspiral phase within a small modulation, typically , which is equivalent to the orbital eccentricity of . This amount of orbital eccentricity is as small as that observed in the nonspinning BH case with a low mass ratio kst2010 ().

### iv.2 Global properties of the disk

The mass of the remnant disk reflects the significance of NS tidal disruption in a clear way because the disk formation is a result of tidal disruption. A massive disk is formed if tidal disruption of the NS occurs far outside the ISCO. If the mass shedding starts in the vicinity of or inside the ISCO, only a small portion of the mass is left outside the AH. The material is not left outside the AH when the mass shedding does not occur before the BH swallows the NS, and the merger of a BH-NS binary may be indistinguishable from that of a BH-BH binary except for very small tidal corrections to the inspiral. Thus, the mass of a remnant disk is a reliable indicator of the degree of tidal disruption.

Figure 5 plots the time evolution of the rest mass located outside the AH, , for and 3 with different nondimensional BH spin parameters , 0.5, 0, and . In both plots, and HB EOS are adopted. We note that the results revised from the previous work kst2010 () are plotted for . The dependence of on for HB EOS found here is similar to those for other EOSs. We set the time origin to be an approximate merger time . These plots indicate that the mass of the material left outside the AH relaxes to a quasisteady value for –4 ms, and the relaxed value increases monotonically as the BH spin increases from retrograde to prograde. This is consistent with the decrease of the BH ISCO radius with the increase of its spin, as described in Sec. IV.1. In particular, the remnant disk mass at ms after the merger is for all the EOSs with and , as shown in Table 5, and for , irrespective of the EOS. The formation of such a massive disk may be encouraging for the BH-NS binary merger hypothesis of a short-hard GRB. For the cases, by contrast, massive accretion disks of are not expected to be formed as merger remnants even for unless the EOS is extremely stiff (the NS radius is 15 km). This fact indicates that the retrograde BH spin is unfavorable for producing a central engine of a short-hard GRB.

The prograde BH spin enhances the disk formation dramatically for a BH-NS binary with a high mass ratio, for which the disk mass is very low when the BH is nonspinning. We plot the time evolution of for and 5 with different EOSs in Fig. 6. In both plots, and are adopted. Figure 6 clearly shows that a massive accretion disk is formed for and 5 if the BH has a prograde spin of . Namely, the formation of a massive accretion disk is universal for the merger of a BH-NS binary with a mass ratio of as far as and (equivalently, ). Note that a heavy BH of is predicted to be realistic as an astrophysical consequence of the stellar evolution with solar metallicity mcclintockremillard () (see, e.g., Ref. bbfrvvh2010 () for a population synthesis study) and hence as a possible progenitor of the short-hard GRB.

For more quantitative discussion, we plot the disk mass estimated at ms after the merger for all the models with and for models with as a function of the NS compactness, , in Fig. 7. Numerical values of are shown in Table 5, as well as other quantities associated with the merger remnants. For any fixed value of , a negative correlation between and is found to hold in Fig. 7. This correlation indicates that the NS with a larger compactness is less subject to tidal deformation and disruption than the NS with a smaller compactness for any fixed value of . This correlation is expected from the nature of a tidal force as a finite-size effect, as found in the study of nonspinning BH-NS binaries kst2010 (). On the other hand, Fig. 7 again shows that the prograde BH spin increases the disk mass for any fixed value of . A remarkable fact is that the disk mass does not decrease steeply to a value of as the compactness increases for binaries with . We expect that the coalescence of a BH-NS binary with may always produce a remnant disk of within a plausible range of the NS compactness, , although it is possible only if for and for or .

The dependence of the disk mass on the NS compactness is different for different values of the mass ratio. We plot in Fig. 8 the disk mass as a function of the NS compactness as in Fig. 7, but for and 0.5. This figure shows that the disk mass depends more strongly on when the mass ratio, , is larger. The disk mass is larger for smaller values of when the EOS is soft and , except for HB-Q2M135a75 and HB-Q3M135a75, for which the disk masses depend only weakly on . This dependence on is expected from the comparison between the mass-shedding radius, , and the ISCO radius, ,

 rshedrISCO∝C−1Q−2/3, (33)

where we assume Newtonian gravity for simplicity. This relation states that a larger amount of mass can escape the capture by the BH and can form an accretion disk when is small because the mass shedding sets in at relatively more distant orbit. However, the disk mass may be larger for larger values of when the EOS is stiff as for and . This should be ascribed to the redistribution process of the specific angular momentum of the NS to the disrupted material and to subsequent behavior of the material (such as collision of the fluid elements in spiral arms). This feature suggests that a binary with a larger value of , say , possibly form a massive remnant disk of and could be a progenitor of a short-hard GRB if the EOS is stiff and the BH has a large spin .

To clarify the dependence of the disk mass on the BH spin, we plot the disk mass as a function of in Fig. 9. The EOS (and, equivalently, ) is the same for each plot. Again, we find a monotonic and steep increase of the disk mass as the increase of for the fixed EOS and mass ratio. The enhancement of the disk mass by a prograde spin is more dramatic for the compact NS (for the soft EOS). For example, the difference in the disk mass between the cases of and is only by a factor of when and 2H EOS is adopted. This low amplification is natural because tidal disruption of a large NS occurs at an orbit far enough from the ISCO for a substantial amount of the disrupted material to escape the capture by the BH irrespective of and because at such a large orbital separation the spin-orbit coupling effect is relatively weak. On the other hand, a few-orders-of-magnitude amplification of the disk mass is seen when and HB EOS is adopted.

Finally, we comment on a possible unbound outflow. To estimate the rest mass of unbound material, we compute

 Mub≡∫r>rAHρ∗H(−ut−1)d3x, (34)

where is a step function. Here, the material with should be considered to have an unbound orbit. We find that can be larger than at ms after the merger for the stiff EOS like 2H and H, and . However, does not approach a constant value and rather continues to decrease. Therefore, it is unclear whether estimated at 10 ms after the merger can really become unbound or not, and we do not show the precise values of . When the EOS is not stiff, is negligible within the accuracy of our simulations.