# A statistical standard siren measurement of the Hubble constant from the LIGO/Virgo gravitational wave compact object merger GW190814 and Dark Energy Survey galaxies

A. Palmese, J. Devicente, M.E.S. Pereira, J. Annis, W. Hartley, K. Herner, M. Soares-Santos, M. Crocce, D. Huterer, I. Magaña Hernandez, et al.

► To cite this version:

A. Palmese, J. Devicente, M.E.S. Pereira, J. Annis, W. Hartley, et al., A statistical standard siren measurement of the Hubble constant from the LIGO/Virgo gravitational wave compact object merger GW190814 and Dark Energy Survey galaxies. The Astrophysical Journal Letters, 2020, 900 (2), pp.L33. 10.3847/2041-8213/abaeff. hal-02899247

HAL Id: hal-02899247

https://hal.science/hal-02899247v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CrossMark

# A Statistical Standard Siren Measurement of the Hubble Constant from the LIGO/Virgo Gravitational Wave Compact Object Merger GW190814 and Dark Energy Survey Galaxies

A. Palmese<sup>1,2</sup>, J. de Vicente<sup>3</sup>, M. E. S. Pereira<sup>4</sup>, J. Annis<sup>1</sup>, W. Hartley<sup>5,6,7</sup>, K. Herner<sup>1</sup>, M. Soares-Santos<sup>4</sup>, M. Crocce<sup>8,9</sup>, D. Huterer<sup>10</sup>, I. Magaña Hernandez<sup>11</sup>, A. Garcia<sup>4</sup>, J. Garcia-Bellido<sup>12</sup>, J. Gschwend<sup>13,14</sup>, D. E. Holz<sup>2</sup>, R. Kessler<sup>2,15</sup>

O. Lahav<sup>6</sup>, R. Morgan<sup>16</sup>, C. Nicolaou<sup>6</sup>, C. Conselice<sup>17</sup>, R. J. Foley<sup>18</sup>, M. S. S. Gill<sup>19</sup>, T. M. C. Abbott<sup>20</sup>, M. Aguena<sup>13,21</sup>, S. Allam<sup>1</sup>, S. Avila<sup>12</sup>, K. Bechtol<sup>16</sup>, E. Bertin<sup>22,23</sup>, S. Bhargava<sup>24</sup>, D. Brooks<sup>6</sup>, E. Buckley-Geer<sup>1</sup>, D. L. Burke<sup>19,25</sup>

M. Carrasco Kind<sup>26,27</sup>, J. Carretero<sup>28</sup>, F. J. Castander<sup>8,9</sup>, C. Chang<sup>2,15</sup>, M. Costanzi<sup>29,30</sup>, L. N. da Costa<sup>13,14</sup>, T. M. Davis<sup>31</sup>, S. Desai<sup>32</sup>

<sup>32</sup> H. T. Diehl<sup>1</sup>, P. Doel<sup>6</sup>, A. Drlica-Wagner<sup>1,2,15</sup>, J. Estrada<sup>1</sup>, S. Everett<sup>18</sup>, A. E. Evrard<sup>10,33</sup>, E. Fernandez<sup>28</sup>, D. A. Finley<sup>1</sup>, B. Flaugher<sup>1</sup>, P. Fosalba<sup>8,9</sup>, J. Frieman<sup>1,2</sup>, E. Gaztanaga<sup>8,9</sup>, D. W. Gerdes<sup>10,33</sup>, D. Gruen<sup>19,25,34</sup>

R. A. Gruendl<sup>26,27</sup>, G. Gutierrez<sup>1</sup>, S. R. Hinton<sup>31</sup>, D. L. Hollowood<sup>18</sup>, K. Honscheid<sup>35,36</sup>, D. J. James<sup>37</sup>, S. Kent<sup>1,2</sup>, E. Krause<sup>38</sup>, K. Kuehn<sup>39,40</sup>, H. Lin<sup>1</sup>, M. A. G. Maia<sup>13,14</sup>, M. March<sup>41</sup>, J. L. Marshall<sup>42</sup>, P. Melchior<sup>43</sup>, F. Menanteau<sup>26,27</sup>

R. Miquel<sup>28,44</sup>, R. L. C. Ogando<sup>13,14</sup>, F. Paz-Chinchón<sup>27,45</sup>, A. A. Plazas<sup>43</sup>, A. Roodman<sup>19,25</sup>, M. Sako<sup>41</sup>, E. Sanchez<sup>3</sup>, V. Scarpine<sup>1</sup>

M. A. Troxel<sup>49</sup>, D. L. Tucker<sup>1</sup>, A. R. Walker<sup>20</sup>, W. Wester<sup>1</sup>, R. D. Wilkinson<sup>24</sup>, and J. Zuntz<sup>50</sup>

(DES Collaboration)

<sup>1</sup> Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA; palmese@fnal.gov

<sup>2</sup> Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

<sup>3</sup> Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

<sup>4</sup> Department of Physics, Brandeis University, Waltham, MA 02453, USA

<sup>5</sup> Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ernest Ansermet, CH-1211 Geneva, Switzerland

<sup>6</sup> Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

<sup>7</sup> Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland

<sup>8</sup> Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain

<sup>9</sup> Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain

<sup>10</sup> Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

<sup>11</sup> University of Wisconsin-Milwaukee, Milwaukee, WI 53201, USA

<sup>12</sup> Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

<sup>13</sup> Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

<sup>14</sup> Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

<sup>15</sup> Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA

<sup>16</sup> Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue Madison, WI 53706-1390, USA

<sup>17</sup> University of Nottingham, School of Physics and Astronomy, Nottingham NG7 2RD, UK

<sup>18</sup> Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

<sup>19</sup> SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

<sup>20</sup> Cerro Tololo Inter-American Observatory, NSF's National Optical-Infrared Astronomy Research Laboratory, Casilla 603, La Serena, Chile

<sup>21</sup> Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, S ao Paulo, SP, 05314-970, Brazil

<sup>22</sup> CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

<sup>23</sup> Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

<sup>24</sup> Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK

<sup>25</sup> Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA

<sup>26</sup> Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

<sup>27</sup> National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, IL 61801, USA

<sup>28</sup> Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain

<sup>29</sup> INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy

<sup>30</sup> Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

<sup>31</sup> School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

<sup>32</sup> Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

<sup>33</sup> Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

<sup>34</sup> Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA

<sup>35</sup> Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

<sup>36</sup> Department of Physics, The Ohio State University, Columbus, OH 43210, USA

<sup>37</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

<sup>38</sup> Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA

<sup>39</sup> Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

<sup>40</sup> Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

<sup>41</sup> Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

<sup>42</sup> George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA

<sup>43</sup> Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

<sup>44</sup> Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

<sup>45</sup> Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

<sup>46</sup> Austin Peay State University, Dept. Physics, Engineering and Astronomy, P.O. Box 4608 Clarksville, TN 37044, USA

<sup>47</sup> School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

<sup>48</sup> Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
<sup>49</sup>Department of Physics, Duke University Durham, NC 27708, USA

<sup>50</sup>Institute for Astronomy, University of Edinburgh, Edinburgh EH9 3HJ, UK

Received 2020 June 25; revised 2020 August 12; accepted 2020 August 13; published 2020 September 7

Abstract

We present a measurement of the Hubble constant $H_0$ using the gravitational wave (GW) event GW190814, which resulted from the coalescence of a 23 $M_\odot$ black hole with a 2.6 $M_\odot$ compact object, as a standard siren. No compelling electromagnetic counterpart has been identified for this event; thus our analysis accounts for thousands of potential host galaxies within a statistical framework. The redshift information is obtained from the photometric redshift (photo-z) catalog from the Dark Energy Survey. The luminosity distance is provided by the LIGO/Virgo gravitational wave sky map. Since this GW event has the second-smallest localization volume after GW170817, GW190814 is likely to provide the best constraint on cosmology from a single standard siren without identifying an electromagnetic counterpart. Our analysis uses photo-z probability distribution functions and corrects for photo-z biases. We also reanalyze the binary black hole GW170814 within this updated framework. We explore how our findings impact the $H_0$ constraints from GW170817, the only GW merger associated with a unique host galaxy. From a combination of GW190814, GW170814, and GW170817, our analysis yields $H_0 = 72.0_{-8.2}^{+12} \text{ km s}^{-1} \text{ Mpc}^{-1}$ (68% highest-density interval, HDI) for a prior in $H_0$ uniform between [20and140] km s$^{-1}$ Mpc$^{-1}$. The addition of GW190814 and GW170814 to GW170817 improves the 68% HDI from GW170817 alone by ~18%, showing how well-localized mergers without counterparts can provide a significant contribution to standard siren measurements, provided that a complete galaxy catalog is available at the location of the event.

