The recent discovery by LIGO/Virgo of a merging binary having a $\sim 23 M_\odot$ black hole and a $\sim 2.6 M_\odot$ compact companion has triggered a debate regarding the nature of the secondary, which falls into the so-called mass gap. Here we explore some consequences of the assumption that the secondary was a neutron star (NS). We show with concrete examples of heretofore viable equations of state (EOSs) that rapid uniform rotation is certainly adequate to explain the existence of a stable $\sim 2.6 M_\odot$ NS for moderately stiff EOS but may not be adequate for soft EOSs. Moreover, rotation may not even be necessary for sufficiently stiff ones. In particular, several soft EOSs favored by GW170817 and with maximum spherical masses of $\sim 2.1 M_\odot$ \textit{cannot} be invoked to explain this object as a uniformly rotating NS. Absolute upper mass limits on the maximum spherical NS derived from GW170817 immediately suggest that this unknown compact companion can be a slowly or even a nonrotating NS. Stiff EOSs that can achieve this are neither rejected nor favored by GW170817, and they are in accord with the results of NICER.

A wide range of mechanisms predict present-day s-wave dark matter (DM) annihilation cross-sections that are orders of magnitude below current experimental sensitivity. We explore the capability of DM density spikes around the Milky Way's supermassive black hole to probe such faint signals of DM annihilations, considering a range of possible spike and halo distributions. As an exemplar of a theory with a suppressed s-wave annihilation cross-section, we consider a hidden sector axion portal model of DM. In this model, the leading contribution to the annihilation cross-section in the early universe is p-wave, while s-wave annihilations occur at higher order in the coupling constant. We provide a unified treatment of DM freezeout in this model including both s- and p-wave annihilations and analytically determine the photon spectrum for the dominant DM annihilation process in the universe today. We find that Fermi and H.E.S.S. observations of the Galactic Center offer excellent sensitivity to this model over a wide range of parameter space, with prospects depending sensitively on the properties of the DM spike as well as the central halo.

Recently, we have shown that dynamically stable ergostar solutions (equilibrium neutron stars that contain an ergoregion) with a compressible and causal equation of state exist [A. Tsokaros, M. Ruiz, L. Sun, S. L. Shapiro, and K. Ury\=u, Phys. Rev. Lett. 123, 231103 (2019)]. These stars are hypermassive, differentially rotating, and highly compact. In this work, we make a systematic study of equilibrium models in order to locate the position of ergostars in parameter space. We adopt four equations of state that differ in the matching density of a maximally stiff core. By constructing a large number of models both with uniform and differential rotation of different degrees, we identify the parameters for which ergostars appear. We find that the most favorable conditions for the appearance of dynamically stable ergostars are a significant finite density close to the surface of the star (i.e., similar to self-bound quark stars) and a small degree of differential rotation.

Binary neutron star (NSNS) mergers can be sources of gravitational waves coincident with electromagnetic counterpart emission. To solidify their role as multimessenger sources, we present fully 3D, general relativistic, magnetohydrodynamic simulations of spinning NSNSs initially on quasicircular orbits that merge and undergo delayed collapse to a black hole (BH). The NSNSs consist of two identical stars modeled as $\Gamma=2$ polytropes with spin $\chi_{NS}= 0.36$ aligned along the direction of the total orbital angular momentum $L$. Each star is initially threaded by a dynamical unimportant interior dipole B-field. The field is extended into the exterior where a nearly force-free magnetosphere resembles that of a pulsar. The magnetic dipole moment $\mu$ is either aligned or perpendicular to $L$ and has the same initial magnitude for each orientation. For comparison, we also impose symmetry across the orbital plane in one case where $\mu$ in both stars is aligned along $L$. We find that the lifetime of the transient hypermassive neutron star remnant, the jet launching time, and the ejecta are very sensitive to the B-field orientation. By contrast, the physical properties of the BH + disk remnant, such as the mass and spin of the BH, the accretion rate, and the electromagnetic luminosity, are roughly independent of the initial B-field orientation. In addition, we find imposing symmetry across the orbital plane does not play a significant role in the final outcome of the mergers. Our results show that an incipient jet emerges only when the seed B-field has a sufficiently large-scale poloidal component aligned to $L$. The lifetime [$\Delta t\gtrsim 140(M_{NS}/1.625M_\odot)\rm ms$] and Poynting luminosities [$L_{EM}\simeq 10^{52}$erg/s] of the jet, when it forms, are consistent with typical short gamma-ray bursts, as well as with the Blandford--Znajek mechanism for launching jets.

Can one distinguish a binary black hole undergoing a merger from a binary neutron star if the individual compact companions have masses that fall inside the so-called mass gap of $3-5\ M_\odot$? For neutron stars, achieving such masses typically requires extreme compactness and in this work we present initial data and evolutions of binary neutron stars initially in quasiequilibrium circular orbits having a compactness $C=0.336$. These are the most compact, nonvacuum, quasiequilibrium binary objects that have been constructed and evolved to date, including boson stars. The compactness achieved is only slightly smaller than the maximum possible imposed by causality, $C_{\rm max}=0.355$, which requires the sound speed to be less than the speed of light. By comparing the emitted gravitational waveforms from the late inspiral to merger and postmerger phases between such a binary neutron star vs a binary black hole of the same total mass we identify concrete measurements that serve to distinguish them. With that level of compactness, the binary neutron stars exhibit no tidal disruption up until merger, whereupon a prompt collapse is initiated even before a common core forms. Within the accuracy of our simulations the black hole remnants from both binaries exhibit ringdown radiation that is not distinguishable from a perturbed Kerr spacetime. However, their inspiral leads to phase differences of the order of $\sim 5$ rad over an $\sim 81$ km separation (1.7 orbits) while typical neutron stars exhibit phase differences of $\geq 20$ rad. Although a difference of $\sim 5$ rad can be measured by current gravitational wave laser interferometers (e.g., aLIGO/Virgo), uncertainties in the individual masses and spins will likely prevent distinguishing such compact, massive neutron stars from black holes.

Multi-messenger astrophysics is a fast-growing, interdisciplinary field that combines data, which vary in volume and speed of data processing, from many different instruments that probe the Universe using different cosmic messengers: electromagnetic waves, cosmic rays, gravitational waves and neutrinos. In this Expert Recommendation, we review the key challenges of real-time observations of gravitational wave sources and their electromagnetic and astroparticle counterparts, and make a number of recommendations to maximize their potential for scientific discovery. These recommendations refer to the design of scalable and computationally efficient machine learning algorithms; the cyber-infrastructure to numerically simulate astrophysical sources, and to process and interpret multi-messenger astrophysics data; the management of gravitational wave detections to trigger real-time alerts for electromagnetic and astroparticle follow-ups; a vision to harness future developments of machine learning and cyber-infrastructure resources to cope with the big-data requirements; and the need to build a community of experts to realize the goals of multi-messenger astrophysics.

We construct the first dynamically stable ergostars (equilibrium neutron stars that contain an ergoregion) for a compressible, causal equation of state. We demonstrate their stability by evolving both strict and perturbed equilibrium configurations in full general relativity for over a hundred dynamical timescales ($\gtrsim 30$ rotational periods) and observing their stationary behavior. This stability is in contrast to earlier models which prove radially unstable to collapse. Our solutions are highly differentially rotating hypermassive neutron stars with a corresponding spherical compaction of $C=0.3$. Such ergostars can provide new insights into the geometry of spacetimes around highly compact, rotating objects and on the equation of state at supranuclear densities. Ergostars may form as remnants of extreme binary neutron star mergers and possibly provide another mechanism for powering short gamma-ray bursts.

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.

We perform long-term simulations of spinning binary neutron stars, with our highest dimensionless spin being $\chi \sim 0.32$. To assess the importance of spin during the inspiral we vary the spin, and also use two equations of state, one that consists of plain nuclear matter and produces compact stars (SLy), and a hybrid one that contains both nuclear and quark matter and leads to larger stars (ALF2). Using high resolution that has grid spacing $\Delta x\sim 98$ m on the finest refinement level, we find that the effects of spin in the phase evolution of a binary system can be larger than the one that comes from tidal forces. Our calculations demonstrate explicitly that although tidal effects are dominant for small spins ($\lesssim 0.1$), this is no longer true when the spins are larger, but still much smaller than the Keplerian limit.

Events GW170817 and GRB 170817A provide the best confirmation so far that compact binary mergers where at least one of the companions is a neutron star (NS) can be the progenitors of short gamma-ray bursts (sGRBs). An open question for GW170817 remains the values and impact of the initial NS spins. The initial spins could possibly affect the remnant black hole (BH) mass and spin, the remnant disk and the formation and lifetime of a jet and its luminosity. Here we summarize our general relativistic magnetohydrodynamic simulations of spinning, NS binaries undergoing merger and delayed collapse to a BH. The binaries consist of two identical NSs, modeled as $\Gamma=2$ polytropes, in quasicircular orbit, each with spins $\chi_{\rm{NS}}=-0.053,\,0,\,0.24$, or $0.36$. The stars are endowed initially with a dipolar magnetic field extending from the interior into the exterior, as in a radio pulsar. Following merger, the redistribution of angular momentum by magnetic braking and magnetic turbulent viscosity in the hypermassive neutron star (HMNS) remnant, along with the loss of angular momentum due to gravitational radiation, induce the formation of a massive, nearly uniformly rotating inner core surrounded by a magnetized Keplerian disk-like envelope. The HMNS eventually collapses to a BH, with spin $a/M_{\rm BH} \simeq 0.78$ independent of the initial spin of the NSs, surrounded by a magnetized accretion disk. The larger the initial NS spin the heavier the disk. After $\Delta t\sim 3000-4000 M \sim 45-60(M_{\rm NS}/1.625M_\odot)\rm ms$ following merger, a mildly relativistic jet is launched. The lifetime of the jet [$\Delta t\sim 100-140(M_{\rm NS}/1.625M_\odot)\rm ms$] and its outgoing Poynting luminosity [$L_{\rm EM}\sim 10^{51.5\pm 1}\rm erg/s$] are consistent with typical sGRBs, as well as with the Blandford--Znajek mechanism for launching jets and their associated Poynting luminosities.

Fragmentation of highly differentially rotating massive stars that undergo collapse has been suggested as a possible channel for binary black hole formation. Such a scenario could explain the formation of the new population of massive black holes detected by the LIGO/VIRGO gravitational wave laser interferometers. We probe that scenario by performing general relativistic magnetohydrodynamic simulations of differentially rotating massive stars supported by thermal radiation pressure plus a gas pressure perturbation. The stars are initially threaded by a dynamically weak, poloidal magnetic field confined to the stellar interior. We find that magnetic braking and turbulent viscous damping via magnetic winding and the magnetorotational instability in the bulk of the star redistribute angular momentum, damp differential rotation and induce the formation of a massive and nearly uniformly rotating inner core surrounded by a Keplerian envelope. The core + disk configuration evolves on a secular timescale and remains in quasi-stationary equilibrium until the termination of our simulations. Our results suggest that the high degree of differential rotation required for $m=2$ seed density perturbations to trigger gas fragmentation and binary black hole formation is likely to be suppressed during the normal lifetime of the star prior to evolving to the point of dynamical instability to collapse. Other cataclysmic events, such as stellar mergers leading to collapse, may therefore be necessary to reestablish sufficient differential rotation and density perturbations to drive nonaxisymmetric modes leading to binary black hole formation.

As of today, we have directly detected exactly one source in both gravitational waves (GWs) and electromagnetic (EM) radiation, the binary neutron star merger GW170817, its associated gamma-ray burst GRB170817A, and the subsequent kilonova SSS17a/AT 2017gfo. Within ten years, we will detect hundreds of events, including new classes of events such as neutron-star-black-hole mergers, core-collapse supernovae, and almost certainly something completely unexpected. As we build this sample, we will explore exotic astrophysical topics ranging from nucleosynthesis, stellar evolution, general relativity, high-energy astrophysics, nuclear matter, to cosmology. The discovery potential is extraordinary, and investments in this area will yield major scientific breakthroughs. Here we outline some of the most exciting scientific questions that can be answered by combining GW and EM observations.

We present a new initial data formulation to solve the full set of Einstein equations for spacetimes that contain a black hole under general conditions. The method can be used to construct complete initial data for spacetimes (the full metric) that contain a black hole. Contrary to most current studies the formulation requires minimal assumptions. For example, rather than imposing the form of the spatial conformal metric we impose 3 gauge conditions adapted to the coordinates describing the system under consideration. For stationary, axisymmetric spacetimes our method yields Kerr-Schild black holes in vacuum and rotating equilibrium neutron stars. We demonstrate the power of our new method by solving for the first time the whole system of Einstein equations for a nonaxisymmetric, self-gravitating torus in the presence of a black hole. The black hole has dimensionless spin $J_{\rm bh}/M_{\rm bh}^2=0.9918$, a rotation axis tilted at a $30^\circ$ angle with respect to the angular momentum of the disk, and a mass of $\sim 1/5$ of the disk.

This report provides an overview of recent work that harnesses the Big Data Revolution and Large Scale Computing to address grand computational challenges in Multi-Messenger Astrophysics, with a particular emphasis on real-time discovery campaigns. Acknowledging the transdisciplinary nature of Multi-Messenger Astrophysics, this document has been prepared by members of the physics, astronomy, computer science, data science, software and cyberinfrastructure communities who attended the NSF-, DOE- and NVIDIA-funded "Deep Learning for Multi-Messenger Astrophysics: Real-time Discovery at Scale" workshop, hosted at the National Center for Supercomputing Applications, October 17-19, 2018. Highlights of this report include unanimous agreement that it is critical to accelerate the development and deployment of novel, signal-processing algorithms that use the synergy between artificial intelligence (AI) and high performance computing to maximize the potential for scientific discovery with Multi-Messenger Astrophysics. We discuss key aspects to realize this endeavor, namely (i) the design and exploitation of scalable and computationally efficient AI algorithms for Multi-Messenger Astrophysics; (ii) cyberinfrastructure requirements to numerically simulate astrophysical sources, and to process and interpret Multi-Messenger Astrophysics data; (iii) management of gravitational wave detections and triggers to enable electromagnetic and astro-particle follow-ups; (iv) a vision to harness future developments of machine and deep learning and cyberinfrastructure resources to cope with the scale of discovery in the Big Data Era; (v) and the need to build a community that brings domain experts together with data scientists on equal footing to maximize and accelerate discovery in the nascent field of Multi-Messenger Astrophysics.

For isentropic fluids, dynamical evolution of a binary system conserves the baryonic mass and circulation; therefore, sequences of constant rest mass and constant circulation are of particular importance. In this work, we present the extension of our Compact Object CALculator (\cocal{}) code to compute such quasiequilibria and compare them with the well-known corotating and irrotational sequences, the latter being the simplest, zero-circulation case. The circulation as a measure of the spin for a neutron star in a binary system has the advantage of being exactly calculable since it is a local quantity. To assess the different measures of spin, such as the angular velocity of the star, the quasilocal, dimensionless spin parameter $J/M^2$, or the circulation $\mathcal{C}$, we first compute sequences of single, uniformly rotating stars and describe how the different spin diagnostics are related to each other. The connection to spinning binary systems is accomplished through the concept of circulation and the use of the constant rotational velocity formulation. Finally, we explore a modification of the latter formulation that naturally leads to differentially rotating binary systems.

