Identifying the hosts of binary black hole and neutron star-black hole mergers
with next-generation gravitational-wave detectors

The LIGO-Virgo-KAGRA (LVK) Collaboration’s third observing run has increased the number of compact binary mergers detected in gravitational waves to almost one hundred (Abbott et al., 2021). Robust population analyses of these events have made great strides in understanding various aspects of the overall population, including the black hole mass function, merger rate, spin population, and evolution with redshift (Abbott et al., 2023). However, the formation channels for these mergers remain poorly understood. Plausible theories include isolated evolution of main sequence stars in binaries (e.g., de Mink & Mandel (2016); Woosley (2016); Belczynski et al. (2016); Stevenson et al. (2017)), dynamical exchanges in dense star clusters (e.g., Rodriguez et al. (2016); Askar et al. (2017); Di Carlo et al. (2019); Gerosa & Fishbach (2021)), black hole encounters in gaseous environments such as active galactic nuclei (AGN) (e.g., Bartos et al. (2017); Secunda et al. (2019); Yang et al. (2019); Tagawa et al. (2020)), among others (see Mapelli (2021) for a thorough review).

Associating a merger with its host galaxy can help build a much more complete picture of its natal state and evolution, aiding our understanding of how these compact binaries were formed. Only one GW event, the multi-messenger binary neutron star (BNS) merger GW170817, has been associated with its host galaxy (Abbott et al., 2017a, b). While BNS mergers like GW170817 are expected to produce electromagnetic (EM) waves along with GWs, allowing for relatively straightforward host identification, all binary black hole (BBH) and most neutron star–black hole (NSBH) mergers are likely to be EM-dark (Foucart, 2020; Biscoveanu et al., 2022; Chen et al., 2021; Fragione, 2021). With current-generation GW detectors, pinning down a host galaxy for a merger without an EM counterpart is nonviable: at design sensitivity with a four-detector network, the best localized BBHs will have 90% localization volumes of $\mathcal{O}(10^{3})$ Mpc³ (Petrov et al., 2021), insufficient to associate a merger with its host.

xg GW detectors such as Cosmic Explorer (CE) (Evans et al., 2021) and Einstein Telescope (ET) (Punturo et al., 2010) promise to unveil almost every BBH merger out to redshifts greater than 10 (Vitale & Evans, 2017; Regimbau et al., 2017; Chen et al., 2021). This data will be immensely powerful for population studies, but there will also be events at low redshifts which are detected at such high signal-to-noise ratios that they will be localized to miniscule volumes, allowing for identification of the host galaxy with GWs alone (Chen & Holz, 2016; Vitale & Whittle, 2018; Borhanian et al., 2020). We show in this work that these next generation (XG) detectors will prove to be powerful engines for finding host galaxies of compact binary coalescences such as BBH mergers.

The host of a CBC provides rich information about how the progenitors came to be and how they evolved. For example, some theories predict that the dense, gaseous environments of AGN disks can catalyze the formation and inspiral of BBHs (Bartos et al., 2017; Secunda et al., 2019; Yang et al., 2019). In some cases, such a BBH merger may produce an EM counterpart (McKernan et al., 2019; Graham et al., 2020). However, it can be difficult to disentangle such a counterpart with usual AGN variability or flares, making it challenging to verify these predictions (Palmese et al., 2021). \Acgw-only localizations of BBH mergers to AGN hosts would provide strong evidence for the existence of this channel.

Even for individual hosts not as unique as AGN, building up demographics of hosts can be powerful, providing insight into characteristics such as late-type vs early-type galaxies, host stellar mass, and star formation rate (SFR). One useful parameterization connecting these host observables to binary formation and evolution is the delay time distribution (DTD), which describes the time between the zero-age main sequence formation of the stars in the binary and their eventual merger as compact objects.

