We present SphericalNR, a new framework for the publicly available Einstein Toolkit that numerically solves the Einstein field equations coupled to the equations of general relativistic magnetohydrodynamics (GRMHD) in a 3+1 split of spacetime in spherical coordinates without symmetry assumptions. The spacetime evolution is performed using reference-metric versions of either the Baumgarte-Shapiro-Shibata-Nakamura equations or the fully covariant and conformal Z4 system with constraint damping. We have developed a reference-metric version of the Valencia formulation of GRMHD with a vector potential method, guaranteeing the absence of magnetic monopoles during the evolution. In our framework, every dynamical field (both spacetime and matter) is evolved using its components in an orthonormal basis with respect to the spherical reference-metric. Furthermore, all geometric information about the spherical coordinate system is encoded in source terms appearing in the evolution equations. This allows for the straightforward extension of Cartesian high-resolution shock-capturing finite volume codes to use spherical coordinates with our framework. To this end, we have adapted GRHydro, a Cartesian finite volume GRMHD code already available in the Einstein Toolkit, to use spherical coordinates. We present the full evolution equations of the framework, as well as details of its implementation in the Einstein Toolkit. We validate SphericalNR by demonstrating it passes a variety of challenging code tests in static and dynamical spacetimes.

We study effects of heating by dark matter (DM) annihilation on black hole gas accretion. We observe that, for reasonable assumptions about DM densities in spikes around supermassive black holes, as well as DM masses and annihilation cross-sections within the standard WIMP model, heating by DM annihilation may have an appreciable effect on the accretion onto Sgr A$^*$ in the Galactic center. Motivated by this observation we study the effects of such heating on Bondi accretion, i.e. spherically symmetric, steady-state Newtonian accretion onto a black hole. We consider different adiabatic indices for the gas, and different power-law exponents for the DM density profile. We find that typical transonic solutions with heating have a significantly reduced accretion rate. However, for many plausible parameters, transonic solutions do not exist, suggesting a breakdown of the underlying assumptions of steady-state Bondi accretion. Our findings indicate that heating by DM annihilation may play an important role in the accretion onto supermassive black holes at the center of galaxies, and may help explain the low accretion rate observed for Sgr A$^*$.

The "direct collapse" scenario has emerged as a promising evolutionary track for the formation of supermassive black holes early in the Universe. In an idealized version of such a scenario, a uniformly rotating supermassive star spinning at the mass-shedding (Keplerian) limit collapses gravitationally after it reaches a critical configuration. Under the assumption that the gas is dominated by radiation pressure, this critical configuration is characterized by unique values of the dimensionless parameters $J/M^2$ and $R_p/M$, where $J$ is the angular momentum, $R_p$ the polar radius, and $M$ the mass. Motivated by a previous perturbative treatment we adopt a fully nonlinear approach to evaluate the effects of gas pressure on these dimensionless parameters for a large range of masses. We find that gas pressure has a significant effect on the critical configuration even for stellar masses as large as $M \simeq 10^6 M_{\odot}$. We also calibrate two approximate treatments of the gas pressure perturbation in a comparison with the exact treatment, and find that one commonly used approximation in particular results in increasing deviations from the exact treatment as the mass decreases, and the effects of gas pressure increase. The other approximation, however, proves to be quite robust for all masses $M \gtrsim 10^4 M_{\odot}$.

We investigate the dynamical stability of relativistic, differentially rotating, quasi-toroidal models of neutron stars through hydrodynamical simulations in full general relativity. We find that all quasi-toroidal configurations studied in this work are dynamically unstable against the growth of non-axisymmetric modes. Both one-arm and bar mode instabilities grow during their evolution. We find that very high rest mass configurations collapse to form black holes. Our calculations suggest that configurations whose rest mass is less than the binary neutron star threshold mass for prompt collapse to black hole transition dynamically to spheroidal, differentially rotating stars that are dynamically stable, but secularly unstable. Our study shows that the existence of extreme quasi-toroidal neutron star equilibrium solutions does not imply that long-lived binary neutron star merger remnants can be much more massive than previously found. Finally, we find models that are initially supra-Kerr ($J/M^2>1$) and undergo catastrophic collapse on a dynamical timescale, in contrast to what was found in earlier works. However, cosmic censorship is respected in all of our cases. Our work explicitly demonstrates that exceeding the Kerr bound in rotating neutron star models does not imply dynamical stability.

The discovery of quasars at increasingly large cosmological redshifts may favor "direct collapse" as the most promising evolutionary route to the formation of supermassive black holes. In this scenario, supermassive black holes form when their progenitors - supermassive stars - become unstable to gravitational collapse. For uniformly rotating stars supported by pure radiation pressure and spinning at the mass-shedding limit, the critical configuration at the onset of collapse is characterized by universal values of the dimensionless spin and radius parameters $J/M^2$ and $R/M$, independent of mass $M$. We consider perturbative effects of gas pressure, magnetic fields, dark matter and dark energy on these parameters, and thereby determine the domain of validity of this universality. We obtain leading-order corrections for the critical parameters and establish their scaling with the relevant physical parameters. We compare two different approaches to approximate the effects of gas pressure, which plays the most important role, find identical results for the above dimensionless parameters, and also find good agreement with recent numerical results.

Numerical relativity codes that do not make assumptions on spatial symmetries most commonly adopt Cartesian coordinates. While these coordinates have many attractive features, spherical coordinates are much better suited to take advantage of approximate symmetries in a number of astrophysical objects, including single stars, black holes and accretion disks. While the appearance of coordinate singularities often spoils numerical relativity simulations in spherical coordinates, especially in the absence of any symmetry assumptions, it has recently been demonstrated that these problems can be avoided if the coordinate singularities are handled analytically. This is possible with the help of a reference-metric version of the Baumgarte-Shapiro-Shibata-Nakamura formulation together with a proper rescaling of tensorial quantities. In this paper we report on an implementation of this formalism in the Einstein Toolkit. We adapt the Einstein Toolkit infrastructure, originally designed for Cartesian coordinates, to handle spherical coordinates, by providing appropriate boundary conditions at both inner and outer boundaries. We perform numerical simulations for a disturbed Kerr black hole, extract the gravitational wave signal, and demonstrate that the noise in these signals is orders of magnitude smaller when computed on spherical grids rather than Cartesian grids. With the public release of our new Einstein Toolkit thorns, our methods for numerical relativity in spherical coordinates will become available to the entire numerical relativity community.