Black hole-neutron star (BHNS) mergers are one of the most promising targets for multimessenger astronomy. Using general relativistic magnetohydrodynamic simulations of BHNS undergoing merger we showed that a magnetically--driven jet can be launched by the remnant if the NS is endowed with a dipole B field extending from the interior into the exterior as in a radio pulsar. These self-consistent studies considered a BHNS system with mass ratio $q=3:1$, BH spin $a/M_{BH}=0.75$ aligned with the total orbital angular momentum (OAM), and a NS that is irrotational, threaded by an aligned B field, and modeled by an $\Gamma$--law equation of state with $\Gamma=2$. Here, as a crucial step in establishing BHNS systems as viable progenitors of central engines that power short gamma--ray bursts (sGRBs) and thereby solidify their role as multimessenger sources, we survey different BHNS configurations that differ in BH spin ($a/M_{BH} =-0.5,\,0,\,0.5,\,0.75$), in the mass ratio ($q=3:1$ and $q=5:1$), and in the orientation of the B field (aligned and tilted by $90^\circ$ with respect to the OAM). We find that by $\Delta t\sim 3500M-4000M \sim 88(M_{NS}/1.4M_\odot){\rm ms}-100(M_{NS}/1.4M_\odot)\rm ms$ after the peak gravitational wave signal a jet is launched in the cases where the initial BH spin is $a/M_{BH}= 0.5$ or $0.75$. The lifetime of the jets[$\Delta t\sim 0.5(M_{NS}/1.4M_\odot){\rm s-0.7}(M_{NS}/1.4M_\odot)\rm s$] and their Poynting luminosities [$L_{jet}\sim 10^{51\pm 1}\rm erg/s$] are consistent with sGRBs, as well as with the Blandford--Znajek mechanism. By the time we terminate our simulations, we do not observe either an outflow or a large-scale B field collimation in the other configurations we simulate. These results suggest that future multimessenger detections from BHNSs are more likely produced by binaries with highly spinning BH companions and small tilt-angle B fields.

Collapsing supermassive stars (SMSs) with masses $M \gtrsim 10^{4-6}M_\odot$ have long been speculated to be the seeds that can grow and become supermassive black holes (SMBHs). We previously performed GRMHD simulations of marginally stable magnetized $\Gamma = 4/3$ polytropes uniformly rotating at the mass-shedding limit to model the direct collapse of SMSs. These configurations are supported entirely by thermal radiation pressure and model SMSs with $M \gtrsim 10^{6}M_\odot$. We found that around $90\%$ of the initial stellar mass forms a spinning black hole (BH) surrounded by a massive, hot, magnetized torus, which eventually launches an incipient jet. Here we perform GRMHD simulations of $\Gamma \gtrsim 4/3$, polytropes to account for the perturbative role of gas pressure in SMSs with $M \lesssim 10^{6}M_\odot$. We also consider different initial stellar rotation profiles. The stars are initially seeded with a dynamically weak dipole magnetic field that is either confined to the stellar interior or extended from its interior into the stellar exterior. We find that the mass of the BH remnant is $90\%-99\%$ of the initial stellar mass, depending sharply on $\Gamma -4/3$ as well as on the initial stellar rotation profile. After $t\sim 250-550M\approx 1-2\times 10^3(M/10^6M_\odot)$s following the BH formation, a jet is launched and it lasts for $\sim 10^4-10^5(M/10^6M_\odot)$s, consistent with the duration of long gamma-ray bursts. Our results suggest that the Blandford-Znajek mechanism powers the jet. They are also in agreement with our proposed universal model that estimates accretion rates and luminosities that characterize magnetized BH-disk remnant systems that launch a jet. This model helps explain why the outgoing luminosities for vastly different BH-disk formation scenarios all reside within a narrow range ($\sim 10^{52 \pm 1} \rm erg/s$), roughly independent of $M$.

We adopt the fluid conduction approximation to study the evolution of spherical star clusters and self-interacting dark matter (SIDM) halos. We also explore the formation and dynamical impact of density cusps that arise in both systems due to the presence of a massive, central black hole. The large N-body, self-gravitating systems we treat are "weakly-collisional": the mean free time between star or SIDM particle collisions is much longer than their characteristic crossing (dynamical) time scale, but shorter than the system lifetime. The fluid conduction model reliably tracks the "gravothermal catastrophe" in star clusters and SIDM halos without black holes. For a star cluster with a massive, central black hole, this approximation reproduces the familiar Bahcall-Wolf quasistatic density cusp for the stars bound to the black hole and shows how the cusp halts the "gravothermal catastrophe" and causes the cluster to re-expand. An SIDM halo with an initial black hole central density spike that matches onto to an exterior NFW profile relaxes to a core-halo structure with a central density cusp determined by the velocity dependence of the SIDM interaction cross section. The success and relative simplicity of the fluid conduction approach in evolving such "weakly-collisional", quasiequilibrium Newtonian systems motivates its extension to relativistic systems. We present a general relativistic extension here.

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.

Recent numerical simulations in general relativistic magnetohydrodynamics (GRMHD) provide useful constraints for the interpretation of the GW170817 discovery. Combining the observed data with these simulations leads to a bound on the maximum mass of a cold, spherical neutron star (the TOV limit): ${M_{\rm max}^{\rm sph}}\lesssim 2.74/\beta$, where $\beta$ is the ratio of the maximum mass of a uniformly rotating neutron star (the supramassive limit) over the maximum mass of a nonrotating star. Causality arguments allow $\beta$ to be as high as $1.27$, while most realistic candidate equations of state predict $\beta$ to be closer to $1.2$, yielding ${M_{\rm max}^{\rm sph}}$ in the range $2.16-2.28 M_\odot$. A minimal set of assumptions based on these simulations distinguishes this analysis from previous ones, but leads to a similar estimate. There are caveats, however, and they are enumerated and discussed. The caveats can be removed by further simulations and analysis to firm up the basic argument.

We perform magnetohydrodynamic simulations in full general relativity of disk accretion onto nonspinning black hole binaries with mass ratio 36:29. We survey different disk models which differ in their scale height, total size and magnetic field to quantify the robustness of previous simulations on the initial disk model. Scaling our simulations to LIGO GW150914 we find that such systems could explain possible gravitational wave and electromagnetic counterparts such as the Fermi GBM hard X-ray signal reported 0.4s after GW150915 ended. Scaling our simulations to supermassive binary black holes, we find that observable flow properties such as accretion rate periodicities, the emergence of jets throughout inspiral, merger and post-merger, disk temperatures, thermal frequencies, and the time-delay between merger and the boost in jet outflows that we reported in earlier studies display only modest dependence on the initial disk model we consider here.

Inspiraling and merging binary neutron stars are not only important source of gravitational waves, but also promising candidates for coincident electromagnetic counterparts. These systems are thought to be progenitors of short gamma-ray bursts (sGRBs). We have shown previously that binary neutron star mergers that undergo {\it delayed} collapse to a black hole surrounded by a {\it weighty} magnetized accretion disk can drive magnetically-powered jets. We now perform magnetohydrodynamic simulations in full general relativity of binary neutron stars mergers that undergo {\it prompt} collapse to explore the possibility of jet formation from black hole-{\it light} accretion disk remnants. We find that after $t-t_{\rm BH}\sim 26(M_{\rm NS}/1.8M_\odot)$ms [$M_{\rm NS}$ is the ADM mass] following prompt black hole formation, there is no evidence of mass outflow or magnetic field collimation. The rapid formation of the black hole following merger prevents magnetic energy from approaching force-free values above the magnetic poles, which is required for the launching of a jet by the usual Blandford--Znajek mechanism. Detection of gravitational waves in coincidence with sGRBs may provide constraints on the nuclear equation of state (EOS): the fate of an NSNS merger--delayed or prompt collapse, and hence the appearance or nonappearance of an sGRB--depends on a critical value of the total mass of the binary, and this value is sensitive to the EOS.

We perform MHD simulations in full GR of uniformly rotating stars that are marginally unstable to collapse. Our simulations model the direct collapse of supermassive stars (SMSs) to seed black holes (BHs) that can grow to become the supermassive BHs at the centers of quasars and AGNs. They also crudely model the collapse of massive Pop III stars to BHs, which could power a fraction of distant, long gamma-ray bursts (GRBs). The initial stellar models we adopt are $\Gamma = 4/3$ polytropes seeded with a dynamically unimportant dipole magnetic field (B field). We treat initial B-field configurations either confined to the stellar interior or extending out from the interior into the stellar exterior. The BH formed following collapse has mass $M_{BH} \simeq 0.9M$ (where $M$ is the mass of the initial star) and spin $a_{BH}/M_{BH}\simeq 0.7$. A massive, hot, magnetized torus surrounds the remnant BH. At $\Delta t\sim 400-550M\approx 2000-2700(M/10^6M_\odot)$s following the gravitational wave (GW) peak amplitude, an incipient jet is launched. The disk lifetime is $\Delta t\sim 10^5(M/10^6M_\odot)$s, and the jet luminosity is $L_{EM}\sim 10^{51-52}$ ergs/s. If $\gtrsim 1-10\%$ of this power is converted into gamma rays, SWIFT and FERMI could potentially detect these events out to large redshifts $z\sim 20$. Thus, SMSs could be sources of ultra-long GRBs and massive Pop III stars could be the progenitors that power a fraction of the long GRBs observed at redshift $z \sim 5-8$. GWs are copiously emitted during the collapse, and peak at $\sim 15(10^6 M_{\odot}/M)\rm mHz$ ($\sim 0.15(10^4 M_{\odot}/M)\rm Hz$), i.e., in the LISA (DECIGO/BBO) band; optimally oriented SMSs could be detectable by LISA (DECIGO/BBO) at $z \lesssim 3$ ($z \lesssim 11$). Hence $10^4 M_{\odot}$ SMSs collapsing at $z\sim 10$ are promising multimessenger sources of coincident gravitational and electromagnetic waves.

We have performed magnetohydrodynamic simulations in general relativity of binary neutron star and binary black hole-neutron star mergers, as well as the magnetorotational collapse of supermassive stars. In many cases the outcome is a spinnng black hole (BH) immersed in a magnetized disk, with a jet emanating from the poles of the BH. While their formation scenarios differ and their BH masses, as well as their disk masses, densities, and magnetic field strengths, vary by orders of magnitude, these features conspire to generate jet Poynting luminosities that all lie in the same, narrow range of $\sim 10^{52\pm1}~{\rm erg ~s^{-1}}$. A similar result applies to their BH accretion rates upon jet launch, which is $\sim 0.1-10~{\rm M_{\odot}~s^{-1}}$. We provide a simple model that explains these unanticipated findings. Interestingly, these luminosities reside in the same narrow range characterizing the observed luminosity distributions of over 400 short and long GRBs with distances inferred from spectroscopic redshifts or host galaxies. This result, together with the GRB lifetimes predicted by the model, supports the belief that a compact binary merger is the progenitor of an SGRB, while a massive, stellar magnetorotational collapse is the progenitor of an LGRB.

Targets for ground-based gravitational wave interferometers include continuous, quasiperiodic sources of gravitational radiation, such as isolated, spinning neutron stars. In this work we perform evolution simulations of uniformly rotating, triaxially deformed stars, the compressible analogues in general relativity of incompressible, Newtonian Jacobi ellipsoids. We investigate their stability and gravitational wave emission. We employ five models, both normal and supramassive, and track their evolution with different grid setups and resolutions, as well as with two different evolution codes. We find that all models are dynamically stable and produce a strain that is approximately one-tenth the average value of a merging binary system. We track their secular evolution and find that all our stars evolve towards axisymmetry, maintaining their uniform rotation, kinetic energy, and angular momentum profiles while losing their triaxiality.

We perform magnetohydrodynamic simulations in full general relativity (GRMHD) of quasi-circular, equal-mass, binary neutron stars that undergo merger. The initial stars are irrotational, $n=1$ polytropes and are magnetized. We explore two types of magnetic-field geometries: one where each star is endowed with a dipole magnetic field extending from the interior into the exterior, as in a pulsar, and the other where the dipole field is initially confined to the interior. In both cases the adopted magnetic fields are initially dynamically unimportant. The merger outcome is a hypermassive neutron star that undergoes delayed collapse to a black hole (spin parameter $a/M_{\rm BH} \sim 0.74$) immersed in a magnetized accretion disk. About $4000M \sim 60(M_{\rm NS}/1.625M_\odot)$ ms following merger, the region above the black hole poles becomes strongly magnetized, and a collimated, mildly relativistic outflow --- an incipient jet --- is launched. The lifetime of the accretion disk, which likely equals the lifetime of the jet, is $\Delta t \sim 0.1 (M_{\rm NS}/1.625M_\odot)$ s. In contrast to black hole--neutron star mergers, we find that incipient jets are launched even when the initial magnetic field is confined to the interior of the stars.

We reinvestigate the effect of annihilations on the distribution of collisionless dark matter (DM) in a spherical density spike around a massive black hole. We first construct a very simple, pedagogic, analytic model for an isotropic phase space distribution function that accounts for annihilation and reproduces the "weak cusp" found by Vasiliev for DM deep within the spike and away from its boundaries. The DM density in the cusp varies as $r^{-1/2}$ for $s$-wave annihilation, where $r$ is the distance from the central black hole, and is not a flat "plateau" profile. We then extend this model by incorporating a loss cone that accounts for the capture of DM particles by the hole. The loss cone is implemented by a boundary condition that removes capture orbits, resulting in an anisotropic distribution function. Finally, we evolve an initial spike distribution function by integrating the Boltzmann equation to show how the weak cusp grows and its density decreases with time. We treat two cases, one for $s$-wave and the other for $p$-wave DM annihilation, adopting parameters characteristic of the Milky Way nuclear core and typical WIMP models for DM. The cusp density profile for $p$-wave annihilation is weaker, varying like $\sim r^{-0.34}$, but is still not a flat plateau.

We perform general-relativistic hydrodynamical simulations of dynamical capture binary neutron star mergers, emphasizing the role played by the neutron star spin. Dynamical capture mergers may take place in globular clusters, as well as other dense stellar systems, where most neutron stars have large spins. We find significant variability in the merger outcome as a function of initial neutron star spin. For cases where the spin is aligned with the orbital angular momentum, the additional centrifugal support in the remnant hypermassive neutron star can prevent the prompt collapse to a black hole, while for antialigned cases the decreased total angular momentum can facilitate the collapse to a black hole. We show that even moderate spins can significantly increase the amount of ejected material, including the amount unbound with velocities greater than half the speed of light, leading to brighter electromagnetic signatures associated with kilonovae and interaction of the ejecta with the interstellar medium. Furthermore, we find that the initial neutron star spin can strongly affect the already rich phenomenology in the postmerger gravitational wave signatures that arise from the oscillation modes of the hypermassive neutron star. In several of our simulations, the resulting hypermassive neutron star develops the one-arm ($m=1$) spiral instability, the most pronounced cases being those with small but non-negligible neutron star spins. For long-lived hypermassive neutron stars, the presence of this instability leads to improved prospects for detecting these events through gravitational waves, and thus may give information about the neutron star equation of state.

Using general-relativistic hydrodynamical simulations, we show that merging binary neutron stars can form hypermassive neutrons stars that undergo the one-arm spiral instability. We study the particular case of a dynamical capture merger where the stars have a small spin, as may arise in globular clusters, and focus on an equal-mass scenario where the spins are aligned with the orbital angular momentum. We find that this instability develops when post-merger fluid vortices lead to the generation of a toroidal remnant - a configuration whose maximum density occurs in a ring around the center-of-mass - with high vorticity along its rotation axis. The instability quickly saturates on a timescale of $\sim 10$ ms, with the $m=1$ azimuthal density multipole mode dominating over higher modes. The instability also leaves a characteristic imprint on the post-merger gravitational wave signal that could be detectable if the instability persists in long-lived remnants.

In the extreme violence of merger and mass accretion, compact objects like black holes and neutron stars are thought to launch some of the most luminous outbursts of electromagnetic and gravitational wave energy in the Universe. Modeling these systems realistically is a central problem in theoretical astrophysics, but has proven extremely challenging, requiring the development of numerical relativity codes that solve Einstein's equations for the spacetime, coupled to the equations of general relativistic (ideal) magnetohydrodynamics (GRMHD) for the magnetized fluids. Over the past decade, the Illinois Numerical Relativity (ILNR) Group's dynamical spacetime GRMHD code has proven itself as a robust and reliable tool for theoretical modeling of such GRMHD phenomena. However, the code was written "by experts and for experts" of the code, with a steep learning curve that would severely hinder community adoption if it were open-sourced. Here we present IllinoisGRMHD, which is an open-source, highly-extensible rewrite of the original closed-source GRMHD code of the ILNR Group. Reducing the learning curve was the primary focus of this rewrite, with the goal of facilitating community involvement in the code's use and development, as well as the minimization of human effort in generating new science. IllinoisGRMHD also saves computer time, generating roundoff-precision identical output to the original code on adaptive-mesh grids, but nearly twice as fast at scales of hundreds to thousands of cores.

