The grand challenges of contemporary fundamental physics---dark matter, dark energy, vacuum energy, inflation and early universe cosmology, singularities and the hierarchy problem---all involve gravity as a key component. And of all gravitational phenomena, black holes stand out in their elegant simplicity, while harbouring some of the most remarkable predictions of General Relativity: event horizons, singularities and ergoregions. The hitherto invisible landscape of the gravitational Universe is being unveiled before our eyes: the historical direct detection of gravitational waves by the LIGO-Virgo collaboration marks the dawn of a new era of scientific exploration. Gravitational-wave astronomy will allow us to test models of black hole formation, growth and evolution, as well as models of gravitational-wave generation and propagation. It will provide evidence for event horizons and ergoregions, test the theory of General Relativity itself, and may reveal the existence of new fundamental fields. The synthesis of these results has the potential to radically reshape our understanding of the cosmos and of the laws of Nature. The purpose of this work is to present a concise, yet comprehensive overview of the state of the art in the relevant fields of research, summarize important open problems, and lay out a roadmap for future progress.

The recent detection of gravitational waves and electromagnetic counterparts emitted during and after the collision of two neutron stars marks a breakthrough in the field of multi-messenger astronomy. Numerical relativity simulations are the only tool to describe the binary's merger dynamics in the regime when speeds are largest and gravity is strongest. In this work we report state-of-the-art binary neutron star simulations for irrotational (non-spinning) and spinning configurations. The main use of these simulations is to model the gravitational-wave signal. Key numerical requirements are the understanding of the convergence properties of the numerical data and a detailed error budget. The simulations have been performed on different HPC clusters, they use multiple grid resolutions, and are based on eccentricity reduced quasi-circular initial data. We obtain convergent waveforms with phase errors of 0.5-1.5 rad accumulated over approximately 12 orbits to merger. The waveforms have been used for the construction of a phenomenological waveform model which has been applied for the analysis of the recent binary neutron star detection. Additionally, we show that the data can also be used to test other state-of-the-art semi-analytical waveform models.

We present the first set of numerical relativity simulations of binary neutron mergers that include spin precession effects and are evolved with multiple resolutions. Our simulations employ consistent initial data in general relativity with different spin configurations and dimensionless spin magnitudes $\sim 0.1$. They start at a gravitational-wave frequency of $\sim392$~Hz and cover more than $1$ precession period and about 15 orbits up to merger. We discuss the spin precession dynamics by analyzing coordinate trajectories, quasi-local spin measurements, and energetics, by comparing spin aligned, antialigned, and irrotational configurations. Gravitational waveforms from different spin configuration are compared by calculating the mismatch between pairs of waveforms in the late inspiral. We find that precession effects are not distinguishable from nonprecessing configurations with aligned spins for approximately face-on binaries, while the latter are distinguishable from a nonspinning configurations. Spin precession effects are instead clearly visible for approximately edge-on binaries. For the parameters considered here, precession does not significantly affect the characteristic postmerger gravitational-wave frequencies nor the mass ejection. Our results pave the way for the modeling of spin precession effects in the gravitational waveform from binary neutron star events.

We present new (3+1)D numerical relativity simulations of the binary neutron star (BNS) merger and postmerger phase. We focus on a previously inaccessible region of the binary parameter space spanning the binary's mass-ratio $q\sim1.00-1.75$ for different total masses and equations of state, and up to $q\sim2$ for a stiff BNS system. We study the mass-ratio effect on the gravitational waves (GWs) and on the possible electromagnetic emission associated to dynamical mass ejecta. We compute waveforms, spectra, and spectrograms of the GW strain including all the multipoles up to $l=4$. The mass-ratio has a specific imprint on the GW multipoles in the late-inspiral-merger signal, and it affects qualitatively the spectra of the merger remnant. The multipole effect is also studied by considering the dependency of the GW spectrograms on the source's sky location. Unequal mass BNSs produce more ejecta than equal mass systems with ejecta masses and kinetic energies depending almost linearly on $q$. We estimate luminosity peaks and light curves of macronovae events associated to the mergers using a simple approach. For $q\sim2$ the luminosity peak is delayed for several days and can be up to four times larger than for the $q=1$ cases. The macronova emission associated with the $q\sim2$ BNS is more persistent in time and could be observed for weeks instead of few days ($q=1$) in the near infrared. Finally, we estimate the flux of possible radio flares produced by the interaction of relativistic outflows with the surrounding medium. Also in this case a large $q$ can significantly enhance the emission and delay the peak luminosity. Overall, our results indicate that BNS merger with large mass ratio have EM signatures distinct from the equal mass case and more similar to black hole - neutron star binaries.

