# Multiband Gravitational Wave Cosmography with Dark Sirens

Brian C. Seymour  0000-0002-7865-1052<EMAIL_ADDRESS>Hang Yu
0000-0002-6011-6190 Yanbei Chen  0000-0002-9730-9463 TAPIR, Walter Burke
Institute for Theoretical Physics, California Institute of Technology,
Pasadena, CA 91125, USA

###### Abstract

Gravitational waves might help resolve the tension between early and late
Universe measurements of the Hubble constant, and this possibility can be
enhanced with a gravitational wave detector in the decihertz band as we will
demonstrate in this study. Such a detector is particularly suitable for the
multiband observation of stellar-mass black hole binaries between space and
ground, which would significantly improve the source localization accuracy
thanks to a long baseline for timing triangulation, hence promoting the "dark
siren" cosmology. Proposed decihertz concepts include DECIGO/B-DECIGO, TianGO,
and others. We consider here the prospects of multiband observation of dark
siren binaries with a variety of network configurations. We find that a
multiband observation can uniquely identify a black hole binary to a single
galaxy to a cosmological distance, and thus a dark siren behaves as if it had
an electromagnetic counterpart. Considering only fully localized dark sirens,
we use a Fisher matrix approach to estimate the error in the Hubble constant
and matter density parameter. We find that a decihertz detector substantially
improves our ability to measure cosmological parameters because it enables
host galaxies to be identified out to a larger distance without the
systematics from statistical techniques based on comparing the population
distribution.

##  I Introduction

The Hubble-Lemaître constant $H_{0}$ describes the current expansion rate of
the universe. Currently, there is substantial deviation between Planck
measurements of the cosmic microwave background fluctuations [1] and SH0ES
measurements of Type 1a supernova with the distance ladder [2, 3]. Notably,
the Hubble tension between these early and late universe measurements differs
by at least $4\sigma$ [4, 5]. Moreover, the tension has occurred since the
first Planck results [6] and strengthened with time. It is important to
validate whether such Hubble tension truly exists or whether it is due to
astrophysical systematics because it could signify violation from the
$\Lambda$CDM concordance model [5, 7, 8]. One signature for departure from the
concordance model would be apparent redshift evolution $H_{0}$ [9].

If the Hubble tension is proven robust with further measurements, there are a
number of possible explanations with new physics. An extensive discussion of
these possibilities is given in a recent review [10], and a comparison between
many of such theories is given in [11]. These proposals can be classified
generally either as early-time modification of the sound horizon or late-time
modification of the Hubble expansion.

Let us first discuss new physics before recombination which would lower the
value of Hubble constant as measured by Planck. Early dark energy adds an
additional scalar field which acts like a cosmological constant and ends after
recombination [12, 13]. The time of matter-radiation equality can be shifted
by adding additional relativistic degrees of freedom with dark radiation [14,
15, 16] or neutrino self-interactions [17, 18]. Finally, there are proposals
that the photon mass cannot be perfectly measured due to the lifetime of the
universe from Heisenberg’s uncertainty principle, and this translates to
uncertainty on the Hubble constant [19, 20].

There are also a number of ways to create a smooth late time deformation in
$H(z)$ with unchanged CMB physics. These include phantom dark energy [21, 22],
running vacuum model [23, 24], phenomenologically emergent dark energy [25],
vacuum phase transition [26, 27], and a phase transition in dark energy [28,
29, 30]. Many of these modify the equation of state parameter of the dark
energy or change how $\Omega_{\Lambda}(z)$ evolves with redshift. Another way
to change the evolution of $H(z)$ is by introducing additional interactions.
These class includes well-known beyond GR theories such as Brans-Dicke gravity
[31], $f(R)$ gravity [32], and Galileon gravity [33]. Additionally, this
includes interacting dark energy [34, 35, 36, 37] where dark energy and dark
matter interact, and decaying dark matter [38, 39, 40] where dark matter
decays into an unknown dark radiation. Finally, the homogeneous and isotropic
assumption of $\Lambda$CDM can be broken with chameleon dark energy [41, 42,
43], cosmic voids [44], and inhomogeneous causal horizons [45].

The detection of gravitational waves (GW) can provide an independent late
universe measurement of the Hubble constant. By measuring the expansion rate
in the late universe, GW could be used as an independent measurement of the
Hubble constant from SH0ES. Furthermore, with a distribution of GW events at
low redshifts ($z\sim 0.0-0.5$), an anomalous evolution of the expansion rate
could be observed. In particular, the luminosity distance of the source can be
obtained from the measured gravitational waveform [46]. A Hubble constant
measurement can be readily attained from a standard siren: a binary neutron
star (BNS) merger with a coincident EM counterpart [46, 47]. With the optical
measurement of the redshift from EM followup and luminosity distance
measurement from the GW detector, one can directly measure the Hubble
constant. Indeed, the Hubble constant was measured with the BNS GW170817 [48]
and its corresponding EM counterpart [49, 50, 51].

However, only a small number of GW events are expected to be bright BNS
mergers with EM counterparts. The majority of observed GW events are binary
black hole (BBH) events without EM counterparts, which are thus known as dark
sirens. Notably, many BBH events have already been detected and cataloged [52,
53, 54, 55, 56, 57]. Dark sirens can measure the Hubble constant by
statistical techniques using galaxy catalogues [46, 58, 59, 60, 61, 62, 63,
64] and features in the mass distribution [65, 66, 67, 68, 69, 70]. These
statistical techniques can be further extended with realistic galaxy
clustering which provide improvements in identifying the redshift due to
galaxy density correlations [71, 72, 73, 74]. These statistical techniques
have been applied to the GWTC-3 catalog, and the Hubble constant is measured
as $H_{0}=68^{+13}_{-12}\,\,\text{km}\,\text{s}^{-1}\text{Mpc}^{-1}$ using
only dark sirens [75] at 68% credible level. By combining the statistical
method with the only standard siren GW170817, the Hubble constant is measured
as $H_{0}=68^{+8}_{-6}\,\,\text{km}\,\text{s}^{-1}\text{Mpc}^{-1}$. For
reference, GW170817 alone gives a Hubble constant value of
$H_{0}=69^{+17}_{-8}\,\,\text{km}\,\text{s}^{-1}\text{Mpc}^{-1}$ [48]. We need
to bear in mind that the statistical dark siren approach relies fundamentally
on population models so there is additional systematic uncertainties [75, 70].
In contrast, the Planck Hubble constant measurement was
$H_{0}=67.4^{+0.5}_{-0.5}\,\,\text{km}\,\text{s}^{-1}\text{Mpc}^{-1}$ [1] and
the SH0ES measurement was
$H_{0}=72.5^{+1.0}_{-1.0}\,\,\text{km}\,\text{s}^{-1}\text{Mpc}^{-1}$ [76]
which corresponds to the $4\sigma$ tension [4].

One new potential class of detector is one in the decihertz range (0.01 - 1
Hz), and such a detector may aid in measuring the Hubble constant. This
detector would lie in between the millihertz LISA band [77] and the 10 - 1000
Hz ground band. A decihertz detector has many advantages for measuring the
Hubble constant. First, it would provide early warning for BNS mergers which
would help guarantee EM identification [78, 79]. Second, a joint decihertz
detection would improve the parameter estimation for stellar mass BBH by
measuring their waves several years before they enter into the ground band
[78, 80]. Since statistical approaches to dark sirens are degraded by having
too many galaxies inside of the localization volume, having a better angular
localization will significantly help measure the cosmological parameters.
Furthermore, the fascinating possibility of a multiband detection exists where
a decihertz detector observes a BBH inspiral and then the ground based
detectors measure the merger and ringdown. A decihertz multiband detection has
been found to substantially improve parameter estimation accuracy [78]. By
combining decihertz and ground detectors, the detector network can uniquely
localize a BBH to a its host galaxy without any EM counterpart. While a ground
network can do this on its own [81], the addition of a decihertz detector will
significantly increase the range at which the BBH can be localized. In this
way, a multiband detection of a BBH can behave like a standard siren.

Right now, there are a number of existing and proposed gravitational wave
detectors. Advanced LIGO [82], Advanced Virgo [83], and KAGRA [84] are
operating ground based gravitational wave detectors and are second-generation
(2G) detectors. Following the 2G detectors, LIGO Voyager aims to maximize the
reach of existing LIGO observatory facilities by adding cryogenic operation,
heavier silicon test masses, and improved quantum squeezing [85, 86]. Einstein
Telescope [87] and Cosmic Explorer [88] are the 3rd generation of ground-based
detectors with planned arm lengths of 10 km and 40 km respectively which aim
to begin observation in the mid 2030s. 3G detectors can break the distance-
inclination degeneracy using higher order spherical harmonic modes which would
improve Hubble constant measurement [89].

At frequencies below $\sim 1$ Hz, detecting gravitational waves may best be
carried out in space due to technical challenges [90, 91]. LISA [77, 92],
TianQin [93], and Taiji [94, 95], are proposed space based detectors which
focus on the $\sim 10^{-3}-10^{-1}$ Hz bands. LISA can measure the Hubble
constant with dark sirens [96] with accuracy of 5% and may be able to measure
it with EMRIs [97] to an accuracy of 1% to 3%, though it is likely that ground
detectors will surpass this by the time it operational. Ref. [98] studied
measuring the Hubble constant measurement with TianQin and LISA/Einstein
Telescope. In the far future, there are proposals for a microhertz GW detector
[99]. At very low frequencies, it may be possible for a pulsar timing array to
measure the effect of a super massive black hole binary [100]. Furthermore,
there are a number of space based plans for a decihertz detector in the
$0.01-1$ Hz band. The Japanese detector DECIGO is an ambitious prospect that
consists of three clusters of interferometers with a 1000km arm length [101,
102, 103]. Big Bang Observer is concept like DECIGO by the European Space
Agency [104]. Previous work found that Big Bang Observer alone would provide
precision cosmological tests by measuring and localizing nearly every GW event
in the universe [105]. Recently, Ref. [106] studied the capabilities of DECIGO
and other decihertz detectors to measure the Hubble constant. B-DECIGO is a
planned pathfinder mission of DECIGO with a single interferometer and a 100 km
arm length [101, 102]. Finally, TianGO is a space based decihertz concept
which is designed with nearer-term technology [78, 107].

