text:
# Critical Behaviour in the Spin Fluctuations and Superfluid Density of La2-xSrxCuO4

## Abstract

We study the doping dependence of low frequency spin fluctuations and the zero-temperature superfluid density of La2-xSrxCuO4 using the muon spin relaxation ($`\mu `$SR) and ac-susceptibility techniques. Superconductivity is found to coexist with low frequency spin fluctuations over a large region of the superconducting phase diagram. The characteristic temperature of spin fluctuations detected by $`\mu `$SR decreases with increasing $`x`$ and vanishes above a critical doping $`x_c0.19`$. This value of $`x_c`$ coincides with the doping at which the normal state pseudogap extrapolates to zero. The superfluid density behaves in the opposite way to the low frequency spin fluctuations. It increases with $`x`$ and becomes nearly doping-independent for $`x>x_c`$. These results are consistent with predictions involving quantum criticality at $`x_c`$.

The importance of quantum criticality and the interplay between magnetism and superconductivity in elucidating the pairing mechanism in high-$`T_c`$ superconductors (HTS) have been discussed in many theoretical and experimental works since their discovery . Spin fluctuations and quantum criticality have been extensively studied in heavy fermion compounds, where the appearance of superconductivity as the Neel temperature is suppressed with pressure, has provided strong evidence both for a quantum critical point and for a pairing mechanism mediated by spin fluctuations . Unlike heavy fermion compounds the properties of HTS are usually tuned by varying the carrier concentration (doping). Fundamental physical quantities such as the superconducting condensation energy  and the superfluid density  show distinct changes at a critical doping $`x_c`$ which is slightly above the optimal carrier concentration defined as the point at which the superconducting transition temperature, $`T_c`$, is maximum. Recently, Tallon and Loram  have analysed the doping dependence of the normal state pseudogap from the specific heat and a set of other experimental data and suggested that the normal state pseudogap vanishes abruptly at $`x_c`$, as predicted for quantum criticality . It is therefore important to search for a critical behaviour in HTS using direct measurements of the spin or charge fluctuations.

In this Letter we present a systematic study of the doping dependence of low frequency spin fluctuations and the superfluid density of La2-xSrxCuO4. From the muon depolarisation and the in-plane penetration depth determined both by $`\mu `$SR and ac-susceptibility, we identify a critical doping $`x_c`$ at which qualitative changes occur in both the spin excitations and the superfluid density. For $`x<x_c`$ there are significant low frequency magnetic fluctuations whereas the superfluid density changes in the opposite manner. Above $`x_c`$ these fluctuations are not detected and the superfluid density becomes nearly doping independent as expected for conventional BCS superconductors. The correlation between these two quantities provides evidence for the presence of a quantum critical point, in support of theoretical predictions made by several groups .

Zero-field $`\mu `$SR is a sensitive modern technique for studying low frequency spin dynamics. It is a local probe and can detect very small (1$`G`$) internal magnetic fields . The data obtained by $`\mu `$SR give a direct measure of the low energy spin dynamics in the sample . From the time dependence of the depolarisation of muons, one can obtain important information about the characteristic energy or temperature scale of spin fluctuations and the spin-glass transition temperature.

La2-xSrxCuO4 (LSCO) is one of the simplest HTS systems. It has a simple crystal structure and its carrier concentration can be accurately controlled. Undoped La2CuO4 is an insulator with long range antiferromagnetic order. Substitution of Sr2+ for La3+ introduces holes into the CuO2 planes and at $`x0.02`$ the system exhibits a short range ordered, antiferromagnetic or spin glass phase . Superconductivity emerges near $`x=0.05`$ and the superconducting transition temperature approximately follows a parabolic doping dependence and vanishes at $`x0.30`$ .

The samples we studied were high quality polycrystalline La2-xSrxCuO4 ($`x=0.030.24`$) synthesised in Cambridge. The phase purity of these samples, as characterised by powder x-ray diffraction and micro-Raman spectroscopy  as well as extensive transport  and thermodynamic measurements  was found to be better than $`1\%`$. Their lattice parameters and $`T_c`$ values were also in good agreement with published data obtained on powders as well as single crystals .

Zero-field (ZF) and transverse-field (TF) $`\mu `$SR experiments were performed at the pulsed muon source, ISIS Facility, Rutherford Appleton Laboratory and spectra were collected over the temperature range from $`40mK`$ to $`150K`$. In a $`\mu `$SR experiment, $`100\%`$ spin polarised positive muons are implanted into a specimen and evolve in their local magnetic environment. The muon decays with a life time $`2.2\mu s`$, emitting a positron preferentially in the direction of the muon spin at the time of decay. By accumulating time histograms of such positrons one may deduce the muon depolarisation rate as a function of time after implantation. The pulsed muon facility at ISIS allows the muon polarisation function to be measured for up to 6-7 muon life times, and thus small changes in the depolarisation function can be detected. This property of the ISIS muon facility has allowed us to study the temperature dependence of slow spin fluctuations in a HTS. The muon is expected to reside in the most electronegative site of the lattice. In the case of LSCO this is 1.0Å away from the apical O-2 site of the CuO2 plane . As will be discussed later, this assignment is supported by the strong enhancement in the muon depolarisation rate with Zn doping since Zn is known to substitute in the CuO2 planes. We therefore believe that the results reported here are dominated by the magnetic correlations in the CuO2 planes.

The zero temperature superfluid density, $`\rho ^s(0)`$, was determined from measurements of the magnetic penetration depth $`\lambda _{ab}`$ for current flow in the CuO2 planes, ($`\rho ^s`$$``$$`\lambda _{ab}^2`$). $`\lambda _{ab}`$ was obtained for all samples using the low-field $`ac`$ -susceptibility technique (typically at $`1G`$ and $`333Hz`$) for grain-aligned powders . Details of the technique can be found in refs . For each doping concentration of LSCO, the result presented here is typical of those obtained for 2-4 samples prepared independently by aligning powders from the same polycrystalline pellet. Four pieces were cut out from each aligned sample and measured. In other words 8-16 samples were investigated for each doping content. The values of $`\lambda _{ab}^2(0)`$ for samples with $`x`$$`0.15`$ were also confirmed by standard TF-cooled $`\mu `$SR experiments performed on unaligned powders at 400 Gauss. In a TF-cooled $`\mu `$SR experiment, the field distribution of a flux-line lattice produced by an external field is probed by muons. The depolarisation rate of the initial muon spin is proportional to $`\lambda _{ab}^2(0)`$. For $`x<0.15`$, $`\lambda _{ab}^2(0)`$ obtained from TF-$`\mu `$SR is less accurate than that obtained from the $`ac`$-susceptibility measurements because of a strong enhancement in the muon depolarisation rate below $`0.3`$ $`T_c`$  probably due to the rapid increase in slow spin fluctuations with underdoping as discussed below.

Figure 1 shows the typical time dependence of the ZF muon asymmetry for $`x=0.08`$ ($`T_c=21K`$) for two characteristic temperatures. In all samples the high temperature form of the depolarisation is Gaussian, consistent with dipolar interactions between the muons and their near neighbour nuclear moments. This was verified by applying a $`50G`$ longitudinal field, which completely suppressed the depolarisation. The onset of dynamical relaxation processes becomes apparent at low temperatures by a deviation of the depolarisation from a Gaussian behaviour.

The depolarisation rate can be well described by a stretched exponential function $`G_z(t)=A_1exp((\gamma t)^\beta )+A_2`$, where $`\gamma `$ is the dynamic muon spin depolarisation rate and $`A_2`$ accounts for a small time independent background arising from muons stopping in the silver backing plate. As can be seen in Fig. 2 the exponent $`\beta `$ starts from the value 2 at high temperatures, corresponding to Gaussian relaxation, but decreases steadily at lower temperatures, suggesting the development of a low frequency component in the spectral weight of the spin dynamics. We take the temperature at which $`\beta `$ first drops below the value 2 as one indicator of the energy scale of the magnetic correlations in a given sample. This temperature is denoted $`T_f`$. At low temperatures (Fig. 2) the values of $`\beta `$ decrease towards the value 0.5. At lower temperatures still the form of the relaxation function changes (Fig. 1a): there is an initial rapid drop of asymmetry, followed by a long-time tail with a slower relaxation. This is very charactersitic of the behaviour found in spin glass samples below $`T_g`$ , where the initial rapid drop is ascribed to the static distribution of random local fields while the long-time behaviour results from dynamical processes. In this regime the data were fitted to the form $`G_z(t)=A_1exp(\gamma _1t)+A_2exp((\gamma _2t)^\beta )+A_3`$. We identify the spin glass temperature $`T_g`$ as the temperature at which the value of $`\beta `$ (Fig. 2) reaches the value 0.5. This ”root exponential” form for the relaxation function is a common feature of spin glasses , and in the present samples this temperature coincided with the cross-over in behaviour (as shown in Fig. 1) and with a peak in the longitudinal relaxation rate. All the samples with $`x=0.030.125`$ followed the same behaviour.

Our data indicate that the spin glass phase persists beyond $`x=0.125`$. In fact the onset of the spin glass phase for $`x=0.125`$ occurs at a higher temperature than that for $`x=0.10`$. This may be due to the formation of strongly correlated antiferromagnetic stripe domains in this range of doping . For $`x=0.15`$ and $`0.17`$, $`T_g`$ becomes very small ($`<45mK`$) and $`T_f`$ is approximately equal to $`8K`$ and $`2K`$, respectively. For $`x0.20`$, there are no changes in the depolarisation function to the lowest temperature measured ($`40mK`$). The presence of a finite $`T_g`$ for $`x>0.125`$ is verified by the Zn doping studies discussed below.

Figure 3 shows the doping dependences of $`T_g`$ and $`T_f`$ together with $`T_c`$. Although the freezing of spins occurs at very low temperatures, low frequency spin fluctuations appear at significantly higher temperatures. Both $`T_g`$ and $`T_f`$ are found to decrease with increasing doping and tend to zero at $`x_c0.19`$. Their behaviour resembles that of the normal state gap  determined from measurements above $`T_c`$, which also tends to vanish at the same $`x_c`$. This is better demonstrated in the semi-log plot of $`T_g`$ and $`T_f`$ shown in the inset of the figure. The similarity in the doping dependence of $`T_g`$ and $`T_f`$ suggests that both quantities may have the same scaling behaviour. The exponential doping dependence of $`T_g`$ has also been found in Y1-yCayBa2Cu3O6.02 for doping of up to 0.09 holes per planar Cu atom . This suggests that the trend of $`T_g`$ shown in Fig. 3 is common to all high-$`T_c`$ materials.

The present $`\mu `$SR data strongly suggest that there is a certain carrier concentration at which two magnetic energy scales associated with the superconducting state vanish. Although they are both much smaller than the normal state gap, there is a similarity in that they all decrease with doping and extrapolate to zero at $`x_c0.19`$  suggestive of a quantum critical point .

