A Standard Siren Measurement of the Hubble Constant from GW170817 without the Electromagnetic Counterpart

The first multimessenger detection of a binary neutron star (BNS) merger, GW170817, by LIGO (Aasi et al. 2015) and Virgo (Acernese et al. 2015) enabled the first standard siren measurement of the Hubble constant, H₀, ushering in the era of gravitational-wave (GW) cosmology (Abbott et al. 2017a, 2017c, 2017d). This H₀ measurement combined the luminosity distance to the source, as measured from the GW signal (Schutz 1986), with the known redshift of the host galaxy, NGC 4993. NGC 4993 was identified as the unique host galaxy following the discovery of an optical transient located only ∼10 arcsec from NGC 4993 (Abbott et al. 2017d; Coulter et al. 2017; Soares-Santos et al. 2017). The probability of a chance coincidence between the GW signal and the optical transient was estimated to be ≲0.5% (Soares-Santos et al. 2017), and the probability of a chance association between the optical transient and NGC 4993 is ≲0.004% (Abbott et al. 2017a).

The original proposal by Schutz (1986) to measure the Hubble constant with GW detections of compact binary mergers did not involve electromagnetic counterparts. Instead, Schutz considered bright galaxies in the GW localization region as potential hosts to the merger. Each galaxy provides a redshift that, when combined with the GW-measured luminosity distance, gives a separate estimate of H₀. The final H₀ measurement from a single event is the sum of all contributions from the individual galaxies. The first detailed exploration of this method on simulated data, and with the first use of a galaxy catalog (Aihara et al. 2011), was by Del Pozzo (2012). An up-to-date forecast incorporating realistic detection rates, galaxy peculiar velocities, large-scale structure, and additional considerations can be found in Chen et al. (2018). We refer to this approach of measuring H₀ as the "statistical" method (Schutz 1986; MacLeod & Hogan 2008; Petiteau et al. 2011; Chen et al. 2018), compared with the "counterpart" method in which an electromagnetic (EM) counterpart provides a unique host galaxy association. In the limit where the GW event is so well-localized that there is only one potential host galaxy in the GW localization error box (Chen & Holz 2016), the statistical method reduces to the counterpart method. In the opposite limit, where the GW event is poorly localized, there are so many potential host galaxies that the distinct peaks from individual galaxies are washed out, and the H₀ measurement is uninformative (Chen et al. 2018).

The statistical approach may be the only way to do standard siren science with binary black holes, because they are not expected to have EM counterparts. We emphasize that although the statistical measurements for a given event are inferior to the counterpart case, combining many of these measurements leads to increasingly precise constraints (Schutz 1986; Del Pozzo 2012; Chen et al. 2018; Nair et al. 2018). In ground-based gravitational wave detector networks, the rate of detection of binary black holes is significantly higher than that for neutron stars (Abbott et al. 2016, 2017b, 2017c), although the higher rate is not expected to compensate for the inferior constraints (Chen et al. 2018). Nonetheless, the black hole systems can be observed to much higher redshifts, potentially providing constraints on the evolution history of the universe out past the turnover between dark matter and dark energy domination (Del Pozzo 2012; Dominik et al. 2015; Belczynski et al. 2016; Fishbach et al. 2018). Because these systems are farther away, however, it will be a greater challenge to supply a sufficiently complete galaxy catalog.

In this paper we carry out a measurement of H₀ using the GW data from GW170817 and a catalog of potential host galaxies within the GW localization region. In other words, we explore how tight the H₀ measurement from GW170817 would have been if an EM counterpart had not been detected or if a unique host galaxy had not been identified. We present our methods in Section 2, a discussion of the galaxy selection in Section 3, a discussion of the GW constraints in Section 4, results in Section 5, and conclude in Section 6.

We follow the statistical framework presented in Chen et al. (2018; see also Del Pozzo 2012; Gray et al. 2019). We include the role of GW selection effects, galaxy catalog incompleteness, galaxy luminosities, and redshift uncertainties in our analysis. The posterior on H₀ given the GW and EM data, x_GW and x_EM, is:

$$\begin{eqnarray}\begin{array}{rcl}p({H}_{0}| {x}_{\mathrm{GW}},{x}_{\mathrm{EM}}) & = & \displaystyle \frac{{p}_{0}({H}_{0})}{\beta ({H}_{0})}\displaystyle \int p({x}_{\mathrm{GW}}| {\hat{D}}_{L}(z,{H}_{0}),{\rm{\Omega }})\\ & & \times \ p({x}_{\mathrm{EM}}| z,{\rm{\Omega }}){p}_{0}(z,{\rm{\Omega }})d{\rm{\Omega }}{dz},\end{array}\end{eqnarray} \tag{ 1 }$$

where ${\hat{D}}_{L}(z,{H}_{0})$ is the luminosity distance of a source at redshift z for a given H₀ (fixing other cosmological parameters to the Planck values; Ade et al. 2016)⁸⁶ , Ω is the sky position, and β(H₀) is a normalization term to ensure that the likelihood normalizes to 1 when integrated over all detectable GW and EM data sets (Mandel et al. 2016). The term p₀(H₀) represents the prior on the Hubble constant. A detailed derivation of Equation (1), including the role of the normalization term β(H₀), is provided in Appendix A.