Unified Astronomy Thesaurus concepts: Cosmology (343); Cosmological parameters (339); Gravitational waves (678); Hubble constant (758); Redshift surveys (1378); Surveys (1671); Gravitational wave astronomy (675)

1. Introduction

The first detection of gravitational waves (GW) from a binary black hole merger (GW150914; Abbott et al. 2016), and only 2 years later the first detection of a binary neutron star merger (GW170817; Abbott et al. 2017) with an associated electromagnetic counterpart (Arcavi et al. 2017; Coulter et al. 2017; LIGO Scientific Collaboration et al. 2017; Lipunov et al. 2017; Soares-Santos et al. 2017; Tanvir et al. 2017; Valenti et al. 2017), have generated tremendous excitement among the astrophysics community. One of the most appealing applications of gravitational wave detections is in cosmological analyses. The GW signal from mergers of neutron star and black hole binary systems are in fact absolute distance indicators, and can be used as “standard sirens,” as first proposed in Schutz (1986). The luminosity distance to the source can be inferred from the gravitational wave signal, and if a redshift measurement is also available (e.g., through identification of the host galaxy), we can measure the present rate of expansion of the Universe $H_0$ via the distance–redshift relation.

New and independent measurements of the Hubble constant are of great interest to cosmology. Measurements obtained with type Ia Supernovae (SN Ia) and inferred from the cosmic microwave background (CMB; e.g., Planck Collaboration et al. 2018; Freedman et al. 2019; Riess et al. 2019) have each reached the 1%–2% precision level, although the latest measurements from SN Ia (Riess et al. 2019) and CMB (Planck Collaboration et al. 2018) disagree at the 4.4σ level. Independent, precise, and accurate measurements of $H_0$ could help clarify whether the tension arises from beyond ΛCDM physics or unknown systematics (Freedman 2017; Mörtsell & Dhawan 2018; Verde et al. 2019).

The standard siren methodology was applied for the first time to the binary neutron star (BNS) merger GW170817 (Abbott et al. 2017a), thanks to the association of its electromagnetic (EM) counterpart to the nearby host galaxy NGC 4993 (e.g., Palmese et al. 2017). However, no other

compelling counterparts to GW events have been identified to date (e.g., Andreoni et al. 2020; Garcia et al. 2020; Morgan et al. 2020; Vieira et al. 2020). Events without counterparts can also be used for cosmological analyses, using a statistical approach first proposed in Schutz (1986). If a complete catalog of potential host galaxies exists within the event localization region, their redshift distribution can provide the redshift information needed to infer cosmological parameters from the distance–redshift relation. We refer to this method as the “dark” or “statistical” standard siren method, as opposed to the “bright” or “counterpart” case. Del Pozzo (2012) and Chen et al. (2018) provide a Bayesian framework that enables a measurement of $H_0$ with the statistical approach. This framework is implemented in Fishbach et al. (2019) using GW170817, and in Soares-Santos et al. (2019), where we measure the Hubble constant using a binary black hole (GW170814; Abbott et al. 2017b) for the first time. LIGO Scientific Collaboration & Virgo Collaboration (2019b) combine dark sirens from all of the binary black holes detected during O1 and O2, but find that the method could only bring a ~7% improvement to the GW170817 bright siren constraint due to the lack of complete public galaxy catalogs in the regions of interest and to the poor localization of most of the events. Eventually, a large sample of events, combined with wide-field ongoing and upcoming galaxy surveys (e.g., LSST, DESI, 4MOST), will enable precise cosmological measurements with dark and bright standard sirens (e.g., Holz & Hughes 2005; MacLeod & Hogan 2008; Nissanke et al. 2010; Del Pozzo 2012; Nissanke et al. 2013; Nishizawa 2017; Chen et al. 2018; Vitale & Chen 2018; Nair et al. 2018; Feeney et al. 2019; Mortlock et al. 2019; Palmese et al. 2019; Yu et al. 2020). In particular, Chen et al. (2018) forecast a 2% precision within 5 yr for bright standard sirens detected by LIGO/Virgo. Nair et al. (2018) predict a ~7% measurement with 25 binary black hole (BBH) events from next-generation GW detectors, while Yu et al. (2020) anticipate that a 1%–4% measurement will be possible with LIGO/Virgo BBHs if these mergers
originate in groups and clusters of galaxies. Interesting cosmological constraints can also be placed using gravitational wave compact binary mergers and large galaxy surveys by measuring their peculiar velocity power spectrum, along with their overdensity and cross-correlation power spectra, as proposed in Palmese & Kim (2020).

With the promises of gravitational wave cosmology in mind, the Dark Energy Survey (DES) collaboration and external collaborators launched the DES gravitational waves (DESGW) program. As part of this program, we seek optical emission from LIGO/Virgo events using the Dark Energy Camera (DECam; Flaugher et al. 2015). Our searches resulted in the independent discovery of the kilonova associated with GW170817 (Soares-Santos et al. 2017), and in some of the most stringent limits for optical emission from a binary black hole (Cowperthwaite et al. 2016; Soares-Santos et al. 2016; Doctor et al. 2019) and other compact object (Morgan et al. 2020) mergers.