Figure 1: Comparison of detectors with a sample waveform. We plot the
instrumental sensitivities for TianGO (orange), LIGO Voyager (blue), LISA
(purple), aLIGO design sensitivity (red), Cosmic Explorer 2 (brown), Einstein
Telescope D (pink), TianQin (gray), and B-DECIGO (yellow). We also show a
sample TianGO waveform for a typical BBH merger (black) at $z=0.3$,
$\mathcal{M}_{c}=25M_{\odot}$, $q=1.05$, and $T_{\text{obs}}=5\text{ yr}$
assuming observed by a TianGO-like detector. On the top axis, we give the time
until merger.

For this analysis, we study how well we can measure the expansion rate of the
universe by measuring BBH with future ground detectors and decihertz concepts.
We consider two representative decihertz detectors: (i) TianGO in the LIGO
Voyager era, and (ii) B-DECIGO in the ET/CE era. TianGO is chosen because it
represents a possible near term decihertz detector. In such a timescale, it
would be operational in late 2020s/early 2030s and be working with the LIGO
Voyager network. B-DECIGO is a longer term prospect, which would be
operational in the late 2030s.

We forecast how well a dark siren can be localized with the Fisher matrix
formalism [108, 109] with both detector setups. If such a comoving volume
contains only one galaxy, we consider the dark siren to be localized. We
consider the case where localized events will have measured redshift due to
either spectroscopic follow-up or from a complete galaxy catalog. We find that
adding a decihertz detector to the network improves the range at which a dark
siren can be localized. We then constrain the Hubble constant and matter
density parameter by stacking the localized dark siren events together with
the BBH merger rate inferred by LIGO/Virgo/KAGRA [110]. We assume that the
Hubble constant and matter density are the Planck values and fix all other
cosmological parameters. Our study motivates how a decihertz detector can
complement the cosmological measurement capabilities of ground based
detectors.

The rest of the paper is organized as follows. In Sec. II, we describe the
observed strain in a space based detector, and we use the Fisher matrix
formalism to forecast the measurement uncertainties with a multiband
detection. In Sec. III, we describe how we stack localized events together and
the forecast dark siren constraints on the Hubble constant and matter density
parameter for various detector setups. We then conclude this work in Sec. IV.
Finally, App. A delves into the space-based waveform specifics, and App. B
justifies the conservative approach of considering only localized dark sirens.
Throughout the work, we use $G=c=1$.

##  II Measurement of a Binary Black Hole

###  II.1 TianGO Waveform

Let us first model the waveform in a space detector. TianGO is orbiting the
sun at an inclination of $60^{\circ}$, similar to the orbit of LISA [111].
Thus, there are two coordinate frames for the geometry of TianGO. We denote
the ecliptic frame to have basis
$(\boldsymbol{\hat{\bar{x}}},\boldsymbol{\hat{\bar{y}}},\boldsymbol{\hat{\bar{z}}})$
where $\boldsymbol{\hat{\bar{z}}}$ is normal to the orbit of the earth. The
frame with $(\boldsymbol{\hat{x}},\boldsymbol{\hat{y}},\boldsymbol{\hat{z}})$
is fixed on the center of TianGO with
$(\boldsymbol{\hat{x}},\boldsymbol{\hat{y}})$ oriented along its two arms. We
denote $\boldsymbol{\hat{N}}$ as the line of sight vector and
$\boldsymbol{\hat{L}}$ is the direction of binary angular momentum. We can
write the waveform as [80]

$\tilde{h}(f)=\Lambda(f)e^{-i\left[\Phi_{P}(f)+\Phi_{D}(f)\right]}\tilde{h}_{c}(f)\,,$
(1)

where $\tilde{h}_{c}(f)$ is the carrier waveform, $\Lambda(f)$ is the
amplitude in Eq. (2), $\Phi_{P}(f)$ is the polarization phase in Eq. (3), and
$\Phi_{D}(f)$ is the phase modulation due to Doppler effect in Eq. (4). The
carrier waveform is independent of the antenna patterns and only depends on
the intrinsic parameters $(\mathcal{M}_{z},q,D_{L},t_{c},\phi_{c})$ where
$\mathcal{M}_{z}=(1+z)\mathcal{M}_{c}$ is the detector frame chirp mass, $q$
is the mass ratio, $D_{L}$ is the luminosity distance, and $t_{c},\phi_{c}$
are the time and phase of coalescence. Because we wish to model the
gravitational waveform over the frequencies in both TianGO and Voyager, the
carrier waveform is modeled with a phenomenological waveform that combines
inspiral, merger and ringdown. Specifically, we use a `IMRPhenomD` waveform
[112, 113].

The notable difference for a space-based detector compared to a ground one is
that the orientation and location change with time. Thus, the amplitude and
polarization phase which characterize the antenna patterns acquire a frequency
dependence and are derived in [114, 115] for a space based detector. We write
them as

$\displaystyle\Lambda(f)$
$\displaystyle=\left[A_{+}^{2}F_{+}^{2}(f)+A_{\times}^{2}F_{\times}^{2}(f)\right]^{1/2}\,,$
(2) $\displaystyle\Phi_{P}$
$\displaystyle=\arctan\left[\frac{-A_{\times}F_{\times}(f)}{A_{+}F_{+}(f)}\right]\,.$
(3)

$F_{+,\times}(\phi_{S},\theta_{S},\psi_{S})$ are the detector beam pattern
coefficient where $(\phi_{S},\theta_{S})$ are the direction of
$\boldsymbol{\hat{N}}$ in the TianGO corotating frame and the barred ones
denote quantities in the ecliptic frame, and $\psi_{S}$ is the polarization
phase. The polarization amplitudes are
$A_{+}=1+(\boldsymbol{\hat{L}}\cdot\boldsymbol{\hat{N}})^{2}$ and
$A_{\times}=2\boldsymbol{\hat{L}}\cdot\boldsymbol{\hat{N}}$. Additionally,
there is a phase modulation due to the Doppler effect induced by the orbital
motion of the detector (which we have assumed to be a heliocentric one),

$\displaystyle\Phi_{D}(f)$ $\displaystyle=2\pi f\tau\,,$ (4)
$\displaystyle=2\pi
fR_{\text{AU}}\sin\bar{\theta}_{S}\cos\left(\bar{\phi}_{t}(f)-\bar{\phi}_{S}\right)\,,$
(5)

where $\tau=-\boldsymbol{d}\cdot\boldsymbol{\hat{N}}$, $\boldsymbol{d}$ is the
vector from barycenter to detector, $R_{\rm AU}$ is one AU, and
$\bar{\phi}_{t}(f)$ is the azimuthal location of the solar orbit of the
detector. The explicit expressions for $F_{+,\times}$,
$\boldsymbol{\hat{L}}\cdot\boldsymbol{\hat{N}}$, $\bar{\phi}_{t}(f)$ are given
in App. A. The ground waveforms are the same as Eq. (1), but they are
approximated as $f\rightarrow\infty$ for $\Lambda(f),\Phi_{P}(f),\Phi_{D}(f)$
since the antenna patterns are nearly constant while it is in band.

In Fig. 1, we plot a sample TianGO BBH waveform, along with the sensitivity of
some gravitational wave detectors. This waveform terminates on the left side
because of the 5 year observation time. It exhibits amplitude modulation
around $f\sim 2\cdot 10^{-2}$ Hz because TianGO’s orientation
$\boldsymbol{\hat{N}}$ is changing with a period of a year.

###  II.2 Parameter Estimation Background

Let us now describe how we use the Fisher analysis to estimate parameter
uncertainties. The Fisher matrix formalism provides a useful approximation to
parameter estimation in the high SNR limit [108, 109, 116]. We consider a
binary with parameters $\boldsymbol{\theta}^{a}$ and

$\boldsymbol{\theta}^{a}=\left(\ln\mathcal{M}_{z},q,\ln
D_{L},t_{c},\phi_{c},\bar{\phi}_{S},\bar{\theta}_{S},\bar{\phi}_{L},\bar{\theta}_{L}\right)\,.$
(6)

The variance for a specific parameter $\boldsymbol{\theta}^{a}$ is found on
the diagonal of the inverse of the Fisher matrix

$\Delta\boldsymbol{\theta}^{a}=\sqrt{\left(\Gamma^{-1}\right)_{aa}}\,,$ (7)

where the Fisher information matrix is defined as

$\Gamma_{ab}\equiv\left(\frac{\partial\tilde{h}}{\partial\boldsymbol{\theta}_{a}}\Big{|}\frac{\partial\tilde{h}}{\partial\boldsymbol{\theta}_{b}}\right)\,,$
(8)

and the waveform template $\tilde{h}(f,\boldsymbol{\theta})$ is a function of
frequency $f$ and parameters $\boldsymbol{\theta}$. The inner product between
two signals $\tilde{h}(f),\tilde{g}(f)$ is defined as

$\left(\tilde{g}\big{|}\tilde{h}\right)=4\,\text{Re}\int_{0}^{\infty}\frac{\tilde{g}^{\ast}(f)\tilde{h}(f)}{S_{n}(f)}df$
(9)

where $S_{n}(f)$ is the detector noise spectral density. In the case of a
network of detectors, we sum the individual Fisher matrix for each detector
$d$

$\left(\Gamma_{ab}\right)^{\text{net}}=\sum_{d}\Gamma_{ab}^{d}\,.$ (10)

###  II.3 Results from Parameter Estimation

To understand how a decihertz detector can enhance the parameter estimation of
a BBH, we examine the results obtained using TianGO with the HLI Voyager
network. The luminosity distance is defined by

$D_{L}(z)=\frac{1+z}{H_{0}}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}$ (11)

where

$E(z)\equiv\sqrt{\Omega_{m}\left(1+z\right)^{3}+\Omega_{\Lambda}}\,.$ (12)

For precision tests of cosmology, we are mostly interested in the luminosity
distance accuracy and volume localization. The size of the solid angle ellipse
$\Delta\Omega$ can be expressed by [114]

$\Delta\Omega=2\pi\sin\bar{\theta}_{S}\sqrt{\Sigma_{\bar{\phi}_{S}\bar{\phi}_{S}}\Sigma_{\bar{\theta}_{S}\bar{\theta}_{S}}-\left(\Sigma_{\bar{\theta}_{S}\bar{\theta}_{S}}\right)^{2}}\,.$
(13)

The uncertainty in comoving volume can be related to the angular uncertainty
by Eq. (28) of Ref. [117]

$\Delta V_{\mathrm{C}}=\frac{D_{L}^{2}}{(1+z)^{2}}\Delta\Omega\Delta D_{C}\,,$
(14)

where the comoving distance equals $D_{C}=D_{L}/\left(1+z\right)$. Using a
change of variables, the comoving volume uncertainty can be rewritten as

$\Delta
V_{\mathrm{C}}=\frac{D_{L}^{2}}{(1+z)^{3}+D_{L}H(z)\left(1+z\right)}\Delta\Omega\Delta
D_{L}\,,$ (15)

