# Coalescing Binary Neutron Stars

## Table of Contents

1. Introduction 3

1.1 Binary Star Coalescence 3

1.2 Laser Interferometers 5

1.3 General Relativity 5

2. Astrophysical Motivation and Applications 6

2.1 Gravitational Wave Astronomy 6

2.2 Gamma-Ray Bursts 7

2.3 The R-Process Problem 8

3. Calculating Gravitational Radiation Waveforms 8

3.1 The Inspiral Waveform 9

3.2 The Coalescence Waveform 9

3.3 Phase Errors in the Inspiral Waveform 10

4. Hydrodynamic Instabilities and Coalescence 11

4.1 The Stability of Binary Equilibrium Configurations 12

4.2 Mass Transfer and the Dependence on the Mass Ratio 16

4.3 Neutron Star Physics 19

5. The Stability of Compact Binaries in General Relativity 20

5.1 The ISCO in Relativistic Close Binaries 20

5.2 Binary-Induced Collapse Instability 23

5.3 The Final Fate of Mergers 25

5.4 Numerical Relativity and Future Prospects 26

6. Nonsynchronized Binaries 28

6.1 Irrotational Equilibrium Sequences 28

6.2 Coalescence of Nonsynchronized Binaries 29

References 32

## 1 Introduction

Binary neutron stars are among the most promising sources of gravitational waves for future detection by laser interferometers such as LIGO (Abramovici et al1992), VIRGO (Bradaschia et al1990), TAMA (Kuroda et al1997) and GEO (Hough 1992; Danzmann 1998). Binary neutron stars are known to exist and for some of the systems in our own galaxy (like the relativistic binary radio pulsars PSR B1913+16 and PSR B1534+12), general relativistic (hereafter GR) effects in the binary orbit have been measured to high precision (Taylor & Weisberg 1989; Stairs et al1998). With the construction of laser interferometers well underway, it is of growing urgency that we be able to predict theoretically the gravitational waveform emitted during the inspiral and the final coalescence of the two stars. Relativisitic binary systems, like binary neutron stars (NS) and binary black holes (BH) pose a fundamental challenge to theorists, as the two-body problem is one of the outstanding unsolved problems in classical GR.

### 1.1 Binary Star Coalescence

The coalescence and merging of two stars into a single object is the almost inevitable end-point of close binary evolution. Dissipation mechanisms such as friction in common gaseous envelopes, tidal dissipation, magnetic breaking, or the emission of gravitational radiation, are always present and cause the orbits of close binary systems to decay. Examples of the coalescence process for Newtonian systems that are of great current interest include the formation of blue stragglers in globular clusters from mergers of main-sequence star binaries, and the nuclear explosion or gravitational collapse of white dwarf mergers with total masses above the Chandrasekhar limit (for other examples and discussions, see, e.g., Bailyn 1993; Chen & Leonard 1993; Iben, Tutukov, & Yungelson 1996; Rasio 1995).

For most close binary systems the terminal stage of orbital decay is always hydrodynamic in nature, with the final merging of the two stars taking place on a time scale comparable to the orbital period. In many systems this is because mass transfer from one star to the other can lead to a rapid shrinking of the binary separation, which in turns accelerates the mass transfer rate, leading to an instability (for a recent discussion and references, see Soberman, Phinney, & van den Heuvel 1997). In addition to mass transfer instabilities, global hydrodynamic instabilities can drive a close binary system to rapid coalescence once the tidal interaction between the two stars becomes sufficiently strong. The existence of these global instabilities for close binary equilibrium configurations containing a compressible fluid, and their particular importance for binary NS systems, was demonstrated for the first time by the authors (Rasio & Shapiro 1992, 1994, 1995; hereafter RS1–3) using numerical hydrodynamic calculations.

Instabilities in close binary systems can also be studied using analytic methods. The classical analytic work for close binaries containing an incompressible fluid (Chandrasekhar 1969) was extended to compressible fluids in the work of Lai, Rasio, & Shapiro (1993a,b, 1994a,b,c, hereafter LRS1–5). This analytic study confirmed the existence of dynamical and secular instabilities for sufficiently close binaries containing polytropes (idealized stellar models obeying an equation of state of the form $`P=K\rho ^\mathrm{\Gamma }`$, where $`P`$ is pressure, $`\rho `$ is the rest-mass density, $`K`$ is a constant, and $`\mathrm{\Gamma }`$ is the adiabatic exponent related to the polytropic index $`n`$ according to $`\mathrm{\Gamma }=1+1/n`$). Although these simplified analytic studies can give much physical insight into difficult questions of global fluid instabilities, fully numerical calculations remain essential for establishing the stability limits of close binaries accurately and for following the nonlinear evolution of unstable systems all the way to complete coalescence. Given the absence of any underlying symmetry in the problem, these calculations must be done in 3 spatial dimensions plus time and therefore require supercomputers.

A number of different groups have now performed such calculations, using a variety of numerical methods and focusing on different aspects of the problem. Nakamura and collaborators (see Nakamura 1994 and references therein) were the first to perform 3D hydrodynamic calculations of binary NS coalescence, using a traditional Eulerian finite-difference code. Instead, RS used the Lagrangian method SPH (Smoothed Particle Hydrodynamics). They focused on determining the stability properties of initial binary models in strict hydrostatic equilibrium and calculating the emission of gravitational waves from the coalescence of unstable binaries. Many of the results of RS were later independently confirmed by New & Tohline (1997), who used completely different numerical methods but also focused on stability questions, and by Zhuge, Centrella, & McMillan (1994, 1996), who also used SPH. Zhuge et al(1996) also explored in detail the dependence of the gravitational wave signals on the initial NS spins. Davies et al(1994) and Ruffert et al(1996, 1997) have incorporated a treatment of the nuclear physics in their hydrodynamic calculations (done using SPH and PPM codes, respectively), motivated by cosmological models of gamma-ray bursts (see Sec. 2.2).

In GR, strong-field gravity between the masses in a binary system is alone sufficient to drive a close circular orbit unstable. In close NS binaries, GR effects combine nonlinearly with Newtonian tidal effects so that close binary configurations can become dynamically unstable earlier during the inspiral phase (i.e., at larger binary separation and lower orbital frequency) than predicted by Newtonian hydrodynamics alone. The combined effects of relativity and hydrodynamics on the stability of close compact binaries have only very recently begun to be studied. Preliminary results have been obtained using both analytic approximations (basically, post-Newtonian generalizations of LRS; see Lai 1996; Taniguchi & Nakamura 1996; Lai & Wiseman 1997; Lombardi, Rasio, & Shapiro 1997; Taniguchi & Shibata 1997; Shibata & Taniguchi 1997), as well as numerical hydrodynamics calculations in 3D incorporating simplified treatments of relativistic effects (Shibata 1996; Baumgarte et al1997; Baumgarte et al1998a,b; Mathews & Wilson 1997; Shibata, Baumgarte, & Shapiro 1998; Wang, Swesty, & Calder 1998). Several groups, including a NASA Grand Challenge team (Seidel 1998; Swesty & Saylor 1997), are working on a fully relativistic calculation of the final coalescence, combining the techniques of numerical relativity and numerical hydrodynamics in 3D.

### 1.2 Laser Interferometers

It is useful to recall some of the vital statistics of the LIGO/VIRGO/GEO/TAMA network now under construction (see Thorne 1996 for an excellent review and references). It consists of earth-based, kilometer-scale laser interferometers most sensitive to waves in the $`1010^3`$Hz band. The expected rms noise level has an amplitude $`h_{rms}<10^{22}`$. The most promising sources for such detectors are NS–NS, NS–BH and BH–BH coalescing binaries. The event rates are highly uncertain but astronomers (e.g., Phinney 1991; Narayan, Piran & Shemi 1991) estimate that in the case of NS–NS binaries, which are observed in our own galaxy as binary radio pulsars, the rate may be roughly $`3\mathrm{y}\mathrm{r}^1(\mathrm{distance}/200\mathrm{Mpc})^3`$. For binaries containing black holes, the typical BH mass range in the frequency range of interest is $`2300M_{}`$. For typical NS–NS binaries, the total inspiral timescale across the detectable frequency band is approximately 15 mins. During this time the number of cycles of gravitational waves, $`𝒩_{cyc}`$, is approximately 16,000.

