First measurement of the Hubble constant from a dark standard siren using the Dark Energy Survey galaxies and the LIGO/Virgo binary-black-hole merger GW170814

TL;DR

This paper reports the first measurement of the Hubble constant using a dark standard siren from a binary black hole merger, combining gravitational wave data with galaxy surveys to provide a novel, statistical approach to cosmology.

Contribution

It introduces a new method for estimating the Hubble constant using black-hole mergers without electromagnetic counterparts, leveraging galaxy catalogs for statistical inference.

Findings

Measured H_0 = 75.2^{+39.5}_{-32.4} km/s/Mpc with a single event.

Results are consistent with SN Ia and CMB measurements.

Demonstrates potential for improved constraints with more events.

Abstract

We present a multi-messenger measurement of the Hubble constant H_0 using the binary-black-hole merger GW170814 as a standard siren, combined with a photometric redshift catalog from the Dark Energy Survey (DES). The luminosity distance is obtained from the gravitational wave signal detected by the LIGO/Virgo Collaboration (LVC) on 2017 August 14, and the redshift information is provided by the DES Year 3 data. Black-hole mergers such as GW170814 are expected to lack bright electromagnetic emission to uniquely identify their host galaxies and build an object-by-object Hubble diagram. However, they are suitable for a statistical measurement, provided that a galaxy catalog of adequate depth and redshift completion is available. Here we present the first Hubble parameter measurement using a black-hole merger. Our analysis results in $H_0 = 75.2^{+39.5}_{-32.4}~{\rm km~s^{-1}~Mpc^{-1}}$, which is consistent with both SN Ia and CMB measurements of the Hubble constant. The quoted 68% credible region comprises 60% of the uniform prior range [20,140] ${\rm km~s^{-1}~Mpc^{-1}}$, and it depends on the assumed prior range. If we take a broader prior of [10,220] ${\rm km~s^{-1}~Mpc^{-1}}$, we find $H_0 = 78^{+ 96}_{-24}~{\rm km~s^{-1}~Mpc^{-1}}$ ($57\%$ of the prior range). Although a weak constraint on the Hubble constant from a single event is expected using the dark siren method, a multifold increase in the LVC event rate is anticipated in the coming years and combinations of many sirens will lead to improved constraints on $H_0$.

Tables

Table 1. Table 1: Hubble constant estimate from this work. All H 0 subscript 𝐻 0 H_{0} values and errors are in km ​ s − 1 ​ Mpc − 1 km superscript s 1 superscript Mpc 1 ~{}{\rm km~{}s^{-1}~{}Mpc^{-1}} . The uncertainty from the flat prior only is derived by assuming the same H 0 subscript 𝐻 0 H_{0} maximum found in the analysis. Quoted uncertainties represent 68% confidence level around the maximum of the posterior, and they are statistical only. The last column quantity ( σ H 0 / H 0 subscript 𝜎 subscript 𝐻 0 subscript 𝐻 0 \sigma_{H_{0}}/H_{0} prior) corresponds to 68% times the prior width divided by H 0 subscript 𝐻 0 H_{0} .

Prior	$H_0$	$+{\sigma }_{H_0}$	$-{\sigma }_{H_0}$	${\sigma }_{H_0}/H_0$	${\sigma }_{H_0}/H_0$ prior
$[20,140]$	75	40	32	47.8%	54.3%

Equations

$$p(H_0|d_{\mathrm{GW}},d_{\mathrm{EM}})\propto p(d_{\mathrm{GW}},d_{\mathrm{EM}}|H_0)p(H_0).$$

$$\begin{array}{r}p(d_{\mathrm{GW}},d_{\mathrm{EM}}|\{z_j,{\hat{\Omega }}_j\},H_0)\propto \sum _i{w}_i\int \mathrm{d}{d}_L\mathrm{d}{\hat{\Omega }}_{\mathrm{GW}}p(d_{\mathrm{GW}}|d_L,{\hat{\Omega }}_{\mathrm{GW}})\\ \times p(d_{\mathrm{EM}}|\{z_j,{\hat{\Omega }}_j\}){\delta }_D(d_L-d_L(z_i,H_0)){\delta }_D({\hat{\Omega }}_{\mathrm{GW}}-{\hat{\Omega }}_i),\end{array}$$

$$\begin{array}{rl}p(z_i,{\hat{\Omega }}_i)\mathrm{d}{z}_i\mathrm{d}{\hat{\Omega }}_i& \propto \frac{1}{V_{\max }}\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{z}_i\mathrm{d}{\hat{\Omega }}_i}\mathrm{d}{z}_i\mathrm{d}{\hat{\Omega }}_i\propto \frac{1}{V_{\max }}\frac{r^2(z_i)}{H(z_i)}\mathrm{d}{z}_i\mathrm{d}{\hat{\Omega }}_i,\end{array}$$

$$p(d_{\mathrm{EM}}|\{z_j\})=\prod _{i}p(z_{\mathrm{obs},i}|z_i)=\prod _i\mathcal{N}(z_{\mathrm{obs},i},{\sigma }_{z,i};z_i).$$

