LSU Digital Commons LSU Digital Commons A Gravitational-wave Measurement of the Hubble Constant A Gravitational-wave Measurement of the Hubble Constant following the Second Observing Run of Advanced LIGO and Virgo following the Second Observing Run of Advanced LIGO and Virgo

This paper presents the gravitational-wave measurement of the Hubble constant ( H 0 ) using the detections from the ﬁ rst and second observing runs of the Advanced LIGO and Virgo detector network. The presence of the transient electromagnetic counterpart of the binary neutron star GW170817 led to the ﬁ rst standard-siren measurement of H 0 . Here we additionally use binary black hole detections in conjunction with galaxy catalogs and report a joint measurement. Our updated measurement is H 0 = -+ 69 816 km s − 1 Mpc − 1 ( 68.3% of the highest density posterior interval with a ﬂ at-in-log prior ) which is an improvement by a factor of 1.04 ( about 4% ) over the GW170817-only value of -+ 69 817 km s − 1 Mpc − 1 . A signi ﬁ cant additional contribution currently comes from GW170814, a loud and well-localized detection from a part of the sky thoroughly covered by the Dark Energy Survey. With numerous detections anticipated over the upcoming years, an exhaustive understanding of other systematic effects are also going to become increasingly important. These results establish the path to cosmology using gravitational-wave observations with and without transient electromagnetic counterparts.


