Cosmological Inference using Gravitational Wave Standard Sirens: A Mock Data Challenge

TL;DR

This paper develops and tests a method for estimating the Hubble constant using gravitational wave standard sirens and galaxy catalogs, accounting for selection effects and catalog incompleteness, through mock data analyses.

Contribution

It introduces a comprehensive approach to estimate H_0 with both counterpart and galaxy catalog methods, including treatment of selection effects and catalog incompleteness.

Findings

Achieves 4.4% measurement precision with 250 detections and 50% complete catalogs.

Demonstrates unbiased H_0 recovery in simulated data.

Shows catalog incompleteness impacts measurement accuracy.

Abstract

The observation of binary neutron star merger GW170817, along with its optical counterpart, provided the first constraint on the Hubble constant $H_0$ using gravitational wave standard sirens. When no counterpart is identified, a galaxy catalog can be used to provide the necessary redshift information. However, the true host might not be contained in a catalog which is not complete out to the limit of gravitational-wave detectability. These electromagnetic and gravitational-wave selection effects must be accounted for. We describe and implement a method to estimate $H_0$ using both the counterpart and the galaxy catalog standard siren methods. We perform a series of mock data analyses using binary neutron star mergers to confirm our ability to recover an unbiased estimate of $H_0$. Our simulations used a simplified universe with no redshift uncertainties or galaxy clustering, but with different magnitude-limited catalogs and assumed host galaxy properties, to test our treatment of both selection effects. We explore how the incompleteness of catalogs affects the final measurement of $H_0$, as well as the effect of weighting each galaxy's likelihood of being a host by its luminosity. In our most realistic simulation, where the simulated catalog is about three times denser than the density of galaxies in the local universe, we find that a 4.4\% measurement precision can be reached using galaxy catalogs with 50\% completeness and $\sim 250$ binary neutron star detections with sensitivity similar to that of Advanced LIGO's second observing run.

Figures

Tables

Table 1. Table 1: A summary of the parameters present in the methodology.

Parameter	Definition
$H_0$	The Hubble constant.
$N_{\text{det}}$	The number of events detected during the observation period.
$x_{\text{GW}}$	The GW data associated with some GW source, $s$.
$D_{\text{GW}}$	Denotes that a GW signal was detected, i.e. that $x_{\text{GW}}$ passed some detection statistic threshold ${\rho }_{\text{th}}$.
$g$	Denotes that a galaxy is ($G$), or is not ($\overline{G}$), contained within the galaxy catalog.
$x_{\text{EM}}$	The EM data associated with some EM counterpart.
$D_{\text{EM}}$	Denotes that an EM counterpart was detected, i.e. that $x_{\text{EM}}$ passed some threshold.

Table 2. Table 2: A summary of the main results. We quote the peak value and the 68.3% highest density error region for the posterior probability on H 0 subscript 𝐻 0 H_{0} for each of the MDAs combining all 249 events. The fractional uncertainty is defined as the half-width of the 68.3% highest density probability interval divided by the simulated value of H 0 = 70 subscript 𝐻 0 70 H_{0}=70 km s − 1 ​ Mpc − 1 superscript km s 1 superscript Mpc 1 \mbox{km s}^{-1}\,\mbox{Mpc}^{-1} .

MDA	Host galaxy preference	Completeness777The completeness is calculated as a number completeness using Eq. (12) for MDAs 1 and 2, and as a luminosity completeness using Eq. (15) for MDA 3, out to a fiducial distance of 115 Mpc, such that it is indicative of the fraction of host galaxies which are inside the galaxy catalog in both cases.	$m_{\text{th}}$	Analysis assumption	$H_0$ (${\text{km s}}^{-1}{\text{Mpc}}^{-1}$)	Fractional uncertainty
0	Known host	-	-	direct counterpart	${69.08}_{-0.80}^{+0.79}$	$1.13$%
1	equal weights	100%	-	unweighted catalog	${68.91}_{-1.22}^{+1.36}$	$1.84$%
2a	equal weights	75%	19.5	unweighted catalog	${69.97}_{-1.50}^{+1.59}$	$2.21$%
2b	equal weights	50%	18	unweighted catalog	${70.14}_{-1.67}^{+1.80}$	$2.48$%
2c	equal weights	25%	16	unweighted catalog	${70.14}_{-2.18}^{+2.29}$	$3.20$%
3a	luminosity weighted	50%	14	weighted catalog	${70.83}_{-2.72}^{+3.55}$	$4.48$%
3b	luminosity weighted	50%	14	unweighted catalog	${69.50}_{-3.24}^{+4.20}$	$5.31$%

Table 3. Table 3: A summary of the parameters present in the network diagram, Fig. 9 .

Parameter	Definition
$H_0$	The Hubble constant
$x_{\text{GW}}$	The GW data associated with some GW source, $s$.
$D_{\text{GW}}$	Denotes that a GW signal was detected, i.e. that $x_{\text{GW}}$ passed some detection statistic threshold ${\rho }_{\text{th}}$.
$g$	Denotes that a galaxy is ($G$), or is not ($\overline{G}$), contained within the galaxy catalog.
$s$	Denotes that a GW signal was emitted.
$M$	Absolute magnitude.
$z$	Redshift.
$\Omega$	Sky location (right ascension and declination).
$d_L$	Luminosity distance.
$m$	Apparent magnitude.
$m_{\text{th}}$	Apparent magnitude threshold of the galaxy catalog.
${\rho }_{\text{th}}$	SNR threshold of the detector network.
Cosmological model	The cosmological model assumed for the analysis. Typically a Friedmann-Lemaître-Robertson-Walker universe.
Schechter parameters	The parameters which characterize the assumed absolute magnitude distribution of galaxies in the universe.
GW population	The assumed underlying population of GW sources.
Source frame parameters	Source frame parameters of a GW source, e.g. component masses, spins, inclination and polarization.
Detector frame parameters	As above, but redshifted into the detector frame.
Detector configuration	The network set up, including which detectors are included in the search and their noise floors.

Equations

$$H(z)=H_0\sqrt{{\Omega }_{\text{m}}(1+z{)}^3+{\Omega }_{\text{k}}(1+z{)}^2+{\Omega }_{\Lambda }},$$

$$d_L(z)$$

$$d_L(z)\approx cz/H_0,$$

$${\tilde{h}}_{+}(f)$$

$${\tilde{h}}_{\times }(f)$$

$$p(H_0|\{x_{\text{GW}}\},\{D_{\text{GW}}\})\propto p(H_0)p(N_{\text{det}}|H_0)\prod _i^{N_{\text{det}}}p({x_{\text{GW}}}_i|{D_{\text{GW}}}_i,H_0)$$

$$p(x_{\text{GW}}|D_{\text{GW}},H_0)$$

$$=\frac{p(x_{\text{GW}}|H_0)}{p(D_{\text{GW}}|H_0)},$$

$$p(D_{\text{GW}}|H_0)={\int }_{\rho >{\rho }_{\text{th}}}^{\infty }p(x_{\text{GW}}|H_0)d{x}_{\text{GW}}.$$

$$p(x_{\text{GW}}|D_{\text{GW}},H_0)$$

$$=p(x_{\text{GW}}|G,D_{\text{GW}},H_0)p(G|D_{\text{GW}},H_0)$$

$$+p(x_{\text{GW}}|\bar{G},D_{\text{GW}},H_0)p(\bar{G}|D_{\text{GW}},H_0).$$

