# Newtonian hydrodynamics of the coalescence of black holes with neutron stars II: Tidally locked binaries with a soft equation of state.

## 1 Introduction

Angular momentum losses to gravitational radiation are expected to lead to the coalescence of binary systems containing black holes, and/or neutron stars (when the initial binary separation is small enough for the decay to take place in less than the Hubble time). This type of evolution has been suggested in a variety of contexts as possibly giving rise to observable events, such as gamma—ray bursts (GRBs) and bursts of gravitational waves (see e.g. Thorne 1995). Additionally, it could help explain the observed abundances of heavy elements in our galaxy (Lattimer & Schramm 1974; 1976) if the star is tidally disrupted in the encounter (see Wheeler 1971). Study of such events could also provide constraints on the equation of state at supra–nuclear densities.

After the recent measurement of redshifts to their afterglows \[Metzger et al. 1997, Kulkarni et al. 1998, Djorgovski et al. 1998\], it is now generally believed that GRBs originate at cosmological distances. The calculated event rates \[Lattimer & Schramm 1976, Narayan, Piran & Shemi 1991, Tutukov & Yungelson 1993, Lipunov, Postnov & Prokhorov 1997, Portegies Zwart & Yungelson 1998, Bełczyński & Bulik 1999\] for merging compact binaries are compatible with the observed frequency of GRBs (on the order of one per day). The preferred model for the production of a GRB invokes a relativistic fireball from a compact ‘central engine’ that would produce the observable $`\gamma `$–rays through internal shocks (Mészáros & Rees 1992; 1993). This model requires the presence of a relatively baryon–free line of sight from the central engine to the observer along which the fireball can expand at ultrarelativistic speeds. Additionally, the short–timescale variations seen in many bursts (often in the millisecond range) probably arise within the central engine \[Sari & Piran 1997\].

The coalescence of binary neutron star systems or black hole–neutron star binaries was suggested as a mechanism capable of powering the gamma–ray bursts, either during the binary merger itself or through the formation of a dense accretion disk which could survive long enough to accomodate the variable timescales of GRBs (Paczyński 1986; Goodman 1986; Goodman, Dar & Nussinov 1987; Eichler et al. 1988; Paczyński 1991; Mészáros & Rees 1992; Woosley 1993; Jaroszyński 1993; 1996; Witt et al. 1994; Wilson, Mathews & Marronetti 1996; Ruffert, Janka & Schäfer 1996; Lee & Kluźniak 1997; Ruffert, Janka, Takahashi & Shäfer 1997; Katz 1997; Kluźniak & Lee 1998; Ruffert & Janka 1998; Popham, Woosley & Fryer 1998; McFadden & Woosley 1998; Ruffert & Janka 1999). The enormous amount of gravitational energy that would be liberated in such an event could account for the energetics of the observed GRBs and neutrino–antineutrino annihilation may power the necessary relativistic fireball.

In previous work (Lee & Kluźniak 1995; 1998 (hereafter Paper I); Kluźniak & Lee 1998), we have studied the coalescence of a neutron star with a stellar–mass black hole for a stiff ($`\mathrm{\Gamma }=3`$) polytropic equation of state and a range of mass ratios. We found that the neutron star was not entirely disrupted, but rather remained in orbit (with a greatly reduced mass) about the black hole after a quick episode of mass transfer. Thus the duration of the coalescence process would be extended from a few milliseconds to possibly several tens of milliseconds. The observed outcome seemed favorable for the production of a GRB since in every case we found a baryon–free axis in the system, along the axis of rotation.

In the present paper, we investigate the coalesence of a black hole–neutron star binary for a soft equation of state (with an adiabatic index $`\mathrm{\Gamma }=5/3`$) and a range of mass ratios. Our initial conditions are as in Paper I in that they correspond to tidally locked binaries. Complete tidal locking is not realistically expected \[Bildsten & Cutler 1992\], but it can be considered as an extreme case of angular momentum distribution in the system. In the future we will explore configurations with varying degrees of tidal locking.