We present a new method to measure or constrain p-wave-suppressed cross sections for dark matter (DM) annihilations inside the steep density spikes induced by supermassive black holes. We demonstrate that the high DM densities, together with the increased velocity dispersion, within such spikes combine to make thermal p-wave annihilation cross-sections potentially visible in gamma-ray observations of the Galactic center (GC). The resulting DM signal is a bright central point source with emission originating from DM annihilations in the absence of a detectable spatially-extended signal from the halo. We define two simple reference theories of DM with a thermal p-wave annihilation cross-section and establish new limits on the combined particle and astrophysical parameter space of these models, demonstrating that Fermi is currently sensitive to thermal p-wave DM over a wide range of possible scenarios for the DM distribution in the GC.

We perform magnetohydrodynamic simulations in full general relativity (GRMHD) of a binary black hole-neutron star on a quasicircular orbit that undergoes merger. The binary mass ratio is 3:1, the black hole initial spin parameter $a/m=0.75$ ($m$ is the black hole Christodoulou mass) aligned with the orbital angular momentum, and the neutron star is an irrotational $\Gamma=2$ polytrope. About two orbits prior to merger (at time $t=t_B$), we seed the neutron star with a dynamically weak interior dipole magnetic field that extends into the stellar exterior. At $t=t_B$ the exterior has a low-density atmosphere with constant plasma parameter $\beta\equiv P_{\rm gas}/P_{\rm mag}$. Varying $\beta$ at $t_B$ in the exterior from $0.1$ to $0.01$, we find that at a time $\sim 4000{\rm M} \sim 100(M_{\rm NS}/1.4M_\odot)$ms [M is the total (ADM) mass] following the onset of accretion of tidally disrupted debris, magnetic winding above the remnant black hole poles builds up the magnetic field sufficiently to launch a mildly relativistic, collimated outflow - an incipient jet. The duration of the accretion and the lifetime of the jet is $\Delta t\sim 0.5(M_{\rm NS}/1.4M_\odot)$s. Our simulations furnish the first explicit examples in GRMHD which show that a jet can emerge following a black hole - neutron star merger.

We report results from simulations in general relativity of magnetized disks accreting onto merging black hole binaries, starting from relaxed disk initial data. The simulations feature an effective, rapid radiative cooling scheme as a limiting case of future treatments with radiative transfer. Here we evolve the systems after binary-disk decoupling through inspiral and merger, and analyze the dependence on the binary mass ratio with $q\equiv m_{\rm bh}/M_{\rm BH}=1,1/2,$ and $1/4$. We find that the luminosity associated with local cooling is larger than the luminosity associated with matter kinetic outflows, while the electromagnetic (Poynting) luminosity associated with bulk transport of magnetic field energy is the smallest. The cooling luminosity around merger is only marginally smaller than that of a single, non-spinning black hole. Incipient jets are launched independently of the mass ratio, while the same initial disk accreting on a single non-spinning black hole does not lead to a jet, as expected. For all mass ratios we see a transient behavior in the collimated, magnetized outflows lasting $2-5 ( M/10^8M_\odot ) \rm days$ after merger: the outflows become increasingly magnetically dominated and accelerated to higher velocities, boosting the Poynting luminosity. These sudden changes can alter the electromagnetic emission across the jet and potentially help distinguish mergers of black holes in AGNs from single accreting black holes based on jet morphology alone.

Robust gauge conditions are critically important to the stability and accuracy of numerical relativity (NR) simulations involving compact objects. Most of the NR community use the highly robust---though decade-old---moving-puncture (MP) gauge conditions for such simulations. It has been argued that in binary black hole (BBH) evolutions adopting this gauge, noise generated near adaptive-mesh-refinement (AMR) boundaries does not converge away cleanly with increasing resolution, severely limiting gravitational waveform accuracy at computationally feasible resolutions. We link this noise to a sharp (short-wavelength), initial outgoing gauge wave crossing into progressively lower resolution AMR grids, and present improvements to the standard MP gauge conditions that focus on stretching, smoothing, and more rapidly settling this outgoing wave. Our best gauge choice greatly reduces gravitational waveform noise during inspiral, yielding less fluctuation in convergence order and $\sim 40%$ lower waveform phase and amplitude errors at typical resolutions. Noise in other physical quantities of interest is also reduced, and constraint violations drop by more than an order of magnitude. We expect these improvements will carry over to simulations of all types of compact binary systems, as well as other $N$+1 formulations of gravity for which MP-like gauge conditions can be chosen.

If the supermassive black hole Sgr A* at the center of the Milky Way grew adiabatically from an initial seed embedded in an NFW dark matter (DM) halo, then the DM profile near the hole has steepened into a spike. We calculate the dramatic enhancement to the gamma ray flux from the Galactic center (GC) from such a spike if the 1-3 GeV excess observed in Fermi data is due to DM annihilations. We find that for the parameter values favored in recent fits, the point source-like flux from the spike is 35 times greater than the flux from the inner $1^\circ$ of the halo, far exceeding all Fermi point source detections near the GC. We consider the dependence of the spike signal on astrophysical and particle parameters and conclude that if the GC excess is due to DM, then a canonical adiabatic spike is disfavored by the data. We discuss alternative Galactic histories that predict different spike signals, including: (i) the nonadiabatic growth of the black hole, possibly associated with halo and/or black hole mergers, (ii) gravitational interaction of DM with baryons in the dense core, such as heating by stars, or (iii) DM self-interactions. We emphasize that the spike signal is sensitive to a different combination of particle parameters than the halo signal, and that the inclusion of a spike component to any DM signal in future analyses would provide novel information about both the history of the GC and the particle physics of DM annihilations.

Adopting our new method for matching general relativistic, ideal magnetohydrodynamics to its force-free limit, we perform the first systematic simulations of force-free pulsar magnetospheres in general relativity. We endow the neutron star with a general relativistic dipole magnetic field, model the interior with ideal magnetohydrodynamics, and adopt force-free electrodynamics in the exterior. Comparing the spin-down luminosity to its corresponding Minkowski value, we find that general relativistic effects give rise to a modest enhancement: the maximum enhancement for $n=1$ polytropes is $\sim 23\%$. Evolving a rapidly rotating $n=0.5$ polytrope we find an even greater enhancement of $\sim 35\%$. Using our simulation data, we derive fitting formulas for the pulsar spin-down luminosity as a function of the neutron star compaction, angular speed, and dipole magnetic moment. We expect stiffer equations of state and more rapidly spinning neutron stars to lead to even larger enhancements in the spin-down luminosity.

We adopt the conduction fluid approximation to model the steady-state distribution of matter around a massive black hole at the center of a weakly collisional cluster of particles. By "`weakly collisional" we mean a cluster in which the mean free time between particle collisions is much longer than the characteristic particle crossing (dynamical) time scale, but shorter than the cluster lifetime. When applied to a star cluster, we reproduce the familiar Bahcall-Wolf power-law cusp solution for the stars bound to the black hole. Here the star density scales with radius as $r^{-7/4}$ and the velocity dispersion as $r^{-1/2}$ throughout most of the gravitational well of the black hole. When applied to a relaxed, self-interacting dark matter (SIDM) halo with a velocity-dependent cross section $\sigma \sim v^{-a}$, the gas again forms a power-law cusp, but now the SIDM density scales as $r^{-\beta}$, where $\beta = (a+3)/4$, while its velocity dispersion again varies as $r^{-1/2}$. Results are obtained first in Newtonian theory and then in full general relativity. Although the conduction fluid model is a simplification, it provides a reasonable first approximation to the matter profiles and is much easier to implement than a full Fokker-Planck treatment or an $N$-body simulation of the Boltzmann equation with collisional perturbations.

The Numerical INJection Analysis (NINJA) project is a collaborative effort between members of the numerical relativity and gravitational-wave astrophysics communities. The purpose of NINJA is to study the ability to detect gravitational waves emitted from merging binary black holes and recover their parameters with next-generation gravitational-wave observatories. We report here on the results of the second NINJA project, NINJA-2, which employs 60 complete binary black hole hybrid waveforms consisting of a numerical portion modelling the late inspiral, merger, and ringdown stitched to a post-Newtonian portion modelling the early inspiral. In a "blind injection challenge" similar to that conducted in recent LIGO and Virgo science runs, we added 7 hybrid waveforms to two months of data recolored to predictions of Advanced LIGO and Advanced Virgo sensitivity curves during their first observing runs. The resulting data was analyzed by gravitational-wave detection algorithms and 6 of the waveforms were recovered with false alarm rates smaller than 1 in a thousand years. Parameter estimation algorithms were run on each of these waveforms to explore the ability to constrain the masses, component angular momenta and sky position of these waveforms. We also perform a large-scale monte-carlo study to assess the ability to recover each of the 60 hybrid waveforms with early Advanced LIGO and Advanced Virgo sensitivity curves. Our results predict that early Advanced LIGO and Advanced Virgo will have a volume-weighted average sensitive distance of 300Mpc (1Gpc) for $10M_{\odot}+10M_{\odot}$ ($50M_{\odot}+50M_{\odot}$) binary black hole coalescences. We demonstrate that neglecting the component angular momenta in the waveform models used in matched-filtering will result in a reduction in sensitivity for systems with large component angular momenta. [Abstract abridged for ArXiv, full version in PDF]

We report on simulations in general relativity of magnetized disks onto black hole binaries. We vary the binary mass ratio from 1:1 to 1:10 and evolve the systems when they orbit near the binary-disk decoupling radius. We compare (surface) density profiles, accretion rates (relative to a single, non-spinning black hole), variability, effective $\alpha$-stress levels and luminosities as functions of the mass ratio. We treat the disks in two limiting regimes: rapid radiative cooling and no radiative cooling. The magnetic field lines clearly reveal jets emerging from both black hole horizons and merging into one common jet at large distances. The magnetic fields give rise to much stronger shock heating than the pure hydrodynamic flows, completely alter the disk structure, and boost accretion rates and luminosities. Accretion streams near the horizons are among the densest structures; in fact, the 1:10 no-cooling evolution results in a refilling of the cavity. The typical effective temperature in the bulk of the disk is $\sim 10^5 (M/10^8 M_\odot)^{-1/4} (L/L_{\rm edd})^{1/4} {\rm K}$ yielding characteristic thermal frequencies $\sim 10^{15} (M/10^8 M_\odot)^{-1/4} (L/L_{\rm edd})^{1/4}(1+z)^{-1}{\rm Hz} $. These systems are thus promising targets for many extragalactic optical surveys, such as LSST, WFIRST, and PanSTARRS.

We present a new computational method for smoothly matching general relativistic ideal magnetohydrodynamics (MHD) to its force-free limit. The method is based on a flux-conservative formalism for MHD and its force-free limit, and a vector potential formulation for the induction equation to maintain the zero divergence constraint for the magnetic field. The force-free formulation we adopt evolves the magnetic field and the Poynting vector, instead of the magnetic and electric fields. We show that our force-free code passes a robust suite of tests, performed both in 1D flat spacetime and in 3D black hole spacetimes. We also demonstrate that our matching technique successfully reproduces the aligned rotator force-free solution. Our new techniques are suitable for studying electromagnetic effects and predicting electromagnetic signals arising in many different curved spacetime scenarios. For example, we can treat spinning neutron stars, either in isolation or in compact binaries, that have MHD interiors and force-free magnetospheres.

We perform the first general relativistic force-free simulations of neutron star (NS) magnetospheres in orbit about spinning and non-spinning black holes. We find promising precursor electromagnetic emission: typical Poynting luminosities at, e.g., an orbital separation of 6.6 times the NS radius are L ~ 6 x 10^{42} erg/s for a 1.4 solar-mass NS with a 10^{13}G polar magnetic field. The Poynting flux peaks within a broad beam of ~40 degrees in the azimuthal direction and within ~60 degrees from the orbital plane, establishing a possible lighthouse effect. Our calculations, though preliminary, preview more detailed simulations of these systems that we plan to perform in the future.

We consider a geometrically thin, Keplerian disk in the orbital plane of a binary black hole (BHBH) consisting of a spinning primary and low-mass secondary (mass ratio q < 1). To account for the principle effects of general relativity (GR), we propose a modification of the standard Newtonian evolution equation for the (orbit-averaged) time-varying disk surface density. In our modified equation the viscous torque in the disk is treated in full GR, while the tidal torque is handled in the Newtonian limit. Our GR-hybrid treatment is reasonable because the tidal torque is concentrated near the orbital radius of the secondary and is most important prior to binary-disk decoupling, when the orbital separation is large and resides in the weak-field regime. The tidal torque on the disk diminishes during late merger and vanishes altogether following merger. By contrast, the viscous torque drives the flow into the strong-field region and onto the primary during all epochs. Following binary coalescence, the viscous torque alone governs the time-dependent accretion onto the remnant, as well as the temporal behavior, strength and spectrum of the aftermath electromagnetic radiation from the disk. We solve our GR-hybrid equation for a representative BHBH-disk system, identify several observable EM signatures of the merger, and compare results obtained for the gas and EM radiation with those found with the Newtonian prescription.

TwoPunctures is perhaps the most widely-adopted code for generating binary black hole "puncture" initial data and interpolating these (spectral) data onto evolution grids. In typical usage, the bulk of this code's run time is spent in its spectral interpolation routine. We announce a new publicly-available spectral interpolation routine that improves the performance of the original interpolation routine by a factor of ~100, yielding results consistent with the original spectral interpolation routine to roundoff precision. This note serves as a guide for installing this routine both in the original standalone TwoPunctures code and the Einstein Toolkit supported version of this code.

We present new numerical techniques we developed for launching the first parameter study of {\it magnetized} black hole--neutron star (BHNS) mergers, varying the magnetic fields seeded in the initial neutron star. We found that magnetic fields have a negligible impact on the gravitational waveforms and bulk dynamics of the system during merger, regardless of magnetic field strength or BH spin. In a recent simulation, we seeded the remnant disk from an unmagnetized BHNS merger simulation with large-scale, purely poloidal magnetic fields, which are otherwise absent in the full simulation. The outcome appears to be a viable sGRB central engine.

We present results from the first fully general relativistic, magnetohydrodynamic (GRMHD) simulations of an equal-mass black hole binary (BHBH) in a magnetized, circumbinary accretion disk. We simulate both the pre and post-decoupling phases of a BHBH-disk system and both "cooling" and "no-cooling" gas flows. Prior to decoupling, the competition between the binary tidal torques and the effective viscous torques due to MHD turbulence depletes the disk interior to the binary orbit. However, it also induces a two-stream accretion flow and mildly relativistic polar outflows from the BHs. Following decoupling, but before gas fills the low-density "hollow" surrounding the remnant, the accretion rate is reduced, while there is a prompt electromagnetic (EM) luminosity enhancement following merger due to shock heating and accretion onto the spinning BH remnant. This investigation, though preliminary, previews more detailed GRMHD simulations we plan to perform in anticipation of future, simultaneous detections of gravitational and EM radiation from a merging BHBH-disk system.

Black hole--neutron star (BHNS) binary mergers can form disks in which magnetorotational instability (MRI)-induced turbulence may drive accretion onto the remnant BH, supporting relativistic jets and providing the engine for a short-hard gamma-ray burst (SGRB). Our earlier study of magnetized BHNSs showed that NS tidal disruption winds the magnetic field into a toroidal configuration, with poloidal fields so weak that capturing MRI with full-disk simulations would require $\sim 10^8$ CPU-hours. In that study we imposed equatorial symmetry, suppressing poloidal magnetic fields that might be generated from plasma crossing the orbital plane. Here we show that initial conditions that break this symmetry (i.e., {\it tilted} poloidal magnetic fields in the NS) generate much stronger poloidal fields in the disk, indicating that asymmetric initial conditions may be necessary for establishing BHNS mergers as SGRB progenitors via fully general relativistic MHD simulations. We demonstrate that BHNS mergers may form an SGRB engine under the right conditions by seeding the remnant disk from an unmagnetized BHNS simulation with purely poloidal fields dynamically unimportant initially, but strong enough to resolve MRI. Magnetic turbulence occurs in the disk, driving accretion and supporting Poynting-dominated jet outflows sufficient to power an SGRB.