As first emphasized by Schutz (1986), the GW signal from a compact binary coalescence allows for a direct measurement of the distance to the source, as well as its sky location. This measurement is represented in the GW likelihood term, $p({x}_{\mathrm{GW}}| {D}_{L},{\rm{\Omega }})$, which is the probability of the GW data in the presence of signal from a compact binary with parameters D_L and Ω marginalized over the other parameters of the signal (including the inclination angle, component masses, spins, and/or tides). The corresponding posterior $p({D}_{L},{\rm{\Omega }}| {x}_{\mathrm{GW}})\propto p({x}_{\mathrm{GW}}| {D}_{L},{\rm{\Omega }}){p}_{0}({D}_{L},{\rm{\Omega }})$ is summarized in the three-dimensional sky map, which provides a fit to the posterior samples provided by the GW parameter estimation pipeline LALInference (Veitch et al. 2015; Singer et al. 2016a, 2016b). For this analysis, we use the publicly released three-dimensional sky map from Abbott et al. (2019; see Sections 4 and 5). To get the likelihood from the posterior probability, we must first divide out the default "volumetric" distance prior, ${p}_{0}({D}_{L},{\rm{\Omega }})\propto {D}_{L}^{2}$.

Meanwhile, the EM likelihood term $p({x}_{\mathrm{EM}}| z,{\rm{\Omega }})$ is the probability of the electromagnetic data in the presence of signal from a compact binary with parameters z and Ω. In the absence of an EM counterpart and/or a host galaxy identification, we assume the measurement $p({x}_{\mathrm{EM}}| z,{\rm{\Omega }})$ is completely uninformative, and set:

$$\begin{eqnarray}&&p({x}_{\mathrm{EM}}| z,{\rm{\Omega }})\propto 1.\end{eqnarray} \tag{ 2 }$$

In this case, the redshift information enters only through the prior term, ${p}_{0}(z,{\rm{\Omega }})$, which we take to be a galaxy catalog. The detection of an electromagnetic counterpart typically results in $p({x}_{\mathrm{EM}}| z,{\rm{\Omega }})$ being strongly peaked around some $\hat{{\rm{\Omega }}}$ allowing the identification of a host galaxy. We note that in some cases an optical transient may be identified, but it may not be possible to uniquely identify the associated host galaxy. In these circumstances one could perform a pencil-beam survey of the region surrounding the transient (e.g., at distances of ≲100 kpc from the line of sight to the transient), and sharply reduce the relevant localization volume (Chen et al. 2018). This reduces the number of potential host galaxies, and thereby improves the measurement.

The galaxy catalog term p₀(z, Ω) is given by:

$$\begin{eqnarray}&&{p}_{0}(z,{\rm{\Omega }})={{fp}}_{\mathrm{cat}}(z,{\rm{\Omega }})+(1-f){p}_{\mathrm{miss}}(z,{\rm{\Omega }}),\end{eqnarray} \tag{ 3 }$$

where ${p}_{\mathrm{cat}}$ is a catalog of known galaxies, p_miss represents the distribution of missing galaxies, and f denotes the overall completeness fraction of the catalog. The contribution from the known galaxies is:

$$\begin{eqnarray}&&{p}_{\mathrm{cat}}(z,{\rm{\Omega }})=\displaystyle \sum _{i}^{{N}_{\mathrm{gal}}}{w}_{i}N({\bar{z}}_{i},{\sigma }_{z};z)N({\bar{{\rm{\Omega }}}}_{i},{\sigma }_{{\rm{\Omega }}};{\rm{\Omega }}),\end{eqnarray} \tag{ 4 }$$

where ${\bar{z}}_{i},{\bar{{\rm{\Omega }}}}_{i}$ denotes the (peculiar velocity-corrected) "Hubble" redshifts and sky coordinates of all galaxies in the catalog, and $N(\mu ,\sigma ;x)$ denotes the normal probability density function with mean μ and standard deviation σ evaluated at x. To account for peculiar velocity uncertainties, which can be significant for nearby sources, we assume that the true Hubble velocity is normally distributed about the measured Hubble velocity with an uncertainty of cσ_z (Scolnic et al. 2018). On the other hand, the uncertainty on the sky coordinates of galaxies in the catalog is negligible for our purposes, so we always approximate $N(\bar{{\rm{\Omega }}},{\sigma }_{{\rm{\Omega }}};{\rm{\Omega }})$ by $\delta (\bar{{\rm{\Omega }}}-{\rm{\Omega }})$.

The weights w_i can be chosen to reflect the a priori belief that a galaxy could host a GW source. For example, setting all weights to ${w}_{i}=\tfrac{1}{{N}_{\mathrm{gal}}}$ corresponds to equal probability for each galaxy to host a gravitational wave source. In general, because we might expect that the BNS rate is traced by some combination of stellar mass and/or star formation rate (Phinney 1991; Leibler & Berger 2010; Fong & Berger 2013; Chruslinska et al. 2018), we may assign unequal weights to galaxies based on these (or any other relevant observable) quantities, ensuring that the weights sum to unity. In the following, we use a galaxy's B-band luminosity as a proxy for its star formation rate, and its K-band luminosity as a proxy for its total stellar mass (Bell et al. 2003; Singer et al. 2016a); these are very rough estimates, but serve to demonstrate the method. In these cases, we apply weights proportional to B-band or K-band luminosity, ${w}_{i}\propto {L}_{B}^{i}$ or ${w}_{i}\propto {L}_{K}^{i}$, and explore the dependence of the result on those of weightings.