We report on a new open-source, user-friendly numerical relativity code package called SENR/NRPy+. Our code extends previous implementations of the BSSN reference-metric formulation to a much broader class of curvilinear coordinate systems, making it ideally-suited to modeling physical configurations with approximate or exact symmetries. In the context of modeling black hole dynamics, it is orders of magnitude more efficient than other widely used open-source numerical relativity codes. NRPy+ provides a Python-based interface in which equations are written in natural tensorial form and output at arbitrary finite difference order as highly efficient C code, putting complex tensorial equations at the scientist's fingertips without the need for an expensive software license. SENR provides the algorithmic framework that combines the C codes generated by NRPy+ into a functioning numerical relativity code. We validate against two other established, state-of-the-art codes, and achieve excellent agreement. For the first time we demonstrate--in the context of puncture, trumpet, and dual black hole evolutions--nearly exponential convergence of constraint violation and gravitational waveform errors to zero as the order of spatial finite difference derivatives is increased, while holding the spherical-like coordinate grids fixed at moderate resolution. Such behavior outside the horizons is remarkable, as numerical errors do not converge to zero inside horizons, and all points along the polar axis are coordinate singularities. The formulation addresses such coordinate singularities via cell-centered grids and a simple change of basis that analytically regularizes tensor components with respect to the coordinates. Future plans include extending this formulation to allow dynamical coordinate grids and bispherical-like distribution of points to efficiently capture orbiting compact binary dynamics.

In this paper we report on what we believe is the first successful implementation of relativistic hydrodynamics, coupled to dynamical spacetimes, in spherical polar coordinates without symmetry assumptions. We employ a high-resolution shock-capturing scheme, which requires that the equations be cast in flux-conservative form. One example of such a form is the :Valencia" formulation, which has been adopted in numerous applications, in particular in Cartesian coordinates. Here we generalize this formulation to allow for a reference-metric approach, which provides a natural framework for calculations in curvilinear coordinates. In spherical polar coordinates, for example, it allows for an analytical treatment of the singular r and sin(\theta) terms that appear in the equations. We experiment with different versions of our generalized Valencia formulation in numerical implementations of relativistic hydrodynamics for both fixed and dynamical spacetimes. We consider a number of different tests -- non-rotating and rotating relativistic stars, as well as gravitational collapse to a black hole -- to demonstrate that our formulation provides a promising approach to performing fully relativistic astrophysics simulations in spherical polar coordinates.

We perform hydrodynamical simulations of neutron-star mergers for a large sample of temperature-dependent, nuclear equations of state, and determine the threshold mass above which the merger remnant promptly collapses to form a black hole. We find that, depending on the equation of state, the threshold mass is larger than the maximum mass of a non-rotating star in isolation by between 30 and 70 per cent. Our simulations also show that the ratio between the threshold mass and maximum mass is tightly correlated with the compactness of the non-rotating maximum-mass configuration. We speculate on how this relation can be used to derive constraints on neutron-star properties from future observations.

In the absence of symmetry assumptions most numerical relativity simulations adopt Cartesian coordinates. While Cartesian coordinates have some desirable properties, spherical polar coordinates appear better suited for certain applications, including gravitational collapse and supernova simulations. Development of numerical relativity codes in spherical polar coordinates has been hampered by the need to handle the coordinate singularities at the origin and on the axis, for example by careful regularization of the appropriate variables. Assuming spherical symmetry and adopting a covariant version of the BSSN equations, Montero and Cordero-Carri\'on recently demonstrated that such a regularization is not necessary when a partially implicit Runge-Kutta (PIRK) method is used for the time evolution of the gravitational fields. Here we report on an implementation of the BSSN equations in spherical polar coordinates without any symmetry assumptions. Using a PIRK method we obtain stable simulations in three spatial dimensions without the need to regularize the origin or the axis. We perform and discuss a number of tests to assess the stability, accuracy and convergence of the code, namely weak gravitational waves, "hydro-without-hydro" evolutions of spherical and rotating relativistic stars in equilibrium, and single black holes.

We apply a von Zeipel gravity darkening model to corotating binaries to obtain a simple, analytical expression for the emergent radiative flux from a tidally distorted primary orbiting a point-mass secondary. We adopt a simple Roche model to determine the envelope structure of the primary, assumed massive and centrally condensed, and use the results to calculate the flux. As for single rotating stars, gravity darkening reduces the flux along the stellar equator of the primary, but, unlike for rotating stars, we find that gravity brightening enhances the flux in a region around the stellar poles. We identify a critical limiting separation beyond which hydrostatic equilibrium no longer is possible, whereby the flux vanishes at the point on the stellar equator of the primary facing the companion. For equal-mass binaries, the total luminosity is reduced by about 13 % when this limiting separation is reached.

White dwarf-neutron star binaries generate detectable gravitational radiation. We construct Newtonian equilibrium models of corotational white dwarf-neutron star (WDNS) binaries in circular orbit and find that these models terminate at the Roche limit. At this point the binary will undergo either stable mass transfer (SMT) and evolve on a secular time scale, or unstable mass transfer (UMT), which results in the tidal disruption of the WD. The path a given binary will follow depends primarily on its mass ratio. We analyze the fate of known WDNS binaries and use population synthesis results to estimate the number of LISA-resolved galactic binaries that will undergo either SMT or UMT. We model the quasistationary SMT epoch by solving a set of simple ordinary differential equations and compute the corresponding gravitational waveforms. Finally, we discuss in general terms the possible fate of binaries that undergo UMT and construct approximate Newtonian equilibrium configurations of merged WDNS remnants. We use these configurations to assess plausible outcomes of our future, fully relativistic simulations of these systems. If sufficient WD debris lands on the NS, the remnant may collapse, whereby the gravitational waves from the inspiral, merger, and collapse phases will sweep from LISA through LIGO frequency bands. If the debris forms a disk about the NS, it may fragment and form planets.

Most numerical models of binary stars - in particular neutron stars in compact binaries - assume the companions to be either corotational or irrotational. Either one of these assumptions leads to a significant simplification in the hydrodynamic equations of stationary equilibrium. In this paper we develop a new formalism for the construction of binary stars with circulation intermediate between corotational and irrotational. Generalizing the equations for irrotational flow we cast the Euler equation, which is an algebraic equation in the case of corotational or irrotational fluid flow, as an elliptic equation for a new auxiliary quantity. We also suggest a parameterized decomposition of the fluid flow that allows for a variation of the stellar circulation.