where $H(z)=H_{0}E(z)$.

Systematic errors beyond the detector sensitivity can degrade the accuracy of
the luminosity distance. The first of which is the gravitational lensing which
changes the luminosity distance. We use the fit from [118]

$\frac{\left(\Delta
D_{L}\right)_{\text{lens}}}{D_{L}}=0.066\left[\frac{1-(1+z)^{-0.25}}{0.25}\right]^{1.8}\,.$
(16)

Once a particular galaxy is identified, the peculiar velocity adds uncertainty
to the amount of cosmological redshift. The measured redshift is the sum of
the cosmological and Doppler redshift. We can express the peculiar velocity
systematic error as [119]

$\frac{\left(\Delta
D_{L}\right)_{\text{pv}}}{D_{L}}=\Big{|}1-\frac{\left(1+z\right)^{2}}{D_{L}H(z)}\Big{|}\sigma_{v}\,,$
(17)

where we have assumed $\sigma_{v}=200\text{ km s}^{-1}/c$. The relative
magnitude of this effect decreases rapidly with distance since the
cosmological redshift increases while the RMS peculiar velocity is
approximately constant.

Figure 2 gives the measurement accuracy for luminosity distance, angular
resolution, and spatial localization. We considered a binary of
$\mathcal{M}_{c}=25M_{\odot}$, $q=1.05$, a trailing angle between earth and
TianGO of $t_{a}=5^{\circ}$, and a $5$ year observation. The measurement
accuracy strongly depends upon inclination $\iota$ of the binary, in addition
to orientation of the detector network at merger. Therefore, we randomize over
$(\bar{\phi}_{S},\bar{\theta}_{S},\bar{\phi}_{L},\bar{\theta}_{L})$ in the
figure. The line represents the median measurement accuracy while the shaded
region contains $80\%$ of possible systems. While we use a 5 year observing
time for TianGO, the TianGO’s parameter estimation isn’t particularly
sensitive to the observing time as long as it’s above $\sim 1\text{ week}$ as
most of the SNR comes from frequencies above $0.1\text{ Hz}$ (see Fig. 1).

In the top part of Fig. 2, we show the fractional uncertainty in the
luminosity distance $\Delta D_{L}/D_{L}$ versus redshift. One can see that the
addition of TianGO doesn’t significantly improve the ability to measure the
luminosity distance compared with the HLI network. Most of the SNR from the
event comes from the ground network, so the addition of TianGO improves the
luminosity distance measurement by a factor of only 1.5111Note that we have
published a previous paper where we found that TianGO improved the luminosity
distance measurement of the HLI Voyager network (Fig. 3 and Fig. 13 of [78]).
There was an error in the space waveform code.. We also plot the lensing and
peculiar velocity systematic errors here. We see that the systematic error due
to peculiar velocity is only large enough to affect our measurement for very
close events. Meanwhile, the effect of lensing is negligible and can be
ignored in the future sections about cosmology.

In the middle panel of Fig. 2, we give the angular resolution $\Delta\Omega$
versus redshift. We see an angular resolution improvement by a factor of 20
for the addition of TianGO to the HLI Voyager network. The long baseline
between earth and TianGO is responsible for this upgraded sky localization
sensitivity.

Finally, let us describe the comoving volume localization in the bottom panel
of Fig. 2. We plot the comoving volume localization from Eq. (15), and find
that adding TianGO improves the comoving volume localization by a factor 30.
We use a comoving galaxy density of $n_{\text{gal}}=0.01\text{
gal}/\text{Mpc}^{3}$ [81]. This corresponds to the number density which are
about 25% as bright as the Milky Way. This is because the majority of the GW
are expected to come from galaxies at least this luminous [59]. If
$n_{\text{gal}}\Delta V_{C}<1$, we say the galaxy was localized. Using this
criterion, we find that HLI Voyager can localize galaxies up to $z\sim 0.15$,
while TianGO + HLI Voyager can localize them up to $z\sim 0.30$. Note that
error bands are large and asymmetric because the line-of-sight direction and
detector configuration greatly affect the measurement accuracy. For example, a
gravitational wave that is face-on to the ecliptic plane would be poorly
localized by TianGO since the Doppler term does not give any information,
while a GW coming edge-on will measure $\Delta\Omega$ well (and
correspondingly $\Delta V_{C}\sim D_{L}^{2}\Delta D_{L}\Delta\Omega$).

Figure 2:  Measurement accuracy for luminosity distance, angular resolution
and comoving volume localization versus redshift. We plot these measurement
uncertainties for TianGO + HLI Voyager (red), HLI Voyager (blue), and TianGO
(orange). Because we randomize over the angular extrinsic parameters
$(\bar{\phi}_{S},\bar{\theta}_{S},\bar{\phi}_{L},\bar{\theta}_{L})$, we plot
both the median measurement with the line and the shaded region where $80\%$
of binaries lie. We use $\mathcal{M}=25M_{\odot}$, $q=1.05$, $t_{a}=5^{\circ}$
and $T_{\text{obs}}=5\text{ yr}$. We use a galaxy number density per comoving
volume of $n_{\text{gal}}=0.01\text{ gal}/\text{Mpc}^{3}$ to convert comoving
volume localization to estimate our ability to identify the GW source.

###  II.4 Event Rate

To infer cosmological parameters, we stack all dark siren events that the
network can localize. Let us now estimate how many dark sirens can be
localized. First, the merger rate density $\mathcal{R}(z)$ describes the
number of mergers in a comoving volume per year. We model it with a power law
model and choose with $\kappa=2.7$ so that it corresponds to the Madau-
Dickinson star formation rate [120]

$\mathcal{R}(z)=\mathcal{R}_{0}\left(1+z\right)^{\kappa}\,.$ (18)

Since this is the source frame merger rate density, an additional factor of
$1/(1+z)$ is needed to convert time from the source frame to the detector
frame. Therefore, we write the detector-frame merger rate of sources with
$z<z_{m}$ as

$R_{\text{obs}}(z_{m})=\int_{0}^{z_{m}}\mathcal{R}(z^{\prime})\frac{1}{1+z^{\prime}}\frac{dV_{c}}{dz^{\prime}}dz^{\prime}\,,$
(19)

where

$\frac{dV_{c}}{dz}=\frac{4\pi}{H_{0}}\frac{d_{c}^{2}(z)}{E(z)}\,.$ (20)

We use the BBH merger rate $\mathcal{R}_{0}=20\text{ Gpc}^{-3}\text{yr}^{-1}$
and $\kappa=2.7$ which consistent with GWTC-3 [110].

In Fig. 3, we give the number of detections per year which can be fully
localized for HLI Voyager with and without TianGO. We see that TianGO will
nearly double the range at which a BBH can be localized to a single host. This
corresponds to an order of magnitude increase in localization rate.
Furthermore, since the localizations occur at higher redshift, we can probe
cosmological parameters beyond just the Hubble constant.

Figure 3: On the top, we plot the probability of an event being localized as a
function of redshift for HLI Voyager (green) and TianGO + HLI Voyager (blue).
We see that adding TianGO to the HLI Voyager network would nearly double the
range at which we can localize a dark siren event. On the bottom, we plot the
expected number of localizations in the comoving volume sphere. We use the
merger rate equal to the star formation rate (red dashed) from Eq. (19). We
find that the number of yearly localizations will increase by a factor of 10
by adding TianGO. This figure assumes the same binary parameters as Fig. 2,
but also uniformly samples the observation time $T_{\text{obs}}\in[0,5]\text{
yr}$.

##  III Cosmological Constraints

Figure 4: The confidence intervals for Hubble constant $H_{0}$ and matter
density parameter $\Omega_{m}$ from HLI Voyager (blue), TianGO + HLI Voyager
(red), Planck 2018 (yellow) [1], 2 CE 2’s + ET-D (cyan), and B-DECIGO + 2 CE
2’s + ET-D (pink). GW constraints come from Eq. (21) using only fully
localized BBH events during a five year observation. We use chirp mass
$\mathcal{M}_{c}=25M_{\odot}$ and merger rate density at the star formation
rate. One can see that adding TianGO to the Voyager network would improve
error in the measurement of the Hubble constant and the matter density
parameter. Moreover, the 3G ground network sees a similar improvement with the
addition of B-DECIGO assuming it is in a heliocentric orbit. We include the
forecasted cosmology constraints for other detector configurations in Tab. 1.
Figure 5:  Constraints on the expansion rate as a function of redshift for
various forecasted and current measurements at the 68% CL. We plot the
forecasted constraints on HLI Voyager (blue) and TianGO + HLI Voyager (red).
We also plot current expansion rate constraints from Planck 2018 (yellow) [1]
and from GW170817 (pink) [48]. We produce this plot assuming the Planck
parameters as the true values when computing the Fisher matrix, and
incorporate only uncertainty on $(H_{0},\Omega_{m})$ for the shaded regions.
Notice that three Voyagers can measure the expansion rate relatively
accurately below $z\sim 0.1$. Furthermore, adding the decihertz detector
TianGO enhances the ability for the expansion rate to be measured.

Given a set of gravitational wave observations, we wish to compute the
consistent values of the cosmology. Others have studied how to measure the
Hubble constant with dark sirens using statistical inference [59, 121, 70].
Currently, statistical methods are used because the LVK’s best localized BBHs
have comoving volume resolution of $\Delta V_{c}\sim 10^{5}\text{Mpc}^{3}$
[75] which has thousands of galaxies inside. Since our sources are well
localized, we can directly measure the redshift of each dark siren event from
the uniquely identified host galaxy. We demonstrate this in 2D with a mock
simulation in App. B that the likelihood function breaks down to the
particularly simple answer for well localized sources. We stress that our
approach of using the localization condition of $n_{\text{gal}}\Delta V_{C}<1$
is a conservative approach. This doesn’t require a catalogue since optical
telescopes can measure the redshift of the galaxy after the event.
Furthermore, galaxy clustering can improve the cosmology constraints [71].
Additionally, more massive galaxies are statistically more likely to be the
source of the GW, so this would further improve the ability to localize a GW
in the Bayesian approach. Under the localization assumption, a dark siren
(BBH) will behave like a bright one (i.e., BNS) for cosmology.

Let us now describe how to compute confidence intervals on the cosmology with
a set of dark siren observations. For a set of cosmological parameters
$\boldsymbol{H}=(H_{0},\Omega_{m},...)$, we can compute their confidence
intervals with a Fisher matrix

