# Calibrating systematic errors in the distance determination with the
luminosity-distance space large scale structure of dark sirens and its
potential applications

Pengjie Zhang1,2,3, Hai Yu1,2
1Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai, 200240, China
2Key Laboratory for Particle Astrophysics and Cosmology (MOE)/Shanghai Key
Laboratory for Particle Physics and Cosmology, China
3Tsung-Dao Lee Institute, Shanghai Jiao Tong University , Shanghai 200240,
China
Email to<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The cosmological luminosity-distance can be measured from gravitational wave
(GW) standard sirens, free of astronomical distance ladders and the associated
systematics. However, it may still contain systematics arising from various
astrophysical, cosmological and experimental sources. With the large amount of
dark standard sirens of upcoming third generation GW experiments, such
potential systematic bias can be diagnosed and corrected by statistical tools
of the large scale structure of the universe. We estimate that, by cross-
correlating the dark siren luminosity-distance space distribution and galaxy
redshift space distribution, multiplicative error $m$ in the luminosity
distance measurement can be constrained with $1\sigma$ uncertainty
$\sigma_{m}\sim 0.1$. This is already able to distinguish some binary black
hole origin scenarios unambiguously. Significantly better constraints and
therefore more applications may be achieved by more advanced GW experiments.

###### keywords:

Cosmology: large-scale structure of Universe, gravitational waves

##  1 Introduction

About 50 gravitational wave (GW) events of binary compact star coalescence
have been detected by advanced LIGO and advanced Virgo collaboration
(LVC)(Abbott et al., 2016, 2017b, 2017c, 2017d, 2019a; Abbott et al., 2020b;
Abbott et al., 2020c; The LIGO Scientific Collaboration et al., 2020a; Abbott
et al., 2020a; The LIGO Scientific Collaboration et al., 2020b). The total
number of detections will b dramatically increased by three orders of
magnitude by the third generation GW experiments such as the proposed Einstein
Telescope (ET111https://www.et-gw.eu/) and Cosmic Explorer (CE, Dwyer et al.
(2015); Abbott et al. (2017a); Reitze et al. (2019)). These experiments will
be able to detect binary black holes (BBHs) within our horizon (Dwyer et al.,
2015). Given the estimated BBH merger rate $R=53.2^{+58.5}_{-28.8}$ Gpc-3 yr-1
($90\%$ credibility, Abbott et al. (2019b)), $\sim 10^{5}$ BBH detections per
year are achievable. The angular localization accuracy will improve to $\la 1$
deg2, and the measurement accuracy of luminosity-distance through the standard
siren method(Schutz, 1986; Abbott et al., 2017e) will reach $\sigma_{D}/D\la
0.1$ (Zhao & Wen, 2018). The detection of higher spherical harmonics beyond
the quadruple mode (GW190814, Abbott et al. (2020c) & GW190412, The LIGO
Scientific Collaboration et al. (2020a)) and the analysis using the full
waveforms instead of the restricted waveforms(Arun et al., 2007; Van Den
Broeck & Sengupta, 2007; Ajith & Bose, 2009; Roy et al., 2019) will lead to
further significant improvement in the angular resolution and distance
determination.

Binary neutron star (BNS) GW events have weaker amplitude, so the redshift
coverage, angular resolution and distance determination accuracy by ET and CE
will be significantly worse than that of BBH GW events (Zhao & Wen, 2018;
Reitze et al., 2019). Nevertheless, more ambitious GW experiments such as Big
Bang Observer (BBO, Cutler & Harms (2006)) are capable of achieving $1\%$
distance measurement and arcminute angular localization. Since BNS event rate
is an order of magnitude higher (Abbott et al., 2017d, 2019a; Abbott et al.,
2020b), eventually we would expect $\sim 10^{6}$ GW events (mostly BNS events)
per year of known 3D positions all over the observable universe.

These upcoming 3D mappings of the universe through GW events will enable the
measurement and application of a new type of large scale structure (LSS) of
the universe, namely the luminosity-distance space LSS (Zhang, 2018a, b;
Namikawa, 2020; Libanore et al., 2020). Most of these GW events will not have
electromagnetic (EM) counterparts or host galaxies identified, so the usual
standard siren cosmology techniques (Schutz, 1986; Holz & Hughes, 2005; Cutler
& Holz, 2009; Abbott et al., 2017e; Soares-Santos et al., 2019; The LIGO
Scientific Collaboration et al., 2019) are not applicable. Nevertheless, these
dark standard sirens construct an excellent data set of LSS. Due to the fact
that gravitational waves are transparent while EM waves are not, the observed
GW event spatial distribution is immune to galactic and extragalactic
contaminations. Its spatial selection function is expected to be more
homogeneous and well understood. Together with the potentially full coverage
of cosmic volume within the horizon, the luminosity-distance space LSS is
promising to provide clean information on all linear LSS modes of our
universe. Furthermore, this LSS is spatially correlated with the usual LSS in
redshift space and various 2D projections (e.g. weak lensing), making them
highly complementary (Zhang, 2018a, b; Namikawa et al., 2016; Oguri, 2016;
Nair et al., 2018; Mukherjee & Wandelt, 2018; Soares-Santos et al., 2019; The
LIGO Scientific Collaboration et al., 2019).