The inspiral and merger of a binary neutron star (NSNS) can lead to the formation of a hypermassive neutron star (HMNS). As the HMNS loses thermal pressure due to neutrino cooling and/or centrifugal support due to gravitational wave (GW) emission, and/or magnetic breaking of differential rotation it will collapse to a black hole. To assess the importance of shock-induced thermal pressure and cooling, we adopt an idealized equation of state and perform NSNS simulations in full GR through late inspiral, merger, and HMNS formation, accounting for cooling. We show that thermal pressure contributes significantly to the support of the HMNS against collapse and that thermal cooling accelerates its "delayed" collapse. Our simulations demonstrate explicitly that cooling can induce the catastrophic collapse of a hot hypermassive neutron star formed following the merger of binary neutron stars. Thus, cooling physics is important to include in NSNS merger calculations to accurately determine the lifetime of the HMNS remnant and to extract information about the NS equation of state, cooling mechanisms, bar instabilities and B-fields from the GWs emitted during the transient phase prior to BH formation.

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.

Tidal torques acting on a gaseous accretion disk around a binary black hole can create a gap in the disk near the orbital radius. At late times, when the binary inspiral timescale due to gravitational wave emission becomes shorter than the viscous timescale in the disk, the binary decouples from the disk and eventually merges. Prior to decoupling the balance between tidal and viscous torques drives the disk to a quasistationary equilibrium state, perturbed slightly by small amplitude, spiral density waves emanating from the edges of the gap. We consider a black hole binary with a companion of smaller mass and construct a simple Newtonian model for a geometrically thin, Keplerian disk in the orbital plane of the binary. We solve the disk evolution equations in steady state to determine the quasistationary, (orbit-averaged) surface density profile prior to decoupling. We use our solution, which is analytic up to simple quadratures, to compute the electromagnetic flux and approximate radiation spectrum during this epoch. A single nondimensional parameter Td/Tvis, equal to the ratio of the tidal to viscous torque at the orbital radius, determines the disk structure, including the surface density profile, the extent of the gap, the existence of an inner disk, and the accretion rate. The solution reduces to the Shakura-Sunyaev profile for a stationary accretion disk around a single black hole in the limit of small Td/Tvis. Our solution may be useful for choosing physical parameters and setting up quasistationary disk initial data for detailed numerical simulations that begin prior to decoupling and track the subsequent evolution of a black hole binary-disk system.

As a neutron star (NS) is tidally disrupted by a black hole (BH) companion at the end of a BH-NS binary inspiral, its magnetic fields will be stretched and amplified. If sufficiently strong, these magnetic fields may impact the gravitational waveforms, merger evolution and mass of the remnant disk. Formation of highly-collimated magnetic field lines in the disk+spinning BH remnant may launch relativistic jets, providing the engine for a short-hard GRB. We analyze this scenario through fully general relativistic, magnetohydrodynamic (GRMHD) BHNS simulations from inspiral through merger and disk formation. Different initial magnetic field configurations and strengths are chosen for the NS interior for both nonspinning and moderately spinning (a/M=0.75) BHs aligned with the orbital angular momentum. Only strong interior (Bmax~10^17 G) initial magnetic fields in the NS significantly influence merger dynamics, enhancing the remnant disk mass by 100% and 40% in the nonspinning and spinning BH cases, respectively. However, detecting the imprint of even a strong magnetic field may be challenging for Advanced LIGO. Though there is no evidence of mass outflows or magnetic field collimation during the preliminary simulations we have performed, higher resolution, coupled with longer disk evolutions and different initial magnetic field configurations, may be required to definitively assess the possibility of BHNS binaries as short-hard GRB progenitors.

We present fully general relativistic (GR) simulations of binary white dwarf-neutron star (WDNS) inspiral and merger. The initial binary is in a circular orbit at the Roche critical separation. The goal is to determine the ultimate fate of such systems. We focus on binaries whose total mass exceeds the maximum mass (Mmax) a cold, degenerate EOS can support against gravitational collapse. The time and length scales span many orders of magnitude, making fully general relativistic hydrodynamic (GRHD) simulations computationally prohibitive. For this reason, we model the WD as a "pseudo-white dwarf" (pWD) as in our binary WDNS head-on collisions study [PRD83:064002,2011]. Our GRHD simulations of a pWDNS system with a 0.98-solar-mass WD and a 1.4-solar-mass NS show that the merger remnant is a spinning Thorne-Zytkow-like Object (TZlO) surrounded by a massive disk. The final total rest mass exceeds Mmax, but the remnant does not collapse promptly. To assess whether the object will ultimately collapse after cooling, we introduce radiative thermal cooling. We first apply our cooling algorithm to TZlOs formed in WDNS head-on collisions, and show that these objects collapse and form black holes on the cooling time scale, as expected. However, when we cool the spinning TZlO formed in the merger of a circular-orbit WDNS binary, the remnant does not collapse, demonstrating that differential rotational support is sufficient to prevent collapse. Given that the final total mass exceeds Mmax, magnetic fields and/or viscosity may redistribute angular momentum and ultimately lead to delayed collapse to a BH. We infer that the merger of realistic massive WDNS binaries likely will lead to the formation of spinning TZlOs that undergo delayed collapse.

We recently developed a new general relativistic magnetohydrodynamic code with adaptive mesh refinement that evolves the electromagnetic (EM) vector potential (A) instead of the magnetic fields directly. Evolving A enables one to use any interpolation scheme on refinement level boundaries and still guarantee that the magnetic field remains divergenceless. As in classical EM, a gauge choice must be made when evolving A, and we chose a straightforward "algebraic" gauge condition to simplify the A evolution equation. However, magnetized black hole-neutron star (BHNS) simulations in this gauge exhibit unphysical behavior, including the spurious appearance of strong magnetic fields on refinement level boundaries. This spurious behavior is exacerbated when matter crosses refinement boundaries during tidal disruption of the NS. Applying Kreiss-Oliger dissipation to the evolution of the magnetic vector potential A slightly weakens this spurious magnetic effect, but with undesired consequences. We demonstrate via an eigenvalue analysis and a numerical study that zero-speed modes in the algebraic gauge, coupled with the frequency filtering that occurs on refinement level boundaries, are responsible for the creation of spurious magnetic fields. We show that the EM Lorenz gauge exhibits no zero-speed modes, and as a consequence, spurious magnetic effects are quickly propagated away, allowing for long-term, stable magnetized BHNS evolutions. Our study demonstrates how the EM gauge degree of freedom can be chosen to one's advantage, and that for magnetized BHNS simulations the Lorenz gauge constitutes a major improvement over the algebraic gauge.

Simultaneous gravitational and electromagnetic wave observations of merging black hole binaries (BHBHs) can provide unique opportunities to study gravitation physics, accretion and cosmology. Here we perform fully general relativistic, hydrodynamic simulations of equal-mass, nonspinning BHBHs coalescing in a circumbinary disk. We evolve the metric using the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of Einstein's field equations with standard moving puncture gauge conditions. We handle the hydrodynamics via a high-resolution shock-capturing (HRSC) scheme. We track the inspiral starting from a binary separation of 10M, where M is the total binary mass. We take the disks to have an inner radius at R_in~15M to account for the hollow created by the binary torques. Our disks extend to R=65M and have an initial scale height of H/R=0.03-0.11. The gas is governed by a Gamma-law EOS, with Gamma equal to 5/3, 4/3, and 1.1. Disks are allowed to relax in the "early inspiral" epoch to provide quasistationary realistic initial data. We then evolve the metric and matter during the "late inspiral and merger" epoch. The later simulations are designed to track BHBH inspiral following disk-binary decoupling, through merger and ringdown, terminating before viscosity has time to fill the hollow about the remnant. We compute the gas flow and accretion rate and estimate the electromagnetic luminosity due to bremsstrahlung and synchrotron emission as a perturbation for optically thin disks. The synchrotron component of the luminosity peaks in the infrared band and should be detectable by WFIRST and possibly the LSST for a 10^8 M_sun binary embedded in a disk with a density n~10^12/cm^3 at z=1, beginning with a maximum value of $L~10^46 n_12^2 M_8^3 erg/s at decoupling, and decreasing steadily over a timescale of ~100 M_8 hours to a value of L~10^45 n_12^2 M_8^3 erg/s at merger.

Theories of gravity other than general relativity (GR) can explain the observed cosmic acceleration without a cosmological constant. One such class of theories of gravity is f(R). Metric f(R) theories have been proven to be equivalent to Brans-Dicke (BD) scalar-tensor gravity without a kinetic term. Using this equivalence and a 3+1 decomposition of the theory it has been shown that metric f(R) gravity admits a well-posed initial value problem. However, it has not been proven that the 3+1 evolution equations of metric f(R) gravity preserve the (hamiltonian and momentum) constraints. In this paper we show that this is indeed the case. In addition, we show that the mathematical form of the constraint propagation equations in BD-equilavent f(R) gravity and in f(R) gravity in both the Jordan and Einstein frames, is exactly the same as in the standard ADM 3+1 decomposition of GR. Finally, we point out that current numerical relativity codes can incorporate the 3+1 evolution equations of metric f(R) gravity by modifying the stress-energy tensor and adding an additional scalar field evolution equation. We hope that this work will serve as a starting point for relativists to develop fully dynamical codes for valid f(R) models.

We simulate head-on collisions from rest at large separation of binary white dwarf -- neutron stars (WDNSs) in full general relativity. Our study serves as a prelude to our analysis of the circular binary WDNS problem. We focus on compact binaries whose total mass exceeds the maximum mass that a cold degenerate star can support, and our goal is to determine the fate of such systems. A fully general relativistic hydrodynamic computation of a realistic WDNS head-on collision is prohibitive due to the large range of dynamical time scales and length scales involved. For this reason, we construct an equation of state (EOS) which captures the main physical features of NSs while, at the same time, scales down the size of WDs. We call these scaled-down WD models "pseudo-WDs (pWDs)". Using pWDs, we can study these systems via a sequence of simulations where the size of the pWD gradually increases toward the realistic case. We perform two sets of simulations; One set studies the effects of the NS mass on the final outcome, when the pWD is kept fixed. The other set studies the effect of the pWD compaction on the final outcome, when the pWD mass and the NS are kept fixed. All simulations show that 14%-18% of the initial total rest mass escapes to infinity. All remnant masses still exceed the maximum rest mass that our cold EOS can support (1.92 solar masses), but no case leads to prompt collapse to a black hole. This outcome arises because the final configurations are hot. All cases settle into spherical, quasiequilibrium configurations consisting of a cold NS core surrounded by a hot mantle, resembling Thorne-Zytkow objects. Extrapolating our results to realistic WD compactions, we predict that the likely outcome of a head-on collision of a realistic, massive WDNS system will be the formation of a quasiequilibrium Thorne-Zytkow-like object.

We have written and tested a new general relativistic magnetohydrodynamics (GRMHD) code, capable of evolving MHD fluids in dynamical spacetimes with adaptive-mesh refinement (AMR). Our code solves the Einstein-Maxwell-MHD system of coupled equations in full 3+1 dimensions, evolving the metric via the Baumgarte-Shapiro Shibata-Nakamura (BSSN) formalism and the MHD and magnetic induction equations via a conservative, high-resolution shock-capturing scheme. The induction equations are recast as an evolution equation for the magnetic vector potential, which exists on a grid that is staggered with respect to the hydrodynamic and metric variables. The divergenceless constraint div(B)=0 is enforced by the curl of the vector potential. Our MHD scheme is fully compatible with AMR, so that fluids at AMR refinement boundaries maintain div(B)=0. In simulations with uniform grid spacing, our MHD scheme is numerically equivalent to a commonly used, staggered-mesh constrained-transport scheme. We present code validation test results, both in Minkowski and curved spacetimes. They include magnetized shocks, nonlinear Alfv\'en waves, cylindrical explosions, cylindrical rotating disks, magnetized Bondi tests, and the collapse of a magnetized rotating star. Some of the more stringent tests involve black holes. We find good agreement between analytic and numerical solutions in these tests, and achieve convergence at the expected order.

We propose a new radial coordinate to write the Kerr metric in puncture form. Unlike the quasi-radial coordinate introduced previously, the horizon radius remains finite in our radial coordinate in the extreme Kerr limit a/M -> 1. This significantly improves the accuracy of the evolution of black holes with spins close to the extreme Kerr limit. We are able to evolve accurately both stationary and boosted black holes with spins as high as a/M=0.99 using initial data constructed in these new puncture coordinates. Initial data of compact binaries with rapidly spinning black holes can be constructed using our proposed new puncture metric for the background conformal metric. Our simulations for single black holes suggest that such initial data can be evolved successfully by the moving puncture technique.

Merging supermassive black hole-black hole (BHBH) binaries produced in galaxy mergers are promising sources of detectable gravitational waves. If such a merger takes place in a gaseous environment, there is a possibility of a simultaneous detection of electromagnetic and gravitational radiation, as the stirring, shock heating and accretion of the gas may produce variability and enhancements in the electromagnetic flux. Such a simultaneous detection can provide a wealth of opportunities to study gravitational physics, accretion physics, and cosmology. We investigate this scenario by performing fully general relativistic, hydrodynamic simulations of merging, equal-mass, nonspinning BHBH binaries embedded in gas clouds. We evolve the metric using the BSSN formulation with standard moving puncture gauge conditions and handle the hydrodynamics via a high-resolution shock-capturing (HRSC) scheme. We consider both "binary Bondi accretion" in which the binary is at rest relative to the ambient gas cloud, as well as "binary Bondi-Hoyle-Lyttleton accretion" in which the binary moves relative to the gas cloud. The gas cloud is assumed to be homogeneous far from the binary and governed by a \Gamma-law equation of state. We vary \Gamma between 4/3 and 5/3. For each simulation, we compute the gas flow and accretion rate and estimate the electromagnetic luminosity due to bremsstrahlung and synchrotron emission. We find evidence for significant enhancements in both the accretion rate and luminosity over values for a single black hole of the same mass as the binary. We estimate that this luminosity enhancement should be detectable by LSST for a 10^6 M_sun binary in a hot gas cloud of density n~10/cm^3 and temperature T~10^6 K at z=1, reaching a maximum of L~3x10^43 erg/s, with the emission peaking in the visible band.

Tidal torques from a binary black hole (BHBH) empty out the central regions in any circumbinary gaseous accretion disk. The balance between tidal torques and viscosity maintain the inner edge of the disk at a radius r ~ 1.5a -- 2a, where a is the binary semimajor axis. Eventually, the inspiraling binary decouples from disk and merges, leaving behind a central hollow ("donut hole") in the disk orbiting the remnant black hole. We present a simple, time-dependent, Newtonian calculation that follows the secular (viscous) evolution of the disk as it fills up the hollow down to the black hole innermost stable circular orbit and then relaxes to stationary equilibrium. We use our model to calculate the electromagnetic radiation ("afterglow") spectrum emitted during this transient accretion epoch. Observing the temporal increase in the total electromagnetic flux and the hardening of the spectrum as the donut hole fills may help confirm a BHBH merger detected by a gravitational wave interferometer. We show how the very existence of the initial hollow can lead to super-Eddington accretion during this secular phase if the rate is not very far below Eddington prior to decoupling. Our model, though highly idealized, may be useful in establishing some of the key parameters, thermal emission features and scalings that characterize this transient. It can serve as a guide in the design and calibration of future radiation-magnetohydrodynamic simulations in general relativity.

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.

