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.

We present results from a numerical code implementing a new method to solve the master equations describing the evolution of linear perturbations in a spherically symmetric but inhomogeneous background. This method can be used to simulate several configurations of physical interest, such as relativistic corrections to structure formation, the lensing of gravitational waves and the evolution of perturbations in a cosmological void model. This paper focuses on the latter problem, i.e. structure formation in a Hubble scale void in the linear regime. This is considerably more complicated than linear perturbations of a homogeneous and isotropic background because the inhomogeneous background leads to coupling between density perturbations and rotational modes of the spacetime geometry, as well as gravitational waves. Previous analyses of this problem ignored this coupling in the hope that the approximation does not affect the overall dynamics of structure formation in such models. We show that for a giga-parsec void, the evolution of the density contrast is well approximated by the previously studied decoupled evolution only for very large-scale modes. However, the evolution of the gravitational potentials within the void is inaccurate at more than the $10\%$ level, and is even worse on small scales.

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 present a new three-dimensional general-relativistic hydrodynamic evolution scheme coupled to dynamical spacetime evolutions which is capable of efficiently simulating stellar collapse, isolated neutron stars, black hole formation, and binary neutron star coalescence. We make use of a set of adapted curvi-linear grids (multipatches) coupled with flux-conservative cell-centered adaptive mesh refinement. This allows us to significantly enlarge our computational domains while still maintaining high resolution in the gravitational-wave extraction zone, the exterior layers of a star, or the region of mass ejection in merging neutron stars. The fluid is evolved with a high-resolution shock capturing finite volume scheme, while the spacetime geometry is evolved using fourth-order finite differences. We employ a multi-rate Runge-Kutta time integration scheme for efficiency, evolving the fluid with second-order and the spacetime geometry with fourth-order integration, respectively. We validate our code by a number of benchmark problems: a rotating stellar collapse model, an excited neutron star, neutron star collapse to a black hole, and binary neutron star coalescence. The test problems, especially the latter, greatly benefit from higher resolution in the gravitational-wave extraction zone, causally disconnected outer boundaries, and application of Cauchy-characteristic gravitational-wave extraction. We show that we are able to extract convergent gravitational-wave modes up to (l,m)=(6,6). This study paves the way for more realistic and detailed studies of compact objects and stellar collapse in full three dimensions and in large computational domains. The multipatch infrastructure and the improvements to mesh refinement and hydrodynamics codes discussed in this paper will be made available as part of the open-source Einstein Toolkit.

Using accurate numerical relativity simulations of (nonspinning) black-hole binaries with mass ratios 1:1, 2:1 and 3:1 we compute the gauge invariant relation between the (reduced) binding energy $E$ and the (reduced) angular momentum $j$ of the system. We show that the relation $E(j)$ is an accurate diagnostic of the dynamics of a black-hole binary in a highly relativistic regime. By comparing the numerical-relativity $E^{\rm NR} (j)$ curve with the predictions of several analytic approximation schemes, we find that, while the usual, non-resummed post-Newtonian-expanded $E^{\rm PN} (j)$ relation exhibits large and growing deviations from $E^{\rm NR} (j)$, the prediction of the effective one-body formalism, based purely on known analytical results (without any calibration to numerical relativity), agrees strikingly well with the numerical-relativity results.

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.

The accurate modeling of gravitational radiation is a key issue for gravitational wave astronomy. As simulation codes reach higher accuracy, systematic errors inherent in current numerical relativity wave-extraction methods become evident, and may lead to a wrong astrophysical interpretation of the data. In this paper, we give a detailed description of the Cauchy-characteristic extraction technique applied to binary black hole inspiral and merger evolutions to obtain gravitational waveforms that are defined unambiguously, that is, at future null infinity. By this method we remove finite-radius approximations and the need to extrapolate data from the near zone. Further, we demonstrate that the method is free of gauge effects and thus is affected only by numerical error. Various consistency checks reveal that energy and angular momentum are conserved to high precision and agree very well with extrapolated data. In addition, we revisit the computation of the gravitational recoil and find that finite radius extrapolation very well approximates the result at $\scri$. However, the (non-convergent) systematic differences to extrapolated data are of the same order of magnitude as the (convergent) discretisation error of the Cauchy evolution hence highlighting the need for correct wave-extraction.

Binary black-hole systems with spins aligned or anti-aligned to the orbital angular momentum provide the natural ground to start detailed studies of the influence of strong-field spin effects on gravitational wave observations of coalescing binaries. Furthermore, such systems may be the preferred end-state of the inspiral of generic supermassive binary black-hole systems. In view of this, we have computed the inspiral and merger of a large set of binary systems of equal-mass black holes with spins parallel to the orbital angular momentum but otherwise arbitrary. Our attention is particularly focused on the gravitational-wave emission so as to quantify how much spin effects contribute to the signal-to-noise ratio, to the horizon distances, and to the relative event rates for the representative ranges in masses and detectors. As expected, the signal-to-noise ratio increases with the projection of the total black hole spin in the direction of the orbital momentum. We find that equal-spin binaries with maximum spin aligned with the orbital angular momentum are more than "three times as loud" as the corresponding binaries with anti-aligned spins, thus corresponding to event rates up to 30 times larger. We also consider the waveform mismatch between the different spinning configurations and find that, within our numerical accuracy, binaries with opposite spins S_1=-S_2 cannot be distinguished whereas binaries with spin S_1=S_2 have clearly distinct gravitational-wave emissions. Finally, we derive a simple expression for the energy radiated in gravitational waves and find that the binaries always have efficiencies E_rad/M > 3.6%, which can become as large as E_rad/M = 10% for maximally spinning binaries with spins aligned with the orbital angular momentum.