To calculate the term ${p}_{\mathrm{miss}}$ in Equation (3), we assume that on large scales, the distribution of galaxies, ${p}_{0}(z,{\rm{\Omega }})$ is uniform in comoving volume. Let ${p}_{\mathrm{vol}}(z,{\rm{\Omega }})$ denote the cosmologically homogeneous and isotropic distribution normalized over the volume contained in the range ${z}_{\min }\lt z\lt {z}_{\max }$ considered in our analysis. (The result does not depend on the choice of z_min or z_max provided that the interval encompasses all possible redshifts of the source for all allowed values of H₀.) Assuming all galaxies are weighted equally, the distribution of missing galaxies is written as:

$$\begin{eqnarray}&&{p}_{\mathrm{miss}}(z,{\rm{\Omega }})=\displaystyle \frac{\left[1-{P}_{\mathrm{complete}}(z)\right]{p}_{\mathrm{vol}}(z,{\rm{\Omega }})}{(1-f)},\end{eqnarray} \tag{ 5 }$$

where ${P}_{\mathrm{complete}}(z)$ is the probability that a galaxy at redshift z is in the catalog, and the completeness fraction f is given by:

$$\begin{eqnarray}&&f={\int }_{{z}_{\min }}^{{z}_{\max }}{P}_{\mathrm{complete}}(z){p}_{\mathrm{vol}}(z,{\rm{\Omega }})\,{dz}\,d{\rm{\Omega }}.\end{eqnarray} \tag{ 6 }$$

We can similarly add galaxy weightings to an incomplete catalog by computing the luminosity distribution of the "missing galaxies" as a function of redshift, $p(L| z,\mathrm{missing})$. We calculate this distribution by assuming that the luminosities of the missing galaxies together with those in the catalog are distributed according to the Schechter function. Then, the weights of the missing galaxies are given by:

$$\begin{eqnarray}&&w(z)\propto \int {Lp}(L| z,\mathrm{missing})\,{dL},\end{eqnarray} \tag{ 7 }$$

and, weighting each missing galaxy by its luminosity, Equation (5) becomes:

$$\begin{eqnarray}&&{p}_{\mathrm{miss}}(z,{\rm{\Omega }})\propto w(z){p}_{\mathrm{miss}}(z,{\rm{\Omega }}).\end{eqnarray} \tag{ 8 }$$

Note that we have assumed that the coverage of the catalog is uniform over the sky and so P_complete is independent of Ω. (This is true over the relevant sky area for the present analysis, but the method can be easily generalized to add an Ω-dependence.) Alternate approaches of taking into account the incompleteness of galaxy catalogs are being explored in Gray et al. (2019). However, in the present case of a single nearby source where the catalog is largely complete, the differences in results from the various approaches are small, and in particular well within the statistical uncertainties. In Section 3, the completeness function P_complete is estimated for the galaxy catalog used in the analysis.

To demonstrate the statistical method, we apply the analysis described above to 249 simulated BNS GW detections from the First Two Years (F2Y) catalog (Singer et al. 2014) and the MICE simulated galaxy catalog (Carretero et al. 2015, 2017; Crocce et al. 2015; Fosalba et al. 2015a, 2015b; Hoffmann et al. 2015). We assign each BNS detection from the F2Y 2016 scenario (roughly corresponding to O2) to a galaxy in the MICE catalog with a redshift that matches the injected distance and assumed H₀ value (H₀ = 70 km s⁻¹ Mpc⁻¹). For each event, we rotate the sky coordinates of the galaxies in the catalog so that the sky position of the host galaxy matches the true sky position of the BNS injection. We then carry out the statistical method using the three-dimensional sky map for each mock BNS and the galaxies in MICE, assuming no peculiar velocities or incompleteness, and assigning weights to the galaxies in MICE so that the redshift distribution matches the injected redshift distribution of the F2Y data set, p(z) ∝ z². This last step is necessary in order to ensure that the selection effects are incorporated consistently between the injections and the likelihood. The results are shown in Figure 1. Even in the best-case scenario of perfectly known galaxy redshifts and a complete catalog, the H₀ posteriors from most individual events are nearly flat over the prior range. Combining the 249 individual events, the final H₀ posterior is ${H}_{0}\,={70.1}_{-1.9}^{+1.9}$ km s⁻¹ Mpc⁻¹ (68.3% credible interval), corresponding to a convergence rate of $\sim 40 \% /\sqrt{N}$ where N is the number of events, consistent with Chen et al. (2018).  For counterpart standard sirens the convergence is $\sim 15 \% /\sqrt{N}$ (Chen et al. 2018). As is visible in the figure, we confirm that the method is unbiased, with the result for large numbers of detections approaching the true value of H₀. We note that most of the simulated detections in the F2Y data set have much larger localization volumes than GW170817, which was an unusually loud, nearby source that was detected while all three detectors were operational. Therefore, we expect the statistical H₀ measurement from GW170817 to be unusually informative compared to an average event. We quantify this expectation in Section 5.

Zoom In Zoom Out Reset image size