To further understand the physical properties near $`x_c`$ we have performed detail measurements of the absolute values of the in-plane penetration depth for samples prepared from the same batches as those measured by ZF-$`\mu `$SR. The lower panel of Fig. 3 shows the doping dependence of $`\lambda _{ab}^2(0)`$. In contrast to $`T_f`$ and $`T_g`$, $`\lambda _{ab}^2(0)`$ increases with doping and becomes nearly doping independent above $`x_c`$. We note the inverse relationship between $`\lambda _{ab}^2(0)`$ and the magnetic effects determined by $`\mu `$SR. The suppression of $`\lambda _{ab}^2`$ in the underdoped region has been reported from the early days of HTS . Recent penetration depth measurements in LSCO and HgBa2CuO4+x showed that not only the in-plane but also $`c`$-axis zero temperature superfluid density tend to saturate above $`x_c`$ . The present results allow us to identify a correlation between $`\lambda _{ab}^2(0)`$ and the magnetic effects determining $`T_g`$ and $`T_f`$.

Distinct changes at $`x_c`$ are not limited to the properties presented here, but have also been found in many other quantities of HTS. As mentioned already, the superconducting condensation energy and specific heat jump are maximum at $`x_c`$ and drop quickly with underdoping  furthermore the temperature dependence of $`\lambda _{ab}^2`$ only obeys the $`d`$-wave weak-coupling BCS formula for $`x>x_c`$ .

As a further check of the validity of $`x_c`$ in the context of the magnetic scales discussed here we performed ZF- $`\mu `$SR measurements for La2-xSrxCu1-yZnyO4 with $`x=0.15`$, 0.20, and $`y=0.01`$, $`0.02`$. Substitution with Zn slows down the spin fluctuations and greatly enhances the muon depolarisation rate at low temperatures . As a consequence of this, both $`T_g`$ and $`T_f`$, are significantly enhanced. For $`x(Sr)=0.15`$ and $`y=0.01`$ and 0.02, we find that $`T_g`$ $`=2.5K`$ and $`3.5K`$ and $`T_f`$ $`=15K`$ and $`25K`$, respectively. However, for $`x=0.20`$ the muon depolarization rate is unchanged by the Zn substitution. This result also suggests that there is indeed a critical point in the slightly overdoped region. Furthermore, the enhancement of the depolarisation rate by Zn doping confirms that it is the spin fluctuations on the CuO2 planes that are probed by $`\mu `$SR.

In conclusion, our results confirm that spin glass freezing is present in the superconducting state of underdoped LSCO , but show more clearly that it extends slightly above optimal doping. We have also provided evidence for very slow spin fluctuations in the same doping region. The doping level at which these characteristic features disappear is very close to the special point $`x=0.19`$  suggesting a connection between the magnetic energy scales identified in this work, the normal state energy gap and the anomalous properties of underdoped HTS. Our comparison with the measured superfluid density implies that $`low`$ $`frequency`$ spin fluctuations compete with superconductivity in the cuprates. The special doping level $`x=0.19`$ could therefore mark the point at which these effects disappear and the remaining high frequency magnetic excitations are favourable for superconductivity in the cuprates. The competition between quasi-static magnetic correlations with superconductivity is reminiscent of the way the Neel temperature of heavy fermion compounds goes to zero at a quantum critical point where the magnetic excitations are then most favourable for superconductivity.

We are grateful to A.D. Taylor of the ISIS Facility, Rutherford Appleton Laboratory for the allocation of muon beam time and W.Y. Liang for encouragement and support. C.P. thanks C. Bernhard, J.I. Budnick, S. Chakravarty, S.M. Hayden, S. A. Kivelson, W.Y. Liang, P.B. Littlewood, J.W. Loram, Ch. Niedermayer and C.M. Varma for useful discussions, and The Royal Society and Trinity College, Cambridge for financial support.


format: text

text:
# Signatures of Exo-solar Planets in Dust Debris Disks

## 1. INTRODUCTION

A significant limitation to unambiguous planet detection in circumstellar disks is the contaminating thermal emission from dust in the target systems. We propose to turn this hazard into an advantage for detecting and characterizing planetary systems.The major goal of the present Letter is to explore how the presence of the planet(s) impacts the disk via particular resonances responsible for the specific asymmetric features in the dust, and how these features would enable us to derive the major parameters of the yet invisible planetary system. The major obstacle in providing reliable, high-resolution numerical simulations – the particle-number limitation – is solved by our recently elaborated very efficient approach to numerical modeling of distributions of test particles in an external gravitational field (Ozernoy, Gorkavyi, & Taidakova 1998, 2000; Gorkavyi et al. 2000; Taidakova & Gorkavyi 1999).

## 2.  RESONANCE STRUCTURES IN DUSTY CIRCUMSTELLAR DISKS

Both the Poynting-Robertson (P-R) and stellar wind drags tend to induce dust inflow toward the star. As the dust passes by the planets in its infall, it interacts with them by accumulating in the outer planetary resonances, which, as we demonstrate below, could serve as an efficient means of planet detection. The existence of a resonance ring associated with the Earth’s orbit was predicted by Jackson & Zook (1989); such a ring was indeed discovered in the IRAS (Dermott et al. 1994) and COBE (Reach et al. 1995) data. The resonant design in the dust depends upon the capture probability, the particle lifetime in the resonance, the amplitude of libration, and the typical eccentricity of resonant particles. The total mass of dust contained in a particular resonance is the product of the probability for the dust particle being captured into the resonance, $`W_j`$; the lifetime spent in the resonance, $`\tau _j`$ (in a sufficiently dense environment, this time is due to particle collisions); and the total dust flow per unit time in the vicinity of the resonance, $`\mathrm{d}M/\mathrm{d}t`$, i.e. $`M_jW_j\tau _j(\mathrm{d}M/\mathrm{d}t)`$. Having introduced, for convenience, the mass of ‘background’ dust, $`M_b(1W_j)\tau (\mathrm{d}M/\mathrm{d}t)`$, one gets $`M_j/M_bW_j\left(1W_j\right)^1\left(\tau _j/\tau \right)`$, where $`\tau `$ is the minimal timescale of the P-R drag ($`\tau _{PR}`$), gravitational scattering ($`\tau _{gs}`$), and collisional ($`\tau _{coll}`$) characteristic times. Three particular cases are well distinguishable:

1: $`\tau \tau _{PR}<\tau _{gs},\tau _{coll}`$. Actually, this collisionless case is realized, for example, for Neptune in the Solar system. Assuming that $`W_j1W_j`$, one gets $`M_j/M_b\tau _j/\tau _{PR}`$. To get a discernable resonant structure, one needs to have $`\tau _j\tau _{PR}`$.

2: $`\tau \tau _{gs}<\tau _{coll},\tau _{PR}`$. This case is realized for a rather massive planet like Jupiter. If collisions are essential (as they are for Vega’s disk, see Backman 1998), then $`\tau _j\tau _{coll}`$ and one gets $`M_j/M_b\tau _{coll}/\tau _{gs}\genfrac{}{}{0pt}{}{_>}{^{}}1`$. The resonant structure has a high-contrast if $`\tau _{gs}\tau _{coll}`$.

3: $`\tau \tau _{coll}<\tau _{gs},\tau _{PR}`$. In this collisional case one has $`\tau _j\tau _{coll}`$ giving $`M_j/M_bW_j/(1W_j)`$, which results in a high-contrast resonant structure whenever $`W_j1`$.

It is worth emphasizing that if collisions of resonant particles are important (cases 2 and 3), the resonant pattern can nevertheless exist (cf. Wyatt et al. 1999) as long as a source of particle replenishment is available.

## 3. MODELLING CIRCUMSTELLAR DUST DISKS

3.1. The Procedure

The simplest planetary system under consideration, consisting of just one planet (orbiting a star of mass $`M`$ at a circular orbit of radius $`a`$) possesses only mean motion resonances, without most of secular resonances (except the Kozai resonances). We have computed about 300 model disks in order to explore the effects of various resonances adopting the following parameters: (1) the mass of the planet, $`m_{pl}`$; (2) the particle size coupled with stellar luminosity and mass, as described by the parameter $`\beta 0.3L_{}/(M_{}r)`$, which is the ratio of the stellar light pressure and the gravitational force applied to the dust grain of density 2 g $`\mathrm{cm}^3`$ and radius $`r`$ (in $`\mu `$m), with $`L_{}`$ and $`M_{}`$ the star’s luminosity and mass in solar units; and (3) lifetime of dust particles, $`\tau `$. Below, we describe some of these models.

3.2. Typical Resonant Structures Exemplified by the

Resonances 1:1, 2:1, 3:1, and 3:2

In the absence of a planet, the surface density of dust between the star and the pericenter of dust sources is constant, while the number density of dust varies as $`r^1`$, $`r`$ being heliocentric distance (Gorkavyi et al. 1997). The presence of a planet induces dramatic changes in the dust distribution. This is illustrated in Fig. 1, which shows the geometry of dust distribution in different mean motion resonances with the planet.

In the resonance 1:1, the orbits responsible for the dust design are of three kinds: (1) a horseshoe-like orbit that embraces the two Lagrange points L4 and L5 (see Fig. 1a); (2) an orbit with a small libration around L4 (Fig. 1b, upper feature); and (3) an orbit with a substantial libration around L5 (Fig. 1b, lower feature).

The dust structure in the resonance 2:1 looks similarly: the inner orbital loop precesses relative to L4 and L5 in the same way as the epicycle of particles in the 1:1 resonance. That loop can form a horse-like arc along the planet’s orbit (Fig. 1c) or can librate near either L4 or L5 (Fig. 1d illustrates the case of a large libration).

In the resonance 3:1, the dust particles librate around the L4 with a smaller amplitude (Fig. 1e); no horse-like libration has been observed. Finally, Fig. 1f illustrates the resonance 3:2 with an appreciable libration, which results in a symmetric four-clump pattern. When libration is small, the resonance contains a two-clump pattern.

In all cases mentioned above, the vicinity of the planet is free of dust. Each resonant dust structure is accompanied by other specific caverns in dust distribution (see Fig. 1).

3.3. Dependence of Structure on Location of Dust Sources

The scattered cometary population produces dust whose distribution in orbital elements tends to be rather smooth. A well-pronounced resonant structure in the dust can only be formed if there is a planet that massive enough to capture some dust particles into resonances and to eject the remaining non-resonant dust particles out (Case 2 in Sec. 2).

Dust sources like the Kuiper belt objects (KBOs; hereafter ‘kuiperoids’), of which a substantial part is in resonances, can produce the contrast resonant structures if (1) the drifting dust particles have a long lifetime and eventually are captured into resonances (the situation is similar to the formation of a resonant dust belt near the Neptunian orbit in the Solar system) or (2) both $`\beta `$ and the lifetime of dust particles are small enough, so that the particles produced by the resonant kuiperoids cannot leave the resonances and thereby they form resonant structures (Case 3 in Sec. 2). Orbits of resonant objects look similar to dust orbits shown in Fig. 1.

3.4. Dependence of Structure on Planetary Mass