Black hole-neutron star (BHNS) binary mergers are candidate engines for generating both short-hard gamma-ray bursts (SGRBs) and detectable gravitational waves. Using our most recent conformal thin-sandwich BHNS initial data and our fully general relativistic hydrodynamics code, which is now AMR-capable, we are able to efficiently and accurately simulate these binaries from large separations through inspiral, merger, and ringdown. We evolve the metric using the BSSN formulation with the standard moving puncture gauge conditions and handle the hydrodynamics with a high-resolution shock-capturing scheme. We explore the effects of BH spin (aligned and anti-aligned with the orbital angular momentum) by evolving three sets of initial data with BH:NS mass ratio q=3: the data sets are nearly identical, except the BH spin is varied between a/M = -0.5 (anti-aligned), 0.0, and 0.75. The number of orbits before merger increases with a/M, as expected. We also study the nonspinning BH case in more detail, varying q between 1, 3, and 5. We calculate gravitational waveforms for the cases we simulate and compare them to binary black-hole waveforms. Only a small disk (< 0.01 M_sun) forms for the anti-aligned spin case (a/M = -0.5) and for the most extreme mass ratio case (q=5). By contrast, a massive (M_disk is about 0.2 M_sun), hot disk forms in the rapidly spinning (a/M = 0.75) aligned BH case. Such a disk could drive a SGRB,possibly by, e.g., producing a copious flux of neutrino-antineutino pairs.

Black hole-neutron star (BHNS) binaries are expected to be among the leading sources of gravitational waves observable by ground-based detectors, and may be the progenitors of short-hard gamma ray bursts (SGRBs) as well. Here, we discuss our new fully general relativistic calculations of merging BHNS binaries, which use high-accuracy, low-eccentricity, conformal thin-sandwich configurations as initial data. Our evolutions are performed using the moving puncture method and include a fully relativistic, high-resolution shock-capturing hydrodynamics treatment. Focusing on systems in which the neutron star is irrotational and the black hole is nonspinning with a 3:1 mass ratio, we investigate the inspiral, merger, and disk formation in the system. We find that the vast majority of material is promptly accreted and no more than 3% of the neutron star's rest mass is ejected into a tenuous, gravitationally bound disk. We find similar results for mass ratios of 2:1 and 1:1, even when we reduce the NS compaction in the 2:1 mass ratio case. These ambient disks reach temperatures suitable for triggering SGRBs, but their masses may be too small to produce the required total energy output. We measure gravitational waveforms and compute the effective strain in frequency space, finding measurable differences between our waveforms and those produced by binary black hole mergers within the advanced LIGO band. These differences appear at frequencies corresponding to the emission that occurs when the NS is tidally disrupted and accreted by the black hole. The resulting information about the radius of the neutron star may be used to constrain the neutron star equation of state.

We construct new models of black hole-neutron star binaries in quasiequilibrium circular orbits by solving Einstein's constraint equations in the conformal thin-sandwich decomposition together with the relativistic equations of hydrostationary equilibrium. We adopt maximal slicing, assume spatial conformal flatness, and impose equilibrium boundary conditions on an excision surface (i.e., the apparent horizon) to model the black hole. In our previous treatment we adopted a "leading-order" approximation for a parameter related to the black-hole spin in these boundary conditions to construct approximately nonspinning black holes. Here we improve on the models by computing the black hole's quasilocal spin angular momentum and setting it to zero. As before, we adopt a polytropic equation of state with adiabatic index Gamma=2 and assume the neutron star to be irrotational. In addition to recomputing several sequences for comparison with our earlier results, we study a wider range of neutron star masses and binary mass ratios. To locate the innermost stable circular orbit we search for turning points along both the binding energy and total angular momentum curves for these sequences. Unlike for our previous approximate boundary condition, these two minima now coincide. We also identify the formation of cusps on the neutron star surface, indicating the onset of tidal disruption. Comparing these two critical binary separations for different mass ratios and neutron star compactions we distinguish those regions that will lead to a tidal disruption of the neutron star from those that will result in the plunge into the black hole of a neutron star more or less intact, albeit distorted by tidal forces.

We construct quasiequilibrium sequences of black hole-neutron star binaries in general relativity. We solve Einstein's constraint equations in the conformal thin-sandwich formalism, subject to black hole boundary conditions imposed on the surface of an excised sphere, together with the relativistic equations of hydrostatic equilibrium. In contrast to our previous calculations we adopt a flat spatial background geometry and do not assume extreme mass ratios. We adopt a Gamma=2 polytropic equation of state and focus on irrotational neutron star configurations as well as approximately nonspinning black holes. We present numerical results for ratios of the black hole's irreducible mass to the neutron star's ADM mass in isolation of M_{irr}^{BH}/M_{ADM,0}^{NS} = 1, 2, 3, 5, and 10. We consider neutron stars of baryon rest mass M_B^{NS}/M_B^{max} = 83% and 56%, where M_B^{max} is the maximum allowed rest mass of a spherical star in isolation for our equation of state. For these sequences, we locate the onset of tidal disruption and, in cases with sufficiently large mass ratios and neutron star compactions, the innermost stable circular orbit. We compare with previous results for black hole-neutron star binaries and find excellent agreement with third-order post-Newtonian results, especially for large binary separations. We also use our results to estimate the energy spectrum of the outgoing gravitational radiation emitted during the inspiral phase for these binaries.

We present new sequences of general relativistic, quasiequilibrium black hole-neutron star binaries. We solve for the gravitational field in the conformal thin-sandwich decomposition of Einstein's field equations, coupled to the equations of relativistic hydrostatic equilibrium for a perfect fluid. We account for the black hole by solving these equations in the background metric of a Schwarzschild black hole whose mass is much larger than that of the neutron star. The background metric is treated in Kerr-Schild as well as isotropic coordinates. For the neutron star, we assume a polytropic equation of state with adiabatic index Gamma=2, and solve for both irrotational and corotational configurations. By comparing the results of irrotational and synchronized configurations with the same background metric, we conclude that the effect of the rotation on the location of tidal break-up is only on the order of a few percent. The different choices in the background also lead to differences of order a few percent, which may be an indication of the level to which these configurations approximate quasiequilibrium.

Many of the recent numerical simulations of binary black holes in vacuum adopt the moving puncture approach. This successful approach avoids the need to impose numerical excision of the black hole interior and is easy to implement. Here we wish to explore how well the same approach can be applied to moving black hole punctures in the presence of relativistic hydrodynamic matter. First, we evolve single black hole punctures in vacuum to calibrate our BSSN (Baumgarte-Shapiro-Shibata-Nakamura) implementation and to confirm that the numerical solution for the exterior spacetime is invariant to any ``junk'' (i.e., constraint-violating) initial data employed in the black hole interior. Then we focus on relativistic Bondi accretion onto a moving puncture Schwarzschild black hole as a numerical testbed for our high-resolution shock-capturing relativistic hydrodynamics scheme. We find that the hydrodynamical equations can be evolved successfully in the interior without imposing numerical excision. These results help motivate the adoption of the moving puncture approach to treat the binary black hole-neutron star problem using conformal thin-sandwich initial data.

