key: cord-0167713-btazcjl6
authors: Vasylyev, Sergiy; Filippenko, Alex
title: A Measurement of the Hubble Constant using Gravitational Waves from the Neutron-Star Black-Hole Candidate GW190814
date: 2020-07-22
journal: nan
DOI: nan
sha: 5324e54dd951a3d3b41892d8919ebd3307f760b3
doc_id: 167713
cord_uid: btazcjl6

We present a test of the statistical method introduced by Bernard F. Shutz in 1986 using only gravitational waves to infer the Hubble constant (H$_0$) from GW190814, the first first high-probability neutron-star--black-hole (NS-BH) merger candidate detected by the Laser Interferometer Gravitational Wave Observatory (LIGO) and the Virgo interferometer. We apply a baseline test of this method to the binary neutron star (BNS) merger GW170817 and find H$_0 = 70^{+35.0}_{-18.0}$ km s$^{-1}$ Mpc$^{-1}$ (maximum a posteriori and 68.3% highest density posterior interval) for a galaxy $B$-band luminosity threshold of $L_B geq 0.001 L_B^*$ with a correction for catalog incompleteness. Repeating the calculation for GW190814, we obtain H$_0 = 67^{+41.0}_{-26.0}$ km s$^{-1}$ Mpc$^{-1}$ and H$_0 = 71^{+34.0}_{-30.0}$ km s$^{-1}$ Mpc$^{-1}$ for $L_B geq 0.001 L_B^*$ and $L_B geq 0.626 L_B^*$, respectively. Combining the posteriors for both events yields H$_0 = 70^{+29.0}_{-18.0}$ km s$^{-1}$ Mpc$^{-1}$, demonstrating the improvement on constraints when using multiple gravitational wave events. We also confirm the results of other works that adopt this method, showing that increasing the $L_B$ threshold enhances the posterior structure and slightly shifts the distribution's peak to higher H$_0$ values. We repeat the joint inference using the low-spin PhenomPNRT and combined (SEOBNRv4PHM + IMRPhenomPv3HM) posterior samples for GW170817 and GW190814, respectively, achieving a tighter constraint of H$_0 = 69^{+29.0}_{-14.0}$ km s$^{-1}$ Mpc$^{-1}$.

Gravitational wave (GW) and electromagnetic (EM) follow-up observations of black hole (BH) and neutron star (NS) mergers provide a novel method of probing dense astrophysical environments and enable a unique channel for cosmology. The August 2017 discovery of NS-NS merger GW170817 by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and the Virgo interferometer, as well as the subsequent multimessenger observations by the coordination of thousands of astronomers, has brought forth a new era of astronomy (Abbott et al. 2017a) . Two weeks after this discovery, the LIGO detectors in Hanford, Washington and Livingston, Louisiana, as well as the newly added Virgo detector in Italy, ceased operation for upgrades, marking the end of Observing Run 2 (O2). The analysis of GW170817 covered the entire EM spectrum, giving insights into the generation of short-duration gamma-ray bursts (sGRBs; Goldstein et al. 2017) , general relativity in the strong-gravity regime (Abbott et al. 2019a) , the production of r-process elements (Kasen et al. 2017) , constraints on the NS equation of state (Annala et al. 2018) , and an independent measurement of the Hubble constant H 0 (Abbott et al. 2017b ).

Observing Run 3 (O3) began in late April of 2019, promising a higher rate of GW detections owing to sensitivity upgrades to the three detectors. The search volume increased by ∼ 100% from O2 to O3 according to Table 2 in KAGRA Collaboration, LIGO Scientific Collaboration, and Virgo Collaboration et al. (2018) , thereby allowing for a significantly deeper survey. The end of O3 was initially planned for April 30, 2020 in order to conduct additional upgrades, but O3 was ended prematurely on March 27, 2020 because of the COVID-19 pandemic. In summary, O3 produced 36 BH-BH mergers, 1 NS-NS merger, and the first-ever candidate NS-BH merger (GW190814) with at least 90% confidence. In this paper, we will focus on the candidate NS-BH merger GW190814, as its exquisite localization (see Section 2.3) makes it a prime object for the statistical inference of H 0 described below. Hereinafter, we will refer to GW190814 as a NS-BH merger, though we acknowledge that this event has not been explicitly shown to be a NS-BH; it likely contains compact objects with a mass range falling within the criterion for a NS-BH merger. 1