In this article, we measure $H_0$ using the gravitational wave event GW190814 (LIGO Scientific Collaboration & Virgo Collaboration 2019a; Abbott et al. 2020), which resulted from the inspiral and merger of a $23 M_\odot$ black hole with a $2.6 M_\odot$ compact object at $241_{-45}^{+41}$ Mpc (90% credible interval; Abbott et al. 2020). The secondary component was either the lightest black hole or the heaviest neutron star ever found in a binary system. This event is particularly interesting for a dark standard siren analysis because counterpart searches have not identified a convincing counterpart thus far (e.g., Andreoni et al. 2020; Ackley et al. 2020; Gomez et al. 2019; Morgan et al. 2020; Vieira et al. 2020; Watson et al. 2020) and because its localization volume is the second smallest after GW170817. Additionally, the event localization region falls within the DES footprint, making DES galaxy catalogs an ideal sample for a measurement of $H_0$. Abbott et al. (2020), which became publicly available almost concurrently with this work, provides a brief description of a dark siren analysis for this event using an inhomogeneous galaxy sample, mostly incomplete at the redshift range of interest. In this work, we provide a detailed standard siren analysis for GW190814 that takes advantage of a more complete galaxy catalog. Compared to other previous works on dark standard sirens, we also improve upon the treatment of the photometric redshifts (photo-z's) by taking into account their full probability distribution function (PDF) rather than a Gaussian approximation of the former, and by correcting for photo-z biases in the data. This methodology ensures a more accurate recovery of the true redshift distribution of galaxies, as it has been extensively studied in the literature, particularly within the context of weak gravitational lensing analyses (e.g., Lima et al. 2008; Oyaizu et al. 2008; Cunha et al. 2009). We also reanalyze GW170814 using DES galaxies within this updated framework and provide a combined $H_0$ measurement using three events.

We describe our data set in Section 2 and the methods used in Section 3. Our results and discussion follow in Section 4, and the conclusions are in Section 5. We assume a flat $\Lambda$CDM cosmology with $\Omega_m = 0.3$ and $H_0$ values in the 20–140 km s$^{-1}$ Mpc$^{-1}$ range. When not otherwise stated, quoted error bars represent the 68% credible interval (CI).

## 2. Data

### 2.1. The LIGO/Virgo GW Data

The gravitational wave data used in this work come from the publicly available sky maps. First, we use the LALInference sky map (LIGO Scientific Collaboration & Virgo Collaboration 2019a)<fnref n="51" /> released in 2019 August. The area of the sky enclosing 90% of the localization probability is 23 deg², and it is shown in Figure 1 together with the GW170814 90% CI contours and the DES footprint (in blue). Angular and distance probabilities are provided in HEALPix (Górski et al. 2005) pixels, where the distance probability distribution is approximated by a Gaussian along each line of sight. The maximum probability pixel is centered at R.A., decl. = (12.832, −25.241) deg. Marginalized over the whole sky, the luminosity distance has a mean value of 267 Mpc and a standard deviation of 52 Mpc in the Gaussian approximation. The resulting 90% CI comoving volume is $9.2 \times 10^4$ Mpc$^3$, which is two orders of magnitude larger than the GW170817 volume, but two orders of magnitude smaller than GW170814. We derive the redshift range of interest for this analysis by considering the highest value of $H_0$ in our prior range (i.e., 140 km s$^{-1}$ Mpc$^{-1}$) with a flat $\Lambda$CDM cosmology having $\Omega_m = 0.3$, and the high limit of the luminosity distance at 90 and 99.7% CI. These two CI limits correspond to $z = 0.15$ and 0.18, implying that our analysis is sensitive to galaxies at redshifts below these values.

After submission of this manuscript, we added results using the sky map released by the LVC in 2020 July<fnref n="52" /> (hereafter called the *higher modes map*). This map provides more precise parameter estimation thanks to the higher-order modes detected in the offline analysis by the LVC (Abbott et al. 2020). As a result, the 90% CI localization area becomes 18.5 deg², and the luminosity distance is $241 \pm 26$ Mpc (mean and standard deviation). The updated 90% CI comoving volume is about a third of the one from the LALInference map, $3.2 \times 10^4$ Mpc$^3$.

### 2.2. The DES Data

The DES<fnref n="53" />(The Dark Energy Survey Collaboration 2005, 2016) is an optical near-infrared survey that imaged 5000 deg² of the South Galactic Cap over six years in the *grizY* bands. DES used the ~3 deg² DECam, mounted on the Blanco 4-m telescope at the Cerro Tololo Inter-American Observatory (CTIO) in Chile. In our analysis, we use data from the first 3 yr of observations (“Y3”; DES Collaboration et al. 2018; 2013 September–2016 February).

Our reduced data products are from the DES Data Management (DESDM) pipeline (Morganson et al. 2018), which includes calibration of the single-epoch images, background subtraction, co-addition of the background-subtracted images, and cutting the co-added images into tiles. The source catalog was created with SOURCE EXTRACTOR (SEXTRACTOR, Bertin & Arnouts 1996), which detects objects on the *riz* co-added images. For galaxies, the median $10\sigma$ limiting magnitudes in the Y3 data are $g = 24.33, r = 24.08, i = 23.44, z = 22.69$, and $Y = 21.44$ mag (DES Collaboration et al. 2018). The data were carefully selected to produce a Y3 “gold” catalog, as

<footnote n="51">⁵¹ https://gracedb.ligo.org/superevents/S190814bv/view/</footnote>

<footnote n="52">⁵² https://dcc.ligo.org/LIGO-P2000230/public</footnote>

<footnote n="53">⁵³ www.darkenergysurvey.org</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** *Left:* LIGO/Virgo GW dark standard sirens analyzed in this paper, where the shaded regions represent the 90% CI localization from the sky maps. For GW190814, the localization comes from the LALInference sky map. The blue contour represents the DES footprint. *Top right:* difference between several inferred redshift distributions $dN/dz$ and the distribution $(dN/dz)_{com}$ for a uniform-number-density case with $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$. This difference is shown to highlight the presence of overdensities and underdensities in the $dN/dz$ distribution. The redshifts have been inferred with the DNF method using DES photometry. The redshift range shown is the one relevant for an $H_0$ prior extending out to 140 km s$^{-1}$ Mpc$^{-1}$. The dotted line represents the distribution obtained when considering the mean photo-z value of each galaxy as provided by DNF. The pink curve is given by Monte Carlo (MC) sampling the photo-z PDF for each galaxy, assuming that this is a Gaussian centered on the mean value (from the dotted line) with standard deviation given by the $1\sigma$ uncertainty provided by DNF. The blue solid curve is a sum of the full, non-Gaussian, redshift PDFs of each galaxy, and is used in our final results. *Bottom right:* DNF photo-z bias $\Delta z = \langle z_{sp} \rangle - \langle z_{ph} \rangle$ in bins of photo-z measured on a validation sample of ~140,000 galaxies with spectroscopic data in the DES footprint. The photo-z estimator used to compute the bias is the mean DNF PDF value.</figcaption></figure>

described in I. Sevilla-Noarbe et al. 2020, in preparation. In this work, we use the gold catalog, and the photometry is derived through the Single Object Fitting (SOF) pipeline that relies on the ngmix code.<fnref n="54" /> The SOF fluxes have been used to compute photometric redshifts using the Directional Neighborhood Fitting (DNF; De Vicente et al. 2016) method, as described in more detail in the next subsection.

Galaxy properties for the Y3 sample (A. Palmese et al. 2020, in preparation) are derived using the SOF photometry and by fixing the redshift to the DNF photo-z mean value, or spectroscopic redshift, where available. In particular, stellar mass and absolute magnitudes are derived through a broadband Spectral Energy Distribution (SED) fitting of galaxy magnitudes with LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006). Estimates of the galaxy properties used here from DES data alone have been tested and studied in several DES works (Palmese et al. 2016; Etherington et al. 2017; Palmese et al. 2020). We add a 0.05 systematic uncertainty in quadrature to the magnitudes, to account for systematic uncertainties in magnitude estimation and model variance. The templates used for the galaxies’ SED fitting are the simple stellar populations (SSP) from Bruzual & Charlot (2003), with three metallicities (0.2 Z⊙, Z⊙, and 2.5 Z⊙), a Chabrier (2003) Initial Mass Function (IMF), and a Milky Way (Allen 1976) extinction law with five different values between 0 and 0.5 for the $E(B-V)$ reddening. The star formation history (SFH) chosen is exponentially declining with the age of the galaxy t as $e^{-t/\tau}$ with $\tau = 0.1, 0.3, 1, 2, 3, 5, 10, 15$, and 30 Gyr.

As discussed in DES Collaboration et al. (2018), the Y3 gold catalog source list is 95% complete for galaxies within our apparent magnitude limit of $r < 23.35$. By converting the apparent magnitude source completeness into a completeness in redshift intervals using the DNF photo-z’s, in Soares-Santos et al. (2019) we found that our sample with $r < 23.35$ is >93% complete within $z < 0.26$. As a result, our catalog is highly complete over the redshift range of interest; thus we do not apply completeness corrections. However, the sample is still magnitude limited, implying that intrinsically fainter galaxies will be more easily observed at the lowest redshifts. We overcome this problem by defining a volume-limited sample, obtained by applying a luminosity cut. Following Pozzetti et al. (2010) and Hartley et al. (2013), we identify galaxies that are bright enough to be complete and representative of the real galaxy population. In order to facilitate a comparison with GW170814, we adopt the same luminosity cut used in Soares-Santos et al. (2019) because such a limit is valid in the range $0 < z < 0.26$, which includes the redshifts of interest here. This cut is made at −17.2 (computed assuming $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, and it is rescaled accordingly for other $H_0$ values) in r-band absolute magnitude and ~3.8 × 10⁸ M⊙ in stellar mass. As in Soares-Santos et al. (2019), the sample remaining after the cut contains >77% of the stellar mass in the volume, assuming the galaxies follow a Schechter stellar mass function with the best fit values from Weigel et al. (2016).

In the next subsection and in particular for Figure 2, we have used the LALInference map. As the higher modes map has a smaller spatial localization and smaller luminosity distance uncertainties, the ~1800 galaxy sample encompassing the 90%

<footnote n="54">⁵⁴ https://github.com/esheldon/ngmix</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Hubble constant posterior distributions for GW190814 (blue curves) and GW170814 (gray curves), comparing results using the Gaussian redshift PDF versus the full PDFs. The latter case does not depend on the choice of redshift cut, as long as this contains the range of z supported by the $H_0$ prior. The 68% HDI of all PDFs is shown by the dashed lines. Posterior s have been rescaled for visualization purposes. The sky map used here for GW190814 is the LALInference map released in 2019 August.</figcaption></figure>

CI volume using it is a strict subset of the ~3,800 galaxy sample selected using the LALInference.

### 2.2.1. The Redshift Catalog

The DNF method applied to Y3 data provides redshift information for each galaxy in the form of a PDF, from which we compute a mean redshift, and half of the central 68th percentile width. Our method requires having individual galaxies' PDFs, as will be clear in Section 3. We use the full PDF for each galaxy and contrast the results with a Gaussian approximation of the PDF using mean redshifts with the $1\sigma$ values. Where available, we use public spectroscopic redshifts to complement the sample.

The final galaxy redshift distributions are shown in the top-right panel of Figure 1. There are 2,684 galaxies in the GW190814 90% probability region, supporting $H_0$ values out to $140 \text{ km s}^{-1} \text{ Mpc}^{-1}$ ($z \sim 0.135$). Of these, ~700 have available spectroscopic redshifts from 2dF, 2dFLens, and 6dF (Colless et al. 2001; Jones et al. 2009). The top-right plot in Figure 1 shows the difference between several inferred redshift distributions $dN/dz$ and the distribution $(dN/dz)_{\text{com}}$ for a uniform-number-density case (i.e., a distribution which is uniform in comoving volume) with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. The different lines show the results from different redshift estimators. The dotted line represents the $dN/dz$ given by the mean photo-z value of each galaxy, as provided by DNF. The pink curve is given by Monte Carlo (MC) sampling a Gaussian approximation of the photo-z PDF for each galaxy. In this case, we assume that the PDF is a Gaussian with width given by the $1\sigma$ uncertainty provided by DNF. The blue solid curve is a sum of the full, non-Gaussian, redshift PDFs of each galaxy.