$$p(d_{\rm GW}|d_{L},\hat{\Omega})\propto p(\hat{\Omega})\frac{1}{\sqrt{2\pi}\sigma(\hat{\Omega})}{\rm exp}\Big{[}-\frac{\big{(}d_{L}-\mu(\hat{\Omega})\big{)}^{2}}{2\sigma^{2}(\hat{\Omega})}\Big{]}N(\hat{\Omega})\,,$$

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

$$p(H_0|d_{\mathrm{GW}},d_{\mathrm{EM}})\propto \frac{p(H_0)}{V[d_{L,\mathrm{GW}}^{\max }(H_0)]}\sum _i\frac{1}{{\mathcal{Z}}_i}\int \mathrm{d}{z}_{i}p(d_{\mathrm{GW}}|d_L(z_i,H_0),{\hat{\Omega }}_i)p(d_{\mathrm{EM}}|z_i)\frac{r^2(z_i)}{H(z_i)},$$

$$\begin{array}{r}p(H_0|\{d_{\mathrm{GW},j}\},d_{\mathrm{EM}})\propto p(H_0)\int {\mathrm{d}}^N{z}_k{\mathrm{d}}^N{\hat{\Omega }}_{k}p(z_k,{\hat{\Omega }}_k)\times \\ p(d_{EM}|\{z_k,{\hat{\Omega }}_k\})\left[\prod _{j}p(d_{GW,j}|\{z_k,{\hat{\Omega }}_k\})\right].\end{array}$$

Full text

catalogs — cosmology: observations — gravitational waves — surveys

\AuthorCollaborationLimit

=3000

\reportnum

DES-2018-0389 \reportnumFERMILAB-PUB-18-629-AE

1 Introduction

Unlike most extragalactic distance observables, mergers of neutron star and black hole binary systems are absolute distance indicators. Often referred to as “standard sirens”, they emit gravitational waves (GW) from which the luminosity distance can be inferred without relying on any calibration with respect to another source: the rate of change in frequency gives the system’s size and thus the intrinsic amplitude, which is compared against the observed signal amplitude to obtain the distance to the source. If redshifts are associated with those sirens (in the simplest case, the host galaxy is identified and its redshift is obtained via spectroscopic follow up), a measurement of the present rate of expansion of the Universe $H_{0}$ can be achieved via the distance–redshift relation. The use of gravitational wave sources as cosmological probes was first proposed by Schutz (1986), and recently revisited in several works (e.g. Holz & Hughes 2005).

For dark energy research, the possibility of measuring $H_{0}$ directly and independently from other methods is of great interest. Local measurements obtained from type Ia Supernovae (SN Ia) and other distance indicators, as well as the predicted value inferred from the cosmic microwave background at $z\sim 1100$, have achieved remarkable precision of $1-2.5\%$ (e.g. Riess et al. 2018; Planck Collaboration et al. 2018). They disagree, however, by more than $3\sigma$ and interpreting this tension as evidence for beyond-$\Lambda$CDM dark energy or new physics at the early universe requires new measurements of great precision and accuracy (Freedman 2017; Mörtsell & Dhawan 2018). Those measurements are one of the greatest challenges faced by current experiments in cosmology because the observables are subject to correlated systematic effects arising from their complex astrophysics. As estimates become more precise, this challenge becomes more severe and the need for novel independent methods becomes more pressing. Those methods, however, are few and hard to come by. One possibility is standard sirens, which remained elusive for almost 30 years, until the detection of the first gravitational wave event (GW150914; Abbott et al. 2016). The first standard siren-based $H_{0}$ measurement (Abbott et al. 2017a) came with the discovery of the binary–neutron–star (BNS) merger GW170817 (Abbott et al. 2017) and its associated electromagnetic counterpart (LIGO Scientific Collaboration et al. 2017; Soares-Santos et al. 2017; Arcavi et al. 2017; Coulter et al. 2017; Lipunov et al. 2017; Tanvir et al. 2017; Valenti et al. 2017). Several studies have developed methodologies to infer cosmological parameters from standard sirens and establish their constraining power (Schutz, 1986; Holz & Hughes, 2005; MacLeod & Hogan, 2008; Nissanke et al., 2010; Del Pozzo, 2012; Nissanke et al., 2013; Nishizawa, 2017; Chen et al., 2018; Feeney et al., 2018; Vitale & Chen, 2018; Mortlock et al., 2018). Chen et al. (2018) predict that we will be able to constrain $H_{0}$ with 2% precision within 5 years with standard sirens detected by LIGO/Virgo, while Nair et al. (2018) predict a $\sim 7\%$ measurement with just 25 binary black hole (BBH) events from the Einstein telescope.

Anticipating that the LIGO/Virgo Collaboration (LVC) network of gravitational wave detectors would eventually achieve sensitivity sufficient to enable standard siren–based measurements, the Dark Energy Survey (DES) collaboration and external collaborators launched in 2015 the DES gravitational waves (DESGW) program. DESGW uses DECam to search for optical emission associated with LVC detected mergers and pursues cosmological measurements with standard sirens. In particular, the multi-messenger shared discovery of the neutron–star merger GW170817 and of its optical kilonova, resulted in a measurement of $H_{0}$ (Abbott et al. 2017a) that inaugurated the era of siren-based cosmology. We have also performed the most comprehensive searches for optical emission to black hole events, including GW150914 (Soares-Santos et al. 2016), GW151226 (Cowperthwaite et al. 2016), and GW170814 (Doctor et al. 2018). These events are expected to be dark, although the possibility of optical emission has yet to be observationally excluded.