To measure H₀ statistically with GW170817, we use version 2.3 of the GLADE galaxy catalog to construct our redshift prior in Equation (3) (Dálya et al. 2018). GLADE provides galaxy redshifts in the heliocentric frame, corrected for peculiar motions using the peculiar velocity catalog of Carrick et al. (2015). For galaxies that are also listed in the group catalog of Kourkchi & Tully (2017), as identified by a common Principal Galaxy Catalog identifier, we apply an additional correction to correct their velocities to the radial velocity of the group. We assume the group velocity is given by the unweighted mean of the velocities of all member galaxies, although we note that for the dominant group containing NGC 4993, careful group modeling has been done (Hjorth et al. 2017). Finally, we correct all heliocentric velocities to the reference frame of the cosmic microwave background (Hinshaw et al. 2009) and assign a 200 km s⁻¹ Gaussian uncertainty to the "Hubble velocity" of each galaxy in the catalog (corrected by all peculiar motions; Carrick et al. 2015; Scolnic et al. 2018).

GLADE also provides luminosity information for galaxies, listing apparent magnitudes in the B-, J-, H-, and K-bands. We use the reported B-band luminosities to characterize the completeness of the catalog (a small fraction of galaxies do not have B-band apparent magnitudes reported in the catalog; we remove these galaxies from our analysis, assuming that their magnitudes are below our adopted luminosity cutoff). Following Gehrels et al. (2016) and Arcavi et al. (2017), we adopt B-band Schechter function parameters ${\phi }^{\star }=5.5\times {10}^{-3}{h}_{0.7}^{3}$ Mpc${}^{-3}$, ${\alpha }_{B}=-1.07$, ${L}_{B}^{\star }=2.45\times {10}^{10}{h}_{0.7}^{-2}{L}_{B,\odot }$ throughout. The corresponding characteristic absolute magnitude is ${M}_{B}^{\star }=-20.47+5{\mathrm{log}}_{10}{h}_{0.7}$. We will also consider the K-band magnitudes reported in GLADE when applying galaxy weights, and we use the K-band Schechter function parameters of ${\alpha }_{K}=-1.02$, ${M}_{K}^{\star }=-23.55\,+5{\mathrm{log}}_{10}{h}_{0.7}$ (Lu et al. 2016).

Figure 2 summarizes the completeness of GLADE as a function of redshift. We find that GLADE is complete up to redshifts z ∼ 0.06 for galaxies brighter than $\sim 0.626{h}_{0.7}^{-2}\ {L}_{B}^{\star }$, corresponding to about 0.66 of the Milky Way luminosity for h_0.7 = 1. Galaxies brighter than $\sim 0.626{L}_{B}^{\star }$ make up half of the total luminosity for the given Schechter function parameters. We find that for $z\lesssim 0.03$, GLADE is complete for galaxies down to 2.5 times dimmer, or $\sim 0.25{L}_{B}^{\star }$, corresponding to ${M}_{B}=-18.96+5{\mathrm{log}}_{10}{h}_{0.7}$ (see also Figure 2 of Arcavi et al. 2017). Such galaxies make up 75% of the total B-band luminosity. If we consider galaxies down to $\sim 0.05{L}_{B}^{\star }$ (${M}_{B}=-17.22+5{\mathrm{log}}_{10}{h}_{0.7}$), GLADE is $\sim 70 \%$ complete at z ∼ 0.03, and even if we consider galaxies down to $\sim 0.01{L}_{B}^{\star }$ (${M}_{B}=-15.47+5{\mathrm{log}}_{10}{h}_{0.7}$), including 99% of the total B-band luminosity, GLADE is ≳80% complete for z ≲ 0.01, and ∼40% complete at z ∼ 0.03. In the K-band, we find that with our assumed K-band Schechter function parameters, GLADE is complete up to z ∼ 0.045 for galaxies with ${L}_{K}\gt 0.36{L}_{K}^{\star }$, which contain 70% of the total K-band luminosity, and up to z ∼ 0.03 for galaxies with ${L}_{K}\gt 0.1{L}_{K}^{\star }$, which contain 90% of the total luminosity. For galaxies brighter than ${L}_{K}=0.005{L}_{K}^{\star }$, which make up more than 99% of the total K-band luminosity, GLADE is ∼70% complete at z = 0.01.

Zoom In Zoom Out Reset image size

From the GW data alone, GW170817 is the best-localized GW event to date. The original analysis by the LIGO-Virgo collaboration (Abbott et al. 2017c) reported a 90% localization area of 28 deg² and a 90% localization volume of 380 Mpc³ (assuming Planck cosmology; Ade et al. 2016), while the most recent analysis (Abbott et al. 2019) improves the 90% localization area to 16 deg² and the 90% volume to 215 Mpc³. We use this updated three-dimensional sky map (Singer et al. 2016a, 2016b) from Abbott et al. (2019) throughout.⁸⁷ Figure 3 shows the two-dimensional sky map together with the galaxies in the GLADE catalog within the localization region. Figure 4 shows that, although there are a total of 408 galaxies within the 99% localization area (see Figure 3), most of the galaxies with high sky-map probability come in a few distinct groups: a dominant group at z ∼ 0.01 regardless of the assumed luminosity threshold, followed by a secondary group at z ∼ 0.006 containing only moderately faint galaxies. Therefore, there are only a few distinct redshifts that can possibly correspond to the measured distance of GW170817, and we expect that combining the galaxy catalog with the GW localization will yield an informative measurement of the Hubble constant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The three-dimensional sky map also provides an approximation to the luminosity distance posterior along each line of sight. As usual, the distance to GW170817 is determined directly from the gravitational waves, and is calibrated by general relativity (Schutz 1986). No distance ladder is required.

