A Gravitational-wave Measurement of the Hubble Constant Following the Second Observing Run of Advanced LIGO and Virgo

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,

$$\begin{eqnarray}&&\,\begin{array}{l}p({x}_{\mathrm{GW}},{x}_{\mathrm{EM}}| {D}_{\mathrm{GW}},{D}_{\mathrm{EM}},{H}_{0})\\ =\displaystyle \frac{p({x}_{\mathrm{GW}}| {H}_{0})p({x}_{\mathrm{EM}}| {H}_{0})}{p({D}_{\mathrm{EM}}| {D}_{\mathrm{GW}},{H}_{0})p({D}_{\mathrm{GW}}| {H}_{0})}.\end{array}\,\end{eqnarray} \tag{ 4 }$$

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₀) to be a constant that does not depend on H₀. This leads to

$$\begin{eqnarray}&&p({x}_{\mathrm{GW}},{x}_{\mathrm{EM}}| {D}_{\mathrm{GW}},{D}_{\mathrm{EM}},{H}_{0})\approx \displaystyle \frac{p({x}_{\mathrm{GW}}| {H}_{0})p({x}_{\mathrm{EM}}| {H}_{0})}{p({D}_{\mathrm{GW}}| {H}_{0})}.\end{eqnarray} \tag{ 5 }$$

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 $\bar{G})$:

$$\begin{eqnarray}&&\begin{array}{l}p({x}_{\mathrm{GW}}| {D}_{\mathrm{GW}},{H}_{0})\,=\,p({x}_{\mathrm{GW}}| G,{D}_{\mathrm{GW}},{H}_{0})p(G| {D}_{\mathrm{GW}},{H}_{0})\\ \,+\,p({x}_{\mathrm{GW}}| \bar{G},{D}_{\mathrm{GW}},{H}_{0})p(\bar{G}| {D}_{\mathrm{GW}},{H}_{0}).\end{array}\end{eqnarray} \tag{ 6 }$$

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 $\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_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:

$$\begin{eqnarray}&&\,\begin{array}{l}p({x}_{\mathrm{GW}}| G,{D}_{\mathrm{GW}},{H}_{0})\\ =\,\displaystyle \frac{\displaystyle \sum _{j=1}^{{N}_{\mathrm{gal}}}\int p({x}_{\mathrm{GW}}| {z}_{j},{{\rm{\Omega }}}_{j},{H}_{0})p(s| M({z}_{j},{m}_{j},{H}_{0}))p({z}_{j}){{dz}}_{j}}{\sum _{j=1}^{{N}_{\mathrm{gal}}}\int p({D}_{\mathrm{GW}}| {z}_{j},{{\rm{\Omega }}}_{j},{H}_{0})p(s| M({z}_{j},{m}_{j},{H}_{0}))p({z}_{j}){{dz}}_{j}}.\end{array}\,\end{eqnarray} \tag{ 7 }$$

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₀) is the absolute magnitude (for the given H₀), and p(s∣M(z_j , m_j , H₀)) 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₀)) as

$$\begin{eqnarray}p(s| M({z}_{j},{m}_{j},{H}_{0}))\propto \left\{\begin{array}{ll}\mathrm{constant} & \mathrm{if}\ \mathrm{unweighted}\\ L({M}_{j}({H}_{0})) & \mathrm{if}\ \mathrm{luminosity}\,\mathrm{weighted}.\end{array}\right.\end{eqnarray} \tag{ 8 }$$

The likelihood when the host galaxy is not in the catalog, $p({x}_{\mathrm{GW}}| \bar{G},{D}_{\mathrm{GW}},{H}_{0})$, is a ratio of marginalized integrals:

$$\begin{eqnarray}&&\begin{array}{l}p({x}_{\mathrm{GW}}| \bar{G},{D}_{\mathrm{GW}},{H}_{0})\\ =\displaystyle \frac{{\iiint }_{z({m}_{\mathrm{th}},M,{H}_{0})}^{\infty }p({x}_{\mathrm{GW}}| z,{\rm{\Omega }},{H}_{0})p(z)p({\rm{\Omega }})p(M| {H}_{0})p(s| M){dzd}{\rm{\Omega }}{dM}}{{\iiint }_{z({m}_{\mathrm{th}},M,{H}_{0})}^{\infty }p({D}_{\mathrm{GW}}| z,{\rm{\Omega }},{H}_{0})p(z)p({\rm{\Omega }})p(M| {H}_{0})p(s| M){dzd}{\rm{\Omega }}{dM}}.\end{array}\,\end{eqnarray} \tag{ 9 }$$

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_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,

$$\begin{eqnarray}&&\,\begin{array}{l}p(G| {D}_{\mathrm{GW}},{H}_{0})\\ =\displaystyle \frac{{\iiint }_{0}^{z({m}_{\mathrm{th}},M,{H}_{0})}p({D}_{\mathrm{GW}}| z,{\rm{\Omega }},{H}_{0})p(z)p({\rm{\Omega }})p(M| {H}_{0})p(s| M){dzd}{\rm{\Omega }}{dM}}{{\iiint }_{0}^{\infty }p({D}_{\mathrm{GW}}| z,{\rm{\Omega }},{H}_{0})p(z)p({\rm{\Omega }})p(M| {H}_{0})p(s| M){dzd}{\rm{\Omega }}{dM}},\\ \mathrm{and}\ p(\bar{G}| {D}_{\mathrm{GW}},{H}_{0})=1-p(G| {D}_{\mathrm{GW}},{H}_{0}).\end{array}\,\end{eqnarray} \tag{ 10 }$$

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

$$\begin{eqnarray}&&p(z)\propto \displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{{\rm{c}}}(z)}{{dz}}R(z).\end{eqnarray} \tag{ 11 }$$