Information about the last stages of a binary neutron star inspiral and the final merger can be extracted from quasi-equilibrium configurations and dynamical evolutions. In this article, we construct quasi-equilibrium configurations for different spins, eccentricities, mass ratios, compactnesses, and equations of state. For this purpose we employ the SGRID code, which allows us to construct such data in previously inaccessible regions of the parameter space. In particular, we consider spinning neutron stars in isolation and in binary systems; we incorporate new methods to produce highly eccentric and eccentricity reduced data; we present the possibility of computing data for significantly unequal-mass binaries; and we create equal-mass binaries with individual compactness up to 0.23. As a proof of principle, we explore the dynamical evolution of three new configurations. First, we simulate a $q=2.06$ mass ratio which is the highest mass ratio for a binary neutron star evolved in numerical relativity to date. We find that mass transfer from the companion star sets in a few revolutions before merger and a rest mass of $\sim10^{-2}M_\odot$ is transferred between the two stars. This configuration also ejects a large amount of material during merger, imparting a substantial kick to the remnant. Second, we simulate the first merger of a precessing binary neutron star. We present the dominant modes of the gravitational waves for the precessing simulation, where a clear imprint of the precession is visible in the (2,1) mode. Finally, we quantify the effect of an eccentricity reduction procedure on the gravitational waveform. The procedure improves the waveform quality and should be employed in future precision studies, but also other errors, notably truncation errors, need to be reduced in order for the improvement due to eccentricity reduction to be effective. [abridged]

We study equal and unequal-mass neutron star mergers by means of new numerical relativity simulations in which the general relativistic hydrodynamics solver employs an algorithm that guarantees mass conservation across the refinement levels of the computational mesh. We consider eight binary configurations with total mass $M=2.7\,M_\odot$, mass-ratios $q=1$ and $q=1.16$, and four different equation of states (EOSs), and one configuration with a stiff EOS, $M=2.5M_\odot$ and $q=1.5$. We focus on the post-merger dynamics and study the merger remnant, dynamical ejecta and the postmerger gravitational wave spectrum. Although most of the merger remnants form a hypermassive neutron star collapsing to a black hole+disk system on dynamical timescales, stiff EOSs can eventually produce a stable massive neutron star. Ejecta are mostly emitted around the orbital plane; favored by large mass ratios and softer EOS. The postmerger wave spectrum is mainly characterized by non-axisymmetric oscillations of the remnant. The stiff EOS configuration consisting of a $1.5M_\odot$ and a $1.0M_\odot$ neutron star shows a rather peculiar dynamics. During merger the companion star is very deformed; about~$\sim0.03M_\odot$ of rest-mass becomes unbound from the tidal tail due torque; and the merger remnant forms stable neutron star surrounded by a massive accretion disk $\sim0.3M_\odot$. Similar configurations might be particularly interesting for electromagnetic counterparts. Comparing results obtained with and without the conservative mesh refinement algorithm, we find that post-merger simulations can be affected by systematic errors if mass conservation is not enforced in the mesh refinement strategy. However, mass conservation also depends on grid details and on the artificial atmosphere setup. [abridged]

Binary neutron stars in circular orbits can be modeled as helically symmetric, i.e., stationary in a rotating frame. This symmetry gives rise to a first integral of the Euler equation, often employed for constructing equilibrium solutions via iteration. For eccentric orbits, however, the lack of helical symmetry has prevented the use of this method, and the numerical relativity community has often resorted to constructing initial data by superimposing boosted spherical stars without solving the Euler equation. The spuriously excited neutron star oscillations seen in evolutions of such data arise because such configurations lack the appropriate tidal deformations and are stationary in a linearly comoving---rather than rotating---frame. We consider eccentric configurations at apoapsis that are instantaneously stationary in a rotating frame. We extend the notion of helical symmetry to eccentric orbits, by approximating the elliptical orbit of each companion as instantaneously circular, using the ellipse's inscribed circle. The two inscribed helical symmetry vectors give rise to approximate instantaneous first integrals of the Euler equation throughout each companion. We use these integrals as the basis of a self-consistent iteration of the Einstein constraints to construct conformal thin-sandwich initial data for eccentric binaries. We find that the spurious stellar oscillations are reduced by at least an order of magnitude, compared with those found in evolutions of superposed initial data. The tidally induced oscillations, however, are physical and qualitatively similar to earlier evolutions. Finally, we show how to incorporate radial velocity due to radiation reaction in our inscribed helical symmetry vectors, which would allow one to obtain truly non-eccentric initial data when our eccentricity parameter $e$ is set to zero.