Introduction
Gravitational waves (GWs) from compact binary coalescences allow for the direct measurement of the luminosity distance to their source. This makes them standard-distance indicators, and in conjunction with an identified host galaxy or a set of possible host galaxies, they can be used as standard sirens to construct a redshift-distance relationship and measure cosmological parameters like the Hubble constant (H 0 ; Schutz 1986;Holz & Hughes 2005;MacLeod & Hogan 2008;Nissanke et al. 2010;Sathyaprakash et al. 2010). The GW signature from the binary neutron star (BNS) merger GW170817, along with its coincident electromagnetic (EM) transient associated with the host galaxy NGC4993, led to a first standard-siren measurement of H 0 ). This measurement is independent of other state-of-the-art measurements of H 0 , and in particular, independent of the cosmic distance ladder used to calibrate standardizable sources like SNe Ia (Type Ia supernovae). The importance of an independent measurement of H 0 is worth highlighting. With the Planck 2018 data release (Planck Collaboration et al. 2020), and the recalibration of supernovae using Large Magellanic Cloud Cepheids (Riess et al. 2019), the tension between early universe measurements of H 0 from Planck and local measurements from the Supernovae H0 for the Equation of State (SH0ES) project has risen to the 4.4σ level. Independent measurements using cosmological baryon acoustic oscillations to calibrate SNe Ia via the inverse distance ladder (Macaulay et al. 2019) and gravitational lensing of quasars in the nearby universe (H0LiCOW Collaboration, Birrer et al. 2019) favor to some degree the early universe Planck and the local SH0ES measurements, respectively. A complementary measurement of H 0 from the multimessenger GW astronomy sector would help clarify whether the current tension is a statistical anomaly or evidence for new physics beyond the Lambda cold dark matter (ΛCDM) model of cosmology. 196 The GW standard-siren measurement in Abbott et al. (2017a) is broadly consistent with other measurements. By combining information from multiple detections, one can improve the accuracy reaching about 1% with ( )  100 detections in the coming years (Dalal et al. 2006;Nissanke et al. 2013;Chen et al. 2018;Vitale & Chen 2018;Feeney et al. 2019;Mortlock et al. 2019).
An unambiguous identification of the host galaxy is unlikely for all BNS detections; only a crude estimate of the sky position may be available. Moreover there are sources such as 195 Deceased, 2018 July.
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. 196 Cosmological parameters can potentially be inferred from GW observations alone by estimating the redshift using the known physics of neutron stars (Messenger & Read 2012) or their astrophysical mass distribution (Finn 1994;; however these methods are not expected to find an application in context of the current generation of advanced ground-based detectors. Prospects of GW cosmology using robust astrophysical features of black hole mass distributions have recently been explored in Farr et al. (2019). binary black hole (BBH) mergers with no expected EM counterparts. Even in the absence of an EM counterpart, the method outlined in Schutz (1986) can be used: with a set of potential host galaxies identified in a galaxy catalog for each detection, one can build up information by a process of statistical cross-correlation. The method was demonstrated on a set of simulations in Del Pozzo (2012), where a 5% estimate on H 0 was obtained from ( )  100 detections in an idealized situation of nearby events and complete galaxy catalogs; similar results with projections for third-generation detectors have been obtained in Nair et al. (2018). It has further been shown in Chen et al. (2018) that the main benefit of the galaxycatalog method would be for the case of multiple well-localized sources.
An understanding of GW selection effects Chen et al. 2018;Mandel et al. 2019;Vitale et al. 2020) and features of galaxy catalogs, such as their incompleteness and measurement uncertainties, is necessary for an accurate measurement of H 0 . Prescriptions to handle incomplete galaxy catalogs have been outlined in Chen et al. (2018, and Gray et al. (2020), and an extensive study of selection effects including galaxy catalog completeness has been performed in Gray et al. (2020). The simulations in Gray et al. (2020) demonstrate that the current method is able to handle galaxy catalog completeness down to ∼25% for ( )  100 detections in the nearby universe (luminosity distance d L  115 Mpc) with idealized galaxy catalogs, without introducing additional systematic features. The galaxy-catalog method was used in Fishbach et al. (2019) to infer H 0 from GW170817 without its optical counterpart. An estimate of H 0 from GW170814 and the photometric redshift catalog from the Dark Energy Survey (DES) Year 3 data was also obtained in Soares-Santos et al. (2019).
In this paper we report the first joint GW estimate of H 0 from detections during O1 and O2, the first and second observing runs of the Advanced LIGO and Virgo detector network. For our final result, along with the BNS GW170817, we choose the six O1 and O2 BBH detections, which satisfy our selection criterion of a network signal-to-noise ratio (S/N) > 12 in at least one search pipeline for GW detections. The detections for which we have significant galaxy catalog support are GW150914, GW151226, GW170608, and GW170814. For these events, we expect the inferred H 0 posteriors to be driven by the additional EM information. The two remaining BBHs that satisfy the selection criterion, GW170104 and GW170809, need to be retained in the analysis for consistency with the assumed population model. These events for which a significant fraction of potential host galaxies are missing in the associated galaxy catalog, can also potentially contribute to the H 0 measurement, since there is information about cosmology present in the observed source distribution Farr et al. 2019). The distribution of the observed source parameters, including the observed luminosity distance distribution, is driven by H 0 in addition to the intrinsic astrophysical source distributions. With a known underlying astrophysical distribution of sources, H 0 would be measurable solely from the observed distance distribution of BBH detections (Chen & Holz 2014). This latter contribution is not expected to provide significant information at this stage given how uncertain the inferred astrophysical distribution parameters such as the power-law index of the binary mass distribution, or the index characterizing the evolution of the binary merger rate with redshift are, even when H 0 is held fixed (Abbott et al. 2019a). In an ideal situation, one should jointly estimate these astrophysical population parameters along with H 0 . We choose to fix the astrophysical population to a fiducial distribution instead, and perform our analysis with different choices for the mass distribution and binary merger rate evolution with redshift in order to quantify possible systematic effects resulting from this assumption. We set aside a more thorough treatment involving the marginalization over the unknown astrophysical distribution for future work.
The main result of our analysis-a posterior distribution on H 0 -is dominated by the contribution from GW170817 with its optical counterpart, with a modest improvement from the inclusion of the O1 and O2 BBH. These results, possibly refined and marginalized over the aforementioned assumptions, can be used as a prior for future GW estimates of H 0 . The analysis performed in this paper thus serves as a precursor of future analyses for the third and subsequent observing runs of the advanced detector network.
The rest of this paper is arranged as follows. We describe our method in Section 2. We summarize the GW detections we use in our analysis and the corresponding EM data in Section 3. Our main results are presented in Section 4, with a more detailed discussion and a study of possible systematic effects in Section 5. We conclude in Section 6 and highlight some future directions and prospects.
Throughout this paper we assume a ΛCDM cosmology and use the best-fit Planck 2015 values of Ω m = 0.308 and Ω Λ = 0.692, respectively for the fractional matter and dark energy densities in the present epoch (Planck Collaboration 2016). Although these parameters enter the redshift-distance relationship central to the method for Bayesian inference of H 0 , we have verified that our results are robust with regards to a variation of their values within the current measurement uncertainties.