Here, V_c(z) is the comoving volume as a function of redshift and the factor (1 + z)⁻¹ 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)=\mathrm{constant}$ throughout (other than in Section 5, where we consider an alternative redshift-dependent rate model). The prior on the GW sky location p(Ω) is taken to be uniform across the sky. The term p(M∣H₀) 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:

$$\begin{eqnarray}&&p(M| {H}_{0})\propto {10}^{-0.4(\alpha +1)(M-{M}^{* }({H}_{0}))}\exp \left[-{10}^{-0.4(M-{M}^{* }({H}_{0}))}\right].\end{eqnarray} \tag{ 12 }$$

Following Gehrels et al. (2016), we use B-band luminosity function parameters α = − 1.07 for the slope of the Schechter function and ${M}^{* }({H}_{0})=-19.7+5{\mathrm{log}}_{10}h$ for its characteristic absolute magnitude (with h ≡ H₀/100 km s⁻¹ Mpc⁻¹) throughout, unless otherwise specified. ¹⁹⁸ For the upper limits of integration over M, we choose the magnitude of the dimmest galaxies to be $-12.2+5{\mathrm{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₀) 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).

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₀), which requires removing the priors used for the parameter estimation. In particular, this means that the ${d}_{{\rm{L}}}^{2}$ prior must be removed before the priors on z and H₀ 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 source-frame mass distribution, while the data x_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_GW∣H₀) 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₁, m₂∣ μ), 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₀,

$$\begin{eqnarray}&&{m}_{\mathrm{1,2}}=\displaystyle \frac{{m}_{1,2}^{\det }}{1+z({d}_{{\rm{L}}},{H}_{0})}.\end{eqnarray} \tag{ 13 }$$

We marginalize out the source-frame mass dependence using the prior p(m₁, m₂∣ μ). The GW likelihood is a function of d_L and Ω as before, but now is also a function of H₀. 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₀, leading to additional constraining power (see Section 5.2 for details).

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² to as high as 1666 deg². A summary of the relevant parameters of all the GW detections are given in Table 1.

Table 1. Relevant Parameters of the O1 and O2 Detections: Network S/N for the Search Pipeline (PyCBC/GstLAL) in Which the Signal is the Loudest, 90% Sky Localization Region ΔΩ (deg²), Luminosity Distance d_L (Mpc, Median with 90% Credible Intervals), and Estimated Redshift z_event (median with a 90% Range Assuming Planck 2015 Cosmology) from Abbott et al. (2019b)

Event	S/N	ΔΩ/deg2	dL/Mpc	zevent	V/Mpc3	Galaxy Catalog	Number of Galaxies	mth	p(G∣zevent, DGW)
GW150914	24.4	182	${440}_{-170}^{+150}$	${0.09}_{-0.03}^{+0.03}$	3.5 × 106	GLADE	3910	17.92	0.42
GW151012	10.0	1523	${1080}_{-490}^{+550}$	${0.21}_{-0.09}^{+0.09}$	5.8 × 108	GLADE	78195	17.97	0.01
GW151226	13.1	1033	${450}_{-190}^{+180}$	${0.09}_{-0.04}^{+0.04}$	2.4 × 107	GLADE	27677	17.93	0.41
GW170104	13.0	921	${990}_{-430}^{+440}$	${0.20}_{-0.08}^{+0.08}$	2.4 × 108	GLADE	42221	17.76	0.01
GW170608	15.4	392	${320}_{-110}^{+120}$	${0.07}_{-0.02}^{+0.02}$	3.4 × 106	GLADE	6267	17.84	0.60
GW170729	10.8	1041	${2840}_{-1360}^{+1400}$	${0.49}_{-0.21}^{+0.19}$	8.7 × 109	GLADE	77727	17.82	< 0.01
GW170809	12.4	308	${1030}_{-390}^{+320}$	${0.20}_{-0.07}^{+0.05}$	9.1 × 107	GLADE	18749	17.62	< 0.01
GW170814	16.3	87	${600}_{-220}^{+150}$	${0.12}_{-0.04}^{+0.03}$	4.0 × 106	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 × 107	GLADE	1059	17.51	< 0.01
GW170823	11.5	1666	${1940}_{-900}^{+970}$	${0.35}_{-0.15}^{+0.15}$	3.5 × 109	GLADE	117680	17.98	< 0.01

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.

Download table as:  ASCIITypeset image

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

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 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 (Abbott et al. 2018b; Drlica-Wagner et al. 2018) 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 the SDSS (including curated versions such as the Gravitational Wave Events in Sloan (GWENS) catalog, M. Rahman et al. 2019, in preparation), could be used to improve the current analysis for events that 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), ¹⁹⁹ 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

3.2.2. DES Year 1

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₀ 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₀. We additionally show as the vertical lines in Figure 1 the median redshift for each event z_event (calculated assuming a 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.

Zoom In Zoom Out Reset image size

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 $0.05{L}_{g}^{* }$ within the corresponding event's 99% sky localization region defined as

$$\begin{eqnarray}&&{p}_{\mathrm{cat}}(z)\equiv \int p({x}_{\mathrm{GW}}| {\rm{\Omega }}){p}_{0}(z,{\rm{\Omega }})d{\rm{\Omega }},\end{eqnarray} \tag{ 14 }$$

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₀(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₀ 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.

Zoom In Zoom Out Reset image size

3.4.1. Redshift Uncertainties

3.4.2. Source-frame Luminosities

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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_⊙, consistent with all the considered BBHs for all values of H₀ within the prior range. Reducing this cutoff to a slightly restrictive ${M}_{\max }=50\,{M}_{\odot }$ (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₀ 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)=\mathrm{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)³ 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)⁻³ 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₀, 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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-of-catalog 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.

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 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 in the current measurement of H₀. 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.

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