Inferring Host-galaxy Properties of LIGO–Virgo–KAGRA’s Black Holes

The host galaxies of astrophysical electromagnetic (EM) transients have offered interesting insights into the astrophysics of their progenitors. For instance, the host-galaxy population of Type Ia supernovae and core-collapse supernovae are different (e.g., Foley & Mandel 2013); this is expected given that both these classes have different progenitors. Similarly, long and short gamma-ray bursts (GRBs) are known to trace different host-galaxy properties (Fong et al. 2010, 2022; Levesque 2014; Lyman et al. 2017; Nugent et al. 2022). The inferred star formation histories (SFHs) of short GRB host galaxies have also been used to constrain their delay-time distribution, finding a preference for low delay times but also evidence for a high delay-time tail (Zevin et al. 2022). Interesting conclusions about fast radio burst (FRB) progenitors have also been drawn using their host galaxies (Bhandari et al. 2022). The follow-up and cataloging of host galaxies of supernovae, short GRBs, and FRBs continues to be an important research direction.

In a similar vein, information about the host galaxies of gravitational-wave (GW) transients could shed light on many interesting aspects of their formation and evolution. For instance, Safarzadeh & Berger (2019) propose to use the population of host galaxies of binary neutron stars (BNSs) to measure their underlying delay-time distribution, whereas Adhikari et al. (2020) find that detections of localized BNSs would be sufficient to provide interesting constraints on BNS progenitor models, including the relation to star formation and delay-time distributions. Any information about probable host galaxies could also be used to reduce the number of candidate galaxies searched to localize a GW transient, and in designing optimal weighting schemes in the statistical method to measure cosmological parameters (Schutz 1986; Del Pozzo 2012; Chen et al. 2018; Fishbach et al. 2019; Abbott et al. 2021a, 2023a).

A number of works have made assumptions on the formation and evolution of both galaxies and compact binaries to forward model the distribution of GW transient host galaxies (O’Shaughnessy et al. 2010; Lamberts et al. 2016; Artale et al. 2019; Mapelli et al. 2019; Toffano et al. 2019; Santoliquido et al. 2022; Rauf et al. 2023; Srinivasan et al. 2023). For example, Santoliquido et al. (2022) found that, contingent on their modeling assumptions, binary black holes (BBHs) prefer merging in high-mass galaxies. However, unlike their EM cousins, the host galaxies of GW transients cannot usually be identified, owing to poor pointing accuracy of GW detectors. The notable exception to this is the BNS merger GW170817 (Abbott et al. 2017a), whose EM emission helped localize the source to the galaxy NGC 4993 (Abbott et al. 2017b). It comes as no surprise that although the LIGO–Virgo–KAGRA (LVK; Acernese et al. 2015; LIGO Scientific Collaboration et al. 2015; Akutsu et al. 2021) network of detectors has reported the detection of ∼90 BBHs (The LIGO Scientific Collaboration et al. 2021), none of them have been confidently associated with a host galaxy. ⁹ Even with upgrades to the current detector network with the inclusion of LIGO-India (Saleem et al. 2022), or with more sensitive next-generation detectors such as Cosmic Explorer (Evans et al. 2021) and the Einstein Telescope (Punturo et al. 2010), only a small fraction of BBHs are expected to be sufficiently localized such that the host galaxy can be confidently identified (Chen & Holz 2016; Nishizawa 2017; Abbott et al. 2020a; Borhanian et al. 2020; Branchesi et al. 2023; Gupta et al. 2023).

Thinking about the detected BBHs as a population instead of individual events has allowed measurements of the underlying mass and spin distribution of black holes in binaries, as well the merger rate of binaries as a function of redshift, z (Abbott et al. 2019, 2021b, 2023b). Specifically, the redshift evolution of the merger rate (Fishbach et al. 2018, 2021; Abbott et al. 2021b, 2023b) has yielded constraints on the delay-time distribution and formation metallicity of BBHs (Fishbach & Kalogera 2021; Fishbach & van Son 2023; Turbang et al. 2023), and the properties of globular clusters (GCs; Fishbach & Fragione 2023). In this work, we show that the redshift evolution of the merger rate inferred from the latest GW transient catalog GWTC-3 (The LIGO Scientific Collaboration et al. 2021) can already begin to constrain the host-galaxy parameter space of detected BBHs. This results directly from the observation that galaxies weighted by different physical properties have different evolutions for their (weighted) number density (see, e.g., Faber et al. 2007). For instance, there are more star-forming galaxies per unit volume at z = 1 as compared to z ≈ 0, while there is a lesser number of massive galaxies at higher redshifts.