The connection between host observables and the DTD can be seen in e.g., the SFR of hosts. If most hosts are found to have low SFR, the DTD should have more support for longer delay times, since the star formation which produced the progenitor stars could have been quenched by the time the binary’s merger occurs. On the other hand, star forming hosts lend support to shorter delay times, since the host galaxy will not have evolved substantially between the formation of the individual stars in the binary and their merger.

Since the lifetime of the massive stars that result in GW-detectable binaries is short ($\lesssim$ 10 Myr), the DTD of CBCs essentially depends on the separation of the binary at the start of the inspiral, and the duration of the subsequent inspiral. The DTD is typically parameterized by a power law, $dN/dt\propto t^{\alpha}$, where $N$ is the number of mergers, $t$ is the time, and $\alpha$ is the power law index. This parameterization is based on the distribution of initial separations of massive binaries and inspiral physics from general relativity, and has also been confirmed by simulations (Peters, 1964; Piran, 1992; Zevin et al., 2022). A minimum delay time $t_{\rm min}$ is often used to account for the lifetime of the massive progenitors and the time required to form the resulting compact object binary (Safarzadeh & Berger, 2019; Zevin et al., 2022).

A number of studies have found that the characterization of CBC hosts can lead to robust constraints on the DTD. Adhikari et al. (2020) show from galaxy simulations that $\mathcal{O}(100)$ events with hosts will constrain the DTD, parameterized as $dN/dt\propto(t-t_{\rm min})^{\alpha}$, to Gyr precision in the minimum time delay $t_{\rm min}$. McCarthy et al. (2020) predict a similar constraint using an observed galaxy catalog, finding that $\mathcal{O}(100)$ events will constrain $\alpha$ and $t_{\rm min}$ to 10% precision. Using a combination of population synthesis codes and simulated galaxy catalogs, Artale et al. (2019, 2020) connect the SFR of a given host galaxy to a merger’s time delay, finding that the DTD can vary greatly depending on the classification of the host galaxy. Other works such as those from Lamberts et al. (2016), Cao et al. (2018), Safarzadeh et al. (2019b, a); Safarzadeh & Berger (2019), and Fishbach & Kalogera (2021) have explored various aspects of hosts of binary mergers, including constraints on the DTD and hosts of BBHs. Evaluating these predictions with observations of mergers and their hosts will be a major component of XG science.

Finding CBC hosts with known redshifts is also useful for cosmology, particularly in measuring the Hubble constant $H_{0}$. Most cosmological studies with BBHs employ either a “spectral siren” method by fitting the black hole mass function to the BBH population (Farr et al., 2019; You et al., 2021; Ezquiaga & Holz, 2022) or a “dark siren” method by using galaxy catalogs and the GW localization (Oguri, 2016; Fishbach et al., 2019; Soares-Santos et al., 2019; Abbott et al., 2021; Alfradique et al., 2024). By uniquely identifying host galaxies using XG detectors, BBHs can be used as “golden dark sirens” to directly measure $H_{0}$ (Nishizawa, 2017; Gupta et al., 2023).

Previous studies have explored the localization capabilities of XG detectors and shown their promise. Vitale & Whittle (2018) showed with full Bayesian parameter estimation that certain events, such as the “gold-plated” GW150914, can be localized to within 1% uncertainty in distance and to a 90% sky area smaller than $10^{-2}$ deg² for a network comprised of CE and ET. Using Fisher matrix analyses, Borhanian & Sathyaprakash (2022), Pieroni et al. (2022), and Gupta et al. (2023) show that various XG detector networks will be able to localize hundreds to thousands of BBHs per year to within 1% in distance and 0.1 deg² on the sky. In particular, Borhanian & Sathyaprakash (2022) and Gupta (2023) point out that XG localizations of BBHs and NSBHs respectively could be powerful tools for resolving the Hubble tension (see Di Valentino et al. (2021) for a review). Chen et al. (2022) investigate the feasibility of identifying hosts of lensed BBHs using XG GW detectors and next-generation galaxy surveys, finding that enough lensed BBHs and hosts may be detected to begin to distinguish BBH formation channels.