$$p(x_{\text{GW}},$$

$$=\frac{p(x_{\text{GW}},x_{\text{EM}}|H_0)p(D_{\text{GW}},D_{\text{EM}}|x_{\text{GW}},x_{\text{EM}},H_0)}{p(D_{\text{GW}},D_{\text{EM}}|H_0)},$$

$$=\frac{p(x_{\text{GW}}|H_0)p(x_{\text{EM}}|H_0)}{p(D_{\text{EM}}|D_{\text{GW}},H_0)p(D_{\text{GW}}|H_0)}.$$

$$\varphi (L)dL=n^{\ast }{\left(\frac{L}{L^{\ast }}\right)}^{\alpha }{e}^{-L/L^{\ast }}\frac{dL}{L^{\ast }},$$

$$f_{\text{completeness}}(d_L)=\frac{{\sum }_j^{\text{MDA2}}\Theta (d_L-d_{L_j})}{{\sum }_k^{\text{MDA1}}\Theta (d_L-d_{L_k})},$$

$$p(L)=\beta \frac{L}{⟨L⟩}p(L)+(1-\beta )x(L)$$

$$x(L)=\dfrac{p(L)}{1-\beta}\Bigg{[}1-\beta\dfrac{L}{\langle L\rangle}\Bigg{]}\,,$$

$$f_{\text{completeness}}(d_L)=\frac{{\sum }_j^{\text{MDA3}}{L}_j\Theta (d_L-d_{L_j})}{{\sum }_k^{\text{complete}}{L}_k\Theta (d_L-d_{L_k})},$$

$$p(z,\Omega ,M,m|s,H_0)=p(m|z,\Omega ,M,s,H_0)p(z,\Omega ,M|s,H_0),$$

$$p(z,\Omega ,M,m|s,H_0)$$

$$=\delta (m-m(z,M,H_0))\frac{p(s|z)p(z)}{p(s)}p(\Omega )\frac{p(s|M,H_0)p(M|H_0)}{p(s|H_0)}.$$

$$p(s|M,H_0)$$

$$p(s|z)$$

$$p(x_{\text{GW}}|G,D_{\text{GW}},s,H_0)$$

$$p(z,\Omega ,M,m|G,s,H_0)$$

$$p(x_{\text{GW}}|G,D_{\text{GW}},s,H_0)$$

$$p($$

$$p(G|D_{\text{GW}},s,H_0)$$

$$=\text{iiiint}\Theta [m_{\text{th}}-m]\frac{p(D_{\text{GW}}|z,\Omega ,M,m,s,H_0)p(z,\Omega ,M,m|s,H_0)}{p(D_{\text{GW}}|s,H_0)}dzd\Omega dMdm,$$

$$=\frac{1}{p(D_{\text{GW}}|s,H_0)}\text{iiiint}\Theta [m_{\text{th}}-m]p(D_{\text{GW}}|z,\Omega ,s,H_0)p(z,\Omega ,M,m|s,H_0)dzd\Omega dMdm.$$

$$p(G|D_{\text{GW}},s,H_0)$$

Full text

I Introduction

The idea that gravitational waves (GW) detections can be used for the inference of cosmological parameters, such as the Hubble constant ($H_{0}$), was first proposed over three decades ago by Bernard Schutz Schutz (1986). The key to this process is that GW signals from compact binary coalescences act as standard sirens, in the sense that they provide a self-calibrated luminosity distance to the source. This can be obtained directly from the GW signal, and is therefore entirely independent of the cosmic distance ladder Holz and Hughes (2005); Dalal et al. (2006); Nissanke et al. (2010, 2013); Chen et al. (2018); Feeney et al. (2019); Di Valentino et al. (2018); Mortlock et al. (2019); Farr et al. (2019). With the addition of redshift information for each source we then have the required input for cosmological inference.

At the time of writing, the current percent level state-of-the-art electromagnetic (EM) measurements of $H_{0}$ are in tension with each other. The Planck experiment uses measurements of cosmic microwave background (CMB) anisotropies and provides a value of $H_{0}=67.4\pm 0.5$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ Aghanim et al. (2018). The Supernovae, $H_{0}$, for the Equation of State of Dark energy (SH0ES) experiment measures distances to Type Ia supernovae standard candles making use of the cosmic distance ladder, and gives $H_{0}=74.03\pm 1.42$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ Riess et al. (2019). These two independent measurements of $H_{0}$ are in tension at the level of $\sim 4.4$-$\sigma$ Riess et al. (2019). While the early-universe Planck measurements are also favored by measurements using supernovae calibrated with Baryon Acoustic Oscillations Macaulay et al. (2019), and the SH0ES results agree with local gravitational lensing measurements by the H0LiCOW Collaboration Birrer et al. (2019), calibration of supernovae using the Tip of the Red Giant Branch yields $H_{0}$ midway between the two Freedman et al. (2019).

This indicates the possibility that at least one of these measurements is subject to unknown systematics, or it could be an indication of new physics causing the discrepancy between the local measurements and the non-local (early universe) CMB based measurement. This makes a GW standard siren measurement of $H_{0}$ particularly interesting, as this will provide an alternative local constraint on $H_{0}$. In this manner, the use of GWs as standard sirens may allow us to arbitrate the current situation, indicating either a bias in the current measurements, or pointing towards new physics.

The detection of the binary neutron star (BNS) event GW170817 Abbott et al. (2017a), together with its optical counterpart (Abbott et al., 2017b; Soares-Santos et al., 2017) led to the first standard siren measurement of $H_{0}$ Abbott et al. (2017c). The counterpart associated with GW170817 allowed for the identification of its host galaxy, NGC4993, and hence a direct measurement of its redshift, which in turn resulted in the inferred value $H_{0}$=$70^{+12}_{-8}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$. Future counterpart standard siren measurements are expected to constrain $H_{0}$ to the percent level Dalal et al. (2006); Nissanke et al. (2010, 2013); Chen et al. (2018); Feeney et al. (2019).

Central to the aims of this paper is the case where an EM counterpart is not observed, and how $H_{0}$ inference can still be performed. In particular, the method proposed by Schutz in 1986 Schutz (1986); Del Pozzo (2012) allows the use of galaxy catalogs to provide redshift information for potential host galaxies within the event’s GW sky-localization. The idea is that, by marginalizing over the possible discrete values of redshift for each GW detection we account for uncertainty as to which galaxy is the true host. By combining the information from many GW events, the contributions from the true host galaxies will grow since they will all share the same true $H_{0}$. Contributions from the others will statistically average out, leading to a constraint on $H_{0}$ and possibly other cosmological parameters.

Over the course of the first observing run (O1) and the second observing run (O2) a total of 11 GW events were detected by the advanced LIGO and Virgo detectors: 10 are binary black hole (BBH) events and one is the above-mentioned BNS event GW170817 Abbott et al. (2019a). The “galaxy catalog” method has been independently applied to both the BNS event GW170817 (without assuming NGC4993 is the host) Fishbach et al. (2019), and the BBH event GW170814 Soares-Santos et al. (2019) resulting in posterior probability distributions on $H_{0}$ where the posterior from GW170814 was broader than (but consistent with) that obtained from GW170817. The difference in the widths of the $H_{0}$ constraints is an expected result due to the larger localization volume associated with GW170814, and the high number of galaxies it contained. Using the detections from O1 and O2, multiple GW events have been combined to give the latest standard siren measurement of $H_{0}$ Abbott et al. (2019b) using the methodology presented in this paper.

Predictions suggest that it will be possible to constrain $H_{0}$ to less than $2\%$ within 5 years of the start of the third observing run (O3) and to $1\%$ within a decade, though this is dependent on the number of events observed with EM counterparts Chen et al. (2018), and this may change as our understanding of astrophysical rates improves, and would require the detector amplitude calibration error to be measured to better than this precision. Simulations in Chen et al. (2018) and Fishbach et al. (2019), which assume complete catalogs based on realistic large-scale structure simulations, find that for BNSs without counterparts, the convergence is $40\%/\sqrt{N}$. The convergence found there for BBHs is much slower, as BBHs are typically detected at greater distances with larger localization volumes.

The prospects of identifying a transient EM counterpart will certainly increase, and correspondingly, the number of candidate host galaxies in a catalog will decrease, with improved event sky-localizations as future GW observatories join the detector network Abbott et al. (2018). With the Japanese detector KAGRA Abbott et al. (2018) having joined O3 in early 2020, and LIGO-India approved for construction Padma (2019), the next decade of standard siren cosmology is set to be very exciting.

O3 began on April 1st 2019 and consists of 11 months’ worth of data. The sensitivities of the LIGO and Virgo detectors have improved since O2, leading to an increased detection rate of GW candidates111In the first half of O3 the detectors averaged the detection of one GW candidate per week. If all of these candidates are ultimately identified as real GW events, then O3 within its first two months will have exceeded the total number of detections of O1 and O2. Abbott et al. (2019c). This is the first observing run for which there will be 3 detectors operating for the entirety of the run. Having more detectors improves the duty-cycle of the network, i.e. the fraction of run time for which one or more detectors in the network is online, and also increases the rate of three-detector detections, which will likely be better localized on the sky than the two-detector ones. This is important, both in terms of performing EM follow-up for EM counterparts practically Abbott et al. (2019), and for reducing the number of possible host galaxies for events in the case where a counterpart is not observed.

This paper presents the Bayesian framework behind the gwcosmo code, a product of the LIGO and Virgo Collaborations (LVC) which was used to measure $H_{0}$ using detected GW events from O1 and O2 Abbott et al. (2019b). The method detailed in this paper is also expected to be implemented in future LIGO/Virgo/KAGRA standard siren measurements. We present results from a series of mock data analyses which were designed specifically to test this method’s robustness against some of the most common pitfalls, in particular, GW selection effects which affect all $H_{0}$ measurements, and EM selection effects, which are relevant in the context of galaxy catalogs. This method builds upon the Bayesian framework first presented in Del Pozzo (2012) which has subsequently been extended, modified and independently derived by multiple authors Nissanke et al. (2013); Chen et al. (2018); Nair et al. (2018); Fishbach et al. (2019); Soares-Santos et al. (2019). The framework here is broadly equivalent to that in Chen et al. (2018); Fishbach et al. (2019), however the mathematics and implementation differ, most notably in the treatment of EM selection effects. With specific care regarding selection effects we outline methods for constraining $H_{0}$ using both the “galaxy catalog” and “EM counterpart” approaches.

This paper is the first to explicitly test the robustness of a coded implementation of this methodology through use of galaxy catalogs which are incomplete and do not contain all of the GW host galaxies. Additionally, the GW data used in these MDAs were produced using an end-to-end simulation, including searching for “injected” signals in real detector data followed by a full parameter estimation to obtain the GW posterior samples Singer et al. (2014); Berry et al. (2015), making this the most realistic set of simulated GW data to be used to explore GW cosmology to date. The analyses start with the most simplistic scenario, and increase in complexity with each iteration in order to ensure that the gwcosmo code is able to pass each level satisfactorily before moving onto the next.

This paper is structured as follows. Section II presents the Bayesian framework used to estimate the posterior on $H_{0}$. Section III discusses the design and preparation of the MDAs. In Section IV we present our results. We conclude in Section V giving a detailed discussion of results and providing guidance for future work. Some of the details of the Bayesian method have been set aside to be discussed in an Appendix.

II Methodology

The late-time cosmological expansion in a Friedmann-Lemaître-Robertson-Walker universe is characterized by the Hubble-Lemaître parameter as a function of the redshift $z$,

[TABLE]

where $H_{0}$ is the Hubble constant, the rate of expansion in the current epoch, and $\Omega_{\text{m}}$ and $\Omega_{\Lambda}$ are the fractional matter density (including baryonic and cold dark matter) and fractional dark energy density (assumed to be due to a cosmological constant) respectively; $\Omega_{\text{k}}$ is the fractional curvature energy density which is identically zero for a “flat” universe consistent with observations. Additionally, we have the constraint $\Omega_{\text{m}}+\Omega_{\text{k}}+\Omega_{\Lambda}=1$ for all the components contributing to the energy density of universe at the present epoch.

The expansion history of the universe maps to a “redshift-distance relation” associating the redshift $z$ of observable sources to their luminosity distance $d_{L}(z)$ (see e.g. Hogg (1999)) as,

[TABLE]

for a flat universe. From the relation between observed $z$ and $d_{L}$ to sources (EM sources such as variable stars or supernovae, or GW sources), one can measure the cosmological parameters appearing in $H(z)$. With knowledge of the other cosmological parameters $\{\Omega_{\text{m}},\Omega_{\text{k}},\Omega_{\Lambda}\}$ coming from independent observations, the redshift-distance relation can be used to measure $H_{0}$. We would like to note that with prior knowledge on the other cosmological parameters coming from EM observations, the measurement made with GW detections are not strictly independent measurements.

At low redshifts $z\ll 1$, the redshift-distance relation can be approximately described by the linear Hubble relation,

[TABLE]

which contains $H_{0}$ but is independent of the other cosmological parameters. With this approximate linear relation at low redshifts, any measurement of $H_{0}$ with GWs is independent of the values of the other cosmological parameters.

II.1 Standard Sirens

The amplitude of the observed strain is inversely proportional to the luminosity distance to the GW source. For compact binary sources in quasi-circular orbits, the two polarizations of the gravitational wave signal can be written to leading order as a function of frequency $f$ as  Maggiore (2008)

[TABLE]

where $\phi(\mathcal{M}_{z},t)$ is the phase of the signal. The redshifting of the signal is accounted for by using the parameter $\mathcal{M}_{z}\equiv\mathcal{M}(1+z)$, the “redshifted chirp mass,” to describe the signal as observed in the detector. Since $\mathcal{M}_{z}$ appears in both the phase and the amplitude, and in practice is more strongly constrained by $\phi(\mathcal{M}_{z},f)$, the dominant uncertainty on the signal amplitude results from the uncertainties on luminosity distance, $d_{L}$ and inclination angle $\iota$. Each detector sees a linear combination of the two polarizations, $\tilde{h}(t)=F_{+}\tilde{h}_{+}+F_{\times}\tilde{h}_{\times}$, where $F_{+,\times}$ are the antenna response functions of the detector, which vary over the sky position and polarization angle of the source. Given multiple detectors at distant sites it is possible to simultaneously infer the parameters of the source, and therefore find a direct estimate of its luminosity distance Veitch et al. (2015). This makes compact binaries self-calibrated luminosity distance indicators or “standard sirens” unlike EM distance indicators which need to undergo calibration via multiple rungs of the cosmic distance ladder. The redshift of the GW source, also required for cosmological inference, remains degenerate with the source’s mass, contained within $\mathcal{M}_{z}$, and needs to be estimated in alternate ways. The precision of the $d_{L}$ estimate is limited because of correlations with other parameters, particularly the inclination angle $\iota$ Nissanke et al. (2013). In this work we simulate these effects as part of our end-to-end analysis, described in Sec. III.

II.2 Galaxy Information

There are multiple ways in which EM observations can provide complementary redshift222There are ways of obtaining the redshift independent of EM observations, by using known population properties such as the mass distribution Taylor et al. (2012); Farr et al. (2019), or the neutron star equation-of-state Messenger and Read (2012). information. A BNS event may be detected in coincidence with an EM counterpart, which can be associated with the host galaxy to provide a direct measurement of the redshift of the source. More generically, a GW event may not have a detected EM counterpart, in which case one needs to fall back on the method outlined by Schutz Schutz (1986) and use potential host galaxies within the event’s sky localization region for the redshift information for the source. Two possibilities come up: (i) to use available galaxy catalogs, or (ii) to conduct dedicated EM follow-up on the event’s sky region, mapping the galaxies within that area to as great a depth as possible to maximize the redshift information available.

When using galaxy catalogs to provide the prior redshift information, the possibility that the host galaxy lies beyond the reach of the catalog must be taken into account. EM telescopes are flux limited, which means that galaxy catalogs are inherently biased towards containing objects which are brighter and/or nearer-by (although there may be other selection effects due to galaxy color or size, depending on the catalog). These EM selection effects must be accounted for. Carrying out dedicated EM follow-up will, to some degree, mitigate this issue, as it will allow for far deeper coverage over a small section of the sky. For nearby events, the possibility that the host galaxy lies above the telescope’s upper threshold may be negligible. However, the time and resources required for dedicated EM follow-up means that the default approach for GW events observed without counterparts will be to use pre-existing catalogs.

In either case, the uncertainty associated with each galaxy’s redshift must be taken into account, including the redshift error due to the galaxy’s peculiar velocity, $v_{p}$, and, in cases where the redshift is estimated photometrically, a much larger uncertainty due to the photometric algorithm. Peculiar velocities are significant for nearby galaxies. The effect of the peculiar velocity on the measurement of $H_{0}$ may be small if there are a large number of potential host galaxies in the GW event’s sky-localization, but for a small number of galaxies, and for the counterpart case, this effect is particularly noticeable. For GW170817 at a nearby distance of about $40\,\text{Mpc}$, the peculiar velocity contribution was large as $10\%$ of the total observed redshift Abbott et al. (2017c), and different procedures of reconstructing the peculiar velocity field led to residual uncertainties on the redshift of between $2\%$ and $8\%$ Abbott et al. (2017c); Springob et al. (2014); Carrick et al. (2015); Hjorth et al. (2017); Guidorzi et al. (2017). The impact on $H_{0}$ measurement of peculiar velocities and their reconstruction is of topical interest, and has been the subject of several recent studies including Howlett and Davis (2020); Mukherjee et al. (2019); Nicolaou et al. (2019). Photometric redshifts on the other hand are important slightly farther away due to lack of spectroscopic data in galaxy catalogs. The “photo-$z$” are estimated using fitting and machine learning algorithms De Vicente et al. (2016); Sadeh et al. (2016), which often have large $\mathcal{O}(1)$ fractional uncertainties associated with them. While various caveats and subtleties for a realistic measurement have been outlined in Abbott et al. (2019b), the impact of photo-$z$ uncertainties on $H_{0}$ measurement is not precisely quantified in literature yet. Our present mock data analyses ignore these crucial redshift uncertainties altogether, and the impact of their magnitudes, profiles, and other systematic artefacts are left aside for possible future study.

II.3 Bayesian Framework

This section presents an overview of the Bayesian framework of the gwcosmo methodology. Parameters which appear explicitly in this overview are defined in Table 1, while Table 3 in Appendix A.2 provides an extended list of parameter definitions, alongside a network diagram which demonstrates the conditional dependence of these parameters (see Fig. 9).

The posterior probability on $H_{0}$ from $N_{\text{det}}$ GW events is computed as follows:

[TABLE]

where $\{x_{\text{GW}}\}$ is the set of GW data, $D_{\text{GW}}$ indicates that the event was detected as a GW and $p(H_{0})$ is the prior on $H_{0}$. For a given $H_{0}$, the term $p(N_{\text{det}}|H_{0})$ is the probability of detecting $N_{\text{det}}$ events. It depends on the intrinsic astrophysical rate of events in the source frame, $R=\frac{\partial{N}}{\partial V\partial T}$. The total number of expected events is given by $N_{\text{det}}=R\,\langle VT\rangle$, where $\langle VT\rangle$ is the average of the surveyed comoving volume multiplied by the observation time. By choosing a scale-free prior on rate, $p(R)\propto 1/R$, the dependence on $H_{0}$ drops out Fishbach et al. (2018). For simplicity this approximation is made throughout the analysis and therefore $p(N_{\text{det}}|H_{0})$ is absent from further expressions.

The remaining term factorizes into likelihoods for each detected event. Using Bayes’ theorem we can write it as,

[TABLE]

where we set $p(D_{\text{GW}}|x_{\text{GW}},H_{0})=1$, since the analysis is only carried out when the signal-to-noise ratio (SNR), $\rho$, associated with $x_{\text{GW}}$ passes some detection statistic threshold $\rho_{\text{th}}$ – it is a prerequisite that the event has been detected. Calculating $p(D_{\text{GW}}|H_{0})$ requires integrating over all possible realizations of GW events, with a lower integration limit of $\rho_{\text{th}}$:

[TABLE]

For explicit details on the calculation of $p(D_{\text{GW}}|H_{0})$ see Appendix A.5. The term $p(D_{\text{GW}}|H_{0})$ depends on properties of the GW source population (e.g. the mass distribution), but in this work, for simplicity, it is assumed that the population properties are known exactly.

II.3.1 The galaxy catalog method

In the galaxy catalog case, the EM information enters the analysis as a prior, made up of a series of possibly smoothened delta functions333While uncertainties on the galaxy sky-coordinates can be safely ignored, the error on the redshift can be modeled with a Gaussian or a more complicated distribution. at the redshift, right ascension (RA) and declination (dec) of the possible source locations. As we are in the regime where (especially for BBHs) galaxy catalogs cannot be considered complete out to the distances to which GW events are detectable, we have to consider the possibility that the host galaxy is not contained within the galaxy catalog due to being dimmer than the apparent magnitude threshold. In order to do so, we marginalize the likelihood over the case where the host galaxy is, and is not, in the catalog (denoted by $G$ and $\bar{G}$ respectively):

[TABLE]

While theoretically equivalent to and consistent with the methodology presented in Chen et al. (2018); Fishbach et al. (2019), the mathematics and implementation here differ, most notably in the treatment of EM selection effects, and our focus on whether the host galaxy is contained within the galaxy catalog or not, rather than calculating a “completeness fraction” in order to weight the in-catalog and out-of-catalog likelihood contributions. This, alongside the modeling of EM selection effects using an apparent magnitude threshold, which has not been done before, accounts for the main differences between this derivation and those presented in earlier works. The methodology presented here aligns directly with the implementation of the gwcosmo code. We leave the details of this derivation to Appendix A.2.

II.3.2 The counterpart method

The method outlined above is for the galaxy catalog case, in which no EM counterpart is observed, or expected. We also consider the case where we observe an EM counterpart. The main difference is the inclusion of a likelihood term for the EM counterpart data, mirroring that of the GW data.

The likelihood in this case, which is the term within the product in Eq. (6), is given by:

[TABLE]

where $x_{\text{EM}}$ refers to the EM counterpart data and $D_{\text{EM}}$ denotes that the counterpart was detected. In the numerator we have assumed that the GW and EM data are independent of each other and so the joint GW-EM likelihood factors out. $p(D_{\text{GW}},D_{\text{EM}}|x_{\text{GW}},x_{\text{EM}},H_{0})$ is further factorized as $p(D_{\text{EM}}|D_{\text{GW}},x_{\text{GW}},x_{\text{EM}},H_{0})p(D_{\text{GW}}|x_{\text{GW}},x_{\text{EM}},H_{0})$. The first term is equal to 1, as this method is only used when we have observed an EM counterpart, meaning that by definition $x_{\text{EM}}$ has passed some threshold for detectability set by EM telescopes. The second term also goes to 1, due to the same threshold argument as in section II.3.

For simplicity, in this paper we make the assumption that the detection of an EM counterpart is flux-limited and, as in Abbott et al. (2017c), that the detectability of EM counterparts extends well beyond the distance to which BNSs are detectable with O2-like LIGO and Virgo sensitivity. Following this, we make the assumption that the term $p(D_{\text{EM}}|D_{\text{GW}},H_{0})\approx 1$, and leave a more rigorous analysis of the $H_{0}$-dependence of this term for a future study.

In an ideal scenario, the observation of an EM counterpart will allow for the identification of one of the galaxies in the neighboring region as the host of the GW event. In the case where the EM counterpart cannot be unambiguously linked to a host galaxy, this uncertainty can also be taken into account. See Appendix A.4 for more details.

III The Mock Data Analyses

In this section we describe a series of mock data analyses (MDAs) that we use to test our implementation of the Bayesian formalism described in Section II and its ability to infer the posterior on $H_{0}$ under different conditions. For each case, the MDA consists of (i) simulated GW data, and (ii) a corresponding mock galaxy catalog. In all cases, we make several idealized assumptions regarding both the GW and galaxy data. On the GW side, the detection efficiency and the source population properties are assumed to be known exactly. On the galaxy side, the luminosity function and magnitude limit are also assumed to be known exactly in each case, so that the incompleteness correction can be calculated exactly. Further, we neglect the effects of large-scale structure and redshift uncertainties in the mock catalogs.

For each of the MDAs we use an identical set of simulated BNS events from The First Two Years of Electromagnetic Follow-Up with Advanced LIGO and Virgo dataset Singer et al. (2014); Berry et al. (2015)444The set of simulations in Berry et al. (2015) are more realistic with the same injections in (recolored) detector data as opposed to Gaussian noise used in Singer et al. (2014). Correspondingly, the detection criterion is in terms of a false alarm rate (FAR) rather than a threshold on the SNR. This is an important distinction, particularly affecting events marginally close to the detection threshold. We use the simplified set of simulations in Singer et al. (2014) noting potential caveats.. The set of BNS events comes from an end-to-end simulation of approximately 50,000 “injected” events in detector noise corresponding to a sensitivity similar to what was achieved during O2. Only a subset (approximately 500 events) were “detected” by a network of two or three detectors with the GstLAL matched filter based detection pipeline Messick et al. (2017). From the above detections, 249 events were randomly selected (in a way that no selection bias was introduced), and these events underwent full Bayesian parameter estimation using the LALInference software library Veitch et al. (2015) to obtain gravitational wave posterior samples and skymaps. Consistency with the First Two Years parameter estimation results in terms of sky localization areas and 3D volumes was demonstrated in Del Pozzo et al. (2018). It is these 249 events of the First Two Years dataset and the associated GW data which we use for our analysis.

The galaxy catalogs for each iteration of the MDA described below are designed to test a new part of the gwcosmo methodology in a cumulative fashion, starting with GW selection effects, adding in EM selection effects, and finally testing the ability to utilize the information available in the observed brightness of host galaxies, by weighting the galaxies with a function of their intrinsic luminosities.

The starting point for the galaxy catalogs is to take all 50,000 injected events from the First Two Years dataset and simulate a mock universe, which contain a galaxy corresponding to each injected event’s sky location and luminosity distance, where the latter is converted to a redshift using a fiducial “simulated” $H_{0}$ value of $70$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$. The First Two Years data was originally simulated in a universe where GW events followed a $d_{L}^{2}$ distribution, and there was no distinction between the source frame and the (redshifted) detector frame masses. Though not ideal, this data reasonably mimics a low redshift universe ($z\ll 1$) in which the linear Hubble relation of Eq. (3) holds, and galaxies follow a $z^{2}$ distribution. We use the same linear relation for the generation of the MDA universe (i.e. a set of simulated galaxy catalog parameters) for each of the MDAs. It should be emphasized that the Bayesian method for estimation of $H_{0}$ outlined in Section II above is general, and can be extended to realistic scenarios with a non-linear cosmology with $\{\Omega_{\text{m}},\Omega_{\text{k}},\Omega_{\Lambda}\}$ held fixed. So, in particular, the method is applicable for events which are detected at higher distances, where the low redshift approximation breaks down. The restriction to a linear cosmology in this paper comes only due to the use of the MDA dataset. We would like to note that by using a linear cosmology, we are not testing possible effects introduced by the presence of other cosmological parameters. The analysis at large redshifts may, for example, be sensitive to the values (or the assumed prior ranges) of the parameters like $\Omega_{\text{m}}$ and $\Omega_{\Lambda}$.

The first four columns of Table 2 summarize the characteristics of each of the galaxy catalogs created and how they correspond to each MDA. We give a brief description for each of the cases below.

III.1 MDA0: Known Associated Host Galaxies

MDA0 is the simplest version of the MDAs, in which we identify with certainty the host galaxy for each GW event, and is equivalent to the direct counterpart case. As the galaxies are generated with no redshift uncertainties or peculiar velocities, the results will be (very) optimistic. This MDA provides the “best possible” constraint on $H_{0}$ using the 249 events, which then allows for comparison with the other MDAs.

III.2 MDA1: Complete Galaxy Catalog

The MDA1 universe consists of the full set of 50,000 galaxies out to $z\approx 0.1$ ($d_{L}\approx 428\,\text{Mpc}$) in the original First Two Years dataset. This gives a galaxy number density of $\sim$ 1 per $7000\,\text{Mpc}^{3}$, which is $\sim 35$ times sparse compared to the actual density of galaxies in the local universe Dálya et al. (2018). Additional galaxies are generated beyond the edge of the dataset universe, uniformly across the sky and uniformly in comoving volume, thereby extending the universe out to a radius of 2000 Mpc ($z=0.467$ for $H_{0}=70$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$). This means that, even allowing $H_{0}$ to be as large as 200 $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$, the edge of the MDA universe is more than twice the highest redshift associated with the farthest detection (which is at $\sim 270$ Mpc)555For MDA1 and for all subsequent MDAs, it has been tested that the artificial “edge of the universe” has no bearing on the results.. Each of the 249 detected BNS have a unique associated host galaxy contained within the MDA1 catalog. This catalog is thus complete in the sense that it contains every galaxy in the simulated universe. We refer to the MDA universe as MDA1 throughout the rest of the paper, and similarly for the subsequent MDAs.

MDA1 is designed to test our treatment of GW selection effects, by ensuring that given a set of sources and access to a complete catalog, our methodology and analysis produces a result consistent with the simulated value of $H_{0}$.

III.3 MDA2: Incomplete Galaxy Catalog

MDA2 is designed to test our treatment of EM selection effects, by applying an apparent magnitude threshold to the MDA universe, such that a certain fraction of the host galaxies is not contained in it. This is a necessary consideration, given that we are in the regime where GW signals are being detected beyond the distance to which the current galaxy catalogs can be considered to be complete. This has been true for BBHs detections since O1, and is true of BNSs as well in O3.

In order to create the catalog for MDA2, we start with the initial MDA1 universe and assign luminosities to each of the galaxies within it. We assume that the luminosity distribution of the galaxy catalog is known to the observer throughout and follows a Schechter function of the form Schechter (1976)

[TABLE]

where $L$ denotes a given galaxy luminosity and $\phi(L)\,dL$ is the number of galaxies within the luminosity interval $[L,L+dL]$. The characteristic galaxy luminosity is given by $L^{*}=1.2\times 10^{10}\,h^{-2}\,L_{\odot}$ with solar luminosity $L_{\odot}=3.828\times 10^{26}\,\mathrm{W}$, and $h\equiv H_{0}/(100\mbox{km s}^{-1}\,\mbox{Mpc}^{-1})$666We note that the parameter $L^{*}$ of the Schechter luminosity function itself depends on $H_{0}$, which we allow to vary and hence take into account within our formalism., $\alpha=-1.07$ characterizes the exponential drop off of the luminosity function, and $n^{*}$ denotes the number density of objects in the MDA universe (in practice, this only acts as a normalization constant). The integral of the Schechter function diverges at $L\rightarrow 0$, requiring a lower luminosity cutoff for the dimmest galaxies in the universe which we set to $L_{\text{lower}}=0.001L^{*}$. This choice is arbitrary for our purpose here, but small enough to include almost all objects classified as galaxies in real catalogs like GLADE Dálya et al. (2018).

These luminosities are then converted to apparent magnitudes using $m\equiv 25-2.5\log_{10}(L/L^{*})+5\log_{10}(d_{L}/\text{Mpc})$, and an apparent magnitude threshold $m_{\text{th}}$ is applied as a crude characterization of the selection function of an optical telescope observing only objects with $m<m_{\text{th}}$. MDA2 is broken into three sub-MDAs, in order to test our ability to handle different levels of galaxy catalog completeness dictated by different telescope sensitivity thresholds. In each case, the catalog completeness is defined as the ratio of the number of galaxies inside the catalog relative to the number of galaxies inside the MDA universe, out to a reference fiducial distance $d_{L}$,

[TABLE]

where the numerator is a sum over the galaxies contained within the MDA2 catalog out to some reference distance $d_{L}$, and the denominator is a sum over the galaxies in the MDA1 catalog.

Apparent magnitude thresholds of $m_{\text{th}}=19.5$, $18$, and $16$ are chosen for the three sub-MDAs, which correspond to cumulative number completeness fractions of $~{}75\%$, $~{}50\%$ and $~{}25\%$ respectively, evaluated at a distance of $d_{L}=115\,\text{Mpc}$, chosen such that given the luminosity distance distribution of detected BNSs, the completeness fraction for the sub-MDA to this distance is roughly indicative of the percentage of host galaxies which remain inside the galaxy catalog. The left panel of Fig. 1 shows how the completeness of each of the MDA2 catalogs drop off as a function of distance.

III.4 MDA3: Luminosity Weighting

MDA3 is designed to test the effect of weighting the likelihood of any galaxy being host to a GW event as a function of their luminosity. It is probable that the more luminous galaxies are also more likely hosts for compact binary mergers; the luminosity in blue (B-band) is indicative of a galaxy’s star formation rate, for example, while the luminosity in high infrared (K-band) is a tracer of the stellar mass Singer et al. (2016); Leibler and Berger (2010); Fong et al. (2013). The bulk of the host probability is expected to be contained within a smaller number of brighter galaxies, effectively reducing the number of galaxies which need to be considered. Additional information from luminosity is thus expected to improve the constraint on $H_{0}$ by narrowing its posterior probability density.

For MDA3, the probability of a galaxy hosting a GW event is chosen to be proportional to the galaxy’s luminosity. Because the GW events for these MDAs were generated in advance, and we are retroactively simulating the universe in which they exist, generating the MDA3 universe required some care: luminosities have to be assigned to the host galaxies and the non-host galaxies in such a way that our choice of simulated luminosity weighting is correctly represented within the galaxy catalog.

As with MDA2, the luminosity distribution of the galaxies in the universe is assumed to follow the Schechter luminosity function as in Eq. (11) (referred to from now on as $p(L)$). However, the joint probability of a single galaxy having luminosity $L$ and hosting a GW event (which emits a signal, $s$) is $p(L,s)\propto L\,p(L)\,,$ where we assume that the probability of a galaxy of luminosity $L$ hosting a source is proportional to the luminosity itself . All host galaxies thus have luminosities sampled from $L\,p(L)$. In this context, we must consider all galaxies which hosted GW events, not just those from which a signal was detected. With this in mind, the overall luminosity distribution has the following form:

[TABLE]

where $\beta$ is the fraction of host galaxies to total galaxies for the observed time period ($1\geq\beta\geq 0$), $L/\langle L\rangle$ is the normalized luminosity, and $x(L)$ is the unknown luminosity distribution of galaxies which did not host GW events, which we can sample for a given value of $\beta$.

Rearranging to obtain the only unknown, $x(L)$, gives

[TABLE]

and from this we see there is an additional constraint on $\beta$, because the term inside the brackets must be $>0$. The maximum value that $\beta$ can take is given by $\beta_{\text{max}}=\langle L\rangle/L_{\text{max}}$, where $L_{\text{max}}$ is the maximum luminosity from the Schechter function, and $\langle L\rangle$ is the mean. From the Schechter function parameters detailed in section III.3, $\beta_{\text{max}}\approx 0.015$.

The original First Two Years data was generated by simulating $\sim 50,000$ BNS events, of which $\sim 500$ were detected, of which 249 randomly selected detections underwent parameter estimation. The number of “hosting” and “non-hosting” galaxies have to be rescaled to represent this. Thus half of the original galaxies were denoted as hosts (including those associated with the 249 detected GW events). However, in order to satisfy the requirements for $\beta$, a greater density of non-hosting galaxies had to be added to the universe before luminosities could be assigned. Thus for MDA3, the density of galaxies is increased by a factor of 100, with the acknowledgement that this would lead to a broadening of the final posterior. MDA3 has a galaxy density of $\sim 1$ galaxy per $70\,\text{Mpc}^{3}$, which is about $3$ times denser than the actual density of galaxies in the local universe Dálya et al. (2018). This also means that MDA3 is not directly comparable with the previous MDA versions, save MDA0. The galaxies which are hosts are assigned luminosities from $Lp(L)$, and non-hosts from $x(L)$ above.

In order to include EM selection effects, an apparent magnitude cut $m_{\text{th}}$ of 14 is applied, such that the completeness of the galaxy catalog is $\sim 50\%$ out to the same fiducial distance of $115\,\text{Mpc}$ as in MDA2. In this case, completeness is however defined in terms of the fractional luminosity contained in the catalog, rather than in terms of numbers of objects:

[TABLE]

where the numerator is summed over the galaxies inside the MDA3 apparent magnitude-limited catalog, and the denominator is summed over the galaxies in the whole MDA3 universe. This is shown in the right panel of Fig. 1. As the host galaxies are luminosity weighted, the cumulative luminosity completeness is representative of the percentage of BNS event hosts inside the catalog.

IV Results

In this section we summarize the results for the mock data analyses described in Section III. We show the combined posteriors on $H_{0}$ for each MDA, discuss the convergence to the simulated value of $H_{0}$ = $70$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$and calculate the precision of the combined measurement under each set of conditions. In Table 2 we list the measured values of the Hubble constant for the combined 249 event posterior (maximum a-posteriori and 68.3% highest density posterior intervals) all computed with a uniform prior on $H_{0}$ in the range of $[20,200]\ \mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$, as well as the corresponding fractional uncertainties for each of the MDAs.

IV.1 MDA0: Known Associated Host Galaxies

We first consider the simple case where we identify the true host galaxy for every event and determine the resulting 249-event combined $H_{0}$ posterior. Fig. 2 presents the results of this analysis. The likelihoods for each individual GW event are shown (normalized relative to each other but scaled with respect to the combined posterior for clarity) shaded by the event’s optimal SNR in the detector network, as defined in Cutler and Flanagan (1994). In this case, each likelihood is informative, having a clearly-defined peak corresponding to finding the likely values of $H_{0}$ for the known galaxy redshift. Each curve traces the information in the corresponding $d_{L}$ distribution, which is usually unimodal, but in some cases may have two or more peaks Singer et al. (2014); Berry et al. (2015). We see that the peaks of the individual likelihoods do not necessarily correspond to the true value $H_{0}=70\,\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$, but there is always support for it, leading to the combined posterior, which is overlaid in thick purple. This gives us a statistical estimate for the maximum a-posteriori value and $68.3\%$ maximum-density credible interval for $H_{0}$ as $69.08^{+0.79}_{-0.80}\,\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$. The final result combining all the 249 events have converged to a relatively symmetric “Gaussian” distribution Walker (1969).

The result of MDA0 provides us with the best possible $H_{0}$ estimate given the set of GW detections, since this case corresponds to perfect knowledge of the host galaxies. This gives us a benchmark against which other versions of the MDA can be compared. Since this is a best-case scenario, we have the least statistical uncertainty in the final result, making any systematic bias more apparent than for the subsequent MDAs. For the combined result with 249 events, the simulated value is contained within the support of the posterior distribution of $H_{0}$.

MDA0 demonstrates the importance of correctly accounting for GW selection effects. We are biased towards detecting sources which are nearer-by, and which are optimally orientated (closer to face-on). If an analysis is performed without taking into consideration the denominator $p(D_{\text{GW}}|H_{0})$ of Eq. (7), which corrects for this, the posterior density on $H_{0}$ converges to a value different from its simulated value of $70\,\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$. This can be seen in Fig. 2, where the dashed purple line shows the MDA0 combined posterior for all 249 events, neglecting GW selection effects entirely. We leave a detailed exploration of what level of accuracy in the GW selection function is required in order to move beyond 249 BNS-with-counterpart events, and simply note that in this case, it is sufficient enough that any biases which could affect the next stages of the MDA do not arise from the GW selection effects.

IV.2 MDA1: Complete Galaxy Catalog

The next more complex case is MDA1, where we assume no counterpart was observed, and resort to using a galaxy catalog. MDA1 uses a complete galaxy catalog containing all potential hosts – an optimistic scenario, in which EM selection effects do not need to be considered. The results with MDA1 already show a wider posterior distribution on $H_{0}$ ($68.91^{+1.36}_{-1.22}\,\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$) because of lack of certainty of the host galaxy (Fig. 3). The introduction of this uncertainty means that the contributions from each event will be smoothed out, depending on the size of the event’s sky localization and the number of galaxies within it. As can be seen in Fig. 4, there is a far higher proportion of events for which the likelihood is relatively broad and less informative, in comparison to Fig. 2. However, many events clearly have a small number of galaxies in their sky-area, and hence still show clear peaks.

IV.3 MDA2: Incomplete Galaxy Catalog

The next most complex scenario is the case where we have incomplete galaxy catalogs, limited by an apparent magnitude threshold. This gives us the first case where accounting for EM selection effects is important. To investigate this, we consider three galaxy catalogs, with apparent magnitude thresholds of $m_{\text{th}}=19.5$, $18$ and $16$, with respective completeness fractions of $75\%$, $50\%$ and $25\%$ in addition to the complete catalog for MDA1 (see III.3 for details). The combined 249-event posterior distributions on $H_{0}$ are shown in Fig. 5.

As the catalogs become less complete, the combined $H_{0}$ posterior becomes wider. This is because the probability that the host galaxy is inside the catalog decreases. The contribution from the galaxies within the catalog is reduced, and the uninformative contribution from the out-of-catalog term in Eq. (9) increases. This is visible in the individual likelihoods shown in Fig. 6, where instead of decreasing toward zero at high values of $H_{0}$, many of the individual likelihoods tend toward a constant. This is because, in the absence of EM data, and with the linear Hubble relation assumed in this work, the number of unobserved galaxies increases without limit as $d_{L}^{2}$. This is seen mostly for events at high distances (where the host has a lower probability of being in the catalog), or for well-localized events where there is no catalog support at the relevant redshifts within the event’s sky area. However, enough events are detected at low distances, where the catalogs are more complete and so provide informative redshift information, to produce an upper constraint on $H_{0}$.

We estimate $H_{0}=69.97^{+1.59}_{-1.50}$, $70.14^{+1.80}_{-1.67}$, and $70.14^{+2.29}_{-2.18}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ respectively for galaxy catalogs of 75%, 50%, and 25% completeness. See section IV.5 for a more in depth comparison of how galaxy catalog completeness affects posterior width.

Our exercise demonstrates that we need to know (or assess) the completeness of galaxy catalogs, and put in an appropriate out-of-catalog term in the analysis. For any of the MDA2 catalogs, if we assume that the galaxy catalog is complete, when in reality it is not, we get a posterior distribution on $H_{0}$ which is inconsistent with a value of 70 $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$. This is because the well-localized events for which the host is not inside the catalog do not have support for the correct value of $H_{0}$. In real catalogs, galaxy clustering might ensure that there are nearby bright galaxies in the catalog, partially mitigating this bias.

IV.4 MDA3: Luminosity Weighting

Until now we have considered all galaxies in our catalog to be equally likely to host a gravitational-wave source. In MDA3 we analyze the case described in Sec. III.4 where this is no longer true by constructing a galaxy catalog such that the probability of any single galaxy hosting a GW source is directly proportional to its luminosity. MDA3 includes the same EM selection effects as MDA2, in the sense that the catalog is magnitude limited. The completeness of this catalog, defined in terms of luminosity rather than numbers of galaxies, as defined in Eq. (15), is 50% out to 115 Mpc. This is indicative that approximately 50% of the detected GW events have host galaxies inside the catalog.

To investigate the importance of luminosity weighting, MDA3 was analyzed twice under different assumptions, given in Eq. (18). In the first, the analysis was matched to the known properties of the galaxy catalog, such that the probability of any galaxy hosting a GW event was proportional to its luminosity. In the second, we feigned ignorance and ran the analysis with the assumption that each galaxy was equally likely to be host to a GW event (as was true in MDAs 1 and 2). This allows us to determine the effect of ignoring galaxy weighting with this dataset. The combined $H_{0}$ posteriors for both cases are shown in Fig. 7. The estimated values of the Hubble constant are $70.83^{+3.55}_{-2.72}$  $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$(assuming hosts are luminosity weighted) and $69.50^{+4.20}_{-3.24}$  $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$(assuming equal weights). By weighting the host galaxies with the correct function of their luminosities, which happens to be known in this case, the constraint on $H_{0}$ improves — the uncertainty narrows by a factor of 1.2, compared to the case in which equal weights are assumed. Both results are consistent with the fiducial $H_{0}$ value of 70 $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$. In the limit of a far greater number of events, one might expect to see a bias emerge in the case in which the assumptions in the analysis do not match those with which the catalog was simulated. The luminosity weighting of host galaxies, by its very nature, increases the probability that the host galaxy is inside the galaxy catalog; assuming equal weighting gives disproportionate weight to the contribution that comes from beyond the galaxy catalog. However, for the 249 BNS events considered here, the final posteriors are too broad to be able to detect any kind of bias.

IV.5 Comparison between the MDAs

So far we have focused on individual event likelihoods and combined results for all 249 events. Our dataset also allows us to to assess the convergence for the combined Hubble posterior as we add events. We calculate the intermediate combined posteriors as a function of the number of events, and show the resulting convergence in Fig. 8. We plot the fractional $H_{0}$ uncertainty (defined here as the half-width of the 68.3% credible interval divided by $H_{0}$, $\Delta^{68.3\%}_{H0}/2H_{0}$), against the number of events we include in a randomly-selected group. The scatter between realizations of the group is indicated by the error bars, which encompass $68.3\%$ of their range. There is a considerable variation between different realizations, for the incomplete catalogs. For example, of the 100 realizations we used, for 25% completeness and 40 events, there are groups that give $\sim 10\%$ precision, but others that give $\sim 70\%$ precision.

With a sufficiently large number of events, we expect a $1/\sqrt{N}$ scaling of the uncertainty with the number of events Nissanke et al. (2013); Chen et al. (2018). To check whether this behavior is indeed true, we fit the results for each MDA to the expected scaling, obtaining the coefficient of $1/\sqrt{N}$ by maximizing its likelihood given the fractional uncertainties and their variances from the different realizations. The coefficient of the scaling is automatically dominated by the fractional uncertainties at large $N$ where the variances are small. We show this scaling for each MDA as a set of dashed lines in Fig. 8.

It can be seen that for each MDA, the results converge to the expected $1/\sqrt{N}$ scaling. The number of events required before this behavior is reached is dependent on the amount of EM information available on average for each event, in agreement with the results of Chen et al. (2018). The direct counterpart case is always on the trend after $\mathcal{O}(10)$ events, and shows a $\sim 18\%/\sqrt{N}$ convergence, comparable to and consistent with the results in Chen et al. (2018); Feeney et al. (2019). With the most complete galaxy catalogs, if the host galaxy is not directly identified it will take tens of events before this behavior is reached. However, even the least complete catalog (25%) appears to have reached this behavior by the time all 249 events are combined. It should be noted that as the catalogs for MDAs 1 and 2 were not simulated realistically, their low density relative to the density of the universe means that these numbers should not be taken as predictions of how fast $1/\sqrt{N}$ may be reached (except, perhaps, in the counterpart case, although one should bear in mind that even for that case, peculiar velocities and redshift uncertainties have been neglected). Even with a galaxy catalog which is $25\%$ complete, MDA2 gives a result which is only about a factor of $3$ times worse than the counterpart case.

As the density of galaxies in MDA3 was increased by 2 orders of magnitude over MDAs 1 and 2, the final posteriors cannot be directly compared between MDAs. However, by plotting the equivalent convergence figure for MDA3 (including the “known host” case as a reference, see Fig. 8), the impact of increasing the density of galaxies in the universe on the rate at which the posterior converges on the $1/\sqrt{N}$ behavior becomes clear. When there are more host galaxies, the results are overall less precise, and take longer to reach the $1/\sqrt{N}$ trend. As expected, using luminosity-weighting of potential host galaxies as an assumption in the analysis concentrates the probability to a smaller number of galaxies, leading to a more precise result.

IV.6 Limited Robustness Studies

Our results are expected to be sensitive to the luminosity distribution parameters — if one uses values of the Schechter function parameters $\alpha$ and $L^{*}$ in the analysis which are different from the ones used to simulate the galaxy catalogs, one would expect to end up with a bias in the results. With variations of these parameters within their current measurement uncertainties, we have however demonstrated that the resulting variation in the final result is small compared to the statistical uncertainties reached with the current set of MDAs. Furthermore we have also demonstrated that our results are robust against a small $\mathcal{O}(1)$ variation in the value of the telescope sensitivity threshold $m_{\text{th}}$.

V Conclusions and Outlook

The $H_{0}$ measurement using GW standard sirens has been demonstrated with recent events, including both the counterpart method for GW170817 Abbott et al. (2017c), and the galaxy catalog method Soares-Santos et al. (2019); Fishbach et al. (2019). These approaches are combined in the analysis of both BNS and BBH events from the first two observing runs of the advanced detector network Abbott et al. (2019b), using the method described in this paper. Future measurements will rely on a combination of counterpart and catalog methods, as appropriate for each new detected event, with catalog incompleteness playing an important role for the more distant, yet more common, BBHs. This paper outlines a coherent approach that tackles both of these scenarios, including treatment of selection effects in both GWs (due to the limited sensitivity of GW detectors) and EM (due to the flux-limitations of EM observing channels). We performed a series of MDAs to validate our method using up to 249 observed events. For each of the MDAs analyzed, the final posterior on $H_{0}$ is found to be consistent with the value of $\acs{H0}=70\,\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ used to simulate the MDA galaxy catalogs, demonstrating that our method can produce sufficiently unbiased results for treating these numbers of events, in our simulations.

GW selection effects are inherent in every version of the MDA and were corrected for by the term $p(D_{\text{GW}}|H_{0})$ in the denominator of Eq. (7). EM selection effects are addressed in MDAs 2 and 3 by the out-of-catalog terms containing $\bar{G}$ in Eq. (9). In both these MDAs, in spite of having an apparent magnitude-limited galaxy catalog, we are able to accurately infer $H_{0}$ without any bias. MDA2 further demonstrates our ability to account for missing host galaxies down to a level where only $25\%$ of events have hosts inside the catalog. Even in this case, we converge to the correct $H_{0}$ value, to the level of precision which could be reached by 249 events.

MDA3 demonstrates a clear tightening of the posterior distribution when we can assume that GW events trace the galaxy luminosities, compared to the case in which we treat all galaxies as equally likely hosts. The “uniform weights” analysis of MDA3 remains consistent with the simulated $H_{0}$ value. Hence we are unable to conclude whether an incorrect assumption would lead to a biased result, as one might expect. We used only 249 events for our MDAs. With enough events of comparable nature the bias would be detected. Future work will expand these studies to include a larger numbers of simulated GW events, and will be able to discern smaller sources of systematic effects.

Although the galaxy-catalog standard-siren measurement of $H_{0}$ is less precise than the counterpart measurement, it is still able to constrain $H_{0}$, but requires at least an order of magnitude more events in order to reach a comparable accuracy (in the most realistic case of MDA3). These MDAs have validated our method and implementation in simplified scenarios. However future work will be needed to improve on this in several directions, to test its applicability to BBHs (which are detectable out to much farther distances), realistic cosmology, and real galaxy catalogs Abbott et al. (2019b); Chen et al. (2018).

In both the counterpart and galaxy catalog cases, the lack of redshift uncertainties and peculiar velocities implies that the contributions from individual galaxies are a lot more precise than they would be in reality. Moreover, the simulated catalogs in MDAs 1 and 2 have a low density of galaxies compared to the universe, making them more informative than real catalogs. MDA3, with a galaxy density of 1 galaxy per 70 Mpc3, comes closest to the actual density of galaxies in the local universe of $\sim 1$ galaxy per 200 Mpc3 Dálya et al. (2018). In this scenario there is still a clear convergence towards the simulated $H_{0}$ value. In comparison to actual catalogs such as GLADE Dálya et al. (2018), the apparent magnitude threshold of $14$ is very low, so we expect a real catalog-only analysis to fall somewhere between MDAs 2 and 3. We caution the reader that with tens of events, the precision of results can vary by almost an order of magnitude depending on the particular realization of the detected population, before eventually converging to the expected $1/\sqrt{N}$ behaviour Nissanke et al. (2013); Chen et al. (2018). Analyzing more realistic catalogs will also require a sky-varying EM selection function, as the magnitude threshold varies significantly on the sky according to the design of particular surveys.

The galaxy distribution in these simulated catalogs is uniform in comoving volume. Although it has not been studied here, clustering of galaxies is expected to improve the constraint on $H_{0}$ (see, e.g. MacLeod and Hogan (2008); Chen et al. (2018)), since even when the host is not in the catalog, it is likely that there will be observed galaxies nearby.

The Advanced LIGO - Virgo second observing run Abbott et al. (2019a) has confirmed that BBH systems are detected at higher rates than BNSs. Since their greater mass allows them to be observed at much greater distances, where galaxy catalogs are incomplete, the catalog method including EM selection effects is particularly important. With the catalog of GW events expected to expand at an increasing rate in future observing runs, our analysis will evolve to meet the challenges that come with it, and give us the fullest picture of cosmology as revealed by gravitational waves.

Acknowledgements.

We thank members of the LIGO-Virgo Collaboration for valuable discussions pertaining to the writing of this paper, and in particular Nicola Tamanini for a careful reading of the manuscript. AG additionally thanks P. Ajith, Walter Del Pozzo, Anuradha Samajdar, and Chris Van Den Broeck for discussion at various stages of the work. RG, CM and JV are supported by the Science and Technology Research Council (grant No.ST/L000946/1). IMH is supported by the NSF Graduate Research Fellowship Program under grant DGE-17247915. IMH also acknowledges support from NSF Grant No. PHY-1607585. HQ is supported by Science and Technology Facilities Council (grant No.ST/T000147/1). AS thanks Nikhef for its hospitality and support from the Amsterdam Excellence Scholarship (2016-2018). HYC was supported by the Black Hole Initiative at Harvard University, through a grant from the John Templeton Foundation. MF and DEH were supported by NSF grant PHY-1708081. They were also supported by the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation. AG is supported by the research programme of the Netherlands Organisation for Scientific Research (NWO). DEH gratefully acknowledges support from the Marion and Stuart Rice Award. We are grateful for computational resources provided by the Leonard E Parker Center for Gravitation, Cosmology and Astrophysics at the University of Wisconsin-Milwaukee, and those provided by Cardiff University, and funded by an STFC grant supporting UK Involvement in the Operation of Advanced LIGO. This article has been assigned LIGO document number LIGO-P1900017.

Appendix A Detailed methodology

A.1 A note on luminosity weighting and redshift evolution

The probability for a galaxy to host a GW event is not uniform over all the galaxies present in the catalog. Indeed, brighter galaxies are supposed to have an higher star-formation rate and hence have an higher probability to host a GW event. Also galaxies at higher redshifts may be more likely to be hosts, as mergers are expected to be more frequent Madau and Dickinson (2014). Our prior belief for a galaxy at redshift $z$, sky position $\Omega$ and absolute and relative magnitudes $M,m$, to host a GW source $s$ can be expressed as

[TABLE]

where if we assume that $z,\Omega$ and $M$ are conditionally independent given $s,H_{0}$,

[TABLE]

In the last equation we used the explicit relation between apparent magnitude and $z,M$ and $H_{0}$. The probability $p(z)$ is the prior distribution of galaxies in the universe, taken to be uniform in comoving volume-time, $p(\Omega)$ is the prior on galaxy sky location, assumed uniform over the celestial sphere, and $p(M|H_{0})$ is the distribution of absolute magnitudes represented by the Schechter function. In the sections below we will show that the terms $p(s)$ and $p(s|H_{0})$ cancel out with other terms, and so their exact form does not need to be considered. $p(s|M,H_{0})$ can take the form

[TABLE]

We refer to the above equation as luminosity weighting. The term $p(s|z)$ represents the probability for the merger rate to depend on the redshift,

[TABLE]

For the MDAs in this paper with $z\ll 1$, we assume a constant rate model but a more generic model with $p(s|z)\propto(1+z)^{\lambda}$ can be used with detections at higher redshifts. This was the case of  Abbott et al. (2019b), for example, in which a $p(s|z)\propto(1+z)^{3}$ was assumed.

A.2 A detailed breakdown of the galaxy catalog case

This section presents a more detailed look into the galaxy catalog method presented in section II.3.1. The approach is summarized in Fig. 9, a network diagram which shows how each of the parameters of this extended derivation fit together and their dependencies on each other. The parameters which appear in this diagram, and in the following subsections, are defined in Table 3.

The subsections below provide derivations of the individual components of Eq. 9. Note that in the cases where the integration boundaries are not specified, they can be assumed to cover the full parameter space.

A.2.1 Likelihood when host is in catalog: $p(x_{\text{GW}}|G,D_{\text{GW}},H_{0})$

The likelihood in the case where the host galaxy is inside the galaxy catalog, $p(x_{\text{GW}}|G,D_{\text{GW}},H_{0})$, can be obtained from the marginalization over redshift, sky location, absolute magnitude and apparent magnitude. If we assume that the GW data, $x_{\text{GW}}$, is independent of the galaxy catalog $G$, $m$ and $M$ we can write

[TABLE]

The probability density function $p(z,\Omega,M,m|G,s,H_{0})$ is taken as a sum of delta functions with specific $z$, $\Omega$ and $m$ corresponding to the location of each galaxy in the catalog. This can be further factorized as

[TABLE]

where we have assumed again a relation between the apparent magnitude, redshift, $H_{0}$ and absolute magnitude. This allows us to integrate over the absolute magnitude in Eq. 20 and obtain

[TABLE]

Remembering that $p(z,\Omega,m|G)$ represents the distribution of the galaxies in the catalog, we can replace the integral above with a sum over the galaxies.

[TABLE]

where we have factorized $p(z_{i}|s)$ and $p(M(z_{i},m_{i},H_{0})|s)$, together with the term $p(s|z,\Omega,M(z,m,H_{0}),m,G,H_{0})$. Finally expanding the denominator $p(D_{\text{GW}}|G,s,H_{0})$ in the same way, we can recover the likelihood for the “in catalog” part of the galaxy catalog method.

[TABLE]

Notably, in the case the galaxies in the catalogs are provided along with their redshift uncertainties $p(z_{i})$, these can be implemented in the above equations as:

[TABLE]

A.2.2 Probability the host galaxy is in the galaxy catalog: $p(G|D_{\text{GW}},H_{0})$ and $p(\bar{G}|D_{\text{GW}},H_{0})$

The probability that the host galaxy is inside the galaxy catalog, given that a GW signal was detected, can be expressed as

[TABLE]

If we assume that the galaxy catalog is apparent magnitude-limited, such that only galaxies which are observed above some detection threshold are contained within it, we can approximate $p(G|z,\Omega,M,m,D_{\text{GW}},s,H_{0})$ as a Heaviside step around the detection threshold $m=m_{\text{th}}$.

[TABLE]

We now expand $p(z,\Omega,M,m|s,H_{0})$ as in Eq 17 and we obtain

[TABLE]

The term $p(D_{\text{GW}}|s,H_{0})$ can be expanded in a similar way and finally gives the probability for the host galaxy to be in the catalog.

[TABLE]

As the probabilities of being in the catalog and not in the catalog must be complementary, we have,

[TABLE]

A.2.3 Likelihood when host is not in catalog: $p(x_{\text{GW}}|\bar{G},D_{\text{GW}},H_{0})$

We follow an approach similar to the one presented in Appendix A.2.1. We expand

[TABLE]

The prior term, $p(z,\Omega,M,m|s,H_{0})$ can now be expanded as it was in Eq 17. Substituting this in, and utilizing a Heaviside step function to represent the galaxy catalog’s apparent magnitude threshold for $p(\bar{G}|z,\Omega,M,m,s,H_{0})$,

[TABLE]

Expanding the denominator, $p(D_{\text{GW}}|\bar{G},s,H_{0})$, in the same way gives an equivalent term,

[TABLE]

And substituting this back into Eq 31 finally gives,

[TABLE]

A.3 The catalog patch case

While in general the galaxy catalog method derived in A.2 was for use with a galaxy catalog which covers the entire sky, a small modification allows the use of catalogs which only cover a patch of sky, as long as the patch can be specified using limits in right ascension and declination. If we represent the sky area covered by the catalog as $\Omega_{\text{cat}}$, and the area outside the catalog as $\Omega_{\text{rest}}$, such that $\Omega_{\text{cat}}+\Omega_{\text{rest}}$ covers the whole sky, this can be written as follows:

[TABLE]

The first term is equivalent to the regular galaxy catalog case, but with limits on the integral over $\Omega$, while the second term has no $G$ and $\bar{G}$ terms, and covers the rest of the sky from redshift 0 to $\infty$.

A.4 Direct and pencil beam counterpart cases

The “direct” method assumes that the counterpart has been unambiguously linked to the host galaxy of the GW event, such that the redshift and sky location of that galaxy can be taken to be that of the GW event with certainty, see Eq. 10. Instead the numerator is calculated by evaluating the GW likelihood at the delta-function location of the counterpart in $z$ and $\Omega$, and the term in the denominator is evaluated as:

[TABLE]

for priors $p(z)$ and $p(\Omega)$ (note that this is independent of galaxy catalog data).

The “pencil-beam” method makes the assumption that while the sky location of the galaxy associated with the counterpart is that of the GW event, we may not make a direct association to a known galaxy but to a set of potential candidate hosts. We can use the EM constrained sky localization and therefore return to the question of whether the host is within or beyond the galaxy catalog. In this case, the likelihood takes the same form as in the galaxy catalog case, but evaluated along the line of sight of the candidate counterparts.

A.5 GW selection effects

Eq. 8 in section II.3 can be written as:

[TABLE]

where $p(D_{\text{GW}}|x_{\text{GW}},H_{0})$ is a binary quantity which is 1 if the SNR of $x_{\text{GW}}$ passes $\rho_{th}$, and 0 otherwise.

Looking at the individual components of Eq. 9 in their expanded forms (e.g. Eq. 25, 29 and 34), $p(D_{\text{GW}}|H_{0})$ only appears in an expanded form, where it is additionally conditioned on $z$ and $\Omega$. Calculating $p(D_{\text{GW}}|z,\Omega,H_{0})$ requires integrating over all realizations of GW events (detected and not), for a range of $z$, $\Omega$ and $H_{0}$ values, and applying a detection threshold ($\rho_{th}$) which all events must pass in order to be deemed detected.

Practically, Monte-Carlo integration can be used:

[TABLE]

where ${x_{\text{GW}}}_{i}$ corresponds to an event, the parameters of which have been randomly drawn from the prior distributions of parameters which affect an event’s detectability (mass, inclination, polarization, and sky location) and the event’s $\rho_{i}$ is calculated for specific values of $z$ and $H_{0}$.

[TABLE]

which gives a smooth function for $p(D_{\text{GW}}|z,\Omega,H_{0})$, which drops from 1 to 0 over a range of $z$, $\Omega$ and $H_{0}$ values.

For the MDA, we use a $\rho_{th}$ of 8 in each detector (assuming every event was detected by two detectors) and the 2016 PSD from (Singer et al., 2014, Fig. 1), and evaluate $p({D_{\text{GW}}}_{i}|{x_{\text{GW}}}_{i},z,\Omega,H_{0})$ for 5000 samples, such that the integral converges. For this analysis, we assume that the probability of detection is averaged over the course of the entire simulated observation period, such that the dependence of $D_{\text{GW}}$ on $\Omega$ is smeared out over the course of many days. We approximate this to mean that $p(D_{\text{GW}}|z,H_{0})$ is uniform over the sky (ignoring the mild declination dependence which would remain after the rotation of the Earth is taken into account). Fig. 10 shows how the probability of detection behaves as a function of $z$ for different values of $H_{0}$.

A.6 Prior mass distribution

An event’s detectability is dependent on its observed (redshifted) detector-frame mass, $M_{z}$, but priors on the mass refer to their source-frame mass. When calculating $p(D|H_{0})$ the masses are drawn from the priors on source mass, $p(M_{1},M_{2})$ and then converted to observed masses through the equation:

[TABLE]

However, when we use GW data in the form of posterior samples, the prior used to generate those is uniform on the redshifted mass, $M_{z}$ Veitch et al. (2015). Due to the way the MDA GW data was generated, with masses chosen on the detector-frame, rather than the source-frame, this was not something which had to be considered. With real GW data, as the redshift is linked directly to $H_{0}$, it is necessary to take into account the redshifting of the masses explicitly.

In general, when calculating $p(D|H_{0})$ for BBHs, the primary mass $M_{1}$ is drawn from a power-law with slope $\alpha$, between limits $[a,b]\ M_{\odot}$. The secondary mass, $M_{2}$ is drawn from a uniform distribution between $aM_{\odot}$ and $M_{1}$ Abbott et al. (2019c), to give (for $\alpha\neq-1$):

[TABLE]

This is related to the redshifted mass by the Jacobian:

[TABLE]

Substituting in our expression for $p(M_{1},M_{2})$:

[TABLE]

The factor of $(1+z)^{2}$ cancels in the numerator and denominator. As all redshift (and hence $H_{0}$) dependence has been removed, no correction is required for the differing priors. For the case in which $\alpha=-1$, it can be shown that all redshift dependence falls out as well, meaning that as long as the prior mass distribution takes the form of a power-law, no prior correction is required. This will not be the case for all mass distributions.