A standard siren measurement of the Hubble constant from GW170817 without the electromagnetic counterpart

We perform a statistical standard siren analysis of GW170817. Our analysis does not utilize knowledge of NGC 4993 as the unique host galaxy of the optical counterpart to GW170817. Instead, we consider each galaxy within the GW170817 localization region as a potential host; combining the redshift from each galaxy with the distance estimate from GW170817 provides an estimate of the Hubble constant, $H_0$. We then combine the $H_0$ values from all the galaxies to provide a final measurement of $H_0$. We explore the dependence of our results on the thresholds by which galaxies are included in our sample, as well as the impact of weighting the galaxies by stellar mass and star-formation rate. Considering all galaxies brighter than $0.01 L^\star_B$ as equally likely to host a BNS merger, we find $H_0= 76^{+48}_{-23}$ km s$^{-1}$ Mpc$^{-1}$ (maximum a posteriori and 68.3% highest density posterior interval; assuming a flat $H_0$ prior in the range $\left[ 10, 220 \right]$ km s$^{-1}$ Mpc$^{-1}$). Restricting only to galaxies brighter than $0.626 L^\star_B$ tightens the measurement to $H_0= 77^{+37}_{-18}$ km s$^{-1}$ Mpc$^{-1}$. We show that weighting the host galaxies by stellar mass or star-formation rate provides entirely consistent results with potentially tighter constraints. While these statistical estimates are inferior to the value from the counterpart standard siren measurement utilizing NGC 4993 as the unique host, $H_0=76^{+19}_{-13}$ km s$^{-1}$ Mpc$^{-1}$ (determined from the same publicly available data), our analysis is a proof-of-principle demonstration of the statistical approach first proposed by Bernard Schutz over 30 years ago.


INTRODUCTION
The first multi-messenger detection of a binary neutron star (BNS) merger, GW170817, enabled the first standard siren measurement of the Hubble constant, H 0 , ushering in the era of gravitational-wave (GW) cosmology (Abbott et al. 2017a,c,d). This H 0 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 (Coulter et al. 2017;Soares-Santos et al. 2017;Abbott et al. 2017d). 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 be-tween 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 all 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 0 . The final H 0 measurement from a single event is the sum of all contributions from the individual galaxies. A Bayesian formalism for this method was first introduced by Del Pozzo (2012), who applied it to simulated GW data in order to predict the expected H 0 constraints from second-generation detections. An upto-date forecast incorporating realistic detection rates, galaxy peculiar velocities, large-scale structure, and ad-arXiv:1807.05667v1 [astro-ph.CO] 16 Jul 2018 ditional considerations can be found in .
We refer to this approach of measuring H 0 as the "statistical" method (Schutz 1986;MacLeod & Hogan 2008;Petiteau et al. 2011;, 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 0 measurement is uninformative .
The statistical approach may be the only way to do standard siren science with binary black holes, since 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;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(Abbott et al. , 2017b, although the higher rate is not expected to compensate for the inferior constraints . 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 0 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 0 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 §2, a discussion of the galaxy selection in §3, a discussion of the gravitational-wave constraints in §4, results in §5, and conclude in §6.

METHODS
We follow the statistical framework presented in  (see also Del Pozzo 2012; Gray et al. 2018). We include the role of GW selection effects, galaxy catalog incompleteness, galaxy luminosities, and redshift uncertainties in our analysis. The posterior on H 0 given the GW and EM data, x GW and x EM , is: whereD L (z, H 0 ) is the luminosity distance of a source at redshift z for a given H 0 (fixing other cosmological parameters to the Planck values; Ade et al. 2016) 1 , Ω is the sky position, and β(H 0 ) is a normalization term to ensure that the likelihood normalizes to 1 when integrated over all detectable GW and EM datasets (Mandel et al. 2016). The term p 0 (H 0 ) represents the prior on the Hubble constant. A detailed derivation of Equation 1, including the role of the normalization term β(H 0 ), is provided in the Appendix. 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 GW |D L , Ω), 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. The corresponding posterior p(D L , Ω | x GW ) ∝ p(x GW |D L , Ω)p 0 (D L , Ω) is summarized in the 3-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,b). For this analysis, we use the publicly released 3-dimensional sky map from Abbott et al. (2018a) (see §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 , Ω) ∝ D 2 L . Meanwhile, the EM likelihood term p(x EM | z, Ω) 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 EM |z, Ω) is completely uninformative, and set: p(x EM |z, Ω) ∝ 1. (2) In this case, the redshift information enters only through the prior term, p 0 (z, Ω), which we take to be a galaxy catalog. The detection of an electromagnetic counterpart typically results in p(x EM | z, Ω) being strongly peaked around someΩ 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 ). This reduces the number of potential host galaxies, and thereby improves the measurement. The galaxy catalog term p 0 (z, Ω) is given by: where p 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: wherez i ,Ω i denotes the (peculiar velocity-corrected) "Hubble" redshifts and sky coordinates of all galaxies in the catalog, and N (µ, σ; 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. 2017). 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 (Ω, σ Ω ; Ω) by δ(Ω − Ω). The weights w i can be chosen to reflect the a priori belief that a galaxy could host a gravitational-wave source. For example, setting all weights to w i = 1 N gal corresponds to equal probability for each galaxy to host a gravitational wave source. In general, since 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 ∝ L i B or w i ∝ L i K , and explore the dependence of the result on these of weightings.
To calculate the term p miss in Equation 3, we assume that on large scales, the distribution of galaxies, p 0 (z, Ω) is uniform in comoving volume. Let p vol (z, Ω) denote the cosmologically homogeneous and isotropic distribution normalized over the volume contained in the range z min < z < z max considered in our analysis. (The result does not depend on z max provided it is large enough to encompass all possible redshifts of the source for all allowed values of H 0 .) Assuming all galaxies are weighted equally, the distribution of missing galaxies is written as: where P complete (z) is the probability that a galaxy at redshift z is in the catalog, and the completeness fraction f is given by: 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, 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: and, weighting each missing galaxy by its luminosity, Equation 5 becomes: 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.) In the next section, the completeness function P complete is estimated for the galaxy catalog used in the analysis.
Completeness of the GLADE catalog as a function of redshift for galaxies brighter than 0.25L B (solid blue curve), 0.05L B (dashed green curve), and 0.01L B (dotdahsed orange curve), calculated by comparing the redshift distribution of galaxies in GLADE to a distribution that is constant in comoving volume. For galaxies brighter than 0.626L B , GLADE is complete across the entire redshift range shown.