Dark sirens can also be used for cosmology using a statistical method, as first proposed in Schutz (1986). Provided a catalog of potential host galaxies within the event localization region, their redshifts will contribute in a probabilistic way to the measurement of $H_{0}$, depending on the galaxies’ distance and sky position. This approach has been developed within a Bayesian framework by Del Pozzo (2012) and Chen et al. (2018) and implemented in Fishbach et al. (2018) using GW170817, which produced results consistent with the first measurement (Abbott et al. 2017a) where the identified host galaxy, NGC 4993 (e.g., Palmese et al. 2017), was used. Eventually, a large sample of events will enable precise cosmological measurements using the dark siren approach.

In this work, we measure $H_{0}$ using the gravitational wave event GW170814 (Abbott et al. 2017b) as a dark siren. GW170814 resulted from the inspiral and merger of a binary black hole system at a luminosity distance of $540^{+130}_{-210}$ Mpc (median value with 90% credible interval). The masses of the black holes were $30.5^{+5.7}_{-3.0}$ and $25.3^{+2.8}_{-4.2}~{}M_{\odot}$, each. GW170814 is the first BBH detected by a triple network (including LIGO Hanford and Livingston, plus Virgo), and it has the smallest localization volume of any of the BBH events detected by LVC thus far. Therefore the number of potential host galaxies is lower compared to other events, making GW170814 the most appropriate event for this measurement. Additionally, the event localization region falls within the DES footprint, making DES galaxy catalogs a prime sample for measurement of $H_{0}$. With this one event, our goal is to provide a proof of principle measurement, addressing the challenges that are specific to the dark siren method, and establishing its potential to yield precision cosmology results in the near future.

A key component of the measurement is crafting the appropriate galaxy catalog: completeness, as well as precise and accurate photometric redshifts (photo–$z$’s), throughout the entire volume probed are required. The overlap of GW170814’s area with DES allows us to employ galaxy catalogs produced from the first three years of the survey (DES Y3; Abbott et al. 2018). This first dark siren measurement is a step towards incorporating this new cosmological probe into the portfolio of cosmic surveys for dark energy.

A detailed description of the data used in this analysis is provided in §2, followed by a description of our implementation of the method in §3. We present our results and discussion in §4, and our conclusions in §5. Throughout this paper we assume a flat $\Lambda$CDM cosmology with $\Omega_{m}=0.3$ and $H_{0}$ values in the $20-140~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ range. All quoted error bars represent the $68\%$ confidence level (CL), unless otherwise stated.

2 Data

2.1 The LVC sky map

The sky map used in this work is the publicly available LALInference map (LIGO Scientific Collaboration & Virgo Collaboration, 2017)111https://dcc.ligo.org/LIGO-T1700453/public, provided in HEALPix (Górski et al. 2005) pixels. The luminosity distance probability distribution is approximated with a Gaussian in each pixel. The region of interest, enclosing 90% of the localization probability, is 61.66 deg2. The projected sky map and the distribution of luminosity distance mean values from the LVC distance likelihood in each pixel within the region of interest are shown in Figure 1. The probability peak is located at RA, Dec = $(47.523,-44.856)$ deg. At the peak location, the luminosity distance is 504.7 Mpc and the Gaussian width is 91.9 Mpc. Using the limiting values of our $H_{0}$ prior range ([20,140] $~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$) we can convert the 90% and $99.7\%$ distance range into a redshift range ($0.02<z<0.26$ and $z<0.3$, respectively) for this analysis.

2.2 The DES galaxy catalog

The DES222www.darkenergysurvey.org(The Dark Energy Survey Collaboration 2005; Dark Energy Survey Collaboration et al. 2016) is an optical-near-infrared survey that images 5000 ${\rm deg}^{2}$ of the South Galactic Cap in the $grizY$ bands. The survey is being carried out using a $\sim 3$ $\textrm{deg}^{2}$ CCD camera (the DECam, see Flaugher et al. 2015) mounted on the Blanco 4-m telescope at the Cerro Tololo Inter-American Observatory (CTIO) in Chile. The data used here are from the first 3 years of observations (September 2013 – February 2016, Abbott et al. 2018).

The DES Data Management (DESDM) pipeline was used for data reduction (Morganson et al. 2018).The process includes calibration of the single-epoch images, which are co–added after background subtraction and then cut into tiles. The source catalogue was created using Source Extractor (SExtractor, Bertin & Arnouts 1996) to detect objects on the $riz$ co-added images. The median $10\sigma$ limiting magnitudes of Y3 data for galaxies are $g=24.33$, $r=24.08$, $i=23.44$, $z=22.69$, and $Y=21.44$ mag (Abbott et al. 2018). The photometry used in this work is part of a value–added Y3 catalog not released with DR1, and is the result of the Multi-Object Fitting (MOF) pipeline that uses the ngmix code.333https://github.com/esheldon/ngmix Following a procedure similar to Drlica-Wagner et al. (2017) for Year 1 data, the DES collaboration made further selections to produce a high-quality object catalog called the Y3 “gold” catalog. For this sample, redshifts have been computed using the Directional Neighborhood Fitting (DNF; De Vicente et al. 2016), and they are not included in DR1.