We follow the inspiral and merger of equal-mass black holes (BHs) by the moving puncture technique and demonstrate that both the exterior solution and the asymptotic gravitational waveforms are unchanged when the initial interior solution is replaced by constraint-violating ``junk'' initial data. We apply this result to evolve conformal thin-sandwich (CTS) binary BH initial data by filling their excised interiors with arbitrary, but smooth, initial data and evolving with standard puncture gauge choices. The waveforms generated for both puncture and filled-CTS initial data are remarkably similar, and there are only minor differences between irrotational and corotational CTS BH binaries. Even the interior solutions appear to evolve to the same constraint-satisfying solution at late times, independent of the initial data.

We construct quasiequilibrium sequences of black hole-neutron star binaries for arbitrary mass ratios by solving the constraint equations of general relativity in the conformal thin-sandwich decomposition. We model the neutron star as a stationary polytrope satisfying the relativistic equations of hydrodynamics, and account for the black hole by imposing equilibrium boundary conditions on the surface of an excised sphere (the apparent horizon). In this paper we focus on irrotational configurations, meaning that both the neutron star and the black hole are approximately nonspinning in an inertial frame. We present results for a binary with polytropic index n=1, mass ratio M_{irr}^{BH}/M_{B}^{NS}=5 and neutron star compaction M_{ADM,0}^{NS}/R_0=0.0879, where M_{irr}^{BH} is the irreducible mass of the black hole, M_{B}^{NS} the neutron star baryon rest-mass, and M_{ADM,0}^{NS} and R_0 the neutron star Arnowitt-Deser-Misner mass and areal radius in isolation, respectively. Our models represent valid solutions to Einstein's constraint equations and may therefore be employed as initial data for dynamical simulations of black hole-neutron star binaries.

We have calculated the first dynamical evolutions of merging black hole-neutron star binaries that treat the combined spacetime in a nonperturbative general relativistic framework. Using the conformal flatness approximation, we have studied how the location of the tidal disruption radius with respect to the the black hole horizon and innermost stable circular orbit (ISCO) affects the qualitative evolution of the system. Based on simple arguments, we show that for a binary mass ratio q>~0.24, tidal disruption occurs outside the ISCO, while the opposite is true for q<~0.24. When tidal disruption occurs sufficiently far outside the ISCO, mass is transferred unstably from the neutron star to the black hole, resulting in the complete disruption of the neutron star. When tidal disruption occurs slightly within the ISCO, we find that some of the mass forms an extremely hot disk around the black hole. The resulting configurations in this case are excellent candidates for the progenitors of short-hard gamma ray bursts.

This paper provides a brief review of the history of our understanding and knowledge of black holes. Starting with early speculations on ``dark stars'' I discuss the Schwarzschild "black hole" solution to Einstein's field equations and the development of its interpretation from "physically meaningless" to describing the perhaps most exotic and yet "most perfect" macroscopic object in the universe. I describe different astrophysical black hole populations and discuss some of their observational evidence. Finally I close by speculating about future observations of black holes with the new generation of gravitational wave detectors.

The recent localization of some short-hard gamma ray bursts (GRBs) in galaxies with low star formation rates has lent support to the suggestion that these events result from compact object binary mergers. We discuss how new simulations in general relativity are helping to identify the central engine of short-hard GRBs. Motivated by our latest relativistic black hole-neutron star merger calculations, we discuss a scenario in which these events may trigger short-hard GRBs, and compare this model to competing relativistic models involving binary neutron star mergers and the delayed collapse of hypermassive neutron stars. Distinguishing features of these models may help guide future GRB and gravitational wave observations to identify the nature of the sources.

We calculate the first dynamical evolutions of merging black hole-neutron star binaries that construct the combined black hole-neutron star spacetime in a general relativistic framework. We treat the metric in the conformal flatness approximation, and assume that the black hole mass is sufficiently large compared to that of the neutron star so that the black hole remains fixed in space. Using a spheroidal spectral methods solver, we solve the resulting field equations for a neutron star orbiting a Schwarzschild black hole. The matter is evolved using a relativistic, Lagrangian, smoothed particle hydrodynamics (SPH) treatment. We take as our initial data recent quasiequilibrium models for synchronized neutron star polytropes generated as solutions of the conformal thin-sandwich (CTS) decomposition of the Einstein field equations. We are able to construct from these models relaxed SPH configurations whose profiles show good agreement with CTS solutions. Our adiabatic evolution calculations for neutron stars with low compactness show that mass transfer, when it begins while the neutron star orbit is still outside the innermost stable circular orbit, is more unstable than is typically predicted by analytical formalisms. This dynamical mass loss is found to be the driving force in determining the subsequent evolution of the binary orbit and the neutron star, which typically disrupts completely within a few orbital periods. The majority of the mass transferred onto the black hole is accreted promptly; a significant fraction (~30%) of the mass is shed outward as well, some of which will become gravitationally unbound and ejected completely from the system. The remaining portion forms an accretion disk around the black hole, and could provide the energy source for short-duration gamma ray bursts.

We construct numerical models of the newly discovered binary pulsar J0737-3039A, both with a fully relativistic, uniformly rotating, equilibrium code that handles arbitrary spins and in the relativistic, slow-rotation approximation. We compare results for a representative sample of viable nuclear equations of state (EOS) that span three, qualitatively different, classes of models for the description of nuclear matter. A future dynamical measurement of the neutron star's moment of inertia from pulsar timing data will impose significant constraints on the nuclear EOS. Even a moderately accurate measurement (<~ 10 %) may be able to rule out some of these competing classes. Using the measured mass, spin and moment of inertia to identify the optimal model computed from different EOSs, one can determine the pulsar's radius.

We present a formalism for constructing quasi-equilibrium binary black hole initial data suitable for numerical evolution. We construct quasi-equilibrium models by imposing an approximate helical Killing symmetry appropriate for quasi-circular orbits. We use the sum of two Kerr-Schild metrics as our background metric, thereby improving on conformally flat backgrounds that do not accommodate rotating black holes and providing a horizon-penetrating lapse, convenient for implementing black hole excision. We set inner boundary conditions at an excision radius well inside the apparent horizon and construct these boundary conditions to incorporate the quasi-equilibrium condition and recover the solution for isolated black holes in the limit of large separation. We use our formalism both to generate initial data for binary black hole evolutions and to construct a crude quasi-equilibrium, inspiral sequence for binary black holes of fixed irreducible mass.

We present a new numerical method for the construction of quasiequilibrium models of black hole-neutron star binaries. We solve the constraint equations of general relativity, decomposed in the conformal thin-sandwich formalism, together with the Euler equation for the neutron star matter. We take the system to be stationary in a corotating frame and thereby assume the presence of a helical Killing vector. We solve these coupled equations in the background metric of a Kerr-Schild black hole, which accounts for the neutron star's black hole companion. In this paper we adopt a polytropic equation of state for the neutron star matter and assume large black hole--to--neutron star mass ratios. These simplifications allow us to focus on the construction of quasiequilibrium neutron star models in the presence of strong-field, black hole companions. We summarize the results of several code tests, compare with Newtonian models, and locate the onset of tidal disruption in a fully relativistic framework.