We combine the GW distance posterior for GW170817 with the redshift for each potential host galaxy within the localization region. As detailed in Section 2, each galaxy produces a posterior probability for H₀, and we combine these estimates among all the galaxies in the localization region to arrive at a final estimate for H₀. We adopt a flat prior in H₀ over the range 10–220 km s⁻¹ Mpc⁻¹. The results are presented in Figure 5. Because the galaxies are predominantly found in one galaxy group at $z\sim 0.01$, the H₀ posterior shows a clear peak at ${H}_{0}\approx 76$ km s⁻¹ Mpc⁻¹. And because NGC 4993, the true galaxy host of GW170817, is a member of the group at $z\sim 0.01$, we should not be surprised to learn that the peak in H₀ is consistent with the H₀ estimate from the GW170817 standard siren measurement including the counterpart (Abbott et al. 2017a). Because this analysis has been performed on a three-dimensional sky map using an approximation to the distance posteriors, rather than using the full three-dimensional LIGO/Virgo posteriors, the results do not agree precisely with those of Abbott et al. (2017a), and in particular the position of the peak in Figure 5 is at H₀ = 76 km s⁻¹ Mpc⁻¹ instead of H₀ = 70 km s⁻¹ Mpc⁻¹. This is because our three-dimensional sky map approximates the distance posterior along each line of sight by a simple two-parameter Gaussian fit (see Equation (1) of Singer et al. 2016a), which is an imperfect approximation to the true, asymmetric distance posterior (Chen & Holz 2017; Del Pozzo et al. 2018). On the other hand, the analysis in Abbott et al. (2017a) utilizes the full distance posterior along the line of sight to NGC 4993 rather than the Gaussian approximation.

Zoom In Zoom Out Reset image size

Figure 5 shows four different posterior probability distributions, each using a different threshold for the galaxy catalog. In the "assuming counterpart" case, NGC 4993 (which is assumed to be the true host galaxy to GW170817) is given a weight of 1, and all the other galaxies in the localization volume are given a weight of 0. We find ${H}_{0}={76}_{-13}^{+19}$ km s⁻¹ Mpc⁻¹ (maximum a posteriori and 68.3% highest density posterior interval) for our default flat prior, or ${H}_{0}={74}_{-12}^{+18}$ km s⁻¹ Mpc⁻¹ for a flat-in-log prior (the prior choice in Abbott et al. 2017a). This peak is slightly shifted compared to the result presented in Abbott et al. (2017a), ${H}_{0}={70}_{-8}^{+12}$ km s⁻¹ Mpc⁻¹, due to the usage of the Gaussian fit to the distance posterior found in the three-dimensional sky map as discussed above.

The other curves in Figure 5 assume different limiting thresholds for what constitutes a potential host galaxy. For a luminosity threshold of $L\gt 0.626{L}_{B}^{\star }$, we find ${H}_{0}={77}_{-18}^{+37}$ km s⁻¹ Mpc⁻¹.⁸⁸ As the threshold is lowered, additional galaxies fall into the sample, and the H₀ posterior is broadened. For a limiting B-band magnitude of $0.25{L}_{B}^{\star }$, we need to account for the incompleteness of the galaxy catalog at redshifts $z\gtrsim 0.03$, and for $0.01{L}_{B}^{* }$, we need to account for the incompleteness at $z\gtrsim 0.01$, as described in Section 2. The incompleteness correction leads to a slight additional broadening of the H₀ posterior, but the clear peak at H₀ ≈ 76 km s⁻¹ Mpc⁻¹ remains: we find ${H}_{0}={74}_{-24}^{+45}$ km s⁻¹ Mpc⁻¹ for a luminosity threshold of $L\gt 0.25{L}_{B}^{\star }$ ${H}_{0}={76}_{-23}^{+48}$ km s⁻¹ Mpc⁻¹ for a luminosity threshold of $L\gt 0.01{L}_{B}^{\star }$. This peak is the result of the galaxy group at $z\sim 0.01$, of which NGC 4993 is a member.

The curves in Figure 6 weight each galaxy by its B-band luminosity (a proxy for its recent star formation history; right) or its K-band luminosity (a proxy for its stellar mass; left). The peak at ${H}_{0}\approx 76$ km s⁻¹ Mpc⁻¹ becomes more pronounced when galaxies are weighted by their luminosity, as the group containing NGC 4993 consists of many bright, mostly red galaxies. If we assume that the probability of hosting a BNS merger is proportional to a galaxy's B-band luminosity, the posterior on H₀ tightens from ${H}_{0}\in [53,124]$ km s⁻¹ Mpc⁻¹ (68.3% highest density posterior interval) when applying equal weights to all galaxies brighter than 0.01 ${L}_{B}^{\star }$ to ${H}_{0}\in [54,120]$ km s⁻¹ Mpc⁻¹. Applying K-band luminosity weights to galaxies brighter than 0.005 ${L}_{K}^{\star }$, the 68.3% posterior interval tightens from ${H}_{0}\in [61,137]$ km s⁻¹ Mpc⁻¹ to ${H}_{0}\in [57,118]$ km s⁻¹ Mpc⁻¹. Although these results are suggestive that weighting by stellar-mass or star formation rate may lead to faster convergence, the properties of BNS host galaxies are still uncertain, and it is impossible to establish this definitively with a single event. As the source sample increases it is expected to relate to some combination of these quantities, and incorporating these trends will lead to improvements in the statistical H₀ analysis.