Our work is structured as follows. In Section 2 we assume that the measured GWTC-3 merger rate is described by a simple mixture model of galaxies weighted by their star formation rate (SFR) and stellar mass, and infer the mixing fraction between them. We derive an upper limit for the possible contribution of BBHs from galaxies weighted by stellar mass. In Section 3 we constrain the delay-time distribution of BBHs using a physically motivated model. We also seed the SFHs of galaxies taken from UniverseMachine (Behroozi et al. 2019) with the inferred delay-time distribution to study the resulting host-galaxy properties, and present results derived using this procedure. We summarize our results and dwell on directions for future work in Section 4. Throughout we assume cosmological parameters consistent with the Planck 2018 data release (Planck Collaboration et al. 2020) as implemented in the Planck18 class in astropy ¹⁰ (Astropy Collaboration et al. 2013, 2018), and that galaxy formation histories are well emulated by UniverseMachine.

We first demonstrate that it is possible to differentiate between the host-galaxy populations of BBHs using the redshift evolution of the merger rate. If BBHs were hosted by a sample of host galaxies purely weighted by their SFR, the rate of BBHs would increase as a function of redshift, roughly as (1 + z)^2.5, up to z ∼ 1.5, following the evolution of the cosmic star formation rate density (cSFRD; ¹¹ Madau & Dickinson 2014; Madau & Fragos 2017; Behroozi et al. 2019), ¹² which peaks at that redshift. ¹³ On the other hand, if BBHs were hosted by a sample of galaxies purely weighted by their stellar mass, M_*, the BBH number density would decrease as a function of redshift; this is because galaxies build up their mass over cosmic time, and the total mass in stars per unit volume in the local Universe is greater than that at higher redshift. At every redshift where UniverseMachine ¹⁴ has support, we weight every galaxy in direct proportion to its corresponding stellar mass, finding that the weighted number density of the galaxies would scale as (1 + z)^−0.64.

Following Adhikari et al. (2020), we assume for simplicity that BBHs merge in a sample of galaxies either weighted by their stellar mass or SFR, with their relative abundances specified by a mixing fraction α_SM. Such a weighting in stellar mass and SFR physically corresponds to long and short delay times, respectively (see Section 3.2 for a discussion). The effective model for the rate of BBH mergers, R(z), can be written as,

We plot the inferred R(z) along with the stellar-mass- and SFR-weighted expectations in Figure 1. The black traces show the redshift evolution of the merger rate inferred from GWTC-3. By comparing the redshift evolution predicted by Equation (1) for different values of α_SM with the inference from GWTC-3, we ascertain α_SM < 43% [90% confidence limit (CL); uniform prior on α_SM], i.e., a maximum of 43% of the total mergers follow a sample of galaxies weighted by stellar mass, with the rest weighted by SFR. The hypothesis that BBHs are hosted purely by a set of galaxies weighted by their stellar mass (i.e., α_SM = 100%) is ruled out by the redshift evolution of the merger rate at >99.9%. As a corollary, at least 57% host galaxies should belong to a galaxy sample weighted by SFR. As mentioned earlier, this model is rather simple, but it sheds light on the fact that the redshift distribution is a powerful discriminator for the set of host galaxies within which GW events merge. For instance, in the statistical dark-siren analysis to measure the Hubble parameter, individual candidate host galaxies are typically weighted by luminosities in passbands that trace stellar mass or their SFR (see, e.g., Fishbach et al. 2019); our results show that a stellar-mass weighting is inconsistent with the redshift evolution of the BBH merger rate, and applying such a weight during the analysis should be avoided.

Zoom In Zoom Out Reset image size

We now turn to exploring a more physically motivated model for the evolution of binaries.