GALAXY CATALOGS
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 which are also listed in the group catalog of Kourkchi & Tully (2017), as identified by a common PGC 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. 2017).
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 φ = 5.5 × 10 −3 h 3 0.7 Mpc 3 , α B = −1.07, L B = 2.45 × 10 10 h 2 0.7 L B, throughout. The corresponding characteristic absolute magnitude is M B = −20.47 + 5 log 10 h 0.7 . We will also consider the K-band magnitudes reported in GLADE when applying galaxy weights, and we use the Kband Schechter function parameters of α K = −1.02, M K = −23.55 + 5 log 10 h 0.7 (Lu et al. 2016). Figure 1 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 ∼ 0.626h −2 0.7 L B , corresponding to about 0.66 of the Milky Way luminosity for h 0.7 = 1. Galaxies brighter than ∼ 0.626L B make up half of the total luminosity for the given Schechter function parameters. We find that for z 0.03, GLADE is complete for galaxies down to 2.5 times dimmer, or ∼ 0.25L B , corresponding to M B = −18.96 + 5 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 ∼ 0.05L B (M B = −17.22 + 5 log 10 h 0.7 ), GLADE is ∼ 70% complete at z ∼ 0.03, and even if we consider galaxies down to ∼ 0.01L B (M B = −15.47+5 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 Two-dimensional localization region of GW170817 (blue contours) with the sky coordinates of the 408 GLADE galaxies (green crosses) within the 99% localization area and the redshift range 0 < z 0.046 (for an H0 prior range of H0 ∈ [10, 220] km s −1 Mpc −1 ). The light and dark blue contours enclose the 50% and 90% probability regions, respectively, and the shading of the galaxy markers denotes their redshifts, corrected for peculiar and virial motions as described in the text.  Figure 3. Probability distribution of the redshifts of potential hosts to GW170817 weighted by the GW sky map probability, p(z) = p(xGW | Ω)p0(z, Ω)dΩ, compared to a uniform in comoving volume distribution of galaxies, p vol (z). For the orange histogram, we include all galaxies in the catalog brighter than 0.626L B . For galaxies brighter than 0.626L B , the catalog is complete over the redshift range. However, when we lower the luminosity cutoff to 0.25L B (yellow histogram) or 0.005L K (green and blue), we must account for catalog incompleteness at higher redshifts by considering the redshift and luminosity distributions of the missing galaxies (see §2). The yellow (green) histogram additionally weights each galaxy by its B-band (K-band) luminosity. If the ratio p(z)/ p vol (z) were completely flat, we would expect an uninformative H0 measurement in which our posterior recovers our prior. However, in all instances there is a dominant peak at z ∼ 0.01, suggesting that the resulting H0 measurement will be informative. Adding in luminosity weights, especially in the K-band, makes the peak more dominant.
K-band Schechter function parameters, GLADE is complete up to z ∼ 0.045 for galaxies with L K > 0.36L K , which contain 70% of the total K-band luminosity, and up to z ∼ 0.03 for galaxies with L K > 0.1L K , which contain 90% of the total luminosity. For galaxies brighter than L K = 0.005L K , which make up more than 99% of the total K-band luminosity, GLADE is ∼ 70% complete at z = 0.01.

SOURCE LOCALIZATION AND DISTANCE
From the GW data alone, GW170817 is the bestlocalized GW event to date. The original analysis by the LIGO-Virgo collaboration (Abbott et al. 2017c) reported a 90% localization area of 28 deg 2 and a 90% localization volume of 380 Mpc 3 (assuming Planck cosmology; Ade et al. 2016), while the most recent analysis (Abbott et al. 2018a) improves the 90% localization area to 16 deg 2 and the 90% volume to 215 Mpc 3 . We use this updated 3-dimensional sky map (Singer et al. 2016a,b) from Abbott et al. (2018a) throughout. 2 Figure 2 shows the 2-dimensional sky map together with the galaxies in the GLADE catalog within the localization region. Figure 3 shows that, although there are a total of 408 galaxies within the 99% localization area (see Figure 2), 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.
The 3-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 Equal galaxy weights . Posterior probability of H0 under various assumptions regarding the potential host galaxy. We adopt a flat H0 prior in the range H0 ∈ [10, 220] km/s/Mpc. For the dashed orange curve, we assume that only galaxies brighter than 0.626L B can host BNS events, meaning that the galaxy catalog is complete over the relevant redshift range. The solid green curve lowers the luminosity cutoff to 0.25L B , and accounts for the mild incompleteness of the catalog above redshifts z ∼ 0.03. The dotted blue curves incorporate all galaxies brighter than 0.01L , accounting for the incompleteness of faint galaxies at redshifts z 0.01. The dot-dashed pink curve shows the H0 measurement assuming the host galaxy is known to be NGC 4993.
calibrated by general relativity (Schutz 1986). No distance ladder is required.

RESULTS
We combine the GW distance posterior for GW170817 with the redshift for each potential host galaxy within the localization region. As detailed in §2, each galaxy produces a posterior probability for H 0 , and we combine these estimates among all the galaxies in the localization region to arrive at a final estimate for H 0 . We adopt a flat prior in H 0 over the range 10-220 km/s/Mpc. The results are presented in Figure 4. Because the galaxies are predominantly found in one galaxy group at z ∼ 0.01, the H 0 posterior shows a clear peak at H 0 = 76 km/s/Mpc. And because NGC 4993, the true galaxy host of GW170817, is a member of the group at z ∼ 0.01, we should not be surprised to learn that the peak in H 0 is consistent with the H 0 estimate from the GW170817 standard siren measurement including the counterpart (Abbott et al. 2017a). Because this analysis has been performed on the public 3-dimensional sky map, rather than proprietary LIGO/Virgo data, the results do not agree precisely with those of Abbott et al. (2017a), and in particular, the position of the peak in Figure 4 is at H 0 = 76 km/s/Mpc instead of H 0 = 70 km/s/Mpc. This is because the 3-dimensional sky map approximates the distance posterior along each line-of-sight by a simple 2-parameter Gaussian fit (see Eq. 1 of Singer et al. 2016a), which is an imperfect approximation to the true, asymmetric distance posterior Del Pozzo et al. 2018). On the other hand, the analysis in Abbott et al. (2017a) utilizes the full distance posteriors along each line-of-sight. Figure 4 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. This peak is slightly shifted compared to the result presented in Abbott et al. (2017a) due to a different prior choice (flat in H 0 compared to flat in log H 0 ) and the usage of the Gaussian fit to the distance posterior found in the 3-dimensional sky map as discussed above.
The other curves in Figure 4 assume different limiting thresholds for what constitutes a potential host galaxy. As the threshold is lowered, additional galaxies fall into the sample, and the H 0 posterior is broadened. For a limiting B-band magnitude of 0.25L B , we need to account for the incompleteness of the galaxy catalog at redshifts z 0.03, and for 0.01L * B , we need to account for the incompleteness at z 0.01, as described in §2. The incompleteness correction leads to a slight additional broadening of the H 0 posterior, but the clear peak at H 0 = 76 km/s/Mpc remains. This peak is the result of the galaxy group at z ∼ 0.01, of which NGC 4993 is a member. The curves in Figure 5 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 ∼ 76 km/s/Mpc becomes more pronounced when galaxies are weighted by their luminosity, since the group containing NGC 4993 consists of many bright, mostly red galaxies. 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 0 analysis.