We have tested the aforementioned redshift estimators using a spectroscopic sample of 140,000 galaxies matched to DES Y3 objects as a validation sample (Gschwend et al. 2018). This sample does not contain any galaxies used for training DNF. The metric that we choose to evaluate the performance of the photo-z algorithm in recovering the redshift distribution is

$$ \Delta \bar{N}(z)(\%) = \sum_k \frac{|N_{\text{ph}}(z_k) - N_{\text{sp}}(z_k)|}{N_{\text{sp}}(z_k)} \times 100, \quad (1) $$

where the sum extends over the redshift bins of interest, and $N_{\text{ph}}(z_k)$ and $N_{\text{sp}}(z_k)$ are the photometric and spectroscopic sample galaxy counts in the $k$th bin, respectively. In the case of the full photo-z PDF, $N_{\text{ph}}(z_k)$ is computed by summing the contributions of single galaxies' PDFs in each bin $k$. We consider redshift bins of width 0.05 out to $z < 0.6$. Compared to using the mean photo-z values for $N_{\text{ph}}(z_k)$, this metric is improved by more than a factor of two when $N_{\text{ph}}(z_k)$ is summed by integrating the PDF of all galaxies over the $k$th redshift bin. Using MICE simulations (Crocce et al. 2015; Fosalba et al. 2015), we also find that the $\Delta \bar{N}(z)$ metric is significantly reduced (again by a factor of ~2) when using the full PDF versus the mean DNF values. The magnitude of the improvement is equivalent to adding $u$ band to the measurements. This is due to the fact that the full PDF captures the effect of color-redshift degeneracies. On the other hand, these degeneracies are not well captured when only a point estimate such as the mean is taken into account, or when the PDF is approximated by a Gaussian (Buchs et al. 2019; J. Gschwend et al. 2020, in preparation). We conclude that accounting for full PDF is more accurate than using point estimates, as they are able to best reproduce the redshift distribution of galaxies. This conclusion is broadly supported by several works on photometric surveys (e.g., Lima et al. 2008; Oyaizu et al. 2008; Cunha et al. 2009).

Using the spectroscopic validation sample, we estimate the extent of systematic biases that propagate to a bias in the $H_0$ posterior. We find that in the redshift range of interest for this work, $z < 0.3$, the bias $\Delta z = \langle z_{\text{sp}} \rangle - \langle z_{\text{ph}} \rangle$ in bins of photo-z with width 0.1 is of the order ~0.002. The photo-z estimator used to compute this bias is the mean DNF value. We marginalize over the redshift dependent bias shown in the bottom-right panel of Figure 1.

We notice that the redshift distributions show the large “galaxy wall” already seen for GW170814 (see Figure 1 of Soares-Santos et al. 2019), spanning most of the DES footprint and confirmed by spectroscopic data, at $z \sim 0.06$. This overdensity is also present in Bayesian Photometric Redshift (BPZ; Benítez 2000) redshifts, which is a template-based method, proving that the structure is not a result of machine learning training in DNF.

## 3. Method

In this work, we follow the approach of Chen et al. (2018) and slightly adapt it for our purposes, similarly to Soares-Santos et al. (2019). The posterior probability of $H_0$, given the GW data $d_{\text{GW}}$ from a single event detection and EM data $d_{\text{EM}}$ from a galaxy survey, can be written as

$$ p(H_0|d_{\text{GW}}, d_{\text{EM}}) \propto p(d_{\text{GW}}, d_{\text{EM}|H_0})p(H_0), \quad (2) $$

following Bayes' theorem. The joint likelihood $p(d_{\text{GW}}, d_{\text{EM}|H_0})$ can be written as the product of two individual likelihoods, $p(d_{\text{GW}|H_0})$ and $p(d_{\text{EM}|H_0})$, since the GW and EM data are independent. The GW likelihood is marginalized over all variables except for the true luminosity distance $d_L$ and solid angle $\hat{\Omega}_{\text{GW}}$ of the GW source, while the EM likelihood explicitly depends on the true redshift $z_i$ and solid angle $\hat{\Omega}_i$ for each galaxy $i$. We refer to $\hat{\Omega}_i$ as “solid angles,” but these are vectors with the angular position of the source/galaxy as direction, subtending the area (~3 $\times 10^{-3} \text{ deg}^2$) of each
HEALPix pixel. We assume that the GW source is located within one of the galaxies in the galaxy catalog. Therefore, the location and distance of the GW source can be directly related to solid angle $\hat{\Omega}_i$ and the redshift (through the cosmology) of each galaxy. Finally, we marginalize over the choice of galaxy $i$, over the true luminosity distance and over the position of the GW source, and write the joint, marginal likelihood as

$$ p(d_{\text{GW}}, d_{\text{EM}}|\{z_j, \hat{\Omega}_j\}, H_0) \\ \propto \sum_i w_i \int dd_L d\hat{\Omega}_{\text{GW}} p(d_{\text{GW}}|d_L, \hat{\Omega}_{\text{GW}}) \\ \times p(d_{\text{EM}}|\{z_j, \hat{\Omega}_j\}) \delta_D(d_L - d_L(z_i, H_0)) \delta_D(\hat{\Omega}_{\text{GW}} - \hat{\Omega}_i), \quad (3) $$

where $\{z_j, \hat{\Omega}_j\}$ are all of the galaxies' redshift and solid angle, $w_i$ are weights that represent the relative probability that the $i$th galaxy hosts a GW source, and $\delta_D$ is the Dirac delta function. In principle, the $w_i$ weights would be based on galaxy properties, but as we do not know whether GW host galaxy properties such as star formation rate or luminosity differ from the general population, here we take the weights to be uniform across all galaxies. We therefore omit them in the following equations.

Equation (3) needs to be marginalized over the galaxies' redshifts and sky positions, which require a prior $p(z_i, \hat{\Omega}_i)$. To first order, galaxies are uniformly distributed in comoving volume $V$, assuming that our sample is volume limited within $V_{\text{max}}$:

$$ p(z_i, \hat{\Omega}_i) dz_i d\hat{\Omega}_i \propto \frac{1}{V_{\text{max}}} \frac{d^2V}{dz_i d\hat{\Omega}_i} dz_i d\hat{\Omega}_i \\ \propto \frac{1}{V_{\text{max}}} \frac{r^2(z_i)}{H(z_i)} dz_i d\hat{\Omega}_i, \quad (4) $$

where $r$ is the comoving distance to the galaxy.

If we assume that the galaxies' positions $\{\hat{\Omega}_j\}$ are known with exquisite precision, these can easily be marginalized over with a delta function about the observed values, reducing the marginal EM likelihood to $p(d_{\text{EM}}|\{z_j\})$. The assumption on the precision of the galaxies' position is realistic given that spatial probabilities for the GW sources vary significantly over scales that are usually larger than a galaxy's size at the redshifts of interest (in this specific work, they are uniform in each HEALPIX pixel). The marginal EM likelihood is then given by

$$ p(d_{\text{EM}}|\{z_j\}) = \prod_k p(\tilde{z}_k|z_k), \quad (5) $$

where the EM data $d_{\text{EM}}$ here only refers to the measured (photometric or spectroscopic) redshift $\tilde{z}_k$ for each galaxy $k$. In the photo-z Gaussian approximation case mentioned previously, we approximate Equation (5) with a product of Gaussian distributions, $N$, for each galaxy, centered around the observed redshift values $z_{\text{obs},k}$ with a width given by the redshift's uncertainty $\sigma_{z,k}$ for each galaxy $k$.

We compute the marginal GW likelihood $p(d_{\text{GW}}|d_L, \hat{\Omega})$ according to Singer et al. (2016):

$$ p(d_{\text{GW}}|d_L, \hat{\Omega}) \propto p(\hat{\Omega}) \frac{1}{\sqrt{2\pi}\sigma(\hat{\Omega})} \\ \times \exp\left[-\frac{(d_L - \mu(\hat{\Omega}))^2}{2\sigma^2(\hat{\Omega})}\right] N(\hat{\Omega}), \quad (6) $$

with the sky map providing the position probability $p(\hat{\Omega})$, mean $\mu$, normalization $N$, and scale $\sigma$ at each position.