Although much of the current theoretical focus is directed toward LIGO-type experiments, other detectors that may come on-line in the future will also be important. For example, LISA is a proposed space-based, 5 million-kilometer interferometer that will be placed in heliocentric orbit (see, e.g., Danzmann 1998). The relevant frequency band for LISA is $`10^41`$Hz. The most promising sources in this band include short-period, galactic binaries of all types (main sequence binaries; white dwarf-white dwarf binaries, and binaries containing neutron stars and stellar-mass black holes) as well as supermassive BH–BH binaries. The typical black hole mass in a detectable BH–BH binary must be between $`10^310^8M_{}`$, where the upper mass limit is set by the lower bound on the observable frequency.

### 1.3 General Relativity

Solving the binary coalescence problem will ultimately require the full machinery of general relativity. Indeed, many of the key issues cannot even be raised in the context of Newtonian or even post-Newtonian gravitation.

Consider, for example, the recent controversial claim by Wilson, Mathews and Marronetti (Wilson & Mathews 1995; Wilson et al1996; Marronetti et al1998) that massive neutron stars in close binaries can collapse to black holes prior to merger. Catastrophic collapse of equilibrium fluids to black holes is a consequence of the nonlinear nature of Einstein’s field equations and can only be addressed in full GR. Resolving the issue of neutron star collapse prior to merger has huge consequences for predictions of gravitational waveforms from neutron star binaries. In addition, if the neutron stars do undergo collapse, their final coalescence cannot serve as a source of gamma-rays.

There are other aspects of the inspiral problem that require a fully relativistic treatment, even for a qualitative understanding. For example, if the stars do not undergo collapse prior to coalescence, their combined mass is likely to exceed the maximum allowed mass of a cold, rotating star upon merger. In this case the merged remnant must ultimately undergo collapse to a black hole. But it is not clear whether this final collapse proceeds immediately, on a dynamical timescale (ms), or quasi-statically, on a neutrino dissipation timescale (secs). The latter is possible since the merged remnant may be hot, following shock heating, and the thermal component of the pressure may be adequate to keep the star in quasi-equilibrium until neutrinos carry off this thermal energy and with it the thermal pressure support against collapse (Baumgarte & Shapiro 1998). Moreover, it is by no means clear how much angular momentum the rotating remnant will possess at the onset of collapse or what the final fate of the system will be if the angular momentum of the remnant exceeds the maximum allowed value for a Kerr black hole, $`J/M^2=1`$ (We adopt units such that $`G=c=1`$ throughout this paper unless otherwise specified). Will the excess angular momentum be radiated away or ejected via a circumstellar ring or torus? These issues have crucial observational implications and can only be addressed by simulations performed in full GR.

## 2 Astrophysical Motivation and Applications

### 2.1 Gravitational Wave Astronomy

Coalescing compact binaries are very strong sources of gravitational radiation that are expected to become directly detectable with the new generation of laser interferometers now under construction (see Sec. 1). In addition to providing a major new confirmation of Einstein’s theory of general relativity, including the first direct proof of the existence of black holes (Flanagan & Hughes 1998,a,b; Lipunov et al1997), the detection of gravitational waves from coalescing binaries at cosmological distances could provide accurate independent measurements of the Hubble constant and mean density of the Universe (Schutz 1986; Chernoff & Finn 1993; Marković 1993). For a recent review on the detection and sources of gravitational radiation, see Thorne (1996).

Expected rates of NS binary coalescence in the Universe, as well as expected event rates in forthcoming laser interferometers, have now been calculated by many groups. Although there is some disparity between various published results, the estimated rates are generally encouraging. Simple statistical arguments based on the observed local population of binary radio pulsars with probable NS companions lead to an estimate of the rate of NS binary coalescence in the Universe of order $`10^7`$yr$`^1`$Mpc<sup>-3</sup> (Narayan et al1991; Phinney 1991). In contrast, theoretical models of the binary star population in our Galaxy suggest that the NS binary coalescence rate may be much higher, $`>10^6`$yr$`^1`$Mpc<sup>-3</sup> (Tutukov & Yungelson 1993; see also the more recent studies by Portegies Zwart & Spreeuw 1996 and Lipunov et al1998).

Finn & Chernoff (1993) predicted that an advanced LIGO detector could observe as many as 70 NS merger events per year. This number corresponds to a Galactic NS merger rate $`R10^6\mathrm{yr}^1`$ derived from radio pulsar surveys. More recently, however, van den Heuvel & Lorimer (1996) revised this number to $`R0.8\times 10^5\mathrm{yr}^1`$, using the latest galactic pulsar population model of Curran & Lorimer (1995). This value is consistent with the upper limit of $`10^5\mathrm{yr}^1`$ for the Galactic binary NS birth rate derived by Bailes (1996) on the basis of very general statistical considerations about pulsars.

Near the end of the inspiral, when the binary separation becomes comparable to the stellar radii, hydrodynamic effects become important and the character of the waveforms will change. Special purpose narrow-band detectors that can sweep up frequency in real time will be used to try to catch the corresponding final few cycles of gravitational waves (Meers 1988; Strain & Meers 1991; Danzmann 1998). In this terminal phase of the coalescence, the waveforms contain information not just about the effects of GR, but also about the internal structure of the stars and the nuclear equation of state (hereafter EOS) at high density. Extracting this information from observed waveforms, however, requires detailed theoretical knowledge about all relevant hydrodynamic processes. This question is discussed in more detail in Section 4 below.

### 2.2 Gamma-Ray Bursts

Many theoretical models of gamma-ray bursts (GRB) have postulated that the energy source for the bursts could be coalescing compact (NS–NS or NS–BH) binaries at cosmological distances (Paczyński 1986; Eichler et al1989; Narayan, Paczyński, & Piran 1992). The isotropic angular distribution of the bursts detected by the BATSE experiment on the Compton GRO satellite (Meegan et al1992) strongly suggests a cosmological origin, as does the distribution of number versus intensity of the bursts. In addition, the rate of GRBs detected by BATSE, of order one per day, is in rough agreement with theoretical predictions for the rate of NS binary coalescence in the Universe (cf. above). During the past two years, the first X-ray “afterglows,” as well as radio and optical counterparts of several GRBs have been observed after the burst positions were measured accurately with the BeppoSAX satellite (e.g., Costa et al1997). These observations have provided very strong additional evidence for a cosmological origin of the bursts. Most importantly, the recent detections of the optical counterparts of several GRBs at high redshifts ($`Z=0.84`$ for GRB 970508, $`Z=3.4`$ for GRB 971214; see Metzger et al1997 and Kulkarni et al1998) have firmly established that at least some gamma-ray bursts originate at cosmological distances.

To model the gamma-ray emission realistically, the complete hydrodynamic and nuclear evolution during the final merging of the two NS, especially in the outermost, low-density regions of the merger, must be understood in detail. This is far more challenging than understanding the emission of gravitational waves, which is mostly sensitive to the bulk motion of the fluid, but is totally insensitive to nuclear processes taking place in low-density regions. Numerical calculations of NS binary coalescence including some treatment of the nuclear physics have been performed in Newtonian theory by Davies et al(1994; see also Rosswog et al1998a,b) and Ruffert et al(1996, 1997). The most recent results from these calculations indicate that, even under the most favorable conditions, the energy provided by $`\nu \overline{\nu }`$ annihilation during the merger is too small by at least an order of magnitude, and more probably two or three orders of magnitude, to power typical gamma-ray bursts at cosmological distances (Janka & Ruffert 1996). The discrepancy has now become even worse given the higher energies required to power bursts at some of the observed high redshifts ($`10^{54}`$erg for isotropic emission in the case of GRB 971214). However, with sufficient beaming of the gamma ray emission, scenarios in which the merger leads to the formation of a rapidly rotating black hole surrounded by a torus of debris, and where the energy of the burst comes from either the binding energy of the debris, or the spin energy of the black hole, are still viable (Mészáros, Rees, & Wijers 1998).

### 2.3 The R-Process Problem

Recent calculations have raised doubts on the ability of supernovae to produce r-process nuclei in the correct amounts (e.g., Meyer & Brown 1997). Instead, decompressed nuclear matter ejected during binary NS coalescence, or during the tidal disruption of a NS by a BH, may provide a good alternative or supplementary site for the r-process (Symbalisty & Schramm 1992; Eichler et al1989; Rosswog et al1998a,b).

The recent SPH calculations by Rosswog et al(1998b) suggest that the amount of mass ejected during binary NS coalescence may be sufficient for an explanation of the observed r-process abundances. Their preliminary abundance calculations show that practically all the material is subject to r-process conditions. The calculated abundance patterns can reproduce the basic features of the solar r-process abundances very well, including the peak near A=195, which is obtained without any tuning of the initial entropies. Thus, it is possible that all the observed r-process material could be explained by mass ejection during neutron star mergers.