Binary neutron stars (NSNS) are expected to be among the leading sources of gravitational waves observable by ground-based laser interferometers and may be the progenitors of short-hard gamma ray bursts. We present a series of general relativistic NSNS coalescence simulations both for unmagnetized and magnetized stars. We adopt quasiequilibrium initial data for circular, irrotational binaries constructed in the conformal thin-sandwich (CTS) framework. We adopt the BSSN formulation for evolving the metric and a high-resolution shock-capturing scheme to handle the magnetohydrodynamics. Our simulations of unmagnetized binaries confirm the results of Shibata, Taniguchi and Uryu (2003). In cases in which the mergers result in a prompt collapse to a black hole, we are able to use puncture gauge conditions to extend the evolution and determine the mass of the material that forms a disk. We find that the disk mass is less than 2% of the total mass in all cases studied. We then add a small poloidal magnetic field to the initial configurations and study the subsequent evolution. For cases in which the remnant is a hypermassive neutron star, we see measurable differences in both the amplitude and phase of the gravitational waveforms following the merger. For cases in which the remnant is a black hole surrounded by a disk, the disk mass and the gravitational waveforms are about the same as the unmagnetized cases. Magnetic fields substantially affect the long-term, secular evolution of a hypermassive neutron star (driving `delayed collapse') and an accretion disk around a nascent black hole.

Many systems of current interest in relativistic astrophysics require a knowledge of radiative transfer in a magnetized gas flowing in a strongly-curved, dynamical spacetime. Such systems include coalescing compact binaries containing neutron stars or white dwarfs, disks around merging black holes, core collapse supernovae, collapsars, and gamma-ray burst sources. To model these phenomena, all of which involve general relativity, radiation (photon and/or neutrino), and magnetohydrodynamics, we have developed a general relativistic code capable of evolving MHD fluids and radiation in dynamical spacetimes. Our code solves the coupled Einstein-Maxwell-MHD-Radiation system of equations both in axisymmetry and in full 3+1 dimensions. We evolve the metric by integrating the BSSN equations, and use a conservative, high-resolution shock-capturing scheme to evolve both the MHD and radiation moment equations. In this paper, we implement our scheme for optically thick gases and grey-body opacities. Our code gives accurate results in a suite of tests involving radiating shocks and nonlinear waves propagating in Minkowski spacetime. In addition, to test our code's ability to evolve the relativistic radiation-MHD equations in strong-field dynamical spacetimes, we study "thermal Oppenheimer-Snyder collapse" to a black hole, and find good agreement between analytic and numerical solutions.

We study the evolution in axisymmetry of accretion disks formed self-consistently through collapse of magnetized hypermassive neutron stars to black holes. Such stars can arise following the merger of binary neutron stars. They are differentially rotating, dynamically stable, and have rest masses exceeding the mass limit for uniform rotation. However, hypermassive neutron stars are secularly unstable to collapse due to MHD-driven angular momentum transport. The rotating black hole which forms in this process is surrounded by a hot, massive, magnetized torus and a magnetic field collimated along the spin axis. This system is a candidate for the central engine of a short-hard gamma-ray burst (GRB). Our code integrates the coupled Einstein-Maxwell-MHD equations and is used to follow the collapse of magnetized hypermassive neutron star models in full general relativity until the spacetime settles down to a quasi-stationary state. We then employ the Cowling approximation, in which the spacetime is frozen, to track the subsequent evolution of the disk. This approximation allows us to greatly extend the disk evolutions and study the resulting outflows, which may be relevant to the generation of a GRB. We find that outflows are suppressed when a stiff equation of state is assumed for low density disk material and are sensitive to the initial magnetic field configuration.

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.

Black holes are popping up all over the place: in compact binary X-ray sources and GRBs, in quasars, AGNs and the cores of all bulge galaxies, in binary black holes and binary black hole-neutron stars, and maybe even in the LHC! Black holes are strong-field objects governed by Einstein's equations of general relativity. Hence general relativistic, numerical simulations of dynamical phenomena involving black holes may help reveal ways in which black holes can form, grow and be detected in the universe. To convey the state-of-the art, we summarize several representative simulations here, including the collapse of a hypermassive neutron star to a black hole following the merger of a binary neutron star, the magnetorotational collapse of a massive star to a black hole, and the formation and growth of supermassive black hole seeds by relativistic MHD accretion in the early universe.

We perform axisymmetric simulations of the magnetorotational collapse of very massive stars in full general relativity. Our simulations are applicable to the collapse of supermassive stars (M > 10^3M_sun) and to very massive Pop III stars. We model our initial configurations by n=3 polytropes. The ratio of magnetic to rotational kinetic energy in these configurations is chosen to be small (1% and 10%). We find that such magnetic fields do not affect the initial collapse significantly. The core collapses to a black hole, after which black hole excision is employed to continue the evolution long enough for the hole to reach a quasi-stationary state. We find that the black hole mass is M_h = 0.95M and its spin parameter is J_h/M_h^2 = 0.7, with the remaining matter forming a torus around the black hole. We freeze the spacetime metric ("Cowling approximation") and continue to follow the evolution of the torus after the black hole has relaxed to quasi-stationary equilibrium. In the absence of magnetic fields, the torus settles down following ejection of a small amount of matter due to shock heating. When magnetic fields are present, the field lines gradually collimate along the hole's rotation axis. MHD shocks and the MRI generate MHD turbulence in the torus and stochastic accretion onto the central black hole. When the magnetic field is strong, a wind is generated in the torus, and the torus undergoes radial oscillations that drive episodic accretion onto the hole. These oscillations produce long-wavelength gravitational waves potentially detectable by LISA. The final state of the magnetorotational collapse always consists of a central black hole surrounded by a collimated magnetic field and a hot, thick accretion torus. This system is a viable candidate for the central engine of a long-soft gamma-ray burst.

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.

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.

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.

Hypermassive neutron stars (HMNSs) are equilibrium configurations supported against collapse by rapid differential rotation and likely form as transient remnants of binary neutron star mergers. Though HMNSs are dynamically stable, secular effects such as viscosity or magnetic fields tend to bring HMNSs into uniform rotation and thus lead to collapse. We simulate the evolution of magnetized HMNSs in axisymmetry using codes which solve the Einstein-Maxwell-MHD system of equations. We find that magnetic braking and the magnetorotational instability (MRI) both contribute to the eventual collapse of HMNSs to rotating black holes surrounded by massive, hot accretion tori and collimated magnetic fields. Such hot tori radiate strongly in neutrinos, and the resulting neutrino-antineutrino annihilation could power short-hard GRBs.

The merger of two neutron stars usually produces a remnant with a mass significantly above the single (nonrotating) neutron star maximum mass. In some cases, the remnant will be stabilized against collapse by rapid, differential rotation. MHD-driven angular momentum transport eventually leads to the collapse of the remnant's core, resulting in a black hole surrounded by a massive accretion torus. Here we present simulations of this process. The plausibility of generating short duration gamma ray bursts through this scenario is discussed.

We study magnetohydrodynamic (MHD) effects arising in the collapse of magnetized, rotating, massive stellar cores to proto-neutron stars (PNSs). We perform axisymmetric numerical simulations in full general relativity with a hybrid equation of state. The formation and early evolution of a PNS are followed with a grid of 2500 x 2500 zones, which provides better resolution than in previous (Newtonian) studies. We confirm that significant differential rotation results even when the rotation of the progenitor is initially uniform. Consequently, the magnetic field is amplified both by magnetic winding and the magnetorotational instability (MRI). Even if the magnetic energy E_EM is much smaller than the rotational kinetic energy T_rot at the time of PNS formation, the ratio E_EM/T_rot increases to 0.1-0.2 by the magnetic winding. Following PNS formation, MHD outflows lead to losses of rest mass, energy, and angular momentum from the system. The earliest outflow is produced primarily by the increasing magnetic stress caused by magnetic winding. The MRI amplifies the poloidal field and increases the magnetic stress, causing further angular momentum transport and helping to drive the outflow. After the magnetic field saturates, a nearly stationary, collimated magnetic field forms near the rotation axis and a Blandford-Payne type outflow develops along the field lines. These outflows remove angular momentum from the PNS at a rate given by \dot{J} \sim \eta E_EM C_B, where \eta is a constant of order 0.1 and C_B is a typical ratio of poloidal to toroidal field strength. As a result, the rotation period quickly increases for a strongly magnetized PNS until the degree of differential rotation decreases. Our simulations suggest that rapidly rotating, magnetized PNSs may not give rise to rapidly rotating neutron stars.

The capacity to model magnetohydrodynamical (MHD) flows in dynamical, strongly curved spacetimes significantly extends the reach of numerical relativity in addressing many problems at the forefront of theoretical astrophysics. We have developed and tested an evolution code for the coupled Einstein-Maxwell-MHD equations which combines a BSSN solver with a high resolution shock capturing scheme. As one application, we evolve magnetized, differentially rotating neutron stars under the influence of a small seed magnetic field. Of particular significance is the behavior found for hypermassive neutron stars (HMNSs), which have rest masses greater the mass limit allowed by uniform rotation for a given equation of state. The remnant of a binary neutron star merger is likely to be a HMNS. We find that magnetic braking and the magnetorotational instability lead to the collapse of HMNSs and the formation of rotating black holes surrounded by massive, hot accretion tori and collimated magnetic field lines. Such tori radiate strongly in neutrinos, and the resulting neutrino-antineutrino annihilation (possibly in concert with energy extraction by MHD effects) could provide enough energy to power short-hard gamma-ray bursts. To explore the range of outcomes, we also evolve differentially rotating neutron stars with lower masses and angular momenta than the HMNS models. Instead of collapsing, the non-hypermassive models form nearly uniformly rotating central objects which, in cases with significant angular momentum, are surrounded by massive tori.

We present long-term (~10^4 M) axisymmetric simulations of differentially rotating, magnetized neutron stars in the slow-rotation, weak magnetic field limit using a perturbative metric evolution technique. Although this approach yields results comparable to those obtained via nonperturbative (BSSN) evolution techniques, simulations performed with the perturbative metric solver require about 1/4 the computational resources at a given resolution. This computational efficiency enables us to observe and analyze the effects of magnetic braking and the magnetorotational instability (MRI) at very high resolution. Our simulations demonstrate that (1) MRI is not observed unless the fastest-growing mode wavelength is resolved by more than about 10 gridpoints; (2) as resolution is improved, the MRI growth rate converges, but due to the small-scale turbulent nature of MRI, the maximum growth amplitude increases, but does not exhibit convergence, even at the highest resolution; and (3) independent of resolution, magnetic braking drives the star toward uniform rotation as energy is sapped from differential rotation by winding magnetic fields.

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.

We study the effects of magnetic fields on the evolution of differentially rotating neutron stars, which can form in stellar core collapse or binary neutron star coalescence. Magnetic braking and the magnetorotational instability (MRI) both redistribute angular momentum; the outcome of the evolution depends on the star's mass and spin. Simulations are carried out in axisymmetry using our recently developed codes which integrate the coupled Einstein-Maxwell-MHD equations. For initial data, we consider three categories of differentially rotating, equilibrium configurations, which we label normal, hypermassive and ultraspinning. Hypermassive stars have rest masses exceeding the mass limit for uniform rotation. Ultraspinning stars are not hypermassive, but have angular momentum exceeding the maximum for uniform rotation at the same rest mass. We show that a normal star will evolve to a uniformly rotating equilibrium configuration. An ultraspinning star evolves to an equilibrium state consisting of a nearly uniformly rotating central core, surrounded by a differentially rotating torus with constant angular velocity along magnetic field lines, so that differential rotation ceases to wind the magnetic field. In addition, the final state is stable against the MRI, although it has differential rotation. For a hypermassive neutron star, the MHD-driven angular momentum transport leads to catastrophic collapse of the core. The resulting rotating black hole is surrounded by a hot, massive, magnetized torus undergoing quasistationary accretion, and a magnetic field collimated along the spin axis--a promising candidate for the central engine of a short gamma-ray burst. (Abridged)

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.

(abridged) We describe highlights of the results of two observing campaigns in 2004 to investigate the correlation of flare activity in Sgr A* in different wavelength regimes, using a total of nine ground and space-based telescopes. We report the detection of several new near-IR flares during the campaign based on {\it HST} observations. The level of near-IR flare activity can be as low as $\sim0.15$ mJy at 1.6 $\mu$m and continuous up to about 40% of the total observing time. Using the NICMOS instrument on the {\it HST}, the {\it XMM-Newton} and CSO observatories, we also detect simultaneous bright X-ray and near-IR flare in which we observe for the first time correlated substructures as well as simultaneous submillimeter and near-IR flaring. X-ray emission is arising from the population of near-IR-synchrotron-emitting relativistic particles which scatter submillimeter seed photons within the inner 10 Schwarzschild radii of Sgr A* up to X-ray energies. In addition, using the inverse Compton scattering picture, we explain the high energy 20-120 keV emission from the direction toward Sgr A*, and the lack of one-to-one X-ray counterparts to near-IR flares, by the variation of the magnetic field and the spectral index distributions of this population of nonthermal particles. In this picture, the evidence for the variability of submillimeter emission during a near-IR flare is produced by the low-energy component of the population of particles emitting synchrotron near-IR emission. Based on the measurements of the duration of flares in near-IR and submillimeter wavelengths, we argue that the cooling could be due to adiabatic expansion with the implication that flare activity may drive an outflow.

Hypermassive neutron stars (HMNSs) -- equilibrium configurations supported against collapse by rapid differential rotation -- are possible transient remnants of binary neutron star mergers. Using newly developed codes for magnetohydrodynamic simulations in dynamical spacetimes, we are able to track the evolution of a magnetized HMNS in full general relativity for the first time. We find that secular angular momentum transport due to magnetic braking and the magnetorotational instability results in the collapse of an HMNS to a rotating black hole, accompanied by a gravitational wave burst. The nascent black hole is surrounded by a hot, massive torus undergoing quasistationary accretion and a collimated magnetic field. This scenario suggests that HMNS collapse is a possible candidate for the central engine of short gamma-ray bursts.

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.

A hypermassive neutron star (HMNS) is a possible transient formed after the merger of a neutron star binary. In the latest magnetohydrodynamic simulations in full general relativity, we find that a magnetized HMNS undergoes `delayed' collapse to a rotating black hole (BH) as a result of angular momentum transport via magnetic braking and the magnetorotational instability. The outcome is a BH surrounded by a massive, hot torus with a collimated magnetic field. The torus accretes onto the BH at a quasi-steady accretion rate ~10 solar mass/s; the lifetime of the torus is ~10 ms. The torus has a temperature \sim 10^{12} K, leading to copious neutrino-antineutrino thermal radiation. Therefore, the collapse of an HMNS is a promising scenario for generating short-duration gamma-ray bursts and an accompanying burst of gravitational waves and neutrinos.

Numerical relativity is an essential tool for solving Einstein's equations of general relativity for dynamical systems characterized by high velocities and strong gravitational fields. The implementation of new algorithms that can solve these nonlinear equations in 3+1 dimensions has enabled us to tackle many long-standing problems of astrophysical interest, leading to an explosion of important new results. Numerical relativity has been used to simulate the evolution of a diverse array of physical systems, including coalescing black hole and neutron star binaries, rotating and collapsing compact objects (stars, collisionless clusters, and scalar fields), and magnetic and viscous stars, to name a few. Numerical relativity has been exploited to address fundamental points of principle, including critical phenomena and cosmic censorship. It holds great promise as a guide for interpreting observations of gravitational waves and gamma-ray bursts and identifying the sources of such radiation. Highlights of a few recent developments in numerical relativity are sketched in this brief overview.

We consider a gravitational wave oscillating in an initially homogeneous, magnetized fluid. The fluid is perfectly conducting and isentropic, and the magnetic field is initially uniform. We find analytic solutions for the case in which the gravitational wave is linear and unaffected by the background fluid and field. Our solutions show how gravitational waves can excite three magnetohydrodynamic (MHD) modes in the fluid: Alfven waves, and both fast and slow magnetosonic waves. Our analytic solutions are particularly useful for testing numerical codes designed to treat general relativistic MHD in dynamical spacetimes, as we demonstrate in a companion paper.

Many problems at the forefront of theoretical astrophysics require the treatment of magnetized fluids in dynamical, strongly curved spacetimes. Such problems include the origin of gamma-ray bursts, magnetic braking of differential rotation in nascent neutron stars arising from stellar core collapse or binary neutron star merger, the formation of jets and magnetized disks around newborn black holes, etc. To model these phenomena, all of which involve both general relativity (GR) and magnetohydrodynamics (MHD), we have developed a GRMHD code capable of evolving MHD fluids in dynamical spacetimes. Our code solves the Einstein-Maxwell-MHD system of coupled equations in axisymmetry and in full 3+1 dimensions. We evolve the metric by integrating the BSSN equations, and use a conservative, shock-capturing scheme to evolve the MHD equations. Our code gives accurate results in standard MHD code-test problems, including magnetized shocks and magnetized Bondi flow. To test our code's ability to evolve the MHD equations in a dynamical spacetime, we study the perturbations of a homogeneous, magnetized fluid excited by a gravitational plane wave, and we find good agreement between the analytic and numerical solutions.

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.

If supermassive black holes (SMBHs) are the energy sources that power quasars and active galactic nuclei, then QSO SDSS 1148+5251, the quasar with the highest redshift (z_QSO=6.43), hosts a supermassive black hole formed within 0.9 Gyr after the Big Bang. This requirement places constraints on the cosmological formation of SMBHs, believed to grow from smaller initial seeds by a combination of accretion and mergers. We focus on gas accretion onto seeds produced by the collapse of Pop III stars at high redshift. We incorporate the results of recent relativistic, MHD simulations of disk accretion onto Kerr black holes to track the coupled evolution of the masses and spins of the holes. We allow for an additional amplification of ~10^4 in the mass of a typical seed due to mergers, consistent with recent Monte Carlo simulations of hierarchical mergers of cold, dark matter halos containing black hole seeds. We find that the growth of Pop III black hole seeds to ~10^9 M_sun by z_QSO = 6.43 favors MHD accretion disks over standard thin disks. MHD disks tend to drive the holes to a submaximal equilibrium spin rate a/M ~ 0.95 and radiation efficiency epsilon_M ~ 0.2, while standard thin disks drive them to maximal spin (a/M = 1) and efficiency (epsilon_M = 0.42). This small difference in efficiency results in a huge difference in mass amplification by accretion at the Eddington limit. The MHD equilibrium efficiency is consistent with the observed ratio of the QSO plus AGN luminosity density to the local SMBH mass density. Our prototype analysis is designed to stimulate the incorporation of results from relativistic stellar collapse and MHD accretion simulations in future Monte Carlo simulations of hierarchical structure formation to better determine the cosmological role of SMBHs and their mass and spin distributions.

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.

Uniformly rotating, homogeneous, incompressible Maclaurin spheroids that spin sufficiently rapidly are secularly unstable to nonaxisymmetric, bar-mode perturbations when viscosity is present. The intuitive explanation is that energy dissipation by viscosity can drive an unstable spheroid to a stable, triaxial configuration of lower energy - a Jacobi ellipsoid. But what about rapidly rotating compressible stars? Unlike incompressible stars, which contain no internal energy and therefore immediately liberate all the energy dissipated by viscosity, compressible stars have internal energy and can retain the dissipated energy as internal heat. Now compressible stars that rotate sufficiently rapidly and also manage to liberate this dissipated energy very quickly are known to be unstable to bar-mode perturbations, like their incompressible counterparts. But what is the situation for rapidly rotating compressible stars that have very long cooling timescales, so that all the energy dissipated by viscosity is retained as heat, whereby the total energy of the star remains constant on a secular (viscous) evolution timescale? Are such stars also unstable to the nonlinear growth of bar modes, or is the viscous heating sufficient to cause them to expand, drive down the ratio of rotational kinetic to gravitational potential energy T/|W| ~ 1/R, where R is the equatorial radius, and turn off the instability before it gets underway? If the instability still arises in such stars, at what rotation rate do they become unstable, and to what final state do they evolve? We provide answers to these questions in the context of the compressible ellipsoid model for rotating stars. The results should serve as useful guides for numerical simulations in 3+1 dimensions for rotating stars containing viscosity.

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.

Simulations in general relativity show that the outcome of collapse of a marginally unstable, uniformly rotating star spinning at the mass-shedding limit depends critically on the equation of state. For a very stiff equation of state, which is likely to characterize a neutron star, essentially all of the mass and angular momentum of the progenitor are swallowed by the Kerr black hole formed during the collapse, leaving nearly no residual gas to form a disk. For a soft equation of state with an adiabatic index \Gamma - 4/3 << 1, which characterizes a very massive or supermassive star supported predominantly by thermal radiation pressure, as much as 10% of the mass of the progenitor avoids capture and goes into a disk about the central hole. We present a semi-analytic calculation that corroborates these numerical findings and shows how the final outcome of such a collapse may be determined from simple physical considerations. In particular, we employ a simple energy variational principle with an approximate, post-Newtonian energy functional to determine the structure of a uniformly rotating, polytropic star at the onset of collapse as a function of polytropic index n, where \Gamma = 1+1/n. We then use this data to calculate the mass and spin of the final black hole and ambient disk. We show that the fraction of the total mass that remains in the disk falls off sharply as 3-n (equivalently, \Gamma - 4/3) increases.

Viscosity and magnetic fields drive differentially rotating stars toward uniform rotation, which has important consequences in many astrophysical contexts. For example, merging binary neutron stars can form a "hypermassive" remnant, i.e. a differentially rotating star with a mass greater than the maximum allowed by uniform rotation. The removal of the centrifugal support provided by differential rotation can lead to delayed collapse of the remnant to a black hole, accompanied by a burst of gravitational radiation. Both magnetic fields and viscosity alter the structure of differentially rotating stars on secular timescales, making numerical tracking of the evolution difficult. Here, we present the first simulations of rapidly rotating stars with shear viscosity in full general relativity. We self-consistently include viscosity in our relativistic hydrodynamic code to solve the relativistic Navier-Stokes equations both in axisymmetry and in full 3+1 dimensions. In axisymmetry, we follow secular evolution with high resolution over dozens of rotation periods (thousands of M). We find that viscosity operating in a hypermassive star generically leads to the formation of a compact, uniformly rotating core surrounded by a low-density disk. These uniformly rotating cores are often, but not always, unstable to gravitational collapse. In the unstable cases, we follow the collapse and determine the mass and the spin of the final black hole and ambient disk. In all cases studied, the rest mass of the resulting disk is found to be 10-20% of the original star, whether surrounding a uniformly rotating core or a rotating black hole. This study foreshadows more detailed, future simulations of secular processes, including magnetic effects, in relativistic stars.

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 present a new numerical scheme for solving the initial value problem for quasiequilibrium binary neutron stars allowing for arbitrary spins. We construct sequences of circular-orbit binaries of varying separation, keeping the rest mass and circulation constant along each sequence. The spin angular frequency of the stars is shown to vary along the sequence, a result that can be derived analytically in the PPN limit. This spin effect, in addition to leaving an imprint on the gravitational waveform emitted during binary inspiral, is measurable in the electromagnetic signal if one of the stars is a pulsar visible from Earth.

We present a numerical code designed to study astrophysical phenomena involving dynamical spacetimes containing black holes in the presence of relativistic hydrodynamic matter. We present evolutions of the collapse of a fluid star from the onset of collapse to the settling of the resulting black hole to a final stationary state. In order to evolve stably after the black hole forms, we excise a region inside the hole before a singularity is encountered. This excision region is introduced after the appearance of an apparent horizon, but while a significant amount of matter remains outside the hole. We test our code by evolving accurately a vacuum Schwarzschild black hole, a relativistic Bondi accretion flow onto a black hole, Oppenheimer-Snyder dust collapse, and the collapse of nonrotating and rotating stars. These systems are tracked reliably for hundreds of M following excision, where M is the mass of the black hole. We perform these tests both in axisymmetry and in full 3+1 dimensions. We then apply our code to study the effect of the stellar spin parameter J/M^2 on the final outcome of gravitational collapse of rapidly rotating n = 1 polytropes. We find that a black hole forms only if J/M^2<1, in agreement with previous simulations. When J/M^2>1, the collapsing star forms a torus which fragments into nonaxisymmetric clumps, capable of generating appreciable ``splash'' gravitational radiation.

We study the magnetic braking and viscous damping of differential rotation in incompressible, uniform density stars in general relativity. Differentially rotating stars can support significantly more mass in equilibrium than nonrotating or uniformly rotating stars. The remnant of a binary neutron star merger or supernova core collapse may produce such a "hypermassive" neutron star. Although a hypermassive neutron star may be stable on a dynamical timescale, magnetic braking and viscous damping of differential rotation will ultimately alter the equilibrium structure, possibly leading to delayed catastrophic collapse. Here we consider the slow-rotation, weak-magnetic field limit in which E_rot << E_mag << W, where E_rot is the rotational kinetic energy, E_mag is the magnetic energy, and W is the gravitational binding energy of the star. We assume the system to be axisymmetric and solve the MHD equations in both Newtonian gravitation and general relativity. Toroidal magnetic fields are generated whenever the angular velocity varies along the initial poloidal field lines. We find that the toroidal fields and angular velocities oscillate independently along each poloidal field line, which enables us to transform the original 2+1 equations into 1+1 form and solve them along each field line independently. The incoherent oscillations on different field lines stir up turbulent-like motion in tens of Alfven timescales ("phase mixing"). In the presence of viscosity, the stars eventually are driven to uniform rotation, with the energy contained in the initial differential rotation going into heat. Our evolution calculations serve as qualitative guides and benchmarks for future, more realistic MHD simulations in full 3+1 general relativity.

We consider a subset of the physical processes that determine the spin j = a/M of astrophysical black holes. These include: (1) Initial conditions. Recent models suggest that the collapse of supermassive stars are likely to produce black holes with j ~ 0.7. (2) Major mergers. The outcome of a nearly equal mass black hole-black hole merger is not yet known, but we review the current best guesses and analytic bounds. (3) Minor mergers. We recover the result of Blandford & Hughes that accretion of small companions with isotropically distributed orbital angular momenta results in spindown, with j ~ M^{-7/3}. (4) Accretion. We present new results from fully relativistic magnetohydrodynamic accretion simulations. These show that, at least for one sequence of flow models, spin equilibrium (dj/dt = 0) is reached for j ~ 0.9, far less than the canonical value 0.998 of Thorne that was derived in the absence of MHD effects. This equilibrium value may not apply to all accretion flows, particularly thin disks. Nevertheless, it opens the possibility that black holes that have grown primarily through accretion are not maximally rotating.

Differential rotation in stars generates toroidal magnetic fields whenever an initial seed poloidal field is present. The resulting magnetic stresses, along with viscosity, drive the star toward uniform rotation. This magnetic braking has important dynamical consequences in many astrophysical contexts. For example, merging binary neutron stars can form "hypermassive" remnants supported against collapse by differential rotation. The removal of this support by magnetic braking induces radial fluid motion, which can lead to delayed collapse of the remnant to a black hole. We explore the effects of magnetic braking and viscosity on the structure of a differentially rotating, compressible star, generalizing our earlier calculations for incompressible configurations. The star is idealized as a differentially rotating, infinite cylinder supported initially by a polytropic equation of state. The gas is assumed to be infinitely conducting and our calculations are performed in Newtonian gravitation. Though highly idealized, our model allows for the incorporation of magnetic fields, viscosity, compressibility, and shocks with minimal computational resources in a 1+1 dimensional Lagrangian MHD code. Our evolution calculations show that magnetic braking can lead to significant structural changes in a star, including quasistatic contraction of the core and ejection of matter in the outermost regions to form a wind or an ambient disk. These calculations serve as a prelude and a guide to more realistic MHD simulations in full 3+1 general relativity.

We introduce a new numerical scheme for solving the initial value problem for quasiequilibrium binary neutron stars allowing for arbitrary spins. The coupled Einstein field equations and equations of relativistic hydrodynamics are solved in the Wilson-Mathews conformal thin sandwich formalism. We construct sequences of circular-orbit binaries of varying separation, keeping the rest mass and circulation constant along each sequence. Solutions are presented for configurations obeying an n=1 polytropic equation of state and spinning parallel and antiparallel to the orbital angular momentum. We treat stars with moderate compaction ((m/R) = 0.14) and high compaction ((m/R) = 0.19). For all but the highest circulation sequences, the spins of the neutron stars increase as the binary separation decreases. Our zero-circulation cases approximate irrotational sequences, for which the spin angular frequencies of the stars increases by 13% (11%) of the orbital frequency for (m/R) = 0.14 ((m/R) = 0.19) by the time the innermost circular orbit is reached. In addition to leaving an imprint on the inspiral gravitational waveform, this spin effect is measurable in the electromagnetic signal if one of the stars is a pulsar visible from Earth.

There is compelling evidence that supermassive black holes (SMBHs) exist. Yet the origin of these objects, or their seeds, is still unknown. We are performing general relativistic simulations of gravitational collapse to black holes in different scenarios to help reveal how SMBH seeds might arise in the universe. SMBHs with ~ 10^9 solar masses must have formed by z > 6, or within 10^9 yrs after the Big Bang, to power quasars. It may be difficult for gas accretion to build up such a SMBH by this time unless the initial seed black hole already has a substantial mass. One plausible progenitor of a massive seed black hole is a supermassive star (SMS). We have followed the collapse of a SMS to a SMBH by means of 3D hydrodynamic simulations in post-Newtonian gravity and axisymmetric simulations in full general relativity. The initial SMS of arbitrary mass M in these simulations rotates uniformly at the mass--shedding limit and is marginally unstable to radial collapse. The final black hole mass and spin are determined to be M_h/M ~ 0.9 and J_h/M_h^2 ~ 0.75. The remaining mass goes into a disk of mass M_{disk}/M ~ 0.1. This disk arises even though the total spin of the progenitor star, J/M^2 = 0.97, is safely below the Kerr limit. The collapse generates a mild burst of gravitational radiation. Nonaxisymmetric bars or one-armed spirals may arise during the quasi-stationary evolution of a SMS, during its collapse, or in the ambient disk about the hole, and are potential sources of quasi-periodic waves, detectable by LISA.

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.

There is compelling evidence that supermassive black holes exist. Yet the origin of these objects, or their seeds, is still unknown. We discuss several plausible scenarios for forming the seeds of supermassive black holes. These include the catastrophic collapse of supermassive stars, the collapse of relativistic clusters of collisionless particles or stars, the gravothermal evolution of dense clusters of ordinary stars or stellar-mass compact objects, and the gravothermal evolution of self-interacting dark matter halos. Einstein's equations of general relativity are required to describe key facets of these scenarios, and large-scale numerical simulations are performed to solve them.

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.

The collapse of a uniformaly rotating, supermassive star (SMS) to a supermassive black hole (SMBH) has been followed recently by means of hydrodynamic simulations in full general relativity. The initial SMS of arbitrary mass M in these simulations rotates uniformly at the mass--shedding limit and is marginally unstable to radial collapse. The final black hole has mass M_h/M = 0.9 and and spin J_h/M_h^2 = 0.75, approximately. The remaining mass goes into a disk of mass M_disk/M = 0.1, also approximately. Here we show that these black hole and disk parameters can be calculated analytically from the initial stellar density and angular momentum distribution. The analytic calculation thereby corroborates and provides a simple physical explanation for the computational discovery that SMS collapse inevitably terminates in the simultaneous formation of a SMBH and a rather substantial ambient disk. This disk arises even though the total spin of the progenitor star, J/M^2 = 0.97, is safely below the Kerr limit. The calculation performed here applies to any marginally unstable n = 3 polytrope uniformly rotating at the break--up speed, independent of stellar mass or the source of internal pressure. It illustrates how the black hole and disk parameters can be determined for the collapse of other types of stars with different initial density and rotation profiles.

We follow the collapse in axisymmetry of a uniformly rotating, supermassive star (SMS) to a supermassive black hole (SMBH) in full general relativity. The initial SMS of arbitrary mass $M$ is marginally unstable to radial collapse and rotates at the mass-shedding limit. The collapse proceeds homologously early on and results in the appearance of an apparent horizon at the center. Although our integration terminates before final equilibrium is achieved, we determine that the final black hole will contain about 90% of the total mass of the system and have a spin parameter $J/M^2 \sim 0.75$. The remaining gas forms a rotating disk about the nascent hole.

A central supermassive black hole (SMBH) with a mass $10^6-10^9 M_\odot$ appears to be a common feature in nearby galaxies and the likely power source in quasars and active galactic nuclei. We demonstrate that the formation of a central black hole is a natural and inevitable consequence of the gravothermal catastrophe in a self-interacting dark matter (SIDM) halo. Through gravothermal evolution driven by collisional relaxation, an SIDM halo will form a massive inner core whose density and velocity dispersion will increase secularly in time. Eventually, the inner core arrives at a relativistic radial instability and undergoes dynamical collapse to a black hole. The initial mass of the black hole will be $10^{-8}-10^{-6}$ of the total mass of the halo. We show that if at formation the overdensity in the SIDM halo is not too large, SMBHs in the observed mass range can form directly in very massive halos following core collapse. Alternatively, with large overdensities, moderate mass halos undergo core collapse to form central seed black holes of intermediate mass, and these holes can then merge and/or accrete to reach the SMBH range. Forming SMBHs by core collapse in SIDM halos requires no baryons, no prior epoch of star formation and no other mechanism of forming black holes seeds.

We study the evolution of an isolated, spherical halo of self-interacting dark matter (SIDM) in the gravothermal fluid formalism. We show that the thermal relaxation time, $t_r$, of a SIDM halo with a central density and velocity dispersion of a typical dwarf galaxy is significantly shorter than its age. We find a self-similar solution for the evolution of a SIDM halo in the limit where the mean free path between collisions, $\lambda$, is everywhere longer than the gravitational scale height, $H$. Typical halos formed in this long mean free path regime relax to a quasistationary gravothermal density profile characterized by a nearly homogeneous core and a power-law halo where $\rho \propto r^{-2.19}$. We solve the more general time-dependent problem and show that the contracting core evolves to sufficiently high density that $\lambda$ inevitably becomes smaller than $H$ in the innermost region. The core undergoes secular collapse to a singular state (the ``gravothermal catastrophe'') in a time $t_{coll} \approx 290 t_r$, which is longer than the Hubble time for a typical dark matter-dominated galaxy core at the present epoch. Our model calculations are consistent with previous, more detailed, N-body simulations for SIDM, providing a simple physical interpretation of their results and extending them to higher spatial resolution and longer evolution times. At late times, mass loss from the contracting, dense inner core to the ambient halo is significantly moderated, so that the final mass of the inner core may be appreciable when it becomes relativistic and radially unstable to dynamical collapse to a black hole.

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 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.

The evolution of the r-mode instability is likely to be accompanied by secular kinematic effects which will produce differential rotation with large scale drifts of fluid elements, mostly in the azimuthal direction. As first discussed by Rezzolla, Lamb and Shapiro 2000, the interaction of these secular velocity fields with a pre-existing neutron star magnetic field could result in the generation of intense and large scale toroidal fields. Following their derivation in the companion paper, we here discuss the numerical solution of the evolution equations for the magnetic field. The values of the magnetic fields obtained in this way are used to estimate the conditions under which the r-mode instability might be prevented or suppressed. We also assess the impact of the generation of large magnetic fields on the gravitational wave detectability of r-mode unstable neutron stars. Our results indicate that the signal to noise ratio in the detection of gravitational waves from the r-mode instability might be considerably decreased if the latter develops in neutron stars with initial magnetic fields larger than 10^10 G.

The instability of r-mode oscillations in rapidly rotating neutron stars has attracted attention as a potential mechanism for producing high frequency, almost periodic gravitational waves. The analyses carried so far have shown the existence of these modes and have considered damping by shear and bulk viscosity. However, the magnetohydrodynamic coupling of the modes with a stellar magnetic field and its role in the damping of the instability has not been fully investigated yet. Following our introductory paper (Rezzolla, Lamb and Shapiro 2000), we here discuss in more detail the existence of secular higher-order kinematical effects which will produce toroidal fluid drifts. We also define the sets of equations that account for the time evolution of the magnetic fields produced by these secular velocity fields and show that the magnetic fields produced can reach equipartition in less than a year. The full numerical calculations as well as the evaluation of the impact of strong magnetic fields on the onset and evolution of the r-mode instability will be presented in a companion paper.

We investigate the equilibrium and stability of supermassive stars of mass $M \agt 10^5M_{\odot}$ in binary systems. We find that corotating binaries are secularly unstable for close, circular orbits with $r \alt 4R(M/10^6M_{\odot})^{1/6}$ where $r$ is the orbital separation and $R$ the stellar radius. We also show that corotation cannot be achieved for distant orbits with $r \agt 12 R (M/10^6M_{\odot})^{-11/24}$, since the timescale for viscous angular momentum transfer associated with tidal torques is longer than the evolution timescale due to emission of thermal radiation. These facts suggest that the allowed mass range and orbital separation for corotating supermassive binary stars is severely restricted. In particular, for supermassive binary stars of large mass $M \agt 6\times 10^6M_{\odot}$, corotation cannot be achieved, as viscosity is not adequate to mediate the transfer between orbital and spin angular momentum. One possible outcome for binary supermassive stars is the onset of quasi-radial, relativistic instability which drives each star to collapse prior to merger: We discuss alternative outcomes of collapse and possible spin states of the resulting black holes. We estimate the frequency and amplitude of gravitational waves emitted during several inspiral and collapse scenarios.

A newly formed black hole may be directly identified if late-time accretion of material from the base of the ejected envelope generates a luminosity that is observable in the tail of the supernova light curve. In this work we estimate the rate at which events where the black hole ``emerges'' in the supernova light curve can be detected with present capabilities. Our investigation is based on an analytical model of the accretion luminosity at emergence as a function of progenitor mass, coupled to the inferred rate of observed Type II supernovae in nearby galaxies. We find through a parameter survey that under optimistic assumptions the potential rate of observable events can be as high as several per year. However, supernovae which produce black holes are also likely to be low energy explosions and therefore subluminous, as was the case for the best candidate to date, SN1997D. If black hole-forming supernovae are underdetected owing to lower luminosities, the rate of observing black hole emergence is probably not larger than once every few years. We therefore emphasize the importance of dedicated searches for nearby supernovae as well as faint supernovae projects for improving the prospects of observationally certifying the supernova--black hole connection.

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.

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.

Diffferentially rotating stars can support significantly more mass in equilibrium than nonrotating or uniformly rotating stars, according to general relativity. The remnant of a binary neutron star merger may give rise to such a ``hypermassive'' object. While such a star may be dynamically stable against gravitational collapse and bar formation, the radial stabilization due to differential rotation is likely to be temporary. Magnetic braking and viscosity combine to drive the star to uniform rotation, even if the seed magnetic field and the viscosity are small. This process inevitably leads to delayed collapse, which will be accompanied by a delayed gravitational wave burst and, possibly, a gamma-ray burst. We provide a simple, Newtonian, MHD calculation of the braking of differential rotation by magnetic fields and viscosity. The star is idealized as a differentially rotating, infinite cylinder consisting of a homogeneous, incompressible conducting gas. We solve analytically the simplest case in which the gas has no viscosity and the star resides in an exterior vacuum. We treat numerically cases in which the gas has internal viscosity and the star is embedded in an exterior, low-density, conducting medium. Our evolution calculations are presented to stimulate more realistic MHD simulations in full 3+1 general relativity. They serve to identify some of the key physical and numerical parameters, scaling behavior and competing timescales that characterize this important process.

We present a new set of spectroscopic and photometric data extending the observations of SN 1997D to over 400 days after the explosion. These observations confirm the peculiar properties of SN 1997D, such as the very low abundance of 56Co (0.002 Msolar) and the low expansion velocity of the ejecta (~1000 km/s). We discuss the implications of these observations for the character of the progenitor and the nature of the remnant showing that a Crab-like pulsar or an accreting neutron star formed in the explosion of a low mass progenitor should already have produced a detectable luminosity at this epoch, in contrast with photometric data. On the other hand, the explosion of a high mass progenitor with the formation of a black hole is consistent with the available observations. The consequences of this conclusion regarding the nature of the explosion and the prospects of directly identifying the black hole are also addressed.

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.

Thermal emission from a rotating, supermassive star will cause the configuration to contract slowly and spin up. If internal viscosity and magnetic fields are sufficiently weak, the contracting star will rotate differentially. For each of the six initial angular momentum distributions considered, a cooling supermassive star will likely encounter the dynamical bar mode instability, which may trigger the growth of nonaxisymmetric bars. This scenario will lead to the generation of long wavelength gravitational waves, which could be detectable by future space-based laser interferometers like LISA.

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.

The existence of supermassive black holes is supported by a growing body of observations. Supermassive black holes and their formation events are likely candidates for detection by proposed long-wavelength, space-based gravitational wave interferometers like LISA. However, the nature of the progenitors of supermassive black holes is rather uncertain. Supermassive black hole formation scenarios that involve either the stellar dynamical evolution of dense clusters or the hydrodynamical evolution of supermassive stars have been proposed. Each of these formation scenarios is reviewed and the evolution of supermassive stars is then examined in some detail. Supermassive stars that rotate uniformly during their secular cooling phase will spin up to the mass-shedding limit and eventually contract to the point of relativistic collapse. Supermassive stars that rotate differentially as they cool will likely encounter the dynamical bar mode instability prior to the onset of relativistic collapse. A supermassive star that undergoes this bar distortion, prior to or during collapse, may be a strong source of quasiperiodic, long-wavelength gravitational radiation.

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.

White dwarfs and neutron stars are stellar objects with masses comparable to that of our sun. However, as the endpoint stages of stellar evolution, these objects do not sustain any thermonuclear burning and therefore can no longer support the gravitational load of their own mass by generating thermal pressure. Rather, matter in their interiors is compressed to much higher densities than commonly found in normal stars, and pressure is created by degenerate fermion kinetic energy and particle interactions. As a result, white dwarfs and neutron stars offer unique cosmic laboratories for studying matter at very high densities. In this review we discuss the basic properties of condensed matter at extreme densities and summarize the extent to which these properties can be examined by observations of compact objects.

If a black hole formed in a core-collapse supernova is accreting material from the base of the envelope, the accretion luminosity could be observable in the supernova light curve. Here we continue the study of matter fall back onto a black hole in the wake of a supernova and examine realistic supernovae models which allow for an early emergence of the accretion luminosity. Such cases may provide a direct observational identification of the black hole formed in the aftermath of the explosion. Our approach combines analytic estimates and fully relativistic, radiation-hydrodynamic numerical computations. We employ a numerical hydrodynamical scaling technique to accommodate the diverse range of dynamical time scales in a single simulation. We find that while in typical Type II supernovae heating by radioactive decays dominates the late-time light curve, low-energy explosions of more massive stars should provide an important exception where the accretion luminosity will emerge while it is still relatively large. Our main focus is on the only current candidate for such an observation, the very unusual SN1997D. Due to the low energy of the explosion and the very small ($2\times10^{-3} M_\sun$) inferred mass of Co56 in the ejected envelope, we find that accretion should become the dominant source of its luminosity during the year 2000. The total luminosity at emergence is expected to lie in the range $0.5-3\times10^{36} $ ergs/s, potentially detectable with HST. We also discuss the more favorable case of explosions which eject negligible amounts of radioactive isotopes and find that the black hole is likely to emerge a few tens of days after the explosion, with a luminosity of $\sim 10^{37} $\ergss.

We show that r-mode oscillations distort the magnetic fields of neutron stars and that their occurrence is likely to be limited by this interaction. If the field is gtrsim 10^{16} (Omega/Omega_B) G, where Omega and Omega_B are the angular velocities of the star and at which mass shedding occurs, r-mode oscillations cannot occur. Much weaker fields will prevent gravitational radiation from exciting r-mode oscillations or damp them on a relatively short timescale by extracting energy from the modes faster than gravitational wave emission can pump energy into them. For example, a 10^{10} G poloidal magnetic field that threads the star's superconducting core is likely to prevent the ell=2 mode from being excited unless Omega exceeds 0.35 Omega_B. If Omega is larger than 0.35 Omega_B initially, the ell=2 mode may be excited but is likely to decay rapidly once Omega falls below 0.35 Omega_B, which happens in lesssim 15^d if the saturation amplitude is gtrsim 0.1. The r-mode oscillations may play an important role in determining the structure of neutron star magnetic fields.

We consider the effect of a positive cosmological constant on spherical gravitational collapse to a black hole for a few simple, analytic cases. We construct the complete Oppenheimer-Snyder-deSitter (OSdS) spacetime, the generalization of the Oppenheimer-Snyder solution for collapse from rest of a homogeneous dust ball in an exterior vacuum. In OSdS collapse, the cosmological constant may affect the onset of collapse and decelerate the implosion initially, but it plays a diminishing role as the collapse proceeds. We also construct spacetimes in which a collapsing dust ball can bounce, or hover in unstable equilibrium, due to the repulsive force of the cosmological constant. We explore the causal structure of the different spacetimes and identify any cosmological and black hole event horizons which may be present.

If a black hole formed in a core-collapse supernova is accreting material from the base of the envelope, the accretion luminosity could be observable in the supernova light curve. We present results of a fully relativistic numerical investigation of the fallback of matter onto a black hole in a supernova and examine conditions which would be favorable for detection of the black hole. In general, heating by radioactive decays is likely to prevent practical detection of the black hole, but we show that low energy explosions of more massive stars may provide an important exception. We emphasize the particular case of SN1997D in NGC1536, for which we predict that the presence of a black hole could be inferred observationally within the next year.

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.

Recent numerical simulations address a conjecture by Shapiro that when two neutron stars collide head-on from rest at infinity, sufficient thermal pressure may be generated to support the hot remnant in quasi-static equilibrium against collapse prior to neutrino cooling. The conjecture is meant to apply even when the total remnant mass exceeds the maximum mass of a cold neutron star. One set of simulations seems to corroborate the conjecture, while another, involving higher mass progenitors each very close to the maximum mass, does not. In both cases the total mass of the remnant exceeds the maximum mass. We point out numerical subtleties in performing such simulations when the progenitors are near the maximum mass; they can explain why the simulations might have difficulty assessing the conjecture in such high-mass cases.

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 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 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 improve and extend Shapiro's model of a relativistic, compact object which is stable in isolation but is driven dynamically unstable by the tidal field of a binary companion. Our compact object consists of a dense swarm of test particles moving in randomly-oriented, initially circular, relativistic orbits about a nonrotating black hole. The binary companion is a distant, slowly inspiraling point mass. The tidal field of the companion is treated as a small perturbation on the background Schwarzschild geometry near the hole; the resulting metric is determined by solving the perturbation equations of Regge and Wheeler and Zerilli in the quasi-static limit. The perturbed spacetime supports Bekenstein's conjecture that the horizon area of a near-equilibrium black hole is an adiabatic invariant. We follow the evolution of the system and confirm that gravitational collapse can be induced in a compact collisionless cluster by the tidal field of a binary companion.

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'').

