A gravitational-wave measurement of the Hubble constant following the second observing run of Advanced LIGO and Virgo

TL;DR

This paper reports a gravitational-wave based measurement of the Hubble constant using multiple detections from LIGO and Virgo, improving the precision of previous estimates and demonstrating the potential of gravitational-wave cosmology.

Contribution

It provides a joint measurement of the Hubble constant combining binary neutron star and black hole detections, advancing gravitational-wave cosmology methods.

Findings

Measured H_0 = 68.7^{+17.0}_{-7.8} km/s/Mpc

Improved measurement precision by about 4% over previous GW170817-only estimate

Highlights importance of future detections and systematic understanding

Abstract

This paper presents the gravitational-wave measurement of the Hubble constant ($H_0$) using the detections from the first 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 first 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 = 68.7^{+17.0}_{-7.8}$ km/s/Mpc (68.3\% of the highest density posterior interval with a flat-in-log prior) which is an improvement by a factor of 1.04 (about 4\%) over the GW170817-only value of $68.7^{+17.5}_{-8.3}$ km/s/Mpc. A significant 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.

Tables

Table 1. Table 1: Relevant parameters of the O1 and O2 detections: network signal-to-noise ratio (SNR) for the search pipeline (PyCBC/GstLAL) in which the signal is the loudest, 90 % percent 90 90\% sky localization region Δ ​ Ω Δ Ω \Delta\Omega ( deg 2 superscript deg 2 \text{deg}^{2} ), luminosity distance d L subscript 𝑑 L d_{\mathrm{L}} ( Mpc , median with 90 % percent 90 90\% credible intervals), and estimated redshift z event subscript 𝑧 event z_{\text{event}} (median with 90 % percent 90 90\% range assuming Planck 2015 cosmology) from Abbott et al. ( 2019b ) . In the remaining columns we report the corresponding 90 % percent 90 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 subscript 𝑚 th m_{\text{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 ) 𝑝 conditional 𝐺 subscript 𝑧 event subscript 𝐷 GW p(G|z_{\text{event}},D_{\text{GW}}) , also evaluated at the median redshift for each event.

Event	SNR	$\Delta \Omega /$${\text{deg}}^2$	$d_L/$Mpc	$z_{\text{event}}$	$V/{\text{Mpc}}^3$	Galaxy catalog	Number of galaxies	$m_{\text{th}}$	$p\left(G|z_{\text{event}},D_{\text{GW}}\right)$
GW150914	24.4	182	${440}_{-170}^{+150}$	${0.09}_{-0.03}^{+0.03}$	$3.5\times {10}^6$	GLADE	3910	17.92	$0.42$
GW151012	10.0	1523	${1080}_{-490}^{+550}$	${0.21}_{-0.09}^{+0.09}$	$5.8\times {10}^8$	GLADE	78195	17.97	$0.01$
GW151226	13.1	1033	${450}_{-190}^{+180}$	${0.09}_{-0.04}^{+0.04}$	$2.4\times {10}^7$	GLADE	27677	17.93	$0.41$
GW170104	13.0	921	${990}_{-430}^{+440}$	${0.20}_{-0.08}^{+0.08}$	$2.4\times {10}^8$	GLADE	42221	17.76	$0.01$
GW170608	15.4	392	${320}_{-110}^{+120}$	${0.07}_{-0.02}^{+0.02}$	$3.4\times {10}^6$	GLADE	6267	17.84	$0.60$
GW170729	10.8	1041	${2840}_{-1360}^{+1400}$	${0.49}_{-0.21}^{+0.19}$	$8.7\times {10}^9$	GLADE	77727	17.82
GW170809	12.4	308	${1030}_{-390}^{+320}$	${0.20}_{-0.07}^{+0.05}$	$9.1\times {10}^7$	GLADE	18749	17.62
GW170814	16.3	87	${600}_{-220}^{+150}$	${0.12}_{-0.04}^{+0.03}$	$4.0\times {10}^6$	DES-Y1	31554	23.84	$>0.99$
GW170817	33.0	16	${40}_{-15}^{+7}$	${0.01}_{-0.00}^{+0.00}$	$227$	–	–	–	–
GW170818	11.3	39	${1060}_{-380}^{+420}$	${0.21}_{-0.07}^{+0.07}$	$1.5\times {10}^7$	GLADE	1059	17.51
GW170823	11.5	1666	${1940}_{-900}^{+970}$	${0.35}_{-0.15}^{+0.15}$	$3.5\times {10}^9$	GLADE	117680	17.98

Equations

$$p(H_0|\{x_{\text{GW}}\},\{D_{\text{GW}}\})\propto p(H_0)p(N_{\text{det}}|H_0)\prod _i^{N_{\text{det}}}p({x_{\text{GW}}}_i|{D_{\text{GW}}}_i,H_0).$$

$$p(x_{\text{GW}}|D_{\text{GW}},H_0)$$

$$p(D_{\text{GW}}|H_0)=\int p(D_{\text{GW}}|x_{\text{GW}},H_0)p(x_{\text{GW}}|H_0)d{x}_{\text{GW}},$$

$$p(x_{\text{GW}},x_{\text{EM}}|D_{\text{GW}},D_{\text{EM}},H_0)$$

$$p(x_{\text{GW}}|D_{\text{GW}},H_0)$$

$$+p(x_{\text{GW}}|\bar{G},D_{\text{GW}},H_0)p(\bar{G}|D_{\text{GW}},H_0).$$

$$p($$

$$=\frac{{\sum }_{j=1}^{N_{\text{gal}}}\int p(x_{\text{GW}}|z_j,{\Omega }_j,H_0)p(s|M(z_j,m_j,H_0))p(z_j)d{z}_j}{{\sum }_{j=1}^{N_{\text{gal}}}\int p(D_{\text{GW}}|z_j,{\Omega }_j,H_0)p(s|M(z_j,m_j,H_0))p(z_j)d{z}_j}.$$

$$p(s|M(z_j,m_j,H_0))$$

$$p(x_{\text{GW}}|\bar{G},D_{\text{GW}},H_0)=\frac{{\iiint }_{z(m_{\text{th}},M,H_0)}^{\infty }p(x_{\text{GW}}|z,\Omega ,H_0)p(z)p(\Omega )p(M|H_0)p(s|M)dzd\Omega dM}{{\iiint }_{z(m_{\text{th}},M,H_0)}^{\infty }p(D_{\text{GW}}|z,\Omega ,H_0)p(z)p(\Omega )p(M|H_0)p(s|M)dzd\Omega dM}.$$

$$p(G|D_{\text{GW}},H_0)$$

$$p(z)\propto \frac{1}{1+z}\frac{d{V}_{\text{c}}(z)}{dz}R(z).$$

$$p(M|H_0)\propto 10^{-0.4(\alpha +1)(M-M^{\ast }(H_0))}\exp \left[-10^{-0.4(M-M^{\ast }(H_0))}\right].$$

$$m_{1,2}=\frac{m_{1,2}^{\text{det}}}{1+z(d_{\mathrm{L}},H_0)}.$$

$$p_{\text{cat}}(z)\equiv \int p(x_{\text{GW}}|\Omega )p_0(z,\Omega )d\Omega ,$$

Full text

1 Introduction

Gravitational waves from compact binary coalescences allow for a 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}$ (Abbott et al., 2017a). 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 Type Ia supernovae. The importance of an independent measurement of $H_{0}$ is worth highlighting. With the Planck 2018 data release (Aghanim 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 SH0ES project has risen to the $4.4$-$\sigma$ level. Independent measurements using cosmological Baryon Acoustic Oscillations to calibrate Type Ia supernovae 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 multi-messenger GW astronomy sector would help clarify whether the current tension is a statistical anomaly or evidence for new physics beyond the $\Lambda$CDM model of cosmology.111Cosmological 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; Taylor & Gair, 2012); 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).

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 one percent with $\mathcal{O}(100)$ detections in the coming years (Dalal et al., 2006; Nissanke et al., 2013; Chen et al., 2018; Feeney et al., 2019; Vitale & Chen, 2018; 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 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 $\mathcal{O}(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 galaxy-catalog method would be for the case of multiple well-localized sources.

An understanding of GW selection effects (Abbott et al., 2017a; 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), Fishbach et al. (2019), 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 $\sim 25\%$ for $\mathcal{O}(100)$ detections in the nearby universe (luminosity distance $d_{\mathrm{L}}\lesssim 115\,\mbox{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 $\text{SNR}>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 which 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 (Taylor et al., 2012; Taylor & Gair, 2012; 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 measureable 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 BBHs. 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 $\Lambda$CDM cosmology and use the best-fit Planck 2015 values of $\Omega_{m}=0.308$, $\Omega_{\Lambda}=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 Hubble constant ($H_{0}$), we have verified that our results are robust with regards to a variation of their values within the current measurement uncertainties.

2 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_{\text{GW}}\}$ of $N_{\text{det}}$ detections and the associated GW data $\{x_{\text{GW}}\}$:

[TABLE]

Here, $D_{\text{GW}i}$ indicates that the event $i$ was detected as a GW, $p(H_{0})$ is the prior on $H_{0}$, and the term $p(N_{\text{det}}|H_{0})$ is the likelihood of detecting $N_{\text{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 $N_{\text{det}}=f_{\text{det}}(H_{0})\,N_{s}$. If the intrinsic astrophysical merger rate, $R$, which $N_{s}$ is proportional to,222The astrophysical merger rate is defined as $R\equiv{\partial{N_{s}}}/{(\partial V\,\partial T)}$, where $N_{s}$ is the number of sources as in the text, $V$ the comoving sensitive volume, and $T$ the time of observation or survey. is marginalized over with a prior $p(R)\propto R^{-1}$, then $p(N_{\text{det}}|H_{0})=\int p(N_{\text{det}}|H_{0},R)\,p(R)\,dR$ is independent of $f_{\text{det}}(H_{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 factorises 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,

[TABLE]

The denominator, $p(D_{\text{GW}}|H_{0})$, is evaluated as an integral over all possible $x_{\text{GW}}$ (Abbott et al., 2017a; Chen et al., 2018; Mandel et al., 2019):

[TABLE]

where $p(D_{\text{GW}}|x_{\text{GW}},H_{0})=1$ in the case where the signal-to-noise ratio (SNR) of $x_{\text{GW}}$ passes some detection threshold, and 0 in the case where it does not.

2.1 The electromagnetic counterpart case

In the presence of an EM counterpart, there is additional information in the EM data, $x_{\text{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_{\text{EM}}$. Thus, for a single event with an EM counterpart,

[TABLE]

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_{\text{EM}}|D_{\text{GW}},H_{0})$ to be a constant that does not depend on $H_{0}$. This leads to

[TABLE]

2.2 The galaxy-catalog case

In the absence of an EM counterpart, the analogous data comes from galaxy catalogs which 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 $\bar{G}$):

[TABLE]

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_{\text{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 $\bar{G}$ 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_{\text{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 Eq. (6) above correspond to the two terms in Eq. (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 modelled by a single $m_{\text{th}}$ in each direction of the sky is an idealization. However this can in future be extended to include non-trivial selection functions modeling the optical sensitivity of the telescope instead of a sharp cutoff at $m_{\text{th}}$, and to composite catalogs coming from multiple surveys.

The quantities appearing on the right in Eq. (6) can be written out explicitly as follows.

[TABLE]

Here, $N_{\text{gal}}$ is the total number of galaxies in the galaxy catalog, $\Omega_{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

[TABLE]

The likelihood when the host galaxy is not in the catalog, $p(x_{\text{GW}}|\bar{G},D_{\text{GW}},H_{0})$, is a ratio of marginalized integrals: {widetext}

[TABLE]

Here the fact that the terms are conditioned on $\bar{G}$ is incorporated into the redshift limits as a function of the apparent magnitude threshold $m_{\text{th}}$ of the galaxy catalog. Finally, the prior probabilities that a given GW detection has or does not have support in the galaxy catalog are respectively

[TABLE]

In Eqs. (9) and (10), $p(z)$ is the prior on the redshift of host galaxies of GW events, taken to be of the form

[TABLE]

Here $V_{\text{c}}(z)$ is the comoving volume as a function of redshift and the factor $(1+z)^{-1}$ converts the merger rate from the source frame to the detector frame. The merger rate density may in general be a function of redshift; however we set $R(z)=\text{constant}$ throughout (other than in Section 5, where we consider an alternative redshift-dependent rate model). The prior on the GW sky location $p(\Omega)$ is taken to be uniform across the sky. The term $p(M|H_{0})$ is the prior on absolute magnitudes for all the galaxies in the universe (not just those inside the galaxy catalog), which we set to follow the Schechter luminosity function:

[TABLE]

Following Gehrels et al. (2016), we use $B$-band luminosity function parameters $\alpha=-1.07$ for the slope of the Schechter function and $M^{*}(H_{0})=-19.7+5\log_{10}h$ for its characteristic absolute magnitude (with $h\equiv H_{0}/100\,\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$) throughout, unless otherwise specified.333The absolute magnitude is related to the intrinsic luminosity of a galaxy by the relation, $M-M^{*}\equiv-2.5\log_{10}(L/L^{*})$. The parameter $M^{*}$ of the Schechter function itself depends on $H_{0}$, which we take into account. For the upper limits of integration over $M$, we choose the magnitude of the dimmest galaxies to be $-12.2+5\log_{10}h$. 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).

2.3 Posterior samples re-weighting

The GW data $x_{\text{GW}}$ has undergone parameter estimation, providing a set of posterior samples on source parameters, including the luminosity distance $d_{\mathrm{L}}$, sky-coordinates $\Omega$, and the observed component masses in the detector frame $m_{1,2}^{\text{det}}$. We need the likelihood of the GW data $p(x_{\text{GW}}|H_{0})$, which requires removing the priors used for the parameter estimation. In particular, this means that the $d_{\mathrm{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 marginalisation over the source-frame mass distribution, while the data $x_{\text{GW}}$ contains information about the detector-frame masses. In order to be self-consistent, the analysis must be performed using the same prior assumptions on both the individual event data, and the normalizing $p(D_{\text{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}|\,\mu)$ where $\mu$ 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}$,

[TABLE]

We marginalize out the source-frame mass dependence, using the prior $p(m_{1},m_{2}|\,\mu)$. The GW likelihood is a function of $d_{\mathrm{L}}$ and $\Omega$ 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).

3 Data

3.1 Gravitational-wave data

The GW searches performed during the first and the second observation runs of Advanced LIGO and Virgo have led to the identification of ten BBH and one BNS mergers (Abbott et al., 2019b). The BNS event GW170817, well-localized and at a nearby distance of $40^{+10}_{-10}$ Mpc, helped discover the electromagnetic transient from the merger, and was subsequently associated with host galaxy NGC4993. The BBHs span a large range of distances from $320^{+120}_{-110}$ Mpc to $2840^{+1400}_{-1360}$ Mpc and are distributed over the sky with $90\%$ credible regions as low as $39$ $\text{deg}^{2}$ to as high as $1666$ $\text{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 SNR $>12$ in at least one of the two pipelines for modeled searches, namely PyCBC and GstLAL. The sensitivity of the results to this choice is explored in Section 5 by reducing the network SNR threshold to $11$. In each case, this SNR threshold is also used in the treatment of GW selection effects as described in Eq. (3) above. In practice a detection is claimed not solely on the basis of the SNR, 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 SNR 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 SNR and select a subset of loud events from the detected population, which is what we do for our analysis. We note that placing an SNR threshold lower than or close to that set by the detection piplines would again be problematic without modifications to the current method of accounting for GW selection effects.

3.2 Galaxy Catalogs

The analysis with BBHs is performed in conjunction with appropriate galaxy catalogs. We use the GLADE catalog (Dálya et al., 2018) as a default, due to its depth and coverage over an extensive region of the sky (see Section 3.2.1). For GW observations that are particularly well-localized, certain galaxy catalogs show a clear improvement in completeness over GLADE within the relevant localization volume of the event. In particular, we use the DES Year 1 (Y1A1 GOLD or simply Y1) catalog (Drlica-Wagner et al., 2018; Abbott et al., 2018b) for the analysis of GW170814 (see Section 3.2.2). GW170818 lies within the footprint of the Sloan Digital Sky Survey (SDSS). While not used in this work, galaxy catalogs based on Data Release 14 (Abolfathi et al., 2018) or Data Release 16 (Ahumada et al., 2020) of SDSS (including curated versions such as the GWENS catalog, Rahman et al., 2019), could be used to improve the current analysis for events which fall in the SDSS footprint.

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),444Available at: https://www.gw-openscience.org/GWTC-1 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.

3.2.1 GLADE

We use the Galaxy List for the Advanced Detector Era (GLADE) version 2.4 galaxy catalog (Dálya et al., 2018),555GLADE is publicly available at: http://glade.elte.hu to construct the observed redshift distributions for the majority of the detected BBHs. The GLADE catalog has an all sky coverage (Fig. 1 of Dálya et al., 2018) since it is constructed from the GWGC (White et al., 2011), 2MPZ (Bilicki et al., 2014), 2MASS XSC (Skrutskie et al., 2006), HyperLEDA (Makarov et al., 2014) and SDSS-DR12Q (Pâris et al., 2017) catalogs. The GLADE catalog is complete (in $B$-band luminosity) out to 37 Mpc and has an estimated completeness of 50% out to 91 Mpc (Fig. 2 of Dálya et al., 2018). At low redshifts ($\lesssim 0.05$), we expect to be dominated by the peculiar velocity field. GLADE reports peculiar-velocity-corrected redshifts following the reconstruction of Carrick et al. (2015). GLADE provides apparent magnitudes in the $B$-band, which we can use directly (i.e. without any photometric transformations) for luminosity weighting of the galaxies.

3.2.2 DES Year 1

The Dark Energy Survey (DES) is a five year survey mapping $\approx$ 300 million galaxies in five filters (grizY) over $5000~{}\text{deg}^{2}$. It is worth noting that the GW170814 sky localization is fully enclosed within the footprint of the DES (Drlica-Wagner et al., 2018; Abbott et al., 2018b) Year 3 (Y3) “gold” catalog. An estimate of $H_{0}$ from the GW170814 distance and the Y3 catalog of the DES has been carried out (Soares-Santos et al., 2019). In this work, we use the publicly available DES-Y1 catalog (Abbott et al., 2018b),666DES-Y1 is available at: https://des.ncsa.illinois.edu/releases/y1a1 to compute the $H_{0}$ posterior for GW170814. 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 which are derived using the Bayesian Photometric Redshift (BPZ) template fitting method (Hoyle et al., 2018). 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 $\sim 5$% of galaxies with inadequate color information. This transformation enables us to apply K corrections to obtain source-frame luminosities (see Sec. 3.4.2 for details). We use the SDSS $g$-band magnitudes, as these are closest to $B$-band, and update the Schechter parameters of our analysis to have $\alpha=-0.89$ and $M^{*}(H_{0})=-19.39+5\log_{10}h$ based on Blanton et al. (2003).

3.3 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_{\text{th}}$. We estimate $m_{\text{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_{\text{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 modelled by a set of limiting magnitude thresholds, can by itself be a non-trivial 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_{\text{th}}$ estimated as described above, and show in Fig.1 the probability of a host galaxy being inside the catalog $p(G|z,D_{\text{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_{\text{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 Eq. (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 Fig. 1 the median redshift for each event $z_{\text{event}}$ (calculated assuming a Planck 2015 cosmology). In the final two columns of Table 1, we report the $m_{\text{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.

3.4 Detailed analysis of DES-Y1

The high completeness fraction of DES-Y1 within the GW170814 sky localization is apparent from Fig. 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_{\text{cat}}(z)/p_{\text{vol}}(z)$ between the probability distribution for the redshifts of potential host galaxies $p_{\text{cat}}(z)$ and of a uniform in comoving volume distribution of galaxies $p_{\text{vol}}(z)$. When computing $p_{\text{cat}}(z)$ we include all galaxies brighter than $0.05L_{g}^{*}$ within the corresponding event’s 99% sky localization region defined as,

[TABLE]

where $p(x_{\text{GW}}|\Omega)$ is the GW likelihood as a function of the sky position $\Omega$ (this effectively weights each galaxy with the 2D skymap probability), and $p_{0}(z,\Omega)$ 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 Fig. 2 we show the distributions $p_{\text{cat}}(z)/p_{\text{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 over-density of galaxies at $z\approx 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 (Hoyle et al., 2018) and machine learning 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.

3.4.1 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 Eq. (7). We model the individual galaxy redshift probability distributions, $p(z_{j})$, as a Gaussian distribution with mean set to the quoted median photometric redshift $z_{\text{photo}}$, and with standard deviation of $\sigma_{z_{\text{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 which 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 Eq. (7) becomes an additional sum over $N_{\text{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.

3.4.2 Source-frame luminosities

For galaxy catalogs which 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); 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); Chilingarian & Zolotukhin (2012). These $K$ corrections are valid out to redshifts of 0.5, and so we apply a hard redshift cut. In order to remain self-consistent, the method is adjusted to include this redshift cut, such that the probability of the host galaxy being in the catalog above $z=0.5$ is zero, and below the cut it is determined according to the apparent magnitude threshold, as before.

In the case where galaxies have high redshifts and large redshift uncertainties, this uncertainty must be accounted for in the calculation of the $K$ corrections. Using only the mean redshift to calculate the $K$ correction would lead to unrealistically bright or faint luminosities at the tail ends of the redshift distribution. In our implementation, the $K$ corrections are calculated in-situ for each redshift Monte Carlo random draw described in 3.4.1, automatically taking the galaxy redshift uncertainties into account.

4 Results

We apply the method described above to obtain a measurement of the Hubble constant using GW standard sirens only. We carry out our analysis with a prior on $H_{0}$ uniform in the interval of $[20,\,140]~{}\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$; we report our final results also using a flat-in-log prior $p(H_{0})\propto H_{0}^{-1}$ in the same interval for ease of comparison with previous studies. We use the marginalized distance likelihood and skymaps constructed from the posterior samples of Abbott et al. (2019b).777For computational convenience, we separately construct a marginalized distance likelihood and a two-dimensional skymap; this approximation will be revisited in the future. For the BBHs, we choose all galaxies in the $99.9\%$ sky region of the corresponding catalog and we further weight the galaxies in proportion to their $B(g)$-band luminosities. In order to calculate the term $p(D_{\text{GW}}|H_{0})$ in the denominator, we use a Monte Carlo integration, sampling parameters which affect an event’s detectability (masses, sky location, inclination angle, and polarisation) from chosen priors. We choose a power-law mass distribution for BBHs with $p(m_{1})\propto m_{1}^{-\alpha}$ and $m_{2}$ uniform in its range with $5M_{\odot}<m_{2}<m_{1}<100M_{\odot}$ in the source frame, and a distribution of merger rates that does not evolve with redshift; for the power-law index $\alpha$, we choose $\alpha=1.6$ (which is supported by Model B of Abbott et al., 2019a). For BNSs, we use a Gaussian mass distribution with a mean of $1.35M_{\odot}$ and a standard deviation of $0.15M_{\odot}$ (Kiziltan et al., 2013).888While this does not significantly affect our current results, we will need to revisit the BNS mass distribution in light of GW190425 (Abbott et al., 2020), which is potentially a BNS merger with heavier components. The remaining GW parameters are marginalized over their natural distributions: uniform in the sky, uniform on the sphere for orientation, uniform in polarization. We use the time-averaged power-spectral-density of detector noise for the corresponding observation run from Abbott et al. (2018a), and for the selection criterion, we use a network SNR threshold of $12$ in at least one search pipeline. The O1 and O2 BBHs which pass this criterion are GW150914, GW151226, GW170104, GW170608, GW170809, and GW170814 (see Table 1).

Our result for these O1 and O2 BBHs is shown in Fig. 3. The detections for which there is considerable support from the galaxies present in the catalog show features of the galaxy catalog in their $H_{0}$ posterior distribution. The GW170814 estimate is qualitatively similar to the result in Soares-Santos 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.

For our final result we combine the contribution of the BBHs with the result from GW170817 obtained using the low spin prior samples from Abbott et al. (2019b) and an estimated Hubble velocity of $v_{H}\equiv cz=3017\pm 166\,\text{km}\text{s}^{-1}$ (where $c$ is the speed of light) for NGC4993 from Abbott et al. (2017a). Our final combined result is shown in Fig. 4, with the posterior distribution plotted assuming a uniform $H_{0}$ prior: we obtain $H_{0}=$ $69.6^{+20.4}_{-8.6}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ ($68.3\%$ highest density posterior interval). To compare with values in the literature, we also use a flat-in-log prior, $p(H_{0})\propto H_{0}^{-1}$, and calculate $H_{0}=$ $68.7^{+17.0}_{-7.8}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$, which corresponds to an improvement by a factor of 1.04 (about 4%) over the GW170817-only value of $68.7^{+17.5}_{-8.3}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$. We also quote the median and symmetric 90% credible interval for this measurement, which is $74.9^{+39.2}_{-13.9}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$.

5 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.

5.1 Selection criterion

We first look into the sensitivity of our results to the GW selection criterion. We reduce the threshold on the network SNR from $12$ chosen in the previous section to $11$ in least one search pipeline. The computation of the GW selection effects is adjusted accordingly. Fig. 5 shows the results with the two sets of assumptions. Reducing the SNR 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-SNR events are introduced, and also because the individual likelihoods change slightly with a different SNR 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-SNR 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.

5.2 Population model

Going back to our default assumption of a network SNR 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 $\alpha=1.6$ (median inferred value using Model B of Abbott et al., 2019a), we choose a shallower flat-in-log mass distribution with $\alpha=1$, and a steeper distribution with $\alpha=2.3$ (both within the support of the inferred range in Abbott et al., 2019a). Our results are shown in the left panel of Fig. 6, and they demonstrate that the systematic differences due to the choice of power-law slope $\alpha$ are insignificant.

As a test case, we vary the upper cut-off for the mass distribution, $M_{\text{max}}$. For our default analysis, $M_{\text{max}}$ was chosen to be $100\,M_{\odot}$, consistent with all the considered BBHs for all values of $H_{0}$ within the prior range. Reducing this cut-off to a slightly restrictive $M_{\text{max}}=50\,M_{\odot}$ (e.g., motivated by the pair instability supernova process, Fowler & Hoyle, 1964), we see a significant difference (right panel of Fig. 6). Lowering $M_{\text{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)=\text{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)\propto(1+z)^{3}$, which traces the star formation rate at low redshifts ($z<2.5$) (Saunders et al., 1990), as well as a merger rate $R(z)\propto(1+z)^{-3}$ which 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 Fig. 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.

5.3 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. 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 Fig. 8.

With uniform luminosity weights, we obtain a result on a joint binary black hole estimate which is close to flat (thin orange line in Fig. 8). This can be understood as follows: 1) The out-of-catalog terms in Eq. (6) take into account the lack of galaxies beyond the apparent magnitude threshold $m_{\text{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.

5.4 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 which incorporate spin and precession effects to different extents (Husa et al., 2016; Khan et al., 2016; Hannam et al., 2014; Pan et al., 2014; Taracchini et al., 2014; Babak et al., 2017). These models are restricted to quasi-circular orbits (i.e., they do not include orbital eccentricity) and 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., 2016a, 2017b), and given the large statistical uncertainties, the ignored effects in waveform modeling are not expected to cause a difference to 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.

5.5 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 ($\lesssim 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.

6 Conclusion and outlook

In this paper we have presented the first measurement of the Hubble constant 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 (Abbott et al., 2017a), along with constraints using galaxy catalogs for the BBH events observed by Advanced LIGO and Virgo in their first and second observing runs. We measure $H_{0}=$ $68.7^{+17.0}_{-7.8}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ ($68.3\%$ highest density posterior interval with a flat-in-log prior). This result is mainly dominated by the information from GW170817 with its counterpart, but does show a modest improvement with the inclusion of the BBHs. The BBHs contribute both from associated galaxy catalogs as well as via their observed luminosity distance distribution. Since the latter contribution is sensitive to the assumptions on the mass distribution and rate evolution, a more thorough treatment requires a marginalization over these unknown population parameters.

The contribution from events without counterparts is dominated by GW170814, for which the associated galaxy catalog is highly complete. This highlights the importance of deeper surveys and of dedicated EM follow-up of sky regions following GW triggers for a better $H_{0}$ measurement. With numerous anticipated detections in the upcoming observing runs with improved detector sensitivities (Abbott et al., 2018a, 2016b, 2016a, 2017c, 2019b, 2019a), these results pave the road towards an era of precision multimessenger cosmology to be performed with a multitude of sources, including both neutron star and black hole mergers, with or without transient EM counterparts.

Acknowledgments

Analyses in this paper made use of NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), Astropy (Robitaille et al., 2013; Price-Whelan et al., 2018), gwcosmo (https://git.ligo.org/lscsoft/gwcosmo), and IPython (Perez & Granger, 2007); plots were produced with Matplotlib (Hunter, 2007), and Seaborn (Waskom et al., 2020). This research made use of the K-corrections calculator service available at http://kcor.sai.msu.ru/.

The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as 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 Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, 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, the Vicepresidència i Conselleria d’Innovació, Recerca i Turisme and the Conselleria d’Educació i Universitat del Govern de les Illes Balears, the Conselleria d’Educació, Investigació, Cultura i Esport de la Generalitat Valenciana, the National Science Centre of Poland, the Swiss National Science Foundation (SNSF), the Russian Foundation for Basic Research, the Russian Science Foundation, the European Commission, 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 Lyon Institute of Origins (LIO), the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, Industry Canada and the Province of Ontario through the Ministry of Economic Development and Innovation, the Natural Science and Engineering Research Council Canada, the Canadian Institute for Advanced Research, the Brazilian Ministry of Science, Technology, Innovations, and Communications, 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 and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, INFN and CNRS for provision of computational resources.