We take into account selection effects arising from the detection of GW events through a normalization factor $\beta(H_0)$, by which we divide the right-hand side of Equation (2). Omitting these effects would bias our inference of $H_0$, in particular because the detection of GW events depends on the Hubble constant, among other variables (e.g., Chen et al. 2018; Mandel et al. 2019; Mortlock et al. 2019). Following Chen et al. (2018) and Mandel et al. (2019), we normalize the likelihood over all possible GW and EM data with a $[\beta(H_0)]^{-1}$ factor. In this work, the galaxy catalog is volume limited and complete beyond the maximum observable distance for the GW events, so that the normalization factor becomes

$$ \beta(H_0) = \frac{V[d_{L,\text{GW}}^{\text{max}}(H_0)]}{V_{\text{max}}(H_0)}, \quad (7) $$

where $V[d_{L,\text{GW}}^{\text{max}}(H_0)]$ is the maximum observable volume for the GW events considered. This volume is computed assuming the distance reach for 30 + 30 $M_{\odot}$ BBH from KAGRA Collaboration et al. (2018), namely 910 and ~1100 Mpc for O2 and O3, respectively. We have verified that the choice of the masses (and, in turn, the exact value of the distance reach) does not have a significant effect on our final results. In fact, the volume dependence on $H_0$ has the most impact on the Hubble constant posterior, while the exact distance value acts similarly to a normalization factor. However, future analyses that aim at precision measurements of cosmological parameters should include a more sophisticated calculation of $\beta(H_0)$, where the detectability of events following a realistic BBH mass function (rather than a Dirac delta centered on a specific mass) and the distributions of the other relevant binary parameters are taken into account. Chen et al. (2018), Fishbach et al. (2019), and Gray et al. (2020) have already performed more general derivations of this term, but we note that Fishbach et al. (2019) show that the volumetric term used in this work is a good approximation for events at distances that are typical for LIGO/Virgo BNSs, which can be extended to the case of GW190814.

We can rewrite Equation (2) as

$$ p(H_0|d_{\text{GW}}, d_{\text{EM}}) \propto \frac{p(H_0)}{V[d_{L,\text{GW}}^{\text{max}}(H_0)]} \\ \times \sum_i \frac{1}{Z_i} \int dz_i p(d_{\text{GW}}|d_L(z_i, H_0), \hat{\Omega}_i) p(d_{\text{EM}}|z_i) \frac{r^2(z_i)}{H(z_i)}, \quad (8) $$

where $Z_i = \int p(d_{\text{EM}}|z_i)r^2(z_i)/H(z_i) dz_i$ are evidence terms that correctly normalize the posterior. We further take into account photometric redshift systematics in the form of a bias $\Delta z$ and marginalize over it. We model the bias prior $p(\Delta z)$ with a Gaussian, where the mean and standard deviation take the
same value, which is the bias evaluated from the DES Y3 validation sample presented in Section 2.2. For simplicity of notation, we remove the subscript *i* for the redshifts inside the integral, since what really differentiates the integral for each galaxy is the different EM likelihood $p_i(d_{EM}|z)$:

$$ p(H_0|d_{GW}, d_{EM}) \propto \frac{p(H_0)}{V[d_{L,GW}^{\max}(H_0)]} \times \sum_i \frac{1}{Z_i} \int dz \, d\Delta z \, p(d_{GW}|d_L(z, H_0), \hat{\Omega}_i) \times p_i(d_{EM}|z, \Delta z) p(\Delta z) \frac{r^2(z)}{H(z)}. \quad (9) $$

The redshift PDF, shifted in redshift by the observed bias $\Delta z$, enters in Equation (9) through $p_i(d_{EM}|z, \Delta z)$.

We note that in each sky map pixel, the marginal GW likelihood of Equation (9) takes the same value for all galaxies in that pixel; therefore the sum only becomes relevant for the EM contribution to the likelihood. That sum, which aside from the bias term is effectively the sum over the photo-z PDFs, can be approximated with the $dN/dz$, as mentioned in Section 2.2. This explains why in that section we focused on validating the redshift distribution.

We also stress that despite having dropped the spatial dependence (on $\hat{\Omega}$) of the EM likelihood, $p_i(d_{EM}|z)$ still differs from galaxy to galaxy and therefore contains such dependence through the galaxy’s position. As a result, the $dN/dz$ interpretation mentioned previously is not uniform across the sky but differs from pixel to pixel.

We can extend this formalism and combine data from a sample of multiple GW events if we assume that the data $\{d_{GW,i}\}$ from each event *i* are independent of each other, and the EM data $d_{EM}$ use a fixed galaxy catalog for all events. If we start again from Equation (2) and add a marginalization over the redshift $z_k$ and position $\hat{\Omega}_k$ for each galaxy *k* out of the N galaxies in the galaxy catalog,

$$ \begin{aligned} p(H_0|\{d_{GW,i}\}, d_{EM}) &\propto p(H_0) \prod_j p(d_{GW,j}, d_{EM}|H_0) \\ &\propto p(H_0)p(d_{EM}|H_0) \prod_j p(d_{GW,j}|H_0) \\ &\propto p(H_0) \int d^N z_k d^N \hat{\Omega}_k p(z_k, \hat{\Omega}_k) \\ &\times p(d_{EM}|H_0, \{z_k, \hat{\Omega}_k\}) \left[ \prod_j p(d_{GW,j}|H_0, \{z_k, \hat{\Omega}_k\}) \right]. \end{aligned} \quad (10) $$

From this equation, the same steps followed between Equations (2) and (9) can be repeated in the same fashion in order to recover the final posterior for the combined sample.

In our analysis, we assume a flat prior on $H_0$ within [20,140] km s⁻¹ Mpc⁻¹, while all the other cosmological parameters are fixed (flat $\Lambda$CDM cosmology with $\Omega_m = 0.3$ and $\Omega_\Lambda = 0.7$). The choice of this prior is mostly dictated by the choices made in previous standard siren measurements, in particular Soares-Santos et al. (2019) and LIGO Scientific Collaboration & Virgo Collaboration (2019b), in order to ensure an easier comparison of the results. Similar to these works, the width of this prior ensures that the result is mostly informed by the LIGO/Virgo and DES data, instead of by external constraints.

In order to avoid confirmation bias, we blind our analysis when estimating the $H_0$ posterior by randomly displacing the

values of the Hubble constant by an unknown amount. We only unblind after establishing confidence in the methodology and data sets with the blinding in place.

## 4. Results and Discussion

The resulting posterior distributions for $H_0$ using GW170814 and GW190814 are shown in Figure 2. We contrast results for both the full redshift PDF used in this work and the Gaussian approximation of the redshift PDF (as done in Soares-Santos et al. 2019 and LIGO Scientific Collaboration & Virgo Collaboration 2019b). The difference in the shape of the posterior is most significant at the lowest $H_0$ values. The final 68% highest-density intervals (HDI; shown in Figure 2 by the dashed lines) are consistent with each other. It is thus unclear at this level of precision whether differences between the PDF choices are statistically significant. However, the Gaussian approximation does not accurately reproduce a realistic photo-z PDF of galaxies, which may contain multiple peaks due to color-redshift degeneracies. Similarly, the sum of these Gaussians does not accurately reproduce the true $dN/dz$. Let us note that the $dN/dz$ is effectively what enters Equation (7) if we approximate it with the sum of the redshift PDFs, as often done in weak lensing (e.g., Benjamin et al. 2013). On the contrary, the sum of redshift PDFs of single galaxies, or equivalently MC draws from those, are optimized to reproduce the true $dN/dz$ in photometric surveys such as DES and CFHTLenS (e.g., Benjamin et al. 2013; Hoyle et al. 2018). Moreover, the dark siren method with the Gaussian approximated PDFs suffers from a dependence on the redshift cut applied, as was already noted in Soares-Santos et al. (2019). This effect is due to the contribution of the Gaussian tails, which becomes non-negligible at the high redshift end—that is, at the high-$H_0$ end of Figure 2, where one marginalizes over thousands of galaxies. Another effect of the Gaussian distributions is the flattening of the posterior, as a consequence of the overdensities being suppressed in the $dN/dz$. The results shown here for the Gaussian approximation case correspond to a redshift cut at the 90% high limit in luminosity distance, assuming the highest $H_0$ value in our prior (cf. Section 2.1). On the other hand, the posterior resulting from the use of full photo-z PDFs is stable against different choices of redshift cuts (changes are at the sub-percent level). For these reasons, we choose to adopt the full PDF result as our fiducial result and thus reanalyze GW170814 in light of these findings. For both events in our fiducial result, we adopt a conservative redshift cut corresponding to the high bound in luminosity distance at >99.7% CI for the highest $H_0$ value considered here.