Coalescing compact binaries with neutron star or black hole components provide the most promising sources of gravitational radiation for detection by the LIGO/VIRGO/GEO/TAMA laser interferometers now under construction. This fact has motivated several different theoretical studies of the inspiral and hydrodynamic merging of compact binaries. Analytic analyses of the inspiral waveforms have been performed in the Post-Newtonian approximation. Analytic and numerical treatments of the coalescence waveforms from binary neutron stars have been performed using Newtonian hydrodynamics and the quadrupole radiation approximation. Numerical simulations of coalescing black hole and neutron star binaries are also underway in full general relativity. Recent results from each of these approaches will be described and their virtues and limitations summarized.

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.

The head-on collision of identical neutron stars from rest at infinity requires a numerical simulation in full general relativity for a complete solution. Undaunted, we provide a relativistic, analytic argument to suggest that during the collision, sufficient thermal pressure is always generated to support the hot remnant in quasi-static stable equilibrium against collapse prior to slow cooling via neutrino emission. Our conclusion is independent of the total mass of the progenitors and holds even if the remnant greatly exceeds the maximum mass of a cold neutron star.

We have performed a series of systematic tests to evaluate the effects of spurious transport in three-dimensional smoothed particle hydrodynamics (SPH) calculations. Our tests investigate (i) particle diffusion, (ii) shock heating, (iii) numerical viscosity, and (iv) angular momentum transport. The results are useful for quantifying the accuracy of the SPH scheme, especially for problems where shear flows or shocks are present, as well as for problems where true hydrodynamic mixing is relevant. We examine the different forms of artificial viscosity (AV) which have been proposed by Monaghan, by Hernquist & Katz, and by Balsara. For each form, our tests suggest a single set of values for the AV parameters $\alpha$ and $\beta$ (coefficients of the linear and quadratic terms) which are appropriate in a large number of situations. We also discuss how these parameters should be adjusted depending on the goals of the particular application. We find that both the Hernquist & Katz and Balsara forms introduce relatively small amounts of numerical viscosity. Furthermore, both Monaghan's and Balsara's AV do well at treating shocks and at limiting the amount of spurious mixing. For these reasons, we endorse the Balsara AV for use in a broad range of applications.