3.1. Redshift Evolution of the Merger Rate, and Constraints on the Delay-time Distribution

For the purposes of this work we model the redshift evolution of the merger rate as depending on the formation rate, R_f , of the progenitor binaries, and the distribution of delay times, p(t_D ), between the formation of the binary progenitors and subsequent merger of the binaries. Assuming these, the merger rate, R, as a function of lookback time, t, can be written as,

Since redshift, z, and lookback time, t, are related for a given cosmology, R(t) can also be equivalently written as R(z). We assume that R_f (t) follows the cSFRD as inferred by UniverseMachine (Behroozi et al. 2019). ¹⁵ For the delay-time distribution model, we employ a power-law prescription parameterized by a power-law index α with a minimum allowed delay time of and maximum allowed delay of 13.5 Gyr. That is,

The latest results from the LVK collaboration (Abbott et al. 2023b) derive constraints on the redshift evolution of the merger rate assuming that it follows R(z) ∼ (1 + z)^κ , and determine the posterior probability distribution on κ: p(κ∣data). The LVK GWTC-3 analysis assumed a flat prior on κ, and so the likelihood satisfied . In order to compare the model expectation from Equation (2) to redshift evolution inference, we need to map the delay-time distribution parameters onto a prediction on κ. Let us denote this model prediction by . Due to the functional form chosen for R(z), can be approximated as,

Since current observations measure the merger rate only out to z = 1, this is a good approximation; however, with more sensitive detectors, the parameters of the delay-time distribution would have to be constrained directly by jointly sampling over them along with hyperparameters of the underlying population model related to mass, spin, etc. ¹⁶ Using the approximation in Equation (4), we can write a likelihood function for α and ,

where the right-hand side is given by the LVK GWTC-3 population analysis.

We construct the joint posterior on α and using the likelihood function above. For this, we assume flat priors on α between [−4, 1], and on between [0, 6] Gyr. We use the Markov Chain Monte Carlo sampler emcee (Foreman-Mackey et al. 2013) to sample from the posterior distribution. Constraints on α and from GWTC-3 are shown in Figure 2. As expected, the constraints on the parameters are stronger than those derived from GWTC-2 (Fishbach & Kalogera 2021) due to the 50% additional events included in GWTC-3. We constrain α < −1.55 and at 90% CL. On the whole, the posteriors for both parameters rail against the lower limit of the prior, illustrating a preference for short delay times in the population. This is consistent with expectations that BBHs follow the cSFRD, and also with other recent work (Fishbach & Kalogera 2021; Fishbach & van Son 2023; Karathanasis et al. 2023; Turbang et al. 2023).

Zoom In Zoom Out Reset image size

3.2. Modeling Galaxy Formation and Evolution, and Constraining Host-galaxy Distributions

To obtain observed galaxy properties across cosmic time, we use the UniverseMachine (Behroozi et al. 2019) set of semianalytical galaxy formation simulations. Instead of directly simulating the gas and star formation during galaxy formation as is done in hydrodynamical simulations such as IllustrisTNG (Pillepich et al. 2018; Nelson et al. 2019) or EAGLE (Crain et al. 2015), UniverseMachine populates galaxies within dark-matter halos in a pure cold dark-matter simulation with a semiempirical model. This model relates the SFR to the halo mass accretion rate, and constrains this with a Monte Carlo scheme using a wide variety of observations over a range of redshifts out to z = 10. More details about UniverseMachine can be found in Behroozi et al. (2019). As a result, one can track the evolution of galaxy properties (e.g., stellar mass, host halo mass, and SFH) across redshift. We use publicly available UniverseMachine mock catalogs ¹⁷ created using the Bolshoi–Planck dark-matter-only simulations (Klypin et al. 2011) with a box size of 250 h⁻¹ Mpc and 2048³ particles. These catalogs allow us to extract galaxy properties in 170 logarithmically spaced redshift bins between z = 0 and z = 13. While we use UniverseMachine for results in this work, the procedure can be repeated with any galaxy formation simulation that tracks galaxy properties across redshift.

For each galaxy at redshift z₀, we calculate the merger rate as,