. This dependence can be revealed from a variety of available runs with different parameters. A low-mass planet, like Earth or Neptune, forms a resonant ring along the planetary orbit associated with high-order resonances, like 3:2, 4:3, etc. (see, e.g. Fig. 1f). For more massive planets, the resonances like 2:1, 3:1, etc. become essential, whereas the innermost resonances overlap (e.g. Wisdom 1980). Since the libration amplitude decreases as $`m_{pl}`$ increases, clumps similar to those seen in Fig.1e (and not arcs like those shown in Fig. 1d) would be typical for $`m_{pl}`$ as large as $`m_J`$, where $`m_J`$ is Jupiter’s mass.

## 4. A LINK TO OBSERVATIONS: IMAGERY OF CIRCUMSTELLAR DISKS

Let us consider some observational signatures of planets embedded in disks. Figures 2a and 2b show the thermal emission from the simulated disks (seen face-on) with an embedded planet, and these are to be compared with the available observational data on Vega and $`ϵ`$ Eri shown in Figures 2c and 2d, respectively (the parent stars are seen almost pole-on).

As a hotter star, Vega might result in a more extensive, compared to $`ϵ`$ Eri, dust production from evaporating comets. We adopt $`\beta =0.3`$ (close to an upper limit of $`\beta =0.5`$ at which the dust particle escapes) and $`\tau _{coll}=2\times 10^4`$ planet revolutions consistent with $`\tau _{coll}10^7`$ yrs (Backman 1998). A planet as massive as $`2m_J`$ induces the resonances $`n:1`$, which results in two dusty clumps on either side of the planet, while the planet ejects the remaining non-resonant dust particles out (Case 2 in Sec. 2). Thus, even if the distribution of the cometary population around Vega might be symmetric, ejection of non-resonant particles by the planet could make the resonant dust pattern highly asymmetric (Fig. 2a), which is reminiscent of the observed two asymmetric clumps near Vega shown in Fig. 2c (Holland et al. 1998).

In contrast to Vega, $`ϵ`$ Eri is a less massive and a colder star with $`L\frac{1}{180}L_{\mathrm{Vega}}`$. Besides, the typical size of dust particles is several times larger than in Vega’s disk (Backman 1998), therefore we adopt $`\beta =0.002`$ for the dust in $`ϵ`$ Eri’s disk. Assuming that kuiperoids serve as the main source of dust in this disk, we find that the dust

distribution might reflect the distribution of the sources located mostly in the resonances such as 2:1 and 3:2. In this case, a comparatively small planetary mass of $`0.2m_J`$ is able to induce, via libration, two arcs, with a pair of clumps at their edges (Fig. 2b). Even if the dust is distributed in equal proportions to both resonances, one could easily get an asymmetry seen in Fig. 2b by assuming that 10% of comets in the resonance 2:1 have a smaller libration amplitude in the lower branch of the pattern. A similar asymmetric structure has been revealed in the submillimeter imagery of the $`ϵ`$ Eri disk shown in Fig. 2d (Greaves et al. 1998).

Our modeling, which indicates that Vega may have a planet of mass $`2m_J`$ at a distance of $`5060`$ AU, and $`ϵ`$ Eri may have a less massive planet of $`m0.2m_J`$ at a similar distance of $`5565`$ AU, is testable. Each resonant feature is stationary in the reference frame corotating with the planet, but it is not so for the observer at Earth. Therefore, if our interpretation of asymmetric clumps as a dynamical resonant structure is correct, the above asymmetric feature revolves around the star (probably, counterclockwise) with an angular velocity of $`1^{}.2`$$`1^{}.6\mathrm{yr}^1`$ (Vega) and $`0^{}.6`$$`0^{}.8\mathrm{yr}^1`$ ($`ϵ`$ Eri) – a prediction that can be tested within several years. If confirmed, the proposed interpretation of the structure in Vega- and $`ϵ`$ Eri–like circumstellar disks seen face-on would make it possible not just to reveal the embedded planet and determine its semimajor axis, but also to constrain its other basic parameters, such as the planet’s mass and even to pinpoint the position of the planet (see Figs. 2a and 2b).

## 5. CONCLUSIONS

As we demonstrated above, highly asymmetric structures in the circumstellar disk and a central ‘hole’ void of dust can serve as indicators of at least one planet embedded in the disk. The numerical simulations described above lead to the following major conclusions:

1. The outermost planet in an exo-planetary system can produce: (a) an asymmetric resonant dust belt consisting of arc(s) and one or more clumps intermittent with one or a few cavities; and (b) a central cavity void of dust.

2. The morphology of these belts and cavities (size, asymmetry, number of clumps, and their pattern) depends upon particular resonances involved and is eventually a function of the planetary mass, location of dust sources, and $`\beta `$.

3. The crucial test of the above picture would be the discovery of revolution of the resonant asymmetric structure around the star. For circumstellar disks in Vega and $`ϵ`$Eri the asymmetric design is expected to revolve, respectively, by $`1^{}.2`$$`1^{}.6`$ and $`0^{}.6`$$`0^{}.8`$ annually.

In sum, our numerical simulations offer a possibility to reveal the presence of the planet(s) in a circumstellar dusty disk seen pole-on and, under further development, could provide a method for determination of planetary parameters using the visible morphology of an outer part of the disk. The proposed approach seems to be sufficiently powerful to provide also the interpretation of disks seen edge-on in scattered light, such as the $`\beta `$ Pic disk revealed by STIS optical observations (Heap et al. 2000).

Acknowledgements. This work has been supported by NASA Grant NAG5-7065 to George Mason University. N.N.G. has ben supported through NAS/NRC Associateship Research program. The authors are thankful to Sally Heap help with programming and helpful discussions. We thank to W.S. Holland, J.S. Greaves and their teams for permission to reproduce their SCUBA images.

REFERENCES

Backman, D.E. 1998, in Exozodiacal Dust Workshop, eds. D.E. Backman et al. (NASA CP-1998-10155; Moffet Field: NASA/ARC), p.13

Dermott, S.F., Jayaraman, S., Xu, Y.L., Gustafson, B.A.S. & Liou, J.C. 1994, Nature 369, 719

Gorkavyi, N.N., Ozernoy, L.M., Mather, J.C. & Taidakova, T. 1997, ApJ 488, 268

Gorkavyi, N.N., Ozernoy, L.M., Mather, J.C. & Taidakova, T., 2000, in NGST Science and Technology Exposition, eds. E.P. Smith & K.S. Long (San Francisco: ASP). ASP Conf. Ser. 207, 462 (astro-ph/9910551)

Greaves, J.S. et al. 1998, ApJ 503, L133

Heap, S.R. et al. 2000, ApJ, in press

Holland, W.S. et al. 1998, Nature, 329, 788

Jackson, A.A. & Zook, H.A. 1989, Nature 337, 629

Ozernoy, L.M., Gorkavyi, N.N., & Taidakova,T. 1998, astro-ph/9812479

Ozernoy, L.M., Gorkavyi, N.N., & Taidakova,T. 2000, Planet. Space Sci. (in press)

Reach, W.T., Franz, B.A., Weiland J.L. et al. 1995, Nature, 374, 521

Taidakova, T.A. & Gorkavyi, N.N. 1999, in The Dynamics of Small Bodies in the Solar Systems: A Major Key to Solar Systems Studies, Eds. B.A. Steves, A.E. Roy, & V.G. Szebehely (NATO ASI Ser. C, 522; Boston: Kluwer), 393

Wisdom, J. 1980, AJ 85, 1122

Wyatt, M.C. et al. 1999, ApJ 527, 918


text:
# Newtonian hydrodynamics of the coalescence of black holes with neutron stars III: Irrotational binaries with a stiff equation of state

## 1 Introduction and motivation

The emission of gravitational waves in a binary system and the accompanying loss of angular momentum produce a decrease in separation, and for certain systems, will inevitably lead to coalescence (if the initial separation is small enough so that the decay will take place in less than the Hubble time). In fact, the binary neutron star systems known—PSR 1913+16 (Hulse & Taylor 1975), PSR 1534+12 (Wolszczan 1991)—have been observed to decay at a rate that matches the prediction of general relativity to high accuracy (Taylor et al. 1992; Stairs et al. 1998). There are no observed black hole–neutron star systems yet, but it is believed that they do in fact exist. Estimates of the event rates can be inferred from the statistics of Hulse–Taylor type systems and from theoretical studies of stellar evolution, and are expected to be about 10-6 to 10-5 per year per galaxy, implying several coalescences per year out to a distance of 1 Gpc could be taking place (Lattimer & Schramm 1976; Narayan, Piran & Shemi 1991; Tutukov & Yungelson 1993; Lipunov, Postnov & Prokhorov 1997; Portegies Zwart & Yungelson 1998; Bethe & Brown 1998; Kalogera 1998; Bełczyński & Bulik 1999).

The astrophysical interest in compact binaries such as the ones treated here is varied. They are primary candidates for detection by the gravitational detectors being currently constructed and expected to begin operation within the next few years, such as LIGO (Abramovici et al. 1992) and VIRGO (Bradaschia et al. 1990). When the binary separation is large (compared with the stellar radius), the system can be effectively considered to consist of point masses, and accurate waveforms can be calculated using the post newtonian approximation (see e.g. Kidder, Will & Wiseman 1992; Cutler et al. 1993; Blanchet et al. 1995). Detection of these signals will yield information about the masses and spins of the system. At small separations, a powerful but short burst of waves containing information about the radii and internal structure of the stars is expected. In particular, it may help constrain the equation of state of matter at high densities. At this stage, hydrodynamics will play an important role in the evolution of the system and must be taken into account to make any sense of possible observations. Clearly, this requires performing accurate, three–dimensional simulations in general relativity. These simulations have not yet been performed, but progress is being made in this direction by several groups (see e.g. Wilson, Mathews & Marronetti 1996; Lombardi, Rasio & Shapiro 1997; Baumgarte et al. 1997; Shibata 1999; Shibata & Uryū 2000).

Theoretical studies of the tidal disruption and coalescence of black holes with neutron stars were carried out many years ago (Wheeler 1971; Lattimer & Schramm 1974, 1976) but the computing resources to perform accurate numerical studies have only recently become available. Numerical work on the final stages of binary evolution and close interactions, including coalescence, has been carried out in the Newtonian case by many authors. These calculations may serve as a useful guide for future work, and also provide benchmarks one can use to compare the results of future general relativistic codes when used in the weak–field limit. The first studies of binary neutron star coalescence were done by Oohara and Nakamura (1989) who focused on the emission of gravitational waves. They later improved their work, using equilibrium conditions for tidally locked systems, and including the effects of radiation reaction in their calculations (Nakamura & Oohara 1989, 1991; Oohara & Nakamura 1990, 1992). Lai, Rasio & Shapiro (1993b, hereafter LRSb) developed a semi–analytical method that allowed them to generalize the study of incompressible rotating ellipsoids done by Chandrasekhar (1969) to compressible ones, using throughout a polytropic equation of state. Rasio & Shapiro (1992, 1994, 1995, hereafter RS92, RS94, RS95) extended this work to dynamical simulations of coalescing binaries using Smooth Particle Hydrodynamics (SPH), exploring the effect of the mass ratio and the stiffness of the equation of state on the global behavior of the system, as well as on the gravitational radiation signal. They clearly showed that hydrodynamics can play a crucial role in the evolution of such systems, making them dynamically unstable for small separations even in the Newtonian regime (Lai, Rasio & Shapiro 1993a, hereafter LRSa). Zhuge, Centrella & McMillan (1994, 1996) also used SPH and concentrated on the gravitational wave spectrum and on the effect of varying the neutron star spin and radius. The thermodynamic evolution of the fluid during and after coalescence has been studied using a physical equation of state by Davies et al. (1994) and Ruffert, Janka & Schäfer (1996) in double neutron star mergers, and recently by Janka et al. (1999) in black hole–neutron star mergers. This has important implications for gamma ray burst models (see below).