Observations combined with theoretical modeling of the light curve of the recently discovered supernova 1997D in NGC 1536 suggest that it might host a black hole formed in the aftermath of the explosion. We consider some observable consequences of a black hole in SN 1997D and estimate the late--time accretion luminosity of the material which falls back onto the hole. We find that this luminosity, decaying with a characteristic power--law dependence on time, may emerge above the emission of the envelope in just a few years. Its detection would thus provide unmistakable evidence for the presence of a black hole.

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 present the first fully relativistic investigation of matter fallback in a supernova. We investigate spherically symmetric supernova fallback using a relativistic radiation hydrodynamics Lagrangian code that handles radiation transport in all regimes. Our goal is to answer the fundamental question: did SN1987A produce a black hole and, if so, when will the hole become detectable ? We compute the light curve, assuming that a black hole has been formed during the explosion, and compare it with the observations. Our preliminary calculations lack radioactive energy input and adopt a very simple chemical composition (pure hydrogen). As a result, our computed models cannot fit the observed data of SN1987A in detail. Nevertheless, we can show that, during the first hours, the accretion flow is self--regulated and the accretion luminosity stays very close to the Eddington limit. The light curve is completely dominated, during the first few weeks, by the emission of the stellar envelope thermal energy, and resembles that obtained in ``standard'' supernova theory. Only long after hydrogen recombination takes place is there even a chance to actually detect radiation from the accreting black hole above the emission of the expanding envelope. The presence of a black hole is thus not inconsistent with observations to date. Because of the exponential decay of the $^{44}$Ti radioactive heating rate, the date of the emergence of the black hole is not very sensitive to the actual parameters of the models and turns out to be about 1000 years. The bulk of the emission then is expected to be in the visible band, but will be unobservable with present instrumentation. We discuss the implications of our results in connection with the possible emergence of a black hole in other supernovae.