## 3 Calculating Gravitational Radiation Waveforms

At present, we do not possess a single, unified prescription for calculating gravitational waveforms over all the regimes and all the corresponding bands of detectable frequencies from such events. Instead, we must be crafty in breaking up the coalescence into several distinct epochs and corresponding frequency bands and employing appropriate theoretical tools to investigate each epoch separately. One of our immediate theoretical goals is to construct a smooth, self-consistent join between the different solutions for the different epochs. Ultimately, we may succeed in formulating a single computational approach that is capable by itself of tracking the entire binary coalescence and merger and determining the waveform over all frequency bands. But for now we must content ourselves with calculating waveforms by any means possible – by any means necessary!

Gravitational waveforms from coalescing compact binaries may be conveniently divided into two main pieces (Cutler et al1993). The inspiral waveform is the low-frequency component emitted early on, before tidal distortions of the stars become important. The coalescence waveform is the high frequency component emitted at the end, during the epoch of distortion, tidal disruption and/or merger. Existing theoretical machinery for handling the separate epochs differs considerably.

### 3.1 The Inspiral Waveform

Most recent calculations of the gravitational radiation waveforms from coalescing binaries have focused on the signal emitted during the last few thousand orbits, as the frequency sweeps upward from about 10 Hz to $`300`$Hz. The waveforms in this regime can be calculated fairly accurately by performing high-order post-Newtonian (hereafter PN) expansions of the equations of motion for two point masses (Lincoln & Will 1990; Junker & Schäfer 1992; Kidder, Will, & Wiseman 1992; Wiseman 1993; Will 1994; Blanchet et al1996). High accuracy is essential here because the observed signals will be matched against theoretical templates. Since the templates must cover $`10^310^4`$ orbits, a phase error as small as $`10^4`$ could in principle prevent detection (Cutler et al1993; Cutler & Flanagan 1994; Finn & Chernoff 1993).

The PN formalism consists of a series expansion in the parameter $`ϵM/rv^2`$, where $`M`$ is the mass of the binary, $`r`$ is the separation and $`v`$ is the orbital velocity. This parameter is small whenever the gravitational field is weak and the velocity is slow. In this formalism, which is essentially analytic, the stars are treated as point masses. The aim of the PN analysis is to compute to $`𝒪[(v/c)^{11}]`$ in order that theoretical waveforms be sufficiently free of systematic errors to be reliable as templates against which the LIGO/VIRGO observational data can be compared (Cutler & Flanagan 1994). For further discussion of the PN formalism and references, see Blanchet & Damour (1992), Kidder, Will & Wiseman (1993), Apostolatos et al(1994), Blanchet et al(1995) and Will & Wiseman (1996).

### 3.2 The Coalescence Waveform

The coalescence waveform is influenced by finite-size effects, like hydrodynamics in the case of neutron stars, and by tidal distortions. For binary neutron stars, many aspects of coalescence can be understood by solving the Newtonian equations of hydrodynamics while treating the gravitational radiation as a perturbation in the quadrupole approximation. Such an analysis is only valid when the two inequalities, $`ϵ1`$ and $`M/R1`$ are both satisfied. Here R is the neutron star radius. Newtonian treatments of the coalescence waveform come in two forms: numerical hydrodynamic simulations in 3D and analytic analyses based on triaxial ellipsoid models of the interacting stars. The ellipsoidal treatments can handle the influence of tidal distortion and internal fluid motions and spin, but not the final merger and coalescence. For a detailed treatment and references, see Chandrasekhar (1969), Carter & Luminet (1985), Kochanek (1992) and LRS. Numerical simulations are required to treat the complicated hydrodynamic interaction with ejection of mass and shock dissipation, which usually accompany the merger (see, e.g., Oohara & Nakamura 1989; RS; Davies et al1994; Zughe et al1994; Ruffert et al1995).

Fully relativistic calculations are required for quantitatively reliable coalescence waveforms. They are also required to determine those qualitative features of the final merger which can only result from strong-field effects (e.g., catastrophic collapse of merging neutron stars to a black hole). These calculations treat Einstein’s equations numerically in 3+1 dimensions without approximation. In the case of neutron stars, the equations of relativistic hydrodynamics must be solved together with Einstein’s field equations. For earlier work in this area, see the articles in Smarr (1979) and Evans, Finn & Hobill (1988); for recent progress see Matzner et al1995 and Wilson & Mathews 1995, and references therein.

### 3.3 Phase Errors in the Inspiral Waveform

Measuring the binary parameters by gravitational wave observations is accomplished by integrating the observed signal against theoretical templates (Cutler et al1993). For this purpose it is necessary that the signal and template remain in phase with each other within a fraction of a cycle ($`\delta 𝒩_{cyc}<0.1`$) as the signal sweeps through the detector’s frequency band. To leading order we may treat the system as a point-mass, nearly-circular Newtonian binary spiraling slowly inward due to the emission of quadrupole gravitational radiation. In this limit the number of cycles spent sweeping through a logarithmic interval of frequency $`f`$ is

$$\left(\frac{d𝒩_{cyc}}{d\mathrm{ln}f}\right)_0=\frac{5}{96\pi }\frac{1}{M_c^{5/3}(\pi f)^{5/3}},$$
(1)   
where the “chirp mass” $`M_c`$ is given by $`M_c\mu ^{3/5}M^{2/5}`$. Here $`\mu `$ is the reduced mass and $`M`$ is the total mass of the binary. It is expected that LIGO/VIRGO measurements will be able to determine the chirp mass to within 0.04 per cent for a NS–NS binary and to within 0.3 per cent for a system containing at least one BH (Thorne 1996).

The PN formalism can be used to determine corrections to eq. (1) arising from PN contributions to the binary orbit. For example, suppose one of the stars has a spin $`𝐒`$ inclined at an angle $`i`$ to the normal direction to the orbital plane. This spin induces a gravitomagnetic field which modifies the orbit of the companion. In addition, the wave emission rate, which determines the inspiral velocity, is augmented above the value due to the familiar time-changing quadrupole mass moment by an additional contribution from the time-changing quadrupole current moment. The result is easily shown to yield a “correction” to the Newtonian binary phase (eq. 1) given by

$$\frac{d𝒩_{cyc}}{d\mathrm{ln}f}=\left(\frac{d𝒩_{cyc}}{d\mathrm{ln}f}\right)_0\left[1+\frac{113}{12}\frac{S}{M^2}x^{3/2}cosi\right],$$
(2)   
where $`x(\pi Mf)^{2/3}M/r`$ and where we have assumed that the mass of the spinning star is much greater than that of the companion. The frequency dependence of the correction term enables us in principle to distinguish this spin contribution from the Newtonian piece. In practice, it turns out that we may need to know independently the value of the spin in order to determine reliably the reduced mass $`\mu `$ (and thereby $`M`$, and the individual masses, since we already know $`M_c`$ from the Newtonian part of eq. 2). If we somehow know that the spin is small, we can determine $`\mu `$ to roughly 1% for NS–NS and NS–BH binaries and 3% for BH–BH binaries (Thorne 1996). Not knowing the value of the spin worsens the accuracy of $`\mu `$ considerably, but this may be improved if wave modulations due to spin-induced Lens-Thirring precession of the orbit are incorporated (Apostolatos et al1994). This example illustrates how the PN formalism may be used to do classical stellar spectroscopy on binary systems containing compact stars.

Models based on Newtonian compressible ellipsoids can be used to analyze finite-size effects that lead to additional corrections to the phase of a NS–NS binary inspiral waveform (LRS3). Consider for definiteness two identical $`1.4M_{}`$ neutron stars, each with radius $`R/M=5`$ and supported by a stiff polytropic equation of state with adiabatic index $`\mathrm{\Gamma }=3`$. Track their orbit as they spiral inward from a separation $`r_i=70R`$ to $`r_f=5R`$, corresponding to a sweep over wave frequency from $`f_i=10`$ Hz to $`f_f=522`$ Hz (recall for Keplerian motion, $`fr^{3/2}`$). To lowest Newtonian order, the total number of wave cycles emitted as the stars sweep through this frequency band is 16,098. If the two stars have zero spin, then the main hydrodynamic correction to the point-mass Newtonian result is due to the static Newtonian quadrupole interaction induced by the tidal field. The change in the number of cycles varies like $`\delta 𝒩_{cyc}^{(I)}r^{5/2}f^{5/3}`$ and therefore arises chiefly at large $`f`$ (small $`r`$). Sweeping through the entire frequency band results in a small change $`\delta 𝒩_{cyc}^{(I)}0.3`$; in the low frequency band from $`10`$Hz to $`300`$Hz, the change is only 0.1. Such a small change probably can be neglected in designing low-$`f`$ wave templates.