The dark siren based LSS does not require EM identification of EM counterparts
or host galaxies. However, all promising applications of the dark siren based
LSS are built on unbiased measurement of the luminosity-distance $D^{\rm
obs}_{L}$, with respect to the cosmological luminosity-distance $D_{L}^{\rm
cos}$. Various experimental, astrophysical, and cosmological issues may cause
systematic bias in the distance measurement, analogous to redshift measurement
in EM observations. (1) The most obvious one is calibration error in the GW
strain amplitude (e.g. Sun et al. (2020)). Although it is expected to be small
(several percent for current experiments), independent verification is useful.
Calibration errors also affect the phase and therefore the inferred distance.
(2) One example of cosmological issues is modification to general relativity,
which may affect the GW propagation and/or the GW generation (e.g. Deffayet &
Menou (2007); Belgacem et al. (2018); Ezquiaga & Zumalacárregui (2018)). If
graviton has mass or if it can propagate in extra dimensions, the resulting
$D_{L}^{\rm obs}$ will be larger than the corresponding luminosity-distance of
EM counterpart, which equals to the luminosity-distance ($D^{\rm cos}$) of the
FRW universe. (3) One example of astrophysical issues is related to the
environment of GW events. This is motivated by the challenge in explaining the
origin of high mass BBH ($\ga 20M_{\odot}$) by stellar evolution models. Chen
et al. (2019) argued that the observed “high mass” BBHs may preferentially
reside within $\la 10$ Schwarzschild radius of the supermassive black hole of
host galaxies. These BBH events then suffer an extra gravitational redshift
$z_{\rm gra}$ and an extra Doppler redshift $z_{\rm Dop}$. The inferred black
hole masses will be overestimated by $(1+z_{\rm gra})(1+z_{\rm Dop})$. In this
scenario, the inferred luminosity-distance is overestimated by the same
factor, $D_{L}^{\rm obs}=D_{L}^{\rm cos}(1+z_{\rm gra})(1+z_{\rm Dop})$. This
effect can lead to a factor of $\sim 2$ bias in $D_{L}$. Other BBH environment
effects may also exist. Gas around the BBH orbit affects the orbital motion
and therefore affects the GW waveform (Chen & Shen, 2019). Although this
effect is likely negligible for high frequency observations such as LIGO, ET
and CE, this demonstrates the need to comprehensive investigation of BBH
environment and its impact on the $D_{L}$ determination.

Given the possible existence of significant systematic error in $D_{L}$,
standard siren cosmology will be significantly impacted. It is the purpose of
this paper to provide a model-independent method of diagnosing and correcting
such systematics. Significant systematic error in $D_{L}$ would cause abnormal
behavior in the angular correlation of GW dark sirens and galaxy
distributions. Similar behavior is well known in the context of LSS. One early
application is the SZ tomography (Zhang & Pen, 2001; Shao et al., 2011), where
cross-correlating wit galaxies help recover the redshift information of the
thermal Sunyaev Zel’dovich (SZ) effect. It is also widely used in the cross
calibration (Newman, 2008) and the self-calibration of galaxy photometric
redshift errors(Schneider et al., 2006; Zhang et al., 2010, 2017). Similar
techniques have been designed to identification of dark siren redshifts
(Namikawa et al., 2016; Oguri, 2016; Nair et al., 2018; Mukherjee & Wandelt,
2018; Zhang, 2018b). The question is to convert the directly measured cross
correlation into unbiased constraint on systematic error in $D_{L}$. We find
that this problem is essentially identical to identify the true redshift
distribution of GW dark sirens. For the later purpose, we have proposed the
$E_{z}$ estimator(Zhang, 2018b). It is essentially the same estimator that one
of the authors proposed in the context of SZ tomography (Zhang & Pen, 2001).
Zhang (2018b) proved that, it is unbiased in recovering the true redshift
distribution, insensitive to model assumptions on dark siren and galaxy
distribution. The same $E_{z}$ estimator, when applying to the current
problem, will not only diagnose the existence of systematic bias in $D_{L}$,
but also constrain it unbiasedly.

This paper is organized as follows. In §2, we briefly describe the $E_{z}$
estimator in the context of calibrating systematic error in $D_{L}^{\rm obs}$.
In §3 we estimate its constraining power. In term of the multiplicative error
$m$, we find that ET and CE, together with pre-existing galaxy surveys, will
be able to constrain with $\sigma_{m}\sim 0.1$. We discuss and summarize in
§4.

##  2 Cross-calibration with the $E_{z}$ estimator

Any error leading to $\langle D_{L}^{\rm obs}-D_{L}^{\rm cos}\rangle\neq 0$
belongs to systematic error. Following the common parametrization in cosmic
shear measurements (Heymans et al., 2006; Massey et al., 2007), we can
classify them into multiplicative errors and additive errors,

$\displaystyle D_{L}^{\rm obs}=D^{\rm cos}_{L}(1+m)+\epsilon\ .$ (1)

$\epsilon$ is the additive error. When $\langle\epsilon\rangle\neq 0$ or when
$\epsilon$ has spatial correlation (analogous to the galaxy intrinsic
alignment in cosmic shear), it is a systematic error. Otherwise it will be the
statistical error in the distance determination. For the current paper, we
will only consider the later case, as would be caused by instrumental noise of
GW experiments. Namely we consider $\epsilon$ as a random variable with zero
mean and pdf $p(\epsilon)$. For brevity, we assume $p(\epsilon)$ to be
Gaussian, with $\sigma_{\epsilon}=\sigma_{D_{L}^{\rm obs}}$. $m$ is the
multiplicative error in the luminosity-distance determination, and $\langle
m\rangle\neq 0$. As discussed in the previous section, it can have various
origins. For the astrophysical origin proposed by Chen et al. (2019),
$1+m=(1+z_{\rm gra})(1+z_{\rm Dop})$. Chen et al. (2019) argued that for BBHs
with mass $\ga 20M_{\odot}$, $m\sim 1$. But BBHs with lower mass may reside in
different environments and are free of this bias. Therefore later we will
consider two fiducial cases of $m=0$ and $m=1$.