Method
We follow and apply the Bayesian analysis described in Gray et al. (2020) to compute the posterior probability density on H 0 , given the set {D GW } of N det detections and the associated GW data {x GW }: Here, D GWi indicates that the event i was detected as a GW, p(H 0 ) is the prior on H 0 , and the term ( | ) p N H det 0 is the likelihood of detecting N det events over the time of observation with the associated detector configuration for the given value of H 0 . The number of detected events is some fraction of the total number of sources N s , and this fraction depends on H 0 ; thus, . If the intrinsic astrophysical merger rate, R, which N s is proportional to, 197 is marginalized over with a prior is independent of ( ) f H det 0 , and thus loses its dependence on H 0 (see, e.g., Fishbach et al. 2018). For simplicity, we make this approximation throughout our analysis. The final term factorizes into the individual likelihoods for each detection. In the following, we write out the expressions for a single GW event i, omitting the subscript i for brevity of notation: The denominator, p(D GW |H 0 ), is evaluated as an integral over all possible x GW Chen et al. 2018;Mandel et al. 2019): where p(D GW |x GW , H 0 ) = 1 in the case where the S/N of x GW passes some detection threshold, and 0 in the case where it does not.

The EM Counterpart Case
In the presence of an EM counterpart, there is additional information in the EM data, x EM , which appears as an EM likelihood term. Together with this, there is an underlying assumption that there has been an EM detection; we denote this assumption as D EM . Thus, for a single event with an EM counterpart, We assume that the detectability of an EM counterpart is dependent on luminosity distance (as opposed to redshift) because it is flux limited. As GW detectability is also a function of luminosity distance, we expect p(D EM |D GW , H 0 ) to be a constant that does not depend on H 0 . This leads to

The Galaxy-catalog Case
In the absence of an EM counterpart, the analogous data comes from galaxy catalogs that provide a set of galaxies and their associated sky locations, redshifts, and apparent magnitudes. As we are in the regime where the detectability of GW sources extends beyond the distance to which current catalogs are complete, the possibility that the GW host galaxy is not contained in the catalog, because it is too faint, has to be taken into account. This is done by marginalizing over the cases where the host is in the catalog (denoted G), and where it is not (denoted¯) G : In the case where the host galaxy is in the catalog, the EM data enters in the form of information available in the galaxy catalog, as also in Fishbach et al. (2019) and Soares-Santos et al. (2019). This assumption is included in our conditioning on G. The EM information is used to modify our priors on galaxy redshift, sky location, and (apparent) magnitude, which are common among all GW events using the same catalog. This differs from the counterpart case where the EM data enters the expression as a likelihood term, x EM , a transient which is informative for only one GW event. In the case where the host galaxy is not in the catalog, the complementary conditionḠ implies that we include the information about the limitations of the EM survey. Following Gray et al. (2020), we model the galaxy catalog as having an apparent magnitude threshold, m th , since galaxy catalogs are flux limited. This magnitude threshold, alongside the intrinsic (absolute) brightness of a galaxy and its luminosity distance, determines the probability that the galaxy is inside or outside the galaxy catalog. The formalism which we adopt from Gray et al. (2020) is closely analogous to the method described in Chen et al. (2018) and Fishbach et al. (2019). The in-catalog and out-of-catalog terms in Equation (6) above correspond to the two terms in Equation (3) of Fishbach et al. (2019). However, the methods start formally differing in the way that the (in)completeness of the galaxy catalog is treated. In the current formalism from Gray et al. (2020), the (in)completeness of the galaxy catalog follows naturally from the parameters of the underlying EM survey(s). The assumption that a galaxy catalog can be modeled by a single m th in each direction of the sky is an idealization. However, this can in the future be extended to include nontrivial selection functions modeling the optical sensitivity of the telescope instead of a sharp cutoff at m th , and to composite catalogs coming from multiple surveys.
The quantities appearing on the right in Equation (6) can be written out explicitly as follows: Here, N gal is the total number of galaxies in the galaxy catalog, Ω j and m j are respectively the sky coordinates and apparent magnitude for galaxy j, and p(z j ) is a distribution representing the redshift of galaxy j. This quantity, p(z j ), which enters as a prior for our analysis, is the posterior distribution on the galaxy redshift obtained following analysis of EM data; details of the profiles chosen for p(z j ) are deferred until Section 3.4.1. The quantity M(z j , m j , H 0 ) is the absolute magnitude (for the given H 0 ), and p(s|M(z j , m j , H 0 )) is the probability of a galaxy with these parameters to host a GW source during the observation time, relative to other galaxies. Formally, s is the statement that a GW has been sourced or emitted (as opposed to being detected); the previous expressions are all implicitly conditioned on the assumption of s. In writing p(s|M), we make the approximation that the probability of a galaxy hosting a source depends only on the intrinsic luminosity of the galaxy, and not on its other parameters or on the properties of the GW source.
In essence, this term allows for weighting galaxies by their luminosities L(M j (H 0 )) as The likelihood when the host galaxy is not in the catalog, unless otherwise specified. 198 For the upper limits of integration over M, we choose the magnitude of the dimmest galaxies to be -+ h 12.2 5 log 10 . The integrals are not sensitive to the choice of their lower limits, i.e., the magnitudes of the brightest galaxies. More complex models for p(M|H 0 ) can be used; in fact, we expect the luminosity distribution of galaxies to also evolve with redshift (Caditz & Petrosian 1989), as well as to depend on galaxy type and color (Madgwick et al. 2002). While the consideration of such dependence is beyond the scope of the current work, we refer the reader to Gray et al. (2020) for a brief discussion on the misspecification of the luminosity function parameters.
Further details and complete derivations for the framework described above are discussed in Gray et al. (2020).