$\tilde{\Gamma}_{ij}=\sum_{\text{event }k}\frac{1}{\left(\Delta
D_{L}(z_{k})\right)^{2}}\frac{\partial D_{L}(z_{k},\boldsymbol{H})}{\partial
H_{i}}\frac{\partial D_{L}(z_{k},\boldsymbol{H})}{\partial H_{j}}\,,$ (21)

where we use the tilde $\tilde{\Gamma}$ to distinguish from the waveform
parameter estimation matrix used in the last section. Then the error in a
cosmological parameter is

$\Delta H_{i}=\sqrt{(\tilde{\Gamma}^{-1})_{ii}}\,.$ (22)

In the nearby universe, the Fisher matrix result reduces to $(\Delta
H_{0}/H_{0})^{2}=(\Delta D_{L}/D_{L})^{2}$

In Fig. 4, we plot the two sigma confidence intervals on the Hubble constant
and matter density parameter using only uniquely localized BBH events. We use
a five year observation period, and randomly pick
$(\bar{\phi}_{S},\bar{\theta}_{S},\bar{\phi}_{L},\bar{\theta}_{L})$. We use
$\mathcal{M}_{c}=25M_{\odot}$, $q=1.05$, a trailing angle of $5^{\circ}$, and
uniformly randomize the time until merger. The luminosity distance of the
events was sampled accordingly by Eq. (19). This corresponds to 2515 events
with $z<0.4$. There were 43 events localized by HLI Voyager alone and 476
events localized by HLI Voyager + TianGO.

The addition of TianGO substantially improves our ability to measure the
cosmology. Fig. 4 shows the improvement of using TianGO for measuring the
Hubble constant and matter density parameter. Because a multiband measurement
increases the distance we can uniquely localize a galaxy, we can measure the
matter density parameter much more accurately. HLI Voyager measures $H_{0}$ to
$1\%$ and $\Omega_{m}$ to $40\%$, and TianGO upgrades $H_{0}$ to $0.3\%$ and
$\Omega_{m}$ to $8\%$, while Planck measured $H_{0}$ to $0.8\%$ and
$\Omega_{m}$ to $2\%$. We also give the uncertainty ellipse for a possible 3G
network consisting of 2 CE2’s and 1 ET-D, and also we combine B-DECIGO with
the 3G network. We can see an improvement in both near-term and long-term
networks by adding a decihertz detector, particularly in the matter density
parameter since its effect is most pronounced at larger redshifts. Using the
covariance matrix containing $(H_{0},\Omega_{m})$, we can see how well the
expansion rate is measured as a function of redshift. In Fig. 5, we plot the
expansion rate $H(z)/(1+z)$ versus redshift where we shade the 68% CL regions.
We can see that gravitational wave detectors are measuring the redshift region
$z\sim 0.2$ well because the localizations are occurring here because most of
localized events are at this redshift. At large redshifts, the cosmic
expansion rate uncertainty grows because the matter density parameter is more
poorly measured. For reference, we also plot the constraints from GW170817 and
Planck 2018. Note that we only show Fig. 5 up to $z=0.6$ since we cannot
measure $\Omega_{\Lambda},\Omega_{k}$ well enough with localized BBH sources.

Finally, we estimate the constraints on the Hubble constant and matter density
parameter for various 2G to 3G detector networks in Tab. 1. Specifically, we
compare the cosmological constraints from localized dark sirens during a 5
year observation period. For the 3G detectors, we consider Cosmic Explorer 2
(CE2), and Einstein Telescope D (ET-D). We see that even with 2 CE2’s and
ET-D, TianGO improves the ability to measure the Hubble constant by a factor
of 2, and the matter density parameter by a factor of 3. This is because we
see a sizable improvement in the number of localized events.

For the long-term multiband case, we use a network consisting of B-DECIGO,
CE2, and ET-D. Because the orbit of B-DECIGO is still under discussion [122],
we placed it in a trailing 5∘ orbit like TianGO. We performed the same
analysis as in Section III. We find that the addition of B-DECIGO can improve
the cosmological measurement capabilities of the 3G detectors.

| $\Delta H_{0}/H_{0}$ | $\Delta\Omega_{m}$ | Localizations / 5 yr | Notes
---|---|---|---|---
| 3 V
---
(+ T)
| $1\times 10^{-2}$
---
($2\times 10^{-3}$)
| $1\times 10^{-1}$
---
($2\times 10^{-2}$)
| 43
---
(476)
| Voyager at Hanford, Livingston, India
---
sites.
| 1 CE2 + 1 ET-D
---
(+ T)
| $2\times 10^{-3}$
---
($6\times 10^{-4}$)
| $1\times 10^{-2}$
---
($3\times 10^{-3}$)
| 382
---
(1930)
| CE2 at Hanford,
---
ET-D at GEO-600 sites.
| 2 CE2 + 1 ET-D
---
(+ T)
| $1\times 10^{-3}$
---
($5\times 10^{-4}$)
| $6\times 10^{-3}$
---
($2\times 10^{-3}$)
| 843
---
(2410)
| CE2 at Hanford, Livingston.
---
ET-D at GEO-600 sites.
| 2 CE2 + 2V
---
(+ T)
| $1\times 10^{-3}$
---
($6\times 10^{-4}$)
| $9\times 10^{-3}$
---
($3\times 10^{-3}$)
| 556
---
(2211)
| CE2 at Virgo, India sites. Voyager at
---
Hanford, Livingston sites.
| 1 CE2 + 1 ET-D
---
(+ B-Decigo)
| $2\times 10^{-3}$
---
($5\times 10^{-4}$)
| $1\times 10^{-2}$
---
($2\times 10^{-3}$)
| 380
---
(4758)
| CE2 at Hanford, ET-D at GEO-600 sites.
---
B-Decigo placed in 5∘ trailing Heliocentric orbit
| 2 CE2 + 1 ET-D
---
(+ B-Decigo)
| $1\times 10^{-3}$
---
($3\times 10^{-4}$)
| $6\times 10^{-3}$
---
($1\times 10^{-3}$)
| 835
---
(5770)
| CE2 at Hanford, Livingston, ET-D at GEO-600 sites.
---
B-Decigo placed in 5∘ trailing Heliocentric orbit

Table 1: Dark siren constraints on the Hubble constant and matter density
parameter for various detector configurations. We use the same methodology as
for this table as in the rest of this paper. We find the Fisher matrix
confidence interval on the cosmological parameters by using only dark sirens
which are completely localized.

##  IV Conclusion

In this paper, we studied how a space-based decihertz detector can enhance the
sensitivity of a ground network for dark siren cosmological measurement. We
construct the case that these detectors will measure a significant number of
‘bright’ dark siren BBH – GW from which we can uniquely localize and uniquely
identify the host galaxy. We then use a Fisher matrix formalism to place
constraints on the cosmological parameters. We estimated how well the Hubble
constant and matter density parameter could be measured by BBH dark sirens
with a five year observation of TianGO plus three LIGO Voyagers. The result is
the multiband detection of dark sirens improves the measurement of the Hubble
constant by about a factor of 3. The larger redshift localized events allows
the matter density parameter to be resolved in the multiband case.

In the future, it would be interesting to extend our analysis to include dark
sirens which are non-uniquely identified, but are still well localized. Since
the fully localized criterion leaves out events with just a small number of
galaxies, information about the cosmology can still be extracted from these
events. Moreover, there are other effects which can improve the sensitivity
further, such as exploiting the clustering of galaxies to improve localization
[71] and weighting the galaxies by luminosity [73].

Measuring the cosmology with gravitational waves is easier when the host
galaxy is uniquely identified. The statistical dark siren approach is
degenerate with parameters such the merger rate evolution with redshift and
the BBH population model (as discussed in the GWTC-3 cosmology paper [75]).
Simultaneously measuring the cosmology and these population parameters can be
done by looking at the distribution of BBH events [123, 74, 70], but would
result in a less sensitive measurement of the cosmological parameters.
Otherwise, if these factors are not jointly measured, this would bias the
measurement of the Hubble constant [124, 70]. Consequently, a multiband
detection of dark sirens with uniquely identified hosts has the potential to
isolate the measurement of cosmological parameters from these population
parameters.

###### Acknowledgements.

B.S. acknowledges support by the National Science Foundation Graduate Research
Fellowship under Grant No. DGE-1745301. H.Y. acknowledges the support of the
Sherman Fairchild Foundation. Y.C. and B.S. acknowledge support from the
Brinson Foundation, the Simons Foundation (Award Number 568762), and by NSF
Grants PHY-2011961, PHY-2011968, PHY–1836809.

##  Appendix A Antenna Patterns of TianGO

The standard formula for the plus and cross antenna patterns of a detector is

$\displaystyle F_{+}$
$\displaystyle=\left(\frac{1+\cos^{2}\theta_{S}}{2}\right)\cos 2\phi_{S}\cos
2\psi_{S}-\cos\theta_{S}\sin 2\phi_{S}\sin 2\psi_{S}\,,$ (23) $\displaystyle
F_{\times}$ $\displaystyle=\left(\frac{1+\cos^{2}\theta_{S}}{2}\right)\cos
2\phi_{S}\cos 2\psi_{S}+\cos\theta_{S}\sin 2\phi_{S}\sin 2\psi_{S}\,.$ (24)

where $(\phi_{S},\theta_{S})$ are in the detector frame. We use the pycbc
detector class to get the ground based antenna patterns [125]. The antenna
patterns of a space detector are more complicated however, because the
detector has changing orientation. This means that the antenna patterns have
time dependence $F_{+,\times}(t)$, which we will use the time frequency
relation to find their frequency dependence.

To find the detector beam pattern coefficients, let us first describe the
geometry of the system. We have two coordinate systems: unbarred coordinates
$(\boldsymbol{\hat{x}},\boldsymbol{\hat{y}},\boldsymbol{\hat{z}})$ which
correspond to the individual detector and barred coordinates
$(\boldsymbol{\hat{\bar{x}}},\boldsymbol{\hat{\bar{y}}},\boldsymbol{\hat{\bar{z}}})$
in the ecliptic frame. The relationship between the orientation of the
detector frame and the ecliptic is