We determine the innermost stable circular orbit (ISCO) of binary neutron stars (BNSs) by performing dynamical simulations in full general relativity. Evolving quasiequilibrium (QE) binaries that begin at different separations, we bracket the location of the ISCO by distinguishing stable circular orbits from unstable plunges. We study Gamma=2 polytropes of varying compactions in both corotational and irrotational equal-mass binaries. For corotatonal binaries we find an ISCO orbital angular frequency somewhat smaller than that determined by applying turning-point methods to QE initial data. For the irrotational binaries the initial data sequences terminate before reaching a turning point, but we find that the ISCO frequency is reached prior to the termination point. Our findings suggest that the ISCO frequency varies with compaction but does not depend strongly on the stellar spin. Since the observed gravitational wave signal undergoes a transition from a nearly periodic ``chirp'' to a burst at roughly twice the ISCO frequency, the measurement of its value by laser interferometers (e.g LIGO) will be important for determining some of the physical properites of the underlying stars

The merger of binary neutron stars is likely to lead to differentially rotating remnants. In this paper, we survey several cold nuclear equations of state (EOSs) and numerically construct models of differentially rotating neutron stars in general relativity. For each EOS we tabulate maximum allowed masses as a function of the degree of differential rotation. We also determine effective polytropic indices, and compare the maximum allowed masses with those for the corresponding polytropes. We consistently find larger mass increases for the polytropes, but even for the nuclear EOSs we typically find maximum masses 50% higher than the corresponding values for nonrotating (TOV) stars. We evaluate our findings for the six observed binary neutron star (pulsar) systems, including the recently discovered binary pulsar J0737-3039. For each EOS we determine whether their merger will automatically lead to prompt collapse to a black hole, or whether the remnant can be supported against collapse by uniform rotation (possibly as a supramassive star) or differential rotation (possibly as a hypermassive star). For hypermassive stars, delayed collapse to a black hole is likely. For the most recent EOSs we survey the merger remnants can all be supported by rotation against prompt collapse, but their actual fate will depend on the nonequilibrium dynamics of the coalescence event. Gravitational wave observations of coalescing binary neutron stars may be able to distinguish these outcomes -- no, delayed or prompt collapse -- and thereby constrain possible EOSs.

We experiment with modifications of the BSSN form of the Einstein field equations (a reformulation of the ADM equations) and demonstrate how these modifications affect the stability of numerical black hole evolution calculations. We use excision to evolve both non-rotating and rotating Kerr-Schild black holes in octant and equatorial symmetry, and without any symmetry assumptions, and obtain accurate and stable simulations for specific angular momenta J/M of up to about 0.9M.

We investigate the dynamical instability of the one-armed spiral m=1 mode in differentially rotating stars by means of 3+1 hydrodynamical simulations in Newtonian gravitation. We find that both a soft equation of state and a high degree of differential rotation in the equilibrium star are necessary to excite a dynamical m=1 mode as the dominant instability at small values of the ratio of rotational kinetic to potential energy, T/|W|. We find that this spiral mode propagates outward from its point of origin near the maximum density at the center to the surface over several central orbital periods. An unstable m=1 mode triggers a secondary m=2 bar mode of smaller amplitude, and the bar mode can excite gravitational waves. As the spiral mode propagates to the surface it weakens, simultaneously damping the emitted gravitational wave signal. This behavior is in contrast to waves triggered by a dynamical m=2 bar instability, which persist for many rotation periods and decay only after a radiation-reaction damping timescale.

We assemble the equations of general relativistic magnetohydrodynamics (MHD) in 3+1 form. These consist of the complete coupled set of Maxwell equations for the electromagnetic field, Einstein's equations for the gravitational field, and the equations of relativistic MHD for a perfectly conducting ideal gas. The adopted form of the equations is suitable for evolving numerically a relativistic MHD fluid in a dynamical spacetime characterized by a strong gravitational field.

We study of the collapse of a magnetized spherical star to a black hole in general relativity theory. The matter and gravitational fields are described by the exact Oppenheimer-Snyder solution for the collapse of a spherical, homogeneous dust ball. We adopt a ``dynamical Cowling approximation'' whereby the matter and the geometry (metric), while highly dynamical, are unaffected by the electromagnetic fields. The matter is assumed to be perfectly conducting and threaded by a dipole magnetic field at the onset of collapse. We determine the subsequent evolution of the magnetic and electric fields without approximation; the fields are determined analytically in the matter interior and numerically in the vacuum exterior. We apply junction conditions to match the electromagnetic fields across the stellar surface. We use this model to experiment with several coordinate gauge choices for handling spacetime evolution characterized by the formation of a black hole and the associated appearance of singularities. These choices range from ``singularity avoiding'' time coordinates to ``horizon penetrating'' time coordinates accompanied by black hole excision. The later choice enables us to integrate the electromagnetic fields arbitrarily far into the future. At late times the longitudinal magnetic field in the exterior has been transformed into a transverse electromagnetic wave; part of the electromagnetic radiation is captured by the hole and the rest propagates outward to large distances. The solution we present for our simple scenario can be used to test codes designed to treat more general evolutions of relativistic MHD fluids flowing in strong gravitational fields in dynamical spacetimes.

Numerical relativity is the most promising tool for theoretically modeling the inspiral and coalescence of neutron star and black hole binaries, which, in turn, are among the most promising sources of gravitational radiation for future detection by gravitational wave observatories. In this article we review numerical relativity approaches to modeling compact binaries. Starting with a brief introduction to the 3+1 decomposition of Einstein's equations, we discuss important components of numerical relativity, including the initial data problem, reformulations of Einstein's equations, coordinate conditions, and strategies for locating and handling black holes on numerical grids. We focus on those approaches which currently seem most relevant for the compact binary problem. We then outline how these methods are used to model binary neutron stars and black holes, and review the current status of inspiral and coalescence simulations.

The merger of binary neutron stars is likely to lead to differentially rotating remnants. In this paper we numerically construct models of differentially rotating neutron stars in general relativity and determine their maximum allowed mass. We model the stars adopting a polytropic equation of state and tabulate maximum allowed masses as a function of differential rotation and stiffness of the equation of state. We also provide a crude argument that yields a qualitative estimate of the effect of stiffness and differential rotation on the maximum allowed mass.