Posterior Samples Re-weighting
The GW data x GW has undergone parameter estimation, providing a set of posterior samples on source parameters, including the luminosity distance d L , sky coordinates Ω, and the observed component masses in the detector frame m 1,2 det . We need the likelihood of the GW data p(x GW |H 0 ), which requires removing the priors used for the parameter estimation. In particular, this means that the d L 2 prior must be removed before the priors on z and H 0 can be applied.
Additionally, a more subtle effect comes into play regarding source-and detector-frame masses. The treatment of GW selection effects requires a marginalization over the sourceframe mass distribution, while the data x GW contains information about the detector-frame masses. In order to be selfconsistent, the analysis must be performed using the same prior assumptions on both the individual event data and the normalizing p(D GW |H 0 ) term. This requires removing the detector-frame mass prior and re-weighting the posterior samples with the source-frame mass prior.
The source-frame mass distribution is p(m 1 , m 2 | μ), where μ denotes the hyper-parameters describing the astrophysical model (concrete details of the assumed model are discussed later in Section 4). To re-weight the likelihood we use the measured detector-frame masses and luminosity distances to compute the corresponding source-frame masses as a function of H 0 , We marginalize out the source-frame mass dependence using the prior p(m 1 , m 2 | μ). The GW likelihood is a function of d L and Ω as before, but now is also a function of H 0 . This effect will be most pronounced for high distance (redshift) events, where the conversion between detector-frame and source-frame masses is most noticeable. Additionally, posterior samples from events with particularly high or low detector-frame masses may become inconsistent with the source-frame mass prior for certain values of H 0 , leading to additional constraining power (see Section 5.2 for details).

GW Data
The GW searches performed during the first and the second observation runs of the Advanced LIGO and Virgo have led to the identification of 10 BBHs and one BNS merger (Abbott et al. 2019b). The BNS event GW170817, well localized and at a nearby distance of -+ 40 10 10 Mpc, helped discover the EM transient from the merger and was subsequently associated with host galaxy NGC4993. The BBHs span a large range of distances from -+ 320 110 120 to -+ 2840 1360 1400 Mpc and are distributed over the sky with 90% credible regions as low as 39 deg 2 to as high as 1666 deg 2 . A summary of the relevant parameters of all the GW detections are given in Table 1.
For the main results presented in Section 4 of this work, we choose the events which meet the selection criterion of network S/N >12 in at least one of the two pipelines for modeled searches, namely, PyCBC and GstLAL. The sensitivity of the 198 The absolute magnitude is related to the intrinsic luminosity of a galaxy by the relation, ( ) -º-M M L L 2.5 log 10 * * . The parameter M * of the Schechter function itself depends on H 0 , which we take into account. results to this choice is explored in Section 5 by reducing the network S/N threshold to 11. In each case, this S/N threshold is also used in the treatment of GW selection effects as described in Equation (3) above. In practice a detection is claimed not solely on the basis of the S/N, but additionally by applying data quality vetoes in order to remove noise transients, and eventually constructing a ranking statistic such as an inverse false alarm rate or a likelihood ratio (Abbott et al. 2019b). A careful treatment should use a threshold on a ranking statistic rather than the S/N as the selection criterion. However a distinction between the two does not cause an appreciable difference if the considered detections are significantly louder than transient noise artifacts (see, e.g., Appendix A.1 of Abbott et al. 2019a). In order to keep the computation of the GW selection effects tractable, one can thus place a more stringent threshold on the S/N and select a subset of loud events from the detected population, which is what we do for our analysis. We note that placing an S/N threshold lower than or close to that set by the detection pipelines would again be problematic without modifications to the current method of accounting for GW selection effects.