Zoom In Zoom Out Reset image size

In order to quantify the degree of information in the GW170817 H₀ posterior compared to an "average" event as expected from the F2Y data set, we consider the difference in the Shannon entropy between the flat prior and the posterior (see Appendix B; Shannon 1948). We compare this measure of information for the statistical GW170817 H₀ posterior to the individual statistical H₀ posteriors from each of the simulated BNS events in Section 2. We find that for a flat prior in the (relatively narrow) range ${H}_{0}\in [40,100]$ km s⁻¹ Mpc⁻¹, the information gained by applying the statistical method to GW170817 is 0.34 bits. Meanwhile, the median information in an individual posterior shown in Figure 1 is only 0.047 bits, so that GW170817 is in the top ∼10% of informative events, even under optimistic assumptions for the simulated detections (i.e., complete galaxy catalogs and perfect redshift measurements). As expected, GW170817 provides an unusually good statistical H₀ constraint.

For the purposes of this calculation, we use the K-band luminosity-weighted ${L}_{K}\gt 0.1{L}_{K}^{\star }$ posterior shown in the right panel of Figure 6 as a representative posterior for the statistical GW170817 H₀ measurement. Over the wider prior ${H}_{0}\in [10,220]$ km s⁻¹ Mpc⁻¹ shown, the information difference between the posterior and the prior is 0.67 bits. The counterpart GW170817 H₀ measurement (dotted–dashed pink curve in Figure 5) has an information gain of 1.55 bits with respect to the wide prior.

We perform a statistical standard siren measurement of the Hubble constant with GW170817. This analysis is the first application of the measurement originally proposed over 30 yr ago by Schutz (1986). We find that the excellent localization of GW170817, together with the large-scale structure causing galaxies to cluster into distinct groups, enables an informative measurement of H₀ even in the absence of a unique host galaxy identification. Including generic and flexible assumptions regarding the luminosities of BNS host galaxies, we find a peak at ${H}_{0}\approx 76$ km s⁻¹ Mpc⁻¹ at ∼2.4–3.7 times the prior probability density. We find the possibility of improved constraints when weighting the potential host galaxies by stellar mass and star formation rate. Including all galaxies brighter than 0.01 ${L}_{B}^{\star }$ (including 99% of the total blue luminosity) we find ${H}_{0}={76}_{-23}^{+48}$ km s⁻¹ Mpc⁻¹, or ${H}_{0}\,={76}_{-21}^{+45}$ km s⁻¹ Mpc⁻¹ when applying B-band luminosity weights (a proxy for star formation rate). Weighting all galaxies brighter than 0.005 ${L}_{K}^{\star }$ by their K-band luminosity (a proxy for stellar mass), we find ${H}_{0}={76}_{-19}^{+41}$ km s⁻¹ Mpc⁻¹.

Although this statistical standard siren measurement of H₀ is less precise than the counterpart measurement of ${H}_{0}\,={76}_{-13}^{+19}$ km s⁻¹ Mpc⁻¹ (for a flat prior and utilizing the distance ansatz in the three-dimensional sky map; see Section 5), it nonetheless shows that interesting constraints on cosmological parameters are possible from GW sources even in the absence of an optical counterpart and an identification of the unique host galaxy (Schutz 1986; Del Pozzo 2012; Chen et al. 2018; Gray et al. 2019). Although detailed studies find that the measurement of cosmological parameters from the counterpart approach is likely to surpass the statistical approach (Chen et al. 2018), the statistical approach offers an important cross-validation of the counterpart standard siren measurements. Furthermore, the statistical approach holds particular promise for binary black hole sources, which are detected at higher rates than BNS systems and are expected to lack electromagnetic counterparts. The inferior quality of the individual H₀ measurements for binary black holes (because of the larger localization volumes) may be compensated for by the improved quantity due to the higher detection rates. The binary black holes will also be detected at much greater distances, and in addition to measuring H₀ may constrain additional cosmological parameters such as the equation of state of the dark energy.