Recent measurements of H 0 reveal at least a 4.2σ tension. The Planck satellite team infers H 0 = 67.4 ± 0.5 km s −1 Mpc −1 using cosmic microwave background (CMB) data and assuming the standard ΛCDM model is correct (Planck Collaboration et al. 2018) . The SH0ES (Supernovae, H 0 , for the Equation of State of Dark Energy) team measured H 0 = 74.03 ± 1.42 km s −1 Mpc −1 with an independent method, combining Cepheid-variable calibrations with luminosity distance measurements of Type Ia supernovae (Riess et al. 2018 . The SH0ES method is sensitive to how well the distance ladder is calibrated. Other recent studies have only enhanced the discrepancy between the two methods mentioned above, for a combined tension of ∼ 6σ (Riess 2020) . For example, the H0LiCOW (H 0 Lenses in COSMOGRAILs Wellspring) and STRIDES (STRonglensing Insights into Dark Energy Survey) teams use measured time delays in the light curves of different images of single strongly lensed quasars to obtain H 0 = 73.3 ± 1.8 km s −1 Mpc −1 and H 0 = 74.2 ± 1.4 km s −1 Mpc −1 , respectively, agreeing with the SH0ES measurement (Wong et al. 2019; Shajib et al. 2020) . Other measurements that broadly agree with the SH0ES value include the tip of the red giant branch (Jang & Lee 2017; Hatt et al. 2018, TRGB1, TRGB2; ) , Mira variables (Huang et al. 2020) , surface brightness fluctuations (SBF), and masers (Verde et al. 2019 ). On the other hand, measurements from big-bang nucleosynthesis and baryon acoustic oscillations (BBN + BAO; Cuceu et al. 2019) , Wilkinson Microwave Anisotropy Probe (WMAP) CMB + BAO (Hinshaw et al. 2009 ), Atacama Cosmology Telescope Polarization camera (ACTPol) + BAO (Louis et al. 2017) , and the South Pole Telescope SZ camera (SPT-SZ) + BAO (Story et al. 2013 ) favor the Planck result.

The above tension between local (SH0ES, etc.) and earlyuniverse (Planck) measurements of H 0 may be the result of underlying systematic effects or astrophysical causes. A new independent method using gravitational-wave sources as "standard sirens," first proposed by Schutz (1986) , could reconcile this discrepancy. This method uses the observed amplitude and the frequency of the gravitational waveform (from which a distance is determined) together with the measured redshift of the source's host galaxy to infer H 0 . Standard sirens do not need a distance ladder and are thereby decoupled from systematic effects that may be introduced in methods relying on it.

The true host galaxy can only be identified with an EM counterpart. GW170817 produced an optical transient powered by radioactive decay, called a kilonova (Kasen et al. 2017) . The source's host galaxy was identified as NGC 4993 at a luminosity distance of 40 +8.0 −14.0 Mpc (Coulter et al. 2017; Soares-Santos et al. 2017; Valenti et al. 2017; Arcavi et al. 2017; Tanvir et al. 2017; Lipunov et al. 2017) . Abbott et al. (2017b) presented the first result using the standard siren method, estimating H 0 = 70.0 +12.0 −8.0 km s −1 Mpc −1 .

BH-BH mergers are not expected to produce an optical transient, unlike NS-NS mergers. Some theoretical models suggest that NS-BH mergers may only produce a counterpart under certain physical conditions (Foucart et al. 2018 (Foucart et al. , 2019 Barbieri et al. 2020) . However, cosmological inferences are still possible without the EM counterpart.

The LIGO-Virgo detectors produce a skymap that constrains the location of a GW source to specific patches in the sky. Using a statistical approach (the catalog method), one can consider every galaxy in the source's localization region as a potential host with some probability. Summing the probability assigned to each galaxy builds the full posterior on H 0 . Chen et al. (2018) show that an H 0 inference from GW sources with optical counterparts will converge faster than for dark sirens. Although using a so-called "dark siren" will not provide as precise a measurement from a single event compared to the case when the host galaxy is known, many more BH-BH mergers (dark sirens) are expected than NS-NS mergers (Baibhav et al. 2019) , providing a useful validation test of the optical-counterpart method. Furthermore, combining several measurements will yield increasingly tighter constraints on H 0 (Chen et al. 2018; Nair et al. 2018) .

The galaxy catalog (statistical) method was first tested on simulated data by Del Pozzo (2012). More recently, Fishbach et al. (2019) used this method to infer H 0 from GW170817 without relying on the EM counterpart. They obtain several estimates for H 0 using various luminosity cuts and weighting schemes to galaxies in the GLADE 2.3 catalog described in Section 2.2. An estimate of H 0 from BH-BH merger GW170814 using a similar statistical method was recently obtained by the Dark Energy Survey ( (Abbott et al. 2019b ). This work also explored the effects of galaxy luminosity weighting on the H 0 posterior shape.

On 2019-08-14, at 21:10:39 UT, LIGO Hanford, LIGO Livingston, and Virgo detected the GW event GW190814 with a false alarm rate (FAR) of approximately 1 per 10 25 years at a luminosity distance of 267 ± 52 Mpc (LIGO Scientific Collaboration & Virgo Collaboration 2019). Although the NS-BH candidate GW190814 did not have an observed EM counterpart, we can still use the gravitational-wave data to produce meaningful results. Analysis by Abbott et al. (2020) showed that the primary and secondary masses of GW190814 are 23 +1.1 −1.0 M and 2.59 +0.08 −0.09 M , respectively. The secondary mass approaches the observational MassGap (3-5 M ), in which there is uncertainty whether the object is either the heaviest neutron star or the lightest black hole ever discovered. This should not affect our results, given our generous prior on the NS mass for this event.

We follow the methodology presented by Chen et al. (2018) , Gray et al. (2019) , and Abbott et al. (2019b) to obtain the H 0 posterior. With this paper we test the statistical method on a new GW-type candidate (NS-BH) and improve the accessibility of the gwcosmo code. See the Appendix for a detailed discussion of the mathematics involved.

Using the publicly available gwcosmo 2 code, we construct a posterior on the Hubble constant using only gravitational waves for the NS-BH merger candidate GW190814. We use the Bayesian framework presented by Chen et al. (2018) and Gray et al. (2019) , which is detailed in the Appendix. We outline our methodology starting with a thorough account of our assumed priors and input parameters used in gwcosmo. Note that we adopt the 02-H0 branch to perform these calculations, as it is the most stable at the time of writing. The preparation and injection of the GLADE 2.0 galaxy catalog is discussed in Section 2.2. The HEALPIX localization skymaps used for all of the calculations are described in Section 2.3. We create a baseline test of our assumptions by comparing to Abbott et al. (2019b) and Fishbach et al. (2019) using GW170817. We then explore parameter space to present multiple H 0 realizations for GW190814.

Our analysis is carried out with both a uniform and "flat log prior" on H 0 over a set of different intervals.We use the definition for the flat log prior p(H 0 ) ∝ H −1 0 (Abbott et al. 2019b ). Below, we describe the options passed to the gwcosmo single posterior script in the gwcosmo code.The italicized items are presented in Table 1 . The uniform component follows p(m 2 ) = constant, while the power-law component takes the form p(m 1 ) ∝ m −α 1 , with the power-law index α = 1.6. Here, m 2 is the secondary (NS) mass and m 1 is the primary (BH) mass.

2. The power spectral density (PSD) parameter is associated with the detector sensitivity during either the O1, O2, or O3 observing runs (we choose O2 for GW170817 and O3 for GW190814).

3. The completeness parameter is defined as the ratio of the number of galaxies in a chosen galaxy catalog to the true number of galaxies in the cosmological volume. We discuss this in more detail in Section 2.2. and in the Appendix. Gray et al. (2019) study the effects of this parameter on the H 0 posterior extensively on simulated merger data in Sections III and IV.

4. Galaxy weighting may be set to either "False" (equal weights) or have B-band luminosity-dependent weights ω i ∝ L i B . A luminosity-dependent weighting scheme follows the assumption that BH and NS merger rates scale with star-formation rates (Fong & Berger 2013) .

5. Luminosity threshold gives the minimum B-band luminosity considered for the calculation; this parameter will be explored in depth in Section 3.

We hold the following parameters constant throughout every calculation.

1. Linear cosmology is set to "False" because we include galaxies with redshift z > 0.1.

3. Basic pdet is set to "False" allowing us to take into account redshifted mass, M z = M (1 + z) ).

4. The uncertainty parameter is set to "True," taking into account the Gaussian uncertainties in redshift for each galaxy (see Section 2.2).

5. The rate evolution parameter is set to "False," which describes a constant merger rate R(z) as it appears in Equation 11 of Abbott et al. (2019b) . We assume the ΛCDM model (Ω m = 0.308, Ω Λ = 0.692). The effect of choosing a merger rate with a dependence on redshift is shown extensively by Abbott et al. (2019b) .

We use the GLADE v2.3 galaxy catalog throughout, adopting the parameters φ * = 1.6 × 10 −2 h 3 Mpc −3 , where h = 0.7 and β = −1.07 for the Schechter B-band luminosity function,

where ρ(x) is the number density of galaxies and the characteristic B-band luminosity L * B corresponds to an absolute magnitude M B = −20.47 (Gehrels et al. 2016) .

The GLADE galaxy catalog contains nearly 3 million galaxies and is complete up to 300 Mpc at L B = 0.626 L * B , corresponding to the median of the luminosity function (Arcavi et al. 2017) . The statistical method is sensitive to the galaxy completeness fraction, f . The Glade catalog's high completeness fraction over the redshifts considered for this study provides a significant advantage over other catalogs (Dlya et al. 2018) . Approximately half of the objects in the catalog have a measured B-band luminosity.

All redshifts in GLADE are corrected for peculiar motions and are in the heliocentric frame (Carrick et al. 2015) . In order to correct for radial group velocities, we crossreference GLADE galaxies by their "PGC ID" with the Principal Galaxy Catalog (PGC) and use the corresponding corrected radial velocities in the heliocentric frame (Kourkchi & Tully 2017) . We then correct these heliocentric velocities to the CMB frame using NASA/IPAC Extragalactic Database (NED) with parameters l apex = 264.14 • , b apex = +48.26 • , and v apex = 371.0 km s −1 (Fixsen et al. 1996) . We define the z rad and z CMB parameters as the radial group velocity and CMB reference frame corrections, respectively. Finally, we assign a 200 km s −1 Gaussian uncertainty to the velocity (cz) for each galaxy.

For our purposes, we only extract the right ascension (degrees), declination (degrees), redshift, apparent B magnitude, and absolute B magnitude (RA, Dec, z, B, B abs , respectively) from the raw catalog and build the pickle 3formatted dictionary that gwcosmo requires. Note that gwcosmo expects the RA and Dec to be in radians. We use the provided B-band magnitude in our galaxy-weighting procedures.

During O2 and O3, LIGO-Virgo released public alerts accompanied by an allsky HEALPix localization skymap in a FITS file format for each GW event. The skymap FITS includes the sky position, probability, and distance for each pixel. The original 90% region for GW170817 (28 deg 2 ) was improved to 16 deg 2 using the LALInference pipeline (Veitch et al. 2015) . For our analysis, we adopt the updated skymap. 4 For GW190814, we use the updated skymap from the GraceDB database, localizing the source to 23 deg 2 and 5 deg 2 in the 90% and 50% confidence regions, respectively 5 . The marginalized distance posterior found in the skymap is a symmetric Gaussian fit to the full, potentially asymmetric posterior sample. According to Fishbach et al. (2019) , this assumption has the effect of moving the peak of the H 0 pos- terior by as much as 9% compared to using a full posterior sample. For our calculations, Equation A5 of the Appendix only includes galaxies in the 99.9% region (assigning weight = 1 to each if equal weights) and assigns a weight of 0 for those outside of the localization region.

We split our analysis into two parts. First, we summarize our results for GW170817 and compare our H 0 inference to previous works using this statistical method Fishbach et al. 2019; Abbott et al. 2019b) . We also explore possible systematic differences and assumptions that may carry into our H 0 calculation for GW190814. We then repeat the procedure over the parameter space introduced in Section 2 for GW190814. All H 0 measurements assume the 68.3% highest density posterior interval. We choose the H 0 prior range as [30, 200] or [40, 140] km s −1 Mpc −1 for ease of comparison to other works. Finally, we repeat the calculation using the full posterior samples released by the LVC on June 25, 2020.

We note that after submission and initial review of this paper, the Dark Energy Survey (DES) Collaboration posted a paper (Palmese et al. 2020) in which a similar analysis of GW190814 is presented, but they did not include posterior samples.

Since the GLADE catalog is 100% complete up to z = 0.03 for galaxies that are above 0.25 L * B , and we consider galaxies well above this redshift threshold and down to 0.001 L * B (M B = −12.96 mag), we set the completeness parameter to "False" in our final calculation unless stated otherwise. We consider a flat log and uniform prior on H 0 on the intervals [30, 200] and [40, 140] km s −1 Mpc −1 . We use the BNS-Gaussian mass distribution centered on µ = 1.35 M and a standard deviation of σ = 0.15 M . Galaxies up to z ≈ 0.5 (z max = 0.5) are allowed, with equal luminosity weights assigned. A 200 km s −1 Gaussian uncertainty is applied to the velocity (cz) of each galaxy. 2662 galaxies fall into the 99% skymap localization corresponding to 34 deg 2 .

We employ a constant rate evolution term discussed in Section 2.1. The resulting H 0 posterior is illustrated in Figure  2 .

Given a uniform prior along the intervals [30, 200] and [40, 140] km s −1 Mpc −1 , we infer a Hubble constant 70 +50.0 −23.0 km s −1 Mpc −1 and 70 +35.0 −18.0 km s −1 Mpc −1 , respectively. Even when accounting for a luminosity cut ≥ 0.626 L * B yielding H 0 = 73 +36.0 −17.0 , our peak is less pronounced and is shifted compared to the H 0 ≥ 74 km s −1 Mpc −1 obtained by Fishbach et al. (2019) . This discrepancy may be caused by differences in both our chosen subset of the GLADE galaxy catalog and our velocity corrections. As a qualitative test, we also calculate H 0 without the radial group velocity and CMB reference frame corrections, yielding 67 +37.0 −19.0 km s −1 Mpc −1 for a uniform prior interval, [40, 140] km s −1 Mpc −1 . Given the sensitivity of the posterior to the injected catalog, we expect a slight deviation if the catalog is handled differently. These systematic differences will carry over into our calculation for the NS-BH merger.

For GW190814, we assume the same B-band Schechter parameters chosen for the GW170817 calculations. Here, we focus only on the [40, 140] km s −1 Mpc −1 prior range on H 0 . We show our results in Table 1 and illustrate the H 0 posterior in Figure 3 , where we plot five realizations for different luminosity considerations. We apply a max-imum redshift limit z max = 0.5 and assume the NSBHuniform mass distribution. 64,735 galaxies are obtained in the 99% localization region. We repeat the calculation for GW190814 to obtain H 0 = 67 +41.0 −26.0 km s −1 Mpc −1 and H 0 = 71 +34.0 −30.0 km s −1 Mpc −1 for L B ≥ 0.001 L * B and L B ≥ 0.626 L * B , respectively. Our tested parameter space is detailed in Table 1 . According to Figure 3 , the posterior peak is more pronounced for stricter luminosity cuts. The peak of our H 0 posterior shifts by ∼ 6% over the range of luminosity cuts, which is in agreement with Fishbach et al. (2019) .

The GW190814 posterior appears significantly flatter compared to GW170817. Given that the GW190814 localization covered a larger volume and included more than ten times as many galaxies as GW170817, we expect the GW190814 posterior to be washed out. The bump in the H 0 = 100-140 km s −1 Mpc −1 range is enhanced when galaxy weighting is set to "True"; it may be an artifact of enhanced GLADE catalog features similar to Figure 2 of Abbott et al. (2019b) .

We then combine the posteriors for GW170817 with S190814 using the gwcosmo combined posterior script, yielding H 0 = 70 +29.0 −18.0 km s −1 Mpc −1 as shown in Figure  4 . Although the peak is centered between the Planck and SH0ES results, we caution that the value is subject to systematics arising from luminosity cuts or weighting. For example, our combined posterior used a luminosity threshold of 0.626 L * B for GW190814, whereas taking the luminosity down to 0.001 L B would produce a peak at H 0 ≤ 70 km s −1 Mpc −1 . The combined posterior demonstrates the ability to further constrain the Hubble constant with multiple GW sources.

Systematic biases in the joint posterior due to varying population parameters of astrophysical sources are expected to be smaller than the statistical uncertainties due to contributions from the galaxy catalog given a high probability that the host galaxy is in the catalog. GW190814 has a median source redshift of z event = 0.053 using the combined waveform model from Abbott et al. (2020) , corresponding to p(G|z event , D w ) > 0.6 (or "high in-catalog probability") with the GLADE catalog according to Figure 1 of Abbott et al. (2019b) . Therefore, we take the contributions from the galaxy catalog to be the dominant source of uncertainties for GW190814.

In light of the LVC data release for GW190814 on June 25, 2020, we apply the full posterior sample to our H 0 calculation (Abbott et al. 2020) . We now account for biases introduced when only using a Gaussian approximation to the distance posterior from the 3D skymap. For GW170817, we use the low-spin PhenomPNRT posterior sample 6 . For GW190814, we use a combined posterior sample consisting of the SEOB-NRv4PHM (EOBNR PHM; Babak et al. 2017; Ossokine et al. 2020) and IMRPhenomPv3HM (Phenom PHM; Khan et al. 2019 Khan et al. , 2020 Waveform Models. 7 We compare the differences in the GW190814 H 0 posterior with and without the use of posterior samples in Figure 5 .

Following the reasoning of Fishbach et al. (2019) , using a full posterior sample (accounting for masses and spins) as opposed to a Gaussian approximation via the 3D skymap can have the effect of shifting the H 0 posterior peak. In our case, we observe a notable shift for both sets of parameters, favoring a higher value for the Hubble constant when using posterior samples. In Figure 6 , we show the combined H 0 posterior following the same procedure used to produce Figure Table 1 . We illustrate the difference between using posterior samples (a combined SEOBNRv4PHM and IMRPhenomPv3HM waveform model) and using a Gaussian approximation to the distance posterior from the 3D skymap. The C (solid green) and E (solid pink) curves are identical to those in Figure 4 . We label C* (dashed green) and E* (dashed purple) to indicate when posterior samples are used. We include the most recent Planck (orange vertical line) and SH0ES (gray vertical line) measurements and their corresponding 1σ uncertainties (shaded regions). Table 1 together with a SEOBNRv4PHM + IMRPhenomPv3HM posterior sample (labeled E* in Figure 5 ). We include the most recent Planck (orange vertical line) and SH0ES (gray vertical line) measurements and their corresponding 1σ uncertainties (shaded regions). ics were identified, including the injected catalog, luminosity weighting, and luminosity thresholds.

The motivation for this method follows the expectation of having significantly more GW-only mergers (dark sirens) than mergers with optical counterparts, providing a valuable test for H 0 inferences using the electromagnetic coun-terpart. As we head into the next generation of gravitationalwave detectors, the standard-siren method will improve the constraints on the Hubble constant. This new independent method may potentially resolve the Hubble tension problem or compel us to reevaluate our cosmological models, specifically ΛCDM.

Here, we summarize the statistical method adopted from Chen et al. (2018) , Gray et al. (2019) , and Abbott et al. (2019b) . The posterior probability on H 0 from N gravitational-wave (GW) events can be computed as

where x GW is the set of GW data and D GW indicates that the detection was made in the form of a GW. Here, p(H 0 ) is the prior on H 0 that we took to be either uniform or flat log over an interval [a,b] . p(N |H 0 ) is the likelihood of detecting N events given a value for H 0 . Using the same prior on the astrophysical rate of events as in Abbott et al. (2019b) , we drop the dependence of this term on H 0 .

For an individual gravitational wave event, the likelihood can be written as

The normalization factor p(D GW |H 0 ) in the denominator of Eq. A2 can be evaluated with the integral

where ρ th is the signal-to-noise ratio (SNR) threshold below which p(D GW |H 0 ) = 0. In our calculations we assume ρ th = 8, the default in gwcosmo.

In our case, we use the galaxy catalog method, in which the likelihood function p(x GW |D GW , H 0 ) can be expanded as 

with G andḠ denoting the cases where the host is in the catalog and where it is not, respectively.

The likelihood when the host galaxy is in the catalog, p(x GW |G, D GW , H 0 ), can be written as

where N gal is the number of galaxies considered in the catalog, Ω(α, δ) is the sky position angle, z i is the redshift of galaxy, s indicates that a GW has been emitted (distinct from detected), and m i and M are respectively the apparent and absolute magnitude of the galaxy. 

where λ is the rate evolution parameter, whose default value is 3.0 in gwcosmo. However, we hold the rate constant throughout all calculations. We define V c (z) as the co-moving volume contained within a redshift z.

Redshift Prior:

dz , if merger rate density = const.

The likelihood when the host galaxy is not in the catalog is defined as 

Equations A10 and A11 are the probabilities that the host is and is not in the galaxy catalog, respectively. The prior on the GW host galaxy sky location, p(Ω), is taken to be uniform across the sky. The prior on the absolute magnitude, p(M |H 0 ), is taken to be proportional to the Schechter luminosity function, the parameters for which are defined in Section 2.2. Here, m th is the apparent magnitude threshold of the flux-limited galaxy catalog. A complete description of the mathematics in gwcosmo is given in the Appendix of Gray et al. (2019) .

GRB Coordinates Network

We thank Benjamin Stahl and Keto Zhang (U.C. Berkeley) for helpful advice on writing this paper, as well as the anonymous referee whose suggestions improved its quality. Ignacio Magaña Hernandez (U.W. Milwaukee) discussed gwcosmo with us. We are also grateful to the lscsoft team for their efforts in providing public access to the gwcosmo code. Generous financial support for this work was provided by Steven Nelson, the Christopher R. Redlich Fund, and the Miller Institute for Basic Research in Science (U.C. Berkeley).