Galaxy Catalogs
The analysis with BBHs is performed in conjunction with appropriate galaxy catalogs. We use the Galaxy List for the Advanced Detector Era (GLADE) catalog (Dálya et al. 2018) as a default, due to its depth and coverage over an extensive region of the sky (see In Table 1 we summarize the galaxy catalogs that we use for our analysis for each of the detections, along with the number of galaxies in the 90% error volume calculated from 3D skymaps constructed from posterior samples associated with the data release of Abbott et al. (2019b), 199 and estimates of the probability that the host galaxy is in the catalog, evaluated at the median redshift for each detection assuming a Planck 2015 cosmology.
In the following, we describe in more detail the galaxy catalogs that we use, quantify the probability that the host galaxy for each event is in the galaxy catalog that is used for its analysis, and discuss the assessment of the completeness over the relevant localization volume for the best localized events. Finally, we quantify the uncertainties associated with the photometric measurement of redshifts in some of these catalogs.  Dálya et al. 2018). At low redshifts ( 0.05), we expect to be dominated by the peculiar-velocity field. GLADE reports peculiar-velocitycorrected redshifts following the reconstruction of Carrick et al. (2015). GLADE provides apparent magnitudes in the B-band, Note. In the remaining columns we report the corresponding 90% 3D localization comoving volumes, the number of galaxies within each volume for public catalogs which we find to be the most complete, and the apparent magnitude threshold, m th , of the galaxy catalog associated with the corresponding sky region (as described in Section 3.3). The final column gives the probability that the host galaxy is inside the galaxy catalog for each event, p(G|z event , D GW ), also evaluated at the median redshift for each event.
which we can use directly (i.e., without any photometric transformations) for luminosity weighting of the galaxies.

DES Year 1
The DES is a 5 yr survey mapping ≈300 million galaxies in five filters (grizY) over 5000 deg 2 . It is worth noting that the GW170814 sky localization is fully enclosed within the footprint of the SDSS Approximately 87% of the probability region for the GW170814 sky localization is enclosed within the DES-Y1 catalog. Analysis with a different catalog provides a parallel measurement of H 0 with GW170814, and (given the catalog differences) can potentially be indicative of systematic effects in the catalogs, such as the treatment of redshift uncertainties (provided that a similar set of galaxies are present in both catalogs).
We select the objects in the DES-Y1 catalog that are classified as high-confidence galaxies using the default classification scheme, "MODEST_CLASS" (Drlica-Wagner et al. 2018; Sevilla-Noarbe et al. 2018). We use the photometric redshifts that are derived using the Bayesian photometric redshift (BPZ) template fitting method . We use the median redshifts provided in the catalog and discard (around 5%) galaxies with redshift errors larger than twice their corresponding quoted median redshift value. Such a choice is not expected to bias our result since the discarded galaxies are highly uninformative.
We convert from the DES grizY magnitudes to the SDSS ugriz system using the photometric transformations provided in the DES-Y1 paper (Drlica-Wagner et al. 2018), which requires discarding a further ∼5% of galaxies with inadequate color information. This transformation enables us to apply K corrections to obtain source-frame luminosities (see Section 3.4.2 for details). We use the SDSS g-band magnitudes, as these are closest to

Probability That the Host Galaxy is in the Catalog
In this work, we assume that we can characterize the completeness of a galaxy catalog using an apparent magnitude threshold (limiting magnitude) m th . We estimate m th by calculating the median value from the apparent magnitude distribution of all the galaxies within the sky localization of each event. For GLADE, this choice allows us to account for some of the larger changes in completeness across the sky, which come from it being a composite catalog, comprised of many surveys of differing depths. Galaxy catalogs are directional, and a more sophisticated analysis would involve calculating the limiting magnitude for a given line of sight.
Obtaining the H 0 posterior distribution would thus require a joint estimate of m th along the lines of sights within an event's sky localization. We leave this for future work. That the completeness of a galaxy catalog is modeled by a set of limiting magnitude thresholds, can by itself be a nontrivial assumption, especially for photometric catalogs, since galaxies may be missing for various reasons other than them being too faint. This will also need to be revisited in the future in a catalog-specific manner.
For now, we use the m th estimated as described above, and show in Figure 1 the probability of a host galaxy being inside the catalog p(G|z, D GW ) as a function of redshift z, for each of the galaxy catalogs under consideration. For GLADE this quantity is calculated for each event using the m th calculated for each event's sky localization. For DES-Y1, the curve is for the patch of sky covering GW170814. These probabilities are calculated using the expressions in Equation (10), but as a function of z, and are therefore independent of the choice of H 0 . We additionally show as the vertical lines in Figure 1 the median redshift for each event z event (calculated assuming a  Table 1. These lines are thick and solid up to the intercept with the galaxy catalog they are used with, and thin and dashed above. Planck 2015 cosmology). In the final two columns of Table 1, we report the m th of the relevant catalog within the sky localization of each event, and the probability of the host galaxy being in the catalog at the median redshift of each event.

