Black-hole–neutron-star mergers at realistic mass ratios: Equation of state and spin orientation effects

Black-hole–neutron-star mergers at realistic mass ratios: Equation of state and spin orientation effects

Abstract

Black-hole–neutron-star mergers resulting in the disruption of the neutron star and the formation of an accretion disk and/or the ejection of unbound material are prime candidates for the joint detection of gravitational-wave and electromagnetic signals when the next generation of gravitational-wave detectors comes online. However, the disruption of the neutron star and the properties of the post-merger remnant are very sensitive to the parameters of the binary (mass ratio, black hole spin, neutron star radius). In this paper, we study the impact of the radius of the neutron star and the alignment of the black hole spin on black-hole–neutron-star mergers within the range of mass ratio currently deemed most likely for field binaries () and for black hole spins large enough for the neutron star to disrupt (). We find that: (i) In this regime, the merger is particularly sensitive to the radius of the neutron star, with remnant masses varying from to for changes of only in the NS radius; (ii) of unbound material can be ejected with kinetic energy , a significant increase compared to low mass ratio, low spin binaries. This ejecta could power detectable post-merger optical and radio afterglows. (iii) Only a small fraction of the Advanced LIGO events in this parameter range have gravitational-wave signals which could offer constraints on the equation of state of the neutron star (at best of the events for a single detector at design sensitivity). (iv) A misaligned black hole spin works against disk formation, with less neutron star material remaining outside of the black hole after merger, and a larger fraction of that material remaining in the tidal tail instead of the forming accretion disk. (v) Large kicks can be given to the final black hole as a result of a precessing BHNS merger, when the disruption of the neutron star occurs just outside or within the innermost stable spherical orbit.

pacs:
04.25.dg, 04.30.-w, 04.40.Dg, 47.75.+f

I Introduction

The next generation of ground-based gravitational-wave detectors (Advanced LIGO, Advanced Virgo and KAGRA Harry (for the LIGO Scientific Collaboration) (2010); Advanced Virgo Baseline Design, The Virgo Collaboration, note VIR–027A–09, May 16 2009 (); Somiya for the KAGRA Collaboration (2012)) will progressively begin taking data over the next decade, opening an entirely new way to observe the universe. One of the sources of gravitational waves that they will detect are compact binary coalescences: binary black holes (BBH), black-hole–neutron-star (BHNS) and binary neutron star (BNS) systems Abadie et al. (2010). In the presence of a neutron star, these gravitational-wave signals could be accompanied by electromagnetic emissions, which would provide better sky localization and additional information about the characteristics and the environment of the binary. The most energetic, and most often discussed potential counterparts are short gamma-ray bursts (SGRBs, see e.g. Lee and Ramirez-Ruiz (2007)), while other possibilities include the x-ray and optical afterglows of SGRBs, optical transients due to the radioactive decay of neutron-rich unbound material, and radio emission from that ejecta as it decelerates in the interstellar medium (see Metzger and Berger (2012) for more details on EM signals emitted by compact binary mergers, and Nissanke et al. (2012) for their detectability). The ejection of a small amount of material at ultrarelativistic speeds as a result of a shock in the region in which two neutron stars first get into contact was also recently proposed as a potential outcome of BNS mergers Kyutoku et al. (2012). Finally, pre-merger electromagnetic transients are also a possibility, and could for example be due to the breaking of the neutron star crust Tsang et al. (2012).