Simulations of binary neutron stars have seen great advances in terms of physical detail and numerical quality. However, the spin of the neutron stars, one of the simplest global parameters of binaries, remains mostly unstudied. We present the first, fully nonlinear general relativistic dynamical evolutions of the last three orbits for constraint satisfying initial data of spinning neutron star binaries, with astrophysically realistic spins aligned and anti-aligned to the orbital angular momentum. The initial data is computed with the constant rotational velocity approach. The dynamics of the systems is analyzed in terms of gauge-invariant binding energy vs. orbital angular momentum curves. By comparing to a binary black hole configuration we can estimate the different tidal and spin contributions to the binding energy for the first time. First results on the gravitational wave forms are presented. The phase evolution during the orbital motion is significantly affected by spin-orbit interactions, leading to delayed or early mergers. Furthermore, a frequency shift in the main emission mode of the hyper massive neutron star is observed. Our results suggest that a detailed modeling of merger waveforms requires the inclusion of spin, even for the moderate magnitudes observed in binary neutron star systems.

Neutron star binaries offer a rich phenomenology in terms of gravitational waves and merger remnants. However, most general relativistic studies have been performed for nearly circular binaries, with the exception of head-on collisions. We present the first numerical relativity investigation of mergers of eccentric equal-mass neutron-star binaries that probes the regime between head-on and circular. In addition to gravitational waves generated by the orbital motion, we find that the signal also contains a strong component due to stellar oscillations (f-modes) induced by tidal forces, extending a classical result for Newtonian binaries. The merger can lead to rather massive disks on the order of 10% of the total initial mass.

We perform a parameter study of non-spinning, equal and unequal mass black hole binaries on generic, eccentric orbits in numerical relativity. The linear momentum considered ranges from that of a circular orbit to ten times that value. We discuss the different manifestations of zoom-whirl behavior in the hyperbolic and the elliptic regime. The hyperbolic data set applies to dynamical capture scenarios (e.g. in globular clusters). Evolutions in the elliptic regime correspond to possible end states of supermassive black hole binaries. We spot zoom-whirl behavior for eccentricities as low as $e\sim0.5$, i.e. within the expected range of eccentricities in massive black hole binaries from galaxy mergers and binaries near galactic centers. The resulting gravitational waveforms reveal a rich structure, which will effectively break degeneracies in parameter space improving parameter estimation.

We present numerical relativity simulations of nine-orbit equal-mass binary neutron star covering the quasicircular late inspiral and merger. The extracted gravitational waveforms are analyzed for convergence and accuracy. Second order convergence is observed up to contact, i.e. about 3-4 cycles to merger; error estimates can be made up to this point. The uncertainties on the phase and the amplitude are dominated by truncation errors and can be minimized to 0.13 rad and less then 1%, respectively, by using several simulations and extrapolating in resolution. In the latter case finite-radius extraction uncertainties become a source of error of the same order and have to be taken into account. The waveforms are tested against accuracy standards for data analysis. The uncertainties on the waveforms are such that accuracy standards are generically not met for signal-to-noise ratios relevant for detection, except for some best cases using extrapolation from several runs. A detailed analysis of the errors is thus imperative for the use of numerical relativity waveforms from binary neutron stars in quantitative studies. The waveforms are compared with the post-Newtonian Taylor T4 approximants both for point-particle and including the analytically known tidal corrections. The T4 approximants accumulate significant phase differences of 2 rad at contact and 4 rad at merger, underestimating the influence of finite size effects. Tidal signatures in the waveforms are thus important at least during the last six orbits of the merger process.

We compare dynamics and waveforms from binary neutron star coalescence as computed by new long-term ($\sim 10 $ orbits) numerical relativity simulations and by the tidal effective-one-body (EOB) model including analytical tidal corrections up to second post-Newtonian order (2PN). The current analytical knowledge encoded in the tidal EOB model is found to be sufficient to reproduce the numerical data up to contact and within their uncertainties. Remarkably, no calibration of any tidal EOB free parameters is required, beside those already fitted to binary black holes data. The inclusion of 2PN tidal corrections minimizes the differences with the numerical data, but it is not possible to significantly distinguish them from the leading-order tidal contribution. The presence of a relevant amplification of tidal effects is likely to be excluded, although it can appear as a consequence of numerical inaccuracies. We conclude that the tidally-completed effective-one-body model provides nowadays the most advanced and accurate tool for modelling gravitational waveforms from binary neutron star inspiral up to contact. This work also points out the importance of extensive tests to assess the uncertainties of the numerical data, and the potential need of new numerical strategies to perform accurate simulations.