As before, the questions motivating our study are: Is the neutron star tidally disrupted by the black hole and does an accretion torus form around the black hole? If so, how long–lived is it? Is the baryon contamination low enough to allow the formation of a relativistic fireball? Is any significant amount of mass dynamically ejected from the system? What is the gravitational radiation signal like, and how does it depend on the equation of state and the initial mass ratio?

In section 2 we present the method we have used to carry out our simulations. This is followed by a presentation of our results in section 3 and a discussion in section 4.

## 2 Numerical method

For the simulations presented in this paper, we have used the method known as Smooth Particle Hydrodynamics (SPH). Our code is three–dimensional and essentially Newtonian. This method has been described often, we refer the reader to Monaghan \[Monaghan 1992\] for a review of the principles of SPH, and to Paper I and Lee \[Lee 1998\] for a detailed description of our own code, including the tree structure used to compute the gravitational field.

We model the neutron star via a polytropic equation of state, $`P=K\rho ^\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }=5/3`$. For the following, we measure distance and mass in units of the radius R and mass $`M_{\mathrm{NS}}`$ of the unperturbed (spherical) neutron star (13.4 km and 1.4 M respectively), except where noted, so that the units of time, density and velocity are

$`\stackrel{~}{t}=1.146\times 10^4\mathrm{s}\times \left({\displaystyle \frac{R}{13.4\text{km}}}\right)^{3/2}\left({\displaystyle \frac{M_{\mathrm{NS}}}{1.4M_{}}}\right)^{1/2}`$ (1)   
$`\stackrel{~}{\rho }=1.14\times 10^{18}\mathrm{kg}\mathrm{m}^3\times \left({\displaystyle \frac{R}{13.4\text{km}}}\right)^3\left({\displaystyle \frac{M_{\mathrm{NS}}}{1.4M_{}}}\right)`$ (2)   
$`\stackrel{~}{v}=0.39c\times \left({\displaystyle \frac{R}{13.4\text{km}}}\right)^{1/2}\left({\displaystyle \frac{M_{\mathrm{NS}}}{1.4M_{}}}\right)^{1/2}.`$ (3)   
and we use corresponding units for derivative quantities such as energy and angular momentum.

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; Paper I). The binary separation is defined henceforth as the distance between the black hole and the center of mass of the SPH particles. During the construction of these configurations, the specific entropies of all particles are maintained constant, i.e. $`K`$=constant in $`P`$=$`K\rho ^\mathrm{\Gamma }`$. The neutron star is modeled with $`N=17,256`$ particles in every case presented in this paper. To ensure uniform spatial resolution, the masses of the particles were made proportional to the Lane–Emden densities on the initial grid.

To carry out a dynamical run, the black hole and every particle are given the azimuthal velocity corresponding to the equilibrium value of the angular frequency $`\mathrm{\Omega }`$ in an inertial frame, with the origin of coordinates at the center of mass of the system. Each SPH particle is assigned a specific internal enery $`u_i=K\rho ^{(\mathrm{\Gamma }1)}/(\mathrm{\Gamma }1)`$, and the equation of state is changed to that of an ideal gas, where $`P=(\mathrm{\Gamma }1)\rho u`$. The specific internal energy of each particle is then evolved according to the first law of thermodynamics, taking into account the contributions from the artificial viscosity present in SPH. During the dynamical runs we calculate the gravitational radiation waveforms in the quadrupole approximation.

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. Using the quadrupole approximation, the rate of energy change for a point–mass binary is given by (see Landau & Lifshitz 1975):

$`{\displaystyle \frac{dE}{dt}}={\displaystyle \frac{32}{5}}{\displaystyle \frac{G^4(M_{\mathrm{NS}}+M_{\mathrm{BH}})(M_{\mathrm{NS}}M_{\mathrm{BH}})^2}{(cr)^5}}`$ (4)   
and the rate of angular momentum loss by

$`{\displaystyle \frac{dJ}{dt}}={\displaystyle \frac{32}{5c^5}}{\displaystyle \frac{G^{7/2}}{r^{7/2}}}M_{\mathrm{BH}}^2M_{\mathrm{NS}}^2\sqrt{M_{\mathrm{BH}}+M_{\mathrm{NS}}}.`$ (5)   
From these equations a radiation reaction acceleration for each component of the binary can be obtained as