Which of these effects occur in practice depends strongly on the parameters of the binary. The exact conditions leading to the emission of a SGRB are not known, and will depend on the physical process responsible for these bursts — most likely either the extraction of black hole rotational energy by the magnetic field of the accretion disk (Blandford-Znajek effect Blandford and Znajek (1977)), magnetically-driven outflows in the accretion disk (Blandford-Payne effect Blandford and Payne (1982)), or the production of ultrarelativistic pairs from the annihilation of pairs, themselves produced by a hot accretion disk surrounding the remnant black hole Lee et al. (2005); Lee and Ramirez-Ruiz (2007). In all cases, the presence of both an accretion disk and a baryon-poor region in which a relativistic jet can be produced appears to be a prerequisite. The amount of matter required in the disk depends on the efficiency of the jet production mechanism and the energy of the burst, with estimates spanning multiple orders of magnitudes (Lee et al. (2005); Giacomazzo et al. (2012). The formation of an accretion disk is a natural result of the merger of BNS systems. For BHNS binaries, however, many initial conditions lead to the direct plunge of the neutron star into the black hole, before the tidal field of the hole can disrupt the neutron star and cause the formation of an accretion disk. Relatively massive black holes () are favored both by current population synthesis models Belczynski et al. (2008, 2010) and by the distribution of black hole masses measured in galactic X-ray binaries Özel et al. (2010). In that regime, tidal disruption only occurs for the most rapidly spinning black holes Foucart et al. (2012).

The ejection of enough unbound material to power detectable electromagnetic transients is not a certainty either. General relativistic simulations of BNS systems have found that only a small amount of mass is unbound by the merger (Hotokezaka et al. (2012), while Newtonian smoothed particle hydrodynamics (SPH) simulations have more optimistic predictions (Korobkin et al. (2012). BHNS systems have more asymmetric mass ratios, and are thus generally more favorable for the ejection of neutron-rich material during the disruption of the neutron star. However, this requires the neutron star to be disrupted in the first place. Numerical simulations have shown that this will only be the case if the mass of the black hole is lower than what population synthesis models predict, or if the spin of the black hole is high. The amount of material ejected in such mergers remains very uncertain. SPH results predict ejected masses of up to  Rosswog (2005). General relativistic simulations for low mass, low spin black holes find little ejected material (Duez et al. (2010), but estimates for rapidly spinning black holes have not been offered yet. We show in this paper that more massive outflows () are likely for high black hole spins. Finally, matter could be ejected due to magnetically-driven Shibata et al. (2011) or neutrino-driven Dessart et al. (2009) winds in the disk.

An important limitation of existing general relativistic simulations of BHNS mergers is the lack of coverage of the range of black hole masses deemed most likely astrophysically. The only simulations considering mass ratios have shown that, even for very large neutron stars () and aligned black hole spins, dimensionless black hole spins are required for tidal disruption to occur Foucart et al. (2012). Simulations for more symmetric mass ratios have however already shown that smaller neutron star radii Duez et al. (2010); Kyutoku et al. (2011) or misaligned black hole spins Foucart et al. (2011) are likely to make tidal disruption harder. In this paper, we begin a more quantitative exploration of these effects. We consider BHNS binaries at mass ratio and with large black hole spins , and vary the radius of the neutron star and the orientation of the black hole spin. We currently limit ourselves to the most basic physical effects which will influence the dynamics of the merger: we treat gravity in a fully general relativistic framework, but use a simple -law equation of state to describe the neutron star matter, and ignore the effects of magnetic fields or neutrino radiation. More realistic equations of state are likely to influence tidal disruption, but to first order the compactness of the neutron star is expected to capture the main physical effects during merger Duez et al. (2010). Magnetic fields are important for the evolution of any post-merger accretion disk, as these disks are unstable to the magnetorotational instability (MRI), but very large internal fields are necessary for magnetic effects to influence the disruption of the neutron star Chawla et al. (2010); Etienne et al. (2012a). Finally, while neutrinos are the main source of cooling in the disk, and are thus crucial to their evolution over a cooling timescale ( for the relatively dense disks considered by Lee et al. Lee et al. (2005), and could be significantly smaller for the lower density disks observed at the end of our simulations), neutrino emission will be negligible before merger (neutron stars are expected to be extremely cool, with ). This was confirmed numerically in the case of BNS Kiuchi et al. (2012). Our simulations will thus capture the physical effects which are important for the evolution of BHNS systems in the last few orbits before their merger, and during the merger itself. They are however not suitable for the long term evolution of the post-merger remnant.

Gravitational waveforms from BHNS and BNS mergers are of particular interest for the constraints that they might offer on the unknown equation of state of the neutron star. Numerical simulations indicate that from the last few orbits of a BNS merger occurring at , constraints of could be obtained on the radius of the neutron star Read et al. (2009). Damour et al. Damour et al. (2012) have shown using Effective One Body waveforms that equations of state effects in the late inspiral could be measured for events of moderate signal-to-noise ratio (). Accurate numerical simulations are however necessary to calibrate such models at high frequency. Numerical results by Bernuzzi et al. Bernuzzi et al. (2012) have confirmed the predictions of Damour et al. Damour et al. (2012) regarding the detectability of these equations of state effects, but existing simulations are not accurate enough to model the waveform with the accuracy required to take full advantage of all of the information that will be available in waveforms detected by Advanced LIGO. For BHNS mergers of nonspinning black holes, tidal effects during the inspiral are too small to be detected directly by Advanced LIGO Hinderer et al. (2010); Damour et al. (2012). But the cutoff in the gravitational-wave spectrum occurring when the neutron star is disrupted by the tidal field of the black hole can be. Semi-analytical models have been developed to attempt to extract that information Vallisneri (2000); Ferrari et al. (2010). Numerical simulations mapping the cutoff frequency across the relevant parameter space are however necessary to better calibrate them.

At low mass ratios () and for nonspinning black holes, Lackey et al. Lackey et al. (2012) have shown from numerical simulations that the combined effects of the tidal interactions during the inspiral and of the high-frequency cutoff of the signal would allow Advanced LIGO to detect variations of in the radius of the neutron star for a favorable event at . This is thus slightly inferior to what can be done for binary neutron star systems located at the same distance from the observer (and, within a fixed volume, we expect many more BNS mergers than BHNS mergers). At higher mass ratio, tidal effects during the inspiral further decrease Hinderer et al. (2010). On the other hand, due to the higher total mass of the system, the amplitude of the signal will be larger. Tidal disruption will also occur at lower frequency, and thus in a more favorable region of the LIGO noise curve. The spin of the black hole can also affect how much tidal distortion occurs during the inspiral (as the neutron star can get closer to a spinning BH before reaching its ISCO, tidal effects can be stronger). How these competing effects will affect our ability to measure the properties of neutron stars in BHNS binaries, or even to differentiate BHNS systems from BBH binaries, remains an important open question.

In this paper, we study the influence of the radius of the neutron star and of the orientation of the black hole spin on the dynamics of the merger of BHNS binaries around the black hole mass currently favored by population synthesis models, focusing on tidal effects during the inspiral, on the initial characteristics of the post-merger remnant (accretion disk and tidal tail formation), and on the properties of the emitted gravitational-wave signal — and in particular the effects of the neutron star radius on that signal, and the conditions under which such effects might be detected by the next generation of gravitational-wave experiments. We stop the simulations after merger, as neglected microphysical effects are likely to become important at later times. We will begin by describing briefly the numerical setup used for our simulations, as well as modifications to the code since the publication of Foucart et al. (2012) (Sec. II). We will then detail the initial configurations evolved (Sec. III), and estimate the accuracy of the results (Sec. IV). Finally, the main physical results are presented in Sec. V.

Ii Numerical Setup

The numerical simulations presented here are performed with the SpEC code SpEC Collaboration (), which evolves Einstein’s equations of general relativity coupled to the relativistic hydrodynamics equations (see Appendix A for details). Einstein’s equations are solved using pseudospectral methods, in the generalized harmonics formulation Lindblom et al. (2006), and excising the black hole interior. The numerical grid on which Einstein’s equations are solved consists at first of 8 spherical shells immediately around the black hole, 8 spherical shells and 1 inner ball in the region close to the neutron star, 24 spherical shells covering the far-field region, and a set of 13 distorted cylindrical shells and filled cylinders connecting them (see Fig. 1). All subdomains are touching but not overlapping. The low, medium and high resolution runs correspond to a total number of points .

Figure 1: Numerical grid before disruption of the neutron star, below the equatorial plane of the binary. The black hole is the excised region on the right (black sphere), while we superpose a linear color scale for the baryon density . The subdomains on the outer edge of the plot connect to spherical shells covering the wave region.

Once the neutron star disrupts, we cannot take advantage of an approximate spherical symmetry around the neutron star, and have to modify the pseudospectral grid. The shells around the black hole are replaced by a set of 264 distorted cubes, in order to allow for high angular resolution as the neutron star falls into the hole. The wave zone is still covered by 24 spherical shells, while the region around the neutron star and the near field region are covered by non-overlapping distorted cubes (see Fig. 2). The resolution is chosen adaptively, by requiring that the relative truncation error for each set of basis functions (measured from the coefficients of the spectral expansion of the evolved variables) is for our 3 resolutions. The actual number of grid points thus vary during the merger, peaking as the neutron star accretes onto the black hole. At medium resolution, we have .

Figure 2: Same as Fig. 1, but during the disruption of the neutron star.

The relativistic hydrodynamics equations are solved on a separate finite difference grid Duez et al. (2008). The grid covers only the region in which matter is present, and expands/contracts at discrete times as needed. Before disruption, we use points for the 3 different resolutions. During merger, we use instead points. For configurations in which the black hole spin is aligned with the orbital angular momentum, we only evolve the region above the orbital plane and impose symmetry conditions on that plane (the number of grid points in the direction orthogonal to the orbital plane is then divided by 2).

Compared to the simulations performed in Foucart et al. (2012), we also fixed an error in the algorithm responsible for communicating source terms between the two numerical grids, which could cause the time stepping algorithm to be effectively of lower order than expected. The correction leads to reduced errors in the observed trajectories and in the phase of the gravitational waveforms (see Sec. IV.2).

Iii Initial Configurations

The initial data for these simulations is constructed as described in Foucart et al.(2008) Foucart et al. (2008). The constraints that Einstein’s equations impose on the initial variables are solved in the extended conformal thin sandwich approximation York (1999), under the assumption that the system is in equilibrium in a frame rotating at angular velocity and contracting with radial velocity . The values of and determine the eccentricity of the system, and are chosen iteratively in order to minimize that eccentricity Pfeiffer et al. (2007). We go through two iterations of that procedure, starting from quasi-circular orbits (i.e. initial data with and chosen so that the initial motion of the neutron star center has no radial component). The residual eccentricities at the end of the iterative procedure are . The free variables in the initial data (conformal metric, extrinsic curvature) are the weighted superposition of an isolated black hole in Kerr-Schild coordinates and of an isolated neutron star in isotropic coordinates, following the method developed by Lovelace et al.(2008) Lovelace et al. (2008) for binary black holes. A more detailed description of the modifications required to apply this method to black-hole–neutron-star systems is given in Foucart et al. (2008) Foucart et al. (2008).

Name
R12i0 0.170 12.2 km 0.0413 20.5 0.004
R13i0 0.156 13.3 km 0.0413 20.3 0.003
R14i0 0.144 14.4 km 0.0413 19.7 0.003
R14i20 0.144 14.4 km 0.0412 0.003
R14i40 0.144 14.4 km 0.0413 0.004
R14i60 0.144 14.4 km 0.0415 0.002
Table 1: Initial configurations studied. All binaries have a mass ratio of 1:7, with the black hole dimensionless spin magnitude being . is the angle between the rotation axis of the black hole and the initial orbital angular momentum of the binary, is the compactness of the neutron star, and the Schwarzschild radius of that neutron star assuming an ADM mass in isolation of . The orbital parameters are the eccentricity , the initial orbital angular velocity times the total mass of the system, , and the number of gravitational-wave cycles before the peak of the wave amplitude, (approximate numbers given for the precessing systems, in which mode mixing makes this variable ill-defined).

Two series of initial configurations are considered in this paper (see Table 1), chosen in order to study separately the effects of the radius of the neutron star and of the orientation of the black hole spin on BHNS mergers at realistic mass ratios. All configurations consider a black hole of mass and spin , where is the Christodolou mass of the black hole, its angular momentum, and the ADM mass in isolation of a neutron star with the same baryon mass as the neutron star in the binary. The neutron star is initially nonspinning. As our initial data does not exactly represent a BHNS binary in quasi-equilibrium, small transients are always observed at the beginning of the simulations. After those transients, the mass and spin of the black hole are slightly modified (by ). All simulations describe the neutron star matter as an ideal fluid with stress-energy tensor , and use a -law equation of state of index :

(1)
(2)

where is the baryon density of the neutron star material, its pressure, its internal energy, its 4-velocity, a free constant and a variable related to the temperature of the fluid ( is the thermal pressure in the fluid). To obtain the physical temperature of the fluid from the variable , we assume that the thermal pressure is the composition of an ideal gas component and a black body component:

(3)

where is the nucleon mass, the Boltzmann constant, and a function of reflecting the fraction of relativistic particles in the gas (see Shibata and Taniguchi (2008); Etienne et al. (2009) for details).

In the first group of simulations, we consider black hole spins aligned with the orbital angular momentum of the system, and modify the radius of the neutron star between and (for ). This is done by modifying the value of the free parameter in the equation of state of the fluid. We only consider this simple variation of the equation of state as we know that, to first order, the radius of the neutron star is the most important contribution to the dependence of BHNS mergers on the equation of state of the fluid Duez et al. (2010), while the tidal deformability determines tidal effects during the inspiral Hinderer et al. (2010) ( being the tidal Love number of the neutron star). With respect to these parameters, the configurations considered here are within the range currently allowed for real neutron stars. They do however fail to reproduce other properties of neutron stars which are not as relevant to this study. For example, all three equations of state have a maximum mass smaller than (which is more important for studies of stellar collapse and of the evolution of hypermassive neutron stars than for tidal disruption and BHNS mergers, in which the maximum density of the fluid only decreases over the course of the evolution), they have a very simplified temperature dependence, and they do not describe the composition of the fluid. The largest neutron star (case R14i0) is identical to the configuration studied in Foucart et al.(2012) Foucart et al. (2012), except that the initial separation is larger than in our previous work. As we will see, these parameters study most of the transition between BHNS mergers resulting in the formation of massive disks and those having nearly no material left out of the black hole a few milliseconds after merger.

As we are considering -law equations of state, it is worth mentioning that any of these simulations actually represent a continuum of systems, as they are invariant through the rescaling

(4)
(5)
(6)

where is the mass scale of the binary, the distance scale, the time scale, and an arbitrary positive constant. A more useful description of the initial conditions would thus use quantities which are also invariant under this rescaling, i.e. the mass ratio , the compactness of the neutron star , and the dimensionless time (in units in which ). Realistic neutron stars of probably have radii in the range  Hebeler et al. (2010), with the most likely values being  Steiner et al. (2010). The values of considered here are thus more likely to be found in neutron stars of ADM mass around or slightly below , while probably unrealistically low for very massive neutron stars ().

The second set of simulations considers variations of the orientation of the black hole spin while maintaining the equation of state fixed (using the larger neutron star with ). In terms of disk formation, this choice of equation of state is clearly optimistic, although not unrealistic, and thus provides an upper bound on the mass remaining outside of the black hole after merger. Moderate misalignments of the spin of the black hole with respect to the angular momentum of the binaries are an expected consequence of the kick that the supernova explosion is likely to impart to the forming neutron star. The actual distribution of the misalignment angle between the spin of the black hole and the orbital angular momentum of the binary is currently unknown, although should probably be favored Belczynski et al. (2008). Here we vary between and , as such misalignments are both physically realistic and covering the range of parameters over which the properties of the final remnant vary significantly (at least for BHNS systems of mass ratio with black hole spins ). Misalignments are also often quoted as the angle between the black hole spin and the total angular momentum of the system. For the systems considered here, , we have .

All simulations begin at a coordinate separation , where is the total ADM mass of the system at infinite separation. This corresponds to an initial orbital angular velocity , or an initial gravitational-wave frequency . Over the course of the simulation, the binary will go through orbits before merging.

Iv Accuracy

The combination of spectral and finite difference methods used in our simulations can make it difficult to obtain strict error estimates: spectral methods are exponentially convergent in regions in which all variables are smooth, but only show polynomial convergence in the presence of discontinuities (such as at the surface of the neutron star or at a shock front). The finite difference methods used to evolve the equations of relativistic hydrodynamics are second order in smooth regions, and first order at the location of shocks. As the region in which we get first order convergence should be of measure zero, we expect at least second-order convergence as we increase the resolution of the finite difference grid. In practice, we generally observe much faster convergence between the 3 resolutions considered here, particularly for quantities evolved on the spectral grid (e.g. trajectories, gravitational-wave signal,…). A conservative estimate of our error would thus be to assume second order convergence between our medium and high resolutions - the actual error being somewhere between that value and the optimistic estimate obtained by simply looking at the difference between those two simulations. The ratio between these pessimistic and optimistic error estimates is for the simulations presented here.

iv.1 Final Remnant

Of the characteristics of the final remnant listed in Table 2, the parameters of the black hole (mass and spin) are the most accurate, with relative errors (i.e. differences of between the medium and high resolutions). Global mass measurements (disk mass, tidal tail mass) are already less accurate, with ( difference measured in the final remnant mass of the medium and high resolution runs for configuration R13i0). Finally, the maximum density within the disk is only order-of-magnitude accurate: the distribution of matter within the disk remains fairly asymmetric and time-dependent at the end of the simulation, and variations of within should still be expected. Measurements of the mass of unbound material and the properties of this ejecta have similar errors, and are discussed in more details in Sec. V.1.3.

iv.2 Waveform Accuracy

Figure 3: Phase error in the dominant (2,2) mode of the gravitational-wave signal for simulation R13i0. We show the phase difference between our standard and high resolutions, both with (dashed red line) and without (solid black line) aligning the waves in phase and time. The dash-dotted green curve shows an estimate of the error in the extrapolation method (obtained by comparing extrapolation using different polynomial orders).

The phase accuracy of the gravitational waveforms in the non-precessing simulations presented here is about a factor of 2 better than in our last set of simulations Foucart et al. (2012), even though the evolutions are orbits longer. This is most likely due to the correction of an error which effectively decreased the order of the time stepping method used in our simulations. Fig. 3 shows the phase difference between the medium and high resolution of the inspiral of simulation R13i0, both without any matching (i.e. by directly computing the phase difference between the output of the 2 resolutions), and after matching the waveforms over one period of the radial oscillation of the orbit, choosing a time and phase shifts minimizing the difference between the two waveforms in the matching region. The first method is the most direct assessment of the effect of numerical errors on the phase of the gravitational-wave signal, which are here of the order of a few tenths of a radian during the inspiral. The second method is more useful when comparing waveforms obtained in simulations starting from different initial conditions, and shows how different the waveforms would look to a gravitational-wave detector. Fig. 3 shows that, for matched waveforms, differences of the order of a few percents of a radian or less cannot be resolved by our numerical simulations.

When considering waveform accuracy, numerical errors due to the discretization of the evolution equations are however only part of the problem. Another potential source of error comes from extracting gravitational waves at finite radii, and then using polynomial extrapolation to obtain the waveform at null infinity Boyle and Mroué (2009). An estimate of the error due to this process can be obtained by comparing the waveforms obtained using different polynomial orders for the extrapolation. Fig. 3 shows that this error is . Phase differences of the same order can also be due to the eccentricity of the binaries, at least for the eccentricities considered here. This can easily be seen from the oscillations in the phase difference between different configurations shown in Fig. 11 (the oscillations in the phase difference between simulations R13i0 and R14i0 are smaller that those between R12i0 and R14i0 because the radial oscillations of the first two cases happen to be nearly in phase at the beginning of the simulation).

Adding all these sources of error, we can thus estimate that differences between numerical waveforms are large enough to be measured at our current accuracy only if, for matched waveforms, we have during inspiral.

V Results

v.1 Non-Precessing Binaries

Inspiral : Tidal Effects

Before the disruption of the neutron star, the main differences between a BHNS inspiral and a BBH inspiral are due to the finite size of the neutron star, and its distortion under the influence of the tidal field of the black hole. The tidal distortion of the neutron star, and in particular its effect on the gravitational-wave signal, has already been studied in the Post-Newtonian(PN) framework. During the early inspiral, Hinderer et al. Hinderer et al. (2010) found that for BNS systems the tidal effects would only be detectable by Advanced LIGO for the most favorable configurations (i.e. the largest neutron stars, see also Pannarale et al. (2011) for similar results considering the tidal effects up to the disruption of the neutron star). Over the last few orbits, Damour et al. Damour et al. (2012) find that tidal parameters would be detectable for BNS mergers of moderate signal-to-noise ratio (). But as these effects are significantly smaller for more asymmetric mass ratios, the detection of tidal effects through gravitational waves is much more difficult for BHNS systems. A more detailed discussion of the detectability of the neutron star equation of state in our mergers is offered in Sec. V.1.5 - but from Fig. 11 alone, where we show the phase difference between our 3 non-precessing simulations, it is easy to see that up to 4 gravitational-wave cycles before the peak of the gravitational-wave signal () the difference between these cases is not resolved numerically. This could however be due either to the fact that tidal effects on the waveform are extremely small, or to a failure of the simulations to capture the tidal distortion of the neutron star properly. Accordingly, we need to test that the neutron star is tidally distorted during inspiral, and that these tidal effects scale as expected. We compute the quadrupole moments of the neutron star

(7)

(where , , is the spatial component of the 4-velocity and is the volume element, so that ), and assume that they are due to first order to the composition of the tidal distortion of the neutron star and of a coordinate boost, acting along orthogonal directions. is then diagonal in the coordinate frame , where are two orthogonal unit vectors in the equatorial plane of the binary and is a unit vector in the direction of the orbital angular momentum (by symmetry, ). The orientation of in the equatorial plane is a priori unknown, and practically determined by solving for the rotation matrix diagonalizing . To first order, we assume that the tidal distortion causes the neutron star to stretch along the direction and contract along and , while the boost causes a contraction along and a stretch along and , i.e.

(8)
(9)
(10)

where is the magnitude of the tidal distortion, and the amplitude of the boost distortion (by construction, . So we can solve any 2 of the 3 equations above for the unknowns and , and the third will automatically be satisfied). From this decomposition, we can also retrieve the lag angle between the tidal bulge and the line connecting the center of the hole and the center of the neutron star (i.e. the angle between and the line connecting the two centers).

These quantities are clearly dependent on the coordinate system chosen. We cannot entirely remove that dependence, but can at least get a reasonable normalization for the quadrupole moments from the quantity

(11)

For all simulations, we find similar boost components : is a function of the binary separation, but does not depend on the equation of state considered (for a true Lorentz boost of a spherical star, we should get ). The lag angle is fairly constant too, with for separations . The tidal component, on the other hand, varies strongly with both the binary separation and the equation of state. In Fig. 4, we show for the three different equations of state considered here, and in the range of binary separations . The tidal distortion goes from being of the same order as the boost effect at to a factor of larger at . Not surprisingly, the larger neutron star is significantly more distorted.

Figure 4: Top: Tidal quadrupole Q normalized by . Bottom: Normalized quadrupole (see eq. 13).

To leading order, we expect the tidal distortion and the tidal field of the black hole () to be related by

(12)

where is the tidal Love number of the neutron star (for polytropes, which at are equivalent to the -law equation of state used in our simulations, was computed by Hinderer Hinderer (2008)). To verify that the tidal effects scale as expected, we thus compute the normalized tidal parameter :

(13)

The superscript refers to values for simulation R14i0, and is the initial separation of the binary. The normalization is the value of for simulation R14i0 at the separation at which we are measuring . The numerical factor is computed in the limit (i.e. with computed for an isolated neutron star), so that . The scalings of , and are chosen so that, as long as the tidal distortion of the neutron star follows the theoretical predictions, we will measure .

In practice, our ability to measure accurately is limited by the fact that the boost and normalization only approximately model the distortion of the NS due to coordinate effects (i.e. the boost, but also any other gauge effect due to the coordinate choices made in the simulation). Measurements of are thus unreliable for . Figure 4 shows that for , the scatter in the measurement of is , while for , it decreases below . All measurements of are compatible with within that scatter, thus showing that the tidal distortion of the neutron star follows approximately the predictions of Ref. Hinderer (2008), even at close separations.

From these computations, we can thus conclude that the tidal distortion of the neutron star during the late inspiral is resolved by our numerical simulations, and scales with the binary separation and the equation of state of the neutron star in the manner expected from theoretical calculations, at least within .

Merger and Disk Formation

Name []
R12i0 0.10 0.06 2 0.923 0.960
R13i0 0.20 0.11 3 0.919 0.950
R14i0 0.30 0.16 21 0.910 0.935
R14i20 0.28 0.15 17 0.909 0.939
R14i40 0.15 0.10 3 0.898 0.959
R14i60 0.03 0.03 0.4 0.862 0.978
Table 2: Properties of the final remnant. is the baryon mass remaining outside of the black hole 5ms after merger. is the baryon mass located at a coordinate radius greater than at the same time. is the maximum density in the disk, the dimensionless spin of the black hole at the end of the simulation, the final Christodolou mass of the black hole and the ADM mass of the system at infinite separation.

An important question when considering BHNS mergers is the form of the post-merger remnant. To first order, this depends on whether the neutron star is disrupted before reaching the innermost stable circular orbit of the black hole or not. In the first case, a large amount of matter can remain outside of the hole after merger in the form of an accretion disk and a tidal tail. In the second case, no matter will remain. For a BHNS merger to be the progenitor of a short gamma-ray burst, the creation of an accretion disk is necessary. Accordingly, SGRBs are only possible if the neutron star disrupts. Stellar disruption is facilitated by low black hole masses, high black hole spins and large neutron stars (see Foucart (2012) for a simple fit to the results of previous numerical simulations). At low mass ratios (), a moderately spinning black hole is generally sufficient to provide disks of . For the more massive black holes considered here, however, this is no longer the case. We have already shown that for spins disk formation is unlikely even for large neutron stars (). The simulations presented here begin to explore how smaller neutron stars fare.

Figure 5: Merger for the non precessing cases R14i0 (left) and R12i0 (right). The top panel shows the system at the time at which half of the neutron star material has been accreted onto the black hole. We show densities down to . The bottom panel shows the remnant 5ms later, plotting densities down to and cutting out the , quadrant. The differences in scale between the 4 figures can be determined knowing that the size of the black hole is always .

Fig. 5 shows snapshot of simulations R14i0 and R12i0, the largest and smallest neutron stars considered here, both in the middle of the stellar disruption, and later. The larger neutron star disrupts far enough from the black hole for a large portion of the matter to be initially ejected into a tidal tail, but the smaller neutron star disrupts just outside of the ISCO of the black hole. The top-right panel of Fig. 5 in particular shows how close the smaller neutron star is to the hole when it disrupts. From this picture, the fact than any material remains outside of the hole after merger is surprising in itself, and an indication of how strong the effects of the black hole spin can be on infalling material.

Figure 6: Baryon mass remaining outside of the black hole as a function of time, for the 3 non-precessing cases R14i0, R13i0 and R12i0. We shift all the curves by the time at which half of the matter has been accreted onto the black hole.
Figure 7: Top: Average surface density 5ms after merger for the 3 non-precessing cases R14i0, R13i0 and R12i0. Bottom: Disk profile for simulation R14i0. Shown are the angular velocity (normalized by the velocity for circular orbits in this metric), the entropy and the sound speed (normalized by ) as a function of the circumferential radius .

Important differences are observed between the final state of these mergers. For the larger neutron star, of the matter remains outside of the black hole after merger. More importantly, about half that material has already formed a thick accretion disk, of about in radius and with peak density . The formation of a hot, thick disk is less obvious for the smaller neutron stars. The amount of material remaining outside of the black hole is by no means negligible ( of the neutron star, see Fig. 6), but the maximum density is about an order of magnitude lower. In fact, if we look at the average surface density as a function of radius (Fig. 7), we see no evidence of an accumulation of higher density material at lower radii (), while that feature is clearly visible for the larger neutron star. From these results, we can also infer that smaller neutron stars would probably be unable to form any long-lived remnant.

Evolutions including all the necessary microphysics (magnetic fields, neutrino cooling) will be necessary to determine how these disks evolve over longer time scales (). We can already see, however, that for all three configurations the material remaining outside of the black hole is hot (), and would be cooled by neutrino emission. In that respect, the differences in density between the remnants could be significant, as they modify the opacity of the disk to neutrino radiation, and thus the efficiency with which the disk can transfer its energy into neutrinos.

The properties of the final remnant are presented in Table 2. Comparing our results for the amount of material remaining outside of the black hole at late times with the predictions of Ref. Foucart (2012), we find good agreement (within of the neutron star mass) for the two smallest stars. The largest star forms a disk heavy enough that we are out of the range in which the predictions of Ref. Foucart (2012) are expected to be valid — and indeed, the disk formed in the simulation is significantly more massive that what Ref. Foucart (2012) would predict. We also find consistency between our simulations and Ref. Foucart (2012) on the neutron star radius below which no matter will remain outside of the black hole after merger (). A more careful examination of the differences between our numerical results and Ref. Foucart (2012) indicates that in the regime of high spin, high black hole masses considered here, the remnant mass probably has a steeper dependence on the radius of the neutron star than what would be guessed from Ref. Foucart (2012), a model fitted mostly to lower mass systems. However, these differences could also be explained by the expected variations in the remnant mass due to the internal structure of the neutron star (i.e. the fact that two neutron stars with the same radius but different internal structure will result in different post-merger disk masses), especially considering the fact that all of the simulations used to fit the model in Foucart (2012) had larger tidal Love number than the neutron stars from simulations R13i0 and R12i0. Overall, the magnitude of the differences between the numerical results and the model are roughly at the expected level. The final black hole masses and spins are also within the expected errors of existing analytical models, i.e. away from the values derived by Pannarale Pannarale (2012).

For all configurations, the disruption and merger of the neutron star occurs over (see Fig. 6). Mass accretion at later times is negligible compared with what is observed in lower mass ratio systems. When a disk forms, its main characteristics are however fairly similar to the lower mass ratio cases — except of course for the aforementioned lower densities and larger disk radii, which are a natural consequence of the higher black hole mass. Fig. 7 shows a few characteristics of the forming accretion disk for the most strongly disrupted case. The surface density peaks at a distance of from the black hole, and the disk extends to about 1. The orbital velocity profile is slightly sub-Keplerian (by about ), while the sound speed is and thus compatible with a thick, thermally supported disk of scale height . This is consistent with the actual scale height of the disk, . Finally, the inner edge of the disk is particularly hot: we plot an estimate of the entropy

(14)

(the effective constant is defined by ), and find for . Within the disk, we still have . As the disk settles down over , we would expect the entropy to exhibit a minimum at the peak of the surface density distribution, as was observed in lower mass ratio systems Duez et al. (2010); Foucart et al. (2011).

Ejecta

The ejection of unbound material by compact binary mergers is a prerequisite for some electromagnetic counterparts, most notably emissions due to the radioactive decay of the neutron-rich ejecta Roberts et al. (2011); Metzger and Berger (2012). This ejecta can be obtained through various physical processes: unbound material in the tidal tail, but also magnetically-driven Shibata et al. (2011) or neutrino-driven Dessart et al. (2009) winds. The study of winds goes beyond the scope of this article, as this requires accurate long-term evolution of the remnant disk and the inclusion of physical processes that are neglected in this work (magnetic fields, neutrino radiation). We will thus limit ourselves to the measurement of unbound material in the tidal tail.

Even the presence of ejected material in the tidal tail can be difficult to assess in general relativistic simulations, particularly at high mass ratios. Maintaining a high enough resolution in both the disk-forming region and the tidal tail is challenging, and in practice matter can only be reliably evolved up to a distance of a few hundred kilometers from the black hole. This is indeed one of the main disadvantage of any grid-based simulations when compared with smoothed particle hydrodynamics methods, which can easily follow the evolution of tidal tails. In a time-independent spacetime and when pressure forces are negligible, it is easy to determine whether material is unbound: if , then the material will escape to infinity (and is the Lorentz factor of the fluid at infinity). This condition is also a fairly good approximation for low-density material far away from the central black hole after a compact binary merger, but becomes more and more inaccurate as one gets close to the black hole, or densities in the tidal tail become higher.

One way to assess whether using the condition to find unbound material is accurate is to follow material over a sufficiently long period of time, and check that doesn’t vary much. In our simulations, however, this only occurs for before the material leaves the numerical grid, which leads to large uncertainties in the amount of unbound material, and its characteristics. A more detailed discussion of these issues will be presented in Deaton et al. (in preparation). Here, we limit ourselves to a discussion of measurements of at relative low radii () and over short timescales, and note the uncertainties due to these approximations.

Figure 8: Mass of the unbound material, as measured by the condition . For each configuration, is the time at which of the neutron star material has been accreted onto the black hole. The dashed vertical line represent the time at which has escaped the grid (the low density tidal tail of simulation R12i0 cannot be followed accurately for more than , at which point we stop measuring the mass of the ejecta).

On Fig. 8, we plot the amount of mass with on the grid (and more than away from the black hole). The most compact neutron star, simulation R12i0, naturally has the least material in a tidal tail: about , of which appear unbound after accretion onto the black hole begins ( after the time at which of the neutron star material has been accreted onto the black hole), with variations of only over the period over which measurements appear reliable. Moving up in stellar radius, simulation R13i0 offers the most reliable measurement of the ejected mass, with a stable value of long before matter starts flowing out of the grid (the total mass of the tidal tail, in this case, is ), and variations of over . From this run, we can also estimate the relative error due to finite numerical resolution in these measurements, and get (at our “medium” resolution, we found ). This is distinct from the uncertainties due to the approximation made when using as a proxy for finding whether material is bound or not, and appears to be the dominant source of error for simulations R12i0 and R13i0. Finally, the largest neutron star shows the most uncertain measurements. Velocities and densities in the ejecta are generally higher: the approximate method takes more time to become accurate, but material remains on the grid for a shorter amount of time. There also is material with flowing directly out of the forming accretion disk, presumably as a result of shocks during disk formation, which makes it impossible to have all of the potential ejecta in the range at any given time. Even by expanding the outer boundary of the grid by compared with the two other runs, we thus find that the measured only appear to settle at the time at which matter starts flowing out of the grid (and boundary effects might influence the properties of the ejecta). It is thus quite likely that the ejected mass is slightly larger than what is observed in the simulation. Overall, adding the two main sources of error (numerical resolution and use of ), we estimate that the ejected masses are

(15)
(16)
(17)
Figure 9: Distribution of asymptotic velocities for the unbound material of simulation R13i0. Two different times of the high-resolution simulation, and one time of the medium resolution simulation are shown.

From the measurements of , we can also determine the distribution of the velocity of the fluid at large distance from the black hole. This is shown in Fig. 9 for configuration R13i0. The distribution peaks at , and most of the ejecta has . The qualitative features of the velocity distribution appear fairly robust when we vary the time at which we measure , and the resolution of the simulation — with more uncertainties for the high-velocity tail of the distribution. The merger with the more compact neutron stars has a very similar velocity profile. The situation is quite different for , where about half of the ejecta initially appear to have . At this point, however, we do not have the ability to follow such material for a long enough period of time to assess the reliability of the velocity estimates of that last configuration, and have to limit ourselves to the observation that the ejecta appears more relativistic for the largest neutron star than in the other cases studied here.

Finally, from these measurements, we can estimate the kinetic energy of the ejecta, which would be available for future emission as it slows down in the interstellar medium. We should note that these results are only order of magnitude estimates, as these energies are sensitive to the high-velocity tail of the velocity distribution, which is poorly constrained in our simulations. Additionally, our energy estimates are particularly unreliable for simulation R14i0, due to the large amount of poorly resolved high-velocity material that is rapidly leaving the grid. We find for simulations (R12i0,R13i0,R14i0) respectively.

Even considering the large uncertainties in these measurements, it is interesting that we consistently find that in this region of the parameter space a few percents of a solar mass can be ejected from the system. This is indeed very different from the results obtained in the limit of low mass, low spin black holes (or for BNS Hotokezaka et al. (2012)), where only a negligible amount of material was found to be unbound. These results indicate that in the case of BHNS binaries, tidal disruption of the neutron star (when it occurs) is likely to be accompanied by the ejection of of neutron-rich material. Such outflows are promising for optical emissions due to the radioactive decay of neutron-rich elements in the ejecta, and the production of heavy elements resulting from r-process nucleosynthesis. These ejecta might even be detectable as a radio afterglow as unbound material decelerates in the interstellar medium Metzger and Berger (2012): the kinetic energy available for radio emission is indeed larger than in supernova explosions. However, the luminosity of the radio afterglow heavily depends on the density of the environment, and BHNS mergers are likely to occur in much lower density environments than supernova explosions. The deceleration of the ejecta in the interstellar medium would then occur over longer timescales, and remain harder to detect. Massive ejecta in BHNS mergers are also limited to high spin configurations, and should not be considered as the norm unless is standard for black holes in compact binaries. And the most energetic ejecta found here, in particular, is a very optimistic scenario for which, in addition of a high black hole spin, we used a very large neutron star. Finally, we should note that the detailed evolution of the tidal tail is likely to depend on details of the equation of state that our simple -law model cannot capture. This problem should thus be revisited carefully with a more realistic modeling of the neutron star fluid.

Gravitational Waveforms

Variations of the gravitational waveform emitted by BHNS mergers with the equation of state of the neutron star are mainly due to two effects: the small tidal distortion of the neutron star during the inspiral, which we discussed in Sec. V.1.1, and the cutoff in the gravitational-wave signal when the neutron star is disrupted (if tidal disruption occurs). In this section, we discuss measurements of these effects in our numerical simulations, while in the next section we will focus on their detectability by the Advanced LIGO detector.

Figure 10: Gravitational-wave emission [(2,-2) mode] for the non-precessing simulations R12i0 and R14i0, extrapolated to infinity and normalized by the ratio of the distance to the observer to the total mass of the binary .
Figure 11: Phase difference between the (2,-2) mode of the gravitational-wave emission of the non-precessing simulations. Both R12i0 and R13i0 are compared with R14i0. The dashed black lines show the matching interval, while the dashed green line shows the time at which the amplitude of the signal peaks for case R14i0 (the case in which the neutron star disrupts at the earliest time).

The effects of tides on the gravitational-wave signal of a BHNS binary before the disruption of the neutron star are expected to be fairly small: Damour et al. Damour et al. (2012) computed the phase difference due to tidal effects in the Fourier transform of the dominant [] mode of the waveform to 2.5PN order in the stationary phase approximation,

(18)

where is the standard PN parameter, is the frequency of the mode of the gravitational strain,

(19)

and is the tidal deformability parameter which for a BHNS binary is

(20)

We note that the sign of Eq. (18) is different from the one given in Damour et al. (2012), due to differences in the convention used to take the Fourier transform of the signal (we use , while the dephasing in the stationary phase approximation was derived with the opposite convention  Cutler and Flanagan (1994)). For and the binary parameters considered here, we get between simulations R12i0 and R14i0 (or, for the phase of the gravitational waveform in the time domain, ). However, these predictions have only been tested on BNS mergers up to , while the tidal disruption of the neutron star in the binaries considered here occurs at . In that regime, the PN expansion no longer appears convergent: the 2.5PN term is of the same order as the leading order term.

Measuring the small phase difference in our simulation is a challenging problem. Fig.10 shows the dominant (2,-2) mode of the gravitational strain for simulations R12i0 and R14i0, after the application of a time and a phase shift chosen to minimize the phase difference within a matching interval spanning 2 cycles of the radial oscillation frequency. Clearly, the waveforms are very similar up to the last cycle before the disruption of the larger neutron star. Fig. 11 shows the phase differences between the waveforms of simulations R12i0, R13i0 and R14i0. They remain under until before the peak of the gravitational-wave signal! Moreover, the differences appear dominated by the influence of the residual eccentricity, not by equation of state effects. In Sec. V.1.1, we showed that the tidal distortion of the neutron star follows fairly closely PN predictions, yet this is clearly insufficient to have a measurable effect on the gravitational-wave signal for most of the evolution.

Figure 12: Phase difference between the (2,-2) mode of the gravitational-wave emission of the non-precessing simulations as a function of the normalized gravitational-wave frequency . The dotted and dot-dashed curves correspond to the 1PN and 2.5PN predictions from Ref. Damour et al. (2012). All curves are matched at . At lower frequency, residual eccentricity in the simulation makes measurements of too noisy to be useful.

A small phase difference between the waveforms can however be hidden by the matching procedure used. Another way to compare the waveforms is to look at the phase as a function of the gravitational-wave frequency. This gets rid of the need to apply an arbitrary time shift to the waveforms. However, the computation of from the numerical waveform is fairly inaccurate, and even a small residual eccentricity can introduce a large noise in the resulting . Computing the difference between two numerical simulations can thus only be done once the evolution of the orbital frequency due to orbital decay becomes fast enough to dominate the effects of eccentricity. Fig. 12 shows measurements of for our numerical simulations. For , there are large oscillations due to the eccentricity of the orbit, and we cannot accurately measure . But in the frequency range , a phase shift of is clearly observed between simulations R12i0 and R13i0, and the same between R13i0 and R14i0. This result lies in between the predictions of the 1PN and 2.5PN approximations. As the 2.5PN results are barely outside of the expected numerical error, these simulations are not accurate enough to improve on the PN predictions for tidal effects at high frequency. But we can confirm that the dephasing obtained from the 2.5PN predictions is at least correct within a factor of up to frequencies (). The 1PN prediction also seem slightly more accurate than the 2.5PN results.

Figure 13: Spectrum of the gravitational-wave signal. is the spectrum of the dominant mode of the gravitational-wave signal as seen by an optimally oriented observer at . The dashed line shows the leading order PN behavior, , with amplitude matched to the numerical results. The Zero-Detuned High Power and Zero-Detuned Low Power noise curves of Advanced LIGO Shoemaker (2010) are also shown, together with 3 potential high-frequency tunings of the detector, at  Shoemaker (2010), and  Nicolas Smith-Lefebvre ().

The effects of the disruption of the neutron star by the tidal field of the black hole are easier to see in the numerical results. Fig. 13 shows the spectrum of the gravitational-wave signal for the three simulations R12i0, R13i0 and R14i0 (for an optimally oriented binary at ). The largest star disrupts earlier, and the gravitational-wave spectrum is cut slightly above . For the smallest star, the cutoff is at about . As opposed to tidal effects, this high-frequency cutoff is very well resolved in the simulations.

Name (km/s) (kHz)
R12i0 0.021 0.16 30 2.1
R13i0 0.017 0.15 45 1.8
R14i0 0.014 0.13 45 1.5
R14i20 0.013 0.13 60
R14i40 0.013 0.11 150
R14i60 0.013 0.09 345
Table 3: Gravitational-wave emission over the course of the simulation as measured at a radius , where is the ADM mass of the system at infinite separation. is the energy contained in the waves, their angular momentum, and the velocity kick given to that black hole. is the cutoff frequency of the gravitational-wave signal defined by Eq. 21.

Other physical quantities which can be extracted from the gravitational-wave signal are summarized in Table 3: the energy and angular momentum radiated, the final kick given to the system, and a cutoff frequency defined (arbitrarily) by the relation

(21)

with (Note that any value of in the range gives nearly identical result as is approximately constant during the inspiral, as shown in Fig. 13). As the binding energy of these systems at is , we see that these binaries will radiate of their energy before merging, and of their angular momentum (the more compact neutron stars naturally radiating more, as they disrupt later). The final kicks remain low (), as is generally observed for non-precessing BHNS binaries.

Detectability of the Neutron Star Radius by Advanced LIGO

Keeping in mind the results of the previous section, we can begin to address another important question: the measurability of finite size effects on the gravitational waveform of BHNS mergers for mass ratios . An earlier analysis of these issues by Lackey et al. Lackey et al. (2012) showed that at lower mass ratios () and for nonspinning black holes Advanced LIGO would be sensitive to differences in the radius of the neutron star of order for an optimally oriented BHNS merger located at . This is due in part to the effects on the waveform of the tidal distortion of the neutron star, and in part to variations in the binary separation at which the neutron star disrupts and the gravitational-wave signal is cut off.

At higher mass ratio, tidal effects are smaller. However, the disruption of the neutron star occurs at a lower frequency and the amplitude of the gravitational-wave signal is larger. It is thus unclear whether finite size effects will be easier or harder to detect. On Fig. 13, we show the spectrum of the gravitational-wave signal as seen by an optimally oriented observer located away from the binary, and compare it with different Advanced LIGO detector’s strain noise spectra (see below). At that distance the differences between the three simulations seem to be marginally measurable. In the rest of this section, we will attempt to quantify this statement more carefully.

To determine whether the difference between two waveforms and can be detected by Advanced LIGO, we use the approximate condition Lindblom et al. (2008)

(22)

where the inner product is defined as

(23)

Here and are the Fourier transforms of two waveforms and , and is the one-sided power spectral density of the detector’s strain noise, defined as

(24)

where is the noise correlation matrix for zero-mean, stationary noise. In this case, we will consider three of the Advanced LIGO guideline noise curves defined in Ref. Shoemaker (2010): the Zero Detuned Low Power spectrum, which is the expected sensitivity of the detector once signal recycling mirrors are in place, the Zero Detuned High Power spectrum, which is the final design sensitivity of Advanced LIGO, and a High Frequency noise curve optimized to take data at . We also consider alternative tunings of the Advanced LIGO detector to and , using noise curves graciously provided to us by Nicolas Smith-Lefebvre Nicolas Smith-Lefebvre (), and generated by the noise simulation package ‘gwinc’ developed by the LIGO collaboration. These noise curves assume an Advanced LIGO detector in the same configuration as the High Frequency model of Ref. Shoemaker (2010), except that the signal-recycling mirror detuning is chosen to tune the detector to higher frequencies.

When taking the inner product , we choose one polarization of and then allow for a time and phase shift in , chosen to maximize the inner product (which is identical to minimizing for ). Finally, we consider only the quadrupolar part of the waveform as measured by an optimally oriented observer:

(25)

In theory, the integration in Eq. (23) should be carried over the entire frequency band of the detector. This is however complicated for our waveforms, as the numerical simulations only cover the high-frequency part of the signal (). To construct a full waveform, the numerical results should be hybridized with some analytical approximation valid at low frequency (PN, Effective One-Body,…). But in the region of parameter space considered here (, ), the error coming from extending these approximants to frequencies is significantly larger than the actual differences expected between waveforms for . Uncertainties in the construction of the hybrid then dominate the measurement of .

We will instead consider three approximations to :

  • Limiting the integration to frequencies , where the waveform is known exactly from the numerical simulations. This neglects tidal effects during the inspiral.

  • Limiting the integration to frequencies , and using PN predictions over the entire frequency range (see Appendix B). In this case, we ignore the effects of the disruption of the neutron star, as well as errors in the PN predictions at high frequencies. As seen in the previous section, the 1PN predictions appear to remain fairly accurate for , and provide at least a good qualitative estimate of the tidal effects during the inspiral.

  • Combining the low-frequency PN predictions with the numerical results at high frequency by matching in the frequency domain the phase difference between two simulations to the predicted PN phase difference. This procedure is detailed in Appendix B.

The matching procedure and the differences between various PN orders each cause uncertainties of in .

NR Only
Zero Det H-P 1.5 0.7 0.9 Zero Det L-P 0.7 0.3 0.4 High Freq() 1.2 0.4 0.9

1PN Only
Zero Det H-P 1.0 0.4 0.6 Zero Det L-P 0.5 0.2 0.3 High Freq() 0.4 0.2 0.3

Hybrid-1PN
Zero Det H-P 2.5 1.2 1.5 Zero Det L-P 1.2 0.6 0.7 High Freq() 1.5 0.6 1.0 High Freq() 3.2 1.4 1.9 High Freq() 2.6 1.3 1.4

Table 4: Detectability of tidal effects in non-precessing BHNS mergers for 3 different Advanced LIGO noise curves from Ref. Shoemaker (2010): Zero-Detuned High Power, Zero-Detuned Low Power and High Frequency tuned at (see also Fig. 13). The quantities , defined in Eq. 22, are the detectability criteria for optimally oriented events at . We consider first the results for the numerical waveform limited to (NR Only), then for the first order PN expansion [Eq. (18)] limited to (1PN Only) and finally for hybrid waveforms matched in spectral space between (Hybrid-1PN). For each case, we compare simulations [R12i0,R14i0] (), [R12i0,R13i0] () and [R13i0,R14i0] (). For the hybrids, we also give results for noise curves tuned to and .

The results for each method are summarized in Table 4. We find that for the Zero Detuned noise curves, the high frequency cutoff is only more important than the low-frequency tidal effects. The tuned High Frequency noise curves, quite naturally, are much more sensitive to the disruption of the neutron star than to tidal effects. For the Zero Detuned High Power noise curve, and using our estimates of the mismatch over the entire LIGO band, the detectability criteria is satisfied if the neutron stars have radii differing by

(26)

where is the distance to the observer. The Low Power noise curve requires the binary to be about twice as close. However, using the low-power detector with high frequency tuning (at ) leads to larger than for the Zero Detuned High Power noise curve by .

Differentiating a binary black hole system from a BHNS binary would only be slightly easier. The phase difference between a point-particle waveform and the waveform of simulation R14i0 is, during inspiral, only times larger than the phase difference between simulations R12i0 and R14i0. And the high-frequency cutoff in the waveform would not help us much, as the most compact star considered here (R12i0) does not disrupt very far from the ISCO.

These results are not very promising. For the Zero Detuned High Power noise curve and the waveforms from simulations R12i0 and R14i0, the condition will only be satisfied for binaries with signal-to-noise ratio , or about of the Advanced LIGO events (assuming that BHNS binaries are equally distributed in volume). A detector tuned at would reach the condition for about twice as many events. The criteria is also an optimistic limit: it does not take into account the fact that all other parameters of the binary are here assumed to be known. Uncertainties in the masses of the objects and the spin of the black hole will significantly affect these results, making it harder to detect equation of state effects. Additionally, both the PN dephasing and the variations in the high-frequency cutoff of the signal are helped by the fact that we are considering a rapidly rotating black hole. For a nonspinning black hole, the neutron star would reach the ISCO at lower frequencies, and tidal effects during the inspiral would be even smaller. And as the neutron star would plunge into the black hole before being disrupted, there is no guarantee that there would be a measurable difference in the cutoff frequency of the waveform (although this question certainly deserves further investigation: the merger would also occur in a more favorable frequency range). Finally, any real data analysis would require the knowledge of the waveform at low frequency, which is not at this point known with enough accuracy for binaries with and . The theoretical detectability conditions considered here are thus certainly too generous.

v.2 Precessing Binaries

Our second set of BHNS mergers considers variations of the orientation of the black hole spin. The starting configuration is the largest neutron star studied in the previous section. The black hole spin is inclined with respect to the orbital angular momentum by . In all cases, the misaligned component of the spin is initially along the line connecting the black hole and neutron star centers. We do not expect this choice to affect the qualitative feature of the merger (disruption, disk formation) Foucart et al. (2011). It should, however, influence the magnitude of the velocity kick given to the final black hole as a result of gravitational-wave emission Lousto and Zlochower (2009). Modifying the initial separation of the binary is effectively identical to a change in the orientation of the black hole spin at constant . Indeed, for a binary in which only one object is spinning, is conserved Apostolatos et al. (1994) and changing the initial separation only modifies the phase of the precession of the binary.

Inspiral: Orbital Precession

Figure 14: Inclination between the initial and current direction of the normal to the orbital plane.

During inspiral, the main difference with the aligned configurations will naturally be the precession of both the black hole spin and the orbital angular momentum around the total angular momentum of the system. A simple coordinate measurement of that precession is presented in Fig. 14, using the angle between the direction of the initial orbital angular momentum and the normal to the orbital plane (defined as , where are the coordinates of the centers of the compact objects, and their velocities).

We see that the amplitude of the precession of the orbital plane vary from for R14i20 to for the most inclined case. If the total angular momentum of the system was constant (i.e. in the absence of gravitational-wave emission), we would expect the amplitude of that precession to be twice the initial angle between the orbital angular momentum and the total angular momentum, i.e. and for simulations R14i20 and R14i60 respectively. Somewhat larger values are expected for radiating systems, as the loss of angular momentum due to gravitational-wave emission is to first order aligned with the current orbital angular momentum of the binary: taking into account the angular momentum lost to gravitational waves listed in Table 3 would correct these estimates to and (which are now overestimates, as part of the radiated angular momentum is emitted during the disruption of the neutron star). Over the course of the binary evolution, the most inclined binary goes through slightly more than half of a precession period, while R14i20 gets close to completing a full precession period (as is defined with respect to the initial orbital plane, we get after a full precession period, and not after half a precession period).

A simple comparison of the measured precession of the orbital plane with the PN predictions from Faye et al. (2006) is also presented on Fig. 14. The 1PN and 2PN curves are obtained by evolving the initial BH spin and orbital angular momentum using the PN formulae from Faye et al. (2006), but assuming that the trajectories of the compact objects are those observed in the simulation (i.e. when computing the relative position and velocity of the objects, we use our results and not a PN evolution of the initial conditions, as the PN equations of motion are quite inaccurate so close to merger). We see that the period over which the orbital plane precesses in the simulations matches the PN predictions very well, while the amplitude of the precession is smaller in the numerical results. Such differences (as well as the additional oscillations in the value of the inclination angle) could however easily be due to the fact that these measurements are certainly not gauge-independent.

A large precession of the orbital plane is quite natural for the high mass ratio, high spin systems considered here. Indeed, the angular momentum of the black hole is while the initial orbital angular momentum is only (and about of the total angular momentum, or of , is radiated in gravitational waves). This means that even for relatively low inclinations of the black hole spin (such as in simulation R14i20), large variations of the orientation of the orbital plane occur. These oscillations and, more importantly, the shift in the phase of the gravitational waveform with respect to aligned spin templates which accompany them, can make the detection of precessing binaries challenging. The higher dimensionality of the parameter space to consider also complicates parameter estimates — although there are also positive effects due to precession: some of the degeneracies existing between the parameters of the system for aligned binaries are broken for precessing systems Damour et al. (2012). As mentioned in the previous section, an inaccurate determination of the parameters of the system increases the error in any measurement of the neutron star equation of state from gravitational waveforms. Obtaining proper constraints on the parameters of a BHNS binary with misaligned spin, which requires reliable templates for precessing systems, is thus a prerequisite to any attempt at constraining the equation of state of neutron stars from BHNS waveforms, at least for the high-spin configurations considered here (unless the spin of the black hole is aligned with the orbital angular momentum of the system by some unknown mechanism during the pre-merger evolution of the binary).

Disruption and Disk Formation

Figure 15: Same as Fig. 5, but for the precessing binaries R14i20 (left) and R14i40 (right).

The effect of spin misalignment on the properties of the remnant of a BHNS merger were first studied in a general relativistic framework in Ref. Foucart et al. (2011). Material on an orbit inclined with respect to the equatorial plane of the black hole reaches the region in which stable orbits no longer exist at a larger separation than material in the equatorial plane on a corotating orbit. Effectively, this means that BHNS mergers with misaligned black hole spins are roughly equivalent to mergers with a lower black hole spin aligned with the orbital angular momentum. Disruption becomes less likely for misaligned configurations, and the mass remaining outside of the black hole at late times is smaller. The smallest radius at which stable orbits exist for black holes with at inclination is equal to the radius of the innermost stable circular orbit of a black hole with aligned spin  Hughes (2001); Fragile et al. (2007). 2 And indeed, the disruption of the neutron star is close to what is expected for such spins: for the low inclination case R14i20, the mass of material remaining outside of the black hole after merger is nearly identical to what was observed in the aligned configuration R14i0 while in simulations R12i40 and R14i60, and of the neutron star mass remain outside of the hole after merger. This is very similar to what our simple fitting model Foucart (2012) would predict ( and of the neutron star matter surviving the merger for ), or what might be inferred from numerical simulations for smaller black hole spins Foucart et al. (2012) (which found of the matter remaining outside of the hole for and none for , all other parameters being identical to the cases considered here). It thus seems fairly likely that, for the purpose of measuring the mass of material remaining outside of the black hole at late times at least, the disruption of the neutron star in BHNS mergers with misaligned black hole spins can be modeled with good accuracy by considering the results of aligned configurations only.

However, there are important qualitative differences between the behavior of aligned and misaligned configurations. Fig. 15 shows snapshot of the two simulations with the lowest inclination angle for the black hole spin (R14i20 and R14i40), at a time at which of the neutron star has been accreted onto the black hole (top) as well as later (bottom). These can be directly compared with Fig. 5 for non-precessing systems. Simulation R14i20 mostly behaves as the aligned case, although the slight precession of the tidal tail with respect to the disk induces small differences in the formation of the disk, and should affect its subsequent evolution as matter falls back from changing directions. At moderate inclinations (R14i40), the changes are more drastic. There is not much of a disk forming: most of the remaining material is in a long tidal tail, which is differentially precessing. For aligned configurations, a disk generally starts to form as the front edge of the accretion flow wraps around the black hole and hits the material from the tidal tail, creating a shock, outflows, and a rapid redistribution of the tidal tail material. Disk-tail interactions are much less visible in inclined simulations (although there are still contacts between the inner and outer edge of the tidal tail as it wraps around the black hole). Existing shocks are however still sufficient to heat the remaining material in simulations R14i20 and R14i40 to an average temperature . Not surprisingly the remnant of the most inclined merger (R14i60), which does not form a disk, is much cooler (). Finally, we find that precession does not prevent the formation of a baryon-poor region along the rotation axis of the black hole: a cone of opening angle is clear of material at densities , at least over the short timescale over which the post-merger remnant is evolved. Close to the black hole, lower density material is not resolved in the simulation. This baryon-poor region was also shown to exist after the merger of precessing BHNS binaries at lower mass ratio Foucart et al. (2011).

The significant asymmetry of this system with respect to the equatorial plane of the black hole, as well as the differential precession between fluid elements, could also have important consequences for the long term evolution of the system. Disk simulations by McKinney et al. McKinney et al. (2013) indicate that coupling between the magnetic field in an accretion disk (and relativistic jet) and the spin of the black hole leads to an alignment of the inner disk with the spin of the black hole, while the outer disk remains misaligned. Similarly, the jet is emitted along the rotation axis of the black hole, but orthogonal to the plane of the outer disk at large distances. How these effects will play out for the more compact disk produced by BHNS mergers is an open question. Another important consequence of spin-orbit misalignment was pointed out by Etienne et al. Etienne et al. (2012b), who showed that asymmetries and motion across the equatorial plane of the black hole contribute to the formation of a larger toroidal field within the disk formed by BHNS mergers (aligned configurations with equatorial symmetry, on the other hand, form disks with mostly poloidal magnetic fields, as the magnetic field comes from the winding of the original field lines frozen within the disrupting neutron star). The toroidal field is amplified by the fastest growing mode of the magnetorotational instability(MRI), the mechanism most likely to allow the magnetic field in the accretion disk to grow up to levels at which jets (and short gamma-ray bursts) could be created. Misaligned black hole spins could thus lead to qualitative difference in the evolution of magnetized disks. They would also make it easier to resolve the growth of the MRI numerically, as the wavelength of the fastest growing MRI mode scales with the magnitude of the toroidal component of the magnetic field Penna et al. (2010). On the other hand, the disks formed from precessing BHNS binaries have lower densities and lower masses. Compared to aligned configurations, a larger fraction of the matter remaining outside of the black hole is sent in the tidal tail. The most inclined configuration studied here (R14i60), which barely disrupts, actually keeps less than of the neutron star material within of the black hole. Nearly all of the remnant mass ( of the neutron star material) is on highly eccentric, differentially precessing orbits.

If the velocity kicks given to neutron stars during supernova explosions are such that a majority of BHNS systems are within the range of inclination of the black hole spin studied here (as might be expected from the results of Belczinsky et al. Belczynski et al. (2008)), the effects of inclination on the orbital evolution of the binary, and consequently on the gravitational-wave signal, would be significant. But the conditions required for the disruption of the neutron star to occur in a significant number of systems would not be dramatically modified from those derived in non-precessing systems: the location of the marginally stable orbit does not change much for inclinations (see e.g. Fragile et al. (2007)), at least for the relatively large spins which we already know are necessary for tidal disruption to occur for mass ratios . The remnant disks would be less massive, but with larger toroidal magnetic fields initially. Considering that even low mass disks () are energetically sufficient to power gamma-ray bursts if the conditions (temperature, magnetic fields) are otherwise right Lee et al. (2005); Giacomazzo et al. (2012), misaligned configurations could in the end prove more favorable than the aligned cases.

Waveforms

Figure 16: Gravitational waveforms as measured by observers who, at , see the binary face-on (solid black line) and edge-on along the line connecting the centers of the two objects (dashed red line). is normalized by the ratio of the distance to the observer to the total mass of the binary . We show both the and polarization, with shifted by . Top: Non-precessing configuration R14i0. The face-on observer is always optimally located, while the edge-on observer receives a weaker signal in (in which higher order modes are however more visible), and no signal at all in . Bottom: Precessing configuration R14i60. The envelope of the signal varies in time as the binary precess, and the optimal orientation is modified accordingly. The merger also occurs at an earlier time, as the orbital hang-up is only due to the aligned component of the black hole spin.

The waveforms from precessing BHNS binaries are significantly different from their non-precessing counterparts, as previously mentioned: the precession of the orbital plane shown in Fig. 14 will cause a modulation of the preferred direction for gravitational-wave emission (see Fig. 16), which has to be properly modeled in order to avoid significantly reducing our sensitivity to any waveform emitted by a precessing system. The modeling of precessing waveforms goes beyond the scope of this article. A more detailed study of the impact of precession on the detectability of BHNS mergers can be found in Brown et al. Brown et al. (2012). Assuming an isotropic distribution of black hole spin, about half of the BHNS binaries within the theoretical range of the next generation of gravitational-wave detectors could be missed given the performance of the current search methods in the regime of high mass ratio, strongly precessing systems. Attempts to reduce the complexity of the problem by studying the waveforms in a preferred frame precessing with the binary are under way Schmidt et al. (2011); Boyle et al. (2011); Schmidt et al. (2012), and could help in the construction of future precessing templates, but the detection of BHNS systems in Advanced LIGO remains an important challenge today (see also Ajith et al. (2012) for an updated template bank for binaries with arbitrary spins). Considering the negligible influence of tidal effects on the waveform before the disruption of the neutron star, these issues are however better studied in the context of black-hole–black-hole binaries, for which longer and more accurate precessing waveforms are available. Results obtained for these systems should be immediately applicable to BHNS inspirals, at least in the range of mass ratios considered in this work (tidal effects are more important for more symmetric mass ratios).

An interesting particularity of BHNS mergers with high black hole spins, precession, and a high enough mass ratio that tidal disruption either does not occur or only occurs very close to the marginally stable orbit, is the possibility for the remnant black hole to receive a significant velocity kick from the merger. This effect is already well-known for binary black hole systems Zlochower et al. (2011); Lousto and Zlochower (2011): binary black holes with parameters similar to the BHNS mergers studied here would receive velocity kicks (and for more favorable parameters). Kicks in BHNS mergers are generally much smaller, with , even in the precessing configurations that we previously studied Foucart et al. (2011). This is due to the fact that most of the kick is received at the time of merger, while in BHNS systems at lower mass ratio the neutron star disrupts before merger, effectively cutting off the asymmetric gravitational-wave emission responsible for the kick. Here, however, the disruption of the neutron star occurs very late in the inspiral — or nearly not at all in the case of the most inclined configuration. The strong distortion of the neutron star as it plunges into the black hole will still reduce gravitational-wave emission at merger, even for configurations in which no matter remains outside of the black hole afterwards, but not nearly as much as for aligned configurations, or for more symmetric mass ratios. Accordingly, we find that much larger kicks are now possible. It should also be noted that the kick obtained from binary mergers is proportional to , with the orbital phase at merger and an unknown phase shift. The only way to measure the maximum kick from a specific configuration is thus to consider a sequence of mergers, spanning a range of phases . All that can be said from a single configuration is that kicks larger than are possible. Without more studies of the dependence of in , we cannot in principle exclude the possibility that kicks for these configurations could be nearly as large as in binary black hole mergers, although maximal values closer to those presented here appear more likely.

Vi Conclusions

We continue our study of BHNS mergers at mass ratio , the regime currently deemed to be the most likely for BHNS binaries in the field. Previous results Foucart et al. (2012) have shown that for such mass ratios, high black hole spins are needed for the neutron star to disrupt, even for large neutron stars (). In this work, we look into the influence of the radius of the neutron star and the orientation of the black hole spin on the merger, focusing on configurations with mass ratio and black hole spin . We show that the transition between configurations for which the disruption of the neutron star causes the formation of massive accretion disks and those where the neutron star just plunges into the black hole is very sensitive to the compactness of the neutron star: while a neutron star of radius leads to the formation of a massive disk () and tidal tail (), a smaller neutron star in an otherwise identical binary forms a much less massive remnant (, ). For neutron stars with radii , we expect no disruption at all. This indicates that in the range of stellar radii currently favored Steiner et al. (2010), a black hole spin is required for a neutron star to disrupt — and thus for post-merger electromagnetic counterparts such as SGRBs and kilonovae to be possible. We also note that for all but the most massive disks a fairly low maximum density is observed, about an order of magnitude lower than for disks of similar masses at lower mass ratios . This is simply the expected geometrical effect: the radius of the disk is roughly proportional to the mass of the final black hole. But that difference could significantly affect the late time evolution of the disk: the opacity of the disk to neutrino radiation will be lower, and the evolution of the magnetic field is likely to be affected as well.

The amount of unbound material, approximated through measurements of the energy of fluid elements in the limit of a time-independent metric, is found to be larger for these high-spin configurations than in previous, lower spin studies of BHNS mergers. The accuracy of these mass measurements is only in our general relativistic simulations — but for all three equations of state studied here, we find ejected mass . The ejecta has a velocity distribution peaking at (except for the larger neutron star, for which we find larger velocities which cannot be accurately measured at this point) and a kinetic energy . It would be a promising setup for a potential “kilonova”, and could even be detected as a radio afterglow. It should however be emphasized that these massive ejecta are only produced in the high spin region of the parameter space, as the neutron star does not disrupt for low black hole spins.

The effect of a misaligned black hole spin is also studied in more details. General relativistic simulations of precessing BHNS binaries had only been performed for one set of binary parameters before this work, for low mass, low spin black holes (,  Foucart et al. (2011)). These high spin configurations allow us to observe the effect of the misalignment of the black hole spin on tidal disruption more accurately. In particular, we confirm that using the radius of the innermost stable spherical orbit  Hughes (2001) as a way to predict the mass remaining outside the black hole at late times is reasonable: it matches the results of an aligned configuration with effective aligned spin and