In the future, it will be preferable to use a hierarchical, fully Bayesian framework for the galaxies photo-zs that could be more easily incorporated into the statistical method used in this work. In fact, we notice that while the DNF PDFs have been carefully calibrated, they do not allow a flexible implementation for Bayesian inference and contain an “implicit prior (e.g., Schmidt et al. 2020) from the training sample that cannot be disentangled. The DES Collaboration plans to provide photo-zs from methods that fit more easily within a Bayesian scheme (along with a Self-Organizing Map method, Buchs et al. 2019; and BPZ), and we plan on implementing these in future analyses.

We note that the GW190814 posterior shows a clear peak between 50 and 80 km s⁻¹ Mpc⁻¹ for the PDF case, corresponding to the large overdensity shown in the top-right panel of
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Hubble constant posterior distribution for GW190814 obtained by marginalizing over ~1800 possible host galaxies from DES (light blue line). The dark blue curve represents the posterior obtained by the combination of GW190814 and GW170814 using DES galaxies. The maximum a posteriori and its 68% CI for the combined result is $H_0 = 78_{-13}^{+57}$ km s⁻¹ Mpc⁻¹ using a flat prior in the range [20,140] km s⁻¹ Mpc⁻¹. The posterior obtained by Nicolaou et al. (2020) for the bright standard siren event GW170817, associated to one galaxy, is shown in gray. A combination of the two dark sirens considered in this work with GW170817 is shown in red and gives $H_0 = 72.0_{-8.2}^{+12}$ km s⁻¹ Mpc⁻¹. The addition of the dark sirens provides a ~18% improvement to the 68% CI from GW170817 alone. The maximum a posteriori is represented by the solid vertical line. Posteriors have been rescaled for visualization purposes. The 68% CI of all PDFs is shown by the dashed lines. The GW190814 sky map used for this plot is the latest map released by the LVC in 2020 July. Constraints from *Planck* (Planck Collaboration et al. 2018) and Riess et al. (2019) (R19) at $1\sigma$ are shown in purple boxes.</figcaption></figure>

Figure 1. On the other hand, both posteriors for GW170814 appear to be flatter. One main difference between the two events is that GW170814 has a significantly larger localization volume (i.e., we marginalize over a larger number of galaxies; ~77,000 versus ~2,700). Posteriors for events that are not well localized are more likely to provide flatter, less informative posteriors because the overdensities are more likely to be washed out in the marginalization over the galaxies, regardless of the photo-z estimator chosen.

The marginalization over the photo-z bias does not have a significant effect at the current level of precision, confirming our expectations. In fact, a redshift bias of $\Delta z \sim 0.002$ is expected to cause a bias on $H_0$ of the level $\Delta H_0 \sim c\Delta z/d_L \sim 2.5$ km s⁻¹ Mpc⁻¹ for GW190814 at $d_L \sim 240$ Mpc. However, the photo-z bias also depends on redshift, and it could have a more complicated effect on the final $H_0$ posterior. It is thus important to take into account these biases as we start combining events.

We note that the posteriors in Figure 2 show significant probability at the high-$H_0$ end. This effect is intrinsic to the dark siren method; thus we do not widen the prior to larger values of $H_0$. In fact, the GW analysis only provides a luminosity distance estimate, which is consistent with arbitrarily large values of $H_0$. These values will be supported in the posterior from a dark siren analysis if there are galaxies at sufficiently large redshifts, which is the case for DES galaxies that extend out to $z > 1$. The advantage of using full PDFs is that one can extend the $H_0$ prior, and thus the redshift range considered, out to larger values without biasing the final result, since the shape of the $H_0$ posterior in a given interval is insensitive to the redshift cut.

Using the PDF stacking method, we then provide an $H_0$ posterior for GW190814 using the higher modes map. This result is shown in light blue in Figure 3. The mode and 68% CI

of the distribution result in $H_0 = 78_{-13}^{+57}$ km s⁻¹ Mpc⁻¹. The differences between this posterior and the dark blue curve in Figure 2, computed with the LALInference map, can be explained as follows. First, the mean value of the luminosity distance marginalized over the whole sky is shifted toward lower values by $\Delta d_L = 26$ Mpc in the latest map, and this results in a shift of the $H_0$ estimate of the order $\sim c z (\Delta d_L) / d_L^2 \sim 8$ km s⁻¹ Mpc⁻¹ toward larger values. Second, the uncertainty on the luminosity distance in the LALInference map is almost double that in the latest map, and a larger distance uncertainty has the effect of smoothing out peaks in the $H_0$ posterior, as shown with simulations in (Soares-Santos et al. 2019). We note that the higher modes map provides more precise distance and localization measurements (Abbott et al. 2020) and thus a more precise $H_0$ measurement. Therefore, we adopt this result as our final constraint.

Our estimate on $H_0$ from GW190814 is consistent with those from Section 6.4 of Abbott et al. (2020), that appeared shortly before submission of this paper. That study uses an independent pipeline and a smaller galaxy sample of ~472 galaxies (in the 90% CI) from the GLADE catalog (Dálya et al. 2018), which is mostly incomplete at the distance of GW190814. Their analysis also finds an overdensity of galaxies around the same redshift as in this work, resulting in an $H_0$ posterior ($H_0 = 75_{-13}^{+59}$ km s⁻¹ Mpc⁻¹) fully consistent with our results.

The final posterior distribution resulting from combining the dark sirens GW190814 with GW170814 is shown in dark blue in Figure 3. The maximum a posteriori and the 68% HDI around it is $H_0 = 77_{-22}^{+41}$ km s⁻¹ Mpc⁻¹ for the combination of both dark sirens (GW190814 with GW170814), using a flat prior in the range [20, 140] km s⁻¹ Mpc⁻¹ in all cases.<fnref n="55" /> We further combine these two dark siren events with the bright standard siren analysis of GW170817 from Nicolaou et al. (2020), which includes a more careful treatment of the peculiar velocity for NGC 4993 compared to Abbott et al. (2017a). We adapt the $H_0$ prior of Nicolaou et al. (2020) to be the same used in this work. The GW170817-only posterior is shown in Figure 3 by the gray shaded region, while the combined result from GW170817 with the dark sirens analyzed in this work (GW170814 plus GW190814) is given by the red solid line. We find that the addition of these dark sirens brings an ~18% improvement to the 68% HDI from GW170817 alone. For reference, we show constraints from *Planck* (Planck Collaboration et al. 2018) and Riess et al. (2019; R19) at $1\sigma$. Our results are broadly consistent with both estimates, which is expected, given the large uncertainty of the standard siren measurements. In the future, this comparison will become more interesting as we combine hundreds of events.

The final $H_0$ results from this and other standard siren analyses are summarized in Table 1. As expected, GW170817 is the most constraining event, followed by GW190814 and GW170814, as a result of the increasing localization volume. Table 1 also shows the value of the 68% HDI from the posterior divided by the 68% CI from the flat prior, labeled as $\sigma_{H_0}/\sigma_{\text{prior}}$.

<footnote n="55">55 Note that the 68% HDI for these multi-modal distributions is composed of two disjoint intervals. In this case, we force the interval to be around the mode in one contiguous interval. The 68% equal-tailed credible interval (i.e., the 16th and 84th quantiles) has a similar width, but we quote the HDI for consistency with previous results.</footnote>
<table><caption>**Table 1**
Hubble Constant Estimates from GW170817 and the Dark Sirens Considered in This Work</caption><thead><tr><th>Event</th><th>Prior</th><th>H<sub>0</sub></th><th>σ<sub>H0</sub>/H<sub>0</sub></th><th>σ<sub>H0</sub>/σ<sub>prior</sub></th><th>Reference</th></tr></thead><tbody><tr><td>GW170817—bright</td><td>[20,140]</td><td>68.8<sup>+17</sup><sub>-76</sub></td><td>18%</td><td>31%</td><td>Adapted from Nicolaou et al. (2020)</td></tr><tr><td>GW170817—dark</td><td>[10,220]</td><td>76<sup>+48</sup><sub>-23</sub></td><td>47%</td><td>50%</td><td>Fishbach et al. (2019)</td></tr><tr><td>GW190814</td><td>[20,140]</td><td>78<sup>+57</sup><sub>-13</sub></td><td>45%</td><td>86%</td><td>This work</td></tr><tr><td>GW170814</td><td>[20,140]</td><td>77<sup>+41</sup><sub>-33</sub></td><td>48%</td><td>90%</td><td>Soares-Santos et al. (2019)</td></tr><tr><td>GW170814</td><td>[20,140]</td><td>72<sup>+51</sup><sub>-20</sub></td><td>49%</td><td>87%</td><td>This work</td></tr><tr><td>GW190814+GW170814</td><td>[20,140]</td><td>77<sup>+41</sup><sub>-22</sub></td><td>44%</td><td>77%</td><td>This work</td></tr><tr><td>All</td><td>[20,140]</td><td>72.0<sup>+12</sup><sub>-8.2</sub></td><td>14%</td><td>25%</td><td>This work</td></tr></tbody></table>