With the presence of $m\neq 0$, the true redshift distribution of dark sirens
within a distance bin ($D_{1}<D^{\rm obs}_{L}<D_{2}$) is

$\displaystyle\bar{n}_{\rm GW}^{\rm cos}(z)$ $\displaystyle=$
$\displaystyle\int_{D_{1}}^{D_{2}}dD_{L}^{\rm obs}\int_{0}^{\infty}dD^{\rm
cos}_{L}\delta_{D}\left(z-z(D^{\rm cos}_{L})\right)$ (2) $\displaystyle\times
p(D^{\rm cos}_{L}|D_{L}^{\rm obs})n_{\rm GW}^{\rm obs}(D_{L}^{\rm obs})\ .$

Here $n_{\rm GW}^{\rm obs}(D_{L}^{\rm obs})$ is the number of dark sirens per
unit $D_{L}^{\rm obs}$ interval. $p(D^{\rm cos}_{L}|D_{L}^{\rm obs})$ is the
PDF of $D^{\rm cos}_{L}$ given $D_{L}^{\rm obs}$. It is determined by the
statistical error ($\epsilon$) probability distribution $p(\epsilon)$ and $m$,
where $\epsilon=D_{L}^{\rm obs}-D^{\rm cos}_{L}(1+m)$. Fluctuations in the
dark siren surface density is then

$\displaystyle\delta_{{\rm GW}}^{\Sigma}(\hat{\theta})$ $\displaystyle=$
$\displaystyle\bar{\Sigma}_{{\rm GW}}^{-1}\left[\int_{0}^{\infty}\delta_{\rm
GW}(z,\hat{\theta})\bar{n}^{\rm cos}_{{\rm GW}}(z)dz\right]\ .$ (3)

Here the mean surface number density

$\displaystyle\bar{\Sigma}_{{\rm GW}}$ $\displaystyle=$
$\displaystyle\int_{D_{1}}^{D_{2}}\bar{n}_{\rm GW}^{\rm obs}(D_{L}^{\rm
obs})dD_{L}^{\rm obs}=\int_{0}^{\infty}\bar{n}^{\rm cos}_{{\rm GW}}(z)dz\ .$
(4)

$\delta_{\rm GW}(z,\hat{\theta})$ is the number overdensity of dark sirens at
cosmological redshift $z$ and angular direction $\hat{\theta}$.

These dark sirens have angular cross-correlation with galaxies. We have the
freedom to weigh the galaxy number distribution $n_{g}(z)$ with an arbitrary
weighting function $W_{g}(z)$. For the weighted galaxy distribution,

$\displaystyle\delta_{\rm g}^{\Sigma}(\hat{\theta})$ $\displaystyle=$
$\displaystyle\bar{\Sigma}^{-1}_{\rm g}\left[\int_{0}^{\infty}\delta_{\rm
g}(z,\hat{\theta})\bar{n}_{\rm g}(z)W_{g}(z)dz\right]\ .$ (5)

Here the mean galaxy surface number density $\bar{\Sigma}_{\rm
g,j}=\int_{0}^{\infty}\bar{n}_{\rm g}(z)W_{g}(z)dz$. The GW-galaxy cross and
galaxy auto angular power spectra at multipole $\ell$ are

$\displaystyle C^{\rm GW-g}(\ell)$ $\displaystyle=$ $\displaystyle\frac{\int
P_{\rm GW-g}(k,z)\bar{n}^{\rm cos}_{\rm GW}\bar{n}_{\rm
g}W_{g}\chi^{-2}\frac{dz}{d\chi}dz}{\bar{\Sigma}_{\rm GW}\bar{\Sigma}_{\rm
g}}\ ,$ (6) $\displaystyle C^{\rm g}(\ell)$ $\displaystyle=$
$\displaystyle\frac{\int P_{\rm g}(k,z)(\bar{n}_{\rm
g}W_{g})^{2}\chi^{-2}\frac{dz}{d\chi}dz}{\bar{\Sigma}_{\rm g}^{2}}\ .$ (7)

Here $\chi$ is the comoving radial distance. $P_{\rm g}(k,z)$ and $P_{\rm
GW-g}(k,z)$ are the 3D galaxy and galaxy-GW host galaxy power spectrum
respectively, with $k=\ell/\chi(z)$ in the above integrals. The above
expressions adopt a flat universe cosmology. But the methodology holds for
curved universe as well.

The true redshift distribution (or $m$) of GW dark standard sirens can be
obtained by maximizing the $E_{z}$ estimator (Zhang, 2018b),

$\displaystyle E_{z}(\ell|W_{g})=\frac{C_{\rm
GW-g}(\ell|W_{g})}{\sqrt{C_{g}(\ell|W_{g})}}\ .$ (8)

Here we highlight the dependence of $E_{z}$, $C_{\rm GW-g}$ and $C_{g}$ on the
weighting function $W_{\rm g}$. By varying $W_{g}$ to find the maximum of
$E_{z}$, unbiased estimation of the true redshift distribution of dark sirens
can be obtained (refer to Zhang (2018b) for the proof). The reason is that
$E_{z}$ reaches maximum when the cross-correlation coefficient between surface
overdensities of GW and weighted galaxies reaches maximum. This is naturally
achieved when the two distributions match each other in redshift. For $W_{g}$
of free functional form, the numerical variation can be challenging. However,
in the context of the current paper, the situation is simpler and we can take
a constrained form of $W_{\rm g}$. Since $m$ is the only unknown parameter,
the form of $W_{g}$ can be fixed up to an unknown fitting parameter $m^{\rm
fit}$ (namely replacing $m$ in Eq. 2 with $m^{\rm fit}$),