where is merger rate of the ith galaxy at redshift z, while is the SFH of the ith galaxy numerically interpolated across the redshift bins. The equation above is the same as Equation (2), except that the quantities pertain to a specific galaxy rather than the total merger rate across all galaxies. When the delay times are short, the merger rate of BBHs will have a redshift evolution very close to that of the SFR density, since in our model the rate of production of binaries is directly proportional to the rate of production of stars. On the other hand when the delay times are long, the galaxies would have evolved significantly, and would hence no longer necessarily be star forming at the merger epoch. These quiescent galaxies would have a very different evolution for their weighted number density as a function of redshift, which would be more similar to a sample of host galaxies weighted by their stellar mass rather than their SFR.

We illustrate this phenomenon in Figure 3, where we plot a two-dimensional BBH-rate-weighted histogram of the stellar mass and SFR of the z = 0.21 host galaxies assuming p(t_D ∣t₀) = δ(t_D − t₀). In this M_*–SFR plane, the distribution of galaxies has two broadly distinct branches: an upper branch corresponding to star-forming galaxies and a lower branch corresponding to quiescent galaxies. Within these two branches M_* and SFR are, on an average, correlated positively. We choose two representative values, t₀ = 0 Gyr and 10 Gyr, corresponding to short and long delay times, respectively. When the delay times are short, the star-forming branch is highly populated, while the quiescent branch is relatively empty—this is because binaries merge at roughly the same redshift that they formed, and their host galaxies do not have much time to quench their star formation (hence stay star forming). Naturally, the redshift evolution of the weighted number density in such a scenario follows the SFR. This scenario changes when the delay times are large—the quiescent branch is now more populated since galaxies have had a long time to evolve and quench their star formation. The level to which the quiescent branch is populated increases with an increase in the delay time.

Zoom In Zoom Out Reset image size

3.3. Results

In order to extract host-galaxy properties from the observed redshift evolution of the merger rate, we convolve the constraints on the delay-time distribution model obtained in Section 3.1 with galaxy SFHs taken from UniverseMachine. For each sample drawn from the joint posterior on α and , we track every galaxy in the catalog and assign to each one a merger rate of binaries according to Equation (2) at every epoch. We obtain a set of mock candidate BBH host galaxies, with complete information about their current stellar mass and SFR, mass assembly and SFH, as well as merger-rate history. Since UniverseMachine models the connection between galaxies and their host halos, we also obtain information about their host halo mass, radius, and peculiar velocity. Corresponding to every posterior sample, one can plot a histogram of galaxy properties weighted appropriately by the merger rate of BBHs in that galaxy; this BBH-rate-weighted histogram will indicate the population of galaxies hosting BBH events.

3.3.1. Stellar Mass, Star Formation Rate, and Mean Stellar Age

Figure 4 shows the median distribution of BBH host galaxies inferred in the M_*−SFR plane using our prescription at redshifts z = 0.21, 0.4, and 0.81. One noticeable broad feature is that the star-forming branch is shifted to higher values of SFR at high redshifts. This is because, consistent with expectations, galaxies at higher redshift on an average form more stars. Notably, there is significant support for BBH host galaxies both in the star-forming branch as well as the quiescent branch, although there is more support in the star-forming branch. This is a direct result of the small but nonzero delay times that are allowed by the constraints in Section 3.1—these delay times are large enough to induce some support for galaxies in the quiescent branch, but not so large as to entirely populate this branch. The median stellar mass of the hosts at z = 0.21 is 1.9 × 10¹⁰ M_⊙, while the median SFR is 0.9 M_⊙ yr⁻¹. The median stellar mass at z = 0.4 and z = 0.81 is 2.1 × 10¹⁰ M_⊙ and 2.1 × 10¹⁰ M_⊙ respectively, while the median SFR is 1.6 M_⊙ yr⁻¹ and 4.5 M_⊙ yr⁻¹ respectively.

Zoom In Zoom Out Reset image size

Figure 5 shows one-dimensional BBH-rate-weighted histograms of M_*, SFR, and mean stellar age, t_age, inferred from our prescription at z = 0.21. Here, t_age is defined on a per galaxy level as,