Suppose instead that each NS has an intrinsic spin. In this case $`\delta 𝒩_{cyc}^{(S)}r^{1/2}f^{1/3}`$ and the change occurs chiefly at low $`f`$ (large $`r`$). Now the quadrupole moments of the stars are induced by spin as well as by tidal fields. The change in the number of wave cycles as the orbit decays to $`r_f`$ is $`\delta 𝒩_{cyc}^{(S)}9/P_{ms}^2`$, where $`P_{ms}`$ is the spin period in msec. Hence for rapidly spinning NS’s with $`P_{ms}<9`$, the effect is potentially important and must be taken into account in theoretical templates.

Unlike many binaries consisting of ordinary stars, NS binaries are not expected to be corotating (synchronous) at close separation, because the viscosities required to achieve synchronous behavior are implausibly large (Kochanek 1992; Bildsten and Cutler 1992; LRS3). Were this otherwise, the resulting corrections on the inspiral waveform phase evolution would be enormous and would dominate the low-$`f`$ phase correction: $`\delta 𝒩_{cyc}^{(SS)}15`$ in orbiting from $`r=r_i`$ to $`r=r_f`$. See Section 6 for a discussion of the final coalescence for nonsynchronous binaries.

## 4 Hydrodynamic Instabilities and Coalescence

Newtonian hydrodynamic calculations in 3D yield considerable insight into the coalescence process. These calculations also serve as benchmarks for future relativistic codes in the weak-field, slow-velocity limit of GR, applicable whenever $`R/M10`$. This section summarizes some of the important physical effects revealed by these Newtonian simulations.

### 4.1 The Stability of Binary Equilibrium Configurations

Hydrostatic equilibrium configurations for binary systems with sufficiently close components can become dynamically unstable (Chandrasekhar 1975; Tassoul 1975). The physical nature of this instability is common to all binary interaction potentials that are sufficiently steeper than $`1/r`$ (see, e.g., Goldstein 1980, §3.6). It is analogous to the familiar instability of test particles in circular orbits sufficiently close to a black hole (Shapiro & Teukolsky 1983, §12.4). Here, however, it is the tidal interaction that is responsible for the steepening of the effective interaction potential between the two stars and for the destabilization of the circular orbit (LRS3). The tidal interaction exists of course already in Newtonian gravity and the instability is therefore present even in the absence of relativistic effects. For sufficiently compact binaries, however, the combined effects of relativity and hydrodynamics lead to an even stronger tendency towards dynamical instability (see §5).

The stability properties of close NS binaries depend sensitively on the NS EOS. Close binaries containing NS with stiff EOS (adiabatic exponent $`\mathrm{\Gamma }>2`$ if $`P=K\rho ^\mathrm{\Gamma }`$, where $`P`$ is pressure and $`\rho `$ is density) are particularly susceptible to a dynamical instability. This is because tidal effects are stronger for stars containing a less compressible fluid (i.e., for larger $`\mathrm{\Gamma }`$). As the dynamical stability limit is approached, the secular orbital decay driven by gravitational wave emission can be dramatically accelerated (LRS2, LRS3). The two stars then plunge rapidly toward each other, and merge together into a single object in just a few rotation periods. This dynamical instability was first identified in RS1, where the evolution of Newtonian binary equilibrium configurations was calculated for two identical polytropes with $`\mathrm{\Gamma }=2`$. It was found that when $`r<3R`$ ($`r`$ is the binary separation and $`R`$ the radius of an unperturbed NS), the orbit becomes unstable to radial perturbations and the two stars undergo rapid coalescence. For $`r>3R`$, the system could be evolved dynamically for many orbital periods without showing any sign of orbital evolution (in the absence of dissipation). Many of the results derived in RS and LRS concerning the stability properties of NS binaries have been confirmed recently in completely independent work by New & Tohline (1997) and by Zhuge, Centrella, & McMillan (1996). New & Tohline (1997) used completely different numerical methods (a combination of a 3D Self-Consistent Field code for constructing equilibrium configurations and a grid-based Eulerian code for following the dynamical evolution of the binaries), while Zhuge et al(1996) used SPH, as did RS.

*

The dynamical evolution of an unstable, initially synchronized (i.e., rigidly rotating) binary containing two identical stars can be described typically as follows (Fig. 1). During the initial, linear stage of the instability, the two stars approach each other and come into contact after about one orbital revolution. In the corotating frame of the binary, the relative velocity remains very subsonic, so that the evolution is adiabatic at this stage. This is in sharp contrast to the case of a head-on collision between two stars on a free-fall, radial orbit, where shocks are very important for the dynamics (RS1). Here the stars are constantly being held back by a (slowly receding) centrifugal barrier, and the merging, although dynamical, is much more gentle. After typically two orbital revolutions the innermost cores of the two stars have merged and the system resembles a single, very elongated ellipsoid. At this point a secondary instability occurs: mass shedding sets in rather abruptly. Material is ejected through the outer Lagrange points of the effective potential and spirals out rapidly. In the final stage, the spiral arms widen and merge together. The relative radial velocities of neighboring arms as they merge are supersonic, leading to some shock-heating and dissipation. As a result, a hot, nearly axisymmetric rotating halo forms around the central dense core. The halo contains about 20% of the total mass and the rotation profile is close to a pseudo-barotrope (Tassoul 1978, §4.3), with the angular velocity decreasing as a power-law $`\mathrm{\Omega }\varpi ^\nu `$ where $`\nu <2`$ and $`\varpi `$ is the distance to the rotation axis (RS1). The core is rotating uniformly near breakup speed and contains about 80% of the mass still in a cold, degenerate state. If the initial NS had masses close to $`1.4M_{}`$, then most recent stiff EOS would predict that the final merged configuration is still stable and will not immediately collapse to a black hole, although it might ultimately collapse to a black hole as it continues to lose angular momentum (see Cook, Shapiro, & Teukolsky 1994).

The emission of gravitational radiation during dynamical coalescence can be calculated perturbatively using the quadrupole approximation (RS1). Both the frequency and amplitude of the emission peak somewhere during the final dynamical coalescence, typically just before the onset of mass shedding. Immediately after the peak, the amplitude drops abruptly as the system evolves towards a more axially symmetric state. For an initially synchronized binary containing two identical polytropes, the properties of the waves near the end of the coalescence depend very sensitively on the stiffness of the EOS (Fig. 2).

When $`\mathrm{\Gamma }<\mathrm{\Gamma }_{crit}`$, with $`\mathrm{\Gamma }_{crit}2.3`$, the final merged configuration is perfectly axisymmetric. Indeed, a Newtonian polytropic fluid with $`\mathrm{\Gamma }<2.3`$ (polytropic index $`n>0.8`$) cannot sustain a nonaxisymmetric, uniformly rotating configuration in equilibrium (see, e.g., Tassoul 1978, §10.3). As a result, the amplitude of the waves drops to zero in just a few periods (RS1). In contrast, when $`\mathrm{\Gamma }>\mathrm{\Gamma }_{crit}`$, the dense central core of the final configuration remains triaxial (its structure is basically that of a compressible Jacobi ellipsoid; cf. LRS1) and therefore it continues to radiate gravitational waves. The amplitude of the waves first drops quickly to a nonzero value and then decays more slowly as gravitational waves continue to carry angular momentum away from the central core (RS2). Because realistic NS EOS have effective $`\mathrm{\Gamma }`$ values precisely in the range 2–3 (LRS3), i.e., close to $`\mathrm{\Gamma }_{crit}2.3`$, a simple determination of the absence or presence of persisting gravitational radiation after the coalescence (i.e., after the peak in the emission) could place a strong constraint on the stiffness of the EOS. General relativity is likely to play an important quantitative role; for example, the critical Newtonian value of polytropic index for the onset of the bar-mode instability is increased to $`n=1.3`$ in GR (Stergioulas & Friedman 1998).

### 4.2 Mass Transfer and the Dependence on the Mass Ratio

Clark & Eardley (1977) suggested that secular, stable mass transfer from one NS to another could last for hundreds of orbital revolutions before the lighter star is tidally disrupted. Such an episode of stable mass transfer would be accompanied by a secular increase of the orbital separation. Thus if stable mass transfer could indeed occur, a characteristic “reversed chirp” would be observed in the gravitational wave signal at the end of the inspiral phase (Jaranowski & Krolak 1992).