CONCLUSION
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 years 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 0 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 = 76 km s −1 Mpc −1 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.
Although this statistical standard siren measurement of H 0 is less precise than the counterpart measurement of H 0 = 76 +19 −13 km s −1 Mpc −1 (for a flat prior and utilizing the distance ansatz in the 3-dimensional sky map; see §5), it nonetheless shows that interesting constraints on cosmological parameters are possible from gravitational-wave sources even in the absence of an optical counterpart and an identification of the unique host galaxy Gray et al. 2018). Although detailed studies find that the measurement of cosmological parameters from the counterpart approach is likely to surpass the statistical approach , the statistical approach offers an important crossvalidation 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 binary neutron star sys- H 0 (km s −1 Mpc −1 ) K-band luminosity weights Figure 5. Posterior probability of H0, weighting all galaxies in the volume by their B-band luminosities, corresponding roughly to weighting by star-formation rate (left), or K-band luminosities, corresponding roughly to weighting by stellar mass (right).
We have applied the necessary completeness correction (see §2). The blue dashed-dot curve shows all galaxies brighter than 0.01L B in B-band (left) or 0.005L K in K-band (right) with equal weights for comparison. Weighting galaxies by their K-band luminosities brings all the curves into very close agreement, because many galaxies in the group at z ∼ 0.01 have brighter than average K-band luminosities (brighter than 1.5L K ) and thus dominate the K-band weighted population and contain the majority of the stellar mass.
tems and are expected to lack electromagnetic counterparts. The inferior quality of the individual H 0 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 0 may con-strain additional cosmological parameters such as the equation of state of the dark energy.
We thank Alison Coil and Risa Wechsler for valuable discussions about galaxy properties.