We also plot, for comparison, histograms for each of the parameters weighted equally by M_* and by SFR. The equal weights case should be thought of as the set of all galaxies in the Universe at a given redshift, while the M_*-weighted and the SFR-weighted galaxies track massive and star-forming galaxies, respectively. The inferred histograms are plotted in green, and the width of the green region should be thought of as coming from the uncertainties in measuring the parameters of the delay-time distribution.

Zoom In Zoom Out Reset image size

We see that the distribution of M_* is broadly consistent with the SFR-weighted sample of galaxies, and is inconsistent with an M_*-weighted sample of galaxies following expectations from Figure 1. Most of the inferred histograms cluster around the SFR-weighted histograms, consistent with the railing of α and toward the lower bound on the prior in Figure 2. A similar trend can be seen in the distribution of SFRs as well, which is consistent with an SFR-weighted sample of galaxies, but also has support for galaxies with slightly lower SFRs. The inferred mean stellar age of host galaxies is 4.4 Gyr, consistent with that of star-forming galaxies.

The general picture at z = 0.21 also holds at the higher redshifts which we considered, albeit with larger errors, since the merger rate at higher redshifts is less well constrained. We discuss results for z = 0.4, 0.81 briefly in the Appendix.

3.3.2. Magnitudes

The light that we observe from galaxies encodes signatures of their star formation and mass assembly history. Given information about the initial mass function and SFH of galaxies, one can forward model galaxy spectra using stellar population synthesis models (see, e.g., Bruzual & Charlot 2003; Conroy et al. 2009). These spectra can then be bandpassed appropriately to derive luminosities and magnitudes. For our inferred host galaxies, we derive (absolute) magnitudes in the Johnson–Morgan B band (Johnson & Morgan 1953) and the Two Micron All-Sky Survey Ks band (Skrutskie et al. 2006), ¹⁸ thought to be proxies for the SFR and stellar mass, respectively (Bell et al. 2003; Singer et al. 2016). To simulate mock galaxy spectra, we use the stellar population synthesis model from Bruzual & Charlot (2003) implemented in the bagpipes software package (Carnall et al. 2018). Since we already have simulated SFHs from UniverseMachine for our mock galaxies, we use the “custom” SFH implemented in bagpipes instead of relying on semianalytical prescriptions for the SFH. Also, since UniverseMachine does not model metallicity evolution, we fix the metallicity Z = 0.5 Z_⊙ at all redshifts, where Z_⊙ = 0.02 is solar metallicity. ¹⁹

As before, we plot the BBH-rate-weighted histogram of the inferred magnitudes in the B band and Ks band, along with histograms weighted equally by M_* and by SFR, in Figure 6. In both bands, the equally weighted set of galaxies has a preference for high magnitudes (i.e., low luminosities), since there are significantly more low-luminosity galaxies in the Universe (see, e.g., Schechter 1976). In contrast, the histogram for the inferred weights peaks at lower magnitudes (i.e., higher luminosities). The median inferred magnitude in the B band is −19.6, while that in the Ks band is −21.5. These inferred histograms of the magnitudes can be directly used to calculate luminosity weights in the dark-siren method for measuring the cosmic expansion history with GW events. In addition, they may enable the design of more optimized follow-up strategies to search for possible EM counterparts to BBH mergers.

Zoom In Zoom Out Reset image size

3.3.3. Halo Mass, Halo Radius, and Peculiar Velocity

The formation and evolution of galaxies is intrinsically connected to their host dark-matter halos (see Wechsler & Tinker 2018 for a review); thus, the properties of a galaxy are correlated with its parent dark-matter halo. ²⁰ Specifically, larger halos host larger and more massive galaxies. This correlation also naturally causes massive galaxies to have a greater peculiar velocity, due to the denser environment (see, e.g., Sheth & Diaferio 2001).