Detailed Analysis of DES-Y1
The high completeness fraction of DES-Y1 within the GW170814 sky localization is apparent from Figure 1. The catalog is expected to be more complete than GLADE since it has a limiting magnitude of approximately 23.8 for DES-Y1. We analyze the EM information coming from this catalog in greater detail. It is helpful to have an assessment of the contribution from potential host galaxies as a function of redshift for these events. In order to quantify this contribution, we perform a treatment analogous to Fishbach et al. (2019) and compute the ratio p cat (z)/p vol (z) between the probability distribution for the redshifts of potential host galaxies p cat (z) and of a uniform in comoving volume distribution of galaxies p vol (z). When computing p cat (z), we include all galaxies brighter than L 0.05 g * within the corresponding event's 99% sky localization region defined as where p(x GW |Ω) is the GW likelihood as a function of the sky position Ω (this effectively weights each galaxy with the 2D skymap probability), and p 0 (z, Ω) represents the galaxy catalog contribution, obtained from the distribution of galaxies in the catalog, marginalized over their redshift uncertainties also obtained from the catalog, and weighted by their probability of hosting a GW source (assuming a Planck 2015 cosmology for the required magnitude conversion). We consider weights for each galaxy proportional to their g-band luminosity as well as uniform weights to explore the effects due to this choice.
In Figure 2, we show the distributions p cat (z)/p vol (z) for the DES-Y1 galaxies within the GW170814 sky localization region, for the redshift range 0 < z < 0.5. The unweighted curve traces the over/under-density of galaxies, and then falls off at larger redshift due to incompleteness in the catalog. The luminosity-weighted redshift distribution is driven partially by the overdensity of galaxies at z ≈ 0.4, and partially by bright high-redshift galaxies. The host galaxies for GW170814 are more likely to be located near the higher galaxy density regions in the DES-Y1 catalog-these features in the redshift prior are expected to drive the inferred H 0 posteriors for the corresponding events. Features we see in the DES-Y1 catalog are not as pronounced as the overdensity in the DES-Y3 data seen in Soares-Santos et al. (2019). While the DES-Y3 survey is deeper, and may reveal finer features, a part of the above difference is likely also driven by the difference in the photometric redshift estimation algorithms, namely, template fitting methods such as BPZ ) and machinelearning-based methods such as the ANNz2 algorithm (Sadeh et al. 2016), with the latter used for GW170814 in (Soares-Santos et al. 2019). Only the former of the two has been used for the DES-Y1 catalog and a combination of both for the DES-Y3 catalog. The different selection criteria for choosing galaxies from the two catalogs, such as the stringent redshift cut placed in Soares-Santos et al. (2019) versus a more relaxed redshift prior used in this work, is another potential source of difference between the corresponding redshift distributions.

Redshift Uncertainties
An important source of measurement uncertainty with galaxy catalogs is in the galaxy redshifts, which are often photometric estimates due to lack of spectroscopic measurements out to large redshifts. In our formalism, the uncertainty in redshift is reflected in the term p(z j ) in Equation (7). We model the individual galaxy redshift probability distributions, p(z j ), as a Gaussian distribution with the mean set to the quoted median photometric redshift z photo and with a standard deviation of s z photo (both obtained from the galaxy catalog data). Photometric redshift estimates are often not approximated by Gaussian distributions, and we make this assumption only due to limited information present in some of the public galaxy catalogs that we use. A rigorous quantification of the effect of this assumption is beyond the scope of this work. In practice, we implement the redshift uncertainty by the process of a Monte Carlo integration: the integral over z j in Equation (7) becomes an additional sum over N photo random samples. This number varies between the different galaxies and is always high enough to ensure that the final event likelihood is independent of the set of random draws from each galaxy.