The DNF method applied to Y3 data provides redshift information for each galaxy in the form of a probability distribution function (PDF), from which a mean redshift, and half of the central 68th percentile width are computed. The width of the PDF can be over or under-estimated due to the sampling of the training set and algorithmic details of DNF. This issue is particularly relevant for the redshift range used in this work, which is low compared to that exploited in weak lensing and large scale structure cosmology, for which the DNF method was optimized. We find that the typical uncertainty below redshift $z\sim 0.1$ is underestimated by a factor of 10 when compared to the typical scatter found for the subset of the galaxies with available spectroscopic redshifts (where the standard deviation is $\sigma\gtrsim 0.015$). Thus, we add a minimum uncertainty of $0.015$ for these low–$z$ galaxies. At $0.1<z<0.3$, the uncertainty is well behaved and the average value follows $\bar{\sigma_{z}}(z)\simeq 0.013~{}(1+z)^{3}$, as we find using an empirical fit.

We produce alternative photo–$z$ estimates with another machine learning code, ANNz2 (Sadeh et al. 2016). This allows us to test the impact of the correction applied to the DNF errors on the posterior of the Hubble constant. Photo–$z$ with ANNz2 have previously been validated for cosmological analyses using DES Science Verification data (Bonnett et al. 2016; Leistedt et al. 2016; Abbott et al. 2016) and for the Kilo–Degree Survey (KiDS; Bilicki et al. 2018), and are produced as part of the DES photo–$z$ pipeline (Gschwend et al. 2018). In particular, it provides error estimates through a $k$–nearest neighbor ($k$NN) method, and dedicated redshifts for the purposes of this analysis. We additionally employ a reweighting technique (Lima et al. 2008) specifically for our galaxy sample to further tune our redshifts. We run ANNz2 in randomized regression mode with 50 Boosted Decision Trees (BDTs), using a spectroscopic sample of $245,458$ matching Y3 galaxies out to redshift $z\gtrsim 1$, randomly split into subsamples for training, testing and validation. The training and the reweighting use $griz$ MOF magnitudes. We find that the typical error roughly follows $\sim 0.02~{}(1+z)^{3}$ in the redshift range of interest. The two algorithms, DNF and ANNz2, gave similar results, see section §4.

These redshifts, together with publicly available spectroscopic redshifts from 2dF, 6dF and SPT–GMOS (Colless et al., 2001; Jones et al., 2009; Bayliss et al., 2016) and the DES MOF photometry, are used to estimate galaxy properties (including stellar mass and absolute magnitude) of this sample. This is achieved 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. 2019). We add a 0.05 systematic uncertainty in quadrature to the magnitudes, to account for systematic uncertainties in magnitude estimation and model variance.444This is a regularization to compensate for the synthetic model set grid and the fact that many SED fitting codes do not include a model error function. The value chosen is based on past experience of what gives stable results. The simple stellar population (SSP) templates used for the fitting are Bruzual & Charlot (2003), with three metallicities ($0.2~{}Z_{\odot}$, $Z_{\odot}$ and $2.5~{}Z_{\odot}$), 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 as ${\rm e}^{-t/\tau}$ with $\tau=0.1,0.3,1,2,3,5,10,15$ and $30$ Gyr.

The source list of the Y3 gold catalogue is $95\%$ complete for galaxies within our apparent magnitude limit, $r<23.35$ (Abbott et al. 2018). This value is computed through the recovery rate of sources from the deeper CFHTLenS survey (Erben et al., 2013), and thus includes the correct distribution of surface brightnesses. Nevertheless, extended, low surface brightness galaxies near our flux limit may be preferentially missed by the detection pipeline. We therefore provide an approximate completeness of sources throughout the redshift range of interest. Using DNF mean redshifts we convert the source completeness to $r<23.35$ from Abbott et al. (2018) (Figure 12) into a completeness in redshift intervals, $\Delta z=0.02$. By taking the peak of the magnitude distribution in each bin as roughly our observed magnitude limit at that redshift, we find our sample is $>93\%$ complete across the range $0<z<0.26$. We further determined that the fraction of low redshift, extended galaxies missed by the DES Y3 pipeline is $\sim 1\%$, when compared with the 2MASS extended source catalog (Huchra et al., 2012). For the purpose of this paper, we choose to ignore those ultra-low $z$ sources as most of them are at $z<0.02$ and are not relevant for the present analysis.