The question was later reexamined by Kochanek (1992) and Bildsten & Cutler (1992), who both argued against the possibility of stable mass transfer on the basis that very large mass transfer rates and extreme mass ratios would be required. Moreover, in LRS3 it was pointed out that mass transfer has in fact little importance for most NS binaries (except perhaps those containing a very low-mass NS). This is because for $`\mathrm{\Gamma }>2`$, dynamical instability always arises before the Roche limit along a sequence of binary configurations with decreasing separation $`r`$. Therefore, by the time mass transfer begins, the system is already in a state of dynamical coalescence and it can no longer remain in a nearly circular orbit. Thus stable mass transfer from one NS to another appears impossible.

In RS2 a complete dynamical calculation was presented for a system containing two polytropes with $`\mathrm{\Gamma }=3`$ and a mass ratio $`q=0.85`$. This value corresponds to what was at the time the most likely mass ratio for the binary pulsar PSR B2303+46 (Thorsett et al1993) and represented the largest observed departure from $`q=1`$ in any known binary pulsar with likely NS companion. The latest observations of PSR B2303+46, however, give a most likely mass ratio $`q=1.30/1.34=0.97`$ (Thorsett & Chakrabarty 1998). For comparison, $`q=1.386/1.442=0.96`$ in PSR B1913+16 (Taylor & Weisberg 1989), $`q=1.349/1.363=0.99`$ for PSR B2127+11C (Deich & Kulkarni 1996), and $`q=1.339/1.339=1`$ for PSR B1534+12 (Wolszczan 1991; Thorsett & Chakrabarty 1998). Neutron star masses derived from observations of binary radio pulsars are all consistent with a remarkably narrow underlying Gaussian mass distribution with $`M_{\mathrm{NS}}=1.35\pm 0.04M_{}`$ (Thorsett & Chakrabarty 1998).

However, it cannot be excluded that other binary NS systems (that may not be observable as binary pulsars) could contain stars with significantly different masses. For a system with $`q=0.85`$, RS2 found that the dynamical stability limit is at $`r/R2.95`$, whereas the Roche limit is at $`r/R2.85`$. The dynamical evolution turns out to be dramatically different from that of a system with $`q=1`$. The Roche limit is quickly reached while the system is still in the linear stage of growth of the instability. Dynamical mass transfer from the less massive to the more massive star begins within the first orbital revolution. Because of the proximity of the two components, the fluid acquires very little velocity as it slides down from the inner Lagrange point to the surface of the other star. As a result, relative velocities of fluid particles remain largely subsonic and the coalescence proceeds quasi-adiabatically, just as in the $`q=1`$ case. In fact, the mass transfer appears to have essentially no effect on the dynamical evolution. After about two orbital revolutions the smaller-mass star undergoes complete tidal disruption. Most of its material is quickly spread on top of the more massive star, while a small fraction of the mass is ejected from the outer Lagrange point and forms a single-arm spiral outflow. The more massive star, however, remains little perturbed during the entire evolution and simply becomes the inner core of the merged configuration. This type of dynamical evolution, which is probably typical for the final merging of two NS with slightly different masses, is illustrated in Fig. 3.

The dependence of the peak amplitude $`h_{max}`$ of gravitational waves on the mass ratio $`q`$ appears to be very strong, and nontrivial. In RS2 an approximate scaling $`h_{max}q^2`$ was derived. This is very different from the scaling obtained for a detached binary system with a given binary separation. In particular, for two point masses in a circular orbit with separation $`r`$ the result would be $`h\mathrm{\Omega }^2\mu r^2`$, where $`\mathrm{\Omega }^2=G(M+M^{})/r^3`$ and $`\mu =MM^{}/(M+M^{})`$. At constant $`r`$, this gives $`hq`$. This linear scaling is obeyed (only approximately, because of finite-size effects) by the wave amplitudes of the various systems at the onset of dynamical instability. For determining the maximum amplitude, however, hydrodynamics plays an essential role. In a system with $`q1`$, the more massive star tends to play a far less active role in the hydrodynamics and, as a result, there is a rapid suppression of the radiation efficiency as $`q`$ departs even slightly from unity. For the peak luminosity of gravitational radiation RS found approximately $`L_{max}q^6`$. Again, this is a much steeper dependence than one would expect based on a simple point-mass estimate, which gives $`Lq^2(1+q)`$ at constant $`r`$. The results of RS are all for initially synchronized binaries, but very similar results have been obtained more recently by Zhuge et al(1996) for binaries containing initially nonspinning stars with unequal masses.

Little is known about the stability of the mass transfer from a NS to a BH. The first 3D hydrodynamic calculations of the coalescence process for NS-BH binaries were performed recently by Lee & Kluzniak (1998) using a Newtonian SPH code. For all mass ratios in the range of about $`13`$, they find that, after a brief episode of mass transfer, the system stabilizes with a remnant NS core surviving in orbit around the more massive BH. This is qualitatively similar to the results obtained in RS2 for a NS-NS binary with mass ratio $`0.5`$. However, for NS-BH binaries, even in the case of a very stiff NS EOS, one expects relativistic effects to be very important, since the Roche limit radius and the ISCO radius around the BH are very close to each other for any BH more massive than the NS. Therefore the results of purely Newtonian calculations for BH-NS binaries may not even provide a qualitatively correct picture of the final merging.

### 4.3 Neutron Star Physics

The most important parameter that enters into quantitative estimates of the gravitational wave emission during the final coalescence is the ratio $`M/R`$ for a NS. In particular, for two identical point masses we know that the wave amplitude $`h`$ obeys $`(r_O/M)h(M/R)`$, where $`r_O`$ is the distance to the observer, and the total luminosity $`L(M/R)^5`$. Similarly the wave frequency $`f_{max}`$ during final merging should satisfy approximately $`f_{max}(M/R)^{3/2}`$ since it is roughly twice the Keplerian frequency for two NS in contact (binary separation $`r23R`$). Thus one expects that any quantitative measurement of the emission near maximum should lead to a direct determination of the NS radius $`R`$, assuming that the mass $`M`$ has already been determined from the low-frequency inspiral waveform (Cutler & Flanagan 1994). Most current NS EOS give $`M/R0.1`$, with $`R10\mathrm{km}`$ nearly independent of the mass in the range $`0.8M_{}<M<1.5M_{}`$ (see, e.g., Baym 1991; Cook et al1994; LRS3; Akmal, Pandharipande and Ravenhall 1998).

However, the details of the hydrodynamics also enter into this determination. The importance of hydrodynamic effects introduces an explicit dependence of all wave properties on the EOS (which we represent here by a single dimensionless parameter $`\mathrm{\Gamma }`$), and on the mass ratio $`q`$. If relativistic effects were taken into account for the hydrodynamics itself, an additional, nontrivial dependence on $`M/R`$ would also be present. This can be written conceptually as

$`\left({\displaystyle \frac{r_O}{M}}\right)h_{max}`$ $``$ $`(q,\mathrm{\Gamma },M/R)\times \left({\displaystyle \frac{M}{R}}\right)`$ (3)   
$`{\displaystyle \frac{L_{max}}{L_o}}`$ $``$ $`(q,\mathrm{\Gamma },M/R)\times \left({\displaystyle \frac{M}{R}}\right)^5`$ (4)   
Combining all the results of RS, we can write, in the limit where $`M/R0`$ and for $`q`$ not too far from unity,

$$(q,\mathrm{\Gamma },M/R)2.2q^2(q,\mathrm{\Gamma },M/R)0.5q^6,$$
(5)   
essentially independent of $`\mathrm{\Gamma }`$ in the range $`\mathrm{\Gamma }2`$–3 (RS2).

The results of RS were for the case of synchronized spins. Recently Zhuge et al(1996) have performed calculations for nonsynchronized binaries and obtained very similar results (but see §6 below). For example, for the coalescence of two nonspinning stars with $`q=1`$ they found $`1.92.3`$ and $`0.290.59`$, where the range of values corresponds to varying $`\mathrm{\Gamma }`$ between $`5/3`$ and 3. Note that the calculations of Zhuge et al(1996) included an approximate treatment of PN effects by setting up an initial inspiral trajectory for two NS of mass $`M=1.4M_{}`$ and radius in the range $`R=1015`$km. Varying the radius of the stars in this range appears to leave the coefficients $``$ and $``$ practically unchanged within their approximation. Zhuge et al(1994, 1996) also compute frequency spectra for the gravitational wave emission and discuss various ways of defining precisely the characteristic frequency $`f_{max}`$.