$\displaystyle\boldsymbol{\hat{x}}(t)$ $\displaystyle=-\frac{\sin
2\phi_{\mathrm{t}}}{4}\boldsymbol{\hat{\bar{x}}}+\frac{3+\cos
2\bar{\phi}_{\mathrm{t}}}{4}\boldsymbol{\hat{\bar{y}}}+\frac{\sqrt{3}}{2}\sin\bar{\phi}_{\mathrm{t}}\boldsymbol{\hat{\bar{z}}}\,,$
$\displaystyle\boldsymbol{\hat{y}}(t)$
$\displaystyle=\boldsymbol{\hat{z}}(t)\times\boldsymbol{\hat{x}}(t)\,,$
$\displaystyle\boldsymbol{\hat{z}}(t)$
$\displaystyle=-\frac{\sqrt{3}}{2}\left(\cos\bar{\phi}_{\mathrm{t}}\boldsymbol{\hat{\bar{x}}}+\sin\bar{\phi}_{\mathrm{t}}\boldsymbol{\hat{\bar{y}}}\right)+\frac{1}{2}\boldsymbol{\hat{\bar{z}}}\,,$
(25)

where the phase of TianGO in the ecliptic frame is equal to

$\bar{\phi}_{t}(f)=\frac{2\pi t(f)}{1\text{ yr}}-t_{a}\,,$ (26)

where $t_{a}$ is the trailing angle, and equal to $5^{\circ}$ for TianGO. The
time as a function of frequency is [78]

$t(f)=t_{c}-5\left(8\pi
f\right)^{-8/3}\mathcal{M}_{z}^{-5/3}\left[1+\frac{4}{3}\left(\frac{743}{336}+\frac{\mu}{M}x-\frac{32\pi}{5}x^{3/2}\right)\right]\,,$
(27)

where $\mu$ is the reduced mass and

$x=\left(\pi M_{z}f\right)^{2/3}\,.$ (28)

We can now write $(\phi_{S}(f),\theta_{S}(f),\psi_{S}(f))$ for the TianGO
detector using Eq. (25),

$\displaystyle\cos\theta_{S}(f)$
$\displaystyle=\frac{1}{2}\cos\bar{\theta}_{S}-\frac{\sqrt{3}}{2}\sin\bar{\theta}_{S}\cos\left(\bar{\phi}_{t}(f)-\bar{\phi}_{S}\right)\,,$
(29) $\displaystyle\phi_{S}(f)$
$\displaystyle=\bar{\phi}_{t}(f)+\arctan\left[\frac{\sqrt{3}\cos\bar{\theta}_{S}+\sin\bar{\theta}_{S}\cos\left(\bar{\phi}_{t}(f)-\bar{\phi}_{S}\right)}{2\sin\bar{\theta}_{S}\sin\left(\bar{\phi}_{t}(f)-\bar{\phi}_{S}\right)}\right]\,.$
(30)

The polarization phase of TianGO is

$\tan\psi_{S}(f)=\frac{\boldsymbol{\hat{L}}\cdot\boldsymbol{\hat{z}}-\left(\boldsymbol{\hat{L}}\cdot\boldsymbol{\hat{N}}\right)\left(\boldsymbol{\hat{z}}\cdot\boldsymbol{\hat{N}}\right)}{\boldsymbol{\hat{N}}\cdot\left(\boldsymbol{\hat{L}}\times\boldsymbol{\hat{z}}\right)}$
(31)

where

$\displaystyle\boldsymbol{\hat{N}}\cdot\boldsymbol{\hat{z}}$
$\displaystyle=\cos\theta_{S}(f)\,,$ (32)
$\displaystyle\boldsymbol{\hat{L}}\cdot\boldsymbol{\hat{z}}$
$\displaystyle=\frac{1}{2}\cos\bar{\theta}_{L}-\frac{\sqrt{3}}{2}\sin\bar{\theta}_{L}\cos\left[\bar{\phi}_{t}(f)-\bar{\phi}_{L}\right]\,,$
(33) $\displaystyle\boldsymbol{\hat{L}}\cdot\boldsymbol{\hat{N}}$
$\displaystyle=\cos\bar{\theta}_{L}\cos\bar{\phi}_{S}+\sin\bar{\theta}_{L}\sin\bar{\theta}_{S}\cos\left(\bar{\phi}_{L}-\bar{\phi}_{S}\right)\,,$
(34)
$\displaystyle\boldsymbol{\hat{N}}\cdot\left(\boldsymbol{\hat{L}}\times\boldsymbol{\hat{z}}\right)$
$\displaystyle=\frac{1}{2}\sin\bar{\theta}_{L}\sin\bar{\theta}_{S}\sin\left(\bar{\phi}_{L}-\bar{\phi}_{S}\right)$
$\displaystyle-\frac{\sqrt{3}}{2}$
$\displaystyle\cos\bar{\phi}_{t}(f)\left(\cos\bar{\theta}_{L}\sin\bar{\theta}_{S}\sin\bar{\phi}_{S}-\cos\bar{\theta}_{S}\sin\bar{\theta}_{L}\sin\bar{\phi}_{L}\right)$
$\displaystyle-\frac{\sqrt{3}}{2}$
$\displaystyle\sin\bar{\phi}_{t}(f)\left(\cos\bar{\theta}_{S}\sin\bar{\theta}_{L}\cos\bar{\phi}_{L}-\cos\bar{\theta}_{L}\sin\bar{\theta}_{S}\cos\bar{\phi}_{S}\right)\,.$
(35)

##  Appendix B Consistency of Statistical Method

In the statistical method, we wish to break the $z-D_{L}$ degeneracy by using
a galaxy catalog with the gravitational wave observation. We will use the
method described in a variety of sources [59, 73]. If we wish to constrain the
cosmological parameters $\boldsymbol{H}$ and have gravitational wave data
$d_{\text{GW}}$, then with Bayes theorem, we have

$p(\boldsymbol{H}|d_{\text{GW}})\propto
p(H_{0})p(d_{\text{GW}}|\boldsymbol{H})$ (36)

where

$\displaystyle p(d_{\text{GW}}|\boldsymbol{H})$
$\displaystyle=\frac{1}{\beta(\boldsymbol{H})}\int
p(d_{\text{GW}},D_{L},\phi_{S},\theta_{S},z|\boldsymbol{H})dD_{L}d\phi_{S}d\theta_{S}dz\,,$
(37) $\displaystyle=\frac{1}{\beta(\boldsymbol{H})}\int
p(d_{\text{GW}}|D_{L}(z,\boldsymbol{H}),\phi_{S},\theta_{S})p_{0}(z,\phi_{S},\theta_{S})d\phi_{S}d\theta_{S}dz\,.$
(38)

the first term in the integral is approximated with a multivariate Gaussian
distribution

$p(d_{\text{GW}}|D_{L}(z,\boldsymbol{H}),\phi_{S},\theta_{S})=N(D_{L}(z,\boldsymbol{H})-\hat{D}_{L},\sigma_{D_{L}}^{2})N(\phi_{S}-\hat{\phi}_{S},\sigma_{\phi_{S}}^{2})N(\theta_{S}-\hat{\theta}_{S},\sigma_{\theta_{S}}^{2})\,,$
(39)

where $N(x-\mu,\sigma^{2})$ is the probability density function of the normal
distribution, $(\hat{D}_{L},\hat{\phi}_{S},\hat{\theta}_{S})$ are the true
event parameters, and $(\sigma_{D_{L}},\sigma_{\phi_{S}},\sigma_{\theta_{S}})$
is given by the Fisher matrix analysis in Eq. (10). The second term in the
integral is the galaxy catalog

$p_{0}(z,\phi_{S},\theta_{S}|\boldsymbol{H})=\frac{1}{N_{\text{gal}}}\sum^{N_{\text{gal}}}_{i}N(z-z^{i},\sigma_{z_{i}}^{2})\delta(\phi_{S}-\phi_{S}^{i})\delta(\theta_{S}-\theta_{S}^{i})\,,$
(40)

where $\sigma_{z_{i}}$ is the variance due to the peculiar velocity. The
variables $(z^{i},\phi^{i},\theta^{i})$ are the mean redshift and angular
location of the $i$th galaxy, while unbarred variables are parameters.

The angular uncertainty is negligible and the distribution is replaced with a
Dirac delta function $\delta(\phi_{S}-\bar{\phi}_{S}^{i})$ and similarly for
$\theta_{S}$. Finally, the normalization $\beta(\boldsymbol{H})$ is

$\beta(\boldsymbol{H})=\int_{d_{\rm GW}>d_{\rm GW}^{\rm
th}}p(d_{\text{GW}},D_{L},\phi_{S},\theta_{S},z|\boldsymbol{H})dD_{L}d\phi_{S}d\theta_{S}dz\,dd_{\text{GW}}\,,$
(41)

where

$p(d_{\text{GW}},D_{L},\phi_{S},\theta_{S},z|\boldsymbol{H})=p(d_{\text{GW}}|D_{L}(z,\boldsymbol{H}),\phi_{S},\theta_{S})p_{0}(z,\phi_{S},\theta_{S})\,.$
(42)

and where $d_{\rm GW}^{\rm th}$ is the detection threshold. Note that Eq. (37)
reduces to the Fisher matrix confidence interval Eq. (21) on $\boldsymbol{H}$
if only one galaxy has nonvanishing likelihood. This reduction can be derived
by examining (38) in the case that there is only one galaxy inside the volume.
This happens when all other galaxies in the sum in
$p_{0}(z,\phi_{S},\theta_{S}|\boldsymbol{H})$ do not contribute to the
integral in Eq. (38).

Now, let us demonstrate the statistical method in 2D and examine its
convergence as a function of the number of galaxies inside the localization
region. We assume that $D_{L}=z/H_{0}$ and that the peculiar velocity
uncertainty is subdominant. Thus, we assume the peculiar velocity is a very
sharp Gaussian and absorb it into $\sigma_{D_{L}}$. If we call
$h=(H_{0})/(H_{0})_{\text{true}}$, the likelihood function is

$p(d_{\text{GW}}|h)=\frac{1}{\beta(h)}\frac{1}{N_{\text{gal}}}\sum_{i}N(\hat{D}_{L}-D_{L}^{i}(h),\sigma_{D_{L}}^{2})N(\hat{\phi}_{S}-\phi_{S}^{i},\sigma_{\phi_{S}}^{2})N(\hat{\theta}_{S}-\theta_{S}^{i},\sigma_{\theta_{S}}^{2})$
(43)

where $D_{L}^{i}(h)=z^{i}/H_{0}=z^{i}/[h(H_{0})_{\rm true}]$ and
$\sigma_{z}=\sigma_{D_{L}}(H_{0})_{\text{true}}$. In this 2D case,
$\beta(h)\propto h^{2}$.

If we need to stack events, we generalize Eq. (43) to be the product of the
likelihood function of each event222Technically, there is another factor
$p(N|h)$ in front of the product which depends on the intrinsic astrophysical
merger rate and comoving volume surveyed. It is discussed after Eq. (7) in
Ref. [73].,