The DES Y3 gold catalog is nonetheless an observed magnitude–limited sample. This analysis requires a volume-limited sample, which we obtain by applying a luminosity cut. In order to determine the appropriate cut to create a volume–limited sample, we compute the completeness limits in terms of absolute quantities (luminositiy and stellar mass). We follow the method outlined in Pozzetti et al. (2010) and Hartley et al. (2013). We identify galaxies with observed magnitudes that are bright enough to be complete and representative of the real galaxy population within redshift bins. To compute the $95\%$ completeness limit in (rest–frame) luminosity, we scale the luminosities of this sample to that which they would have if their observed magnitude were equal to the survey completeness limit, and take the 95th percentile of the resulting luminosity distribution. This value corresponds to $-17.2$ in $r$-band absolute magnitude and $\sim 3.8\times 10^{8}~{}{\rm M}_{\odot}$ in stellar mass for the redshift range of interest. We cut the DES catalog at the specific absolute luminosity value mentioned above. We conclude that our volume–limited galaxy sample is complete within the redshift range of interest for galaxies down to stellar masses of $\sim 3.8\times 10^{8}~{}{\rm M}_{\odot}$. In other words, our galaxy catalog contains $\sim 77\%$ of the total stellar mass in the volume considered by assuming that the galaxies follow a Schechter stellar mass function with the best fit values from Weigel et al. (2016).

The final galaxy stellar mass and redshift distributions of galaxies are shown in Figure 1. The stellar mass map clearly shows the presence of large scale structure, including clusters, voids and filaments. We recognize a number of well–known clusters within the volume of interest, including several Abell clusters. A uniform distribution of galaxies in comoving volume $({\rm d}N/{\rm d}z)_{\rm com}$ has been subtracted from the observed galaxies’ redshift distribution in Figure 1 to highlight the overdensities. The $({\rm d}N/{\rm d}z)_{\rm com}$ distribution has been obtained by assuming $H_{0}=70~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ and it contains the same total number of galaxies as the observed ${\rm d}N/{\rm d}z$ over the redshift range shown. We are able to identify a “wall”–like structure around $z\sim 0.06$ spanning most of the area between $35<{\rm RA}<55$ and $-55<{\rm Dec}<-35$, which is spectroscopically confirmed by 2dF, LCRS (Shectman et al., 1996), and especially 6dF. A broader galaxy overdensity is found around $z\sim 0.12$ (also seen in LCRS and 2dF, and composed of several Abell galaxy clusters). This broad peak is also identified in redshift distributions by other photo–$z$ codes, including a template based code, the Bayesian Photometric Redshift (BPZ; Benítez 2000). We have further verified that overdensities at the lowest redshifts ($z\sim 0.06$) are also present in spectroscopic samples outside of the region of interest. This is expected at these low redshifts, where large scale structure projects onto vast areas of the sky. In summary, there are $77,092$ galaxies within the $90\%$ LIGO/Virgo probability volume, and $105,011$ when $99.7\%$ of the distance probability is considered, of which $\sim 6,000$ have spectroscopic redshifts.

3 Method

In order to estimate the posterior probability of $H_{0}$ given GW data $d_{GW}$ from a single event detection, and electromagnetic (EM) data from a galaxy survey, we follow Chen et al. (2018). By applying Bayes’ theorem, one can write the posterior as:

[TABLE]

We assume that all cosmological parameters except for $H_{0}$ are fixed (Flat $\Lambda$CDM cosmology with $\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$). We treat the joint GW and EM likelihood $p(d_{\rm GW},d_{\rm EM}|H_{0})$ as the product of two individual likelihoods (since the processes involved in producing the data from the two experiments are independent) marginalized over all variables except for the true luminosity distance $d_{L}$ and solid angle $\hat{\Omega}_{\rm GW}$ of the GW source, and for the true host galaxy redshift $z_{i}$ and solid angle $\hat{\Omega}_{i}$. Note that the solid angles $\hat{\Omega}$ are vectors with the angular position of the source/galaxy as direction, and they all subtend the same area ($\sim 3\times 10^{-3}~{}{\rm deg^{2}}$) as the sky is pixelized with HEALPix maps in this work. If we assume that the event happened in one of the observed galaxies $i$, then $\hat{\Omega}_{\rm GW}$ and $\hat{\Omega}_{i}$ are related, and so are $d_{L}$ and $z_{i}$ through the cosmology (in this case, $H_{0}$). By marginalizing also over the choice of galaxy $i$, the joint, marginal likelihood can be written as:

[TABLE]

where $\delta_{D}$ is the Dirac delta function, $w_{i}$ are weights that represent the relative probability that different galaxies host a GW source, and $\{z_{j},\hat{\Omega}_{j}\}$ represents all the galaxies’ redshift and solid angle. These weights could be based on some galaxy properties, such as luminosity or star–formation rate, but here we assume they are uniform across all galaxies given our lack of knowledge of GW host galaxy properties.

We also need to marginalize over the galaxies’ redshifts and sky positions, with a reasonable choice of prior $p(z_{i},\Omega_{i})$. If one assumes that the galaxies are uniformly distributed in comoving volume $V$, and volume–limited within $V_{\rm max}$:

[TABLE]

where $r$ is the comoving distance to the galaxy. While this assumption holds on average over sufficiently large volumes, it is possible that future precision cosmology analyses will require taking into account the real clustering of galaxies in this formalism.