$\displaystyle W_{\rm g}(z|m^{\rm fit})=\frac{\bar{n}^{\rm cos}_{\rm
GW}(z|m^{\rm fit})}{\bar{n}_{g}(z)}\ .$ (9)

When $m^{\rm fit}=m^{\rm true}$, the weighted galaxy redshift distribution
$\bar{n}_{\rm g}W_{\rm g}\propto\bar{n}^{\rm cos}_{\rm GW}$. Namely the
weighted galaxy redshift distribution will match that of the GW redshift
distribution exactly. This will result in the maximum $E_{z}$. So the form of
Eq. 9 includes returns unbiased estimation of $m$. This is also numerically
shown in Fig. 1. Therefore we will adopt this constrained form of $W_{\rm g}$.

For brevity, hereafter we neglect the superscript “fit” in $m^{\rm fit}$ where
it does not cause confusion with the fiducial value $m^{\rm fid/true}$. The
bestfit value of $m$ is obtained when $\partial E_{z}/\partial m=0$. Then we
take the data as $\partial E_{z}/\partial m$. The posterior probability of $m$
is

$\displaystyle p(m)$ $\displaystyle\propto$
$\displaystyle\exp\left(-\frac{1}{2}\Delta\chi^{2}\right)\ ,$
$\displaystyle\Delta\chi^{2}$ $\displaystyle=$
$\displaystyle\sum_{\ell}\left(\frac{\partial E_{z}}{\partial
m}C^{-1}\frac{\partial E_{z}}{\partial m}\right)\ .$ (10)

Here $C$ is the covariance matrix of $\partial E_{z}/\partial m$ at multipole
$\ell$, whose full expression is given in Zhang (2018b). Statistical errors in
$\partial E_{z}/\partial m$ of different $\ell$ do not correlate, so we can
sum over contributions from each $\ell$ bin to obtain the total
$\Delta\chi^{2}$.

Figure 1: $E_{z}$ as a function of the fitting parameter $m^{\rm fit}$, for
the luminosity bin $D_{L}H_{0}/c\in(1.19,1.98)$. The fiducial $m$ is adopted
as $m=0$ and $m=1$, respectively. Indeed, $E_{z}(m^{\rm fit})$ peaks at
$m^{\rm fit}=m$. Namely it provides unbiased estimation of $m$. We have
normalized its amplitude with $E_{z}(m^{\rm fit}=m)$.

##  3 Error forecast

We apply the $E_{z}$ estimator to evaluate its performance in constraining
$m$, for third generation GW experiments such as ET and CE. We adopt the
fiducial cosmology as the $\Lambda$CDM cosmology with $\Omega_{m}=0.268$,
$\Omega_{\Lambda}=1-\Omega_{m}$, $\Omega_{b}=0.044$, $h=0.71$,
$\sigma_{8}=0.83$ and $n_{s}=0.96$. We consider the fiducial $m=1$, as would
be expected in the scenario of (Chen et al., 2019). The first question is that
whether the data alone allows model-independent discrimination against the
case of $m=0$. If yes, then we need to verify that the bestfit indeed gives
$m=1$ and need to quantify $\sigma_{m}$.

Implementing the $E_{z}$ method requires deep and wide galaxy surveys with
accurate redshift distribution. By the time of these GW experiments, stage IV
spectroscopic redshifts surveys of $\sim 10^{8}$ galaxies such as CSST (Zhan,
2018; Gong et al., 2019), DESI (DESI Collaboration et al., 2016), Euclid
(Laureijs et al., 2011), PFS (Takada et al., 2014), WFIRST (Spergel et al.,
2015), and intensity mapping surveys (e.g. SKA(Santos et al., 2015) and
SPHEREx (Doré et al., 2014; Dore, 2018),) will likely be completed and satisfy
the need. More advanced surveys with $\sim 10^{9}$ galaxies with spectroscopic
redshifts and deeper redshift coverage may also be available (e.g. Abdalla et
al. (2015); Dodelson et al. (2016); Wang et al. (2019); Ferraro & Wilson
(2019)). We also expect photometric redshift of imaging surveys such as LSST
will improve to the required accuracy for implementing the $E_{z}$ method.
Therefore for a bin width of $\Delta z\sim 0.2$, we expect $N_{g}$ of the
order $10^{7-8}$.

For the GW experiments, we target at third generation experiments like ET and
CE. Unlike more advanced GW experiments of percent-level accuracy in
luminosity distance determination and arcminute angular resolution for both
BBHs and BNSs (Cutler & Harms, 2006), these experiments may only reach $\sim
1$deg2 in the angular resolution for BBHs (Zhao & Wen, 2018). This results
into a suppression of $C_{\rm GW-g}(\ell)$ significant at $\ell\ga 100$,

$\displaystyle C_{\rm GW-g}(\ell)$ $\displaystyle\rightarrow$ $\displaystyle
C_{\rm GW-g}(\ell)S(\ell)\ ,$ $\displaystyle E_{z}(\ell)$
$\displaystyle\rightarrow$ $\displaystyle E_{z}(\ell)S(\ell)\ ,$
$\displaystyle\frac{\partial E_{z}}{\partial m}$ $\displaystyle\rightarrow$
$\displaystyle\frac{\partial E_{z}}{\partial m}S(\ell)\ .$ (11)

