Black Hole Leftovers: The Remnant Population from Binary Black Hole Mergers

The LVC has employed a suite of population models to describe the GWTC-2 dataset. Herein, we use the LVC results from the population model with the highest Bayes factor, which is the Power-law + Peak model (Abbott et al., 2020b). In this model, the population of primary masses (the more massive component in each binary) is parameterized as a mixture model of a power law plus Gaussian (Talbot & Thrane, 2017). The secondary mass is a power law conditioned on the mass of the primary. The spin magnitude and tilt distributions (at a reference GW frequency of 20 Hz²²2GWTC-2 parameter estimates reference all spins to 20 Hz except for GW190521, which was referenced to 11 Hz due to its large mass. Nevertheless, we treat the GW190521 spins as if they were referenced to 20 Hz. Abbott et al. (2020b) has verified that this treatment results in errors smaller than the error between results with different gravitational waveform models.) are assumed to be independent of one another and the masses. The spin magnitudes are described by a Beta distribution, and the cosines of the tilt angles with respect to the orbital angular momentum are parameterized with a mixture of a Gaussian and uniform distribution. Since the distribution of BH spin azimuthal angles (i.e. spin direction in the orbital plane) are not expressly fit by the LVC population analyses, they inherit the prior used in LVC parameter estimation, namely isotropic in the orbital plane. The full details of the model and priors can be found in Abbott et al. (2020b). The hyperparameters (e.g. power law indices, Gaussian means and variances) that describe these distributions are constrained by the GWTC-2 data, and inferences on them are reported as “hyperposterior” samples by the LVC. Each hyperposterior sample $\Lambda_{j}$ implies a population distribution of inspiral parameters $p(\theta|\Lambda_{j})$, which can in turn be converted to remnant parameters.

To incorporate these hyperposterior inferences, we start by generating a sample of inspiral masses and spins $\theta_{i}$ from a reference distribution $p(\theta|\Lambda_{\rm ref})$, where $\Lambda_{\rm ref}$ is the reference population model. For each inspiral sample $\theta_{i}$ and each hyperposterior sample $\Lambda_{j}$, we calculate a weight $w_{ij}=p(\theta_{i}|\Lambda_{j})/p(\theta_{i}|\Lambda_{\rm ref})$. We apply these weights to the remnant parameter samples $\theta_{i}^{\rm rem}=f_{\rm GR}(\theta_{i})$, which are calculated through a function $f_{\rm GR}(\theta)$ that is specified by GR. These weights along with the remnant samples $\theta_{i}^{\rm rem}$ represent draws from the remnant distributions corresponding to each $\Lambda_{j}$. Finally, to estimate and visualize the distributions of the remnants, we apply Gaussian kernel density estimates to the samples and weights, yielding continuous functions over the remnant parameter space.

The mass, spin, and kick velocity of a BBH remnant are fully specified by the properties and configuration of the inspiraling BHs, and can be calculated in GR. In practice, numerical relativity (NR) simulations are needed to compute the remnant parameters for an arbitrary set of inspiral masses and spins. Since these simulations can be computationally expensive, fast surrogate models have been built to interpolate remnant parameters between NR simulations. We use two such surrogates here, namely the NRSur3dq8Remnant and NRSur7dq4Remnant models from the surfinBH package (Varma et al., 2019b, a). We default to the NRSur7dq4Remnant model which calculates the final mass $M_{f}$, final spin vector $\vec{\chi}_{f}$, and kick velocity vector $\vec{v}_{f}$ given the 3D dimensionless spin vectors of the two inspiraling BHs and their mass ratio. The model is trained for mass ratios between 1 and 4 and can reliably extrapolate out to mass ratios of 6. For mergers with $q>6$, we switch to the NRSur3dq8Remnant model that is trained out to $q=8$ but assumes the BH spins are aligned with the orbital angular momentum. In these $q>6$ cases, we “ignore” the in-plane spin components by setting them to zero. While some systematic error is picked up from reverting to the NRSur3dq8Remnant model, we note that $q>6$ mergers are expected to be only a fraction of a percent of all mergers under the Power Law + Peak population.

To show how the BH remnant parameters change with binary inspiral parameters, we display the remnant properties calculated via NRSur7dq4Remnant for a set of example cases in Figure 1. Each line shows the remnant parameters as a function of mass ratio for a fixed spin configuration (defined at $t=100M$ before peak of the GW). In the top panel, we see that the ratio of final mass $M_{f}$ to initial mass $M$ follows a common trend regardless of spin configuration: it is greater than 0.9 and approaches 1 as the mass ratio increases. On the other hand, the final spin and final velocity are quite sensitive to spins and mass ratio (Fishbach et al., 2017). Although $\chi_{f}\sim 0.7$ at $q=1$ for most cases, there can be significant deviations from that depending on the spins and mass ratio. For example, when the inspiral spins are large and anti-aligned with the orbital angular momentum (green, dash-dotted), the final spin decreases to near zero with increasing mass ratio until $q=4$ and then increases. The final velocity is the most sensitive to input configuration, though the kick velocities tend to be highest with large, misaligned component spin magnitudes.

One other feature of these surrogate models is that they are capable of reporting 1-$\sigma$ uncertainties on remnant parameters due to interpolation error (Varma et al., 2019b, a). To account for the surrogate model uncertainty in our population analysis, we dither each calculated remnant parameter using the surrogates’ reported 1-$\sigma$ intervals.