Assuming that we precisely know the galaxies’ positions $\{\hat{\Omega}_{j}\}$ (which is realistic especially in the limit in which spatial probabilities are considered within HEALPix pixels), we can integrate over the galaxies’ positions as delta functions about the observed values. The marginal EM likelihood reduces to $p(d_{\rm EM}|\{z_{j}\})$, which we approximate for simplicity by a product of Gaussian distributions, ${\mathcal{N}}$, for each galaxy, centred around the observed redshift values $z_{{\rm obs},i}$ with a width given by the redshift’s uncertainty $\sigma_{z,i}$ for each galaxy $i$:

[TABLE]

The marginal GW likelihood $p(d_{\rm GW}|d_{L},\Omega)$ can be computed as prescribed in Singer et al. (2016):

[TABLE]

where the position probability, location, normalization and scale (PROB $p(\hat{\Omega})$, DISTMU $\mu$, DISTNORM $N$ and DISTSTD $\sigma$ respectively) of the luminosity distance at each position are provided in the sky map.

We now consider the selection effects of GW events and galaxies introduced by the experiments’ sensitivities and detection pipelines. We follow the approach of Chen et al. (2018) and Mandel et al. (2018), and include a $[\beta(H_{0})]^{-1}$ factor that normalizes the likelihood over all possible GW and EM data. Given that our galaxy catalog is volume–limited out to larger distances than the maximum observable distance for the GW events, this term reduces to:

[TABLE]

where $V[d_{L,{\rm GW}}^{\rm max}(H_{0})]$ is the maximum observable volume for the GW events considered.

Finally, Eq. (1) becomes: {widetext}

[TABLE]

where $\mathcal{Z}_{i}=\int p(d_{\rm EM}|z_{i})r^{2}(z_{i})/H(z_{i})~{}{\rm d}z_{i}$ are evidence terms that arise from integrating out the other galaxy redshifts in each term of the sum. This formalism can be extended to combine data $\{d_{{\rm GW},j}\}$ and $d_{{\rm EM}}$ from a sample of multiple events $j$, assuming that the GW events are independent and that the galaxy catalog is fixed for all events:

[TABLE]

In the following, we assume a flat prior on $H_{0}$ within [20,140] $~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$, unless otherwise stated. This is a very broad prior, covering a range much larger than current estimates of $H_{0}$. This choice was made as a compromise between the following aspects: i) a result which is mostly informed by the LVC and DES data rather than by external constraints, ii) a result which can be compared with the first standard siren estimate, and iii) a complete galaxy sample which contains most of the stellar mass within the localization volume, to minimize the chance of missing the real host galaxy. As explained in more detail in §4, the redshift cut is related to the $H_{0}$ prior range, and in order to explore higher values of $H_{0}$, one needs to include higher redshift galaxies, and make a higher luminosity cut to preserve the volume limited sample.

A blinded analysis has been performed when estimating the $H_{0}$ posterior from the data to avoid confirmation bias. The values of the Hubble constant have been randomly displaced by an unknown amount, and we unblinded after our pipeline was able to reliably reproduce the input cosmology on simulation tests.

4 Results and Discussion

We apply the described methodology to the DES galaxies’ redshifts and the GW170814 LIGO/Virgo sky map to produce a posterior distribution for the Hubble constant. We find that changes in the $H_{0}$ estimate and its uncertainty between using the corrected DNF photo–$z$’s or the ANNz2 outputs are below the percent level. This agreement is expected, since the two methods produce redshift distributions that are consistent with similar uncertainties. We also add a $0.001$ systematic redshift error in quadrature (corresponding to a typical peculiar velocity of $\sim 300~{}{\rm km~{}s^{-1}}$). The effect of this correction on the posterior is negligible because only a few percent of the galaxies have a spectroscopic redshift, and the effect of peculiar velocities on the remaining galaxies is more than an order of magnitude below their typical photo–$z$ error.

Our maximum a posteriori estimate of the Hubble constant is $H_{0}=75^{+40}_{-32}~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ using a flat prior between 20 and 140 ${\rm km~{}s^{-1}~{}Mpc^{-1}}$. The full posterior distribution is shown in Figure 2, and Table 1 summarizes our findings. The presence of a main, though broad, peak, is expected given the large scale structure seen in the observed volume.

As described in section 2.1, the galaxy sample used in these results is selected as described in §2, and covers the LIGO/Virgo 90% credible localization volume. The distance cut is translated into a redshift cut (made on the mean photo-$z$ value of each galaxy) for a given $H_{0}$ prior. This cut ensures that the galaxy catalog is as complete as possible throughout the whole redshift range of interest for the cosmological parameters used, and includes the fainter galaxies observable for a volume–limited sample defined as in §2. In fact, in order to include more distant galaxies, the luminosity cut needs to be brighter to ensure that the sample is still volume–limited, with the risk of missing the true host galaxy. We have explored the impact of the redshift cut on the $H_{0}$ posterior, while keeping the angular selection to be within the 90% credible localization area. The effect of including galaxies out to $99.7\%$ of the distance localization (corresponding to $z\lesssim 0.3$) is most pronounced at high $H_{0}$ values, as shown by the shaded red region in Figure 2. With this less restrictive cut, the credible region shifts to $H_{0}=77^{+41}_{-33}~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$, showing a $\sim 2\%$ change of the maximum. The effect described here arises from tens of thousand of galaxies at the higher redshifts included with the more relaxed distance cut and the ansatz of Gaussianity of the luminosity distance posterior. In fact, these galaxies contribute with a non–negligible probability to the posterior because of the high $d_{L}$ tail shown in the bottom right panel of Figure 1, and they contribute more significantly at high $H_{0}$ values. This few percent effect is insignificant at the current levels of precision, but will need to be explored in the future using a more realistic luminosity distance posterior.