It is now believed that the gamma ray bursts (GRBs) are of cosmological origin (Meegan et al. 1992; see Fishman & Meegan 1995 for a review), with redshifts to several optical afterglows (Mészáros & Rees 1997a) having been measured recently (Metzger et al. 1997; Djorgovski et al. 1998; Kulkarni et al. 1998, 1999). The extreme energetics of these events (e.g. Kulkarni et al. 1999) and their variability (which arises at the source, Sari & Piran 1997) indicates that the “central engine” must involve a compact object of some sort. The preferred model involves the relativistic expansion of a fireball that produces the gamma rays as a result of internal shocks in the ejecta. Beaming of this fireball would reduce the energy requirements somewhat, and recent observations indicate that this could indeed be the case (Harrison et al. 1999; Stanek et al. 1999). What the central engine itself is made of remains to be seen. The bimodality in burst durations (Kouveliotou et al. 1995) with classes of long and short bursts (separated at about 2 seconds) hint that there might be two different processes producing the GRBs. Different mechanisms involving neutron stars and/or black holes have been proposed, such as binary coalescence of neutron stars (Paczyński 1986; Goodman 1986; Eichler et al. 1989; Narayan, Paczyński & Piran 1992), neutron stars with black holes (Paczyński 1991), catastrophic release of rotational energy through intense magnetic fields (Usov 1992; Kluźniak & Ruderman 1998; Spruit 1999) and failed supernova explosions (Woosley 1993). Janka & Ruffert (1996) estimated that the energy release in neutrinos during a binary neutron star coalescence was insufficient to power a GRB. However, after the initial coalescence an accretion torus may be formed around a resulting black hole, and such a structure might produce a GRB (Jarozyński 1993, 1996; Witt et al. 1994; Popham, Woosley & Fryer 1999; MacFadyen & Woosley 1999). The fireball might be powered by neutrino emission from this disc (e.g. Goodman, Dar & Nussinov 1987; Thompson 1994; Ruffert & Janka 1999) or by the energy of rotation of the black hole at the center of the disc (Mészáros & Rees 1997b) if one is present (Blandford & Znajek 1977).

It is not clear that the r–process nuclei can be produced by supernovae explosions alone in the observed abundances (see e.g. Meyer & Brown 1997; Freiburghaus et al. 1999a). Thus an alternative mechanism seems desirable. The ejection of nuclear matter during binary neutron star or black hole–neutron star mergers and its subsequent decompression could provide an adequate environment for this to occur, as has been suggested before (Lattimer & Schramm 1974, 1976; Symbalisty & Schramm 1982; Eichler et al. 1989). Recently, Newtonian hydrodynamical simulations by Rosswog et al. (1999) have addressed this problem, and the evolution of the ejected matter has been studied in greater detail by Freiburghaus, Rosswog & Thielemann (1999b). So far, they have found that indeed binary neutron star mergers seem to be promising candidates for this process to occur. Since the rates of black hole–neutron star mergers are probably comparable (see above), it seems natural to assume that this will also happen in such systems. Since our equation of state is of a simple form, we do not calculate nuclear reactions in any way, but merely wish to determine how much (if any) matter can be ejected during the dynamical coalescence and/or tidal disruption of the neutron star by the black hole, a necessary first step if it is to contribute to the abundances of r–process material.

In previous work (Lee & Kluźniak 1995; Kluźniak & Lee 1998; Lee & Kluźniak 1999a,b—hereafter Papers I and II) we have studied the dynamical coalescence of tidally locked black hole–neutron star systems. In all cases the neutron star was modeled as a polytrope and the dynamical simulation was followed using an ideal gas equation of state, as we have done for this work. We studied the case of a stiff ($`\mathrm{\Gamma }=3`$, Paper I) and a soft ($`\mathrm{\Gamma }=5/3`$, Paper II) equation of state, mainly to gauge what effect this would have on the outcome of the coalescence. We found that the compressibility of the neutron star affected the process greatly, with the binary surviving after a brief episode of mass transfer for the former, and unstable mass transfer leading to complete tidal disruption in the latter. We also confirmed through these dynamical simulations the presence of a dynamical instability for high mass ratios before Roche lobe overflow ocurred for $`\mathrm{\Gamma }=3`$. In the case of $`\mathrm{\Gamma }=5/3`$ mass transfer was initiated before the onset of dynamical instability (see Figure 1 in Paper II). In all cases we found that the outcome of the event was extremely favorable for the production of gamma ray bursts, with a region of low density along the rotation axis of the binary, and in some cases an accretion disc containing a few tenths of a solar mass around the black hole. This scenario seems particularly well suited for the class of short GRBs (Kluźniak & Lee 1998; see Fryer, Woosley & Hartmann 1999 for a variety of scenarios that could give rise to a GRB).

In the present paper, we extend our work to the case of irrotational binaries, which is believed to be a more realistic approximation to the conditions encountered in such systems (see below, section 3). We have chosen to model the neutron star with a stiff equation of state, using $`\mathrm{\Gamma }=3`$ as before, but also studying the case with $`\mathrm{\Gamma }=2.5`$. This may seem like a small difference, but the outcome is significantly altered, as shown below. We have adopted this value because it closely matches the adiabatic index that is believed to be appropriate for matter at nuclear densities (e.g. in the physical equation of state of Lattimer & Swesty 1991).

The paper is organized as follows. In section 2 we present the numerical method we have used, emphasizing the modifications we have implemented in our code. The initial conditions used for the dynamical simulations are presented in section 3, followed by a presentation of the simulations themselves in section 4. The influence of the choice of initial conditions is explored in section 5. A summary and discussion of our work is given in section 6.

## 2 Numerical method

For the calculations presented in this work, we have used the method known as Smooth Particle Hydrodynamics (SPH) (see Monaghan 1992 for a review and Lee 1998 for a description of our own code). The code is essentially the same one that was used for our previous simulations of tidally locked black hole–neutron star binaries (Paper I, Paper II) with some minor modifications, discussed below. Here we will not discuss the code in detail nor present any tests of the method, since it has now become widely known. We limit the discussion of numerical implementation to the new features of the code.

As before, the black hole is modeled as a Newtonian point mass of mass $`M_{BH}`$ with an absorbing boundary at the Schwarzschild radius $`r_{Sch}=2GM_{\mathrm{BH}}/c^2`$. Any SPH particle that crosses this boundary is absorbed by the black hole, whose mass and momentum are adjusted so as to ensure conservation of total mass and total linear momentum in the system.