Here $S(\ell)$ is the suppression caused by uncertainties in the angular
localization. We take a Gaussian form
$S(\ell)=\exp(-\ell^{2}\sigma^{2}_{\theta}/2)$, whose effective area is
$2\pi\sigma_{\theta}^{2}$. So the typical $\sigma_{\theta}\sim
1^{\circ}=\pi/180$ for ET and CT.

These experiments have a typical distance error $\sigma_{D}/D\la 0.1$ (Zhao &
Wen, 2018). Therefore luminosity bin width smaller than this is unnecessary.
Furthermore, the bin width needs to be sufficiently large to include enough GW
events. We then choose the luminosity bins $D_{L}H_{0}/c\in(0.51,1.19)$,
$(1.19,1.98)$, $(1.98,2.83)$. They correspond to $z\in(0.4,0.8)$, $(0.8,1.2)$,
$(1.2,1.6)$, if the distance measurement is exact. The analysis can be carried
out for other choices of luminosity-distance bins. The investigated ones, with
bin width much larger than $\sigma_{D}$, are chosen mainly for the purpose of
sufficient GW events. The number of GW dark sirens in each bin is $\sim
10^{4}$, with too large uncertainty to predict reliably. Therefore we adopt
$N_{\rm GW}=10^{4}$ as the fiducial value. Associated errors can be
conveniently rescaled by $N^{-1/2}_{\rm GW}$ to other values. The bottom line
is that we will work under the situation $N_{\rm GW}\sim 10^{4}\ll N_{\rm
g}\sim 10^{7-8}$.

Fig. 1 shows the dependence of $E_{z}$ on the fitting parameter $m^{\rm fit}$,
for the luminosity bin $D_{L}H_{0}/c\in(1.19,1.98)$. We show the cases of two
fiducial $m=0.0,1.0$ respectively. It is clear that $E_{z}$ reaches the
maximum when $m^{\rm fit}=m$. This verifies the theoretical prediction that
the $E_{z}$ estimator is unbiased. Furthermore, for the same $m^{\rm fit}$,
$E_{z}$ of different fiducial $m=0,1$ shows significant difference.

Figure 2: $\sqrt{\Delta\chi^{2}}$ as a function of $m^{\rm fit}$, which
determines the posterior probability of $m$ ($p(m^{\rm
fit})\propto\exp(-\Delta\chi^{2}/2)$). The fiducial $m=1$. The luminosity
distance range is $(1.19,1.98)c/H_{0}$. We adopt $N_{\rm GW}=10^{4}$, $b_{\rm
GW}=1$, $\sigma_{\ln D}=0.1$. From top to bottom,
$\sigma_{\theta}=0.5^{\circ},1^{\circ},2^{\circ}$ are adopted respectively.
The $1\sigma$ error in $m$ ranges from $0.07$ to $0.24$.

The constraining power of $E_{z}$ on $m$ is then quantified by the measurement
errors in $E_{z}$ and $\partial E_{z}/\partial m$. The measurement error in
$E_{z}$ arises both from that in $C_{\rm GW-g}$ and $C_{\rm g}$. Given the
availability of powerful galaxy surveys listed above, shot noise in $C_{g}$ at
scales of interest ($\ell\sim 100$) is negligible. Namely $C_{\rm g,N}=4\pi
f_{\rm sky}/N_{\rm g}\ll C_{\rm g}$. But that in the GW dark siren
distribution dominates over the cosmic variance. Namely the shot noise power
spectrum $C_{\rm GW,N}=4\pi f_{\rm sky}/N_{\rm GW}\gg C_{\rm GW}$. Therefore
the statistical error $\sigma_{E_{z}}$ in $E_{z}$ of a multipole bin centered
at $\ell$ and with bin width $\Delta\ell$ is dominated by error in $C_{\rm
GW-g}$. We then have

$\displaystyle\left(\frac{E_{z}}{\sigma_{E_{z}}}\right)^{2}_{\ell}$
$\displaystyle\simeq$ $\displaystyle\frac{2\ell\Delta\ell f_{\rm
sky}}{1+(C_{\rm GW}+C_{\rm GW,N})(C_{\rm g}+C_{\rm g,N})/C_{\rm GW-g}^{2}}\ .$

It can be further simplified,

$\displaystyle\left(\frac{E_{z}}{\sigma_{E_{z}}}\right)^{2}_{\ell}$
$\displaystyle\simeq$ $\displaystyle\frac{2\ell\Delta\ell f_{\rm
sky}}{1+C_{\rm GW,N}C_{\rm g}/C_{\rm GW-g}^{2}}$ (12) $\displaystyle\simeq$
$\displaystyle\frac{2\ell\Delta\ell f_{\rm sky}}{C_{\rm GW,N}C_{\rm g}/C_{\rm
GW-g}^{2}}=\frac{\ell\Delta\ell}{2\pi}N_{\rm GW}E_{z}^{2}$
$\displaystyle\simeq$ $\displaystyle
16\left(\frac{\ell\Delta\ell}{10^{4}}\right)\left(\frac{N_{\rm
GW}}{10^{4}}\right)\left(\frac{E_{z}}{10^{-3}}\right)^{2}\ .$

Notice that $N_{\rm GW}$ is the total number of GW events in the sky area of
the corresponding galaxy survey, which only covers a fraction $f_{\rm sky}$ of
the whole sky.

The statistical error in $\partial E_{z}/\partial m$ is $\sigma_{\partial
E_{z}/\partial m}=C^{1/2}$. The full expression of $C$ is given in Zhang
(2018b) . Adopting the same approximation in $E_{z}$, we have

