# Hydrodynamics of black hole–neutron star coalescence.

## 1 Introduction

Theoretical studies of the binary coalescence of a neutron star with a black hole are as venerable as the Texas Symposia (Wheeler 1971 ). Such coalescences have been implicated in the creation of heavy nuclei through the r–process and were suspected of causing the violent outbursts known as gamma–ray bursts (Lattimer & Schramm 1974, 1976 ; they are also candidate sources for the gravitational radiation detectors currently coming of age, as the double neutron star binaries have been for a long time (Clark & Eardley 1977 ).

No binary system composed of a black hole and a neutron star has as yet been identified as such by the astronomers. However, theoretical estimates and studies of binary stellar evolution predict that such systems are created, and that the expected rate of their coalescence may be between about one per one hundred thousand and one per a million years per galaxy (Lattimer & Schramm 1974, 1976 ; Narayan, Piran & Shemi 1991 ; Tutukov & Yungelson 1993 ; Lipunov, Postnov & Prokhorov 1997 ). This rate is about right to explain the r–process and gamma–ray bursts (although it remains to be demonstrated that the phenomena are related to the coalescence), and would give a very satisfactory event rate for future gravity wave detectors (10 to 100 per year out to 200 Mpc).

The expected outcome of the coalescence ranged from the neutron star “plunging” into the black hole (Bardeen, Press & Teukolsky 1972 ), through the disruption of the neutron star and the formation of a transient accretion stream (Wheeler 1971 , Lattimer & Schramm 1976 ) or of an accretion disk/torus (Paczyński 1991 , Jaroszyński 1993 , Witt et al. 1994 ), to rapid growth of the binary separation in a steady Roche-lobe overflow scenario (Blinnikov et al. 1984 , Portegies Zwart 1998 ). It seemed important to investigate the hydrodynamics of the process numerically, to determine which, if any, of these outcomes are likely.

We have initially performed Newtonian simulations treating the neutron star as a stiff polytrope (Lee & Kluźniak 1995 ), motivated by a desire to determine the timescale of the coalescence with the black hole and to investigate the spatial distribution of the matter lost by the star. These questions were of relevance to the theory of gamma–ray bursts, which seems to allow the coalescence to power the burst in the accepted relativistic shock model, provided that the environment into which the fireball expands has a very small number of baryons at least in some directions (Mészáros & Rees 1992, 1993 ), and that the explosive event is highly variable in time (Sari & Piran 1997 ). We found the coalescence of a neutron star and a black hole promising in these respects, at least in our Newtonian model (Kluźniak & Lee 1998 ).

However, the stiff polytrope ($`\mathrm{\Gamma }=3`$), used by us in the work cited above, differed in one crucial respect from a neutron star—it responded to mass loss by shrinking, instead of expanding. For this reason we have repeated the simulation with a softer polytrope, and we report the results below, and elsewhere (Lee & Kluźniak 1999 ).

## 2 Numerical method

For the simulations presented here, we have used the method known as Smooth Particle Hydrodynamics. Our code is three–dimensional and is completely Newtonian, although removal of angular momentum by gravitational radiation was included, as described below. This method has been described often, we refer the reader to Monaghan (1992)  for a review of the principles of SPH, and to Lee (1998)  and Lee & Kluźniak (1998)  for a detailed description of our own code.

Here, we model the neutron star via a polytropic equation of state, $`P=K\rho ^\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }=5/3`$. The unperturbed (spherical) neutron star has a radius R=13.4 km and mass M=1.4M. The black hole (of mass $`M_{\mathrm{BH}}`$) is modeled as a Newtonian point mass, with a potential $`\mathrm{\Phi }_{\mathrm{BH}}(r)=GM_{\mathrm{BH}}/r`$. We model accretion onto the black hole by placing an absorbing boundary at the Schwarzschild radius ($`r_{Sch}=2GM_{\mathrm{BH}}/c^2`$). Any particle that crosses this boundary is absorbed by the black hole and removed from the simulation. The mass and position of the black hole are continously adjusted so as to conserve total mass and total momentum.

Initial conditions corresponding to tidally locked binaries in equilibrium are constructed in the co–rotating frame of the binary for a range of separations r and a given value of the mass ratio $`q=M_{\mathrm{NS}}/M_{\mathrm{BH}}`$ (Rasio & Shapiro 1994 ; Lee & Kluzniak 1998 ). The neutron star is modeled with $`N=17,256`$ particles at the start of the calculation in every case presented here. We have calculated the gravitational radiation signal emitted during the coalescence in the quadrupole approximation, and refer the reader to Lee & Kluźniak (1999)  for details.

We have included a term in the equations of motion that simulates the effect of gravitational radiation reaction on the components of the binary system in the quadrupole approximation (see Landau & Lifshitz 1975 ). This formulation of the gravitational radiation reaction has been used in SPH before (Davies et al. 1994 , Zhuge et al. 1996 , Rosswog et al. 1999 ) in the case of merging neutron stars, and it is usually switched off once the stars come into contact, when the point–mass approximation clearly breaks down. We are assuming then, that the polytrope representing the neutron star can be considered as a point mass for the purposes of including radiation reaction. Clearly, the validity of this assumption must be verified a posteriori when the simulation has run its course. If the neutron star is disrupted during the encounter with the black hole, this radiation reaction must be turned off, since our formula would no longer give meaningful results. We have adopted a switch for this purpose, as follows: if the center of mass of the SPH particles comes within a prescribed distance of the black hole (effectively a tidal disruption radius), then the radiation reaction is turned off. This distance is set to $`r_{tidal}=CR(M_{\mathrm{BH}}/M_{\mathrm{NS}})^{1/3}`$, where C is a constant of order unity.

## 3 Results

### 3.1 Evolution of the binary

To allow comparisons of results for differing equations of state, we have run simulations with the same initial binary mass ratios as previously explored for $`\mathrm{\Gamma }=3`$ (Lee & Kluźniak 1998 ), namely $`q`$=1, $`q`$=0.8 and $`q`$=0.31. Additionally we have examined the case with mass ratio $`q`$=0.1. Equilibrium sequences of tidally locked binaries were constructed for a range of initial separations, terminating at the point where the neutron star overflows its Roche Lobe (at $`r=r_{RL}`$). In Figure 1a

we show the variation of total angular momentum J for one of these sequences as a function of binary separation (solid line). Following Lai, Rasio & Shapiro (1993b) , we have also plotted the variation in J that results from approximating the neutron star as compressible tri–axial ellipsoid (dashed lines) and as a rigid sphere (dotted lines). In all cases, the SPH results are very close to the ellipsoidal approximation until the point of Roche–Lobe overflow. This result is easy to understand if one considers that the softer the equation of state, the more centrally condensed the neutron star is and the less susceptible to tidal deformations arising from the presence of the black hole.

For $`\mathrm{\Gamma }=3`$ (Lee & Kluźniak 1998 ), the variation in angular momentum as a function of binary separation was qualitatively different (for high mass ratios) from our present findings for $`\mathrm{\Gamma }=5/3`$. For $`q`$=1 and $`q`$=0.8, total angular momentum attained a minimum at some critical separation before  Roche–Lobe overflow occurred. This minimum indicated the presence of a dynamical instability, which made the binary decay on an orbital timescale. This purely Newtonian effect arose from the tidal interactions in the system (Lai, Rasio & Shapiro 1993a ). In the present study, we expect all orbits with initial separations $`rr_{RL}`$ to be dynamically stable.

For polytropes, the mass–radius relationship is $`RM^{(\mathrm{\Gamma }2)/(3\mathrm{\Gamma }4)}`$. For $`\mathrm{\Gamma }`$=5/3, this becomes $`RM^{1/3}`$. Thus, the polytrope considered here responds to mass loss by expanding, as do neutron stars modeled with realistic equations of state (Arnett & Bowers 1977 )–the dynamical disruption of the star reported below seems to be related to this effect. For the polytropic index considered in Lee & Kluźniak (1998) , the star was not disrupted (see also Lee & Kluźniak 1995 ; 1997 ; Kluźniak & Lee 1998 ), but we find no evidence in any of our dynamical calculations for a steady mass transfer in the binary, such as the one suggested in the literature (e.g. Blinnikov et al. 1984 ; Portegies Zwart 1998 ).

Using the quadrupole approximation, one can compute the binary separation as a function of time for a point–mass binary, and obtain

$`r=r_i\left(1t/t_0\right)^{1/4},`$ (1)   
with $`t_0^1=256G^3M_{\mathrm{BH}}M_{\mathrm{NS}}(M_{\mathrm{BH}}+M_{\mathrm{NS}})/(5r_i^4c^5)`$. Here $`r_i`$ is the separation at $`t`$=0. For black hole–neutron star binaries studied here, the timescale for orbital decay because of angular momentum loss to gravitational radiation, $`t_0`$, is on the order of the orbital period, $`P`$ (for $`q`$=1, at an initial separation $`r_i`$=2.7R we find $`t_0`$=6.5 ms and $`P`$=2.24 ms).

### 3.2 Run parameters

In Table 1 we present the parameters distinguishing each dynamical run we performed. All times are in milliseconds and all distances in kilometers. The runs are labeled with decreasing mass ratio (increasing black hole mass), from $`q`$=1 down to $`q`$=0.1. All simulations were run for the same length of time, $`t_{final}=22.9`$ ms (this covers on the order of ten initial orbital periods for the mass ratios considered). The initial separation for each dynamical run is given as $`r_i`$, and the separation at which Roche Lobe overflow from the neutron star onto the black hole occurs is given by $`r_{RL}`$.

The fifth column in Table 1 shows the value of $`t_{rad}`$, when radiation reaction is switched off according to the criterion established in section 2. We note here that run E is probably at the limit of what should be inferred from a Newtonian treatment of such a binary system. The black hole is very large compared to the neutron star, and the initial separation ($`r_i=67.87`$ km) is such that the neutron star is well within the innermost stable circular orbit around a Schwarzschild black hole of the mass considered.

### 3.3 Morphology of the mergers

The initial configurations are close to Roche Lobe overflow, and mass transfer from the neutron star onto the black hole starts within one orbital period for all runs presented here. In every run the binary separation (solid lines in in Figure 1) initially decreases due to gravitational radiation reaction. For high mass ratios, (runs A, B) the separation decays faster than what would be expected of a point–mass binary. This is also the case for a stiff equation of state, in black hole–neutron star mergers (Lee & Kluzniak 1998 ) as well as in binary neutron star mergers (Rasio & Shapiro 1994 ), and merely reflects the fact that hydrodynamical effects are playing an important role. For the soft equation of state studied here, there is the added effect of ‘runaway’ mass transfer because of the mass–radius relationship (see section 3.1). For run C, the solid and dashed lines in Figure 1b follow each other very closely, indicating that the orbital decay is primarily driven by angular momentum losses to gravitational radiation. For run E, the orbit decays more slowly than what one would expect for a point–mass binary. This is explained by the fact that there is a large amount of mass transfer (10% of the initial neutron star mass has been accreted by $`t=t_{rad}`$ in this case) in the very early stages of the simulation, substantially altering the mass ratio in the system (the dashed curves in Figure 1b are computed for fixed masses; at constant total mass, note that from equation 1, lowering the mass ratio in the system slows the orbital decay for $`q<0.5`$).

The general behavior of the system is qualitatively similar for every run. Figure 2 shows density contours in the orbital plane (left panels) and in the meridional plane containing the black hole (right panels) for run D at $`t=5.73`$ ms and $`t=t_f=22.9`$ ms.

The neutron star becomes initially elongated along the binary axis and an accretion stream forms, transferring mass to the black hole through the inner Lagrange point. The neutron star responds to mass loss and tidal forces by expanding, and is tidally disrupted. An accretion torus forms around the black hole as the initial accretion stream winds around it. A long tidal tail is formed as the material furthest from the black hole is stripped from the star. Most of the mass transfer occurs in the first two orbital periods and peak accretion rates reach values between 0.5 M/ms and 1.2M/ms (see Figure 3). The mass accretion rate then drops and the disk becomes more and more azimuthally symmetric, reaching a quasi–steady state by the end of the simulations.

We show in Figure 4 the various energies of the system (kinetic, internal, gravitational potential and total) for run D.

The dramatic drop in total internal energy reflects the intense mass accretion that takes place within the first couple of orbits. Figure 4 also shows \[panel (b)\] the total angular momentum of the system for runs A, B, D and E (the only contribution to the total angular momentum not plotted is the spin angular momentum of the black hole, see below). Angular momentum decreases for two reasons. First, if gravitational radiation reaction is still acting on the system, it will decrease approximately according to the quadrupole formula. Second, whenever matter is accreted by the black hole, the corresponding angular momentum is removed from our system. In reality, the angular momentum of the accreted fluid would increase the spin of the black hole. We keep track of this accreted angular momentum and exhibit its value in Table 3 as the Kerr parameter of the black hole. This shows up as a decrease in the total value of J.

### 3.4 Accretion disk structure

In Table 2 we show several parameters pertaining to the final accretion structure around the black hole for every run. The mass that has been accreted by the black hole is denoted by M<sub>acc</sub>. The disk settles down to a fairly azimuthally symmetric structure within a few initial orbital periods (except for the long tidal tail, which always persists as a well–defined structure), and there is a baryon–free axis above and below the black hole in every case (see below). We have calculated the mass of the remnant disk, $`M_{disk}`$, by searching for the amount of matter that has sufficient specific angular momentum j at the end of the simulation to remain in orbit around the black hole (as in Ruffert & Janka 1998 ). This material has $`j>j_{crit}=\sqrt{6}GM_t/c`$, where $`M_t`$ is the total mass of the system. By the end of the simulations, between 70% and 80% of the neutron star has been accreted by the black hole. It is interesting to note that the final accretion rate (at $`t=t_f`$) appears to be rather insensitive to the initial mass ratio, and is between 2.4 M s<sup>-1</sup> and 6.1 M s<sup>-1</sup>. From this final accretion rate we have estimated a typical timescale for the evolution of the accretion disk, $`\tau _{disk}=M_{disk}/\dot{M}_{final}`$. Despite the difference in the initial mass ratios and the typical sizes of the disks, the similar disk masses and final accretion rates make the lifetimes comparable for every run.

We have plotted azimuthally averaged density and internal energy profiles in Figure 5 for run D. The specific internal energy is greater towards the center of the disk, and flattens out at a distance from the black hole roughly corresponding the density maximum, at $`u3\times 10^{18}`$ erg g<sup>-1</sup>, or 3.1 MeV/nucleon, and is largely independent of the initial mass ratio. The inner regions of the disks have specific internal energies that are greater by approximately one order of magnitude.

Additionally, panel (b) in the same figure shows the azimuthally averaged distribution of specific angular momentum j in the orbital plane for all runs. The curves terminate at $`r_{in}=2r_{Sch}`$. Pressure support in the inner regions of the accretion disks makes the rotation curves sub–Keplerian, while the flattening of distribution marks the outer edge of the disk and the presence of the long tidal tail (see Figure 2),which has practically constant specific angular momentum.

The Kerr parameter of the black hole, given by $`a=J_{\mathrm{BH}}c/GM_{\mathrm{BH}}^2`$, is also shown in Table 3. We have calculated it from the amount of angular momentum lost by the fluid via accretion onto the black hole (see Figure 4b), assuming that the black hole is not rotating at $`t=0`$. The specific angular momentum of the black hole is smaller for lower mass ratios simply because the black hole is initially more massive when q is smaller.

It is of crucial importance for the production of GRBs from such a coalescence event that there be a baryon–free axis in the system along which a fireball may expand with ultrarelativistic velocities (Mészáros & Rees 1992, 1993 ). We have calculated the baryon contamination for every run as a function of the half–angle $`\mathrm{\Delta }\theta `$ of a cone directly above the black hole and along the rotation axis of the binary that contains a given amount of mass $`\mathrm{\Delta }M`$. Table 3 shows these angles (in degrees) for $`\mathrm{\Delta }M/\mathrm{M}_{}=1.4\times 10^3,1.4\times 10^4,1.4\times 10^5`$. There is a greater amount of pollution for high mass ratios (the disk is geometrically thicker compared to the size of the black hole), but in all cases only modest angles of collimation are required to avoid contamination. We note here that the values for $`\theta _5`$ are rough estimates at this stage since they are at the limit of our numerical resolution in the region directly above the black hole.

## 4 Discussion

The numerical simulations reported here were quasi-newtonian. An important caveat to keep in mind is that inclusion of general relativistic effects may lead to results qualitatively different from even post-newtonian treatment (Wilson, Mathews & Marronetti 1996 ).

Our results indicate that the outcome of the binary coalescence depends on the nature of the star orbiting the black hole. As reported previously (Lee & Kluźniak 1995 , 1998 , Kluźniak & Lee 1998 ), when we modeled the star as a polytrope with adiabatic index $`\mathrm{\Gamma }=3`$, the coalescence appeared to be an intermittent process in which the core of the polytrope survives the initial encounters and increases its separation from the black hole, thus extending the merger to possibly $`0.1`$ s. For the softer polytrope discussed here ($`\mathrm{\Gamma }=5/3`$), the star is disrupted completely in a few milliseconds and all that remains after the initial mass transfer is an accretion disk, containing no more than $`1/5`$ of the initial mass, and some ejecta. Perhaps the current simulation with $`\mathrm{\Gamma }=5/3`$, is the more realistic one, because for this polytrope $`dM/dR<0`$, as for physical models of neutron stars (e.g. Arnett & Bowers 1977 ).

In agreement with earlier suggestions (Lattimer & Schramm 1974, 1976 ), we have found that some matter will be ejected from the system, in an amount sufficient to account for the abundance of the r–process nuclei (assuming the r–process does indeed occur during the merger).

The binary coalescence of a neutron star with a black hole remains an attractive theoretical source of gamma–ray bursts. The energy requirements for at least one recently observed burst are so severe, if emission is isotropic (e.g. Kulkarni et al. 1998 ), that some degree of beaming seems desirable. According to the simulations presented here, ultrarelativistic flows are possible in the post-merger system only along the rotational axis of the system in a solid angle of about $`0.1`$ steradian. Proper inclusion of neutrino transport may change this angle somewhat. The rather short ($`50`$ ms) accretion timescale of the remnant disk reported, does not include possible interaction between the disk and the black hole (Blandford & Znajek 1977 ). In fact, the appearance in the simulation of a substantial toroidal disk around the black hole is encouraging, as it may allow the black hole spin to be extracted by the Blandford-Znajek mechanism, possibly powering in this manner the gamma–ray burst fireball (Mészáros & Rees 1997 ).

## ACKNOWLEDGMENTS

We gratefully acknowledge support for this work from DGAPA–UNAM and KBN (grant P03D01311).

## References