The neutron star is modeled as a polytrope with a stiff equation of state, so that the pressure is given by $`P=K\rho ^\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }`$ and $`K`$ being constants (see Paper I). Unless otherwise noted, we measure mass and distance in units of the mass and radius of the unperturbed (spherical) neutron star (13.4 km and 1.4 M respectively), 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\mathrm{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\mathrm{km}}}\right)^3\left({\displaystyle \frac{M_{\mathrm{NS}}}{1.4M_{}}}\right)`$ (2)
$`\stackrel{~}{v}=0.39c\times \left({\displaystyle \frac{R}{13.4\mathrm{km}}}\right)^{1/2}\left({\displaystyle \frac{M_{\mathrm{NS}}}{1.4M_{}}}\right)^{1/2}`$ (3)
For the dynamical simulations presented here, we have used $`N38000`$ SPH particles to model the neutron star (except where noted). The initial (spherical) neutron star is constructed by placing the SPH particles on a uniform three–dimensional grid with particle masses proportional to the Lane–Emden density. This ensures that the spatial resolution is approximately uniform throughout the fluid. This isolated star is then allowed to relax for a period of thirty time units (as defined above) by including a damping term linear in the velocities in the equations of motion. The specific entropies of the particles are kept constant during this procedure (i.e. K is constant in the equation of state $`P=K\rho ^\mathrm{\Gamma }`$).

To perform a dynamical run, the black hole and every SPH particle are given the velocity as determined from the corresponding initial condition (see below) in an inertial frame, with the origin of coordinates at the centre of mass of the system. Each SPH particle is assigned a specific internal energy $`u_i=K\rho ^{(\mathrm{\Gamma }1)}/(\mathrm{\Gamma }1)`$ and the equation of state is changed to that of an ideal gas, $`P=(\mathrm{\Gamma }1)\rho u`$. The specific internal energy is then evolved according to the first law of thermodynamics, taking into account the contribution from artifical viscosity (see below). We vary the initial mass ratio $`q=M_{\mathrm{NS}}/M_{\mathrm{BH}}`$ in the binary by adjusting the mass of the black hole only.

### 2.1 Artificial viscosity

Artificial viscosity is used in SPH to handle the presence of shocks and avoid particle interpenetration. The momentum and energy equations are written as:

$`{\displaystyle \frac{d𝒗_i}{dt}}={\displaystyle \underset{j}{}}m_j\left(2{\displaystyle \frac{\sqrt{P_iP_j}}{\rho _i\rho _j}}+\mathrm{\Pi }_{ij}\right)\mathbf{}_iW_{ij}\mathbf{}\mathrm{\Phi }_i𝒂_i^{\mathrm{RR}},`$
and

$`{\displaystyle \frac{du_i}{dt}}={\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}m_j\left(2{\displaystyle \frac{\sqrt{P_iP_j}}{\rho _i\rho _j}}+\mathrm{\Pi }_{ij}\right)(𝒗_i𝒗_j)\mathbf{}_iW_{ij}.`$
Here $`\mathrm{\Pi }_{ij}`$ is the artificial viscosity term and $`𝒗`$, $`P`$, $`\rho `$, $`u`$, $`\mathrm{\Phi }`$, $`𝒂_i^{\mathrm{RR}}`$ and $`W`$ are the velocity, pressure, density, internal energy per unit mass, gravitational potential, radiation reaction acceleration (see below) and interpolation kernel respectively (for the kernel we use the form of Monaghan & Lattanzio (1985), see also Paper I). In previous work we have used the form introduced by Monaghan (1992), containing both shear and bulk viscosity. We have now changed the viscosity prescription and use the form presented by Balsara (1995), namely:

$`\mathrm{\Pi }_{ij}=\left({\displaystyle \frac{P_i}{\rho _i^2}}+{\displaystyle \frac{P_j}{\rho _j^2}}\right)(\alpha \mu _{ij}+\beta \mu _{ij}^2)`$
where

$`\mu _{ij}=\{\begin{array}{cc}\frac{(𝒗_i𝒗_j)(𝒓_i𝒓_j)}{h_{ij}(|𝒓_i𝒓_j|^2/h_{ij}^2)+\eta ^2}\frac{f_i+f_j}{2c_{ij}}\hfill & (𝒗_i𝒗_j)(𝒓_i𝒓_j)<0\hfill \\ 0\hfill & (𝒗_i𝒗_j)(𝒓_i𝒓_j)0\hfill \end{array}`$
and $`f_i`$ is the form function for particle i defined by

$`f_i={\displaystyle \frac{|\mathbf{}𝒗|_i}{|\mathbf{}𝒗|_i+|\mathbf{}\times 𝒗|_i+\eta ^{}c_i/h_i}}.`$
The factor $`\eta ^{}10^4`$ in the denominator prevents numerical divergences. The sound speed at the location of particle i is denoted by $`c_i`$, and $`\alpha `$ and $`\beta `$ are constants. The divergence and curl of the velocity field are evaluated through

$`(\mathbf{}𝒗)_i={\displaystyle \frac{1}{\rho _i}}{\displaystyle \underset{j}{}}m_j(𝒗_j𝒗_i)\mathbf{}_iW_{ij}`$
and

$`(\mathbf{}\times 𝒗)_i={\displaystyle \frac{1}{\rho _i}}{\displaystyle \underset{j}{}}m_j(𝒗_j𝒗_i)\times \mathbf{}_iW_{ij}`$
We use this new form of the viscosity because it vanishes in regions of large vorticity, when $`\mathbf{}\times 𝒗\mathbf{}𝒗`$. Since in our dynamical simulations accretion structures are often formed (where strong shear is present), we wish to limit the effect of artificial viscosity on the evolution of such structures, while at the same time retain the ability to deal with strong shocks (in regions of strong compression where $`\mathbf{}𝒗\mathbf{}\times 𝒗`$, the form function f approaches unity). Lombardi et al. (1999) have performed a study of the effects of different artificial viscosity prescriptions in several problems involving SPH, and they suggest using $`\alpha \beta \mathrm{\Gamma }/2`$, where $`\mathrm{\Gamma }`$ is the adiabatic index. For all the calculations presented here we have used these values for $`\alpha `$ and $`\beta `$.

### 2.2 Gravitational radiation reaction

We have modified our implementation of gravitational radiation reaction. Previously (see Paper II, also Davies et al. 1994) we used the quadrupole formula for angular momentum loss, and applied a back reaction force to each binary component. This approach is only reasonable as long as the neutron star is not tidally disrupted and can be reasonably approximated as a point mass. Usually this radiation reaction force is switched off when the binary separation (defined as the distance between the black hole and the centre of mass of the tidally deformed neutron star) decreases below a tidal disruption radius, given by $`r_{tidal}R(M_{\mathrm{BH}}/M_{\mathrm{NS}})^{1/3}`$. However, there are situations in which the initial episode of mass transfer in the binary does not disrupt the neutron star completely. Instead, a massive core survives while at the same time an accretion disc forms around the black hole (see the results for a tidally locked binary with a mass ratio of unity in Lee & Kluźniak 1995; Kluźniak & Lee 1998 and Paper I). In this case one still has a stable binary system, composed of the black hole and the surviving core, plus an accretion torus around the black hole. Gravitational radiation losses will affect the evolution of the new binary, but the system composed of the accretion torus and the core (the fluid) can hardly be considered to be a point mass. We thus realized that identifying the core was necessary in order to apply a gravitational raditation reaction term to it and the black hole. This needed to be done while the simulation was progressing, in order to update the calculation of angular momentum losses as needed.

In all cases where it was apparent that this identification was a necessity, the system had evolved through Roche lobe overflow into a configuration consisting of a stable binary (a low–mass neutron star around the black hole) and an accretion torus around the black hole (see Figure 4 in Paper I). It was apparent that the core consisted of matter that was gravitationally self–bound, and that the fluid in the accretion torus was in quasi Keplerian orbits around the black hole. Thus identifying the core became an issue of finding where the gravitationally self–bound matter was located. Our calculation of gravitational forces in the code uses a hierarchical tree structure (see Benz et al. 1990), and it simultaneously calculates what the gravitational potential energy is for every particle, $`\mathrm{\Phi }_i=\mathrm{\Phi }_i^{\mathrm{self}}+\mathrm{\Phi }_i^{\mathrm{BH}}`$. The first term arises from the contribution of all SPH particles, and the second from the presence of the black hole, given simply by $`\mathrm{\Phi }_i^{\mathrm{BH}}=GM_{\mathrm{BH}}m_i/|𝒓_{\mathrm{BH}}𝒓_i|`$. There is also kinetic energy in the fluid that needs to be considered, and we take this into account by calculating the kinetic energy of the SPH particles in the instantaneous frame of reference that is co–moving with the centre of mass of the core (to remove the influence of the kinetic energy of orbital motion). Since we do not know which SPH particles will make up the core beforehand, we use the position of the centre of mass of the core from the previous time step in the computation. Since in our initial condition (see below) the star has not overflowed its Roche lobe, the centre of mass of the core at $`t=0`$ coincides with the centre of mass of the entire neutron star. Thus the kinetic energy we calculate is $`K_i=m_i(\delta v)^2/2`$, where $`\delta 𝒗=𝒗_i𝒗_{\mathrm{core}}`$. If the total mechanical energy of an SPH particle, $`E_i=\mathrm{\Phi }_i^{\mathrm{self}}+K_i`$ is negative, then that particle is identified as belonging to the core of the star. We note that in the simulations presented here, we have not considered the internal energy of the fluid for this determination (we have performed several tests including this term, and the effect it has on the global evolution of the system is negligible). We keep track of the total core mass $`M_{\mathrm{core}}`$ and the position of its centre of mass. The radiation reaction formula is then applied to the black hole and to every SPH particle in the core as in Paper II, producing an acceleration

$`𝒂^i={\displaystyle \frac{1}{q(M_{\mathrm{BH}}+M_{\mathrm{core}})}}{\displaystyle \frac{dE}{dt}}{\displaystyle \frac{𝒗_{cm}^{\mathrm{core}}}{(v_{cm}^{\mathrm{core}})^2}},`$
while for the black hole the acceleration is given by

$`𝒂^{\mathrm{BH}}={\displaystyle \frac{q}{M_{\mathrm{BH}}+M_{\mathrm{core}}}}{\displaystyle \frac{dE}{dt}}{\displaystyle \frac{𝒗^{\mathrm{BH}}}{(v^{\mathrm{BH}})^2}}.`$
The energy loss is given by the quadrupole formula (Landau & Lifshitz 1975):

$`{\displaystyle \frac{dE}{dt}}={\displaystyle \frac{32}{5}}{\displaystyle \frac{G^4(M_{\mathrm{core}}+M_{\mathrm{BH}})(M_{\mathrm{core}}M_{\mathrm{BH}})^2}{(c|𝒓_{\mathrm{BH}}𝒓_{\mathrm{core}}|)^5}}.`$
The emission of gravitational radiation is calculated as described in Papers I and II (see also Finn 1989; Rasio & Shapiro 1992).

Of course, even with this new formulation of gravitational radiation reaction, the star can in principle be completely disrupted, so throughout our dynamical calculations, we monitor the separation of the black hole and the core, and compare it to the current tidal disruption radius. For $`\mathrm{\Gamma }=3`$ we switch gravitational radiation off if the distance between the centre of mass of the core and the black hole drops below $`r_{tidal}`$ (this never actually happened for the simulations presented here), and for $`\mathrm{\Gamma }=2.5`$ we do so when (and if) the mass of the core drops below one tenth of the initial neutron star mass (i.e. 0.14 M). We then perform an a posteriori check to ensure that the gravitational radiation reaction terms are not switched off too late or too early in the simulation. We have altered the criterion used to turn off radiation reaction from that used previously (Paper II) because we noticed that there sometimes occurred a slight jump in the slope of total angular momentum as a function of time ($`dJ/dt`$) at the moment when radiation reaction was turned off. The new procedure eliminates this discontinuity, and when radiation reaction is turned off it is no longer an important factor in the evolution of the system. To make sure that this is so, we continuously compute the radiation reaction timescale $`t_{RR}^1=256G^3M_{\mathrm{BH}}M_{\mathrm{core}}(M_{\mathrm{BH}}+M_{\mathrm{core}})/(5r^4c^5)`$ and an estimate of the current orbital period $`t_{orb}=2\pi /\sqrt{G(M_{\mathrm{BH}}+M_{\mathrm{core}})/r^3}`$, where r is the separation between the black hole and the centre of mass of the core. For the typical separations and masses in the black hole–core binary, by the time the core mass has dropped to 0.14 M, the radiation reaction timescale is much longer (by at least an order of magnitude) than the current orbital period.

## 3 Initial conditions

In previous work (Papers I and II) we have used tidally locked binary systems as initial conditions. This was done initially for two main reasons. First, it is relatively easy to set up self–consistent initial conditions for systems in a state of rigid rotation. Since we are interested in equilibrium configurations, the system can be relaxed in the co–rotating frame of the binary and Coriolis forces can be ignored. Second, we wanted to perform accurate comparisons to the results for two neutron stars using essentially the same approach (Newtonian hydrodynamics and a polytropic equation of state) as RS92, RS94, RS95. On the other hand, it has been known for some time that realistically, tidal locking is not expected in such binary systems, because the viscosity inside neutron stars is not large enough (Bildsten & Cutler 1992; Kochanek 1992).

Thus the case of irrotational systems is physically more interesting. Unfortunately, it is also more complicated. In this scenario the shape of the star is fixed in the co–rotating frame, but there are internal motions with zero circulation. Each component appears to be counter–spinning with the orbital angular velocity, and in an inertial frame, the star has effectively no spin. (We note here that, as pointed out by Rasio & Shapiro (1999), for double neutron star systems this creates a vortex sheet at the interface of the stars when they come into contact that is extremely difficult to handle numerically. For black hole–neutron star systems like the ones presented here, this problem does not exist, since there is no contact surface between the neutron star and the black hole, but a vacuum boundary condition.) Constructing equilibrium initial conditions for irrotational systems in black hole–neutron star binaries is a problem that has only recently been addressed (Uryū & Eriguchi 1999). A method for constructing approximate solutions was developped by LRSb. They used an energy variational method to find equilibrium solutions for a variety of binary systems assuming a polytropic equation of state and approximating the stars as compressible tri–axial ellipsoids. The black hole–neutron star binaries we treat in this work correspond to irrotational Roche–Riemann binaries (see Section 8 in LRSb). We have used this approach in order to construct initial conditions for all the dynamical runs presented here, except two (see below).

We build an initial condition by first constructing a spherical star of given radius and mass, as described in section 2. We then use the method of LRSb to calculate the orbital angular velocity of the binary and the semimajor axes of the Roche–Riemann ellipsoid for the appropriate choice of adiabatic index, initial mass ratio and binary separation (see Table 1). The semimajor axes of the fluid configuration can also be calculated from the SPH numerical solution using

$`a_i=\sqrt{{\displaystyle \frac{5I_{ii}}{\kappa _nM_{\mathrm{NS}}}}}`$
where

$`I_{ii}={\displaystyle \underset{j}{}}m_j(x_j^i)^2.`$
The stiffness of the equation of state enters these equations through the parameter $`\kappa _n`$ ($`\kappa _n=0.815`$ for $`\mathrm{\Gamma }=3`$ and $`\kappa _n=0.75`$ for $`\mathrm{\Gamma }=2.5`$). The first and second semimajor axes of the tri–axial ellipsoid lie in the orbital plane, with the first one along the line joining the two binary components. The third axis is oriented perpendicular to this plane (along the axis of rotation). The position of each SPH particle is then re–scaled (independently along each coordinate axis) so that the new fluid configuration has the appropriate semimajor axes. It is this new ellipsoid that is then used as an initial condition for the corresponding dynamical run. The initial velocity is given by the orbital angular velocity (for the azimuthal component) plus the (small) radial velocity corresponding to point–mass inspiral. We show in Figure 1 the variation in total angular momentum as a function of binary separation for irrotational Roche–Riemann binaries (with various mass ratios and adiabatic indices in the equation of state), as calculated using the method of LRSb. The curves exhibit a turning point as the separation is decreased, indicating the presence of an instability in the system. Strictly speaking, the ellipsoidal approximation breaks down close to this turning point, and a full equilibrium solution is necessary (see Uryū & Eriguchi 1999). However, we have chosen the values of the initial separation for our dynamical runs $`r_i`$ to be slightly above the turning point (see Table 1). The ellipsoidal approximation is then still reasonable, and our equilibrium configurations have not yet reached the point where the neutron star will overflow its Roche lobe. When the dynamical simulation is initiated, the separation will decrease due to the emission of gravitational waves, and mass transfer will start promptly. If one were to construct the appropriate Roche–Riemann ellipsoid at a large separation and then allow the system to evolve dynamically through gravitational wave emission, the neutron star would be spun–up to some degree for several reasons. First, since the neutron star cannot react instantly to changes in the tidal potential produced by the black hole, a tidal lag angle will be induced as the separation decreases, so that the tri–axial ellipsoid and the binary axis are no longer aligned (see Lai, Rasio & Shapiro 1994). This tidal lag increases as the separation is decreased. Immediately before Roche lobe overflow occurs, the presence of a dynamical instability can make it quite large (on the order of $`10^{}`$). Second, any amount of viscosity will transfer orbital angular momentum to the spin of the neutron star, contributing to orbital decay. This applies both to a real, physical viscosity, but also to any numerical effect that might be present in our scheme. We have thus tried to minimize the effect of numerical noise and viscosity on the initial evolution of the binary by proceeding as described above, always keeping in mind that our initial condition is not fully self–consistent. We hope to improve on this in future work. The initial separations we have chosen are similar to what we have presented before for the case of tidally locked black hole–neutron star binaries (Paper I, Paper II).

We also include in Table 1 the initial parameters for two runs (A31S and B31S) that used initially spherical neutron stars for the dynamical calculations, with a Keplerian orbital angular velocity. We have performed these runs to gauge the effect non–equilibrium initial conditions will have on the evolution of the system. The results of these runs are presented in section 5.

## 4 Results

We now present the result of dynamical calculations performed using the initial conditions described in section 3.

### 4.1 Morphology of the mergers

The binary separation decreases as a result of angular momentum losses to gravitational radiation, and in every case the neutron star overflows its Roche lobe within one orbital period, initiating mass transfer to the black hole.

Figures 2 and 3 show a sequence of density contour plots in the orbital plane of the binary for runs A50 and B31. The neutron star becomes elongated (a), and an accretion stream forms between it and the black hole. This stream winds around the black hole, colliding with itself and producing a thick accretion torus (b), (c), (d). At the same time as the accretion torus is being formed, a long tidal tail of neutron star matter being ejected through the outer Lagrange point appears. The most striking difference between the two cases is that for $`\mathrm{\Gamma }=3`$ (run A50) the neutron star core clearly survives this first mass transfer episode as a coherent body, while for $`\mathrm{\Gamma }=2.5`$ (run B31) the star is almost completely disrupted, although one can discern a bulge in the tidal tail in panels (c) and (d) in Figure 3. The core (in the case of run A50) and the bulge (in the case of run B31) make a second periastron passage around the black hole–see panels (f,g,h). In the first case, a second, less pronounced stream forms, feeding the accretion disc, as well as a smaller, secondary tidal tail (Figure 2g-h). In the latter case, the final disruption of the core is not so evident, but one can see that the accretion disc has a complex structure, with two partial rings on one side of the disc (Figure 3g-h).

The separation between the black hole and the centre of mass of the core can be seen in Figure 4, where we plot it for all runs. The separation initially decreases at a rate consistent with the quadrupole approximation for a point–mass binary, and subsequently deviates from this as the separation decreases at an even faster rate, particularly for the higher mass ratios ($`q=0.5`$ and $`q=0.31`$). This is a reflection of the importance of the hydrodynamics in the evolution of the system. It was pointed out by LRSa that for such compact binaries (black hole–neutron star or double neutron star systems), a purely Newtonian dynamical instability can make the orbital separation decay on an orbital timescale in the absence of gravitational radiation losses, especially for relatively stiff equations of state such as the ones being considered here. This is precisely what is seen directly in the curves in Figure 4, and also observed for the case of tidally locked binaries with a stiff equation of state (Paper I). The initial episode of mass transfer then makes the separation increase, as we have seen above. In all cases with $`\mathrm{\Gamma }=3`$ the binary is not disrupted, and the surviving core makes successive periastron passages, transferring some mass to the accretion disc or directly to the black hole each time. We note here that for all of these runs, the radiation reaction force is always included in the equations of motion (see Table 1). For $`\mathrm{\Gamma }=2.5`$, the neutron star is completely disrupted after (at most) the second periastron passage, and the accretion disc evolves steadily towards an axisymmetric configuration. Our formulation of radiation reaction is strictly valid only for circular orbits. However, the orbital eccentricity (when the core survives) is small ($`e0.2`$) and thus we believe the errors in the rate of energy and angular momentum loss are not excessive (20% and 10% respectively). In all cases, the primary tidal tail (formed during the first episode of mass transfer) persists as a well–defined large scale structure throughout the simulation. At late times, we see the formation of knots through a sausage instability at approximately regular intervals all along the tail (this has been observed in this kind of simulation before, see e.g. RS94). Their masses are fairly uniform, of order $`2\times 10^3`$. This does not occur for soft equations of state (with $`\mathrm{\Gamma }2`$, see Paper II and RS92, RS95).

For a fixed value of the adiabatic index, varying the initial mass ratio $`q`$ leads to qualitatively similar results, although there are quantitative differences. The most important of these is that for $`\mathrm{\Gamma }=3`$ and lower mass ratios (runs A31 and A20) the amount of matter that is left in orbit around the black hole after each periastron passage is much smaller than for run A50, since most of it is directly accreted by the black hole. This is in part due to the fact that the black hole is physically larger in those cases, and thus the infalling gas is more likely to cross the horizon than to wind around the black hole and produce a larger accretion disc.

We show in Figure 5 the mass accretion rates onto the black hole as a function of time. The peak accretion rates (see Table 2) occur during the initial episode of mass transfer, and are of order 0.05 (equivalent to 0.6 M ms-1), largely independent of the initial mass ratio and the adiabatic index. There are however, important qualitative and quantitative differences at later times. Whenever a torus forms around the black hole, the accretion rate decreases rather smoothly as the simulation progresses (runs A50, B50, B31 and B20). For the cases in which, as stated above, there is practically no such structure (A31 and A20), the accretion rate is episodic and directly reflects the periastron passages of the surviving core. This is not easily seen at our current level of resolution, since in these runs, the amount of mass transferred in these events is small ($`\mathrm{\Delta }M10^3`$) and hence very few SPH particles are actually lost to the black hole (on the order of 50). In Figure 5a the mass transfer episode occurring at $`t140`$ is barely visible.

During our simulations, the total angular momentum in the system, J, is not conserved for two reasons. First, since we have included a gravitational radiation reaction term in the equations of motion, some of it will be lost from the system.

We have plotted J as a function of time in Figure 6 for all runs, and this effect can be seen in the early stages of the simulation, before substantial mass transfer has taken place and hydrodynamic effects have become important (before $`t10`$). After this, the decrease is due primarily to accretion onto the black hole (when matter is accreted through the horizon, the total mass and linear momentum in the system are conserved, but not the total angular momentum),

and we assume that it goes into spinning it up, thus giving it a finite Kerr parameter $`a=J_{\mathrm{BH}}c/GM_{\mathrm{BH}}^2`$, shown in Table 2 (column 9) at the end of the simulation and assuming that at $`t=0`$, $`a=0`$.

Why the neutron star core moves out after the first mass transfer episode can be seen directly in Figure 7, where we plot the specific angular momentum of the core as a function of time. This initially decreases slightly (for the same reasons we have outlined above), and then increases sharply. The oscillations in the curves for $`\mathrm{\Gamma }=3`$ at late times reflect each successive periastron passage and corresponding episode of mass transfer, and show that there is no significant orbital decay (as is also evident from Figure 4a). The slight variations arise because our numerical definition of the core (for gravitational radiation reaction purposes) makes its mass vary somewhat as the binary separation changes, thus altering its total angular momentum content as well. The mass transfer is not stable, and the evolution of the system at this stage is not driven by the loss of angular momentum to gravitational waves as one might initially expect but by hydrodynamic processes. This is the same qualitative result that was obtained for tidally locked binaries with a stiff equation of state (Kluźniak & Lee 1998; Paper I), and the reason is essentially the same. For polytropes, the mass-radius relationship is $`RM^{(\mathrm{\Gamma }2)/(3\mathrm{\Gamma }4)}`$. For $`\mathrm{\Gamma }=3`$ this becomes $`RM^{1/5}`$ and for $`\mathrm{\Gamma }=2.5`$, $`RM^{1/7}`$. Thus the neutron stars presented here respond to mass loss by shrinking slightly, rather than expanding as should be the case realistically (Arnett & Bowers 1977). In the event of conservative mass transfer in a binary system, if the donor is the lower mass component, the orbital separation will increase. Strictly speaking we are not dealing with a conservative system in this case, but the general response of the system is the same. Since the neutron star shrinks upon losing mass, it eventually cuts off the mass transfer stream to the black hole when it no longer fills its Roche lobe, and we are left with a stable binary. The surviving orbit is slightly eccentric, and thus successive periastron passages allow a small amount of mass transfer to take place. The evolution of the binary is driven by two competing effects. On the one hand, gravitational radiation reaction drains angular momentum from the system and thus tends to make the separation decrease (it will also circularize the orbit, but on a longer timescale than we have modeled here). On the other, every time mass transfer occurs, it tends to increase the separation as outlined above. It would appear from our calculations that these effects cancel each other out on a timescale of many orbits, but not in a manner that allows steady mass transfer, as it has been suggested before (Blinnikov et al. 1984; Portegies Zwart 1998). This is simply because the initial episode of mass transfer (due to angular momentum loss to gravitational radiation, or to a hydrodynamic instability) is much too violent to allow the star to react fast enough. For softer equations of state with $`\mathrm{\Gamma }=2`$ and $`\mathrm{\Gamma }=5/3`$, whether the initial condition is that of a tidally locked (Paper II) or an irrotational binary (Lee 2000, in preparation), the mass transfer process itself can be unstable and lead to complete tidal disruption on a dynamical timescale (in this case the mass–radius relationship is such that the star expands upon losing mass). That stable mass transfer in a black hole–neutron star binary is not possible had been pointed out by Bildsten & Cutler (1992) and Kochanek (1992) using essentially the same arguments we have outlined above to explain our numerical results.

The various energies in the system are shown in Figure 8

for runs A50 and B31. The violent episode of mass transfer in the initial stages of the coalescence is reflected in the drop in total internal energy, while the survival of the binary in a slightly eccentric orbit can be seen from the small oscillations in kinetic and gravitational potential energy for run A50. In contrast, the corresponding plot for run B31 shows only a monotonic change in the curves at late times, when the neutron star has been completely disrupted.

### 4.2 Accretion disc structure

In Table 2 we show several parameters pertaining to the final configuration of the system (at $`t=t_f`$). We calculate the disc mass $`M_{disc}`$ by computing the amount of matter which has angular momentum larger than a critical value, necessary for remaining in orbit about the black hole. This matter has $`j>j_{crit}=\sqrt{6}GM_{total}/c^2`$ (as in Paper II). This value is not zero even for the cases in which no appreciable disc has formed (runs A31, A31S and A20 in particular) because we are not excluding matter in the surviving core in this calculation. We are assuming that the neutron star core will eventually be disrupted and part of it will be accreted by the black hole. Of course for these cases this is a very rough estimate since an explosion of the core (see below) could redistribute angular momentum significantly. In the cases where the neutron star is disrupted and a disc has formed, its structure is becoming roughly azimuthally symmetric by the end of the simulation (excluding the presence of the long tidal tail). Based on the disc masses and final accretion rates (admittedly noisy and poorly resolved) we infer a rough estimate for the lifetime of the disc as $`\tau _{disc}\dot{M}_{final}/M_{disc}`$. The values shown in the table are then between 70 ms and 100 ms. The viscosity in our numerical scheme is purely artificial, so these values must be taken with caution, and probably underestimate the true lifetime of the disc.

We show in Figure 9 density contours for runs A50 and B31 at the end of the simulations.

In run A50, a secondary episode of mass transfer can be seen directly, while a few of the knots formed in the initial tidal tail for run B31 can also be seen. The azimuthally averaged density and internal energy profiles at $`t=t_f`$ for run B31 are shown in Figure 10.

The internal energy has a maximum in the inner regions of the disc, and becomes somewhat flatter in the outer regions, at $`u/100010^5`$, equivalent to $`1.37\times 10^{18}`$ erg g-1, or 1.45 MeV/nucleon. The fluctuations seen in the curves reflect the fact that the accretion disc is not yet azimuthally symmetric. The rotation curve is sub–Keplerian, indicating that pressure support is an important factor in the structure of the disc.

It is apparent in Figure 9 that the region directly above and below the black hole is largely free of matter. This is of crucial importance in the context of gamma–ray bursts, since a low baryon loading (of the order of 10-5 M) is necessary to permit the formation and expansion of an ultrarelativistic fireball (Mészáros & Rees 1992, 1993). Figure 11 shows the mass enclosed in a cone of opening angle $`\mathrm{\Delta }\theta `$ directly above and below the black hole, along the rotation axis. All runs show that this axis is clear of matter within about 10 down to approximately $`10^5`$ M, and only modest collimation of an outflow along the rotation axis is required to avoid contamination (see the last three columns in Table 2). This is the first time we have been able to reach such a low mass resolution limit regarding the baryon pollution, and is due to the improved resolution of our dynamical simulations (the runs shown in Papers I,II used at most $`17,000`$ SPH particles). In the event of disc formation, a fraction of the energy of rotation of the black hole could be extracted from the system through the Blandford–Znajek (1977) mechanism. The black hole is spun up during the coalescence because of the angular momentum of the accreted matter, to rates varying between 20% and 40% of the maximal rotation rate (see Table 2).

### 4.3 Ejected mass and r–process

We calculate the amount of mass that may be dynamically ejected from the system as in Paper II, by identifying matter with a positive total energy (kinetic+gravitational potential+internal) at the end of each simulation (see Table 2). This matter can be found in the outer parts of the tidal tails that are formed during the initial disruption of the star (for $`\mathrm{\Gamma }=2.5`$) and also in the secondary tails formed at secondary periastron passages of the core (for $`\mathrm{\Gamma }=3`$, see Figure 12).

Since it is not clear that the fluid in this later category will truly become unbound from the system (because of possible interactions with other tidal tails or with the surviving core itself), we have not included it in what we consider to be ejected mass, shown in Table 2 (this eliminates only about 10% of the ejected mass for $`\mathrm{\Gamma }=3`$). Within our range of initial mass ratios, there appears to be little difference in the amount of ejected mass. This contrasts with the results of Janka et al. (1999), who find that the dynamical ejection of matter is highly suppressed at higher mass ratios. The reason for this discrepancy could be due to their use of the equation of state of Lattimer and Swesty, and to their different implementation of gravitational radiation reaction (they use the formalism of Blanchet, Damour & Schäfer (1990) as opposed to our use of the quadrupole approximation). We do, however, observe a slight increase for the softer equation of state, probably because in this case the star is completely disrupted, whereas for $`\mathrm{\Gamma }=3`$ a substantial amount of mass remains locked in the core. We find that between $`3.1\times 10^2`$ and $`4.6\times 10^2`$ M can be dynamically ejected from the system. This is about the same as what was found for double neutron star mergers by Rosswog et al. (1999a), who performed a Newtonian study and used the physical equation of state of Lattimer & Swesty (1991). A more detailed hydrodynamic study has been carried out by Freiburghaus et al. (1999b) to more accurately determine the production of heavy elements through the r–process in this ejected material. While our equation of state does not allow us to follow such evolution, determining how much material can be expelled to the interstellar medium from a dynamical coalescence is a necessary first step. In this aspect we consider our results to be upper bounds, mainly because the effects of general relativity will probably make the ejection of neutron star material from the potential well of the binary more difficult.

### 4.4 Emission of gravitational waves

As in previous work (see also Finn 1989; RS92) we have calculated the emitted gravitational radiation waveforms during the merger in the quadrupole approximation for every run. We show in Figure 13 the waveforms and luminosities for runs A20 and B20 (the curves are qualitatively similar for runs A50, A31 and B50, B31 respectively), compared to what these results would be for a point–mass binary decaying through angular momentum loss in the quadrupole approximation. Since the neutron star is completely disrupted for the runs with $`\mathrm{\Gamma }=2.5`$, the amplitude of the waveforms quickly drops effectively to zero (we show upper bounds in Table 3), as the accretion torus becomes ever more azimuthally symmetric. The luminosity initially rises as the separation decreases (before any significant amount of mass transfer has taken place, until $`t10`$) and then decays quickly, exhibiting a single, sharp peak. For $`\mathrm{\Gamma }=3`$ the behavior is initially similar, except that now there remains a coherent signal that clearly has the signature of a binary system, since the core remains in orbit around the black hole for the duration of the simulation. The orbit is slightly elliptical, and the successive peaks in the gravitational radiation luminosity (at $`t75,130,185`$) correspond to the later periastron passages of the core. We detail the information extracted from the waveforms in Table 3, where we show the maximum and final amplitudes for the waveforms, the peak luminosity and the total energy radiated away by the system, and the efficiency of gravitational wave emission $`ϵ=\mathrm{\Delta }E/M_{\mathrm{total}}c^2`$. For reference, $`L_{max}=1`$ (in the units given in the table) corresponds to $`3.036\times 10^{55}`$ erg s-1and $`\mathrm{\Delta }E=10`$ is equivalent to $`3.48\times 10^{52}`$ erg. The variations for different values of the mass ratio (at fixed $`\mathrm{\Gamma }`$) arise because we have normalized results to the initial neutron star mass (1.4 M). However, it is clear that there are quantitative differences depending on the adiabatic index (at a fixed initial mass ratio). The stiffer the equation of state, the more important tidal effects are, and these contribute to orbital decay, as we have seen. Thus the binaries with $`\mathrm{\Gamma }=3`$ are disrupted earlier during the encounter, reaching consistently lower values in amplitude for $`h_{max}`$ and $`L_{max}`$. The total energy emitted by the system should not be taken as an absolute value, since it depends on the choice of the origin of time. We can nevertheless perform comparisons between runs with this figure, since our initial separations are always comparable for a given choice of $`q`$. Here too we see that less energy is emitted for the stiffer equation of state, despite the fact that the binary survives and continues to produce gravitational waves. This is simply because the separation and mass ratio have changed so much during the mass transfer episodes that emission is quickly suppressed by more than two orders of magnitude (see Figure 13b).

## 5 Influence of initial conditions on the dynamical evolution of the system

As mentioned before, we have performed two runs, A31S and B31S, in which the initial condition uses a spherical neutron star, relaxed in isolation. The initial orbital angular velocity used for these dynamical runs is that of a point mass binary with the same initial mass ratio and separation. Otherwise all initial parameters are as in runs A31 and B31 respectively (see Table 1). Constructing fully self–consistent initial conditions for close binaries is not an easy task, particularly for non–synchronized binaries like the ones treated in this paper. The complications that arise also depend on the particular numerical method being used, and there is a trade–off between numerical accuracy and physical sense of the results that is always present. All other things being equal, at larger binary separations, tidal effects are smaller (how much smaller depends on the details of the equation of state), and thus taking a spherical star as an initial condition can be quite reasonable. However, one then needs to evolve the system for a large number of dynamical times (e.g. draining angular momentum through gravitational wave emisison) before the actual coalescence takes place. This has two very important drawbacks. First, much CPU time is wasted modeling the system at high resolution while the hydrodynamics has little effect on the evolution. Zhuge et al. (1994, 1996) have performed dynamical simulations of binary neutron star coalescence using SPH, starting with initial conditions at large separations and removing angular momentum from the system through gravitational wave emission as we have done here. Their spatial resolution was severely limited (the neutron star was modeled with approximately 4000 particles at most) because they had to evolve the system for a large number of dynamical times. Second, spurious numerical effects can degrade the initial condition one imposed at the start (through viscous dissipation and angular momentum transport) so that we have little control over the physical state of the system when the final merging occurs. Thus, many of the simulations of binary coalescence reported in the literature for double neutron star systems (e.g. Oohara & Nakamura 1989; Benz et al. 1990; Davies et al. 1994, Ruffert et al. 1996; Rosswog et al. 1999), and for black hole neutron star binaries (Janka et al. 1999), have used spherical neutron stars as an initial condition, regardless of the initial separation, which is nevertheless small enough so that the two problems stated above are avoided. We note here that it is much easier to build even approximate equilibrium solutions to initial conditions if the equation of state is of a simple form like the one we have used here, rather than a physical equation of state like the one used by Janka et al. (1999) and Rosswog et al. (1999) (this is another trade-off one must balance). Here we wish to explore how one aspect (tri–axial ellipsoid vs. spherical star) of the choice of initial conditions will affect the outcome of the coalescence, in a qualitative and quantitative fashion. As noted above (see section 3), we remind the reader that an irrotational Roche–Riemann ellipsoid does not correspond to a fully self–consistent initial condition either.

The main effect of using a spherical star as an initial condition appears as soon as the dynamical simulation starts, since the neutron star must respond the the instantaneous appearance of a tidal field produced by the black hole. Thus a tidal bulge quickly appears, initially along the axis joining the two binary components. This act drains energy from the orbital motion, and makes the separation decay slightly faster than it otherwise would (see Figure 14a).

The appearance of the tidal bulge also induces radial oscillations of the star which damp out after a few dynamical times, and can be appreciated in Figure 14b, where we plot the total internal energy of the neutron star for runs A31, A31S, B31 and B31S. All curves show variation, since runs A31 and B31 do not correspond to full equilibrium solutions either, but their oscillations are somewhat smaller. This initial perturbation is responsible for all subsequent differences in the variables shown in Tables 2 and 3, and are easily understood. Since the separation decreases at a higher rate, the initial encounter with the black hole is more violent, making the peak accretion rates larger by about 8%, and the surviving core (for $`\mathrm{\Gamma }=3`$) less massive by about 20%. The gravitational wave signal (see Figure 15)

is also affected for the same reasons. The peak amplitude in the waveforms and the peak luminosity are slightly higher (because the minimum separation reached during the first episode of mass transfer is smaller) but decay faster, and the final amplitude (when the binary survives) is smaller (largely because the final mass ratio is lower). However, the total energy radiated away through gravitational waves is significantly lower (by about 15%), because the initial encounter is much more brief. Finally, the faster initial decay is also reflected in a slight phase shift in the gravitational wave signal (see Figure 15b).

## 6 Discussion and conclusions

We have presented the results of three–dimensional hydrodynamical simulations of the coalescence of a neutron star with a black hole. We use a stiff equation of state to model the neutron star as a cold polytrope with $`\mathrm{\Gamma }=3`$ and $`\mathrm{\Gamma }=2.5`$, which is believed to be appropriate for matter at nuclear densities. The black hole is modeled as a Newtonian point mass with an absorbing boundary at the Schwarzschild radius. We have explored a range of mass ratios $`q=M_{\mathrm{NS}}/M_{\mathrm{BH}}`$, and follow the dynamical evolution of the system for approximately 23 ms. It is important to keep in mind that our simulations are completely Newtonian, except for our treatment of gravitational radiation reaction on the hydrodynamics (see section 2.2). Since tidal locking is not expected in this type of system (Bildsten & Cutler 1992; Kochanek 1992) we have used initial conditions corresponding to irrotational binaries in equilibrium, approximating the neutron star as a compressible tri–axial ellipsoid (see LRSb). The dynamical simulations are begun when the system is on the verge of initiating mass transfer.

The loss of angular momentum to gravitational wave emission, together with hydrodynamic instabilities, always drives these systems to initiate mass transfer within one orbital period after the start of the dynamical simulations. For $`\mathrm{\Gamma }=3`$ this episode always leads to the survival of the neutron star core and its transfer to a higher orbit. A long tidal tail of matter stripped from the star appears from the outer Lagrange point, and for high mass ratios ($`q=0.5`$ in our simulations) an accretion disc is clearly visible (for lower mass ratios this is not so evident at our current level of resolution). The new binary system is clearly stable, and the slightly elliptical orbit ($`0.1e0.2`$) is small enough to allow secondary mass transfer events to occur at each periastron passage (two or three more occur during our simulations). Since the mass ratio in the new binary differs substantially from the initial one, and the separation is larger, the timescale for decay of the orbit due to gravitational wave emission is lenghtened considerably, and becomes up to two orders of magnitude greater than the orbital period (calculated for point–mass binaries in the quadrupole approximation). Thus the lifetime of the system may be on the order of tenths of one second. For $`\mathrm{\Gamma }=2.5`$, the neutron star is eventually completely disrupted and an accretion disc is formed around the black hole, but over a longer timescale (essentially up to and including the second periastron passage of the core) than in the previous case. A tidal tail of material thrown out from the outer Lagrange point is also formed as described above. These tails survive as well–defined structures throughout the simulations.

The overall morphology of the events is directly reflected in the gravitational radiation waveforms. In the cases where the stellar core survives, the signal continues to exhibit a finite amplitude (albeit greatly reduced) at a frequency corresponding to the new orbital period (approximately 350 Hz vs. 800 Hz at the start of the dynamical simulations, since the frequency of the gravitational radiation signal is twice that of orbital motion for circular orbits). The eccentricity of the orbit produces periodic peaks in the gravitational wave luminosity. This remnant signal is completely absent in the case of $`\mathrm{\Gamma }=2.5`$, where the waveforms (and luminosity) drop apruptly and practically to zero after the star is disrupted and the accretion disc is formed. The dramatic dependence on the stiffness of the equation of state shows the type of information that could be gleaned from an observation of gravitational waves in the near future, and is similar to what has been found for double neutron star systems by Rasio & Shapiro (RS94), where the persistent emission was due to a lack of axisymmetry in the central object left after the coalescence.

In all cases, the fluid contained in the outer parts of the tidal tails (amounting to 10-2–10-1 M) appears to be on outbound trajectories (it has enough mechanical energy to escape the gravitational potential well of the black hole–debris torus system). This probably represents an upper bound, since the effects of general relativity are likely to lower this value. As we have pointed out above, this may have important implications for the abundances of r–process material in the galaxy. We refer the reader to the more detailed analysis performed by Rosswog et al. (1999) and Freiburghaus et al. (1999b) where they have used the equation of state of Lattimer & Swesty (1991) and a detailed thermodynamic and nuclear network calculation.

The use of adequate initial conditions for detailed dynamical simulations is always an important concern. For tidally locked binaries, setting up equilibrium initial conditions can be done in a relatively simple way (Paper I, Paper II). This is not the case if this restriction is lifted, as it has been for this work. We have chosen to use the semi–analytical method of LRSb to construct approximate equilibrium solutions to the problem of a compressible tri–axial ellipsoid in orbit about a point mass companion (see Uryū & Eriguchi 1999 for full equilibrium configurations). We remark that this does not mean that our initial conditions are entirely self–consistent, since at small separations, a dynamical tidal lag angle is induced in the binary (see Lai et al. 1994), and we have necessarily neglected this effect for now. We have carried out two dynamical simulations in this work that use a spherical neutron star relaxed in isolation as an initial condition, and find that there are quantitative differences in the evolution of the system (see section 5). These can be directly traced to the fact that the spherical star must respond instantaneously to the presence of the gravitational field of the black hole when the dynamical simulation is started. This leads to radial oscillations of the star and to a more rapid orbital decay, as compared to the equilibrium case. The evolution is more rapid and the mass transfer episodes are more violent. Thus the maximum accretion rates and the peak in the gravitational radiation waveforms and luminosities are overestimated (see Table 3 and Figure 15).

We can perform a direct comparison between run A31 and the dynamical coalescence presented in Paper I for a tidally locked binary with initial mass ratio $`q=0.31`$ (see section 5.4.2 and Figures 20–23 in Paper I). Both simulations included the effect of gravitational radiation reaction on the system and were carried out with an adiabatic index $`\mathrm{\Gamma }=3`$, thus any differences are due entirely to the initial condition. In both cases there is an initial episode of intense mass transfer, with the surviving core being transferred to a higher orbit. However, in the case of the tidally locked system, it is less violent, and thus the final separation of the binary is smaller than in run A31, and the orbit is more circular. This allows the secondary episodes of mass transfer to be more intense, by about an order of magnitude (they are clearly resolved in the run shown in Paper I, despite the lower numerical resolution). The evolution is reflected in the gravitational radiation signal, with the final amplitude of the waveforms being greater by a factor of two than in run A31, and the gravitational radiation luminosity exhibiting secondary peaks at each periastron passage that are greater by approximately one order of magnitude. Additionally, no resolvable amount of mass ($`M_{ejected}10^4`$) is dynamically ejected in the tidally locked case, in stark contrast with the result of run A31.

The surviving core (for $`\mathrm{\Gamma }=3`$) can be driven below the minimum mass required for stability by the successive episodes of mass transfer. If this occurs, an explosion may take place (Page 1982; Blinnikov et al. 1984; Colpi, Shapiro & Teukolsky 1991; Sumiyoshi et al. 1998). Exploring this process clearly requires a greater level of detail in the physics input than we have at our disposal for this work, and is beyond the scope of this paper. The presence of the black hole would make this even more complicated.

All the accretion discs that form around the black hole during the coalescence are similar in structure, and by the end of our simulations, they are quite close to being azimuthally symmetric. They have masses of a few tenths of a solar mass, with maximum densities and specific internal energies of order 1011g cm-3 and 1019erg g-1 (or 10 MeV/nucelon) respectively. All of the final configurations have a low degree of baryon contamination along the rotation axis, in the regions directly above and below the black hole. It is low enough so that only modest beaming (of approximately 10) of a relativistic fireball along this axis would be required in order to avoid being stopped by the fluid in the vicinity of the black hole (Mészáros & Rees 1992, 1993). This is encouraging as far as current models of gamma ray bursts (GRBs) are concerned, and confirms our previous results (Lee & Kluźniak 1997; Kluźniak & Lee 1998), where we found that these systems appear to be good candidates for the central engines of short GRBs. The recent simulations of black hole–neutron star coalescence performed by Janka et al. (1999) using the equation of state of Lattimer & Swesty (1991) have produced essentially the same results. The coalescence of binary neutron stars may result in a very similar scenario if the central object collapses to a black hole (Ruffert & Janka 1999).

## Acknowledgments

It is a pleasure to acknowledge many helpful discussions with Włodzimierz Kluźniak, Frederic Rasio, Lars Bildsten and Maximilian Ruffert. I thank the Aspen Center for Physics for its hospitatlity. Support for this work was provided by CONACyT (27987E) and DGAPA–UNAM (IN-119998). I thank the anonymous referee for his comments and suggestions for improvements to the text.