$\displaystyle C$ $\displaystyle\simeq$ $\displaystyle\frac{1}{2\ell\Delta\ell
f_{\rm sky}}\frac{4\pi f_{\rm sky}}{N_{\rm GW}}$ (13) $\displaystyle\times$
$\displaystyle\frac{\int
W^{2}_{g,m}P_{g}\bar{n}_{g}^{2}\chi^{-2}(dz/d\chi)dz}{\int
W^{2}_{g}P_{g}\bar{n}_{g}^{2}\chi^{-2}(dz/d\chi)dz}\ $ $\displaystyle=$
$\displaystyle\frac{2\pi}{\ell\Delta\ell N_{\rm GW}}\frac{\int
W^{2}_{g,m}P_{g}\bar{n}_{g}^{2}\chi^{-2}(dz/d\chi)dz}{\int
W^{2}_{g}P_{g}\bar{n}_{g}^{2}\chi^{-2}(dz/d\chi)dz}\ .$

Here $W_{g,m}(z,m)\equiv\partial W_{g}(z,m)/\partial m$.

Fig. 2 shows the resulting $\Delta\chi^{2}$ as a function of the fitting
parameter $m^{\rm fit}$, for the fiducial case of $m=1.0$ and the luminosity
bin $D_{L}H_{0}/c\in(1.19,1.98)$. We adopt $\sigma_{D}/D=0.1$ and show its
dependence on the angular resolution $\sigma_{\theta}$. $\Delta\chi^{2}=1$
determines the $1\sigma$ error $\sigma_{m}$ in the $m$ constraint. The
resulting $\sigma_{m}$ is shown in Table 1. We show its dependence on the
luminosity bin, angular resolution and $\sigma_{D}/D$.

In general, $\sigma_{m}\sim 0.1$ is achievable. This will robustly distinguish
the scenario of BBH origin with $m=1$ proposed in Chen et al. (2019) from the
commonly assumed $m=0$ case. To further check these scenarios, we can split
BBHs of the same luminosity bin into two samples of mass above and below
$20M_{\odot}$, and apply the $E_{z}$ estimator separately. If the scenario
proposed by Chen et al. (2019) is valid, we will find that $m\sim 1$ and
$\sigma_{m}\sim 0.1$ for the sample of $M>20M_{\odot}$, while $m\sim 0$ and
$\sigma_{m}\sim 0.1$ for the sample of $M<20M_{\odot}$. If we apply the
$E_{z}$ estimator to the whole sample, we will find that $m$ significantly
deviates from zero. Furthermore, the minimum $\Delta\chi^{2}$ will be
significantly larger than unity, meaning bad fit of the data with a single $m$
and multiple origins of BBHs.

Table 1: $\sigma_{m}$, the constraining power in $m$, as a function of observed luminosity distance range, angular location accuracy $\sigma_{\theta}$, and distance measurement error $\sigma_{\ln D}$. The fiducial multiplicative error parameter $m=1$. The listed luminosity distance ranges correspond to redshift ranges $(0.4,0.8)$, $(0.8,1.2)$ and $(1.2,1.6)$ respectively, if we interpret them assuming $m=0$. $\sigma_{\theta}$ is the angular resolution, and third generation GW experiments have $\sigma_{\theta}\sim 1^{\circ}$. The quoted values of $\sigma_{m}$ outside of the parentheses have $\sigma_{\ln D}=0.1$ and those inside have $\sigma_{\ln D}=0.2$. The number of dark sirens within the observed luminosity-distance range is adopted as $N_{\rm GW}=10^{4}$.  $D_{L}^{\rm obs}H_{0}/c$ | $\sigma_{\theta}=0.5^{\circ}$ | $\sigma_{\theta}=1^{\circ}$ | $\sigma_{\theta}=2^{\circ}$
---|---|---|---
$(0.51,1.19)$ | $\sigma_{m}=0.07(0.11)$ | $0.11(0.16)$ | $0.19(0.26)$
$(1.19,1.98)$ | $\sigma_{m}=0.07(0.12)$ | $0.12(0.20)$ | $0.24(0.40)$
$(1.98,2.83)$ | $\sigma_{m}=0.08(0.15)$ | $0.13(0.27)$ | $0.35(--)$

There are many factors affecting the constraining power on $m$. Since
$C\propto N_{\rm GW}^{-1}$ and $dE_{z}/dm\propto b_{\rm GW}$, we have

$\displaystyle\sigma_{m}\propto N_{\rm GW}^{-1/2}b^{-1}_{\rm GW}\ .$ (14)

Here $b_{\rm GW}$ is the bias of GW dark sirens. We adopt $b_{\rm GW}=1$ and
$N_{\rm GW}=10^{4}$ as the fiducial values. We remind that, although
$\sigma_{m}$ depends on $b_{\rm GW}$, the bestfit $m$ does not (Zhang, 2018b).
Furthermore, $\sigma_{m}$ does not depend on the galaxy bias $b_{g}$, under
the limit of negligible shot noise in galaxy distribution. The major limiting
factor, from the observational side, is $N_{\rm GW}$. Therefore the
constraining power on $m$ will also significantly depend on the S/N threshold
of GW events, which changes $N_{\rm GW}$ significantly.