Our result agrees well (as expected, due to the large uncertainty) with the latest CMB estimate of the Hubble constant by the Planck Collaboration ($67.36\pm 0.54~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ from TT,TE,EE+lowP+lensing; Planck Collaboration et al. 2018), and with results using distance ladder methods by ShoES ($73.52\pm 1.62~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$; Riess et al. 2016) and by DES ($67.77\pm 1.30~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ from SN+BAO; Macaulay et al. 2018).

For the bright standard siren measurement using GW170817 and its electromagnetic counterpart, Abbott et al. (2017a) found $H_{0}=70.0^{+12.0}_{-8.0}~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ at 68% credible interval. Without an EM counterpart leading to a unique host galaxy redshift, we would have recovered a broader $H_{0}$ posterior since we average over all possible host galaxies in the localization volume. For example, Fishbach et al. (2018) applied the statistical standard siren method to GW170817 and found a larger uncertainty than the counterpart standard siren result: $H_{0}=76^{+48}_{-23}~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ for a uniform prior over the range $[10,220]~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$. For a BBH standard siren measurement, as in this work, the combination of the larger localization volume (implying a significantly greater number of potential host galaxies) and the large photometric redshift uncertainty for each galaxy results in an even broader $H_{0}$ posterior. Therefore, while applying the statistical standard siren method to GW170817 yields a 68% credible region on $H_{0}$ comprising 34% of the prior range (Fishbach et al., 2018), in this work we obtain a 68% credible region on $H_{0}$ that is 60% of the prior range. We note that the prior used in Fishbach et al. (2018) is 1.75 times broader than the prior used in this work; if we adopt the same broader prior of [10,220] for our analysis of GW170814, we find $H_{0}=78^{+96}_{-24}~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$. The dependence of the width of the $H_{0}$ posterior on the prior width is a consequence of the fact that the GW observation, which provides only a luminosity distance estimate, is consistent with arbitrarily large $H_{0}$’s, if there are galaxies at sufficiently large redshifts. If the galaxy catalogue extends to some redshift, $z_{\rm max}$, the posterior would fall off around $H_{0}\approx cz_{\rm max}/d_{L}$, where $d_{L}$ is the typical luminosity distance from the GW posterior. However, this fall off is artificial since there are galaxies at greater redshifts which are not included in the catalogue. These may be accounted for using catalogue incompleteness corrections. We chose the prior range for this analysis rather than a larger one such that we did not need to include such corrections, which simplifies the analysis. However, dark siren measurements will become particularly interesting when multiple events can be combined and this effect becomes irrelevant (Chen et al. 2018).

The analysis in Fishbach et al. (2018) for GW170817 used the GLADE galaxy catalog (Dálya et al. 2018), and accounted for incompleteness at the distance of GW170817. GLADE becomes significantly incomplete at the distance to GW170814. As GW detectors improve in sensitivity, the majority of dark standard sirens will be detected at even greater distances and with larger localization volumes, well beyond the reach of spectroscopic galaxy catalogs. This highlights the need for reliable and complete photometric galaxy catalogs. Surveys such as DES, Pan–STARRS1 (Chambers et al. 2016) and LSST are therefore likely to play an important role in future constraints from BBH standard sirens.

The assumption throughout this work is that even if the event occurred in a galaxy below our luminosity threshold, large scale structure predicts that fainter galaxies follow the clustering pattern of the more luminous galaxies in our sample. We have verified in our simulations that a threshold up to 1 magnitude brighter than the limit used here to place events has a negligible impact over a sample of 100 events, provided that the catalog is volume–limited for the range of redshifts relevant to the measurement.

Since galaxies are biased tracers of the Universe’s dark matter, some theories predict that the origin of the black holes involved in these GW events is primordial, constituting part or all of the dark matter (Bird et al. 2016; Clesse & García-Bellido 2017; García-Bellido 2017; Clesse & García-Bellido 2018). In that case, GW events follow exactly the underlying dark matter distribution (presenting an unbiased tracer). Because of the stellar mass to dark matter halo connection (see Wechsler & Tinker 2018 and references therein) it is reasonable to weight galaxies by their stellar mass in Eq. (2) as $w_{i}\propto M_{\star}$. The impact of this scaling with stellar mass or star–formation rate has been explored in Fishbach et al. (2018). We find that the stellar mass weighting has a negligible effect on the posterior. This is due to the large volume analyzed (over which the stellar masses tend to be averaged out) and to the precision level of this measurement. In other theories, these black hole binaries are produced in very low metallicity galaxies (e.g. Cao et al. 2018; Mapelli et al. 2018), biased relative to the dark matter distribution differently than the luminous galaxies in our catalog. Annis et al. (in prep.) explore the effect of the tracer bias assumptions on the $H_{0}$ posterior for future analyses aiming at precision measurements.