**Note.** All $H_0$ values and errors are in km s$^{-1}$ Mpc$^{-1}$; all priors are flat. The uncertainty from the flat prior only is derived by assuming the same $H_0$ maximum found in the analysis. Quoted uncertainties represent 68% HDI around the maximum of the posterior. The "$\sigma_{H0}/\sigma_{prior}$" column shows the 68% CI from the posterior divided by 68% CI of the prior width.

In this work, we have ignored the effect of peculiar velocities on the $H_0$ inference. While peculiar velocities can have a significant impact on the $H_0$ inference for nearby events, such as GW170817 (Mukherjee et al. 2019; Palmese et al. 2019; Howlett & Davis 2020; Nicolaou et al. 2020), this effect is expected to be negligible at the distance of GW190814 at the current level of precision. For example, a pessimistic bias of 300 km s$^{-1}$ would result in a ~1 km s$^{-1}$ Mpc$^{-1}$ bias on $H_0$. Moreover, typical biases and uncertainties on peculiar velocities are not significant compared to typical photo-z uncertainties. It is, however, possible that at the lowest redshifts the peculiar motion is correlated across large portions of the sky and that this effect will become relevant after we combine many events. Thus we will address this issue in a separate work.

In this analysis, we have not taken into account possible systematics introduced by the Gaussian approximation made in the sky map for the distance posterior (which enters in Equation (6)). As posterior samples become available for more dark siren events, we will be able to assess the extent of this issue in the future.

It is interesting to note that both GW events overlapping with the DES footprint analyzed thus far are correlated with large overdensities. In fact, the luminosity distance of GW190814 corresponds to $z \sim 0.06$ (for $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$), and it is clear from the top-right panel of Figure 1 that an overdensity of galaxies is present at that redshift. Similarly, there is an overdensity at the location of GW170814 at $z \sim 0.14$. This suggests that GW events could follow the underlying large scale structure (LSS), possibly with a large bias factor. Future work combining more events with LSS could confirm this hypothesis and shed light on the formation channels of these binaries (e.g., Raccanelli et al. 2016).

## 5. Conclusions

In this paper, we present a statistical standard siren measurement using the gravitational wave event GW190814, detected by LIGO/Virgo on 2019 August 14. The redshift information used comes from the DES DNF photo-z's. Our work shows the advantages of using a full redshift PDF of the potential host galaxies versus a Gaussian approximation. Full redshift PDFs provide a better representation of the $dN/dz$, which, together with the luminosity distance estimates, is the ultimate ingredient needed in the dark standard siren framework. We also correct for photo-z biases; however, these are subdominant at the current level of precision. Peculiar motions

are also not expected to provide a significant contribution to our results at the distance of GW190814. We run our analysis both with the preliminary GW190814 sky map released shortly after the GW alert, and with the sky map from the refined LVC analysis that uses higher-order spherical harmonic modes from the GW signal. Because the latter provides better constraints on the binary parameters, we choose it as our final result.

For GW190814, the maximum a posteriori value of $H_0$ with its 68% CI is $H_0 = 78_{-13}^{+57}$ km s$^{-1}$ Mpc$^{-1}$, while it is $H_0 = 77_{-22}^{+41}$ km s$^{-1}$ Mpc$^{-1}$ when combined with our re-analysis of GW170814, using a flat prior in the range [20,140] km s$^{-1}$ Mpc$^{-1}$. Finally, a combination of GW190814, GW170814, and GW170817 yields $H_0 = 72.0_{-8.2}^{+12}$ km s$^{-1}$ Mpc$^{-1}$. The addition of GW190814 and GW170814 to GW170817 improves the 68% CI interval by ~18%, showing how well-localized mergers without counterparts can provide a substantial contribution to standard siren measurements, provided that a complete galaxy catalog is available at the location of the event.

In the future, a combination of well-localized GW events with ongoing and upcoming galaxy surveys, such as DESI, 4MOST, and LSST, will be able to provide competitive constraints on the Hubble constant using the dark standard siren method. We remind the reader that this method is not expected to provide constraints that are competitive with current precision measurements of $H_0$ from single events. In fact, the Hubble constant posteriors for single events can present multiple peaks or be asymmetric distributions, as is the case in this work and in previous dark siren analyses. These effects are due to the presence of overdensities or underdensities along the line of sight at the position of the events, and those will differ from event to event. However, after enough events are combined, the analysis will converge to the true value of $H_0$ (e.g., Gray et al. 2020): the background or foreground overdensities will not be consistently present for the same values of $H_0$, and their contribution to the final posterior will eventually be washed out. If we assume that the statistical uncertainty on $H_0$ scales as $1/\sqrt{N}$, where N is the number of events considered, a combination of $\mathcal{O}(100)$ events with a localization similar to GW170814 and GW190814 will provide a $\mathcal{O}(1)\%$ statistical precision on the Hubble constant. We stress that a range of systematic effects and biases that have been neglected so far will become important at that level of precision (e.g., assumptions about the background cosmology; Keeley et al. 2019; Shafieloo et al. 2020) and will need to be addressed in the coming years. An event sample of this size thus will be
available a few years after LIGO/Virgo run at design sensitivity (KAGRA Collaboration et al. 2018).

Funding for the DES Projects has been provided by the DOE and NSF (USA), MEC/MICINN/MINECO (Spain), STFC (UK), HEFCE (UK), NCSA (UIUC), KICP (University of Chicago), CCAPP (Ohio State), MIFPA (Texas A&M), CNPQ, FAPERJ, FINEP (Brazil), DFG (Germany), and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne Lab, University of California Santa Cruz, University of Cambridge, CIEMAT-Madrid, University of Chicago, University College London, DES-Brazil Consortium, University of Edinburgh, ETH Zürich, Fermilab, University of Illinois, ICE (IEEC-CSIC), IFAE Barcelona, Lawrence Berkeley Lab, LMU München and the associated Excellence Cluster Universe, University of Michigan, NOAO, University of Nottingham, Ohio State University, University of Pennsylvania, University of Portsmouth, SLAC National Lab, Stanford University, University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES Data Management System is supported by the NSF under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-88861, FPA2015-68048, and Centro de Excelencia SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509. Research leading to these results has received funding from the ERC under the EU's 7th Framework Programme, including grants ERC 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020.

This manuscript has been authored by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the US Department of Energy, Office of Science, Office of High Energy Physics. The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes.

A.P. thanks Zoheyr Doctor, Maya Fishbach, and Will Farr for very useful discussion on standard sirens. In this work, we made extensive use of TOPCAT (Taylor 2005).

## ORCID iDs

A. Palmese https://orcid.org/0000-0002-6011-0530

J. Annis https://orcid.org/0000-0002-0609-3987

M. Soares-Santos https://orcid.org/0000-0001-6082-8529

J. Gschwend https://orcid.org/0000-0003-3023-8362

D. E. Holz https://orcid.org/0000-0002-0175-5064

R. Kessler https://orcid.org/0000-0003-3221-0419

R. Morgan https://orcid.org/0000-0002-7016-5471

C. Conselice https://orcid.org/0000-0003-1949-7638

S. Avila https://orcid.org/0000-0001-5043-3662

D. Brooks https://orcid.org/0000-0002-8458-5047

S. Desai https://orcid.org/0000-0002-0466-3288

E. Gaztanaga https://orcid.org/0000-0001-9632-0815

D. W. Gerdes https://orcid.org/0000-0001-6942-2736

D. Gruen https://orcid.org/0000-0003-3270-7644

R. A. Gruendl https://orcid.org/0000-0002-4588-6517

S. R. Hinton https://orcid.org/0000-0003-2071-9349

D. L. Hollowood https://orcid.org/0000-0002-9369-4157

K. Kuehn https://orcid.org/0000-0003-0120-0808

J. L. Marshall https://orcid.org/0000-0003-0710-9474

P. Melchior https://orcid.org/0000-0002-8873-5065

F. Menanteau https://orcid.org/0000-0002-1372-2534

R. Miquel https://orcid.org/0000-0002-6610-4836

A. Roodman https://orcid.org/0000-0001-5326-3486

G. Tarle https://orcid.org/0000-0003-1704-0781