In addition to the shape of the population distributions of mergers, the LVC has inferred the total number of mergers occurring per comoving time and volume in the local universe (Abbott et al., 2020b). Integrating this rate density from the beginning of the universe to now yields the local number density of remnants. One complication in this number density calculation is that the rate density of mergers is likely not constant over cosmic time, because the star formation rate evolves over time, and BBH mergers are expected to occur with some delay from the initial formation of the progenitor stars. To incorporate these considerations, we use a procedure analogous to that in Abbott et al. (2017). We convolve the star formation rate from Madau & Fragos (2017) with a $t^{-1}$ distribution of delay times from progenitor formation to BBH merger and a minimum delay time of $t_{\rm min}=10$ Myr. We then assume this convolution is proportional to the merger rate density $R$. That is, $R(t)=A\int_{t_{\rm min}}^{t_{h}}p(t-t^{\prime})\psi(t^{\prime})dt^{\prime}$, where $A$ is a normalization constant, $p(t-t^{\prime})=(t-t^{\prime})^{-1}$, $t_{h}$ is the Hubble time, and $\psi(t)$ is the cosmic star formation rate. Integrating this over cosmic time yields

where $R_{0}$ is the rate density today, which we take to be the local rate inferred by the LVC. For each hyperposterior sample $R_{0}$ we compute $n$ to build a posterior distribution on the number density.

As mentioned in Section 1, there are two primary ways one might detect a BBH remnant: through gravitational waves emitted by hierarchical mergers or though “indirect” electromagnetic detections of BHs. The existence of hierarchical mergers is still uncertain, but the eventual reach of future gravitational-wave detectors could probe these if they are occurring with an appreciable rate (Hall & Evans, 2019; Ng et al., 2020). If hierarchical mergers are rare, a BBH detected with a component coming exclusively from the distributions in Figure 2 could be a smoking gun of a hierarchical merger with a BBH remnant.

On the other hand, if hierarchical mergers are common, then the GWTC-2 data set and LVC population inferences may already include hierarchical mergers. In that case, our results cannot be interpreted as present-day number densities of remnants, because we do not account for remnants merging again. Instead, our resultant number density distributions could be re-normalized to the inferred merger rates and interpreted as the distributions over parameters of remnants that are being created per unit time-volume. Such a re-normalization would make our results directly analogous to the inferences of the inspiraling population: they are inferences on numbers of events per volume and time, not a number density. To properly compute the present day number density of remnants in this strongly hierarchical scenario, we would need to include second generation mergers through a self-consistent hierarchical population model (Kimball et al., 2020a; Doctor et al., 2020). There are tantalizing hints that some events in GWTC-2 may have a hierarchical origin (Kimball et al., 2020b; Zevin et al., 2020), so future work could include this possibility in the number density calculation.

Now we turn to electromagnetic means of detecting BBH remnants. As with hierarchical mergers, the detectability of electromagnetic signals from BBH remnants is contingent on the number of remnants, the remnants’ likelihood of interacting with other bodies, and the statistical power to differentiate a remnant signal from a non-remnant BH signal background. With a local remnant number density of $660_{-240}^{+440}\,\mathrm{Mpc}^{-3}$, one could expect $O(60,000)$ remnants produced in the Milky Way³³3This assumes a conversion factor of $\sim$ 0.01 Milky-Way-equivalent galaxies per cubic Mpc, which is similar to the conversion used in Abadie et al. (2010)., which matches theoretical expectations (e.g. Lamberts et al., 2018; Olejak et al., 2020). This number of remnants, even ignoring that some may be ejected from the galaxy due to kicks, is far below the limits set on primordial black hole number densities by microlensing studies, which probe populations of black holes comprising percents of the total dark matter mass (Lu et al., 2019; Tisserand et al., 2007; García-Bellido, 2017). As such, we do not expect that these remnants can be detected with current instruments with microlensing. Another avenue for detection might be through observations of dynamics of companion stars or X-ray emission from accretion onto a remnant. These detections are possible only under the highly speculative condition that the remnants encounter other bodies after their creation, which could occur in a cluster environment or if the remnant approaches the center of the galaxy due to dynamical friction. The results of our study could serve as a guide for extracting whether BHs detected via electromagnetic means come from a remnant or stellar population.

In this work we have made a number of assumptions which could affect the resultant remnant parameter distributions and number densities. Firstly, we only consider the results of the LVC’s inferences with the Power Law + Peak model. If the inspiral distribution has more complicated features that are not resolved by this model, the resultant remnant parameter distributions may also pick up additional features. However, the LVC’s analyses with other models show a consistent picture of a lower rate of mergers with larger mass, meaning the overall shape of the remnant mass distributions in Figure 2 should be relatively insensitive to other model choices. Likewise, the final spins of mergers tend to be near $\chi_{f}=0.7$ except in select cases, so the spin distribution is also unlikely to change dramatically with different parameterizations. On the other hand, the kick velocities are much more sensitive to the assumed models. Changes to the parameterization of inspiral mass ratio or spins could affect the kick velocities. For example, if the component BH spin angles in the orbital plane do not match our assumption of being uniformly distributed, which could happen in some dynamical environments (Mould & Gerosa, 2020; Yu et al., 2020), certain kick speeds could become (dis-)favored. Changes in spin orientations alone can change the kick velocities significantly, as seen in Figure 1.

In calculating the number density, we assume that the rate of black hole mergers follows the star formation rate of Madau & Fragos (2017) with a $t^{-1}$ delay time and minimum delay time of 10 Myr. Since the star formation history and black hole merger delay times are still relatively uncertain, the overall number density of remnants calculated herein is subject to change with stronger constraints, but unlikely to change by the multiple orders of magnitude needed to affect our conclusions about remnant detectability, assuming the binaries come from a stellar population. For example, if we use the rate density versus redshift profile from Rodriguez & Loeb (2018) for BBHs from globular clusters and normalize to the inferred LVC rate at $z=0$, the inferred present-day number density changes by $<10\%$. We neglect the possibility that the binaries come from primordial black holes, which would entail a different redshift evolution of the merger rate.