We solve Einstein's field equations coupled to relativistic hydrodynamics in full 3+1 general relativity to evolve astrophysical systems characterized by strong gravitational fields. We model rotating, collapsing and binary stars by idealized polytropic equations of state, with neutron stars as the main application. Our scheme is based on the BSSN formulation of the field equations. We assume adiabatic flow, but allow for the formation of shocks. We determine the appearance of black holes by means of an apparent horizon finder. We introduce several new techniques for integrating the coupled Einstein-hydrodynamics system. For example, we choose our fluid variables so that they can be evolved without employing an artificial atmosphere. We also demonstrate the utility of working in a rotating coordinate system for some problems. We use rotating stars to experiment with several gauge choices for the lapse function and shift vector, and find some choices to be superior to others. We demonstrate the ability of our code to follow a rotating star that collapses from large radius to a black hole. Finally, we exploit rotating coordinates to evolve a corotating binary neutron star system in a quasi-equilibrium circular orbit for more than two orbital periods.

We study the gravitational collapse of a rotating supermassive star (SMS) by means of a (3+1) hydrodynamical simulation in a post-Newtonian (PN) approxi- mation of general relativity. This problem is particularly challenging because of the vast dynamical range in space which must be covered in the course of col- lapse. We evolve a uniformly rotating SMS from the onset of radial instability at R_{p}/M=411, where R_{p} is the proper polar radius of the star and M is the total mass-energy, to the point at which the PN approximation breaks down. We introduce a scale factor and a "comoving" coordinate to handle the large varia- tion in radius during the collapse and focus on the central core. Since T/W, the ratio of the rotational kinetic energy to the gravitational binding energy, is nearly proportional to 1/R_{p} throughout the collapse, the imploding star may ultimately exceed the critical value of T/W for dynamical instability to bar-mode formation. However, for stars rotating uniformly at the onset of col- lapse, we do not find any unstable growth of bars prior to the termination of our simulation. We do find that the collapse is likely to form a supermassive black hole (BH) coherently, with almost all of the matter falling into the hole, leaving very little ejected matter to form a disk. In the absence of nonaxi- symmetric bar formation, the collapse of a uniformly rotating SMS does not lead to appreciable quasi-periodic gravitational wave (GW) emission by the time our integrations terminate. However, the coherent nature of the implosion suggests that rotating SMS collapse will be a promising source of GW bursts. We also expect that, following BH formation, long wavelength quasi-periodic waves will result from quasi-normal ringing. These waves may be detectable by LISA.

We show that a reformulation of the ADM equations in general relativity, which has dramatically improved the stability properties of numerical implementations, has a direct analogue in classical electrodynamics. We numerically integrate both the original and the revised versions of Maxwell's equations, and show that their distinct numerical behavior reflects the properties found in linearized general relativity. Our results shed further light on the stability properties of general relativity, illustrate them in a very transparent context, and may provide a useful framework for further improvement of numerical schemes.

We model the adiabatic inspiral of relativistic binary neutron stars in a quasi-equilibrium (QE) approximation, and compute the gravitational wavetrain from the late phase of the inspiral. We compare corotational and irrotational sequences and find a significant difference in the inspiral rate, which is almost entirely caused by differences in the binding energy. We also compare our results with those of a point-mass post-Newtonian calculation. We illustrate how the late inspiral wavetrain computed with our QE numerical scheme can be matched to the subsequent plunge and merger waveform calculated with a fully relativistic hydrodynamics code.

We evolve a scalar field in a fixed Kerr-Schild background geometry to test simple $(3+1)$-dimensional algorithms for singularity excision. We compare both centered and upwind schemes for handling the shift (advection) terms, as well as different approaches for implementing the excision boundary conditions, for both static and boosted black holes. By first determining the scalar field evolution in a static frame with a $(1+1)$-dimensional code, we obtain the solution to very high precision. This solution then provides a useful testbed for simulations in full $(3+1)$ dimensions. We show that some algorithms which are stable for non-boosted black holes become unstable when the boost velocity becomes high.

We study by computational means the dynamical stability against bar-mode deformation of rapidly and differentially rotating stars in a post-Newtonian approximation of general relativity. We vary the compaction of the star $M/R$ (where $M$ is the gravitational mass and $R$ the equatorial circumferential radius) between 0.01 and 0.05 to isolate the influence of relativistic gravitation on the instability. For compactions in this moderate range, the critical value of $\beta = T/W$ for the onset of the dynamical instability (where $T$ is the rotational kinetic energy and W the gravitational binding energy) slightly decreases from $\sim$ 0.26 to $\sim$ 0.25 with increasing compaction for our choice of the differential rotational law. Combined with our earlier findings based on simulations in full general relativity for stars with higher compaction, we conclude that relativistic gravitation enhances the dynamical bar-mode instability, i.e. the onset of instability sets in for smaller values of $\beta$ in relativistic gravity than in Newtonian gravity. We also find that once a triaxial structure forms after the bar-mode perturbation saturates in dynamically unstable stars, the triaxial shape is maintained, at least for several rotational periods. To check the reliability of our numerical integrations, we verify that the general relativistic Kelvin-Helmholtz circulation is well-conserved, in addition to rest-mass and total mass-energy, linear and angular momentum. Conservation of circulation indicates that our code is not seriously affected by numerical viscosity. We determine the amplitude and frequency of the quasi-periodic gravitational waves emitted during the bar formation process using the quadrupole formula.

Newtonian point mass binaries can be brought into arbitrarily close circular orbits. Neutron stars and black holes, however, are extended, relativistic objects. Both finite size and relativistic effects make very close orbits unstable, so that there exists an innermost stable circular orbit (ISCO). We illustrate the physics of the ISCO in a simple model problem, and review different techniques which have been employed to locate the ISCO in black hole and neutron star binaries. We discuss different assumptions and approximations, and speculate on how differences in the results may be explained and resolved.

A quasi-equilibrium (QE) computational scheme was recently developed in general relativity to calculate the complete gravitational wavetrain emitted during the inspiral phase of compact binaries. The QE method exploits the fact that the the gravitational radiation inspiral timescale is much longer than the orbital period everywhere outside the ISCO. Here we demonstrate the validity and advantages of the QE scheme by solving a model problem in relativistic scalar gravitation theory. By adopting scalar gravitation, we are able to numerically track without approximation the damping of a simple, quasi-periodic radiating system (an oscillating spherical matter shell) to final equilibrium, and then use the exact numerical results to calibrate the QE approximation method. In particular, we calculate the emitted gravitational wavetrain three different ways: by integrating the exact coupled dynamical field and matter equations, by using the scalar-wave monopole approximation formula (corresponding to the quadrupole formula in general relativity), and by adopting the QE scheme. We find that the monopole formula works well for weak field cases, but fails when the fields become even moderately strong. By contrast, the QE scheme remains quite reliable for moderately strong fields, and begins to breakdown only for ultra-strong fields. The QE scheme thus provides a promising technique to construct the complete wavetrain from binary inspiral outside the ISCO, where the gravitational fields are strong, but where the computational resources required to follow the system for more than a few orbits by direct numerical integration of the exact equations are prohibitive.