Source-frame Luminosities
For galaxy catalogs that are complete to high redshifts, it is important to note that the observed apparent magnitudes are redshifted from their source frame. Galaxies do not emit uniformly in all wavelengths, and galaxy surveys are only sensitive in certain bands. In order to compare the properties of galaxies at different redshifts, e.g., to apply luminosity weighting, or to apply a luminosity cut to the galaxy sample, the K correction, described in Oke & Sandage (1968) and Hogg et al. (2002), needs to be applied. In particular, the K correction term is redshift dependent and so neglecting it could lead to systematic bias in the estimation of H 0 .
For DES, we compute the K corrections using the calculator described in Chilingarian et al. (2010) and Chilingarian & Zolotukhin (2012). These K corrections are valid out to redshifts Figure 2. Probability distribution for the redshifts of potential host galaxies p cat (z), with redshift uncertainties taken into account, divided out by a uniform in comoving volume distribution p vol (z) of galaxies. When computing p cat (z) we include all galaxies brighter than L 0.05 g * within the corresponding event's 99% sky localization region and weight each galaxy by weights proportional to their g-band luminosity (solid lines) as well as with uniform weights (dotted lines). We show these distributions for the DES-Y1 galaxies within the GW170814 sky localization region. We also show the 90% median estimated redshift range for GW170814 (calculated assuming a Planck 2015 cosmology) for reference. et al. (2019) with analogous peaks in the posterior distribution. The differences in peak locations can be attributed to a difference in the redshift distribution for the DES-Y3 catalog used in Soares-Santos et al. (2019) versus that for the public DES-Y1 catalog used in this work. For the detections for which the galaxy catalogs are relatively empty, we see the features of the assumptions on mass distribution and redshift evolution of binary merger rate that have entered our analysis. The more distant events, such as GW170809, lead to H 0 estimates pushed to the lower end of the prior. This is due to the following reason. In a universe where host galaxies are distributed uniformly in comoving volume, for lower values of H 0 , the expected distribution of detected GW events favors relatively higher luminosity distances. Thus, events at high luminosity distances have more support for smaller values of H 0 , while relatively nearby events, namely, GW150914, GW151226, GW170608, and GW170814, correspondingly, have lower support at smaller values of H 0 . The information present in the observed luminosity distance distribution would thereby potentially contribute to the H 0 measurement, independent of, or even in absence of information in galaxy catalogs, if the underlying distributions of source parameters were known.
In the absence of knowledge of the astrophysical distribution of BBH source parameters, a thorough treatment would involve a marginalization over all possible mass distributions and rate models. The following section discusses the systematic differences that the assumptions on the assumed population model could lead to.

Systematic Effects
In this section, we repeat the previous analysis with a few alternative assumptions for the GW selection criterion, the GW population model, and source host properties. Subsequently, we discuss other sources of systematic effects, a detailed study of which is beyond the scope of this work.

Selection Criterion
We first look into the sensitivity of our results to the GW selection criterion. We reduce the threshold on the network S/N from 12 chosen in the previous section to 11 in least one search pipeline. The computation of the GW selection effects is adjusted accordingly. Figure 5 shows the results with the two sets of assumptions. Reducing the S/N threshold to 11 introduces two additional events in our analysis, namely, GW170818 and GW170823, neither of which have a significant in-catalog contribution. Differences are expected due to the fact that additional low-S/N events are introduced, and also because the individual likelihoods change slightly with a different S/N threshold used in the GW selection term. In the regime of a large number of events, these two effects are expected to cancel, provided that the additional low-S/N events do not have significant in-catalog support. For our result, this difference is not significant; however, the small variation is a reminder that we are still in the regime of a low number of events.

Population Model
Going back to our default assumption of a network S/N threshold of 12, we test the sensitivity to our assumptions regarding the population model, i.e., the mass distribution and the distribution of binary merger rate with redshift. In addition to the power-law mass distribution with α = 1.6 (median inferred value using Model B of Abbott et al. 2019a), we choose a shallower flat-in-log mass distribution with α = 1, and a steeper distribution with α = 2.3 (both within the support of the inferred range in Abbott et al. 2019a). Our results are shown in the left panel of Figure 6, and they demonstrate that the systematic differences due to the choice of power-law slope α are insignificant.
As a test case, we vary the upper cutoff for the mass distribution, M max . For our default analysis, M max was chosen to be 100 M e , consistent with all the considered BBHs for all values of H 0 within the prior range. Reducing this cutoff to a slightly restrictive  = M M 50 max (e.g., motivated by the pair instability supernova process, Fowler & Hoyle 1964), we see a significant difference (right panel of Figure 6). Lowering M max corresponds to a closer GW detection horizon. This systematically leads each event to prefer slightly lower values of H 0 than in the main result, for the reasons outlined in Section 4, namely, the relationship between the predicted event distribution (from our GW selection effects) and the real detected event distribution.
Next, we relax our assumption on the evolution of rate of binary mergers with redshift. A constant merger rate density, ( ) = R z constant, implicit in the previous treatment, assumed that the merger rate traces the comoving volume. In addition, we repeat our analysis using a merger rate R(z) ∝ (1 + z) 3 that traces the star formation rate at low redshifts (z < 2.5) (Saunders et al. 1990), as well as a merger rate R(z) ∝ (1 + z) −3 that could arise if typical delay times are very long, as may be expected from the chemically homogeneous evolution formation channel (Mandel & de Mink 2016). These relaxed assumptions thus cover a large fraction of physically viable and inferred population models (Abbott et al. 2019a). We show our results for the different assumed redshift evolution models in Figure 7. The model in which the merger rate traces star formation shows a significant difference, with a tendency to prefer lower values of H 0 , compared to the other two models which are quite similar. This is the only systematic effect that leads to a significant difference in the results at this time.