Another assumption of our analysis that needs attention concerns the redshift likelihood. As anticipated above for the GW likelihood, this will not, in general, be well approximated by a Gaussian. In the future, we plan to explore the impact of realistic photometric redshift PDFs on the $H_{0}$ posterior, in order to enable precision cosmology with binary–black–hole events. An analysis with the full, asymmetric, GW likelihood will also be required. While an estimate of those effects is needed, tests on off–source lines–of–sight show that our constraint is likely not strongly impacted by the photo–$z$ training sample or systematic failures.

In the past two LVC observing seasons, black–hole mergers outnumbered neutron star events at a rate of approximately 10 to 1. Uncertainties on the expected detection rate are large, but conservative estimates predict $\sim 1$ event per week for the upcoming observing campaign (scheduled to start in April 2019). The majority of these events will have larger localization volumes than GW170814 (Chen & Holz 2016 estimate $\lesssim 1\%$ of BBHs will be localized to better than $10^{4}~{}$Mpc3) and hence provide poorer constraints than those reported here. However, given the high expected event rate for dark sirens, larger event samples will be available in the future. Chen et al. (2018) provide forecasts using a distribution of realistic localization area, finding that the dark siren method will reach $\sim 10\%$ statistical precision on $H_{0}$ by 2026 from BBH mergers only, and 5–10% precision from BNS mergers if none of them have EM counterparts.

5 Conclusions

In this paper, we have performed the first measurement of the Hubble constant using a gravitational wave detection of a binary–black–hole merger as a dark standard siren and the DES galaxies as a sample of potential host galaxies. Our analysis was blinded to avoid confirmation bias. Our main results, discussed in §4, include a measurement of $H_{0}=75^{+40}_{-32}~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$ for a flat prior within [20,140] km s*-1* Mpc*-1*, consistent with previous measurements of $H_{0}$. The $68\%$ confidence interval quoted here is $60\%$ of the uniform prior range, and it depends on the width of the prior. For example, with a broader prior of [10,220] ${\rm km~{}s^{-1}~{}Mpc^{-1}}$, we find $H_{0}=78^{+96}_{-24}~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$. Albeit weak, this measurement is not uninformative and the method becomes more powerful when we combine large numbers of dark sirens (Chen et al., 2018).

Future dark siren measurements will require complete galaxy catalogs. A wide field galaxy catalog with a DES–like depth is currently available only for $\sim 1/8$ of the sky. However, DES can be complemented with other datasets taken with DECam (such as the Blanco Imaging of the Southern Sky, BLISS, and the Dark Energy Camera Legacy Survey, DECals), to cover the whole Southern sky to a good depth ($r\sim 23.4$, 5$\sigma$ depth). An even deeper survey with more precise photo–z’s, such as the Large Synoptic Survey Telescope (LSST; LSST Dark Energy Science Collaboration 2012), would be of great value for further improving these constraints.

At the expected level of precision from hundreds of events ($<10\%$), systematics will play an important role. In future work, we plan to incorporate systematic uncertainties in our simulated data studies, in order to prepare for precision cosmology analyses on real data. We anticipate that some of the main sources of systematics will be photo–$z$ biases and catastrophic outliers, photo–$z$ training sample sample variance, galaxy catalog cuts and galaxy catalog completeness. In order to achieve the full potential of statistical standard siren cosmology, wide and deep galaxy surveys such as DES and LSST are necessary. Overall, our findings show that the synergy between gravitational wave black–hole merger detections and new generation large galaxy surveys will establish a new powerful probe for precision cosmology.

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(U. 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, UC 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 U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific and Industrial Research of India, the Department of Science and Technology, India, the Science & Engineering Research Board (SERB), India, the Ministry of Human Resource Development, India, the Spanish Agencia Estatal de Investigación, the Vicepresidència i Conselleria d’Innovació, Recerca i Turisme and the Conselleria d’Educació i Universitat del Govern de les Illes Balears, the Conselleria d’Educació, Investigació, Cultura i Esport de la Generalitat Valenciana, the National Science Centre of Poland, the Swiss National Science Foundation (SNSF), the Russian Foundation for Basic Research, the Russian Science Foundation, the European Commission, the European Regional Development Funds (ERDF), the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the Lyon Institute of Origins (LIO), the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, Industry Canada and the Province of Ontario through the Ministry of Economic Development and Innovation, the Natural Science and Engineering Research Council Canada, the Canadian Institute for Advanced Research, the Brazilian Ministry of Science, Technology, Innovations, and Communications, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Leverhulme Trust, the Research Corporation, the Ministry of Science and Technology (MOST), Taiwan and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, MPS, INFN, CNRS and the State of Niedersachsen/Germany for provision of computational resources.