We study the dynamical stability against bar-mode deformation of rapidly and differentially rotating stars in the first post-Newtonian approximation of general relativity. We vary the compaction of the star $M/R$ (where $M$ is the gravitational mass and $R$ the equatorial circumferential radius) between 0.01 and 0.05 to isolate the influence of relativistic gravitation on the instability. For compactions in this moderate range, the critical value of $\beta \equiv T/W$ for the onset of the dynamical instability (where $T$ is the rotational kinetic energy and $W$ the gravitational binding energy) slightly decreases from $\sim 0.26$ to $\sim 0.25$ with increasing compaction for our choice of the differential rotational law. Combined with our earlier findings based on simulations in full general relativity for stars with higher compaction, we conclude that relativistic gravitation {\em enhances} the dynamical bar-mode instability, i.e. the onset of instability sets in for smaller values of $\beta$ in relativistic gravity than in Newtonian gravity. We also find that once a triaxial structure forms after the bar-mode perturbation saturates in dynamically unstable stars, the triaxial shape is maintained, at least for several rotational periods. To check the reliability of our numerical integrations, we verify that the general relativistic Kelvin-Helmholtz circulation is well-conserved, in addition to rest-mass energy, total mass-energy, linear and angular momentum. Conservation of circulation indicates that our code is not seriously affected by numerical viscosity. We determine the amplitude and frequency of the quasi-periodic gravitational waves emitted during the bar formation process using the quadrupole formula.

We present a new method for generating the nonlinear gravitational wavetrain from the late inspiral (pre-coalescence) phase of a binary neutron star system by means of a numerical evolution calculation in full general relativity. In a prototype calculation, we produce 214 wave cycles from corotating polytropes, representing the final part of the inspiral phase prior to reaching the ISCO. Our method is based on the inequality that the orbital decay timescale due to gravitational radiation is much longer than an orbital period and the approximation that gravitational radiation has little effect on the structure of the stars. We employ quasi-equilibrium sequences of binaries in circular orbit for the matter source in our field evolution code. We compute the gravity-wave energy flux, and, from this, the inspiral rate, at a discrete set of binary separations. From these data, we construct the gravitational waveform as a continuous wavetrain. Finally, we discuss the limitations of our current calculation, planned improvements, and potential applications of our method to other inspiral scenarios.

We study the dynamical stability against bar-mode deformation of rapidly spinning neutron stars with differential rotation. We perform fully relativistic 3D simulations of compact stars with $M/R \geq 0.1$, where $M$ is the total gravitational mass and $R$ the equatorial circumferential radius. We adopt an adiabatic equation of state with adiabatic index $\Gamma=2$. As in Newtonian theory, we find that stars above a critical value of $\beta \equiv T/W$ (where $T$ is the rotational kinetic energy and $W$ the gravitational binding energy) are dynamically unstable to bar formation. For our adopted choices of stellar compaction and rotation profile, the critical value of $\beta = \beta_{dGR}$ is $\sim 0.24-0.25$, only slightly smaller than the well-known Newtonian value $\sim 0.27$ for incompressible Maclaurin spheroids. The critical value depends only very weakly on the degree of differential rotation for the moderate range we surveyed. All unstable stars form bars on a dynamical timescale. Models with sufficiently large $\beta$ subsequently form spiral arms and eject mass, driving the remnant to a dynamically stable state. Models with moderately large $\beta \gtrsim \beta_{dGR}$ do not develop spiral arms or eject mass but adjust to form dynamically stable ellipsoidal-like configurations. If the bar-mode instability is triggered in supernovae collapse or binary neutron star mergers, it could be a strong and observable source of gravitational waves. We determine characteristic wave amplitudes and frequencies.

We introduce a new method to construct solutions to the constraint equations of general relativity describing binary black holes in quasicircular orbit. Black hole pairs with arbitrary momenta can be constructed with a simple method recently suggested by Brandt and Bruegmann, and quasicircular orbits can then be found by locating a minimum in the binding energy along sequences of constant horizon area. This approach produces binary black holes in a "three-sheeted" manifold structure, as opposed to the "two-sheeted" structure in the conformal-imaging approach adopted earlier by Cook. We focus on locating the innermost stable circular orbit and compare with earlier calculations. Our results confirm those of Cook and imply that the underlying manifold structure has a very small effect on the location of the innermost stable circular orbit.

We perform 3D numerical simulations in full general relativity to study the stability of rapidly rotating, supramassive neutron stars at the mass-shedding limit to dynamical collapse. We adopt an adiabatic equation of state with $\Gamma = 2$ and focus on uniformly rotating stars. We find that the onset of dynamical instability along mass-shedding sequences nearly coincides with the onset of secular instability. Unstable stars collapse to rotating black holes within about one rotation period. We also study the collapse of stable stars which have been destabilized by pressure depletion (e.g. via a phase transition) or mass accretion. In no case do we find evidence for the formation of massive disks or any ejecta around the newly formed Kerr black holes, even though the progenitors are rapidly rotating.

We construct relativistic equilibrium models of differentially rotating neutron stars and show that they can support significantly more mass than their nonrotating or uniformly rotating counterparts. We dynamically evolve such ``hypermassive'' models in full general relativity and show that there do exist configurations which are dynamically stable against radial collapse and bar formation. Our results suggest that the remnant of binary neutron star coalescence may be temporarily stabilized by differential rotation, leading to delayed collapse and a delayed gravitational wave burst.

We launch a fully relativistic study of the formation of supermassive black holes via the collapse of supermassive stars. Here we initiate our investigation by analyzing the secular evolution of supermassive stars up to the onset of dynamical instability and collapse. We focus on the effects of rotation, assumed uniform, and general relativity. We identify the critical configuration at which radial instability sets in and determine its structure in detail. We show that the key nondimensional ratios R/M, T/|W| and J/M^2 (T is the rotational kinetic energy and W is the gravitational potential energy) for this critical configuration are universal numbers, independent of the mass, spin, radius or history of the star. We compare results from an approximate, analytic treatment with a fully relativistic, numerical calculation and find good agreement. We solve analytically for the time evolution of these parameters up to the onset of instability. Cooling by photon radiation drives the evolution, which is accompanied by mass, angular momentum and entropy loss. The critical configuration serves as initial data for a future relativistic, hydrodynamical, 3D simulation of the collapse of an unstable supermassive star. Since this implosion starts from a universal critical configuration, the collapse is also uniquely determined and should produce a universal gravitational waveform. In this paper we briefly speculate on the possible outcome of this collapse and asses to what extent it offers a promising route to forming a supermassive black hole.