Gravitational wave emission from colliding neutron stars (which may resemble coalescing NS binaries in the highly relativistic limit where a very large radial infall velocity develops prior to final merging) have been calculated recently by RS1 and Centrella & McMillan (1993) using SPH, and by Ruffert & Janka (1998) using a grid-based (PPM) code. However, even for the simplest case of head-on (axisymmetric) collisions in the Newtonian limit, the full dependence of the waveforms on the NS EOS and on the mass ratio has yet to be explored.

## 5 The Stability of Compact Binaries in General Relativity

### 5.1 The ISCO in Relativistic Close Binaries

Over the last two years, various efforts have started to calculate the stability limits for NS binaries including both hydrodynamic finite-size (tidal) effects and relativistic effects. Note that, strictly speaking, equilibrium circular orbits do not exist in GR because of the emission of gravitational waves.However, outside the innermost stable circular orbit (ISCO), the timescale for orbital decay by radiation is much longer than the orbital period, so that the binary can be considered to be in “quasiequilibrium”. This fact allows one to neglect both gravitational waves and wave-induced deviations from a circular orbit to a very good approximation outside the ISCO. Accordingly, the stability of quasi-circular orbits can be studied in the framework of GR by truncating the radiation-reaction terms in a PN expansion of the equations of motion (Lincoln & Will 1990; Kidder et al1992; Will 1994). Alternatively, one can solve a subset of the full nonlinear Einstein equations numerically in the $`3+1`$ formalism on time slices with a spatial 3-metric chosen to be conformally flat (Wilson & Mathews 1989, 1995; Wilson et al1996; Baumgarte et al1997). In the spirit of the York-Lichnerowicz conformal decomposition, which separates radiative variables from nonradiative ones, (Lichnerowicz 1944; York 1971) such a choice is believed to effectively minimize the gravitational wave content of space-time. In addition, one can set the time-derivatives of many of the metric functions equal to zero in the comoving frame, forcing the solution to be approximately time-independent in that frame. The field equations then reduce to a set of coupled elliptic equations (for the $`3+1`$ lapse and shift functions and the conformal factor); see §5.1.2 for more detailed discussion.

Several groups are now working on PN generalizations of the semi-analytic Newtonian treatment of LRS based on ellipsoids. Taniguchi & Nakamura (1996) consider NS–BH binaries and adopt a modified version of the pseudo-Newtonian potential of Paczyński & Wiita (1980) to mimic GR effects near the black hole. Lai & Wiseman (1997) concentrate on NS–NS binaries and the dependence of the results on the NS EOS. They add a restricted set of PN orbital terms to the dynamical equations given in Lai & Shapiro (1995) for a binary system containing two NS modeled as Riemann-S ellipsoids (cf. LRS), but they neglect relativistic corrections to the fluid motion, self-gravity and tidal interaction. Lombardi, Rasio, & Shapiro (1997) include PN corrections affecting both the orbital motion and the interior structure of the stars and explore the consequences not only for orbital stability but also for the stability of each NS against collapse. Taniguchi & Shibata 1997 and Shibata & Taniguchi 1997 provide an analytic treatment of incompressible binaries in the PN approximation. The most important result, on which these various studies all seem to agree, is that neither the relativistic effects nor the Newtonian tidal effects can be neglected if one wants to obtain a quantitatively accurate determination of the stability limits. In particular, the critical frequency corresponding to the onset of dynamical instability can be much lower than the value obtained when only one of the two effects is included. This critical frequency for the “last stable circular orbit” is potentially a measurable quantity (with LIGO/VIRGO) and can provide direct information on the NS EOS.

#### 5.1.1 Post-Newtonian Calculations of the ISCO

Lombardi, Rasio & Shapiro (1997, hereafter LRS97) have calculated PN quasi-equilibrium configurations of binary NS obeying a polytropic equation of state. Surfaces of constant density within the stars are approximated as self-similar triaxial ellipsoids, i.e., they adopt the same ellipsoidal figure of equilibrium (EFE) approximation used previously in the Newtonian study of LRS. An energy variational method is used, with the energy functional including terms both for the internal hydrodynamics of the stars and for the external orbital motion. The leading PN corrections to the internal and gravitational energies of the stars are added, and hybrid orbital terms (which are fully relativistic in the test-mass limit and always accurate to first PN order) are implemented.

The EFE treatment, while only approximate, can find an equilibrium configuration in less than a second on a typical workstation. This speed affords a quick means of generating stellar models and quasi-equilibrium sequences. The results help provide a better understanding of both GR calculations and future detections of gravitational wave signals. In addition, while many treatments of binary NS are currently limited to corotating (synchronized) sequences, the EFE approach allows straightforward construction and comparison of both corotating and (the more realistic) irrotational sequences. The irrotational sequences are found to maintain a lower maximum equilibrium mass than their corotating counterparts, although the maximum mass always increases as the orbit decays.

LRS97 use the second order variation of the energy functional to identify the innermost stable circular orbit (ISCO) along their sequences. A minimum of the energy along a sequence of equilibrium configurations with decreasing orbital separation marks the ISCO, inside of which the orbit is dynamically unstable (Fig. 4). It is often assumed that the ISCO frequency of an irrotational sequence does not differ drastically from the frequency determined from corotating calculations. The results of LRS97 help quantify this difference: the ISCO frequency along an irrotational sequence is about 17% larger than the secular ISCO frequency along the corotating sequence when the polytropic index $`n=0.5`$, and 20% larger when $`n=1`$.

Arras & Lombardi (1998) have suggested an alternative analytic approximation scheme for treating binary neutron stars. In place of an energy variational method which uses a trial density function, the 1PN orbit, Euler and continuity equations are explicitly solved. The only assumptions are that the unperturbed star is a polytrope and that the system is in quasi-equilibrium. The EFE approximation is relaxed and the problem is solved order by order in a triple expansion, with separate expansion parameters for GR, rotational, and tidal effects. This technique is the natural PN generalization of the Chandrasekhar-Milne expansion method used to treat Newtonian binaries. This method improves upon the work of LRS97 by also including PN effects for the internal fluid motion, in addition to the orbital motion. Some strong field effects can be accounted for through a hybrid scheme: energy terms which also exist for isolated non-rotating stars can be replaced with an exact expression obtained by integrating the OV equation. One is free to add any 2nd and higher order PN terms when working to 1PN order.

#### 5.1.2 Fully Relativistic Calculations of the ISCO

The first calculations in full relativity of equal mass, polytropic neutron star binaries in quasiequilibrium, synchronized orbits were performed by Baumgarte et al(1997; 1998a,b). They integrated Einstein’s equations together with the relativistic equations of hydrostatic equilibrium, obtaining numerical solutions of the exact initial-value problem and approximate quasiequilibrium evolution models for these binaries. Their numerical method for the coupled set of nonlinear elliptic equations consisted of adaptive multigrid integrations in 3D, using the DAGH software developed by the Binary Black Hole Grand Challenge Alliance to run the code in parallel (see, e.g., Parashar 1997). DAGH (“Distributive Adaptive Grid Hierarchy”) allows for convenient implementation of parallel and adaptive applications.

Baumgarte et alused the resulting models to construct sequences of constant rest-mass at different radii, locating turning points along binding energy equilibrium curves to identify the onset of orbital instability. By this means they identified the ISCO and its angular velocity. They found, in agreement with Newtonian treatments (e.g., LRS), that an ISCO exists only for polytropic indices $`n1.5`$; for softer equations of state, contact is reached prior to the onset of orbital instability.

The results of Baumgarte et alfor the ISCO are summarized in Table 1 for sequences of constant rest mass $`M_0`$ and polytropic index $`n=1`$. Also included are the values of $`J/M^2`$ for each system at the ISCO. For small rest-masses, this value is larger than unity, so that the two stars cannot form a Kerr black hole following coalescence without having to lose additional angular momentum. Note that the masses of models governed by a polytropic equation of state scale with $`K`$ as indicated in the Table. Generalizing these calculations for realistic equations of state is straightforward, but has not yet been performed.

### 5.2 Binary-Induced Collapse Instability