These highly precise GW localizations constrain the possible site of the merger to a small volume in space, where for some events, only a single possible host galaxy will exist. This volume can then be cross-matched with existing galaxy catalogs to identify and characterize the host galaxy. As new EM facilities such as the Vera Rubin Observatory come online, galaxy catalogs will continue to improve (Ivezić et al., 2019). Where catalogs are not complete, targeted EM observations could search for and characterize the host, and reveal its morphology, stellar mass $M_{*}$, SFR, color, velocity dispersion, and local density around the host (Adhikari et al., 2020).

In this paper, we use a Bayesian analysis to determine the 3-dimensional localization volumes of NSBHs and BBHs in the XG era. From these localizations, we estimate the rate of mergers that can be directly associated with their hosts, using population estimates put forth by Abbott et al. (2023) and Biscoveanu et al. (2022). We find that a BBH merger will be localized to a volume smaller than 100 Mpc³, the average volume for one potential host galaxy, every eight days, and an NSBH merger every 1.5 months. In Sec. 2 we describe the simulation and localization procedure. Then, we present the results and show how they can be used to constrain population synthesis models in Sec. 3. We discuss the implications and limitations of our simulation and conclude in Sec. 4.

In our simulations, we assume one CE placed at the current LIGO Livingston site and one ET at the current Virgo site, using the CE and ET-D sensitivity curves from Abbott et al. (2017c). To simulate the BBHs, we use the population models released by Abbott et al. (2023). We choose the maximum likelihood distribution from the binary masses and redshift for the Powerlaw+Peak model; the mass model and power-law redshift model are described in Abbott et al. (2023) in Appendix B.1.b and Eq. 8 respectively. We list the chosen parameters for these distributions in Table 1. We simulate 1500 BBHs events by sampling the primary mass $m_{1}$, mass ratio $q$, and redshift $z$ from this distribution.

For the NSBH simulation, we use the population model defined in Biscoveanu et al. (2022) to simulate 1500 events. Biscoveanu et al. use a truncated power-law (parameterized by power-law index $\alpha$, minimum BH mass $m_{\rm BH,min}$, and maximum BH mass $m_{\rm BH,max}$, as described in Eq. 3 of their paper) to describe the BH mass distribution, and a truncated Gaussian model (described by mean $\mu$, standard deviation $\sigma$, the maximum neutron star mass $m_{\rm{NS,max}}$ and the drawn black hole mass $m_{\rm BH}$, in Eq. 4 of their paper) for the pairing function between the black hole and the neutron star. The redshift distribution is assumed to be uniform in comoving volume and source frame time. This is an appropriate assumption since the NSBH rate evolution is not expected to be significant within the relevant redshift range. Similar to the BBH simulation, we use the population hyperparameters from the maximum likelihood distribution to sample $m_{1}$, $q$, and $z$. The rate of NSBH mergers is taken from the median of the posterior distribution. These parameters are described in Table 2.

For both types of events, we use the IMRPhenomD waveform (Husa et al., 2016; Khan et al., 2016) to simulate the GW signals. The arrival time of the events are sampled uniformly in one year, the right ascensions and declinations are sampled uniformly across the sky, the polarization angle $\psi$ is uniform between 0 and $2\pi$, and the inclination angle $\iota$ is uniform in $\cos(\iota)$. We neglect the effect of spins to reduce computational costs. The use of information from spins has been shown to improve localizations (van der Sluys et al., 2008; Pankow et al., 2018), meaning that a future study involving spins could improve our results. Since events at large distances will have S/Ns too low to be localized to the volume of a single galaxy, we limit our simulation to $z=0.5$ for BBH and $z=0.3$ for NSBH. (This choice is validated by our results, which show that almost all events with single-host localizations are located within 1 Gpc, far below $z=0.3$.) We assume Planck18 cosmology (Aghanim et al., 2020) in the simulations.