We present a new numerical relativity code designed for simulations of compact binaries involving matter. The code is an upgrade of the BAM code to include general relativistic hydrodynamics and implements state-of-the-art high-resolution-shock-capturing schemes on a hierarchy of mesh refined Cartesian grids with moving boxes. We test and validate the code in a series of standard experiments involving single neutron star spacetimes. We present test evolutions of quasi-equilibrium equal-mass irrotational binary neutron star configurations in quasi-circular orbits which describe the late inspiral to merger phases. Neutron star matter is modeled as a zero-temperature fluid; thermal effects can be included by means of a simple ideal-gas prescription. We analyze the impact that the use of different values of damping parameter in the Gamma-driver shift condition has on the dynamics of the system. The use of different reconstruction schemes and their impact in the post-merger dynamics is investigated. We compute and characterize the gravitational radiation emitted by the system. Self-convergence of the waves is tested, and we consistently estimate error-bars on the numerically generated waveforms in the inspiral phase.

We present the first analytical inspiral-merger-ringdown gravitational waveforms from binary black holes (BBHs) with non-precessing spins, that is based on a description of the late-inspiral, merger and ringdown in full general relativity. By matching a post-Newtonian description of the inspiral to a set of numerical-relativity simulations, we obtain a waveform family with a conveniently small number of physical parameters. These waveforms will allow us to detect a larger parameter space of BBH coalescence, including a considerable fraction of precessing binaries in the comparable-mass regime, thus significantly improving the expected detection rates.

Gravitational waveforms from the inspiral and ring-down stages of the binary black hole coalescences can be modelled accurately by approximation/perturbation techniques in general relativity. Recent progress in numerical relativity has enabled us to model also the non-perturbative merger phase of the binary black-hole coalescence problem. This enables us to \emph{coherently} search for all three stages of the coalescence of non-spinning binary black holes using a single template bank. Taking our motivation from these results, we propose a family of template waveforms which can model the inspiral, merger, and ring-down stages of the coalescence of non-spinning binary black holes that follow quasi-circular inspiral. This two-dimensional template family is explicitly parametrized by the physical parameters of the binary. We show that the template family is not only \emph{effectual} in detecting the signals from black hole coalescences, but also \emph{faithful} in estimating the parameters of the binary. We compare the sensitivity of a search (in the context of different ground-based interferometers) using all three stages of the black hole coalescence with other template-based searches which look for individual stages separately. We find that the proposed search is significantly more sensitive than other template-based searches for a substantial mass-range, potentially bringing about remarkable improvement in the event-rate of ground-based interferometers. As part of this work, we also prescribe a general procedure to construct interpolated template banks using non-spinning black hole waveforms produced by numerical relativity.

We study the transition from inspiral to plunge in general relativity by computing gravitational waveforms of non-spinning, equal-mass black-hole binaries. We consider three sequences of simulations, starting with a quasi-circular inspiral completing 1.5, 2.3 and 9.6 orbits, respectively, prior to coalescence of the holes. For each sequence, the binding energy of the system is kept constant and the orbital angular momentum is progressively reduced, producing orbits of increasing eccentricity and eventually a head-on collision. We analyze in detail the radiation of energy and angular momentum in gravitational waves, the contribution of different multipolar components and the final spin of the remnant. We find that the motion transitions from inspiral to plunge when the orbital angular momentum L=L_crit is about 0.8M^2. For L<L_crit the radiated energy drops very rapidly. Orbits with L of about L_crit produce our largest dimensionless Kerr parameter for the remnant, j=J/M^2=0.724. Generalizing a model recently proposed by Buonanno, Kidder and Lehner to eccentric binaries, we conjecture that (1) j=0.724 is the maximal Kerr parameter that can be obtained by any merger of non-spinning holes, and (2) no binary merger (even if the binary members are extremal Kerr black holes with spins aligned to the orbital angular momentum, and the inspiral is highly eccentric) can violate the cosmic censorship conjecture.