We determine the effect of rotation on the luminosity of supermassive stars. We apply the Roche model to calculate analytically the emitted radiation from a uniformly rotating, radiation-dominated supermassive configuration. We find that the luminosity at maximum rotation, when mass at the equator orbits at the Kepler period, is reduced by ~36% below the usual Eddington luminosity from the corresponding nonrotating star. A supermassive star is believed to evolve in a quasistationary manner along such a maximally rotating ``mass-shedding'' sequence before reaching the point of dynamical instability; hence this reduced luminosity determines the evolutionary timescale. Our result therefore implies that the lifetime of a supermassive star prior to dynamical collapse is ~56% longer than the value typically estimated by employing the usual Eddington luminosity.

We report on our numerical implementation of fully relativistic hydrodynamics coupled to Einstein's field equations in three spatial dimensions. We briefly review several steps in our code development, including our recasting of Einstein's equations and several tests which demonstrate its advantages for numerical integrations. We outline our implementation of relativistic hydrodynamics, and present numerical results for the evolution of both stable and unstable Oppenheimer-Volkov equilibrium stars, which represent a very promising first test of our code.

We present a new set of 3.5 Post-Newtonian equations in which Newtonian hydrodynamics is coupled to the nonconservative effects of gravitational radiation emission. Our formalism differs in two significant ways from a similar 3.5 Post-Newtonian approach proposed by Blanchet (1993, 1997). Firstly we concentrate only on the radiation-reaction effects produced by a time-varying mass-current quadrupole $S_{ij}$. Secondly, we adopt a gauge in which the radiation-reaction force densities depend on the fourth time derivative of $S_{ij}$, rather than on the fifth, as in Blanchet's approach. This difference makes our formalism particularly well-suited to numerical implementation and could prove useful in performing fully numerical simulations of the recently discovered $r$-mode instability for rotating neutron stars subject to axial perturbations.

We include matter sources in Einstein's field equations and show that our recently proposed 3+1 evolution scheme can stably evolve strong-field solutions. We insert in our code known matter solutions, namely the Oppenheimer-Volkoff solution for a static star and the Oppenheimer-Snyder solution for homogeneous dust sphere collapse to a black hole, and evolve the gravitational field equations. We find that we can evolve stably static, strong-field stars for arbitrarily long times and can follow dust sphere collapse accurately well past black hole formation. These tests are useful diagnostics for fully self-consistent, stable hydrodynamical simulations in 3+1 general relativity. Moreover, they suggest a successive approximation scheme for determining gravitational waveforms from strong-field sources dominated by longitudinal fields, like binary neutron stars: approximate quasi-equilibrium models can serve as sources for the transverse field equations, which can be evolved without having to re-solve the hydrodynamical equations (``hydro without hydro'').

Many numerical codes now under development to solve Einstein's equations of general relativity in 3+1 dimensional spacetimes employ the standard ADM form of the field equations. This form involves evolution equations for the raw spatial metric and extrinsic curvature tensors. Following Shibata and Nakamura, we modify these equations by factoring out the conformal factor and introducing three ``connection functions''. The evolution equations can then be reduced to wave equations for the conformal metric components, which are coupled to evolution equations for the connection functions. We evolve small amplitude gravitational waves and make a direct comparison of the numerical performance of the modified equations with the standard ADM equations. We find that the modified form exhibits much improved stability.

We perform simulations of relativistic binary stars in post-Newtonian gravity to investigate their dynamical stability prior to merger against gravitational collapse in a tidal field. In general, our equations are only strictly accurate to first post-Newtonian order, but they recover full general relativity for spherical, static stars. We study both corotational and irrotational binary configurations of identical stars in circular orbits. We adopt a soft, adiabatic equation of state with $\Gamma = 1.4$, for which the onset of instability occurs at a sufficiently small value of the compaction $M/R$ that a post-Newtonian approximation is quite accurate. For such a soft equation of state there is no innermost stable circular orbit, so that we can study arbitrarily close binaries. This choice still allows us to study all the qualitative features exhibited by any adiabatic equation of state regarding stability against gravitational collapse. We demonstrate that, independent of the internal stellar velocity profile, the tidal field from a binary companion stabilizes a star against gravitational collapse.

We study neutrino emission from the remnant of an inspiraling binary neutron star following 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. We find 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. We illustrate these findings with a post-Newtonian, ellipsoidal model calculation. Simple arguments suggest that the remnant may form 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.

We perform fully relativistic calculations of binary neutron stars in corotating, circular orbit. While Newtonian gravity allows for a strict equilibrium, a relativistic binary system emits gravitational radiation, causing the system to lose energy and slowly spiral inwards. However, since inspiral occurs on a time scale much longer than the orbital period, we can treat the binary to be in quasiequilibrium. In this approximation, we integrate a subset of the Einstein equations coupled to the relativistic equation of hydrostatic equilibrium to solve the initial value problem for binaries of arbitrary separation. We adopt a polytropic equation of state to determine the structure and maximum mass of neutron stars in close binaries for polytropic indices n=1, 1.5 and 2. We construct sequences of constant rest-mass and locate turning points along energy equilibrium curves to identify the onset of orbital instability. In particular, we locate the innermost stable circular orbit (ISCO) and its angular velocity. We construct the first contact binary systems in full general relativity. These arise whenever the equation of state is sufficiently soft >= 1.5. A radial stability analysis reveals no tendency for neutron stars in close binaries to collapse to black holes prior to merger.

We analyze the stability of relativistic, quasi-equilibrium binary neutron stars in synchronous circular orbit. We explore stability against radial collapse to black holes prior to merger, and against orbital plunge. We apply theorems based on turning points along uniformly rotating sequences of constant angular momentum and rest mass to locate the onset of secular instabilities. We find that inspiraling binary neutron stars are stable against radial collapse to black holes all the way down to the innermost stable circular orbit.

We perform fully relativistic calculations of binary neutron stars in quasi-equilibrium circular orbits. We integrate Einstein's equations together with the relativistic equation of hydrostatic equilibrium to solve the initial value problem for equal-mass binaries of arbitrary separation. We construct sequences of constant rest mass and identify the innermost stable circular orbit and its angular velocity. We find that the quasi-equilibrium maximum allowed mass of a neutron star in a close binary is slightly larger than in isolation.

We present a numerical scheme that solves the initial value problem in full general relativity for a binary neutron star in quasi-equilibrium. While Newtonian gravity allows for a strict equilibrium, a relativistic binary system emits gravitational radiation, causing the system to lose energy and slowly spiral inwards. However, since inspiral occurs on a time scale much longer than the orbital period, we can adopt a quasi-equilibrium approximation. In this approximation, we integrate a subset of the Einstein equations coupled to the equations of relativistic hydrodynamics to solve the initial value problem for binaries of arbitrary separation, down to the innermost stable orbit.