These 1500 BBH events within $z=0.5$ correspond to 1.25 years of observations using the merger rate density at redshift zero $\mathcal{R}_{0}$ from the maximum likelihood, listed in Table 1. For the NSBHs, the 1500 events within $z=0.3$ will be observed in 6.26 years of XG detector operation.

We then produce the three-dimensional localization uncertainties of these events using the localization algorithm first presented in Chen & Holz (2017) and Chen & Holz (2016). The algorithm is a maximum likelihood method which takes as inputs the detector-frame masses, detector arrival time differences, phase differences, and S/Ns at each detector, then computes a three-dimensional sky localization. While it is not a full GW parameter estimation tool (e.g., including the mass and spin estimate) like LALInference (Veitch et al., 2015) or Bilby (Ashton et al., 2019; Romero-Shaw et al., 2020), the localization algorithm only takes $\mathcal{O}(1-10)$ minutes per event, even for XG S/Ns up to $\sim 1000$, and gives localization precisions consistent with those from full parameter estimation, as shown in Chen & Holz (2017). In addition, this method results in realistic three-dimensional localizations, where the localization regions are comprised of “voxels” which have sizes that increase with distance. This is in contrast to Fisher Information studies which assume Gaussian distributions in distance.

Understanding the rate of XG host galaxy identification requires an estimate of the average volume in which only a single potential BBH or NSBH host resides. Following Gehrels et al. (2016) and Chen & Holz (2016), we integrate the Schechter luminosity function (Schechter, 1976) describing the density of galaxies of a given luminosity $L$, written as $\rho(x)dx=\phi^{*}x^{a}e^{-x}dx$, where $x=L/L^{*}$ and $L^{*}$ is a characteristic luminosity describing the change between the exponential and power-law behaviors of $\rho$. We use the same $\phi^{*}$, $a$, and $L^{*}$ values from $B$-band observations of nearby galaxies as in Gehrels et al. (2016). Integrating the luminosity function, we find that 86% of the total galaxy luminosity is contained in galaxies with $L>0.12L^{*}$. Then, integrating the Schechter function to find the average number density of $L>0.12L^{*}$ galaxies, we arrive at $\phi^{*}\Gamma(a+1,0.12)$ = 0.01 Mpc^-3 (where $\Gamma$ is the incomplete gamma function, see Gehrels et al. (2016) for further details). Since luminosity is robustly correlated with stellar mass, we make the assumption that almost all BBH and NSBH mergers occur in galaxies with $L>0.12L^{*}$. Thus, we assume 100 Mpc³ is the average volume that contains only one potential host galaxy.

Of the 1500 simulated BBHs in this study, 56 had 90% probability localization comoving volumes smaller than 100 Mpc³. 26 BBH events were localized to 90% volumes within 10 Mpc³, with the best localized BBH event having a 90% volume of 0.015 Mpc³. This event was at a distance of 353 Mpc and was recovered with a network S/N of 923. The farthest BBH merger localized to within 100 Mpc³ had a luminosity distance of 2464 Mpc, corresponding to a redshift of 0.43, showing the utility of this method beyond just the local universe. 48 of the 1500 simulated NSBHs were localized to within the single-host volume of 100 Mpc³. 20 of these events were localized to within 10 Mpc³. The closest NSBH merger, at a luminosity distance of 63.5 Mpc, was localized to a volume of $3.5\times 10^{-3}$ Mpc³; the farthest NSBH localized to within 100 Mpc³ was at a luminosity distance of 634 Mpc.

For both BBHs and NSBHs, the distribution of 90% volumes is shown in Fig. 1 and the distribution of 90% areas is shown in Fig. 2. These figures also show the trend of increasing area and volume with luminosity distance, as expected.

As mentioned above, these 1500 BBH events (within $z=0.5$) will be detected in 1.25 years of XG observations. With 56 out of 1500 (3.7%) of such events being localized to within a single host galaxy, an average of 45 single-host BBH events will be detected every year; this is once every eight days. 48 of the 1500 (3.2%) of the NSBHs within $z=0.3$ were localized to within 100 Mpc^-3 for a rate of 7.7 per year, or once every month and a half.