$p(\left\\{d_{\text{GW}}\right\\}|h)=\prod_{\text{event
}e}^{N}p(\left(d_{\text{GW}}\right)_{e}|h)\,.$ (44)

If we assume a uniform prior on $h$, then
$p(h|\left\\{d_{\text{GW}}\right\\})\propto
p(\left\\{d_{\text{GW}}\right\\}|h)$. In Fig. 6, we plot the posterior on $h$
for 30 and 300 events. In this figure, we vary the angular resolution of the
events for each curve. We plot the median number of potential host galaxies
for the events. One can see that as events are nearly perfectly localized
$(n\rightarrow 0)$, the posterior on $h$ approaches the Fisher likelihood in
Eq. (21).

Due to the potential systematics possible in such an experiment, we list the
precise choices we used to make the plot. Our distance resolution was $\Delta
D_{L}/D_{L}=0.15z+10^{-2}$ and our angular resolutions varied between
$\Delta\phi_{S}=\frac{z}{1000}\deg$ to $\Delta\phi_{S}=100z\deg$. These scaled
with redshift linearly due to the SNR scaling of parameter measurement, while
the $10^{-2}$ is the same order as the peculiar velocity error (so a few close
events don’t dominate). We uniformly placed $3\times 10^{6}$ galaxies
thoughout the disc in the $z\in[0,2)$ ’redshift window’. For each event, we
randomly picked a galaxy with $z\in[0,1)$. The particular redshift window can
have a systematic effect on the statistical method [124], and we chose our
galaxy disc to be much bigger than the redshift window to avoid artificial
boundary effects.

Figure 6:  We plot the negative log likelihood of the posterior distribution
on the Hubble constant $p(h|\left\\{d_{\text{GW}}\right\\})$ in the 2D
simulation. Each curve with $n$ labeled is the median number of extra galaxies
in the localization region while the Fisher matrix constraint approximation
from Eq. (21) is also plotted (dashed red). Each curve in the plot corresponds
to picking a different angular resolution for the events. This shows that
measurement of sources with poor angular resolution will result in weaker
Hubble constant constraints due to the increased number of galaxies in the
localization region. We also see that with a higher number of events, the
likelihood distribution for $h$ tightens. Finally, we see that the Bayesian
approach reduces to the Fisher information estimate when there is a uniquely
identified galaxy. This is still a conservative estimation on how well we can
measure the cosmology as the information from $n>0$ systems is discarded.