A. R. Walker https://orcid.org/0000-0002-7123-8943

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, PhRvL, 116, 061102

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, PhRvL, 119, 161101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Natur, 551, 85

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, PhRvL, 119, 141101

Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJL, 896, L44

Ackley, K., Amati, L., Barbieri, C., et al. 2020, arXiv:2002.01950

Allen, D. A. 1976, MNRAS, 174, 29P

Andreoni, I., Goldstein, D. A., Kasliwal, M. M., et al. 2020, ApJ, 890, 131

Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017, Natur, 551, 64

Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540

Benítez, N. 2000, ApJ, 536, 571

Benjamin, J., Van Waerbeke, L., Heymans, C., et al. 2013, MNRAS, 431, 1547

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buchs, R., Davis, C., Gruen, D., et al. 2019, MNRAS, 489, 820

Chabrier, G. 2003, PASP, 115, 763

Chen, H.-Y., Fishbach, M., & Holz, D. E. 2018, Natur, 562, 545

Colless, M., Dalton, G., Maddox, S., et al. 2001, MNRAS, 328, 1039

Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Sci, 358, 1556

Cowperthwaite, P. S., Berger, E., Soares-Santos, M., et al. 2016, ApJL, 826, L29

Crocce, M., Castander, F. J., Gaztañaga, E., Fosalba, P., & Carretero, J. 2015, MNRAS, 453, 1513

Cunha, C. E., Lima, M., Oyaizu, H., Frieman, J., & Lin, H. 2009, MNRAS, 396, 2379

Dálya, G., Galgóczi, G., Dobos, L., et al. 2018, MNRAS, 479, 2374

DES Collaboration, Abbott, T. M. C., Abdalla, F. B., et al. 2018, ApJS, 239, 18

De Vicente, J., Sánchez, E., & Sevilla-Noarbe, I. 2016, MNRAS, 459, 3078

Del Pozzo, W. 2012, PhRvD, 86, 043011

Doctor, Z., Kessler, R., Herner, K., et al. 2019, ApJL, 873, L24

Etherington, J., Thomas, D., Maraston, C., et al. 2017, MNRAS, 466, 228

Feeney, S. M., Peiris, H. V., Williamson, A. R., et al. 2019, PhRvL, 122, 061105

Fishbach, M., Gray, R., Magaña Hernandez, I., et al. 2019, ApJL, 871, L13

Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, AJ, 150, 150

Fosalba, P., Crocce, M., Gaztañaga, E., & Castander, F. J. 2015, MNRAS, 448, 2987

Freedman, W. L. 2017, NatAs, 1, 0169

Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34

Garcia, A., Morgan, R., Herner, K., et al. 2020, arXiv:2007.00050

Gomez, S., Hosseinzadeh, G., Cowperthwaite, P. S., et al. 2019, ApJL, 884, L55

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Gray, R., Hernandez, I. M., Qi, H., et al. 2020, PhRvD, 101, 122001

Gschwend, J., Rossel, A. C., Ogando, R. L. C., et al. 2018, A&C, 25, 58

Hartley, W. G., Almaini, O., Mortlock, A., et al. 2013, MNRAS, 431, 3045

Holz, D. E., & Hughes, S. A. 2005, ApJ, 629, 15

Howlett, C., & Davis, T. M. 2020, MNRAS, 492, 3803

Hoyle, B., Gruen, D., Bernstein, G. M., et al. 2018, MNRAS, 478, 592

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Jones, D. H., Read, M. A., Saunders, W., et al. 2009, MNRAS, 399, 683

KAGRA Collaboration, LIGO Scientific Collaboration & Virgo Collaboration 2018, LRR, 21, 3

Keeley, R. E., Shafieloo, A., L'Huillier, B., & Linder, E. V. 2019, MNRAS, 491, 3983
LIGO Scientific Collaboration, Virgo Collaboration, GBM, F., et al. 2017, ApJL, 848, L12

LIGO Scientific Collaboration & Virgo Collaboration 2019a, GCN, 25333, 1

LIGO Scientific Collaboration & Virgo Collaboration 2019b, A gravitational-wave measurement of the Hubble constant following the second observing run of Advanced LIGO and Virgo, arXiv:1908.06060

Lima, M., Cunha, C. E., Oyaizu, H., et al. 2008, MNRAS, 390, 118

Lipunov, V. M., Gorbovskoy, E., Kornilov, V. G., et al. 2017, ApJL, 850, L1

MacLeod, C. L., & Hogan, C. J. 2008, PhRvD, 77, 043512

Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086

Morgan, R., Soares-Santos, M., Annis, J., et al. 2020, Constraints on the Physical Properties of S190814bv through Simulations based on DECam Follow-up Observations by the Dark Energy Survey, arXiv:2006.07385

Morganson, E., Gruendl, R. A., Menanteau, F., et al. 2018, PASP, 130, 074501

Mortlock, D. J., Feeney, S. M., Peiris, H. V., Williamson, A. R., & Nissanke, S. M. 2019, PhRvD, 100, 103523

Mörtsell, E., & Dhawan, S. 2018, JCAP, 9, 025

Mukherjee, S., Lavaux, G., Bouchet, F. R., et al. 2019, arXiv:1909.08627

Nair, R., Bose, S., & Saini, T. D. 2018, PhRvD, 98, 023502

Nicolaou, C., Lahav, O., Lemos, P., Hartley, W., & Braden, J. 2020, MNRAS, 495, 90

Nishizawa, A. 2017, PhRvD, 96, 101303

Nissanke, S., Holz, D. E., Dalal, N., et al. 2013, arXiv:1307.2638

Nissanke, S., Holz, D. E., Hughes, S. A., Dalal, N., & Sievers, J. L. 2010, ApJ, 725, 496

Oyaizu, H., Lima, M., Cunha, C. E., Lin, H., & Frieman, J. 2008, ApJ, 689, 709

Palmese, A., Annis, J., Burgad, J., et al. 2020, MNRAS, 493, 4591

Palmese, A., Graur, O., Annis, J. T., et al. 2019, BAAS, 51, 310

Palmese, A., Hartley, W., Tarsitano, F., et al. 2017, ApJL, 849, L34

Palmese, A., & Kim, A. G. 2020, arXiv:2005.04325

Palmese, A., Lahav, O., Banerji, M., et al. 2016, MNRAS, 463, 1486

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv:1807.06209

Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13

Raccanelli, A., Kovetz, E. D., Bird, S., Cholis, I., & Muñoz, J. B. 2016, PhRvD, 94, 023516

Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85

Schmidt, S. J., Malz, A. I., Soo, J. Y. H., et al. 2020, Evaluation of probabilistic photometric redshift estimation approaches for LSST, arXiv:2001.03621

Schutz, B. F. 1986, Natur, 323, 310

Shafieloo, A., Keeley, R. E., & Linder, E. V. 2020, JCAP, 2020, 019

Singer, L. P., Chen, H.-Y., Holz, D. E., et al. 2016, ApJS, 226, 10

Soares-Santos, M., Holz, D. E., Annis, J., et al. 2017, ApJL, 848, L16

Soares-Santos, M., Kessler, R., Berger, E., et al. 2016, ApJL, 823, L33

Soares-Santos, M., Palmese, A., Hartley, W., et al. 2019, ApJL, 876, L7

Tanvir, N. R., Levan, A. J., González-Fernández, C., et al. 2017, ApJL, 848, L27

Taylor, M. B. 2005, in ASP Conf. Ser. 347, TOPCAT & STIL: Starlink Table/ VOTable Processing Software, ed. P. Shopbell, M. Britton, & R. Ebert (San Francisco, CA: ASP), 29

The Dark Energy Survey Collaboration 2005, in ASP Conf. Ser. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert (San Francisco, CA: ASP), 29

The Dark Energy Survey Collaboration 2016, MNRAS, 460, 1270

Valenti, S., David, S. J., et al. 2017, ApJL, 848, L24

Verde, L., Treu, T., & Riess, A. G. 2019, NatAs, 3, 891

Vieira, N., Ruan, J. J., Haggard, D., et al. 2020, ApJ, 895, 96

Vitale, S., & Chen, H.-Y. 2018, PhRvL, 121, 021303

Watson, A. M., Butler, N. R., Lee, W. H., et al. 2020, MNRAS, 492, 5916

Weigel, A. K., Schawinski, K., & Bruderer, C. 2016, MNRAS, 459, 2150

Yu, J., Wang, Y., Zhao, W., & Lu, Y. 2020, Hunting for the Host Galaxy Groups of Binary Black Holes and the Application in Constraining Hubble Constant, arXiv:2003.06586