Studies of host galaxies and their relationship with compact binary mergers by Adhikari et al. (2020), McCarthy et al. (2020), and Safarzadeh & Berger (2019), as mentioned in Sec. 1, give estimates for the number of mergers with hosts required to constrain the DTD. Though these studies focus on BNS due to the comparative ease of finding hosts using EM counterparts, their methods and results are extendable to EM-dark BBH and NSBH mergers.

Adhikari et al. (2020) and McCarthy et al. (2020) both predict that with only 10 events, which we find can be observed in about two and a half months for BBHs and just over a year for NSBHs, will be enough to make distinctions between formation channels which trace various galaxy properties. In particular, Adhikari et al. (2020) find that 10 events will be enough to distinguish between hosts that are weighted by SFR or stellar mass from a random sample of hosts. With $\sim 30$ events (eight months for BBHs and 3.9 years for NSBHs), we will be able to discern whether hosts tend to be stellar- or halo-mass weighted. Studies have found that halo mass correlates with globular cluster abundance (Harris et al., 2013; Hudson et al., 2014). A halo mass-weighted sample of hosts might then signify that globular clusters are a more likely site of BBH or NSBH formation compared to the galactic field (Adhikari et al., 2020).

Safarzadeh & Berger (2019) predict that with a fixed $t_{\rm min}$, the DTD power law index $\alpha$ will be constrained to within 30% with $\mathcal{O}(300)$ hosts—about seven years of observations of BBHs and 39 years of NSBHs. McCarthy et al. (2020) and Adhikari et al. (2020) are more optimistic, finding that $\mathcal{O}(100)$ events will be enough to constrain $\alpha$ to about ten percent precision. $t_{\rm min}$, despite being degenerate with $\alpha$, will also be constrained, though to a lesser degree. Thus with two years of observations of BBHs and 13 years of observations of NSBHs, their respective DTD slopes may be determined to 10% or tighter, and the minimum delay time to within a Gyr.

Our results are in broad agreement with Fisher matrix calculations from Borhanian & Sathyaprakash (2022), Gupta et al. (2023), and Pieroni et al. (2022). While the above studies do not directly compute 3D volume localizations, agreement in our area localizations lends support to the robustness of the 3D localization results presented in this paper.

Beyond identifying the host galaxy with GWs alone, it is also possible to determine the offset between the merger site and the center of its host galaxy with precise localizations, as can be done for events with EM counterparts (Abbott et al., 2017d), and which have been useful in studies of short gamma-ray bursts (Bloom et al., 2002; Fong & Berger, 2013; Zevin et al., 2020). The best localized BBH event in our sample has a 90% localization area of $5.1\times 10^{-5}$ deg² at a distance of 353 Mpc. Approximating this as a square and taking the root, this area corresponds to an angular resolution of 26 arcseconds. At 353 Mpc, this can resolve an offset of 43 kpc from the center of a potential host galaxy.

In addition to formation channels, constraining cosmology is another incentive for identifying hosts of BBHs and NSBHs. With our results showing that the identification of a unique host galaxy using GWs alone will be commonplace with XG detectors, we can use these BBHs and NSBHs as standard sirens (Schutz, 1986; Holz & Hughes, 2005; Abbott et al., 2017; Chen et al., 2018), since a precise redshift can be directly measured from the uniquely identified host galaxy. Cosmological measurements made from typically EM-dark BBH and NSBH mergers can then be used as a consistency check against those computed from BNS standard sirens, and help to mitigate any possible systematic errors in a given method. For example, methods relying on the identification of an EM counterpart may suffer from biases due to the inclination angle (Chen, 2020; Chen et al., 2023), which will have no impact on events with GW-only localizations.