APPENDIX
In this appendix we derive the H 0 posterior probability distribution function from Equation 1. We write the likelihood for GW and EM data, x GW and x EM , given a value of H 0 as: and factor the numerator as: The H 0 posterior is related to the likelihood in Equation 2 by a prior: This equation is identical to Equation 1 in the main text. The normalization term β(H 0 ) is given by (see Mandel et al. 2016;: where we assume that only data that is above some threshold is detected, and we define: and similarly: where 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 GW det by assuming a network signal-to-noise ratio threshold of 12 for detection, a monochromatic BNS mass distribution of 1.4-1.4 M , zero spins, and isotropic inclination angles.
In practice, for nearby BNS sources, the term β(H 0 ) is insensitive to the precise details of this calculation or to the choice of z h 0.2, and is essentially β(H 0 ) ∼ H 3 0 . This can be seen as follows. In LIGO-Virgo's second observing run, detectable BNS sources were within ∼ 100 Mpc (Abbott et al. 2018b). For H 0 values within our prior range, this corresponds to redshifts z 0.07, which is much smaller than the maximum detectable galaxy redshift, and so we can work in the limit z h → ∞. 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 GW det depends only on D L and Ω, and is independent of z. With these approximations, we apply a different chain rule factorization to Equation 2 and write: p(x GW , dx EM |H 0 )α(H 0 ) = p(x GW |D L , Ω)p 0 (D L , Ω)p(x EM |ẑ(D L , H 0 ), Ω) dD L dΩ, where α(H 0 ) is a normalization term analogous to β(H 0 ). With this factorization, we can follow the steps in Equation 4 to write α(H 0 ) as: α(H 0 ) = P GW det (D L , Ω)p 0 (D L , Ω)P EM det (ẑ(D L , H 0 ), Ω) dΩ dD L , but this is now a constant (independent of H 0 ) because P EM det (z, Ω) = 1. We can then do a change of variables dD L = c/H 0 dz, and if we assume that p 0 (D L , Ω) ∝ D 2 L , we get: H 0 ), Ω)p(x EM |z, Ω)p 0 (z, Ω) dΩ dz.
(Here we have dropped α(H 0 ) because it is a constant.) This is equivalent to Equation 2 if we set β(H 0 ) ∝ H 3 0 .