Since UniverseMachine also tracks host halo properties along with the galaxy properties, we also have host halo information for the host galaxies. Figure 7 shows the rate-weighted histograms of inferred host dark-matter halo mass, M_halo, halo radius, r_halo, and galaxy peculiar velocity, v_pec. We note that the M_*-weighted galaxies have higher M_halo, v_pec, and r_halo, as compared to those weighted by SFR; this is because the former galaxies are larger and hence lie in larger dark-matter halos. Similar to the case of galaxies presented in the previous subsections, we find that the inferred histograms for host dark-matter halos lie closest to the halos associated with SFR-weighted galaxies. The median M_halo and r_halo are 5.4 × 10¹¹ h⁻¹ M_⊙ and 180 pc, respectively. The median v_pec is 150 km s⁻¹, but has a tail extending out to v_pec ≈ 800 km s⁻¹. This means that the peculiar velocities can shift the redshifts of BBHs at most by v_pec/c ≈ 0.003, and thus are unlikely to affect the uncertain estimates of distances derived from BBHs even in next-generation detectors.

Zoom In Zoom Out Reset image size

3.3.4. Large-scale Bias

Galaxies are biased tracers of the underlying matter density in the Universe—they prefer forming at rare, highly dense peaks of the underlying matter density field. Kaiser (1984) showed that the clustering of galaxies, quantified by the two-point correlation ξ_gal(r), can be modeled as a constant enhancement over the clustering of the underlying matter distribution on large scales: . The value of the bias, b_gal, is directly related to the masses of halos that host galaxies (see, e.g., Cooray & Sheth 2002).

Since BBHs also form and evolve in galaxies, they are a biased tracer of the underlying matter density. We denote b_BBH as the large-scale bias of BBHs, akin to the same quantity for galaxies. Indeed, multiple works have proposed using the bias to probe the astrophysical formation channels of GW transients (Namikawa et al. 2016; Raccanelli et al. 2016; Scelfo et al. 2018; Mukherjee et al. 2021; Vijaykumar et al. 2023). Using the halo mass–bias relation of Tinker et al. (2010), implemented in the colossus software package (Diemer 2018), we are able to predict values for the large-scale bias of BBHs from the inferred host-galaxy distributions.

We plot the 90% credible region of the inferred BBH-rate-weighted bias as a function of redshift in Figure 8. Overall, the bias has an increasing trend toward high redshift—this is a natural consequence of hierarchical structure formation in the Universe. We note that the bias at z = 0 is predicted to be in the range 0.86–0.88 at z = 0.1 at 90% CL, while that at z = 0.9 is predicted to be 1.16–1.40 (90% CL). The inferred bias is consistent with that of the SFR-weighted galaxies, but is inconsistent with the M_*-weighted galaxies; this is because the bias is purely a function of M_halo, and the inferred M_halo distribution is inconsistent with the M_*-weighted galaxies (see Figure 7).

Zoom In Zoom Out Reset image size

We have shown that the measured redshift evolution of the BBH merger rate from GWTC-3 constrains the set of host galaxies with which BBHs are associated. We started with a simple mixture model of galaxies weighted either by SFR or stellar mass. We demonstrated that the hosts of merging BBHs do not exclusively track stellar mass at the >99.9% level—at most 43% of the host galaxies come from a stellar-mass-weighted sample. We then used UniverseMachine simulations, along with data-driven constraints on the BBH delay-time distribution, to infer the host-galaxy population. Our results indicate that most BBHs are hosted by actively star-forming galaxies, but there is support for a small fraction of them being hosted by quiescent galaxies. We provide inferred distributions for the halo masses, radii, peculiar velocities, and stellar ages. We also make predictions for the observed Ks- and B-band magnitudes of BBH host galaxies, as well as the large-scale bias associated with the clustering of BBHs. Our results assume that R(z) ∼ (1 + z)^κ is a good description for the redshift evolution of the merger rate. There are some hints of structure in R(z) beyond this power-law prescription (Edelman et al. 2023; Callister & Farr 2024); if these results hold up with newer data, it might necessitate invoking the presence of two (or more) populations of BBHs with distinct delay-time distributions. However, the inferred delay-time distribution would still prefer short delay times, and hence we do not expect our results to change significantly.