$`𝒂^{}={\displaystyle \frac{1}{q(M_{\mathrm{NS}}+M_{\mathrm{BH}})}}{\displaystyle \frac{dE}{dt}}{\displaystyle \frac{𝒗^{}}{(v^{})^2}}`$ (6)   
$`𝒂^{\mathrm{BH}}={\displaystyle \frac{q}{M_{\mathrm{NS}}+M_{\mathrm{BH}}}}{\displaystyle \frac{dE}{dt}}{\displaystyle \frac{𝒗^{\mathrm{BH}}}{(v^{\mathrm{BH}})^2}}`$ (7)   
where $`v^{}`$ is the velocity of the neutron star and $`v^{\mathrm{BH}}`$ that of the black hole.

We have used this formula for our calculations to simulate the effect of gravitational radiation reaction on the system. Clearly, the application of equation (7) to the black hole in our calculations is trivial, since we always treat it as a point mass. For the neutron star, we have chosen to apply the same acceleration to all SPH particles. This value is that of the acceleration at the center of mass of the SPH particles, so that equation (6) now reads:

$`𝒂^i={\displaystyle \frac{1}{q(M_{\mathrm{NS}}+M_{\mathrm{BH}})}}{\displaystyle \frac{dE}{dt}}{\displaystyle \frac{𝒗_{cm}^{}}{(v_{cm}^{})^2}},`$ (8)   
This formulation of the gravitational radiation reaction has been used in SPH simulations by others \[Davies et al. 1994, Zhuge, Centrella & McMillan 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. 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: the radiation reaction is turned off if the center of mass of the SPH particles comes within a prescribed distance of the black hole (effectively a tidal disruption radius). 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

We now describe our results. First, we present the initial conditions that were used to perform the dynamical runs. We then describe the general morphology of the coalescence events, the detailed structure of the accretion disks that form as a result of the tidal disruption of the neutron star, and the gravitational radiation signal.

### 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 (Paper I), 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 1

we show the variation of total angular momentum J in these sequences as a function of binary separation for the four values of the mass ratio (solid lines). Following Lai, Rasio & Shapiro \[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 for the $`\mathrm{\Gamma }=5/3`$ polytrope 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`$ (Paper I), the variation in angular momentum as a function of binary separation was qualitatively different (for high mass ratios) from our present findings. 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.

There is a crucial difference between the two polytropes considered in Paper I and here. For polytropes, the mass–radius relationship is $`RM^{(\mathrm{\Gamma }2)/(3\mathrm{\Gamma }4)}`$. For $`\mathrm{\Gamma }`$=3 this becomes $`RM^{1/5}`$, while for $`\mathrm{\Gamma }`$=5/3, $`RM^{1/3}`$. Thus, the polytrope considered in Paper I responded to mass loss by shrinking. The $`\mathrm{\Gamma }=5/3`$ 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 Paper I, 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 equations (4) and (5) one can compute the binary separation as a function of time for a point–mass binary in the quadrupole approximation, and obtain

$`r=r_i\left(1t/t_0\right)^{1/4},`$ (9)   
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 in this paper, 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.7 we find $`t_0=56.81\times \stackrel{~}{t}`$=6.5 ms) and $`P=19.58\times \stackrel{~}{t}`$=2.24 ms), so one must analyze whether hydrodynamical effects will drive the coalescence on a comparable timescale. We have performed a dynamical simulation for $`q`$=1, with an initial separation on the verge of Roche Lobe overflow, at $`r`$=2.7, without including radiation reaction in the equations of motion. We show in Figure 2a

the binary separation as a function of time for this calculation (solid line) as well as the separation for a point–mass binary decaying in the quadrupole approximation, using equation (9) (dashed line). In the dynamical simulation the orbital separation remains approximately constant, and begins to decay rapidly around $`t`$=110 (in the units defined in equation ), when mass loss from the neutron star becomes important. Clearly at this stage hydrodynamical effects are dominant, but one must include radiation reaction in the early stages of the process. There is an added (practical) benefit derived from including radiation reaction in these calculations. As seen in Figure 2a, it takes a full 15 ms for the orbit to become unstable. Simulating the behavior of the system at high resolution (practically no SPH particles have been accreted at this stage) for such a long time is computationally expensive, whereas accretion in the early stages of the simulation allows us to perform, in general, more calculations at higher resolution. We have thus included radiation reaction in all the runs (A through E) presented in this paper, and adopted a switch as described in section 2.

### 3.2 Run parameters

In Table 1 we present the parameters distinguishing each dynamical run we performed. All times are in units of $`\stackrel{~}{t}`$ (equation ) and all distances in units of R, the unperturbed (spherical) stellar radius. 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}=200`$, equivalent to 22.9 ms (this covers on the order of ten initial orbital periods for the mass ratios considered).

The fifth column in Table 1 shows the value of $`t_{rad}`$, when radiation reaction is switched off according to the criterion described in section 2. In Figure 3 we show density contours in the orbital plane for runs A, C, D and E at times very close to $`t=t_{rad}`$. The corresponding plot for run B is very similar to that for run A. Note that runs C and D differ only in the corresponding value of $`t_{rad}`$. For run D there is little doubt that approximating the neutron star as a point–mass is still reasonable at this stage, while for run C this is clearly not the case. We can then use these two runs to gauge the effect of our simple radiation reaction formulation on the outcome of the coalescence event. We note here that run E is probably beyond 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=5.05`$, equivalent to 67.87 km) is such that the neutron star is within the innermost stable circular orbit of a test particle around a Schwarzschild black hole of the same mass ($`r_{ms}=9.17`$, equivalent to 123.26 km). Thus we present in Appendix A a dynamical run with initial mass ratio $`q=0.1`$ making use of a pseudo–Newtonian potential for the black hole. For the following, we will at times omit a discussion of run B, as it is qualitatively and quantitatively very similar to run A (both of these have a relatively high mass ratio).

### 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, A through E. Once accretion begins, the total number of particles decreases. Since this compromises resolution, we made a modification to the code for run E to avoid the number of particles from dropping below 9,000. This is done simply by splitting a given fraction of the particles $`N_{split}`$ and creating $`2N_{split}`$ particles from them. Total mass and momentum are conserved during this procedure, and it can be shown that the numerical noise introduced into the smoothed density $`\rho `$ by doing this is of the order of the accuracy of the SPH method itself, $`𝒪`$($`h^2`$), where $`h`$ is the smoothing length \[Meglicki, Wickramasinghe & Bicknell 1993\].

In every run the binary separation (solid lines in in Figure 2b) initially decreases due to gravitational radiation reaction. For high mass ratios, (runs A, B) the separation decays faster that 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 (Paper I) 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 runs C and D, the solid and dashed lines in Figure 2b 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 2 are computed for a fixed mass ratio). From the expression for the timescale for orbital decay $`t_0`$ in equation (9), it is apparent that at constant total mass, lowering the mass ratio slows the orbital decay, when $`q<0.5`$.

The general behavior of the system is qualitatively similar for every run. Figures 45 and 6 show density contours in the orbital plane (left columns) and in the meridional plane containing the black hole (right columns) for runs A, D and E respectively, at $`t=50`$ and $`t=t_f=200`$ (equivalent to 5.73 ms and 22.9 ms). The corresponding plots for runs B and C are very similar to those for runs A and D, respectively. 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.04 and 0.1—equivalent to 0.49 and 1.22 M/ms (see Figure 7). 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 8 the various energies of the system (kinetic, internal, gravitational potential and total) for runs A, D and E.

The dramatic drop in total internal energy reflects the intense mass accretion that takes place within the first couple of orbits. Figure 8 also shows \[in panel (d)\] the total angular momentum of the system (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 equation (5). 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 2 as the Kerr parameter of the black hole. This shows up as a decrease in the total value of J. It is clear from runs C and D that the value of $`t_{rad}`$ influences the peak accretion rate and the mass of the black hole (particularly immediately after the first episode of heavy mass transfer). The maximum accretion rate is different for run C and D by about a factor of 1.4. This is easy to understand, since radiation reaction increases angular momentum losses, and hence more matter is accreted per unit time when it is present.

### 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 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 (Figure 9). 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 1999). This material has $`j>j_{crit}=\sqrt{6}GM_t/c`$, where $`M_t`$ is the total mass of the system. The values shown in Table 2 are equivalent to a few tenths of a solar mass, and again the effect of $`t_{rad}`$ can be seen by comparing runs C and D, where the disk masses differ by a factor of 1.2. 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\times 10^4`$ and $`5\times 10^4`$ (equivalent to 2.4 and 6.1 M s<sup>-1</sup> respectively). 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}`$—for reference, $`\tau =100`$, in the units of equation (1) corresponds to $`11.5`$ ms and thus the values of $`\tau _{disk}`$ given in Table 2 are between 47 and 63 ms. Despite the difference in the initial mass ratios and the typical sizes of the disks ($`r_0`$ is the radial distance from the black hole to the density maximum at $`t=t_f`$), 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 10 for runs A, D and E. 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 $`u2\times 10^2`$. This value corresponds to $`2.74\times 10^{18}`$ erg g<sup>-1</sup> or 2.9 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 (d) 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 11),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 2. We have calculated it from the amount of angular momentum lost via accretion onto the black hole (see Figure 8d), assuming that the black hole is not rotating at $`t=0`$. The final 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. The difference in the value of a for runs C and D is important (almost a factor of 2), and again reflects the influence of gravitational radiation reaction (for a larger value of $`t_{rad}`$ the black hole is spun up to a greater degree because of the larger amount of accreted mass).

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 (and below) the black hole and along the rotation axis of the binary that contains a given amount of mass $`\mathrm{\Delta }M`$. Table 2 shows these angles (in degrees) for $`\mathrm{\Delta }M=10^3,10^4,10^5`$ (equivalent to $`1.4\times 10^3,1.4\times 10^4,1.4\times 10^5`$ M respectively). 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. This can be seen by inspection in Figure 9

where we show the enclosed mass as a function of half–angle $`\mathrm{\Delta }\theta `$ for all runs at $`t=t_f`$.

### 3.5 Ejected mass and r–process

To calculate the amount of dynamically ejected mass during the coalescence process, we look for matter that has a positive total energy (kinetic+gravitational potential+internal) at the end of each simulation. Figure 11 shows a large–scale view of the system at $`t=t_f`$ for runs D and E. The thick black line running across the tidal tail in each case divides matter that is bound to the black hole from that which may be on outbound trajectories. This matter comes from the part of the neutron star that was initially furthest from the black hole and was ejected through the outer Lagrange point in the very early stages of mass transfer. We find that a mass between $`4.4\times 10^3`$ M and $`2.2\times 10^2`$ M can potentially be ejected in this fashion (see Table 2). This is very similar to what has recently been calculated for binary neutron star mergers for a variety of initial configurations \[Rosswog et al. 1999\]. Since it is believed that the event rate for binary neutron star mergers is comparable to that of black hole–neutron star mergers, this could prove to be a sizable contribution to the amount of observed r–process material. This result appears to be strongly dependent on the equation of state, since we previously observed no significant amount of matter being ejected for a system with a stiff equation of state (Paper I). For binary neutron star mergers, Rosswog et al. 1999 have obtained the same qualitative result. We also note that there is a significant difference in the amount of ejected mass for runs C and D (approximately a factor of four), due to the difference in the values of $`t_{rad}`$ (in run D the system loses less angular momentum and thus matter escapes with greater ease).

### 3.6 Gravitational radiation waveforms and luminosities

The emission of gravitational radiation is calculated in all our models in the quadrupole approximation (see e.g. Finn 1989; Rasio & Shapiro 1992), and can be obtained directly from the hydrodynamical variables of the system. The calculation of the gravitational radiation luminosity then requires only an additional numerical differentiation. Figure 12 shows the computed waveforms and luminosities, along with what the waveforms would be for a point–mass binary, also calculated in the quadrupole approximation. It is apparent that hydrodynamical effects play an important role particularly for high mass ratios in the early stages of the coalescence process (see panel (a) in Figure 12). When the neutron star is tidally disrupted, the amplitude of the waveform drops abruptly and practically to zero, since a structure that is almost azimuthally symmetric has formed around the black hole. This is in stark contrast to what occurred for a stiff equation of state (Lee 1998; Paper I; Kluźniak & Lee 1998), when the binary system survived the initial episode of mass transfer and a stable binary was the final outcome. In fact, these waveforms resemble more the case of a double neutron star merger with a soft equation of state \[Rasio & Shapiro 1992\], in which the coalescence resulted in a compact, azimuthally symmetric object surrounded by a dense halo and spiral arms. Table 3 shows the maximum amplitude $`h_{max}`$ for an observer located a distance $`r_0`$ away from the system along the axis of rotation, the maximum luminosity $`L_{max}`$ and the total enery $`\mathrm{\Delta }E_{GW}`$ emitted during the event. This last number should be taken only as an order of magnitude estimate since it depends on the choice of the origin of time. The peak luminosities are $`(R/M_{\mathrm{NS}})^5(L_{max}/L_0)=0.37`$ for run A, $`(R/M_{\mathrm{NS}})^5(L_{max}/L_0)=1.50`$ for run D and $`(R/M_{\mathrm{NS}})^5(L_{max}/L_0)=5.90`$ for run E (equivalent to $`1.12\times 10^{55}`$erg s<sup>-1</sup>, $`4.55\times 10^{55}`$erg s<sup>-1</sup> and $`1.79\times 10^{56}`$erg s<sup>-1</sup> respectively). We note that although the waveforms for runs C (not plotted) and D (panel (b) in Figure 12) are very similar, the maximum amplitudes and the luminosity differ by about 1.3% and 3.4% respectively, a small but non–negligible amount. This is again a reflection of the way in which the radiation reaction was formulated, and indicates that a more rigorous treatment of this effect is necessary.

## 4 Summary and Discussion

We have presented results of hydrodynamical simulations of the binary coalescence of a black hole with a neutron star. We have used a polytropic equation of state (with index $`\mathrm{\Gamma }=5/3`$) to model the neutron star, and a Newtonian point mass with an absorbing surface at the Schwarzschild radius to represent the black hole. All our computations are strictly Newtonian, but we have included a term that approximates the effect of gravitational radiation reaction in the system. We have also calculated the emission of gravitational radiation in the quadrupole approximation.

We have found that for every mass ratio investigated ($`M_{\mathrm{NS}}/M_{\mathrm{BH}}`$=1, 0.8, 0.31 and 0.1) the $`\mathrm{\Gamma }=5/3`$ polytrope (‘neutron star’) is entirely disrupted by tidal forces, and a dense accretion torus, containing a few tenths of a solar mass, forms around the black hole. The maximum densities and specific internal energies in the tori are on the order of $`10^{11}\text{g cm}\text{-3}`$ and $`10^{19}\text{erg g}\text{-1}`$ (or 10 MeV/nucleon) respectively (all simulations were run for approximately 22.9 ms). The final accretion rate is between 2 and 6 solar masses per second, and hence the expected lifetime of the torus $`\tau _{disk}=M_{disk}/\dot{M}_{final}`$ is between 40 and 60 milliseconds.

The rotation axis of the system remains free of matter to a degree that would not hinder the production of a relativistic fireball, possibly giving rise to a gamma ray burst. Although the duration of such a bursts would still be too short to power the longest GRBs, the present scenario could well account for the subclass of short bursts \[Kouveliotou et al. 1995\]. A significant amount of matter (between $`10^2`$ and $`10^3`$ solar masses) is dynamically ejected from the system, and could contribute significantly to the observed abundances of r–process material in our galaxy. The gravitational radiation signal is very similar to that of a point–mass binary until the beginning of mass transfer, particularly for low mass ratios. After mass transfer starts, the amplitude of the waveforms drops dramatically on a dynamical timescale when the accretion torus is formed. In every aspect, the results are dramatically different from what occurs for a stiff equation of state.

## Acknowledgments

We gratefully acknowledge financial support from DGAPA–UNAM and KBN grant 2P03D01311. W.L. thanks Craig Markwardt for helpful discussions concerning the effect of a soft equation of state on the system. It is a pleasure to thank the referee for his most helpful comments.

## Appendix A Dynamical evolution for a pseudo–Newtonian potential

Since, as stated in section 3.2, the system with mass ratio $`q=0.1`$ (run E) is at an initial separation that is within the marginally stable orbit for test particles around a Schwarzschild black hole, we have performed an additional run, altering the form of the potential produced by the black hole. We will perform a detailed description of results elsewhere (Lee 1999). Here we only present a comparison to the results of run E. We have chosen the form proposed by Paczyński & Wiita (1980), namely:

$`\mathrm{\Phi }_{\mathrm{BH}}^{PW}(r)=GM_{\mathrm{BH}}/(rr_{Sch}).`$ (10)   
This potential correctly reproduces the positions of the marginally bound and marginally stable orbits for test particles. A few modifications need to be made to the SPH code to accomodate this new potential. First, the absorbing boundary of the black hole is now placed at a distance $`r_{boundary}=1.5r_{Sch}`$ from the position of the black hole; second, the total gravitational force exerted by the neutron star on the black hole is symmetrized so that total linear momentum is conserved in the system.

A tidally locked equilibrium configuration for a given separation $`r`$ can be constructed with these modifications in the same manner as described in section 3.1. For a test particle in orbit about a Schwarzschild black hole, the marginally stable orbit appears at a separation $`r_{ms}=6GM_{\mathrm{BH}}/c^2`$ because the total angular momentum exhibits a minimum at that point. In our case (where the neutron star has a finite size and mass), the turning point in the curve of total angular momentum as a function of binary separation occurs at approximately $`r=9.1R_{\mathrm{NS}}`$, and so we have chosen this value for the initial separation $`r_i`$ to be used in the dynamical simulation. Gravitational radiation reaction has been implemented as described in section 2, with a slight modification to the definition of $`r_{tidal}`$ to account for the increased strength of the gravitational interactions, so that now $`r_{tidal}=CR(M_{\mathrm{BH}}/M_{\mathrm{NS}})^{1/3}+r_{Sch}`$.

Since gravitational interactions are stronger with the modified form of the gravitational potential, the overall encounter is more violent. The neutron star is tidally disrupted into a long tidal tail in a way similar to that exhibited in run E. The accretion episode is very brief, with a peak accretion rate onto the black hole of $`\dot{M}_{max}=0.7`$, equivalent to 8.5 M/ms and thus substantially higher than that for run E (see Table 2).

We followed the dynamical evolution from $`t=0`$ to $`t_f=200`$, and show final density contour plots in the orbital and meridional plane containing the black hole in Figure 13. By the end of the simulation, the fluid has not formed a quasi–static accretion structure around the black hole as for the Newtonian runs, and 99.2% of the initial neutron star mass has been accreted ($`M_{acc}=0.992`$). The thick black line in Figure 13a divides material that is bound to the black hole from that which is on outbound trajectories (see Figure 11 for a comparison with runs D and E). Overall, a smaller amount of mass is left over after the initial episode of heavy mass transfer (approximately an order of magnitude less than for run E), but a larger fraction ($`M_{ejected}=6.8\times 10^3`$, equivalent to $`9.6\times 10^3`$M) may be dynamically ejected from the system. The region above and below the black hole is devoid of matter to an even greater extent than for the Newtonian case as can be seen in Figure 14,

where we plot the enclosed mass $`\mathrm{\Delta }M`$ as a function of the half angle $`\mathrm{\Delta }\theta `$ of a cone directly above (and below) the black hole at $`t=t_f`$, as in Figure 9. Thus in this scenario the production of a relativistic fireball that could give rise to a gamma–ray burst would require an even more modest degree of beaming in order to avoid baryon contamination.