As one step towards a systematic modeling of the electromagnetic (EM) emission from an inspiralling black hole binary we consider a simple scenario in which the binary moves in a uniform magnetic field anchored to a distant circumbinary disc. We study this system by solving the Einstein-Maxwell equations in which the EM fields are chosen with astrophysically consistent strengths. We consider binaries with spins aligned or anti-aligned with the orbital angular momentum and study the dependence of gravitational and EM signals with these spin configurations. Overall we find that the EM radiation in the lowest l=2, m=2 multipole accurately reflects the gravitational one, with identical phase evolutions and amplitudes that differ only by a scaling factor. We also compute the efficiency of the energy emission in EM waves and find that it is given by E^rad_EM/M ~ 10^-15 (M/10^8 M_Sun)^2 (B/10^4 G)^2, hence 13 orders of magnitude smaller than the gravitational energy for realistic magnetic fields. The corresponding luminosity is much smaller than the accretion luminosity if the system is accreting at near the Eddington rate. Most importantly, this EM emission is at frequencies of 10^-4 (10^8 M_Sun/M) Hz, well outside those accessible to astronomical radio observations. As a result, it is unlikely that the EM emission discussed here can be detected directly and simultaneously with the gravitational-wave one. However, indirect processes, driven by changes in the EM fields behavior could yield observable events. In particular if the accretion rate of the circumbinary disc is small and sufficiently stable over the timescale of the final inspiral, then the EM emission may be observable indirectly as it will alter the accretion rate through the magnetic torques exerted by the distorted magnetic field lines.

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 provide a compact analytic formula to compute the spin of the black hole produced by the coalescence of two black holes. The expression, which uses an analytic fit of numerical-relativity data and relies on four assumptions, aims at modelling generic initial spin configurations and mass ratios. A comparison with numerical-relativity simulations already shows very accurate agreements with all of the numerical data available to date, but we also suggest a number of ways in which our predictions can be further improved.

Binary black-hole systems with spins aligned with the orbital angular momentum are of special interest as they may be the preferred end-state of the inspiral of generic supermassive binary black-hole systems. In view of this, we have computed the inspiral and merger of a large set of binary systems of equal-mass black holes with spins aligned with the orbital angular momentum but otherwise arbitrary. By least-square fitting the results of these simulations we have constructed two "spin diagrams" which provide straightforward information about the recoil velocity |v_kick| and the final black-hole spin a_fin in terms of the dimensionless spins a_1 and a_2 of the two initial black holes. Overall they suggest a maximum recoil velocity of |v_kick|=441.94 km/s, and minimum and maximum final spins a_fin=0.3471 and a_fin=0.9591, respectively.

Determining the final spin of a black-hole (BH) binary is a question of key importance in astrophysics. Modelling this quantity in general is made difficult by the fact that it depends on the 7-dimensional space of parameters characterizing the two initial black holes. However, in special cases, when symmetries can be exploited, the description can become simpler. For black-hole binaries with unequal masses but with equal spins which are aligned with the orbital angular momentum, we show that the use of recent simulations and basic but exact constraints derived from the extreme mass-ratio limit allow to model this quantity with a simple analytic expression. Despite the simple dependence, the expression models very accurately all of the available estimates, with errors of a couple of percent at most. We also discuss how to use the fit to predict when a Schwarzschild BH is produced by the merger of two spinning BHs, when the total angular momentum of the spacetime ``flips'' sign, or under what conditions the final BH is ``spun-up'' by the merger. Finally, suggest an extension of the fit to include unequal-spin binaries, thus potentially providing a complete description of the final spin from the coalescence of generic black-hole binaries with spins aligned to the orbital angular momentum.

Binary black-hole systems with spins aligned with the orbital angular momentum are of special interest, as studies indicate that this configuration is preferred in nature. If the spins of the two bodies differ, there can be a prominent beaming of the gravitational radiation during the late plunge, causing a recoil of the final merged black hole. We perform an accurate and systematic study of recoil velocities from a sequence of equal-mass black holes whose spins are aligned with the orbital angular momentum, and whose individual spins range from a = +0.584 to -0.584. In this way we extend and refine the results of a previous study and arrive at a consistent maximum recoil of 448 +- 5 km/s for anti-aligned models as well as to a phenomenological expression for the recoil velocity as a function of spin ratio. This relation highlights a nonlinear behavior, not predicted by the PN estimates, and can be readily employed in astrophysical studies on the evolution of binary black holes in massive galaxies. An essential result of our analysis is the identification of different stages in the waveform, including a transient due to lack of an initial linear momentum in the initial data. Furthermore we are able to identify a pair of terms which are largely responsible for the kick, indicating that an accurate computation can be obtained from modes up to l=3. Finally, we provide accurate measures of the radiated energy and angular momentum, finding these to increase linearly with the spin ratio, and derive simple expressions for the final spin and the radiated angular momentum which can be easily implemented in N-body simulations of compact stellar systems. Our code is calibrated with strict convergence tests and we verify the correctness of our measurements by using multiple independent methods whenever possible.

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.