Unlike previous works which provide distributions of host galaxies based on prescriptions from population synthesis, our framework uses constraints from observed GW data to solve the inverse problem: we infer a set of host galaxies without assumptions for any specific prescription for binary formation and evolution. While we have concentrated on BBHs in this work, our framework will extend trivially to BNSs and NS–BHs, once we are able to measure the redshift evolution of their merger rates. This framework can also be extended to infer the host-galaxy population of astronomical objects for which source localization is poor (e.g., poorly localized GRBs and FRBs). Our inferred distributions of various parameters can be used to design optimal weighting schemes for the dark-siren method to constrain the expansion history of the Universe. In addition, our predictions for host-galaxy properties may aid in the design of optimal follow-up strategies for BBH events. This will be of particular interest once the BNS and NS–BH extensions to our methods have been implemented, since those searches are expected to have transient EM counterparts.

For the main results in this work, we have assumed that binaries followed a physically motivated merger-rate evolution, given by the UniverseMachine cSFRD convolved with a delay-time distribution. This prescription is expected to work well if all binaries form and merge in isolated galactic field environments, and our results should be interpreted in this context. For instance, binaries could also form and merge in dynamical environments such as GCs or AGN, and these populations will have different progenitor formation rates and delay-time distributions, thus requiring changes to our modeling procedure. Our framework can be readily extended to include these channels given some assumption of the progenitor formation rate, although adding new parameters is likely to increase the statistical error on the results. Since compact binary formation is expected to be a strong function of metallicity (Chruślińska 2024, and references therein), another natural extension of our framework would be to include the metallicity-specific SFH from hydrodynamical simulations of galaxy and structure formation while assigning merger rates to individual galaxies. This would also necessitate using a metallicity-specific cSFRD as the progenitor formation rate in constraining delay-time distributions. It is also possible to relax our assumption of specific progenitor formation rates and delay-time distributions, and instead generalize our approach to compare redshift evolution of galaxies and BBHs in a nonparametric fashion, although this would still be subject to any modeling differences inherent in different galaxy formation models.

We have demonstrated the power of the observed properties of the existing GW population to constrain properties of the associated host galaxies. The relation between GW mergers and their parent galaxies offers a rich vein to mine for astrophysical insights. A more precise future measurement of the rate evolution of binaries will help us further constrain the host environments, possibly helping to disentangle the various formation channels that contribute to the total population of merging binaries observable with ground-based GW detectors.

We are grateful to Jacqueline Antwi-Danso for patiently answering all our questions about galaxy spectral energy distributions, and Peter Behroozi for help with interpreting the UniverseMachine catalogs. We also thank Surhud More for useful discussions, and Alexandra Hanselman and Aditya Kumar Sharma for a reading of our draft. All computations for this work were performed on the Alice cluster at ICTS-TIFR.

A.V.’s work is supported by the Department of Atomic Energy, Government of India, under Project No. RTI4001. A.V. also acknowledges support by a Fulbright Program grant under the Fulbright-Nehru Doctoral Research Fellowship, sponsored by the Bureau of Educational and Cultural Affairs of the United States Department of State and administered by the Institute of International Education and the United States-India Educational Foundation. A.V. and M.F. are supported by the Natural Sciences and Engineering Research Council of Canada (NSERC; funding reference number 568580). D.E.H is supported by NSF grants AST-2006645 and PHY-2110507, as well as by the Kavli Institute for Cosmological Physics (KICP) through an endowment from the Kavli Foundation and its founder Fred Kavli.

We thank KICP for support through “The quest for precision gravitational wave cosmology” workshop where this work was initiated, and ICTS for support through the Largest Cosmological Surveys and Big Data Science (code: ICTS/BigDataCosmo2023/5) workshop. This research has made use of the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00.

This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation.

Software: numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration et al. 2013, 2018), jupyter (Kluyver et al. 2016), pandas (McKinney 2010), seaborn (Waskom 2021), bagpipes (Carnall et al. 2018), and colossus (Diemer 2018).

In this section we provide inferred rate-weighted histograms of M_*, SFR, t_age, M_halo, r_halo, and v_pec, at redshifts of 0.40 (Figure 9) and 0.81 (Figure 10). We also show inferred histograms of B-band and Ks-band absolute magnitudes in Figure 11. We reiterate that our prescription is able to infer galaxy properties at any redshift where a UniverseMachine simulation snapshot exists.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size