One potential extension of our study involves the use of GW waveform models which include higher-order modes, especially for BBHs and NSBHs with highly asymmetric masses. While these models are more expensive to compute, they have been shown to tighten parameter estimation constraints on localization (Ajith & Bose, 2009; Kalaghatgi et al., 2020; Borhanian et al., 2020).

In simulating our BBH and NSBH populations, we necessarily assume a merger rate. While the LVK have detected enough BBH to start placing meaningful constraints on their population distribution, the rate of NSBH coalescences remains highly uncertain. We expect upcoming LVK observing runs to significantly reduce this uncertainty. Any future rate estimate will act as a scaling factor on the overall rate of detections; the 3% proportion of NSBH which will be localized to a single host galaxy for events within $z=0.3$ remains valid as long as the inferred NSBH population is similar to the one simulated here.

Our study neglects data analysis issues that have yet to be solved for XG detectors. These include the identification and isolation of individual events from a background of overlapping signals (Janquart et al., 2022; Samajdar et al., 2021), GW waveform accuracy and systematics at very high S/N (Pürrer & Haster, 2020; Huang et al., 2022; Hu & Veitch, 2022), and others mentioned in Couvares et al. (2021). Our analyses focus primarily on the localization of a particular set of GW events which tend to be at low redshift and recovered with high S/N. For these events, some as-yet unsolved XG analysis problems will not have significant effects; for example, they will stand out from any overlapping sources or noise. Additionally, we expect that many of these challenges will be solved by the time that XG detectors begin observations in the mid-2030s.

We note that our estimate of the host galaxy number density stems from integrating the Schechter function down to 0.12 $L^{*}$. If less luminous galaxies such as dwarf galaxies have nontrivial contributions to the compact binary population (e.g., Conselice et al. (2020)), then the average volume in which only a single host galaxy exists will be smaller than fiducial value of 100 Mpc³. In addition, this value is a statement of the average number density; in reality, galaxies tend to be clustered.

For GW localizations where galaxy catalogs are sparse, a minimal amount of telescope time will be required to characterize a potential host. As an example, a robust spectrum of a Milky Way-like galaxy at 2 Gpc is obtainable with an 1800 s exposure on an 8 m class telescope.

Beyond observations of CBC hosts, \Acdtd constraints can also be obtained in other ways. Using synthetic XG detections of BBH mergers, Vitale et al. (2019) simultaneously fit the SFR and the DTD. They indicated a similar or longer observation time with XG detectors compared to our two years (for BBH observations) to reach a similar 10% constraint on a characteristic DTD delay time. Since the two methods rely on different assumptions, it will be useful to compare the two approaches once XG data begins to flow to verify their results.

We note that precise localizations heavily depend on the size of the XG detector network. If only one XG detector operates, the prospects for constraining BBHs and NSBHs to within a single host galaxy with GWs alone become much more bleak. As an example, if the current best-localized BBH in the simulation (at 353 Mpc, localized to 0.015 Mpc³ with the ET + CE network) were localized with just ET, its 90% localization volume would balloon to 19 000 Mpc³. This emphasizes the importance of a global network of detectors in enabling future GW science.

We have shown that with XG detectors, singling out a unique galaxy as the host to a BBH merger with GWs alone will be a regular occurrence. Every eight days on average, an XG network comprised of CE and ET will detect and localize a BBH merger to a comoving volume smaller than 100 Mpc³, the average volume in which only a single probable host galaxy resides. This information will be invaluable in constraining the DTD and in determining whether hosts of these mergers trace SFR, stellar mass, galaxy morphology, or other parameters. Within eight months of observations, we will be able to determine if BBH host galaxies are stellar- or halo-mass weighted for BBHs, and thus whether mergers are more likely to be formed in globular clusters versus in the fields of galaxies. In two years of observations, the BBH DTD parameters may be constrained to within 30%. For NSBHs, for which the detection rates are much lower, a single-host localization will happen every six weeks. Finally, this method will allow for EM-determined redshifts from hosts of EM-dark BBHs and NSBHs, proving useful for standard siren cosmology as well.