The dependences on other factors are more complicated. These include the
angular localization uncertainty $\sigma_{\theta}$, the distance measurement
error $\sigma_{\ln D}$, the range of luminosity distance and the true value of
$m$. One significant dependence is on $\sigma_{\theta}$ (Fig. 2 & Table 1).
$\sigma_{\theta}$ sets the maximum $\ell_{\rm max}\sim 2/\sigma_{\theta}\simeq
115/(\sigma_{\theta}/1^{\circ})$, available to the given GW experiment. This
factor alone would leads to $\sigma_{m}\propto\sigma_{\theta}$. But since shot
noise also increases with increasing $\ell$, the dependence on
$\sigma_{\theta}$ is slightly weaker. The dependence of $\sigma_{m}$ on
$\sigma_{D}/D$ is also significant, but for a different reason. Larger
$\sigma_{D}/D$ means wider (true) redshift distribution in GW dark sirens,
making the effort of constraining $m$ more difficult. $\sigma_{m}$ increases
slowly with increasing redshift/distance (Table 1), for fixed
$\sigma_{\theta}$, $\sigma_{\ln D}$, and $b_{\rm GW}$. The reason is that same
$\ell$ means smaller $k$ at higher redshift/distance, and therefore weaker
spatial correlation signal.

##  4 Discussion and conclusion

We have investigated the possibility of calibrating multiplicative errors in
the dark standard siren distance determination. Such errors cause systematic
deviation of the measured luminosity distance from the true cosmological
luminosity distance. Through cross correlation with pre-existing galaxy
redshift surveys, such systematic error can be constrained and the constraint
can be made model-independent and unbiased by the proposed $E_{z}$ estimator.
The constraining power is mainly limited by the angular resolution and
distance determination uncertainty. For ET/CE-like GW experiments,
$\sigma_{m}\sim 0.1$ is achievable. This constraint will provide unambiguous
distinction of some BBH origin mechanisms (e.g. Chen et al. (2019)). But we
caution that the discriminating power relies on the capability to
appropriately select GW samples. If the selected GW sample is composed of GW
events of mixed $m$, we may only constrain the mean $m$, but lose the power to
idendify those sub-samples with significant $m\neq 0$.

Multiplicate error is not the only form of systematic errors in the distance
determination. Additive error with non-zero average
($\langle\epsilon\rangle\neq 0$) shifts the mean value of $D_{L}^{\rm obs}$,
analogous to $z_{\rm bias}$ in photo-z errors. Weak lensing magnification in
the luminosity-distance measurement and magnification bias in the number
distribution of GW events and galaxies, can also bias the inference of GW
redshifts. Given the weak signal, the impact of weak lensing on $m$
determination is negligible for $\sigma_{m}\sim 0.1$. But it may become non-
nelgigible when $\sigma_{m}$ approaches $1\%$ or below. Other more complicated
form of systematic errors may also exist. Nevertheless, the same dark siren
LSS based calibration technique and the specific $E_{z}$ estimator are
applicable.

The expected $\sigma_{m}\sim 0.1$ is for third generation GW experiments and
the main limiting factors are the angular resolution and the distance
determination statistical errors. Although it is already useful for certain
applications, it does not meet the requirement to diagnose calibration error
in GW strain measurements, or modification of GW $D_{L}$ and EM $D_{L}$ by
modified gravity models still surviving. These applications require
$\sigma_{m}\sim 0.01$. Such constraining accuracy can be achieved by more
advanced surveys such as BBO. Its arcminute resolution means that multipole
modes with $\ell\ga 10^{2-3}$ are accessible. Its percent level distance
determination accuracy makes a finer luminosity bin size useful. Its
sensitivity to BNS events within the horizon will increase the dark siren
number density by an order of magnitude, making $\sim 10^{5}$ GW events per
luminosity-distance bin feasible. With these threefold improvements,
$\sigma_{m}=0.01$ is promising and new applications are possible.

## Data availability

No new data were generated or analyzed in support of this research.