We thank Chihway Chang, Alison Coil, and Risa Wechsler for valuable discussions about galaxy properties. M.F. was supported by the NSF Graduate Research Fellowship Program under grant DGE-1746045. M.F. and D.E.H. were supported by NSF grant PHY-1708081. They were also supported by the Kavli Institute for Cosmological Physics at the University of Chicago through NSF grant PHY-1125897 and an endowment from the Kavli Foundation. D.E.H. also gratefully acknowledges support from a Marion and Stuart Rice Award. The authors 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. We are grateful to the LIGO and Virgo collaborations for releasing sky map data for GW170817 at dcc.ligo.org/LIGO-P1800061/public. This work has made use of CosmoHub. CosmoHub has been developed by the Port d'Informació Cientfica (PIC), maintained through a collaboration of the Institut de Fsica d'Altes Energies (IFAE) and the Centro de Investigaciones Energticas, Medioambientales y Tecnológicas (CIEMAT), and was partially funded by the "Plan Estatal de Investigación Cientfica y Tcnica y de Innovación" program of the Spanish government. We also acknowledge the First Two Years data release (https://www.ligo.org/scientists/first2years/).

In this appendix we derive the H₀ posterior probability distribution function from Equation (1). We write the likelihood for GW and EM data, ${x}_{\mathrm{GW}}$ and ${x}_{\mathrm{EM}}$, given a value of H₀ as:

$$\begin{eqnarray}&&p({x}_{\mathrm{GW}},{x}_{\mathrm{EM}}| {H}_{0})=\displaystyle \frac{\int p({x}_{\mathrm{GW}},{x}_{\mathrm{EM}},{D}_{L},{\rm{\Omega }},z| {H}_{0})\,{{dD}}_{L}\,d{\rm{\Omega }}\,{dz}}{\beta ({H}_{0})},\end{eqnarray} \tag{ 9 }$$

and factor the numerator as:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int p({x}_{\mathrm{GW}},{x}_{\mathrm{EM}},{D}_{L},{\rm{\Omega }},z| {H}_{0})\,{{dD}}_{L}\,d{\rm{\Omega }}\,{dz}=\displaystyle \int \\ \times \ p({x}_{\mathrm{GW}}| {D}_{L},{\rm{\Omega }})p({x}_{\mathrm{EM}}| z,{\rm{\Omega }})p({D}_{L}| z,{H}_{0}){p}_{0}(z,{\rm{\Omega }}){{dD}}_{L}\,d{\rm{\Omega }}{dz}\\ =\displaystyle \int p({x}_{\mathrm{GW}}| {D}_{L},{\rm{\Omega }})p({x}_{\mathrm{EM}}| z,{\rm{\Omega }})\\ \times \ \delta ({D}_{L}-{\hat{D}}_{L}(z,{H}_{0})){p}_{0}(z,{\rm{\Omega }})\,{{dD}}_{L}\,d{\rm{\Omega }}\,{dz}\\ =\displaystyle \int p({x}_{\mathrm{GW}}| {\hat{D}}_{L}(z,{H}_{0}),{\rm{\Omega }})p({x}_{\mathrm{EM}}| z,{\rm{\Omega }}){p}_{0}(z,{\rm{\Omega }})\,d{\rm{\Omega }}\,{dz}.\end{array}\end{eqnarray} \tag{ 10 }$$

The H₀ posterior is related to the likelihood in Equation (10) by a prior:

$$\begin{eqnarray}&&p({H}_{0}| {x}_{\mathrm{GW}},{x}_{\mathrm{EM}})={p}_{0}({H}_{0})p({x}_{\mathrm{GW}},{x}_{\mathrm{EM}}| {H}_{0}).\end{eqnarray} \tag{ 11 }$$

This equation is identical to Equation (1) in the main text. The normalization term $\beta ({H}_{0})$ is given by (see Mandel et al. 2016; Chen et al. 2018):

$$\begin{eqnarray}\begin{array}{rcl}\beta ({H}_{0}) & = & {\displaystyle \int }_{{x}_{\mathrm{GW}}\gt \mathrm{thresh}}{\displaystyle \int }_{{x}_{\mathrm{EM}}\gt \mathrm{thresh}}\displaystyle \int p\\ & & \times \ ({x}_{\mathrm{GW}},{x}_{\mathrm{EM}},{D}_{L},{\rm{\Omega }},z| {H}_{0})\,{{dD}}_{L}\,{dz}\,d{\rm{\Omega }}\,{{dx}}_{\mathrm{GW}}\,{{dx}}_{\mathrm{EM}}\\ & = & \displaystyle \int {P}_{\det }^{\mathrm{GW}}({\hat{D}}_{L}(z,{H}_{0}),{\rm{\Omega }},z){P}_{\det }^{\mathrm{EM}}(z,{\rm{\Omega }}){p}_{0}(z,{\rm{\Omega }})d{\rm{\Omega }}{dz}\\ & = & {\displaystyle \int }_{0}^{{z}_{{\rm{h}}}}\displaystyle \int \int {P}_{\det }^{\mathrm{GW}}({\hat{D}}_{L}(z,{H}_{0}),{\rm{\Omega }},z){p}_{0}(z,{\rm{\Omega }})\,d{\rm{\Omega }}\,{dz},\end{array}\end{eqnarray} \tag{ 12 }$$

where we assume that only data that is above some threshold is detected, and we define:

$$\begin{eqnarray}&&{P}_{\det }^{\mathrm{GW}}({D}_{L},{\rm{\Omega }},z)\equiv {\int }_{{d}_{\mathrm{GW}}\gt \mathrm{thresh}}p({d}_{\mathrm{GW}}| {D}_{L},{\rm{\Omega }},z){{dd}}_{\mathrm{GW}},\end{eqnarray} \tag{ 13 }$$

and similarly:

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\det }^{\mathrm{EM}}(z,{\rm{\Omega }}) & \equiv & {\displaystyle \int }_{{d}_{\mathrm{EM}}\gt \mathrm{thresh}}p({d}_{\mathrm{EM}}| z,{\rm{\Omega }}){{dd}}_{\mathrm{EM}}\\ & = & { \mathcal H }({z}_{{\rm{h}}}-z),\end{array}\end{eqnarray} \tag{ 14 }$$

where ${ \mathcal H }$ is the Heaviside step function. We assume that the EM likelihood is constant with redshift up to a maximum (horizon) redshift, beyond which we assume there are no detectable host galaxies. In the statistical analysis in which the EM likelihood is assumed to be uninformative, z_h is equivalent to the maximum redshift of our galaxy catalog, or z_max defined before Equation (5). We calculate ${P}_{\det }^{\mathrm{GW}}$ by assuming a network signal-to-noise ratio threshold of 12 for detection, a monochromatic BNS mass distribution of 1.4–$1.4\ {M}_{\odot }$, zero spins, and isotropic inclination angles.