A surprising result coming from the numerical $`3+1`$ relativistic calculations of Wilson and collaborators (Wilson, Mathews, & Marronetti 1996; Mathews & Wilson 1997; Marronetti et al1998; hereafter WMM) is the appearance of a “binary-induced collapse instability” of the NS, with the central density of each star increasing by an amount proportional to $`1/r`$. This result, which is based on integrating an approximate subset of the Einstein field equations (assuming a conformally flat 3-metric), was surprising in light of the earlier demonstration by LRS (see, e.g., Fig 15 of LRS1) that in Newtonian gravitation, the tidal field of a companion tends to stabilize a star against radial collapse, lowering the critical value of $`\mathrm{\Gamma }`$ for collapse below $`4/3`$. Indeed, Newtonian tidal effects make the central density in a star decrease by an amount proportional to $`1/r^6`$; cf. Lai 1996). So if correct, the result of WMM thus would have to be a purely relativistic effect. In effect, the maximum stable mass of a NS in a relativistic close binary system would have to be slightly lower than that of a NS in isolation. An initially stable NS close to the maximum mass could then collapse to a black hole well before getting to the final phase of binary coalescence!

The numerical results of WMM have yet to be confirmed independently by other studies. Even if valid, the WMM effect would be of importance only if the NS EOS is very soft and the maximum stable mass for a NS in isolation is not much larger than $`1.4M_{}`$. More significant, the numerical results of WMM have been criticized by many authors on theoretical grounds. Brady & Hughes (1997) show analytically that, in the limit where the NS companion becomes a test particle of mass $`m`$, the central density of the NS remains unchanged to linear order in $`m/R`$, in contrast to what would be expected from the WMM results. LRS97 and Wiseman (1997) argue that there should be no destabilizing relativistic effect to first PN order. In contrast, WMM claim that their effect is at least partially caused by a nonlinear first PN order enhancement of the gravitational potential. But Lombardi et al(1997) also find that, to first PN order, the maximum equilibrium mass of a NS in a binary increases as the binary separation $`r`$ decreases, in agreement with the fully relativistic numerical calculations of Baumgarte et al(1997). Indeed, in a systematic radial stability analysis of their fully relativistic, corotating binary models, Baumgarte et al(1998a) conclude that the configurations are stable against collapse to black holes all the way down to the ISCO. The conclusion that binary neutron stars are stable to collapse to black holes has also been reached by means of analytic “local-asymptotic-rest-frame” calculations by Flanagan (1998) and Thorne (1997).

A direct demonstration casting doubt on the WMM effect, at least for fluid stars, is provided by the numerical simulations of Shibata, Baumgarte & Shapiro (1998). They perform a fully hydrodynamic evolution of relativistic binary stars to investigate their dynamical stability against gravitational collapse prior to merger. While in general their equations are only strictly accurate to first PN order, they retain sufficient nonlinearity to recover full GR in the limit of spherical, static stars. Shibata et alstudy both corotating and irrotational binary configurations of identical stars in circular orbits. A soft, adiabatic equation of state with $`\mathrm{\Gamma }=1.4`$ is adopted, for which the onset of instability occurs at a sufficiently small value of $`M/R`$ that the PN approximation is quite accurate. For such a soft equation of state there is no innermost stable circular orbit, so that one can study arbitrarily close binaries, while still exploring the same qualitative features exhibited by any adiabatic equation of state regarding stability against gravitational collapse. The main new result of is that, independent of the internal stellar velocity profile, the tidal field from a binary companion stabilizes a star against gravitational collapse. Specifically, one finds that neutron stars which reside on the stable branch of the mass vs central density equilibrium curve in isolation rotate about their companions for many orbital periods without undergoing collapse. Only those models which are well along on the unstable branch in isolation undergo collapse in a binary.

To demonstrate a point of principle, however, Shapiro (1998a) constructed a simple model illustrating how a highly relativistic, compact object which is stable in isolation could be driven dynamically unstable by the tidal field of a binary companion. The compact object consists of a test-particle in a relativistic orbit about a black hole while the binary companion is a distant point mass. This strong-field model suggests that first-order PN treatments of binaries, and stability analyses of binary equilibria based on orbit-averaged, mean gravitational fields, may not be adequate to rule out the instability. The main result of this simple demonstration was to provide a word of caution. On the one hand, there is mounting evidence which argues against the WMM effect. However, the possibility that sufficiently massive, highly compact NS in coalescing binaries can collapse to black holes prior to merger will not be completely ruled out until detailed hydrodynamic simulations in full GR, without approximation, are finally carried out.

### 5.3 The Final Fate of Mergers

Fully relativistic numerical simulations are clearly required to obtain quantitatively reliable coalescence waveforms. However, a numerical approach in full GR is also required for deciding between qualitatively different outcomes, even in the case of neutron stars.

Consider, for example, the simple problem of a nearly head-on collision of two identical neutron stars moving close to free-fall velocity at contact (Shapiro 1998b). Assume that each star has a mass larger than $`0.5M_{max}`$, where $`M_{max}`$ is the maximum mass of a cold neutron star. When the two stars collide, two recoil shocks propagate through each of the stars from the point of contact back along the collision axis. This shock serves to convert bulk fluid kinetic energy into thermal energy. The typical temperature is $`kTM/R`$. What happens next? There are two possibilities. One possibility is that after the merged configuration undergoes one or two large-amplitude oscillations on a dynamical timescale ($``$ms), the coalesced star, which now has a mass larger than $`M_{max}`$, collapses immediately to a black hole. Another possibility is that the thermal pressure generated by the recoil shocks is sufficient to hold up the merged star against collapse in a quasi-static, hot equilibrium state until neutrinos carry away the thermal energy on a neutrino diffusion timescale ($``$10s). The two outcomes are both plausible but very different. The implications for gravitational wave, neutrino and possibly gamma-ray bursts from NS–NS collisions are also very different for the two scenarios. Because the outcomes depend critically on the role of time-dependent, nonlinear gravitation, resolving this issue requires a numerical simulation in full GR.

Baumgarte & Shapiro (1998a) have studied the neutrino emission from the remnant of binary NS coalescence. The mass of the merged remnant is likely to exceed the stability limit of a cold, rotating neutron star. However, the angular momentum of the remnant may also approach or even exceed the Kerr limit, $`J/M^2=1`$, so that total collapse may not be possible unless some angular momentum is dissipated. Baumgarte & Shapiro (1998a) show that neutrino emission is very inefficient in decreasing the angular momentum of these merged objects and may even lead to a small increase in $`J/M^2`$. They illustrate these findings with a PN ellipsoidal model calculation. Simple arguments suggest that the remnant may undergo a bar-mode instability on a timescale similar to or shorter than the neutrino emission timescale, in which case the evolution of the remnant will be dominated by the emission of gravitational waves. But the dominant instability may be the newly discovered r-mode (Andersson 1998; Friedman & Morsink 1998), which has the potential to slow down dramatically rapidly rotating, hot neutron stars like the remnant formed by coalescence. The mechanism is the emission of current-quadrupole gravitational waves, which carry off angular momentum. The process itself may be an interesting source of detectable gravitational waves (Owen et al1998).

### 5.4 Numerical Relativity and Future Prospects

Calculations of coalescence waveforms from colliding black holes and neutron stars require the tools of numerical relativity – the art and science of solving Einstein’s equations numerically on a spacetime lattice. Numerical relativity in 3+1 dimensions is in its infancy and is fraught with many technical complications. Always present, of course, are the usual difficulties associated with solving multidimensional, nonlinear, coupled PDE’s. But these difficulties are not unique to relativity; they are also present in hydrodynamics, for example. But numerical relativity must also deal with special problems, like the appearance of singularities in a numerical simulation. Singularities are regions where physical quantities like the curvature (i.e., tidal field) or the matter density blow up to infinity. Singularities are always present inside black holes. Encountering such a singularity causes a numerical simulation to crash, even if the singularity is inside a black hole event horizon and causally disconnected from the outside world. Another special difficulty that confronts numerical relativity is the challenge of determining the asymptotic gravitational waveform which is generated during a strong-field interaction. The asymptotic waveform is just a small perturbation to the background metric and it must be determined in the wave zone far from the strong-field sources. Such a determination presents a problem of dynamic range: one wants to measure the waveform accurately far from the sources, but one must put most of the computational resources (i.e. grid) in the vicinity of those same sources, where most of the nonlinear dynamics occurs, Moreover, to determine the outgoing asymptotic emission, one must wait for the wave train to propagate out into the far zone, but by then, the simulation may be losing accuracy because of the growth of singularities in the strong-field, near zone.