We present results from a fully relativistic investigation of matter fallback in a supernova, assuming that the central compact star is a black hole (BH). Our goal is to answer the following questions: can we infer from the light curve whether a stellar black hole (BH) has formed in the aftermath of a supernova explosion ? Did SN1897A produce a BH and, if so, when will the hole become detectable ? We show that it should not be possible to infer the presence of a BH in the remnant for another several hundred to a thousand years. In the case of SN1987A, we estimate the BH luminosity to be $L\simeq 5\times 10^{34}$ erg s$^{-1}$, which is well below the present day bolometric luminosity of the remnant ($\sim 10^{36}erg/s$; Suntzeff 1997).

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 construct analytic models of incompressible, rigidly rotating stars in PN gravity and study their stability against nonaxisymmetric Jacobi-like bar modes. PN configurations are modeled by homogeneous triaxial ellipsoids and the metric is obtained as a solution of Einstein's equations in 3+1 ADM form. We use an approximate subset of the equations well-suited to numerical integration for strong field, 3D configurations in quasi--equilibrium. These equations are exact at PN order, and admit an analytic solution for homogeneous ellipsoids. In this paper we present this solution, as well as analytic functionals for the conserved global quantities, M, M_0 and J. By using a variational principle we construct sequences of axisymmetric equilibria of constant density and rest mass, i.e. the PN generalization of Maclaurin spheroids, which are compared to other PN and full relativistic sequences presented by previous authors. We then consider nonaxisymmetric ellipsoidal deformations of the configurations, holding J constant and the rotation uniform, and we locate the point at which the bar modes will be driven secularly unstable by a dissipative agent like viscosity. We find that the value of the eccentricity, as well as the ratios \Omega^2/(\pi\rho_0) and T/|W|, defined invariantly, all increase at the onset of instability as the stars become more relativistic. Since higher degrees of rotation are required to trigger a viscosity-driven bar mode as the star's compactness increases, the effect of GR is to weaken the instability, at least to PN order. This behavior is opposite to that found for secular instability via Dedekind-like modes driven unstable by gravitational radiation, supporting the suggestion that in GR, nonaxisymmetric modes driven unstable by viscosity and gravitational radiation may no longer coincide.

We present a simple model illustrating how a highly relativistic, compact object which is stable in isolation can be driven dynamically unstable by the tidal field of a binary companion. Our compact object consists of a test-particle in a relativistic orbit about a black hole; the binary companion is a distant point mass. Our example is presented in light of mounting theoretical opposition to the possibility that sufficiently massive, binary neutron stars inspiraling from large distance can collapse to form black holes prior to merger. Our strong-field model suggests that first order post-Newtonian treatments of binaries, and stability analyses of binary equilibria based on orbit-averaged, mean gravitational fields, may not be adequate to rule out this possibility.

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.

Using an energy variational method, we calculate quasi-equilibrium configurations of binary neutron stars modeled as compressible triaxial ellipsoids obeying a polytropic equation of state. Our energy functional includes terms both for the internal hydrodynamics of the stars and for the external orbital motion. We add the leading post-Newtonian (PN) corrections to the internal and gravitational energies of the stars, and adopt hybrid orbital terms which are fully relativistic in the test-mass limit and always accurate to PN order. The total energy functional is varied to find quasi-equilibrium sequences for both corotating and irrotational binaries in circular orbits. We examine how the orbital frequency at the innermost stable circular orbit depends on the polytropic index n and the compactness parameter GM/Rc^2. We find that, for a given GM/Rc^2, the innermost stable circular orbit along an irrotational sequence is about 17% larger than the innermost secularly stable circular orbit along the corotating sequence when n=0.5, and 20% larger when n=1. We also examine the dependence of the maximum neutron star mass on the orbital frequency and find that, if PN tidal effects can be neglected, the maximum equilibrium mass increases as the orbital separation decreases.

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 explore the evolution of collisionally merged stars in the blue straggler region of the HR diagram. The starting models for our stellar evolution calculations are the results of the smoothed particle hydrodynamics (SPH) simulations of parabolic collisions between main sequence stars performed by Lombardi, Rasio and Shapiro (1996). Since SPH and stellar evolution codes employ different and often contradictory approximations, it is necessary to treat the evolution of these products carefully. The mixture and disparity of the relevant timescales (hydrodynamic, thermal relaxation and nuclear burning) and of the important physical assumptions between the codes makes the combined analysis of the problem challenging, especially during the initial thermal relaxation of the star. In particular, the treatment of convection is important, and semiconvection must be modeled in some detail. The products of seven head-on collisions are evolved through their initial thermal relaxation, and then through the main sequence phase to the base of the giant branch. Their evolutionary tracks are presented. In contrast to the assumptions in previous work, these collision products do not develop substantial convective regions during their thermal relaxation, and therefore are not mixed significantly after the collision.

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.

We report the results of new SPH calculations of parabolic collisions between two main-sequence stars in a globular cluster. Such collisions are directly relevant to the formation of blue stragglers. In particular, we consider parent stars of mass $M/M_{TO}=0.2,0.5,0.75$, and $1$, where $M_{TO}$ is the cluster turnoff mass (typically about $0.8\,M_\odot$). Low-mass stars (with $M=0.2 M_{TO}$ or $0.5 M_{TO}$) are assumed to be fully convective and are therefore modeled as $n=1.5$ polytropes. Stars at the turnoff (with $M= M_{TO}$) are assumed to be mostly radiative and are modeled as $n=3$ polytropes. Intermediate-mass stars (with $M=0.75 M_{TO}$) are modeled as composite polytropes consisting of a radiative core with polytropic index $n=3$ and a convective envelope with $n=1.5$. We focus our study on the question of hydrodynamic mixing of helium and hydrogen, which plays a crucial role in determining the observable characteristics of blue stragglers. In all cases we find that there is negligible hydrodynamic mixing of helium into the outer envelope of the merger remnant. The amount of hydrogen mixed into the core of the merger depends strongly on the entropy profiles of the two colliding stars. In contrast to what has been done in previous studies, blue stragglers formed by direct stellar collisions should not necessarily be assumed to have initially homogeneous composition profiles.

We report the results of new SPH calculations of parabolic collisions between main-sequence (MS) stars. The stars are assumed to be close to the MS turn-off point in a globular cluster and are therefore modeled as $n=3$, $\Gamma=5/3$ polytropes. We find that the high degree of central mass concentration in these stars has a profound effect on the hydrodynamics. In particular, very little hydrodynamic mixing occurs between the dense, helium-rich inner cores and the outer envelopes. As a result, and in contrast to what has been assumed in previous studies, blue stragglers formed by direct stellar collisions are not necessarily expected to have anomalously high helium abundances in their envelopes or to have their cores replenished with fresh hydrogen fuel.

In addition to their possible relevance to gamma-ray bursts, coalescing binary neutron stars have long been recognized as important sources of gravitational radiation that should become detectable with the new generation of laser interferometers such as LIGO. Hydrodynamics plays an essential role near the end of the coalescence when the two stars finally merge together into a single object. The shape of the corresponding burst of gravitational waves provides a direct probe into the interior structure of a neutron star and the nuclear equation of state. The interpretation of the gravitational waveform data will require detailed theoretical models of the complicated three-dimensional hydrodynamic processes involved. Here we review the results of our recent work on this problem, using both approximate quasi-analytic methods and large-scale numerical hydrodynamics calculations on supercomputers. We also discuss briefly the coalescence of white-dwarf binaries, which are also associated with a variety of interesting astrophysical phenomena.

We employ an approximate treatment of dissipative hydrodynamics in three dimensions to study the coalescence of binary neutron stars driven by the emission of gravitational waves. The stars are modeled as compressible ellipsoids obeying a polytropic equation of state; all internal fluid velocities are assumed to be linear functions of the coordinates. The hydrodynamic equations then reduce to a set of coupled ordinary differential equations for the evolution of the principal axes of the ellipsoids, the internal velocity parameters and the binary orbital parameters. Gravitational radiation reaction and viscous dissipation are both incorporated. We set up exact initial binary equilibrium configurations and follow the transition from the quasi-static, secular decay of the orbit at large separation to the rapid dynamical evolution of the configurations just prior to contact. A hydrodynamical instability resulting from tidal interactions significantly accelerates the coalescence at small separation, leading to appreciable radial infall velocity and tidal lag angles near contact. This behavior is reflected in the gravitational waveforms and may be observable by gravitational wave detectors under construction.

We study the secular evolution and gravitational wave signature of a newly-formed, rapidly rotating neutron star. The neutron star may arise from core collapse in a massive star or from the accretion-induced collapse of a white dwarf. After a brief dynamical phase, the nascent neutron star settles into an axisymmetric, secularly unstable equilibrium state. Gravitational radiation drives the star to a nonaxisymmetric, stationary equilibrium configuration via the bar-mode instability. The emitted quasi-periodic gravitational waves have a unique signature: the wave frequency sweeps downward from a few hundred Hertz to zero, while the wave amplitude increases from zero to a maximum and then decays back to zero. Such a wave signal could be detected by broad-band gravitational wave interferometers currently being constructed. We also characterize two other types of gravitational wave signals that could arise in principle from a rapidly rotating, secularly unstable neutron star: a high-frequency ($f\go 1000$ Hz) wave which increases the pattern-speed of the star, and a wave that actually increases the angular momentum of the star.

We present a new numerical study of the equilibrium and stability properties of close binary systems using the smoothed-particle hydrodynamics (SPH) technique. We adopt a simple polytropic equation of state $p=K\rho^\gam$ with $\gam=5/3$ and $K=\,$constant within each star, applicable to low-mass degenerate dwarfs as well as low-mass main-sequence stars. Along a sequence of binary equilibrium configurations for two identical stars, we demonstrate the existence of both secular and dynamical instabilities, confirming directly the results of recent analytic work. We use the SPH method to calculate the nonlinear development of the dynamical instability and to determine the final fate of the system. We find that the two stars merge together into a single, rapidly rotating object in just a few orbital periods. Equilibrium sequences are also constructed for systems containing two nonidentical stars. These sequences terminate at a Roche limit, which we can determine very accurately using SPH. For two low-mass main-sequence stars with mass ratio $q\lo0.4$ we find that the (synchronized) Roche limit configuration is secularly unstable. Degenerate binary configurations remain hydrodynamically stable all the way to the Roche limit for all mass ratios $q\ne1$. Dynamically unstable mass transfer can also lead to the rapid coalescence of a binary system, but the details of the hydrodynamic evolution are quite different. We discuss the implications of our results for double white-dwarf and W Ursae Majoris systems.

We develop a new formalism to study the dynamics of fluid polytropes in three dimensions. The stars are modeled as compressible ellipsoids and the hydrodynamic equations are reduced to a set of ordinary differential equations for the evolution of the principal axes and other global quantities. Both viscous dissipation and the gravitational radiation reaction are incorporated. We establish the validity of our approximations and demonstrate the simplicity and power of the method by rederiving a number of known results concerning the stability and dynamical oscillations of rapidly rotating polytropes. In particular, we present a generalization to compressible fluids of Chandrasekhar's classical results for the secular and dynamical instabilities of incompressible Maclaurin spheroids. We also present several applications of our method to astrophysical problems of great current interest such as the tidal disruption of a star by a massive black hole, the coalescence of compact binaries driven by the emission of gravitational waves, and the development of instabilities in close binary systems.

We have performed a series of three-dimensional hydrodynamic calculations of binary coalescence using the smoothed particle hydrodynamics (SPH) method. The initial conditions are exact polytropic equilibrium configurations with $\gam > 5/3$, on the verge of dynamical instability. We calculate the emission of gravitational radiation in the quadrupole approximation. The fully nonlinear development of the instability is followed until a new equilibrium configuration is reached. We find that the properties of this final configuration depend sensitively on both the compressibility and mass ratio. An {\em axisymmetric} merged configuration is always produced when $\gam\lo2.3$. As a consequence, the emission of gravitational radiation shuts off abruptly right after the onset of dynamical instability. In contrast, {\em triaxial\/} merged configurations are obtained when $\gam\go2.3$, and the system continues to emit gravitational waves after the final coalescence. Systems with mass ratios $q\ne1$ typically become dynamically unstable before the onset of mass transfer. Stable mass transfer from one neutron star to another in a close binary is therefore probably ruled out. The maximum amplitude $h_{max}$ and peak luminosity $L_{max}$ of the gravitational waves emitted during the final coalescence are nearly independent of $\gam$, but depend very sensitively on the mass ratio $q$.

We present a new analytic study of the equilibrium and stability properties of close binary systems containing polytropic components. Our method is based on the use of ellipsoidal trial functions in an energy variational principle. We consider both synchronized and nonsynchronized systems, constructing the compressible generalizations of the classical Darwin and Darwin-Riemann configurations. Our method can be applied to a wide variety of binary models where the stellar masses, radii, spins, entropies, and polytropic indices are all allowed to vary over wide ranges and independently for each component. We find that both secular and dynamical instabilities can develop before a Roche limit or contact is reached along a sequence of models with decreasing binary separation. High incompressibility always makes a given binary system more susceptible to these instabilities, but the dependence on the mass ratio is more complicated. As simple applications, we construct models of double degenerate systems and of low-mass main-sequence-star binaries. We also discuss the orbital evolution of close binary systems under the combined influence of fluid viscosity and secular angular momentum losses from processes like gravitational radiation. We show that the existence of global fluid instabilities can have a profound effect on the terminal evolution of coalescing binaries. The validity of our analytic solutions is examined by means of detailed comparisons with the results of recent numerical fluid calculations in three dimensions.

We study the importance of hydrodynamic effects on the evolution of coalescing binary neutron stars. Using an approximate energy functional constructed from equilibrium solutions for polytropic binary configurations, we incorporate hydrodynamic effects into the calculation of the orbital decay driven by gravitational wave emission. In particular, we follow the transition between the quasi-static, secular decay of the orbit at large separation and the rapid dynamical evolution of configurations approaching contact. We show that a purely Newtonian hydrodynamic instability can significantly accelerate the coalescence at small separation. Such an instability occurs in all close binary configurations containing sufficiently incompressible stars. Calculations are performed for various neutron star masses, radii, spins, and effective polytropic indices. Typically, we find that the radial infall velocity just prior to contact is about 10\% of the tangential orbital velocity. Post-Newtonian effects can move the stability limit to a larger binary separation, and may induce an even larger radial velocity. We also consider the possibility of mass transfer from one neutron star to the other. We show that stable mass transfer is impossible except when the mass of one of the components is very small (less than about 0.4 solar mass) and the viscosity is high enough to maintain corotation.