In practice, for nearby BNS sources, the term $\beta ({H}_{0})$ is insensitive to the precise details of this calculation or to the choice of ${z}_{{\rm{h}}}\gtrsim 0.2$, and is essentially $\beta ({H}_{0})\sim {H}_{0}^{3}$. This can be seen as follows. In LIGO-Virgo's second observing run, detectable BNS sources were within ∼100 Mpc (Abbott et al. 2018). For H₀ values within our prior range, this corresponds to redshifts $z\lesssim 0.07$, which is much smaller than the maximum detectable galaxy redshift, and so we can work in the limit ${z}_{{\rm{h}}}\to \infty$. We furthermore assume that the large-scale distribution of galaxies is uniform in comoving volume and we use the low-redshift, linear approximation ${H}_{0}={cz}/{D}_{L}$. At the low redshifts of detected BNS events, the redshifting of the GW signal in the detectors is negligible, and so we assume that ${P}_{\det }^{\mathrm{GW}}$ depends only on D_L and Ω, and is independent of z. With these approximations, we apply a different chain rule factorization to Equation (10) and write:

$$\begin{eqnarray}&&\begin{array}{l}p({x}_{\mathrm{GW}},{{dx}}_{\mathrm{EM}}| {H}_{0})\alpha ({H}_{0})=\displaystyle \int p({x}_{\mathrm{GW}}| {D}_{L},{\rm{\Omega }})\\ \times \ {p}_{0}({D}_{L},{\rm{\Omega }})p({x}_{\mathrm{EM}}| \hat{z}({D}_{L},{H}_{0}),{\rm{\Omega }})\,{{dD}}_{L}\,d{\rm{\Omega }},\end{array}\end{eqnarray} \tag{ 15 }$$

where $\alpha ({H}_{0})$ is a normalization term analogous to $\beta ({H}_{0})$. With this factorization, we can follow the steps in Equation (12) to write $\alpha ({H}_{0})$ as:

$$\begin{eqnarray}\begin{array}{rcl}\alpha ({H}_{0}) & = & \displaystyle \int {P}_{\det }^{\mathrm{GW}}({D}_{L},{\rm{\Omega }}){p}_{0}({D}_{L},{\rm{\Omega }})\\ & & \times \ {P}_{\det }^{\mathrm{EM}}(\hat{z}({D}_{L},{H}_{0}),{\rm{\Omega }})d{\rm{\Omega }}\,{{dD}}_{L},\end{array}\end{eqnarray} \tag{ 16 }$$

but this is now a constant (independent of H₀) because ${P}_{\det }^{\mathrm{EM}}(z,{\rm{\Omega }})=1$. We can then do a change of variables ${{dD}}_{L}=c/{H}_{0}\,{dz}$, and if we assume that ${p}_{0}({D}_{L},{\rm{\Omega }})\propto {D}_{L}^{2}$, we get:

$$\begin{eqnarray}\begin{array}{rcl}p({x}_{\mathrm{GW}},{x}_{\mathrm{EM}}| {H}_{0}) & \propto & \displaystyle \frac{1}{{H}_{0}^{3}}\displaystyle \int p({x}_{\mathrm{GW}}| {\hat{D}}_{L}(z,{H}_{0}),{\rm{\Omega }})\\ & & \times \ p({x}_{\mathrm{EM}}| z,{\rm{\Omega }}){p}_{0}(z,{\rm{\Omega }})\,d{\rm{\Omega }}\,{dz}.\end{array}\end{eqnarray} \tag{ 17 }$$

(Here we have dropped $\alpha ({H}_{0})$ because it is a constant.) This is equivalent to Equation (10) if we set $\beta ({H}_{0})\propto {H}_{0}^{3}$.

In order to explore whether the large-scale structure in the GW170817 localization volume, and the resulting statistical H₀ posterior, is typical for galaxies at $z\sim 0.01$, we rotate the true GW170817 sky map to different galaxies in the the MICE simulated catalog and repeat the statistical H₀ measurement. We assume that unlike for real galaxies in the GW170817 localization volume, no detailed observations have been carried out to measure the peculiar velocity field and apply group corrections. We therefore use the uncorrected redshifts given in the MICE catalog, which include a peculiar velocity contribution. The distribution of peculiar velocities is approximately described by a Gaussian of width 400 km s⁻¹, and we incorporate this uncertainty in the simulated H₀ measurements. Figure 7 shows the results for 20 realizations of the GW170817 localization volume centered on different galaxies in the MICE catalog. We see that the true GW170817 statistical H₀ measurement (we once again use the K-band luminosity-weighted ${L}_{K}\gt 0.1{L}_{K}^{\star }$ posterior shown in the right panel of Figure 6 as a representative posterior) is fairly typical among the different realizations. Over 50 different realizations, the information given by the difference in Shannon entropy between prior and posterior is ${0.43}_{-0.19}^{+0.43}$ bits (median and symmetric 90% intervals), whereas the information for the true GW170817 measurement is 0.67 bits. If we lower the peculiar velocity uncertainty in the simulations from 400 km s⁻¹ to 200 km s⁻¹, the GW1708170-like posteriors become slightly more informative on average, with a typical information gain of ${0.57}_{-0.27}^{+0.42}$ bits.

Zoom In Zoom Out Reset image size