Abstract
We use 47 gravitational wave sources from the Third LIGO–Virgo–Kamioka Gravitational Wave Detector Gravitational Wave Transient Catalog (GWTC–3) to estimate the Hubble parameter H(z), including its current value, the Hubble constant H0. Each gravitational wave (GW) signal provides the luminosity distance to the source, and we estimate the corresponding redshift using two methods: the redshifted masses and a galaxy catalog. Using the binary black hole (BBH) redshifted masses, we simultaneously infer the source mass distribution and H(z). The source mass distribution displays a peak around 34 M⊙, followed by a drop-off. Assuming this mass scale does not evolve with the redshift results in a H(z) measurement, yielding 
(68% credible interval) when combined with the H0 measurement from GW170817 and its electromagnetic counterpart. This represents an improvement of 17% with respect to the H0 estimate from GWTC–1. The second method associates each GW event with its probable host galaxy in the catalog GLADE+, statistically marginalizing over the redshifts of each event's potential hosts. Assuming a fixed BBH population, we estimate a value of 
with the galaxy catalog method, an improvement of 42% with respect to our GWTC–1 result and 20% with respect to recent H0 studies using GWTC–2 events. However, we show that this result is strongly impacted by assumptions about the BBH source mass distribution; the only event which is not strongly impacted by such assumptions (and is thus informative about H0) is the well-localized event GW190814.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
The discovery of a gravitational wave (GW) signal from a binary neutron star (BNS) merger (Abbott et al. 2017a) and the kilonova emission from its remnant (Coulter et al. 2017; Abbott et al. 2017b) provided the first GW standard siren measurement of the cosmic expansion history (Abbott et al. 2017c). As pointed out by Schutz (1986), the GW signal from a compact binary coalescence directly measures the luminosity distance to the source without any additional distance calibrator, earning these sources the name "standard sirens" (Holz & Hughes 2005). Measuring the cosmic expansion as a function of the cosmological redshift is one of the key avenues with which to explore the constituents of the universe, along with the other canonical probes such as the cosmic microwave background (CMB; Spergel et al. 2003, 2007; Komatsu et al. 2011; Ade et al. 2014, 2016; Aghanim et al. 2020), baryon acoustic oscillations (Eisenstein & Hu 1998, 1997; Dawson et al. 2013; Alam et al. 2017), type Ia supernovae (Riess et al. 1996; Perlmutter et al. 1999; Riess et al. 2016; Freedman 2017; Riess et al. 2019), strong gravitational lensing (Suyu et al. 2010; Wong et al. 2020), and cosmic chronometers (Moresco et al. 2016; Jimenez et al. 2019).
Even though GW sources are excellent distance tracers, using them to study the expansion history also requires measurement of their redshift. The redshift information is usually degenerate with the source masses in the GW signal, as the redshifted masses affect the GW frequency evolution. However, several techniques are proposed to infer the redshift of GW sources and break the mass–redshift degeneracy. For sources with confirmed electromagnetic (EM) counterparts, the host galaxy and its redshift can be determined directly (Holz & Hughes 2005; Dalal et al. 2006; MacLeod & Hogan 2008; Nissanke et al. 2010; Abbott et al. 2017c; Chen et al. 2018; Feeney et al. 2019). For sources without an EM counterpart, alternative techniques to infer the source redshift include comparing the redshifted mass distribution to an astrophysically motivated source mass distribution (Chernoff & Finn 1993; Taylor & Gair 2012; Farr et al. 2019; You et al. 2021; Mastrogiovanni et al. 2021), obtaining statistical redshift information from galaxy catalogs (Schutz 1986; MacLeod & Hogan 2008; Del Pozzo 2012; Nishizawa 2017; Chen et al. 2018; Nair et al. 2018; Fishbach et al. 2019; Soares-Santos et al. 2019; Gray et al. 2020; Yu et al. 2020; Palmese et al. 2020; Borhanian et al. 2020; Finke et al. 2021), comparing the spatial clustering between GW sources and galaxies (Oguri 2016; Bera et al. 2020; Mukherjee et al. 2020, 2021), leveraging external knowledge of the source redshift distribution (Ding et al. 2019; Ye & Fishbach 2021), and exploiting the tidal distortions of neutron stars (Messenger & Read 2012; Chatterjee et al. 2021).
The third LIGO–Virgo–Kamioka Gravitational Wave Detector (KAGRA) GW transient catalog (GWTC–3; Abbott et al. 2021e) contains 90 compact binary coalescence candidate events with at least a 50% probability of being astrophysical in origin. Out of the events from the third observing run, a notable EM counterpart has been claimed only for the high-mass binary black hole (BBH) event GW190521 (Graham et al. 2020). However, the significance of this association has been reassessed by several authors, who found insufficient evidence to claim a confident association (Ashton et al. 2021; De Paolis et al. 2020; Palmese et al. 2021). In this work, we do not include the redshift information from this putative EM counterpart signal. The only standard siren with an EM counterpart in GWTC–3 remains the BNS GW170817.
The remainder of this paper is organized as follows. In Section 2 we discuss the statistical methods used to infer the cosmological parameters with and without galaxy catalog information. In Section 3 we discuss the properties of the GW events and galaxy catalogs used in this paper. In Section 4 we present the results of our analyses, and in Section 5 we discuss their implications for the cosmological parameters. Finally, in Section 6 we present our conclusions.
2. Method
We use two analysis methods 301 : (i) jointly fitting the cosmological parameters and the source population properties of BBHs, without using galaxy catalog information (Mastrogiovanni et al. 2021; Section 2.1), and (ii) fixing the source population properties, and inferring the cosmological parameters using statistical galaxy catalog information (Gray et al. 2020; Section 2.2).
2.1. Hierarchical Inference without Galaxy Surveys
The GW event catalog can be described by two sets of parameters: a set of population hyperparameters Φ that are common to the entire population of GW sources, and a set of intrinsic parameters that are unique for each event. The cosmological population hyperparameters in this work are the cosmological parameters for a flat universe. For the redshift range considered in the analysis, the contribution to the total energy density from radiation and neutrinos is negligible. Hence, we consider the dark energy density today as ΩDE(z) = 1 − Ωm. The cosmological parameters considered are, therefore, the Hubble constant H0, the matter density Ωm, and the dark energy equation-of-state (EOS) parameter w(z) = w0. (Chevallier & Polarski 2001). The dark energy EOS is defined by w = p/ρ, where for the standard ΛCDM we have w = −1. Additionally, the set of hyperparameters Φ includes parameters describing the source mass distribution and the merger rate density as a function of the redshift.
Given a set of Nobs GW detections associated with the data 
, the posterior on Φ can be expressed as (Mandel et al. 2019; Thrane & Talbot 2019; Vitale et al. 2021)
where p(Φ) is a prior on the population parameters including standard cosmological parameters and parameters governing the redshift and mass distributions of GW sources. The variable θ is the set of parameters intrinsic to each GW event, such as spins, masses, redshift etc., p(xi
∣Φ, θ) is the GW likelihood, 
is the probability of detecting a GW event with intrinsic parameters θ and for population hyperparameters Φ, and ppop(θ∣Φ) is a population modeled prior. The denominator in Equation (1) correctly normalizes the posterior and takes into account selection effects (Mandel et al. 2019). We use the hierarchical statistical framework to infer the population parameters Φ and the prior that they induce on the distributions of the GW parameters.
The intrinsic GW parameters that are interesting for cosmology are those that provide information about the redshift z of the source. For sources that are detected at cosmological distances, GWs provide a measurement of the redshifted masses 
and the luminosity distance DL, rather than the redshift and source masses m1, m2, where
The relation between the source mass and redshifted mass can then be used to probe cosmology even in the absence of an explicit EM counterpart (Taylor et al. 2012; Taylor & Gair 2012) provided the source mass scale can be well characterized. This approach is more effective if the source mass distribution displays sharp features (Farr et al. 2019; Ezquiaga & Holz 2021; You et al. 2021; Mastrogiovanni et al. 2021).
2.1.1. Population models
We model the BBH population as BBHs represent the majority of the detected sources. We now give a general overview of the source mass and redshift models used in this paper; see Appendix A for a complete description of the population models.
We describe the underlying distribution in the redshift and source masses as
where C is a normalization factor, p(m1, m2∣Φm ) is the source frame mass distribution, Φm refers to all the population parameters not related to cosmology, the (1 + z) term encodes the clock difference between the source frame and detector frame, and p(z∣H0, w0, Ωm) is the redshift prior, which is taken to be uniform in comoving volume. The term ψ(z∣γ, k, zp) describes the redshift evolution of the merger rate with a parameterization similar to that of Madau & Dickinson (2014), characterized by a low-redshift power-law slope γ, a peak at redshift zp, and a high-redshift power-law slope k after the peak, or
The above rate evolution model is more complex than that in Fishbach et al. (2018), Abbott et al. (2021b), and Abbott et al. (2021f); this is because when we vary the cosmological parameters, the GW observations may be pushed to higher redshift, z > 2, past the peak of the star formation rate (Madau & Dickinson 2014).
The source mass models are factorized as
where the secondary mass is modeled by a power-law distribution between a minimum mass 
and maximum mass m1. For the primary mass, we implement three phenomenological mass models used in Abbott et al. (2019a, 2021b).
The first phenomenological model, the Truncated model, describes the mass distribution as a power law (PL) between a minimum mass 
and a maximum mass 
(Fishbach & Holz 2017). This model was originally proposed in order to characterize the mass spectra of BBHs using two features: a lower hard mass edge representing the gap between NSs and BHs and an upper mass edge representing the imprint of (pulsational) pair-instability supernovae (Fryer et al. 2001; Heger & Woosley 2002; Farmer et al. 2019; Renzo et al. 2020; Umeda et al. 2020). The BBH mass distribution inferred from GWTC–2 was more complicated than the Truncated, and the second and third models are extensions of the Truncated model that contain more complex structures to better fit the mass distribution (Abbott et al. 2021b). These models were proposed to describe more complex morphologies of the mass spectra such as a possible smooth transition between NSs and BHs and the possibility of having massive BHs formed in channels that are not subject to the (pulsational) pair-instability supernovae, e.g., hierarchical mergers (Mapelli 2021). The second model, the Broken Power Law, consists of two PLs attached at a break point (Abbott et al. 2021b). The third model is a superposition of a Truncated and a Gaussian component referred to as Power Law + Peak, in which the primary mass distribution is described by a PL with the addition of a Gaussian peak with mean μg and variance 
(Talbot & Thrane 2018). Using the Broken Power Law and Power Law + Peak models, GWTC–2 revealed an excess of BBH systems with primary masses in the range ∼30–40 M⊙, followed by a drop-off in the merger rate at high masses (Abbott et al. 2021b). This structure in the PL, modeled either as a break or a Gaussian peak, may represent the imprint of (pulsational) pair-instability supernovae.
In a companion paper investigating the GWTC–3 population (Abbott et al. 2021f) we show the evidence for substructures in the BHs primary mass spectrum around ∼ 10M⊙. However we also find that the simpler Power Law + Peak model is still one of the models preferred by the GWTC–3 data. For this reason, and for simplicity in this paper, we only adopt models that are characterized by a single structure (Broken Power Law and Power Law + Peak) to describe the excess of BHs observed around 35M⊙; this also corresponds to the binaries that we can observe at higher redshifts and for which source mass assumptions could be important.
In order to infer H(z) from the BBH population, a crucial assumption is that the source mass distribution is independent of the redshift (Fishbach et al. 2021). In most BBH formation scenarios, we expect some evolution of the mass distribution with the redshift, due to factors such as the metallicity evolution of the universe (Kudritzki & Puls 2000; Belczynski et al. 2010) and the dependence of the delay time between BBH formation and merger on BBH properties (Kushnir et al. 2016; Gallegos-Garcia et al. 2021; van Son et al. 2022). Nevertheless, if mass features, such as the break in the Broken Power Law model or the peak in the Power Law + Peak model, are caused by the pair-instability supernova process, their location is thought to stay constant within a few solar masses across cosmic time (Farmer et al. 2019). The presence of these sharp mass features drives our cosmological constraints (Farr et al. 2019; Mastrogiovanni et al. 2021). Moreover, the BBH mass distribution is expected to evolve only weakly over the range of redshift accessible to current observations, at a level below current statistical uncertainties (Fishbach & Kalogera 2021; van Son et al. 2022). Although the BH mass spectrum at formation may vary with cosmic time, BBH channels typically predict a wide distribution of delay times between formation and merger, which tends to wash out any dependence of the BBH mass on the merger redshift (Mapelli et al. 2019).
In the following we will neglect the selection effect of spin distribution as the detection probability due to their inclusion should not vary by more than a factor of 2 (Ng et al. 2018b). This is indeed a negligible term with respect to the statistical uncertainties on our posteriors (see Section 5) and the dependence of the selection bias on other parameters, such as H0 and γ, for which it nearly follows a power law.
2.2. Statistical Galaxy Catalog Method
We also use the gwcosmo code (Gray et al. 2020) in the pixelated sky scheme (Gray et al. 2022), i.e., using the HEALPix pixelization algorithm (Górski et al. 2005; Zonca et al. 2019), to infer H0 using information from galaxy surveys. This method assumes a fixed source mass distribution, as well as a fixed-rate evolution for the binaries, and estimates H0 from the GW data using galaxy catalogs to provide statistical information about the GW source redshifts.
When including galaxy catalog information, the prior on redshift can be replaced by the distribution of galaxies in the survey. However, Equation (1) needs to be modified in order to take into account completeness corrections. These extra terms account for the impact of incompleteness, i.e., missing galaxies, due to the limited sensitivity of the catalog, in the GW localization volume. In this case, the posterior is given by
where G is the hypothesis that the GW host galaxy is included in the catalog and 
that it is not, and 
expresses their probabilities for ![$g\in [G,\bar{G}]$](https://content.cld.iop.org/journals/0004-637X/949/2/76/revision1/apjac74bbieqn12.gif)
. The term p(Nobs∣H0, Φm
) is the probability of having Nobs detections. We analytically marginalize over this by assuming a uniform in the log rate prior. The notation 
indicates the hypothesis that an event has been detected. The likelihoods 
are built from the GW data and corrected for the selection effects in the case that the host galaxy is, and is not, inside the catalog; see Gray et al. (2020).
We implement an improved version of the analysis presented in Abbott et al. (2021a) that can estimate H0 for any given sky direction covered by the GW localization by dividing the sky into equal-area pixels. In each pixel, the apparent magnitude threshold (mthr) is taken to be the median of the apparent magnitudes of all the galaxies inside that pixel. This assumption is a conservative choice for approximating the impact of catalog completeness: all galaxies with apparent magnitude fainter than the defined threshold are excluded from the analysis. Using this mthr the completeness is assessed, and the H0 likelihoods are calculated in each pixel. In the end, all the pixel likelihoods are combined using weights proportional to the GW posterior probability in each pixel to give the final H0 posterior of each GW event. Pixels with no GW support have zero contribution, so only the pixels within the 99.9% sky area are used.
These improvements are necessary to correct for galaxy catalog incompleteness in the case that the galaxy surveys contained within the catalog are less sensitive in particular sky areas, such as the directions of the galactic plane. Moreover, the analysis can take into account the fact that the GW luminosity distance posterior conditioned on the sky position might significantly change between different sky positions. The combination of galaxy redshift information and luminosity distance estimation changes from pixel to pixel, leading to a more robust estimation of H0.
3. Event and Catalog Selection
3.1. GW Events
For our main result, we select 47 GW events with a network matched filter signal-to-noise ratio (S/N) > 11 and an inverse false-alarm rate (IFAR) higher than 4 yr, taking their maximum across the different search pipelines from the data distributions released with GWTC–1 (Abbott et al. 2019c), GWTC–2 (Abbott et al. 2021c), and GWTC–3 (Abbott et al. 2021e), with no plausible instrumental origin. For events released with GWTC–3, we use the first version of the data distribution associated with LIGO Scientific Collaboration et al. (2021). We have verified that the selection of GW events and results of our analysis do not significantly change with the second version of the data release associated with GWTC–3.
We also consider events identified by different S/N choices to explore possible systematics in the computation of selection effects (see Appendix B). In the remainder of the paper we will shortly refer to these different ensembles by quoting the threshold S/N choices. Of the 47 events with S/N > 11, 42 are BBH detections, 2 are the BNS events GW170817 (Abbott et al. 2017a) and GW190425 (Abbott et al. 2020a), 2 are the NSBH events GW200105 and GW200115 (Abbott et al. 2021d), and 1 is the asymmetric mass binary GW190814 (Abbott et al. 2020b). A visual representation of the population of BBHs that we detected is provided in Figure 1, where we show the distribution of detector frame masses and luminosity distance of the BBHs. We have tabulated all the GW sources used in this analysis in Table 1, mentioning their source properties, sky localization error, 3D localization volume, number of galaxies in the catalog within the localization volume, and the probability that the GW host is present in the GLADE+catalog (Dálya et al. 2018, 2022). Note that, unlike in Abbott et al. 2021e, the estimation of masses and distances are reported using a prior 
and not uniform in comoving volume as we are interested in showing these values using cosmology-agnostic priors. For the events detected during O1, O2, and O3, we use combined posterior samples from the IMRPhenom (Thompson et al. 2020; Pratten et al. 2021) and SEOBNR (Ossokine et al. 2020; Matas et al. 2020) families, while for the two NSBH events GW200105 and GW200115 we use posterior samples generated with low-spin priors (Abbott et al. 2021d).
302
Figure 1. Distribution of the mass and luminosity distance parameters for the 42 BBH events with S/N> 11. The figure is generated using a Planck cosmology with H0 = 67.9 km s−1 Mpc−1 and Ωm = 0.3065 with a D2 L prior and a uniform prior on the detector frame masses ad stacking posterior samples for each event. This figure is only representative of the events reported in Table 4 and does not indicate the population reconstruction. Left: distribution of the primary detector frame masses (blue solid line) and source frame masses (orange dashed line). Middle: same but for the secondary source mass. Right: distribution of the luminosity distance (bottom axis) and redshift (top axis).
Download figure:
Standard image High-resolution imageTable 1. GW Events Considered in This Paper and Properties Reported With Their Median Values and 90% Symmetric Credible Intervals
| Name | S/N | IFAR |
|
| DL | m1 | m2 | z | ΔΩ | ΔVc | Ngal | p(G∣z, H0) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| [yr] | [M⊙] | [M⊙] | [Mpc] | [M⊙] | [M⊙] | [deg2] | [Gpc3] | |||||
| GW150914 | 24.4 | 1 × 107 |
|
|
|
|
|
| 200 | 3 × 10−4 | 3000(3 × 104) | 50-5(99-72) |
| GW151012 | 10.0 | 100 |
|
|
|
|
|
| 2000 | 0.05 | 6000(5 × 105) | 5-0(73-10) |
| GW151226 | 13.1 | 1 × 107 |
|
|
|
|
|
| 1000 | 0.003 | 1 × 104(1 × 105) | 58-4(100-70) |
| GW170104 | 13.0 | 1 × 107 |
|
|
|
|
|
| 900 | 0.02 | 8000(4 × 105) | 8-0(76-17) |
| GW170608 | 15.4 | 1 × 107 |
|
|
|
|
|
| 400 | 4 × 10−4 | 4000(2 × 104) | 69-19(100-87) |
| GW170729 | 10.8 | 6 |
|
|
|
|
|
| 1000 | 0.3 | 60(8 × 104) | 0-0(12-0) |
| GW170809 | 12.4 | 1 × 107 |
|
|
|
|
|
| 300 | 0.004 | 2000(1 × 105) | 4-0(70-21) |
| GW170814 | 16.3 | 1 × 107 |
|
|
|
|
|
| 90 | 2 × 10−4 | 1000(1 × 104) | 27-1(91-60) |
| GW170817 | 33.0 | 1 × 107 |
|
|
|
|
|
| 20 | 2 × 10−8 | 10(6) | 100-100(100-100) |
| GW170818 | 11.3 | 2 × 104 |
|
|
|
|
|
| 40 | 6 × 10−4 | 80(1 × 104) | 2-0(55-9) |
| GW170823 | 11.5 | 1 × 107 |
|
|
|
|
|
| 2000 | 0.2 | 700(3 × 105) | 0-0(35-0) |
| GW190408_181802 | 14.8 | 1 × 105 |
|
|
|
|
|
| 100 | 0.003 | 80(2 × 104) | 0-0(39-3) |
| GW190412 | 19.7 | 1 × 105 |
|
|
|
|
|
| 20 | 3 × 10−5 | 100(4000) | 7-0(76-51) |
| GW190413_134308 | 10.3 | 6 |
|
|
|
|
|
| 500 | 0.2 | 0(70) | 0-0(0-0) |
| GW190421_213856 | 10.5 | 400 |
|
|
|
|
|
| 1000 | 0.2 | 1(2 × 104) | 0-0(5-0) |
| GW190425 | 12.9 | 30 |
|
|
|
|
|
| 8000 | 0.002 | 3 × 104(6 × 104) | 95-63(100-100) |
| GW190503_185404 | 12.8 | 1 × 105 |
|
|
|
|
|
| 90 | 0.005 | 100(4 × 104) | 0-0(46-1) |
| GW190512_180714 | 12.4 | 1 × 105 |
|
|
|
|
|
| 200 | 0.008 | 100(3 × 104) | 0-0(41-2) |
| GW190513_205428 | 12.9 | 8 × 104 |
|
|
|
|
|
| 500 | 0.04 | 20(3 × 104) | 0-0(18-0) |
| GW190517_055101 | 11.3 | 3000 |
|
|
|
|
|
| 400 | 0.1 | 90(7 × 104) | 0-0(28-0) |
| GW190519_153544 | 13.9 | 1 × 105 |
|
|
|
|
|
| 700 | 0.2 | 9(2 × 104) | 0-0(7-0) |
| GW190521 | 14.4 | 5000 |
|
|
|
|
|
| 900 | 0.4 | 4(900) | 0-0(0-0) |
| GW190521_074359 | 24.7 | 1 × 105 |
|
|
|
|
|
| 500 | 0.01 | 700(3 × 104) | 1-0(59-8) |
| GW190602_175927 | 12.6 | 1 × 105 |
|
|
|
|
|
| 700 | 0.2 | 2(2 × 104) | 0-0(6-0) |
| GW190620_030421 | 10.9 | 90 |
|
|
|
|
|
| 6000 | 2 | 200(1 × 105) | 0-0(6-0) |
| GW190630_185205 | 15.2 | 1 × 105 |
|
|
|
|
|
| 1000 | 0.03 | 9000(5 × 105) | 8-0(76-12) |
| GW190701_203306 | 11.9 | 200 |
|
|
|
|
|
| 40 | 0.003 | 9(6000) | 0-0(17- 0) |
| GW190706_222641 | 12.7 | 2 × 104 |
|
|
|
|
|
| 600 | 0.3 | 0(10) | 0-0(0-0) |
| GW190707_093326 | 13.2 | 1 × 105 |
|
|
|
|
|
| 1000 | 0.02 | 6000(2 × 105) | 20-0(86-26) |
| GW190708_232457 | 13.1 | 3000 |
|
|
|
|
|
| 1 × 104 | 0.2 | 1 × 105(4 × 106) | 10-0(79-24) |
| GW190720_000836 | 11.6 | 1 × 105 |
|
|
|
|
|
| 500 | 0.02 | 5000(2 × 105) | 12-0(81-11) |
| GW190727_060333 | 12.1 | 1 × 105 |
|
|
|
|
|
| 700 | 0.2 | 0(5000) | 0-0(2-0) |
| GW190728_064510 | 13.4 | 1 × 105 |
|
|
|
|
|
| 400 | 0.004 | 2000(1 × 105) | 10-0(79-30) |
| GW190814 | 22.2 | 1 × 105 |
|
|
|
|
|
| 20 | 9 × 10−7 | 90(200) | 76-57(100-100) |
| GW190828_063405 | 16.6 | 1 × 105 |
|
|
|
|
|
| 500 | 0.03 | 30(6 × 104) | 0-0(19-0) |
| GW190828_065509 | 11.1 | 3 × 104 |
|
|
|
|
|
| 600 | 0.03 | 100(5 × 104) | 0-0(32-1) |
| GW190910_112807 | 13.4 | 300 |
|
|
|
|
|
| 9000 | 0.9 | 6000(1 × 106) | 0-0(41-0) |
| GW190915_235702 | 13.1 | 1 × 105 |
|
|
|
|
|
| 300 | 0.02 | 700(1 × 105) | 0-0(35-1) |
| GW190924_021846 | 13.0 | 1 × 105 |
|
|
|
|
|
| 400 | 0.002 | 5000(6 × 104) | 34-1(93-55) |
| GW190929_012149 | 10.3 | 6 |
|
|
|
|
|
| 2000 | 1 | 0(1 × 104) | 0-0(2-0) |
| GW190930_133541 | 10.1 | 80 |
|
|
|
|
|
| 2000 | 0.02 | 7000(2 × 105) | 18-0(85-31) |
| GW191105_143521 | 10.7 | 80 |
|
|
|
|
|
| 800 | 0.02 | 3000(3 × 105) | 1-0(61-10) |
| GW191109_010717 | 15.8 | 6000 |
|
|
|
|
|
| 1000 | 0.3 | 5000(5 × 105) | 3-0(66-0) |
| GW191127_050227 | 10.3 | 4 |
|
|
|
|
|
| 1000 | 1 | 2(10000) | 0-0(3-0) |
| GW191129_134029 | 13.3 | 1 × 105 |
|
|
|
|
|
| 800 | 0.006 | 9000(2 × 105) | 13-0(82-35) |
| GW191204_171526 | 17.1 | 1 × 105 |
|
|
|
|
|
| 300 | 0.001 | 3000(5 × 104) | 21-0(87-52) |
| GW191215_223052 | 10.9 | 1 × 105 |
|
|
|
|
|
| 600 | 0.07 | 100(8 × 104) | 0-0(28-0) |
| GW191216_213338 | 18.6 | 1 × 105 |
|
|
|
|
|
| 500 | 4 × 10−4 | 4000(3 × 104) | 67-19(100-86) |
| GW191222_033537 | 12.0 | 1 × 105 |
|
|
|
|
|
| 2000 | 0.7 | 20(8 × 104) | 0-0(9-0) |
| GW191230_180458 | 10.3 | 20 |
|
|
|
|
|
| 800 | 0.2 | 0(40) | 0-0(0-0) |
| GW200105_162426 | 13.9 | 5 |
|
|
|
|
|
| 7000 | 0.007 | 8 × 104(3 × 105) | 80-26(100-90) |
| GW200112_155838 | 17.6 | 1 × 105 |
|
|
|
|
|
| 4000 | 0.09 | 8000(1 × 106) | 0-0(51-7) |
| GW200115_042309 | 11.5 | 1 × 105 |
|
|
|
|
|
| 700 | 0.001 | 7000(2 × 104) | 77-17(100-85) |
| GW200128_022011 | 10.1 | 200 |
|
|
|
|
|
| 2000 | 1 | 7(3 × 104) | 0-0(4-0) |
| GW200129_065458 | 26.5 | 1 × 105 |
|
|
|
|
|
| 80 | 5 × 10−4 | 400(2 × 104) | 5-0(70-27) |
| GW200202_154313 | 11.3 | 1 × 105 |
|
|
|
|
|
| 200 | 3 × 10−4 | 2000(1 × 104) | 57-8(69-17) |
| GW200208_130117 | 10.8 | 3000 |
|
|
|
|
|
| 40 | 0.004 | 0(2000) | 0-0(11-0) |
| GW200209_085452 | 10.0 | 20 |
|
|
|
|
|
| 900 | 0.4 | 4(7000) | 0-0(3-0) |
| GW200219_094415 | 10.8 | 1000 |
|
|
|
|
|
| 700 | 0.2 | 7(3000) | 0-0(2-0) |
| GW200224_222234 | 19.2 | 1 × 105 |
|
|
|
|
|
| 50 | 0.002 | 30(1 × 104) | 0-0(32-1) |
| GW200225_060421 | 13.1 | 9 × 104 |
|
|
|
|
|
| 500 | 0.01 | 2000(2 × 105) | 3-0(68-9) |
| GW200302_015811 | 10.6 | 9 |
|
|
|
|
|
| 7000 | 0.8 | 7000(1 × 106) | 0-0(50-0) |
| GW200311_115853 | 17.6 | 1 × 105 |
|
|
|
|
|
| 40 | 3 × 10−4 | 90(1 × 104) | 0-0(56-16) |
| GW200316_215756 | 10.1 | 1 × 105 |
|
|
|
|
|
| 200 | 0.005 | 500(2 × 104) | 2-0(59-8) |
Note. First, second, and third columns: GW event label, detected S/N (highest among the different pipelines; see Abbott et al. 2021e), and inverse false-alarm rate. Fourth, fifth, sixth columns: estimated primary and secondary detector frame masses and luminosity distance. Seventh, eighth, and ninth columns: primary and secondary source frame masses and redshift assuming a reference cosmology. Tenth and eleventh columns: sky localization area and 3D localization comoving volume. The twelfth column lists the number of galaxies in GLADE+ inside the localization volume for each event with observed K band (BJ band in parenthesis), while the thirteenth column reports the probability p(G∣z, H0) that GLADE+ detected the GW event host galaxy using the K-band (BJ in parenthesis) luminosity function calculated for a fiducial flat ΛCDM cosmology with H0 = 67.9 km s−1 Mpc−1 and Ωm = 0.3065. The lower and upper bounds on these probabilities are derived from the completeness fractions at the boundaries of the 90% localization volume. We report these values using a distance prior proportional to DL 2 and a uniform in the detector frame masses prior. The events in bold are the ones with S/N > 11 entering all the analyses, while the others are the ones used only for the systematic studies in Appendix B. For events released with GWTC–3, we use the first data release associated with LIGO Scientific Collaboration et al. (2021).
3.2. Description of the GLADE+ Galaxy Catalog
In one of the analyses that we perform, the redshift information is taken from galaxy surveys for all of the events apart from GW170817, for which we assume the redshift information from its EM counterpart. For the analysis taking into account galaxy surveys we use the GLADE+ (Dálya et al. 2018, 2022) all-sky galaxy catalog, which is a revised version of the first GLADE catalog (Dálya et al. 2018) containing about 22 million galaxies. GLADE+ incorporates six different galaxy catalogs and surveys, namely, the Gravitational Wave Galaxy Catalogue (GWGC; White et al. 2011), HyperLEDA (Makarov et al. 2014), the Two Micron All-Sky Survey Extended Source Catalog (2MASS XSC; Skrutskie et al. 2006), the 2MASS Photometric Redshift Catalog (2MPZ; Bilicki et al. 2014), the WISExSCOS Photometric Redshift Catalogue (WISExSCOSPZ; Bilicki et al. 2016), and the Sloan Digital Sky Survey (SDSS) quasar catalog from the 16th data release (DR16Q; Lyke et al. 2020) and covers the full sky with a completeness of about 20% up to 800 Mpc.
303
Most of the galaxies in the GLADE+ catalog have a redshift measurement obtained photometrically using an artificial neural network algorithm (Collister & Lahav 2004) with a relative error 
(Bilicki et al. 2016). The peculiar velocity corrections are implemented for galaxies up to redshifts z < 0.05 using a Bayesian technique (Mukherjee et al. 2021c) that can capture both linear and nonlinear components of the velocity field.
For our main results, we use all galaxies with measured Ks -band (denoted as K band henceforth) luminosity reported in the Vega system, and we assign a probability for each galaxy to be the host of a GW event that is proportional to this luminosity (luminosity weighting). We also explore possible systematics in our results by not using luminosity weighting and by using BJ -band observations; see Section 5.2 for more details. We choose these two bands as we have found that there is a good match between the galaxy luminosity functions and the galaxy number density of the GLADE+ catalog, in particular for the K band; see Appendix C for more details. The K-band galaxies in the GLADE+ catalog are the same as the ones in the GLADE catalog. The galaxies in the BJ band are present only in the GLADE+ catalog.
Figure 2 presents a series of sky maps showing the directional dependence of the K-band apparent magnitude threshold for the GLADE+ galaxies, in superposition with the sky localizations of the GW events included in our analysis. Outside of the galactic plane, mthr ∼ 13.5 on average for the K band, while within the galactic plane region the apparent magnitude threshold is significantly lower (i.e., brighter).
Figure 2. Sky maps showing the GLADE+ K-band apparent magnitude threshold, mthr, generated by dividing the sky into 3.35 deg2 pixels; this is the resolution used for all the events but GW190814. A mask is applied that removes from the figures all pixels with mthr < 12.5 in order to improve the figure readability. Also shown are the 90% CL sky localizations for the GW events considered in this paper.
Download figure:
Standard image High-resolution imageWe assume that the K-band absolute magnitude distribution for GLADE+ galaxies is well described by a Schechter function with parameters (reported for H0 = 100 km s−1 Mpc−1) M*,K = − 23.39 and αK
= −1.09 (Kochanek et al. 2001), while for the BJ
band we use 
and 
(Norberg et al. 2002). We set a bright cutoff high enough to include all the bright galaxies supported by the Schechter function: 
and 
. Further, we consider all the galaxies no fainter than 
. These choices correspond to all galaxies with luminosity L > 0.017L*,K and 
, where L* is the characteristic galaxy luminosity of the Schechter luminosity function. To calculate the rest-frame absolute magnitudes of the galaxies, for a given cosmology, we apply color and evolution corrections as reported in Kochanek et al. (2001) for the K band and Norberg et al. (2002) for the BJ
band.
For all events, apart from GW190814, we carry out the analysis using a pixel size of 3.35 deg2, while for GW190814 we use a pixel size of 0.2 deg2 as the sky localization for this event was 10 times smaller than most of the others. These values have been chosen taking into account the average number of galaxies per square degree reported in GLADE+. Considering the bright and faint limits of the Schechter function assumed with a median apparent magnitude threshold mthr,K = 13.5 for the K band and 
for the BJ
band, we find that there are ∼25 galaxies per square degree in GLADE+ reported in the K band and 500 galaxies per square degree reported in the BJ
band considered in the analysis. Note that the actual galaxy density per square degree is higher outside the galactic plane and in the region of GW190814.
For any given redshift, the completeness fraction, which is the probability that the galaxy catalog contains the host galaxy of the GW event, P(G∣z, H0), is defined as the fraction of galaxies with absolute magnitudes brighter than the absolute magnitude threshold (calculated from mthr), namely,
Here ϕ(L) is the assumed galaxy luminosity function, 
and 
are the minimum and maximum luminosity corresponding to 
and 
, and Lthr is the threshold luminosity for detection, calculated from mthr.
Figure 3 shows the completeness fraction of the GLADE+ catalog, in the K band and BJ band, respectively, as a function of the redshift and for different values of mthr, assuming a fiducial cosmology with H0 = 67.9 km s−1 Mpc−1 and Ωm = 0.3065 (Ade et al. 2016). As can be seen from the figure, the GLADE+ catalog is less complete in the K band than in the BJ band, but we decide to use the K-band data for our main results as they are better described by the Schechter function assumed in our analysis; see Appendix C.
Figure 3. Top: completeness fraction of GLADE+ in the K band, indicating the probability that the catalog contains the host galaxy of a GW event, as a function of the redshift for H0 = 67.9 km s−1 Mpc−1 and Ωm = 0.3065. The different lines are calculated for a given K-band apparent magnitude mthr threshold indicated in the legend. The legend also indicates the fraction of the sky, computed by dividing the sky in equally sized pixels of 3.35 deg2, for which the apparent magnitude threshold is brighter than the one reported in the legend. The fraction of pixels with no galaxies is ∼ 5% for the Bj and K bands. Bottom: same as the top panel but for the BJ band.
Download figure:
Standard image High-resolution imageFor GLADE+, the systematic uncertainties of the photometric redshift reconstruction are inside the statistical errors associated to each galaxy (Bilicki et al. 2014) with a very small percentage of outliers. We do not consider deeper galaxy surveys with a restricted sky area footprint, such as the DES Y1 survey (Drlica-Wagner et al. 2018) or DESI Legacy Imaging Survey (Zou et al. 2019) that are supposed to be complete up to redshift z ∼ 1, as we decide to employ the same all-sky galaxy catalog for all of our events. Moreover, color corrections and photometric redshift reconstruction might need particular attention with deeper galaxy surveys.
4. Results
The first analysis that we present in Section 4.1 will focus on the impact of the BBH population source masses on inference of the cosmological parameters, using the formalism discussed in Section 2.1. This analysis uses no galaxy catalog information; instead constraints on the cosmological parameters will be inferred from the mass scale set by the source mass distribution. We use only the BBHs detections as they are the majority of the sources entering in our cosmological analysis and because a joint description of NSBH, BNS, and BBHs is uncertain.
The second analysis in Section 4.2 fixes the source mass distribution and uses redshift information derived from galaxy catalogs, based on the formalism discussed in Section 2.2. Unless stated otherwise, all the figures are generated with a uniform prior on H0, credible intervals are reported as maximum posterior and 68.3% highest density intervals. We use a flat-in-log prior on H0 only when quoting results combined with GW170817 and its EM counterpart.
4.1. Implications of Population Assumptions for Cosmology
We jointly estimate population-related GW parameters and cosmological parameters using BBH events, as these are the majority of the GW events observed to date. We use the 42 BBH detected events with S/N > 11. We exclude from this analysis GW190814 (Abbott et al. 2021b) given the current uncertainty on the nature of the secondary object in this system.
We consider two cosmological models: (i) a flat w0CDM model with wide priors on the Hubble constant H0, matter density Ωm and dark energy EoS w0 parameter, and (ii) a flat ΛCDM universe with a fixed value of Ωm = 0.3065 (Ade et al. 2016) and dark energy EoS parameter w0 = −1 and with a restricted prior in the H0-tension region (H0 ∈ [65, 77] km s−1 Mpc−1). We refer to model (i) as the w0CDM model and to model (ii) as the H0-tension model. We also adopt wide priors on the hyperparameters of the GW source mass distribution and its merger rate evolution as described in Appendix A. For all the phenomenological mass models assumed we obtain posteriors on the source mass distribution and merger rate parameters, which are compatible with previous population studies (Abbott et al. 2021b, 2021f) and the latest studies with O3b data (Abbott et al. 2021e). See Appendix B for more details.
As evident from values of Bayes factor reported in Table 2, we do not find any preference of the data in supporting any one of the cosmological models (w0CDM model or the H0-tension model) considered in the analysis. As we will see later, this is because the posteriors on Ωm and w0 are not constrained by the GW observations, and the error on the H0 estimation extends beyond the H0-tension region.
Table 2. Logarithm of the Bayes Factor Comparing Runs That Adopt the Same Source Mass Model but Different Cosmologies: Wide Priors (For a General w0CDM Cosmology) vs. Restricted Priors (in the H0-tension Region)
| Mass model |
|
|---|---|
| Truncated | 0.2 |
| Power Law + Peak | − 0.3 |
| Broken Power Law | − 0.4 |
Download table as: ASCIITypeset image
In Table 3 we report the Bayes factors computed between different mass models, for the case of wide priors on the w0CDM cosmological parameters. Consistent with Abbott et al. (2021b, 2021f), we find that, even if we allow the cosmological parameters to vary with wide priors, the Truncated model is still strongly disfavored with respect to the Power Law + Peak and Broken Power Law models by a factor ∼100. This result is consistent with the fact that, as indicated in Figure 1, the source mass distribution contains more structure than a simple Truncated model. As motivated in Abbott et al. (2021b), this comparatively poor fit for the Truncated model is due to the inability of this model to capture a moderate fraction of detected events with high masses, while predicting a large fraction of detected events with lower masses. Using the reduced set of signals with S/N > 11, we do not find any compelling evidence to prefer the Power Law + Peak model over the Broken Power Law model.
Table 3. Logarithm of the Bayes Factor Between the Different Mass Models and the Power Law + Peak Model Preferred by the Data, for the Case of a w0CDM Cosmology with Wide Priors
| Mass model |
|
|---|---|
| Truncated | − 1.9 |
| Power Law + Peak | 0.0 |
| Broken Power Law | − 0.5 |
Download table as: ASCIITypeset image
The marginal posterior distributions that we obtain for the cosmological parameters H0, Ωm, and w0 are shown in Figure 4 for each phenomenological mass model. As anticipated by our Bayes factor results, we find that with the current BBH GW events we cannot constrain the values of these three cosmological parameters, as we obtain broad and uninformative posteriors.
Figure 4. Top panel: marginal posterior distribution for H0. Middle panel: marginal posterior distribution for Ωm. Bottom panel: marginal posterior distribution for w0. In each panel the different lines indicate the three phenomenological mass models. The solid orange line identifies the preferred Power Law + Peak model. The pink shaded areas identify the 68% CI of the cosmological parameters inferred from measurements from the CMB (Ade et al. 2016; apart for w0 that is reported at 95% CI), and the green shaded area in the top panel shows the value of the Hubble constant measured in the local universe (Riess et al. 2019).
Download figure:
Standard image High-resolution imageWith the Power Law + Peak we estimate 
, while for the Broken Power Law model we estimate 
. These constraints on H0, as we will see later, arise from the ability of these models to fit an excess of BBHs with masses around 35 M⊙, which sets a scale for the redshift distribution of BBHs.
We discuss this effect further using the Power Law + Peak model. Figure 5 shows the joint posterior distribution between the cosmological parameters and the parameters μg and 
defined in Equation (A11), which govern the position of the BBH Gaussian excess and the upper end of the source primary mass distribution, respectively.
Figure 5. Posterior probability density for H0 and population parameters 
and γ, governing the position of the Gaussian peak, the upper end of the mass distribution, and the merger rate evolution in the Power Law + Peak mass model. The solid and dashed black lines indicate the 50% and 90% CL contours.
Download figure:
Standard image High-resolution imageThe presence of a peak in the BBH source mass distribution allows us to set a characteristic source mass scale, which informs H(z) and allows us to exclude higher values of H0. Marginalizing over the cosmological parameters, we obtain a central value of 
for the peak position of the Gaussian BBH excess.
We also find that the estimation of population parameters related to the mass and merger rates of BBHs is not strongly dependent on the choice of cosmological priors between the w0CDM model and the H0-tension model. The only parameter showing a nonnegligible impact is μg, which is more precisely measured when using the H0-tension prior, with a posterior standard deviation that is 65% of the value when using the w0CDM model. This is due to the degeneracy observable between μg and H0 in Figure 5.
On the other hand, the disfavored Truncated model shows support at higher H0. This result is due to the fact that the Truncated model is not able to adequately fit the presence of massive binaries while producing an excess of BBHs with masses ∼ 40 M⊙ in the detector frame. For this reason, higher H0 values are more supported as those values place events at higher redshifts, thus reducing their source masses.
When we combine the H0 posteriors from the three mass models with the H0 inferred from the bright standard siren GW170817 (see Figure 6), we find a value of 
for the Power Law + Peak model and 
for the Broken Power Law model. These results represent an improvement of 17% and 12%, respectively, compared with the H0 value reported in Abbott et al. (2021a) that made use of GW170817 and six BBH detections from O2, with redshift information inferred from galaxy catalogs. For the Truncated model, we obtain 
. These results are obtained assuming a redshift-independent mass distribution. Considering a redshift dependence of the mass distribution can degrade the constraints.
Figure 6. Posterior distributions for H0 obtained by combining the H0 posteriors from the 42 BBH detections and the H0 posterior inferred from the bright standard siren GW170817. The pink and green shaded areas identify the 68% CI constraints on H0 inferred from the CMB anisotropies (Ade et al. 2016) and in the local universe from SH0ES (Riess et al. 2019), respectively.
Download figure:
Standard image High-resolution image4.2. Results Using Galaxy Catalog Information
We now discuss constraints on H0 when we fix the source population model but employ galaxy surveys to infer statistical redshift information using the pixelated gwcosmo code (Gray et al. 2020). Our analysis incorporates 47 GW events, comprising 42 BBH detections, GW190814, the 2 BNS events GW170817 and GW190425, and the 2 NSBH events GW200105 and GW200115. We include all galaxies of the GLADE+ catalog that lie inside the 99.9% estimated sky area of each event. We use the GLADE+ K-band data in this analysis, adopting luminosity weights for each galaxy. For a more in-depth discussion about the impact of our BH population assumptions and choice of photometric bands, see Section 5.2.
To describe the distribution of BH primary masses, we use a Power Law + Peak source mass model where we fix the population parameters to the median values obtained in the joint cosmological and population analysis described in Section 4.1. For the rate evolution we adopt γ = 4.59, k = 2.86, and zp
= 2.47, while for the Power Law + Peak model we use α = 3.78, β = 0.81, 
, 
, δm
= 4.8 M⊙
μg = 32.27 M⊙, σg = 3.88 M⊙, and λg = 0.03. For the NS source mass model we consider a uniform distribution between 
and 
, consistent with (Abbott et al. 2021f). We evaluate the GW selection effects using LIGO and Virgo sensitivities during the O1, O2, and O3 runs.
In Figure 7 we show the posteriors for all of the GW events considered in this analysis for the K band. For many of the O3 events, the H0 inference is dominated by the likelihood based on the hypothesis that the host galaxy is not in the catalog (referred to as out-of-catalog). The out-of-catalog term dominates for sources that are localized at redshifts at which the GLADE+ galaxy catalog has a low completeness fraction (see Figure 3). This is the case for most of the GW sources that are BBHs observed at large luminosity distances. Another interesting trend observed in Figure 7 is that, for lower values of H0, the in-catalog likelihood terms tend to dominate because for low H0 values the GW events are placed at smaller redshifts where the galaxy catalog is more complete, as shown in Figure 3.
Figure 7. Plots reporting the results on the H0 inference for each event using the GLADE+
K band and luminosity weighting. Two panels are shown for each event. Top panels: hierarchical likelihood under the hypothesis, G, that the host galaxy is in the catalog (blue solid lines) and under the hypothesis, 
, that the host galaxy is not in the catalog (pink dashed lines). The different lines shown in each panel correspond to the different pixels within the sky localization area for each event. Bottom panels: the blue solid line shows the posterior obtained by summing the terms corresponding to the in-catalog and out-of-catalog hypotheses. The orange dashed line shows the posterior obtained by assuming a galaxy catalog with null completeness. In this case the H0 inference comes entirely from the population assumptions. This plot is intended to show which event is informative on the H0 value and whether the information is coming from population assumptions or galaxy catalog contribution.
Download figure:
Standard image High-resolution imageFor most of these events, the number of galaxies present in the sky localization volume is large enough that the redshift information is still dominated by population assumptions (Section 5.2). GW190814 is the only event for which there is a sufficiently small number of galaxies in its sky localization area of about 18 deg2. Its small area makes this event partially more informative on the value of H0 in comparison to the other GW events. We can see in Figure 7 that, out of all the GW events, the most informative posterior on H0 (compared to the zero galaxy catalog completeness posterior) is from GW190814, provided that the luminosity weighting scheme is applied. We have verified that the H0 posterior with the K band and using luminosity weights does not depend on the faint end magnitude limit used for the analysis. For this event, we infer an H0 constraint of 
(MAP and HDI). We quote the maximum a posteriori probability (MAP) and the corresponding highest density interval (HDI) values in the analysis.
Figure 8 shows the redshift distribution of galaxies in the 90% CI sky area of GW190814 (top panel) and the galaxy catalog completeness (bottom panel), compared to the predicted distribution for a prior that is uniform in the comoving volume. We observe that for the K band the H0 support results from an excess of galaxies, with respect to the uniform in the comoving volume prior, around z ∼ 0.051. Switching off the luminosity weighting assumption decreases the contribution of this excess of galaxies as the completeness is estimated to be lower. The same excess is not visible in the BJ band as more galaxies are reported in this band and some luminous galaxies with measured K-band apparent magnitudes do not have measured apparent magnitudes for the BJ band.
Figure 8. Top panels: distribution of galaxies observed in the GLADE+ K (left panels) and BJ (right panels) bands as a function of the redshift z with and without galaxy luminosity weights, compared with the GW190814 redshift localization in its 90% CI sky area, assuming a cosmology with H0 = 67.9 km s−1 Mpc−1 and Ωm = 0.3065 (green line) and the predicted redshift distribution for a prior that is uniform in the comoving volume (blue dashed line). The redshift distribution is not intended as representative of the reconstructed galaxy distribution as it is calculated by stacking each galaxy redshift localization assumed as a normal distribution. This procedure only serves to give a rough idea where the H0 contribution is coming from. Bottom panels: completeness calculated using the K and BJ band with and without application of the luminosity weighing scheme.
Download figure:
Standard image High-resolution imageDespite the cases where there is a significant in-catalog contribution, the final H0 result is nevertheless dominated by the BBHs population assumptions that are contributing to the out-of-catalog likelihood terms (when the galaxy catalog is not complete) and in the in-catalog terms when a large number of galaxies is present in the GW sky localization volume.
In Figure 9 we show the combined H0 posterior inferred from all of the GW events and for several different scenarios. By using all of the dark sirens, together with K-band galaxy information from GLADE+, we obtain a value of 
. This value is strongly dominated by the BH population assumptions, as can be seen in Figure 9. The H0 value obtained from population assumptions alone (empty catalog case in Figure 9) is 
. When we combine the galaxy catalog measurement with the result from the bright standard siren GW170817, we obtain 
. This value represents an improvement of 42% with respect to the corresponding result obtained with GWTC–1 (Abbott et al. 2021a), and an improvement of 44% with respect to the result of 
obtained using only the GW170817 event (Abbott et al. 2021a).
Figure 9. Hubble constant posterior for several cases. Gray dotted line: posterior obtained using all dark standard sirens without any galaxy catalog information and fixing the BBH population model. Orange dashed line: posterior using all dark standard sirens with GLADE+ K-band galaxy catalog information and fixed population assumptions. Black solid line: posterior from GW170817 and its EM counterpart. Blue solid line: posterior combining dark standard sirens and GLADE+ K-band catalog information (orange dashed line) with GW170817 and its EM counterpart (black solid line). The pink and green shaded areas identify the 68% CI constraints on H0 inferred from the CMB anisotropies (Ade et al. 2016) and in the local universe from SH0ES (Riess et al. 2019), respectively.
Download figure:
Standard image High-resolution image5. Discussion
5.1. Considerations for The BBH-based Population Analysis
We have shown how population assumptions on BBH formation dominate the inference on cosmological parameters and, in particular, we have seen how the presence of an excess of BBHs with primary masses between 30 M⊙ and 40 M⊙ (Farr et al. 2019) sets a scale for the BBH redshifts, thus allowing for a weak constraint on H0.
In the BBH-based population analysis without GW170817, the (H0, Ωm, w0) parameters are not constrained. In Figure 10 we portray these constraints on the expansion rate of the universe, H(z). The best constraint that we obtain on the expansion rate of the universe has a value, and uncertainty (median and symmetric 90% CI), of 
at redshift z = 0 if we include the bright siren GW170817, and 
at redshift z ∼ 0.01 without it.
Figure 10. Evolution of the Hubble parameter predicted from the most preferred mass model Power Law + Peak (blue lines). The yellow shaded area indicates the 90% CL contours identified by the uniform priors on H0, Ωm, and w0 while the blue shaded area indicates the 90% CL contours from the posterior of the preferred mass model. The dashed lines indicate the median of the prior and posterior for H(z).
Download figure:
Standard image High-resolution image5.2. Considerations on The Catalog Analysis
As already discussed in Section 4.2, the H0 inference is dominated by the population assumptions of the underlying BH mass distribution. In particular, as shown in Figure 5, the population parameter that is most strongly correlated with the value of H0 is the position of the BHs excess μg.
In Figure 11 we show the H0 posterior computed with different choices of μg and fixing the remaining parameters. The values of μg are spaced by ∼ 2.5M⊙, which is roughly the uncertainty identified in Section 2.1. It can be seen that the value of μg has a strong impact on the inference of H0. For values of μg higher than the median value of 32.55M⊙, the posterior supports low H0 values as we need GW events to be at a lower redshift to explain the excess of BHs at higher masses. On the other hand, for μg < 32.55 M⊙, the posterior supports higher H0 value in order to place events at higher redshift compatible with an excess of BHs at lower masses.
Figure 11. Systematic effects on the inference of the Hubble constant due to the choice of different values for the mean μg of the Gaussian component in the source mass model, and other population model parameters (upper panel) and different choices for the luminosity band and weighting scheme adopted for the GLADE+ galaxy catalog (lower panel). The pink and green shaded areas identify the 68% CI constraints on H0 inferred from the CMB anisotropies (Ade et al. 2016) and in the local universe from SH0ES (Riess et al. 2019), respectively.
Download figure:
Standard image High-resolution imageWe also explored the effect of raising the maximum mass of the black holes to 
. As can be seen in Figure 11, raising 
does not have a significant effect on the H0 posterior as only few events are present at these masses. One more parameter for which we explored the effect of its variation is the γ parameter in the rate evolution model. In the same plot one can see the H0 posterior for γ = 2.59. This parameter has a stronger effect on the H0 posterior, making the posterior less informative and at the same time moving its peak to higher values.
The galaxy catalog brings additional information only for GW190814, due to the much better sky localization (∼18 deg2) for this event; this has the effect of providing more support for the H0-tension region.
In Figure 12, we show how population assumptions impact the hierarchical likelihood calculation as a function of H0, for the hypotheses that the host galaxy is (or is not) inside the catalog. Population assumptions strongly impact the out-of-catalog term of the likelihood, which is the dominant contribution to the H0 posterior when the event is localized in an area where the galaxy catalog has a low completeness fraction, which happens for most of the GW events that we consider in this analysis. On the other hand, population assumptions are less important for events with a small localization in a region of the galaxy catalog that is complete. In these cases (for example GW190814), the posterior is dominated by the in-catalog likelihood terms and hence exhibits a weak dependence on the population assumptions.
Figure 12. Plots showing, for each event, the hierarchical likelihood as a function of H0 marginalized over sky localization as a function of different population assumptions and using the GLADE+ K band. The y-axes of each panel start from zero and have the same normalization. Two panels are shown for each event. Top panels: hierarchical likelihood under the hypothesis that the host galaxy is not in the catalog, for several different population assumptions. Bottom panels: hierarchical likelihood under the hypothesis that the host galaxy is in the catalog, again for several different population assumptions. This shows how population and galaxy catalog treatment affect the information on the H0 value for each event.
Download figure:
Standard image High-resolution imageThe aforementioned discussion also explains the difference between our results and those found in Finke et al. (2021) using GWTC–2 events. In that work the authors found a weak dependence of their posterior on the source population parameters. However, in their case only a few events, above a given completeness threshold of 70% for the main result, were used to explore systematic effects due to the source population. Moreover, in exploring these systematics they varied the population assumptions only within the range of uncertainties reported in Abbott et al. (2021b), which already assumes a fixed cosmological model with a value of H0 consistent with the Planck results (Ade et al. 2016). Consequently, the results obtained by them are primarily driven by the Planck cosmological parameters.
Finally, we also explore the systematics introduced by choices related to the galaxy catalog data. In Figure 11 we also show the H0 posteriors obtained with the GLADE+ catalog, but using K-band galaxies without luminosity weighting and BJ -band galaxies with luminosity weighting. In both cases, the H0 posterior is not significantly affected by this choice, and it is, again, dominated by the population assumptions.
However, the impact of using luminosity weights is not negligible. For instance, in the case of GW190814 (see Figure 12), removing the extra luminosity weight suppresses the H0 posterior peak around 70 km s−1 Mpc−1. This arises because the luminous galaxies shown in Figure 8, observed in the K band, are now contributing to the GW event redshift localization with the same probability as the other 200+ galaxies included in the GW localization volume. We have verified that the H0 posterior with the K band and using luminosity weights do not depend on the faint end magnitude limit used for the analysis.
5.3. Additional Systematic Uncertainties
As we have seen, the presence of a known source mass scale can be used to measure cosmological parameters (Farr et al. 2019; Mastrogiovanni et al. 2021). However, if the source mass distribution is mismodeled, then the cosmological inference will be biased. With the current set of events, this effect contributes the dominant source of systematic uncertainty in the measurement of H(z). In the population-based method, a key assumption is that the source mass distribution does not evolve with the redshift (Fishbach et al. 2021); any evolution is degenerate with the cosmological inference. Many BBH formation scenarios predict mild evolution in the mass distribution (Mapelli et al. 2019; Weatherford et al. 2021; Gallegos-Garcia et al. 2021; Mapelli et al. 2022; van Son et al. 2022), but given our broad statistical uncertainties, we expect this evolution to only weakly affect our results, and we do not attempt to calibrate the mass models to theory. In the galaxy-catalog-based method, we fix the source mass distribution. In this case, the choice of the peak location μg, associated with excess BHs in the distribution of source masses, is a main source of systematic uncertainty. If the μg parameter is assumed to be lower (or higher) than its true value, this can lead to a higher (or lower) inferred value of H0. The impact of this bias can be reduced if the support from the in-catalog part of the statistical galaxy catalog method is informative, which is mainly possible for sources with better three-dimensional localization error and with a more complete galaxy catalog. In the future, as more GW detectors join the network, resulting in more events with better sky localizations and galaxy catalogs that are more complete at high redshift, the impact of the source population on the galaxy catalog method can be mitigated.
One of the additional sources of contamination, when seeking to infer the true luminosity distance DL (and hence the true source masses) of a GW source in the absence of an EM counterpart, is the possible lensing of the GW signal due to the intervening matter distribution (Schneider et al. 1992; Bartelmann 2010; Nakamura 1998; Wang et al. 1996; Takahashi & Nakamura 2003; Dai et al. 2017; Broadhurst et al. 2018; Diego 2019). In the geometric optics limit, lensing modifies the GW signal by a magnification factor μ that only changes the amplitude of the GW strain, leading to a measured luminosity distance given by 
. In the strong lensing limit, when the value of μ is large, the inferred luminosity distance to the source 
may be substantially lower than the true luminosity distance, i.e., introducing a bias in the measurement of the luminosity distance. However, for the GW detections considered here, the probability of observing such a strongly lensed event is less than one percent (Ng et al. 2018a; Oguri 2018; Mukherjee et al. 2021a), even for a broad range of astrophysical time-delay scenarios (Mukherjee et al. 2021b). Moreover, searches from O1+O2 (Hannuksela et al. 2019) and more recent LIGO–Virgo analysis of the O3a data did not reveal any signs of strong lensing (Abbott et al. 2021f). Furthermore, searches for a stochastic GW background also provide a model-independent bound on the lensing event rate of BBHs (Mukherjee et al. 2021a; Buscicchio et al. 2020; Abbott et al. 2021f), which is again consistent with a low probability of contamination due to strong lensing in the GW sources considered here. In this paper we, therefore, ignore any possible impact of strong lensing on our cosmological parameter estimates.
Apart from strong lensing, weak lensing of GW sources can also be a potential source of contamination. However, due to the effects of sky averaging, weak lensing should not produce a bias in the inferred values of the cosmological parameters, but will introduce additional variance on the luminosity distance of individual sirens at the level of a few percent (Holz & Wald 1998; Hirata et al. 2010) and it is subdominant in comparison to the intrinsic measurement error (of about 20%) on the luminosity distance even for the best source in our current GW catalog. Hence, we also ignore the uncertainty due to weak lensing in this paper. However, in a future analysis, we could include the contribution from weak lensing and its impact on the distance measurement (Holz & Wald 1998; Cutler & Holz 2009; Hirata et al. 2010; Namikawa et al. 2016; Mukherjee et al. 2020).
6. Conclusions
Using the 47 GW events with detected S/N > 11 reported in the third LIGO–Virgo–KAGRA Gravitational Wave Transient Catalog (Abbott et al. 2021e), we have inferred constraints on the cosmological parameters adopting two different approaches: hierarchical inference without galaxy surveys and the statistical galaxy catalog method. We present for the first time analysis that constrains jointly the properties of the population of BBHs and the parameters of the cosmological model, and we have shown the crucial correlation that exists between the two sectors.
We have shown that an excess population of BHs in the mass range 30–40 M⊙, pointed out by Abbott et al. (2021b), is robust to the choice of assumed cosmological model parameters. While our constraints on the present-day matter density, Ωm, and dark energy EoS, w0, parameters are weak, we have measured the Hubble constant to be 
at 68% CL from combining dark sirens with information from the bright siren GW170817 and its EM counterpart (Abbott et al. 2019b). This result represents an improvement of 17% with respect to the H0 value reported from analysis of O1 and O2 data (Abbott et al. 2021a) that made use of galaxy catalogs alone to infer statistical redshift information. In our analysis we also obtain weak constraints on the expansion history as a function of the redshift.
In addition we provide a constraint on the value of H0 adopting a fixed Power Law + Peak population model of BBHs and using statistical redshift information inferred from the GLADE+ galaxy catalog. This analysis obtained, for the K band, 
, which represents an improvement of 42% with respect to Abbott et al. (2021a) alone, and an improvement of 20% with respect to recent H0 studies using GWTC–2 events (Finke et al. 2021). Most of the constraining power in our H0 inference comes from the event GW170817 using its EM counterpart. Combining the above result with information from GW170817 we obtain 
. The most informative dark siren in the GWTC–3 catalog is GW190814, which alone provides an estimate of 
, provided that the luminosity weighting scheme is applied.
A summary of the different H0 values obtained using different data sets and model assumptions can be seen in Table 4. The table is divided into two parts. The first part summarizes the values that we infer for H0 when fixing the population model to the most favorable one and then varying the luminosity band from GLADE+ used in our analysis. We used both the BJ band and the K band, and the results are very similar (see also Figure 11). The second part of the table summarizes the results obtained by marginalizing over the population parameters and using no galaxy catalog information.
Table 4. Values of the Hubble Constant Obtained in this Study Using Different Data Sets and Analysis Methods
| Description | Galaxy catalog | BBH mass model |
|
|
|---|---|---|---|---|
| [km s−1 Mpc−1] | [km s−1 Mpc−1] | |||
| No galaxy catalog, marginalizing over population model, 42 events | - | Truncated |
|
|
| - | Power Law + Peak |
|
| |
| - | Broken Power Law |
|
| |
| Using galaxy catalog, fixed population model, 47 events | GLADE+ K band | Power Law + Peak |
|
|
| GLADE+ BJ band | Power Law + Peak |
|
|
Note. The columns are in order: short description of the sources used in the study with S/N > 11, galaxy catalog used (where appropriate), BBH mass model used, and the 68.3% CL H0 value. The last two columns report the median and symmetric 90% CI H0 values. The values in the parentheses are those obtained after combining with the GW170817 EM counterpart posterior.
Download table as: ASCIITypeset image
Although we have improved our previously reported constraints on the value of H0 using these 47 GW events, our results are still dominated by the systematic effects induced by the assumptions made about the GW source population. The choice of mass scale set by μg, the mass at which the excess of BHs is centered, plays a crucial role in constraining the value of H0.
In the future, with significantly more both bright and dark sirens it will be possible to make robust measurements of H0 and other cosmological parameters. On the one hand, measurement of bright sirens will greatly help with inferring the redshift from direct observations of EM counterparts (Holz & Hughes 2005; Dalal et al. 2006; Nissanke et al. 2010; Chen et al. 2018; Feeney et al. 2019). On the other hand, for dark sirens, the application of cross-correlation techniques to infer the clustering redshift of GW sources (Mukherjee et al. 2021; Bera et al. 2020) using spectroscopic galaxy surveys (Diaz & Mukherjee 2022), the pair-instability supernovae mass scale of black holes (Farr et al. 2019; Mastrogiovanni et al. 2021), the redshift distribution of the GW sources (Ding et al. 2019; Ye & Fishbach 2021), and the tidal distortion of neutron stars (Messenger & Read 2012; Chatterjee et al. 2021) will enable robust measurement of the cosmic expansion history. With the aid of more observations and further development of analysis techniques, we will be able to reduce the current systematics and proceed toward accurate and precision gravitational wave cosmology.
This material is based upon work supported by NSF's LIGO Laboratory, which is a major facility fully funded by the National Science Foundation. The authors also gratefully acknowledge the support of the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS), and the Netherlands Organization for Scientific Research (NWO), for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific and Industrial Research of India, the Department of Science and Technology, India, the Science & Engineering Research Board (SERB), India, the Ministry of Human Resource Development, India, the Spanish Agencia Estatal de Investigación (AEI), the Spanish Ministerio de Ciencia e Innovación and Ministerio de Universidades, the Conselleria de Fons Europeus, Universitat i Cultura and the Direcció General de Política Universitaria i Recerca del Govern de les Illes Balears, the Conselleria d'Innovació Universitats, Ciència i Societat Digital de la Generalitat Valenciana and the CERCA Programme Generalitat de Catalunya, Spain, the National Science Centre of Poland and the European Union–European Regional Development Fund, Foundation for Polish Science (FNP), the Swiss National Science Foundation (SNSF), the Russian Foundation for Basic Research, the Russian Science Foundation, the European Commission, the European Social Funds (ESF), the European Regional Development Funds (ERDF), the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the French Lyon Institute of Origins (LIO), the Belgian Fonds de la Recherche Scientifique (FRS-FNRS), Actions de Recherche Concertées (ARC) and Fonds Wetenschappelijk Onderzoek–Vlaanderen (FWO), Belgium, the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, the Natural Science and Engineering Research Council Canada, Canadian Foundation for Innovation (CFI), the Brazilian Ministry of Science, Technology, and Innovations, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Leverhulme Trust, the Research Corporation, the Ministry of Science and Technology (MOST), Taiwan, the United States Department of Energy, and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, INFN, and CNRS for provision of computational resources.
This work was supported by MEXT, JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, JSPS Grant-in-Aid for Scientific Research on Innovative Areas 2905: JP17H06358, JP17H06361, and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grant-in-Aid for Scientific Research (S) 17H06133 and 20H05639, JSPS Grant-in-Aid for Transformative Research Areas (A) 20A203: JP20H05854, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, National Research Foundation (NRF) and Computing Infrastructure Project of KISTI-GSDC in Korea, Academia Sinica (AS), AS Grid Center (ASGC), and the Ministry of Science and Technology (MoST) in Taiwan under grants including AS-CDA-105-M06, Advanced Technology Center (ATC) of NAOJ, Mechanical Engineering Center of KEK.
We would like to thank all of the essential workers who put their health at risk during the COVID-19 pandemic, without whom we would not have been able to complete this work.
Appendix A: Population Prior Models
A.1. Models for Background Cosmologies
We use a flat ΛCDM cosmological model with dark energy density as a function of the redshift z described by Linder (2003)
where ρΛ is the dark energy density, and w0 is a phenomenological parameter. If the dark energy density is constant during the cosmic expansion, as it is in the standard cosmological model, then w0 = −1. The luminosity distance is calculated as
where Ωm and ΩΛ,0 are the present-day dimensionless matter and dark energy densities, respectively, and ΩΛ,0 = 1 − Ωm.
We consider two sets of priors for the cosmological background model, which are indicated in Table 5. For the first set of priors, we adopt a flat ΛCDM cosmology that restricts the Hubble constant to only the range compatible with the H0 tension, while for the second set we adopt more general, wide priors.
Table 5. Summary of the Priors on the Cosmological Parameters for the Two Sets of Priors Considered
| Restricted priors (H0 tension) | ||
|---|---|---|
| Parameter | Description | Prior |
| H0 | Hubble constant expressed in km s−1 Mpc−1 in the H0-tension region |
 (65, 77) |
| Ωm | Present-day matter density of the universe fixed to the mean value inferred from measurements of the CMB in Ade et al. (2016) | 0.3065 |
| w0 | Dark energy equation-of-state parameter fixed to the value that corresponds to constant density | -1 |
| Wide priors | ||
| Parameter | Description | Prior |
| H0 | Hubble constant expressed in km s−1 Mpc−1 |
 (10.0, 200.0) |
| Ωm | Present-day matter density of the universe |
 (0.0, 1.0) |
| w0 | Dark energy equation-of-state parameter |
 (− 3.0, 0.0) |
Download table as: ASCIITypeset image
A.2. Merger Rate and Redshift Distribution Priors
We model the binary merger rate using a phenomenological model introduced with the form of Madau & Dickinson (2014), motivated by the fact that the binary formation rate might follow the star formation rate. The parameterization that we use Callister et al. (2020) for the merger rate in the detector frame is
where R0 is the binary merger rate today, Vc is the comoving volume, γ and k are the slopes of the two power-law regimes before and after a turning point zp, and Λc are a set of parameters describing the cosmological expansion. The extra 1/(1 + z) factor encodes the clock difference in the source and detector, while the factor (1 + zp) ensures that today the merger rate is R(z = 0) = R0. The redshift prior, normalized over the redshift, can be expressed as
where C is the normalization factor calculated from Equation (A3).
The priors that we use on the merger rate hyperparameters are indicated in Table 6. The prior ranges that we model are wide enough to include the effect of a possible time delay between the formation and the merger of the binary.
Table 6. Summary of the Prior Hyperparameters Used for the Merger Rate Evolution Models Adopted in This Paper
| Parameter | Description | Prior |
|---|---|---|
| R0 | BBH merger rate today in Gpc−3 yr−1 |
 (0.0, 100.0) |
| γ | Slope of the power-law regime for the rate evolution before the point zp |
 (0.0, 12.0) |
| k | Slope of the power-law regime for the rate evolution after the point zp |
 (0.0, 6.0) |
| zp | Redshift turning point between the power-law regimes with γ and k |
 (0.0, 4.0) |
Download table as: ASCIITypeset image
A.3. Phenomenological Mass Priors
The three phenomenological mass models that we implement are a superposition of two probability density distributions and are compatible with the phenomenological priors used for BBHs in Abbott et al. (2021b, 2021f), although the prior ranges on the population parameters are different. The first is a truncated power law 
described by slope α, and lower and upper bounds 
at which there is a hard cutoff
The second is a Gaussian distribution with mean μ and standard deviation σ,
The source mass priors for the BBHs population that we consider are factorized as
where π(m1∣Φm
) is the distribution of the primary mass component while π(m2∣m1, Φm
) is the distribution of the secondary mass component given the primary. For all of the mass models, the secondary mass component m2 is described with a truncated power law with slope β between a minimum mass 
and a maximum mass m1
while the primary mass is described with several models discussed in the following paragraphs.
For some of the phenomenological models, we also apply a smoothing factor to the lower end of the mass distribution
where S is a sigmoid-like window function that adds a tapering of the lower end of the mass distribution. See Equations (B6) and (B7) of Abbott et al. (2021b) for the explicit expression for the window function.
The three phenomenological mass models are highlighted in the following list. In Table 7, we report the prior ranges used for the population hyperparameters.
- 1.Truncated model: It describes the distribution of the primary mass m1 with a truncated power law with slope − α between a minimum mass
and a maximum mass . - 2.Power Law + Peak model: It describes the primary mass component as a superposition of a truncated PL, with slope − α between a minimum mass and a maximum mass , plus a Gaussian component with mean μg and standard deviation σg,
- 3.Broken Power Law model: It describes the distribution of m1 as a PL between a minimum mass and a maximum mass . The Broken Power Law model is characterized by two PL slopes, α1 and α2, and by a break point between the two regimes at , where b is a number ∈ [0, 1]. The broken PL model is
Table 7. Summary of the Priors Used for the Population Hyperparameters for the Three Phenomenological Mass Models
| Truncated | ||
|---|---|---|
| Parameter | Description | Prior |
| α | Spectral index for the PL of the primary mass distribution |
 (1.5, 12.0) |
| β | Spectral index for the PL of the mass ratio distribution |
 (− 4.0, 12.0) |
| Minimum mass of the PL component of the primary mass distribution |
 (2.0 M⊙, 10.0 M⊙) |
| Maximum mass of the PL component of the primary mass distribution |
 (50.0M⊙, 200.0M⊙) |
| Power Law + Peak | ||
| Parameter | Description | Prior |
| α | Spectral index for the PL of the primary mass distribution |
 (1.5, 12.0) |
| β | Spectral index for the PL of the mass ratio distribution |
 (− 4.0, 12.0) |
| Minimum mass of the PL component of the primary mass distribution |
 (2.0M⊙, 10.0M⊙) |
| Maximum mass of the PL component of the primary mass distribution |
 (50.0M⊙, 200.0M⊙) |
| λg | Fraction of the model in the Gaussian component |
 (0.0, 1.0) |
| μg | Mean of the Gaussian component in the primary mass distribution |
 (20.0M⊙, 50.0M⊙) |
| σg | Width of the Gaussian component in the primary mass distribution |
 (0.4M⊙, 10.0M⊙) |
| δm | Range of mass tapering at the lower end of the mass distribution |
 (0.0M⊙, 10.0M⊙) |
| Broken Power Law | ||
| Parameter | Description | Prior |
| α1 | PL slope of the primary mass distribution for masses below mbreak |
 (1.5, 12.0) |
| α2 | PL slope for the primary mass distribution for masses above mbreak |
 (1.5, 12.0) |
| β | Spectral index for the PL of the mass ratio distribution |
 (− 4.0, 12.0) |
| Minimum mass of the PL component of the primary mass distribution |
 (2.0M⊙, 50.0M⊙) |
| mmax | Maximum mass of the primary mass distribution |
 (50.0M⊙, 200.0M⊙) |
| b | Fraction of the way between  and  at which the primary mass distribution breaks |
 (0.0, 1.0) |
| δm | Range of mass tapering on the lower end of the mass distribution |
 (0.0M⊙, 10.0M⊙) |
Download table as: ASCIITypeset image
Appendix B: Full Results of the Population Analysis and Effect of Different S/N Cuts
We provide extra details on the joint inference of cosmological and population parameters using BBHs.
In Figure 13 we show the corner plots for the posterior associated with the Power Law + Peak model (the one implemented for the galaxy catalog analysis) adopting wide priors on the cosmological parameters and for the population of BBHs. As was also shown in Figure 5, the main mass-related population parameters correlating with the cosmological parameters, and in particular H0, are the position of the Gaussian component (BBH excess in the mass distribution) and the higher end of the source mass distribution.
Figure 13. Corner plots for the preferred Power Law + Peak model parameters and cosmological parameters, fitted to BBHs with S/N > 11.
Download figure:
Standard image High-resolution imageThe other parameter that correlates with the estimation of H0 is the rate evolution parameter γ. This parameter models a power-law-increasing merger rate with the redshift. We find that higher values of γ support lower values of H0, which is due to the fact that lowering H0 will place events at lower redshifts that are incompatible with the observed mass distribution; therefore γ tries to correct this by supporting higher redshifts. However, the posterior on the rate evolution is well within the statistical uncertainties given in Abbott et al. (2021f).
We run additional systematic studies using different S/N cuts for the selection of the BBHs to include in our analysis. We explore a higher S/N cut of 12 (more pure) that selects a sample of 35 events. We also explore a lower S/N cut of 10, allowing 60 events with a IFAR > 0.5 yr and no plausible instrumental origin. For this set of events, GW190426_152155 and GW190531_023648 are excluded as their secondary mass extends to lower masses in the NS region.
The marginal H0 posterior for all the mass models is shown in Figure 14, where we show that including more events always produces a posterior on the H0 within the statistical uncertainties of other selection criteria. In all of these cases, the excess of BBHs around 35M⊙ is present for all the S/N cuts, and it is responsible for the preference observed in the H0 posterior.
Figure 14. Marginal posterior probability distributions for H0 with using BBH events with three different S/N cuts and the three mass models.
Download figure:
Standard image High-resolution imageAppendix C: Schechter Luminosity Function Studies
In this appendix we show comparisons of the Schechter luminosity function (LF) for K- and BJ -band galaxies reported in GLADE+. Wrong assumptions on the LF, or incorrect description of the selection biases, could potentially introduce a bias in the inferred H0. Indeed, one of the key assumptions that we have made to construct our completeness corrections is that the galaxy catalog is magnitude limited, i.e., galaxies are not detected only because they are too faint. This cannot be the case if another selection bias (based on, e.g., colors or spectral features) were present, or if color and evolution corrections are not implemented properly.
In Figure 15 we show a comparison of the assumed LF and the number density of galaxies per comoving volume present in GLADE+. In the case that the galaxy catalog can be correctly described as magnitude-limited, we expect that the distribution of the GLADE+ galaxies distribution will match the assumed LF at its bright end, and then will start to decrease when we reach (and exceed) the corresponding absolute magnitude threshold. Galaxies in the K band are well described by this behavior and missing galaxies can be explained by the impact of the apparent magnitude threshold, while for the BJ band there seems to be some additional missing galaxies at low redshift. This observed behavior motivated our decision to present our main results using the K-band magnitudes compiled in the GLADE+ catalog.
Figure 15. Left: comparison of the assumed BJ -band galaxy LF (pink dashed line) and the differential number density of galaxies in different redshift bins (solid colored lines). The vertical dashed line indicates the median absolute magnitude threshold for galaxy detection, as computed in our code. Right: same comparison, but for the K-band LF. In each panel the solid lines are calculated based on the median apparent magnitude threshold, computed over all sky directions, then collecting and making a histogram of all galaxies in the GLADE+ catalog that lie in the appropriate redshift bin and are brighter than the apparent magnitude threshold.
Download figure:
Standard image High-resolution image
Footnotes
- 301
The data distribution files for these methods are available in Zenodo under a Creative Commons (Attribution) license:10.5281/zenodo.5645777.
- 302
Events in the analysis showing differences in posterior samples with different waveforms are GW191109_010717 and GW200129_065458 (Abbott et al. 2021e). These differences are mostly related to the effective and precession spin parameters, which are not considered in this analysis.
- 303
Links to the different constituent galaxy catalogs and surveys in GLADE+ are as follows: GWGC- http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=GWGC, HyperLEDA-http://leda.univ-lyon1.fr/, 2MASS XSC-https://wise2.ipac.caltech.edu/staff/jarrett/2mass/XSC/, 2MPZ-http://ssa.roe.ac.uk/TWOMPZ.html, WISExSCOS-http://ssa.roe.ac.uk/WISExSCOS.html, SDSS-DR16Q-https://www.sdss.org/dr16/algorithms/qso_catalog/.