Luminosity Weighting
The results in the previous section assumed a weighting of galaxies by their luminosities in the B-band, which are indicative of galaxies' star formation rates. Luminosities in the infrared (such as the K-band) are more indicative of the total masses of galaxies; however, infrared luminosities are not present in catalogs like DES-Y1. In order to quantify the difference likely to be caused by alternate ways of weighting the galaxies, we repeat our analysis with no luminosity weighting. These results are shown in Figure 8.
With uniform luminosity weights, we obtain a result on a joint BBH estimate that is close to flat (thin orange line in Figure 8). This can be understood as follows: (1) The out-ofcatalog terms in Equation (6) take into account the lack of galaxies beyond the apparent magnitude threshold m th of the catalog in a uniform way. When galaxies are unweighted, the probability of the host galaxy being inside the catalog is reduced compared to the luminosity-weighted case. More weight is given to the uniform out-of-catalog contribution, and the events become less informative.
(2) The catalogs used in this analysis contain high numbers of low-luminosity galaxies. The contribution from these more evenly distributed dim galaxies, in the unweighted case, again reduces the informativeness of each event and flattens the final result. This is also in agreement with our expectations from Fishbach et al. (2019) and Gray et al. (2020), where weighting by luminosities enhance the features in the posterior distribution coming from the galaxy catalog.

Waveform Models
The posterior samples of Abbott et al. (2019b) used for the results in this paper have been obtained combining the results of gravitational waveform models that incorporate spin and precession effects to different extents (Hannam et al. 2014;Pan et al. 2014;Taracchini et al. 2014;Husa et al. 2016;Khan et al. 2016;Babak et al. 2017). These models are restricted to quasicircular orbits (i.e., they do not include orbital eccentricity) and Figure 5. Sensitivity of the results to the selection criterion on the S/N. The final joint result as well as the contribution from all the BBHs that satisfy the selection criterion are shown for a threshold network S/N of 12 (thick orange) and 11 (thin orange); the variation is shown as a shaded band. The GW170817 counterpart result (gray) is added to guide the eye. Six BBHs (GW150914, GW151226, GW170104, GW170608, GW170809, and GW170814) pass the selection criterion with an S/N >12. Two additional BBHs (GW170818 and GW170823) are included with an S/N >11. These results assume the default m −1.6 power-law distribution on masses and a non-evolving rate model. neglect higher-order harmonics. Systematic differences in GW parameter estimation results with the employed waveform models constitute only a small fraction of the total uncertainty budget (see, e.g., Abbott et al. 2016aAbbott et al. , 2017b, and given the large statistical uncertainties, the ignored effects in waveform modeling are not expected to cause a difference in the current measurement of H 0 . However cumulative systematic effects arising from limitations of waveform models will become increasingly important as the statistical uncertainties become smaller, and in particular, features that can lead to biases in the GW estimation of distance will need to be incorporated.

Detector Calibration
An independent effect to be considered is the calibration of the GW detectors. Currently, the GW parameter estimation results are marginalized over the detector calibration uncertainties ( 4% in amplitude in O1 and O2), which accounts for both the statistical uncertainty and the systematic error correlated between detections (Abbott et al. 2019b). Both the statistical uncertainty and the systematic error in GW detector calibration are much smaller than the other measurement uncertainties, and thus negligible for H 0 estimates from a handful of detections that we have now or expect in the near future (Cahillane et al. 2017). However, the impact of correlated systematic calibration errors between detections will become relatively more important in the long term, with an increasing number of detections driving down the statistical uncertainties, and an improved understanding of other systematic effects that possibly govern our current uncertainty budget. Further quantitative study of the effect of correlated calibration uncertainties is ongoing.

Conclusion and Outlook
In this paper we have presented the first measurement of H 0 using multiple GW observations. Our result reanalyzes and combines the posterior probability distribution obtained from the BNS event GW170817 using the redshift of the host galaxy inferred from the observed EM counterpart , along with constraints using galaxy catalogs for the BBH events observed by the Advanced LIGO and Virgo in   Figure 8. Sensitivity of the results to luminosity weighting. We show how the results vary when we weight the galaxies in the catalog by their B(g)-band luminosity (thick solid) as well as with constant (uniform) weights (thin solid), both assuming a power-law index for the mass distribution, α = 1.6, and constant intrinsic astrophysical merger rate, ( ) = R z constant.