Arguably the most outstanding problem in numerical relativity is the coalescence of binary black holes. The late stages of the merger can only be solved by numerical means. To advance this effort, the National Science Foundation recently funded a “Grand Challenge Alliance” of numerical relativists and computer scientists at various institutions in the United States. At present, no code can integrate two black holes in binary orbit for as long as a few periods, let alone long enough to get a gravitational wave out to, say, 10 per cent accuracy. That is because the multiple complications described above all conspired to make the integration of two black holes increasingly divergent at late times, well before the radiation content could be reliably determined. Most recently, however, the Grand Challenge Alliance has reported several promising developments (for updates, see their web site at http://www.npac.syr.edu/projects/bh). New formulations of Einstein’s field equations have been proposed (Choquet-Bruhat & York 1995; Bona et al1995; van Putten & Eardley (1996); Friedrich 1996; Anderson, Choquet-Bruhat & York 1998) that cast them is a flux-conservative, first order, hyberbolic form where the only nonzero characteristic speed is that of light. As a result of this new formulation, it may be possible to “cut-out” the interior regions of the black holes from the numerical grid and install boundary conditions at the hole horizons (“horizon boundary conditions”). Removing the black hole interiors is crucial since that is where the spacetime singularities reside, and they are the main sources of the computational inaccuracies. So now there is renewed confidence that the binary black hole problem can be solved.

The binary neutron star coalescence problem is both easier and more difficult than the binary black hole problem. It is easier in that there are no singularities and no horizons to contend with numerically. It is more difficult in that one cannot work with the vacuum Einstein equations, but must solve the the equations of relativistic hydrodynamics in conjunction with the field equations. The 3+1 ADM equations may prove adequate to solve the binary neutron star problem. This would be convenient since some of the new hyperbolic formulations require taking derivatives of the original ADM equations, and these may introduce inaccuracies if matter sources are present. A modified set of ADM equations has recently been proposed by Shibata & Nakamura (1995; see also Baumgarte & Shapiro 1998b) which casts the system into a more appealing mathematical form and which exhibits improved stability in tests of gravitational wave propagation. This modified set may prove to be an effective compromise for dealing with the binary neutron star problem.

As discussed previously, there are several independent efforts underway to tackle NS binary coalescence in full GR, including a NASA-sponsored Grand Challenge project (for updates, see the web sites at http://jean-luc.ncsa.uiuc.edu/nsngc and http://wugrav.wustl.edu/Relativ/nsgc.html). It is conceivable that the binary NS problem will be solved before the binary BH problem, at least for the evolutionary phase prior to merger and shock heating. However, any progress in solving either one of these problems will likely serve to advance the other effort as well, given the overlap of numerical algorithms and software.

## 6 Nonsynchronized binaries

### 6.1 Irrotational Equilibrium Sequences

It is very likely that the synchronization time in close NS binaries always remains longer than the orbital decay time due to gravitational radiation (Kochanek 1992; Bildsten & Cutler 1992). In particular, Bildsten & Cutler (1992) show with simple dimensional arguments that one would need an implausibly small value of the effective viscous time, approaching $`t_{visc}R/c`$, in order to reach complete synchronization just before final merging.

In the opposite limiting regime where viscosity is completely negligible, the fluid circulation in the binary system is conserved during the orbital decay and the stars behave approximately as Darwin-Riemann ellipsoids (Kochanek 1992; LRS3). Of particular importance are the irrotational Darwin-Riemann configurations, obtained when two initially nonspinning (or, in reality, slowly spinning) NS evolve in the absence of significant viscosity. Compared to synchronized systems, these irrotational configurations exhibit smaller deviations from point-mass Keplerian behavior at small $`r`$. However, as shown in LRS3, irrotational configurations for binary NS with $`\mathrm{\Gamma }>2`$ can still become dynamically unstable near contact. Thus the final coalescence of two NS in a nonsynchronized binary system can still be driven entirely by hydrodynamic instabilities.

Sequences of Newtonian equilibrium configurations for irrotational binaries were computed by LRS (see especially LRS3) using an enegy variational method and modeling the stars explicitly as compressible Darwin-Riemann ellipsoids. LRS showed that a dynamical instability can occur in all close binary configurations, whether synchronized or not, provided that the system contains sufficiently incompressible stars. For binary systems containing two nonspinning NS with a stiff EOS, the hydrodynamic instability can significantly accelerate the coalescence at small separation, with the radial infall velocity just prior to contact reaching typically about 10% of the tangential orbital velocity.

Using a self-consistent field method, Uryu & Eriguchi (1998) have calculated the first exact 3D equilibrium solutions for irrotational equal-mass binaries with polytropic components in Newtonian gravity. They find that a dynamical instability is reached before contact when the polytropic index $`n<0.7`$, i.e., when $`\mathrm{\Gamma }>2.4`$, in reasonable agreement with the approximate results of LRS. When PN effects are taken into account, however, it is found that dynamical instability sets in before contact for even softer EOS (see LRS97).

Fully relativistic generalizations of the calculations by Uryu & Eriguchi (1998) are currently being performed by several groups. Bonazzola, Gourgoulhon, & Marck (1998) report the first relativistic results from calculations of irrotational equilibrium sequences with constant baryon number. They solve the Einstein field equations numerically in the Wilson-Mathews approximation (cf. §5.2). The velocity field inside the stars is computed by solving an elliptical equation for the velocity scalar potential. Their most significant result is that, although the central NS density decreases much less with the binary separation than in the corotating case, it still decreases. Thus, no tendency is found for the stars to individually collapse to black holes prior to final merging.

### 6.2 Coalescence of Nonsynchronized Binaries

For nonsynchronized binaries, the final hydrodynamic coalescence of the two stars can be very complicated (Fig. 5), leading to significant differences in the gravitational wave emission (Fig. 6) compared to the synchronized case, and an additional dependence of the gravitational radiation waveforms on the stellar spins (not included in eqs. 3–5).

Consider for example the case of an irrotational system (containing two initially nonspinning stars). Because the two stars appear to be counter-spinning in the corotating frame of the binary, a vortex sheet (where the tangential velocity jumps discontinuously by $`\mathrm{\Delta }v=|v_+v_{}|\mathrm{\Omega }r`$) appears when the stellar surfaces come into contact. Such a vortex sheet is Kelvin-Helmholtz unstable on all wavelengths and the hydrodynamics is therefore extremely difficult to model accurately given the limited spatial resolution of 3D calculations, even in the Newtonian limit. The breaking of the vortex sheet generates a large turbulent viscosity so that the final configuration may no longer be irrotational. In numerical simulations, however, vorticity is generated mostly through spurious shear viscosity introduced by the spatial discretization (see, e.g., Lombardi et al1998 for a detailed study of spurious viscosity in SPH simulations). The late-time decay of the gravitational waves seen in Fig. 6 may be dominated by this spurious viscosity.

An additional difficulty is that nonsynchronized configurations evolving rapidly by gravitational radiation emission tend to develop small but significant tidal lags, with the long axes of the two components becoming misaligned (LRS5). This is a purely dynamical effect, present even if the viscosity is zero, but its magnitude depends on the entire previous evolution of the system. Thus the construction of initial conditions for hydrodynamic calculations of nonsynchronized binary coalescence must incorporate the gravitational radiation reaction self-consistently. Instead, previous hydrodynamic calculations of nonsynchronized binary coalescence (Shibata et al1992; Davies et al1994; Zhuge et al1994, 1996; Ruffert et al1997) used very crude initial conditions consisting of two spherical stars placed on an inspiral trajectory calculated for two point masses. The SPH calculation illustrated in Figs. 5 and 6 (performed by the authors) used the ellipsoidal approximation of LRS to construct a more realistic (but still not exact) initial condition for an irrotational system at the onset of dynamical instability. Fully relativistic, self-consistent calculations for the coalescence of nonsynchronized NS binaries have yet to be attempted.

It is a pleasure to thank Thomas Baumgarte and Dong Lai for several useful discussions. F.A.R. has been supported in part by NSF Grant AST-9618116 and by a Sloan Research Fellowship. S.L.S. has been supported in part by NSF Grant AST 96-18524 and NSF Binary Black Hole Grand Challenge Grant NSF PHY/ASC 93-18152/ASC (ARPA supplemented), and by NASA Grant NAG5-7152. This work was supported by the National Computational Science Alliance under Grants AST970022N (F.A.R.), and AST 970023N and PHY 970014N (S.L.S.), and utilized the NCSA SGI/Cray POWER CHALLENGE array and the NCSA SGI/Cray Origin2000. F.A.R. also thanks the Aspen Center for Physics, and the Theoretical Astrophysics Division of the Harvard-Smithsonian Center for Astrophysics for hospitality.     

## References