Acknowledgement.— This work was supported by the National Science Foundation
of China 11621303 & 11653003\. We thank Xian Chen for valuable information on
the BBH origin mechanisms.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Physical Review Letters, 116, 061102
  * Abbott et al. (2017a) Abbott B. P., et al., 2017a, Classical and Quantum Gravity, 34, 044001
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, Physical Review Letters, 118, 221101
  * Abbott et al. (2017c) Abbott B. P., et al., 2017c, Physical Review Letters, 119, 141101
  * Abbott et al. (2017d) Abbott B. P., et al., 2017d, Physical Review Letters, 119, 161101
  * Abbott et al. (2017e) Abbott B. P., et al., 2017e, Nature, 551, 85
  * Abbott et al. (2019a) Abbott B. P., et al., 2019a, Physical Review X, 9, 031040
  * Abbott et al. (2019b) Abbott B. P., et al., 2019b, ApJ, 882, L24
  * Abbott et al. (2020a) Abbott R., et al., 2020a, arXiv e-prints, p. arXiv:2010.14527
  * Abbott et al. (2020b) Abbott B. P., Abbott R., Abbott T. D., Abraham S., et al. 2020b, ApJ, 892, L3
  * Abbott et al. (2020c) Abbott R., Abbott T. D., Abraham S., Acernese F., et al. LIGO Scientific Collaboration Virgo Collaboration 2020c, ApJ, 896, L44
  * Abdalla et al. (2015) Abdalla F. B., et al., 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), p. 17
  * Ajith & Bose (2009) Ajith P., Bose S., 2009, Phys. Rev. D, 79, 084032
  * Arun et al. (2007) Arun K. G., Iyer B. R., Sathyaprakash B. S., Sinha S., van den Broeck C., 2007, Phys. Rev. D, 76, 104016
  * Belgacem et al. (2018) Belgacem E., Dirian Y., Foffa S., Maggiore M., 2018, Phys. Rev. D, 97, 104066
  * Chen & Shen (2019) Chen X., Shen Z.-F., 2019, arXiv e-prints, p. arXiv:1906.11055
  * Chen et al. (2019) Chen X., Li S., Cao Z., 2019, MNRAS, 485, L141
  * Cutler & Harms (2006) Cutler C., Harms J., 2006, Phys. Rev. D, 73, 042001
  * Cutler & Holz (2009) Cutler C., Holz D. E., 2009, Phys. Rev. D, 80, 104009
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, preprint,  (arXiv:1611.00036)
  * Deffayet & Menou (2007) Deffayet C., Menou K., 2007, ApJ, 668, L143
  * Dodelson et al. (2016) Dodelson S., Heitmann K., Hirata C., Honscheid K., Roodman A., Seljak U., Slosar A., Trodden M., 2016, preprint,  (arXiv:1604.07626)
  * Dore (2018) Dore O., 2018, in 42nd COSPAR Scientific Assembly. pp E1.16–15–18
  * Doré et al. (2014) Doré O., et al., 2014, arXiv e-prints, p. arXiv:1412.4872
  * Dwyer et al. (2015) Dwyer S., Sigg D., Ballmer S. W., Barsotti L., Mavalvala N., Evans M., 2015, Phys. Rev. D, 91, 082001
  * Ezquiaga & Zumalacárregui (2018) Ezquiaga J. M., Zumalacárregui M., 2018, Frontiers in Astronomy and Space Sciences, 5, 44
  * Ferraro & Wilson (2019) Ferraro S., Wilson M. J., 2019, BAAS, 51, 72
  * Gong et al. (2019) Gong Y., et al., 2019, ApJ, 883, 203
  * Heymans et al. (2006) Heymans C., et al., 2006, MNRAS, 368, 1323
  * Holz & Hughes (2005) Holz D. E., Hughes S. A., 2005, ApJ, 629, 15
  * Laureijs et al. (2011) Laureijs R., et al., 2011, preprint,  (arXiv:1110.3193)
  * Libanore et al. (2020) Libanore S., et al., 2020, arXiv e-prints, p. arXiv:2007.06905
  * Massey et al. (2007) Massey R., et al., 2007, MNRAS, 376, 13
  * Mukherjee & Wandelt (2018) Mukherjee S., Wandelt B. D., 2018, preprint,  (arXiv:1808.06615)
  * Nair et al. (2018) Nair R., Bose S., Saini T. D., 2018, Phys. Rev. D, 98, 023502
  * Namikawa (2020) Namikawa T., 2020, arXiv e-prints, p. arXiv:2007.04359
  * Namikawa et al. (2016) Namikawa T., Nishizawa A., Taruya A., 2016, Physical Review Letters, 116, 121302
  * Newman (2008) Newman J. A., 2008, ApJ, 684, 88
  * Oguri (2016) Oguri M., 2016, Phys. Rev. D, 93, 083511
  * Reitze et al. (2019) Reitze D., LIGO Laboratory: California Institute of Technology LIGO Laboratory: Massachusetts Institute of Technology LIGO Hanford Observatory LIGO Livingston Observatory 2019, BAAS, 51, 141
  * Roy et al. (2019) Roy S., Sengupta A. S., Arun K. G., 2019, arXiv e-prints, p. arXiv:1910.04565
  * Santos et al. (2015) Santos M., et al., 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), p. 19
  * Schneider et al. (2006) Schneider M., Knox L., Zhan H., Connolly A., 2006, ApJ, 651, 14
  * Schutz (1986) Schutz B. F., 1986, Nature, 323, 310
  * Shao et al. (2011) Shao J., Zhang P., Lin W., Jing Y., 2011, ApJ, 730, 127
  * Soares-Santos et al. (2019) Soares-Santos M., et al., 2019, ApJ, 876, L7
  * Spergel et al. (2015) Spergel D., et al., 2015, preprint,  (arXiv:1503.03757)
  * Sun et al. (2020) Sun L., et al., 2020, arXiv e-prints, p. arXiv:2005.02531
  * Takada et al. (2014) Takada M., et al., 2014, PASJ, 66, R1
  * The LIGO Scientific Collaboration et al. (2019) The LIGO Scientific Collaboration the Virgo Collaboration Abbott B. P., Abbott R., et al. 2019, arXiv e-prints, p. arXiv:1908.06060
  * The LIGO Scientific Collaboration et al. (2020a) The LIGO Scientific Collaboration the Virgo Collaboration Abbott R., et al. 2020a, arXiv e-prints, p. arXiv:2004.08342
  * The LIGO Scientific Collaboration et al. (2020b) The LIGO Scientific Collaboration et al., 2020b, arXiv e-prints, p. arXiv:2010.14533
  * Van Den Broeck & Sengupta (2007) Van Den Broeck C., Sengupta A. S., 2007, Classical and Quantum Gravity, 24, 1089
  * Wang et al. (2019) Wang Y., et al., 2019, BAAS, 51, 508
  * Zhan (2018) Zhan H., 2018, in 42nd COSPAR Scientific Assembly. pp E1.16–4–18
  * Zhang (2018a) Zhang P., 2018a, preprint,  (arXiv:1810.11915)
  * Zhang (2018b) Zhang P., 2018b, arXiv e-prints, p. arXiv:1811.07136
  * Zhang & Pen (2001) Zhang P., Pen U.-L., 2001, ApJ, 549, 18
  * Zhang et al. (2010) Zhang P., Pen U.-L., Bernstein G., 2010, MNRAS, 405, 359
  * Zhang et al. (2017) Zhang L., Yu Y., Zhang P., 2017, ApJ, 848, 44
  * Zhao & Wen (2018) Zhao W., Wen L., 2018, Phys. Rev. D, 97, 064031