We study the inspiral, merger and ringdown of unequal mass black hole binaries by analyzing a catalogue of numerical simulations for seven different values of the mass ratio (from q=M2/M1=1 to q=4). We compare numerical and Post-Newtonian results by projecting the waveforms onto spin-weighted spherical harmonics, characterized by angular indices (l,m). We find that the Post-Newtonian equations predict remarkably well the relation between the wave amplitude and the orbital frequency for each (l,m), and that the convergence of the Post-Newtonian series to the numerical results is non-monotonic. To leading order the total energy emitted in the merger phase scales like eta^2 and the spin of the final black hole scales like eta, where eta=q/(1+q)^2 is the symmetric mass ratio. We study the multipolar distribution of the radiation, finding that odd-l multipoles are suppressed in the equal mass limit. Higher multipoles carry a larger fraction of the total energy as q increases. We introduce and compare three different definitions for the ringdown starting time. Applying linear estimation methods (the so-called Prony methods) to the ringdown phase, we find resolution-dependent time variations in the fitted parameters of the final black hole. By cross-correlating information from different multipoles we show that ringdown fits can be used to obtain precise estimates of the mass and spin of the final black hole, which are in remarkable agreement with energy and angular momentum balance calculations.

Recent progress in numerical relativity has enabled us to model the non-perturbative merger phase of the binary black-hole coalescence problem. Based on these results, we propose a phenomenological family of waveforms which can model the inspiral, merger, and ring-down stages of black hole coalescence. We also construct a template bank using this family of waveforms and discuss its implementation in the search for signatures of gravitational waves produced by black-hole coalescences in the data of ground-based interferometers. This template bank might enable us to extend the present inspiral searches to higher-mass binary black-hole systems, i.e., systems with total mass greater than about 80 solar masses, thereby increasing the reach of the current generation of ground-based detectors.

When unequal-mass black holes merge, the final black hole receives a ``kick'' due to the asymmetric loss of linear momentum in the gravitational radiation emitted during the merger. The magnitude of this kick has important astrophysical consequences. Recent breakthroughs in numerical relativity allow us to perform the largest parameter study undertaken to date in numerical simulations of binary black hole inspirals. We study non-spinning black-hole binaries with mass ratios from $q=M_1/M_2=1$ to $q =0.25$ ($\eta = q/(1 + q)^2$ from 0.25 to 0.16). We accurately calculate the velocity of the kick to within 6%, and the final spin of the black holes to within 2%. A maximum kick of $175.2\pm11$ km s$^{-1}$ is achieved for $\eta = 0.195 \pm 0.005$.

We study the coalescence of non-spinning binary black holes from near the innermost stable circular orbit down to the final single rotating black hole. We use a technique that combines the full numerical approach to solve Einstein equations, applied in the truly non-linear regime, and linearized perturbation theory around the final distorted single black hole at later times. We compute the plunge waveforms which present a non negligible signal lasting for $t\sim 100M$ showing early non-linear ringing, and we obtain estimates for the total gravitational energy and angular momentum radiated.

We study the stability of three-dimensional numerical evolutions of the Einstein equations, comparing the standard ADM formulation to variations on a family of formulations that separate out the conformal and traceless parts of the system. We develop an implementation of the conformal-traceless (CT) approach that has improved stability properties in evolving weak and strong gravitational fields, and for both vacuum and spacetimes with active coupling to matter sources. Cases studied include weak and strong gravitational wave packets, black holes, boson stars and neutron stars. We show under what conditions the CT approach gives better results in 3D numerical evolutions compared to the ADM formulation. In particular, we show that our implementation of the CT approach gives more long term stable evolutions than ADM in all the cases studied, but is less accurate in the short term for the range of resolutions used in our 3D simulations.

We present the first results in a new program intended to make the best use of all available technologies to provide an effective understanding of waves from inspiralling black hole binaries in time for imminent observations. In particular, we address the problem of combining the close-limit approximation describing ringing black holes and full numerical relativity, required for essentially nonlinear interactions. We demonstrate the effectiveness of our approach using general methods for a model problem, the head-on collision of black holes. Our method allows a more direct physical understanding of these collisions indicating clearly when non-linear methods are important. The success of this method supports our expectation that this unified approach will be able to provide astrophysically relevant results for black hole binaries in time to assist gravitational wave observations.

The standard approach to the numerical evolution of black hole data using the ADM formulation with maximal slicing and vanishing shift is extended to non-symmetric black hole data containing black holes with linear momentum and spin by using a time-independent conformal rescaling based on the puncture representation of the black holes. We give an example for a concrete three dimensional numerical implementation. The main result of the simulations is that this approach allows for the first time to evolve through a brief period of the merger phase of the black hole inspiral.