## References

  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Riess  _et al._ [2021a] A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM, Astrophys. J. Lett. 908, L6 (2021a), arXiv:2012.08534 [astro-ph.CO] .
  * Riess  _et al._ [2016] A. G. Riess _et al._ , A 2.4% Determination of the Local Value of the Hubble Constant, Astrophys. J. 826, 56 (2016), arXiv:1604.01424 [astro-ph.CO] .
  * Verde  _et al._ [2019] L. Verde, T. Treu, and A. G. Riess, Tensions between the Early and the Late Universe, Nature Astron. 3, 891 (2019), arXiv:1907.10625 [astro-ph.CO] .
  * Camarena and Marra [2021] D. Camarena and V. Marra, On the use of the local prior on the absolute magnitude of Type Ia supernovae in cosmological inference, Mon. Not. Roy. Astron. Soc. 504, 5164 (2021), arXiv:2101.08641 [astro-ph.CO] .
  * Ade  _et al._ [2014] P. A. R. Ade _et al._ (Planck), Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Vagnozzi [2020] S. Vagnozzi, New physics in light of the $H_{0}$ tension: An alternative view, Phys. Rev. D 102, 023518 (2020), arXiv:1907.07569 [astro-ph.CO] .
  * Jackson [2007] N. Jackson, The Hubble Constant, Living Rev. Rel. 10, 4 (2007), arXiv:0709.3924 [astro-ph] .
  * Krishnan  _et al._ [2021] C. Krishnan, E. O. Colgáin, M. M. Sheikh-Jabbari, and T. Yang, Running Hubble Tension and a H0 Diagnostic, Phys. Rev. D 103, 103509 (2021), arXiv:2011.02858 [astro-ph.CO] .
  * Perivolaropoulos and Skara [2022] L. Perivolaropoulos and F. Skara, Challenges for $\Lambda$CDM: An update, New Astron. Rev. 95, 101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Schöneberg  _et al._ [2022] N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, The H0 Olympics: A fair ranking of proposed models, Phys. Rept. 984, 1 (2022), arXiv:2107.10291 [astro-ph.CO] .
  * Poulin  _et al._ [2019] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Early Dark Energy Can Resolve The Hubble Tension, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Agrawal  _et al._ [2019] P. Agrawal, F.-Y. Cyr-Racine, D. Pinner, and L. Randall, Rock ’n’ Roll Solutions to the Hubble Tension, (2019), arXiv:1904.01016 [astro-ph.CO] .
  * Buen-Abad  _et al._ [2018] M. A. Buen-Abad, M. Schmaltz, J. Lesgourgues, and T. Brinckmann, Interacting Dark Sector and Precision Cosmology, JCAP 01, 008, arXiv:1708.09406 [astro-ph.CO] .
  * Lancaster  _et al._ [2017] L. Lancaster, F.-Y. Cyr-Racine, L. Knox, and Z. Pan, A tale of two modes: Neutrino free-streaming in the early universe, JCAP 07, 033, arXiv:1704.06657 [astro-ph.CO] .
  * Green  _et al._ [2019] D. Green _et al._ , Messengers from the Early Universe: Cosmic Neutrinos and Other Light Relics, Bull. Am. Astron. Soc. 51, 159 (2019), arXiv:1903.04763 [astro-ph.CO] .
  * Cyr-Racine and Sigurdson [2014] F.-Y. Cyr-Racine and K. Sigurdson, Limits on Neutrino-Neutrino Scattering in the Early Universe, Phys. Rev. D 90, 123533 (2014), arXiv:1306.1536 [astro-ph.CO] .
  * Kreisch  _et al._ [2020] C. D. Kreisch, F.-Y. Cyr-Racine, and O. Doré, Neutrino puzzle: Anomalies, interactions, and cosmological tensions, Phys. Rev. D 101, 123505 (2020), arXiv:1902.00534 [astro-ph.CO] .
  * Capozziello  _et al._ [2020] S. Capozziello, M. Benetti, and A. D. A. M. Spallicci, Addressing the cosmological $H_{0}$ tension by the Heisenberg uncertainty, Found. Phys. 50, 893 (2020), arXiv:2007.00462 [gr-qc] .
  * Spallicci  _et al._ [2022] A. D. A. M. Spallicci, M. Benetti, and S. Capozziello, The Heisenberg Limit at Cosmological Scales, Found. Phys. 52, 23 (2022), arXiv:2112.07359 [physics.gen-ph] .
  * Alestas  _et al._ [2020] G. Alestas, L. Kazantzidis, and L. Perivolaropoulos, $H_{0}$ tension, phantom dark energy, and cosmological parameter degeneracies, Phys. Rev. D 101, 123516 (2020), arXiv:2004.08363 [astro-ph.CO] .
  * Di Valentino  _et al._ [2016] E. Di Valentino, A. Melchiorri, and J. Silk, Reconciling Planck with the local value of $H_{0}$ in extended parameter space, Phys. Lett. B 761, 242 (2016), arXiv:1606.00634 [astro-ph.CO] .
  * Banerjee  _et al._ [2021] S. Banerjee, D. Benisty, and E. I. Guendelman, Running Dark Energy and Dark Matter from Dynamical Spacetime, Bulg. J. Phys. 48, 117 (2021), arXiv:1910.03933 [gr-qc] .
  * Basilakos  _et al._ [2016] S. Basilakos, N. E. Mavromatos, and J. Solà, Starobinsky-like inflation and running vacuum in the context of Supergravity, Universe 2, 14 (2016), arXiv:1505.04434 [gr-qc] .
  * Li and Shafieloo [2019] X. Li and A. Shafieloo, A Simple Phenomenological Emergent Dark Energy Model can Resolve the Hubble Tension, Astrophys. J. Lett. 883, L3 (2019), arXiv:1906.08275 [astro-ph.CO] .
  * Di Valentino  _et al._ [2021] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, In the realm of the Hubble tension—a review of solutions, Class. Quant. Grav. 38, 153001 (2021), arXiv:2103.01183 [astro-ph.CO] .
  * Caldwell  _et al._ [2006] R. R. Caldwell, W. Komp, L. Parker, and D. A. T. Vanzella, A Sudden gravitational transition, Phys. Rev. D 73, 023513 (2006), arXiv:astro-ph/0507622 .
  * Banihashemi  _et al._ [2019] A. Banihashemi, N. Khosravi, and A. H. Shirazi, Ginzburg-Landau Theory of Dark Energy: A Framework to Study Both Temporal and Spatial Cosmological Tensions Simultaneously, Phys. Rev. D 99, 083509 (2019), arXiv:1810.11007 [astro-ph.CO] .
  * Banihashemi  _et al._ [2020] A. Banihashemi, N. Khosravi, and A. H. Shirazi, Phase transition in the dark sector as a proposal to lessen cosmological tensions, Phys. Rev. D 101, 123521 (2020), arXiv:1808.02472 [astro-ph.CO] .
  * Benevento  _et al._ [2020] G. Benevento, W. Hu, and M. Raveri, Can Late Dark Energy Transitions Raise the Hubble constant?, Phys. Rev. D 101, 103517 (2020), arXiv:2002.11707 [astro-ph.CO] .
  * Solà Peracaula  _et al._ [2019] J. Solà Peracaula, A. Gomez-Valent, J. de Cruz Pérez, and C. Moreno-Pulido, Brans–Dicke Gravity with a Cosmological Constant Smoothes Out $\Lambda$CDM Tensions, Astrophys. J. Lett. 886, L6 (2019), arXiv:1909.02554 [astro-ph.CO] .
  * D’Agostino and Nunes [2020] R. D’Agostino and R. C. Nunes, Measurements of $H_{0}$ in modified gravity theories: The role of lensed quasars in the late-time Universe, Phys. Rev. D 101, 103505 (2020), arXiv:2002.06381 [astro-ph.CO] .
  * Frusciante  _et al._ [2020] N. Frusciante, S. Peirone, L. Atayde, and A. De Felice, Phenomenology of the generalized cubic covariant Galileon model and cosmological bounds, Phys. Rev. D 101, 064001 (2020), arXiv:1912.07586 [astro-ph.CO] .
  * Wang and Meng [2005] P. Wang and X.-H. Meng, Can vacuum decay in our universe?, Class. Quant. Grav. 22, 283 (2005), arXiv:astro-ph/0408495 .
  * Di Valentino  _et al._ [2017] E. Di Valentino, A. Melchiorri, and O. Mena, Can interacting dark energy solve the $H_{0}$ tension?, Phys. Rev. D 96, 043503 (2017), arXiv:1704.08342 [astro-ph.CO] .
  * Yang  _et al._ [2020] W. Yang, S. Pan, R. C. Nunes, and D. F. Mota, Dark calling Dark: Interaction in the dark sector in presence of neutrino properties after Planck CMB final release, JCAP 04, 008, arXiv:1910.08821 [astro-ph.CO] .
  * Aljaf  _et al._ [2021] M. Aljaf, D. Gregoris, and M. Khurshudyan, Constraints on interacting dark energy models through cosmic chronometers and Gaussian process, Eur. Phys. J. C 81, 544 (2021), arXiv:2005.01891 [astro-ph.CO] .
  * Ichiki  _et al._ [2004] K. Ichiki, M. Oguri, and K. Takahashi, WMAP constraints on decaying cold dark matter, Phys. Rev. Lett. 93, 071302 (2004), arXiv:astro-ph/0403164 .
  * Bjaelde  _et al._ [2012] O. E. Bjaelde, S. Das, and A. Moss, Origin of Delta N_eff as a Result of an Interaction between Dark Radiation and Dark Matter, JCAP 10, 017, arXiv:1205.0553 [astro-ph.CO] .
  * Berezhiani  _et al._ [2015] Z. Berezhiani, A. D. Dolgov, and I. I. Tkachev, Reconciling Planck results with low redshift astronomical measurements, Phys. Rev. D 92, 061303 (2015), arXiv:1505.03644 [astro-ph.CO] .
  * Khoury and Weltman [2004a] J. Khoury and A. Weltman, Chameleon cosmology, Phys. Rev. D 69, 044026 (2004a), arXiv:astro-ph/0309411 .
  * Khoury and Weltman [2004b] J. Khoury and A. Weltman, Chameleon fields: Awaiting surprises for tests of gravity in space, Phys. Rev. Lett. 93, 171104 (2004b), arXiv:astro-ph/0309300 .
  * Cai  _et al._ [2021] R.-G. Cai, Z.-K. Guo, L. Li, S.-J. Wang, and W.-W. Yu, Chameleon dark energy can resolve the Hubble tension, Phys. Rev. D 103, 121302 (2021), arXiv:2102.02020 [astro-ph.CO] .
  * Lombriser [2020] L. Lombriser, Consistency of the local Hubble constant with the cosmic microwave background, Phys. Lett. B 803, 135303 (2020), arXiv:1906.12347 [astro-ph.CO] .
  * Fosalba and Gaztanaga [2020] P. Fosalba and E. Gaztanaga, Explaining Cosmological Anisotropy: Evidence for Causal Horizons from CMB data 10.1093/mnras/stab1193 (2020), arXiv:2011.00910 [astro-ph.CO] .
  * Schutz [1986] B. F. Schutz, Determining the Hubble Constant from Gravitational Wave Observations, Nature 323, 310 (1986).
  * Holz and Hughes [2005] D. E. Holz and S. A. Hughes, Using gravitational-wave standard sirens, Astrophys. J. 629, 15 (2005), arXiv:astro-ph/0504616 .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), A gravitational-wave standard siren measurement of the Hubble constant, Nature 551, 85 (2017a), arXiv:1710.05835 [astro-ph.CO] .
  * Nicholl  _et al._ [2017] M. Nicholl _et al._ , The Electromagnetic Counterpart of the Binary Neutron Star Merger LIGO/VIRGO GW170817. III. Optical and UV Spectra of a Blue Kilonova From Fast Polar Ejecta, Astrophys. J. Lett. 848, L18 (2017), arXiv:1710.05456 [astro-ph.HE] .
  * Coulter  _et al._ [2017] D. A. Coulter _et al._ , Swope Supernova Survey 2017a (SSS17a), the Optical Counterpart to a Gravitational Wave Source, Science 358, 1556 (2017), arXiv:1710.05452 [astro-ph.HE] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017b), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, (2021b), arXiv:2111.03606 [gr-qc] .
  * Venumadhav  _et al._ [2020] T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, New binary black hole mergers in the second observing run of Advanced LIGO and Advanced Virgo, Phys. Rev. D 101, 083030 (2020), arXiv:1904.07214 [astro-ph.HE] .
  * Olsen  _et al._ [2022] S. Olsen, T. Venumadhav, J. Mushkin, J. Roulet, B. Zackay, and M. Zaldarriaga, New binary black hole mergers in the LIGO–Virgo O3a data, (2022), arXiv:2201.02252 [astro-ph.HE] .
  * Nitz  _et al._ [2019] A. H. Nitz, C. Capano, A. B. Nielsen, S. Reyes, R. White, D. A. Brown, and B. Krishnan, 1-OGC: The First Open Gravitational-wave Catalog of Binary Mergers from Analysis of Public Advanced LIGO Data, ApJ 872, 195 (2019), arXiv:1811.01921 [gr-qc] .
  * Nitz  _et al._ [2020] A. H. Nitz, T. Dent, G. S. Davies, S. Kumar, C. D. Capano, I. Harry, S. Mozzon, L. Nuttall, A. Lundgren, and M. Tápai, 2-OGC: Open Gravitational-wave Catalog of Binary Mergers from Analysis of Public Advanced LIGO and Virgo Data, ApJ 891, 123 (2020), arXiv:1910.05331 [astro-ph.HE] .
  * Del Pozzo [2012] W. Del Pozzo, Inference of the cosmological parameters from gravitational waves: application to second generation interferometers, Phys. Rev. D 86, 043011 (2012), arXiv:1108.1317 [astro-ph.CO] .
  * Chen  _et al._ [2018] H.-Y. Chen, M. Fishbach, and D. E. Holz, A two per cent Hubble constant measurement from standard sirens within five years, Nature 562, 545 (2018), arXiv:1712.06531 [astro-ph.CO] .
  * Fishbach  _et al._ [2019] M. Fishbach, R. Gray, I. Magaña Hernandez, H. Qi, A. Sur, F. Acernese, L. Aiello, A. Allocca, and et al., A Standard Siren Measurement of the Hubble Constant from GW170817 without the Electromagnetic Counterpart, ApJ 871, L13 (2019), arXiv:1807.05667 [astro-ph.CO] .
  * Scelfo  _et al._ [2020] G. Scelfo, L. Boco, A. Lapi, and M. Viel, Exploring galaxies-gravitational waves cross-correlations as an astrophysical probe, J. Cosmology Astropart. Phys 2020, 045 (2020), arXiv:2007.08534 [astro-ph.CO] .
  * Finke  _et al._ [2021] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore, and M. Mancarella, Cosmology with LIGO/Virgo dark sirens: Hubble parameter and modified gravitational wave propagation, J. Cosmology Astropart. Phys 2021, 026 (2021), arXiv:2101.12660 [astro-ph.CO] .
  * Cigarrán Díaz and Mukherjee [2022] C. Cigarrán Díaz and S. Mukherjee, Mapping the cosmic expansion history from LIGO-Virgo-KAGRA in synergy with DESI and SPHEREx, MNRAS 511, 2782 (2022), arXiv:2107.12787 [astro-ph.CO] .
  * Mukherjee  _et al._ [2022] S. Mukherjee, A. Krolewski, B. D. Wandelt, and J. Silk, Cross-correlating dark sirens and galaxies: measurement of $H_{0}$ from GWTC-3 of LIGO-Virgo-KAGRA, arXiv e-prints , arXiv:2203.03643 (2022), arXiv:2203.03643 [astro-ph.CO] .
  * Taylor  _et al._ [2012] S. R. Taylor, J. R. Gair, and I. Mandel, Cosmology using advanced gravitational-wave detectors alone, Phys. Rev. D 85, 023535 (2012), arXiv:1108.5161 [gr-qc] .
  * Farr  _et al._ [2019] W. M. Farr, M. Fishbach, J. Ye, and D. E. Holz, A Future Percent-level Measurement of the Hubble Expansion at Redshift 0.8 with Advanced LIGO, ApJ 883, L42 (2019), arXiv:1908.09084 [astro-ph.CO] .
  * Mastrogiovanni  _et al._ [2021] S. Mastrogiovanni, K. Leyde, C. Karathanasis, E. Chassande-Mottin, D. A. Steer, J. Gair, A. Ghosh, R. Gray, S. Mukherjee, and S. Rinaldi, On the importance of source population models for gravitational-wave cosmology, Phys. Rev. D 104, 062009 (2021), arXiv:2103.14663 [gr-qc] .
  * Ezquiaga and Holz [2021] J. M. Ezquiaga and D. E. Holz, Jumping the Gap: Searching for LIGO’s Biggest Black Holes, Astrophys. J. Lett. 909, L23 (2021), arXiv:2006.02211 [astro-ph.HE] .
  * Ezquiaga and Holz [2022] J. M. Ezquiaga and D. E. Holz, Spectral Sirens: Cosmology from the Full Mass Distribution of Compact Binaries, Phys. Rev. Lett. 129, 061102 (2022), arXiv:2202.08240 [astro-ph.CO] .
  * Yu  _et al._ [2022] H. Yu, B. Seymour, Y. Wang, and Y. Chen, Uncertainty and bias of cosmology and astrophysical population model from statistical dark sirens, arXiv e-prints , arXiv:2206.09984 (2022), arXiv:2206.09984 [astro-ph.CO] .
  * MacLeod and Hogan [2008] C. L. MacLeod and C. J. Hogan, Precision of Hubble constant derived using black hole binary absolute distances and statistical redshift information, Phys. Rev. D 77, 043512 (2008), arXiv:0712.0618 [astro-ph] .
  * Nair  _et al._ [2018] R. Nair, S. Bose, and T. D. Saini, Measuring the Hubble constant: Gravitational wave observations meet galaxy clustering, Phys. Rev. D 98, 023502 (2018), arXiv:1804.06085 [astro-ph.CO] .
  * Gray  _et al._ [2020] R. Gray _et al._ , Cosmological inference using gravitational wave standard sirens: A mock data analysis, Phys. Rev. D 101, 122001 (2020), arXiv:1908.06050 [gr-qc] .
  * Mukherjee  _et al._ [2021] S. Mukherjee, B. D. Wandelt, S. M. Nissanke, and A. Silvestri, Accurate precision cosmology with redshift unknown gravitational wave sources, Phys. Rev. D 103, 043520 (2021), arXiv:2007.02943 [astro-ph.CO] .
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Constraints on the cosmic expansion history from GWTC-3, (2021c), arXiv:2111.03604 [astro-ph.CO] .
  * Riess  _et al._ [2021b] A. G. Riess _et al._ , A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km/s/Mpc Uncertainty from the Hubble Space Telescope and the SH0ES Team, (2021b), arXiv:2112.04510 [astro-ph.CO] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Kuns  _et al._ [2020] K. A. Kuns, H. Yu, Y. Chen, and R. X. Adhikari, Astrophysics and cosmology with a decihertz gravitational-wave detector: TianGO, Phys. Rev. D 102, 043001 (2020), arXiv:1908.06004 [gr-qc] .
  * Sedda  _et al._ [2020] M. A. Sedda _et al._ , The missing link in gravitational-wave astronomy: discoveries waiting in the decihertz range, Class. Quant. Grav. 37, 215011 (2020), arXiv:1908.11375 [gr-qc] .
  * Yu and Chen [2021] H. Yu and Y. Chen, Direct determination of supermassive black hole properties with gravitational-wave radiation from surrounding stellar-mass black hole binaries, Phys. Rev. Lett. 126, 021101 (2021), arXiv:2009.02579 [gr-qc] .
  * Chen and Holz [2016] H.-Y. Chen and D. E. Holz, Finding the One: Identifying the Host Galaxies of Gravitational-Wave Sources, (2016), arXiv:1612.01471 [astro-ph.HE] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Akutsu  _et al._ [2019] T. Akutsu _et al._ (KAGRA), KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector, Nature Astron. 3, 35 (2019), arXiv:1811.08079 [gr-qc] .
  * Adhikari  _et al._ [2017] R. Adhikari, N. Smith, A. Brooks, and et al., Ligo voyager upgrade concept,  (2017), LIGO Document T1400226.
  * Adhikari  _et al._ [2020] R. X. Adhikari _et al._ (LIGO), A cryogenic silicon interferometer for gravitational-wave detection, Class. Quant. Grav. 37, 165003 (2020), arXiv:2001.11173 [astro-ph.IM] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The Einstein Telescope: A third-generation gravitational wave observatory, Class. Quant. Grav. 27, 194002 (2010).
  * Reitze  _et al._ [2019] D. Reitze _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Borhanian  _et al._ [2020] S. Borhanian, A. Dhani, A. Gupta, K. G. Arun, and B. S. Sathyaprakash, Dark Sirens to Resolve the Hubble–Lemaître Tension, Astrophys. J. Lett. 905, L28 (2020), arXiv:2007.02883 [astro-ph.CO] .
  * Hall  _et al._ [2021] E. D. Hall _et al._ , Gravitational-wave physics with Cosmic Explorer: Limits to low-frequency sensitivity, Phys. Rev. D 103, 122004 (2021), arXiv:2012.03608 [gr-qc] .
  * Harms  _et al._ [2013] J. Harms, B. J. J. Slagmolen, R. X. Adhikari, M. C. Miller, M. Evans, Y. Chen, H. Müller, and M. Ando, Low-Frequency Terrestrial Gravitational-Wave Detectors, Phys. Rev. D 88, 122003 (2013), arXiv:1308.2074 [gr-qc] .
  * Auclair  _et al._ [2022] P. Auclair _et al._ (LISA Cosmology Working Group), Cosmology with the Laser Interferometer Space Antenna, (2022), arXiv:2204.05434 [astro-ph.CO] .
  * Luo  _et al._ [2016] J. Luo _et al._ (TianQin), TianQin: a space-borne gravitational wave detector, Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Hu and Wu [2017] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity, Natl. Sci. Rev. 4, 685 (2017).
  * Ruan  _et al._ [2020] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Taiji program: Gravitational-wave sources, Int. J. Mod. Phys. A 35, 2050075 (2020), arXiv:1807.09495 [gr-qc] .
  * Del Pozzo  _et al._ [2018] W. Del Pozzo, A. Sesana, and A. Klein, Stellar binary black holes in the LISA band: a new class of standard sirens, Mon. Not. Roy. Astron. Soc. 475, 3485 (2018), arXiv:1703.01300 [astro-ph.CO] .
  * Laghi  _et al._ [2021] D. Laghi, N. Tamanini, W. Del Pozzo, A. Sesana, J. Gair, S. Babak, and D. Izquierdo-Villalba, Gravitational-wave cosmology with extreme mass-ratio inspirals, Mon. Not. Roy. Astron. Soc. 508, 4512 (2021), arXiv:2102.01708 [astro-ph.CO] .
  * Zhu  _et al._ [2022] L.-G. Zhu, L.-H. Xie, Y.-M. Hu, S. Liu, E.-K. Li, N. R. Napolitano, B.-T. Tang, J.-d. Zhang, and J. Mei, Constraining the Hubble constant to a precision of about 1% using multi-band dark standard siren detections, Sci. China Phys. Mech. Astron. 65, 259811 (2022), arXiv:2110.05224 [astro-ph.CO] .
  * Sesana  _et al._ [2021] A. Sesana _et al._ , Unveiling the gravitational universe at $\mu$-Hz frequencies, Exper. Astron. 51, 1333 (2021), arXiv:1908.11391 [astro-ph.IM] .
  * Spallicci [2013] A. D. A. M. Spallicci, On the complementarity of pulsar timing and space laser interferometry for the individual detection of supermassive black hole binaries, Astrophys. J. 764, 187 (2013), arXiv:1107.5984 [gr-qc] .
  * Kawamura  _et al._ [2020] S. Kawamura _et al._ , Current status of space gravitational wave antenna DECIGO and BDECIGO, (2020), arXiv:2006.13545 [gr-qc] .
  * Sato  _et al._ [2017] S. Sato _et al._ , The status of DECIGO, J. Phys. Conf. Ser. 840, 012010 (2017).
  * Seto  _et al._ [2001] N. Seto, S. Kawamura, and T. Nakamura, Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space, Phys. Rev. Lett. 87, 221103 (2001), arXiv:astro-ph/0108011 .
  * Crowder and Cornish [2005] J. Crowder and N. J. Cornish, Beyond LISA: Exploring future gravitational wave missions, Phys. Rev. D 72, 083005 (2005), arXiv:gr-qc/0506015 .
  * Cutler and Holz [2009] C. Cutler and D. E. Holz, Ultra-high precision cosmology from gravitational waves, Phys. Rev. D 80, 104009 (2009), arXiv:0906.3752 [astro-ph.CO] .
  * Liu  _et al._ [2022] M. Liu, C. Liu, Y.-M. Hu, L. Shao, and Y. Kang, Dark-siren Cosmology with Decihertz Gravitational-wave Detectors, (2022), arXiv:2205.06991 [astro-ph.CO] .
  * Kuns [2019] K. A. Kuns,  _Future Networks of Gravitational Wave Detectors: Quantum Noise and Space Detectors_ , Ph.D. thesis, University of California, Santa Barbara (2019).
  * Finn [1992] L. S. Finn, Detection, measurement and gravitational radiation, Phys. Rev. D 46, 5236 (1992), arXiv:gr-qc/9209010 .
  * Cutler and Flanagan [1994] C. Cutler and E. E. Flanagan, Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral wave form?, Phys. Rev. D 49, 2658 (1994), arXiv:gr-qc/9402014 .
  * Abbott  _et al._ [2021d] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), The population of merging compact binaries inferred using gravitational waves through GWTC-3, (2021d), arXiv:2111.03634 [astro-ph.HE] .
  * Dhurandhar  _et al._ [2005] S. V. Dhurandhar, K. R. Nayak, S. Koshti, and J. Y. Vinet, Fundamentals of the LISA stable flight formation, Classical and Quantum Gravity 22, 481 (2005), arXiv:gr-qc/0410093 [astro-ph] .
  * Khan  _et al._ [2016] S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Husa  _et al._ [2016] S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal, Phys. Rev. D 93, 044006 (2016), arXiv:1508.07250 [gr-qc] .
  * Cutler [1998] C. Cutler, Angular resolution of the LISA gravitational wave detector, Phys. Rev. D 57, 7089 (1998), arXiv:gr-qc/9703068 .
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * Vallisneri [2008] M. Vallisneri, Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects, Phys. Rev. D 77, 042001 (2008), arXiv:gr-qc/0703086 .
  * Hogg [1999] D. W. Hogg, Distance measures in cosmology, (1999), arXiv:astro-ph/9905116 .
  * Hirata  _et al._ [2010] C. M. Hirata, D. E. Holz, and C. Cutler, Reducing the weak lensing noise for the gravitational wave Hubble diagram using the non-Gaussianity of the magnification distribution, Phys. Rev. D 81, 124046 (2010), arXiv:1004.3988 [astro-ph.CO] .
  * Gordon  _et al._ [2007] C. Gordon, K. Land, and A. Slosar, Cosmological Constraints from Type Ia Supernovae Peculiar Velocity Measurements, Phys. Rev. Lett. 99, 081301 (2007), arXiv:0705.1718 [astro-ph] .
  * Madau and Dickinson [2014] P. Madau and M. Dickinson, Cosmic Star Formation History, Ann. Rev. Astron. Astrophys. 52, 415 (2014), arXiv:1403.0007 [astro-ph.CO] .
  * Abbott  _et al._ [2021e] B. P. Abbott _et al._ (LIGO Scientific, Virgo), A Gravitational-wave Measurement of the Hubble Constant Following the Second Observing Run of Advanced LIGO and Virgo, Astrophys. J. 909, 218 (2021e), arXiv:1908.06060 [astro-ph.CO] .
  * Kawamura  _et al._ [2019] S. Kawamura _et al._ , Space gravitational-wave antennas DECIGO and B-DECIGO, Int. J. Mod. Phys. D 28, 1845001 (2019).
  * You  _et al._ [2021] Z.-Q. You, X.-J. Zhu, G. Ashton, E. Thrane, and Z.-H. Zhu, Standard-siren Cosmology Using Gravitational Waves from Binary Black Holes, ApJ 908, 215 (2021), arXiv:2004.00036 [astro-ph.CO] .
  * Trott and Huterer [2021] E. Trott and D. Huterer, Challenges for the statistical gravitational-wave method to measure the Hubble constant, (2021), arXiv:2112.00241 [astro-ph.CO] .
  * Biwer  _et al._ [2019] C. M. Biwer, C. D. Capano, S. De, M. Cabero, D. A. Brown, A. H. Nitz, and V. Raymond, PyCBC Inference: A Python-based parameter estimation toolkit for compact binary coalescence signals, Publ. Astron. Soc. Pac. 131, 024503 (2019), arXiv:1807.10312 [astro-ph.IM] .