# Clustering of Gravitational Wave and Supernovae events:
a multitracer analysis
in Luminosity Distance Space

S. Libanore,  M.C. Artale  D. Karagiannis  M. Liguori
N. Bartolo  Y. Bouffanais  M. Mapelli  S. Matarrese

###### Abstract

We study the clustering of Gravitational Wave (GW) merger events and
Supernovae IA (SN), as cosmic tracers in Luminosity Distance Space. We modify
the publicly available CAMB code to numerically evaluate auto- and cross-
power spectra for the different sources, including Luminosity Distance Space
distortion effects generated by peculiar velocities and lensing convergence.
We perform a multitracer Fisher analysis to forecast expected constraints on
cosmological and GW bias coefficients, using outputs from hydrodinamical
N-body simulations to determine the bias fiducial model and considering future
observations from the Vera Rubin Observatory and Einstein Telescope (ET), both
single and in a 3 detector network configuration. We find that adding SN to
the GW merger dataset considerably improves the forecast, mostly by breaking
significant parameter degeneracies, with final constraints comparable to those
obtainable from a Euclid-like survey. GW merger bias is forecasted to be
detectable with good significance even in the single ET case.

##  1 Introduction

The last six years have witnessed the first three observing runs of Advanced
LIGO and Virgo, leading to the detection of more than 50 gravitational wave
(GW) events [2, 1, 3] from compact binary mergers. This sample will
conspicuously grow in the next two observing runs [15]. Furthermore, third
generation GW detectors, such as Einstein Telescope111http://www.et-gw.eu
(ET), will lead to a dramatic increase of the detection rate, observing
$\sim{10^{5}}$ GW events per year with accurate luminosity distance
determination [35]. The sky localization precision will be in a range between
tenths to a hundred deg2 for a single detector and it will significantly
increase using networks of third generation detectors (both ET-like or Cosmic
Explorer222https://cosmicexplorer.org/-like). This scenario will allow us to
perform statistical studies of the clustering properties of GW mergers. Since
these are tracers of the Large Scale Structure (LSS) of the Universe, this
will enable us to constrain both cosmological and merger bias parameters, both
cross-correlating with galaxy maps [42, 22, 12] and by using GW surveys alone
[40, 53, 33, 48].

For analogous studies of 3D galaxy clustering, we use redshifts as distance
indicators; however, for GW events we directly measure luminosity distances,
not redshifts. If we want to identify the redshifts of GW sources, we need to
rely on either electromagnetic counterparts or statistical methods, based on
cross-correlations with theoretical models or other LSS tracers. These
approaches present some drawbacks: both rely on extra-observations or external
datasets and are difficult or impossible to apply for high redshift sources;
moreover, the statistical analysis can depend on specific assumptions. The
natural alternative to redshift identification, which we consider in this
work, is that of directly using luminosity distance as our distance indicator
in the 3D clustering analysis, as originally proposed in [40, 53]. In this
way, no external datasets or assumptions are needed and all mergers, including
very high redshift ones, can be included in the analysis.

In a previous work [33], we produced cosmological and merger bias parameter
forecasts, considering third generation interferometers and more futuristic
scenarios. We underlined that the main advantages of this kind of survey are
the large volumes probed and the capability of constraining the merger bias at
high statistical significance, which can provide interesting information about
the physical nature and properties of mergers themselves. In this paper we
significantly expand our previous study, in three main directions.

Firstly, we include new tracers. In our previous analysis, we considered only
GW events produced by Double Neutron Star (DNS) and Double Black Hole (DBH)
mergers; the same procedure can be used to analyse Black Hole-Neutron Star
mergers (BHNS) as well. Moreover, another type of observable for which a
Luminosity Distance Space (LDS) clustering analysis is natural exists, namely
Supernovae IA (SN). Here, we therefore include SN in the analysis.

Secondly, we improve our treatment of Luminosity Distance Space Distortions
(LDSD). Differently from redshift space distortions (RSD), LDSD display a
leading order contribution coming from the Line-of-Sight (LoS) derivative of
the lensing convergence. In [33], following an heuristic argument presented in
[53], we neglected this contribution in our power spectrum calculation.
Lensing of GW is an active field of study, see e.g. [19]) for a treatment of
lensing effects in the GW signal and e.g. [37] for a weak lensing analysis. As
for LDSD, the lensing term was explicitly derived and discussed in [39], where
the resulting merger power spectrum was evaluated using a Limber
approximation. Following that calculation, here we implement a modified
version of the CAMB code, which evaluates this term in full sky, so that we
can propagate and assess its impact on the final parameter forecasts. We find
that the lensing LDSD contribution only affects the large scales and has a
very small impact on the final parameter forecasts. These findings are in line
with both the power spectrum computation333Note that for this analysis we need
to consider only the monopole term. As shown in [39], the lensing contribution
has a larger effect and a distinctive signature on higher order multipoles. of
[39] and the general argument raised in [53].

Finally, we go beyond the single tracer approach of our previous work and
develop a multitracer analysis. We show that SN play a significant role in
breaking degeneracies arising in a GW-only analysis, thus leading to important
improvements in the expected constraints. As in our previous work, we obtain
our fiducial values for merger bias parameters by building an Halo Occupation
Distribution (HOD), based on full hydrodynamical simulations [5, 6]; the
fiducial model for SN is instead based on the literature [37]. We also study
how the inclusion of a prior on one of such parameters affects the results.
Our full multitracer analysis is performed assuming that GW observations are
made using Einstein Telescope either in single (ET) or in three detector
network configuration (ET$\times$3). For Supernovae IA instead, the Vera Rubin
Observatory444https://www.lsst.org (VRO), previously known as Large Synoptic
Survey Telescope (LSST), is considered.

This paper is structured as follows: in sec. 2 we describe the number
distributions and bias models of the tracers (mergers and SN) used in the
analysis. In sec. 3 we present the observational effects that must be taken
into account when mapping the tracers in LDS, namely the Luminosity Distance
Space Distortions and the lensing (other effects are neglected in this work
being subdominant and still not well described at the theoretical level). Sec.
4 presents the formalism we used to perform the single and multitracer
analysis, while sec. 5 describes the constraints we obtained both on
cosmological and bias parameters.

##  2 Tracer distributions and bias

Working in Luminosity Distance Space can present significant advantages when
studying the spatial distribution of sources such as GW produced by compact
binary mergers and SN. Since we assume General Relativity in our work, the
luminosity distance $D_{L}$ represents the same kind of measurement for both
merger and SN observations. In some scenarios involving dark energy or
modified gravity, however, the luminosity distance extracted from the
gravitational wave signal $D_{L}^{GW}$ turns out to be different from the
electromagnetic one $D_{L}^{EM}$, both in the background (see e.g. [8, 7]) and
when considering fluctuations (compare with [24, 21]), making the radial
coordinates $D_{L}^{GW}$ and $D_{L}^{EM}$ different. In these scenarios, our
analysis would be self-consistent except for the multitracer case with
combined mergers and SN (sec. 5.2). As shown by [25], the combination of the
two tracers would still be possible using other techniques and it should allow
constraining the dark energy or modified gravity models themselves.

To study the constraining power of future LDS surveys, first and foremost the
distribution and clustering of the tracers have to be modelled. On the one
hand, binary coalescences and merger signals are very well described (see e.g.
[16]): from their amplitude and frequency variation it is possible to extract
a full set of intrinsic and observational parameters, including their
luminosity distance $D_{L}$ [47, 51]. On the other hand, the explosion
mechanism which produces the SN determines the light curve shape (see e.g.
[32]) and the absolute magnitude of the peak, allowing us to estimate $D_{L}$
[30] after a proper calibration procedure.

In sec. 2.1, we model the distribution and clustering properties of GW
referring to compact binary mergers by means of numerical simulations,
suitably processed to account for observational selection effects. In sec.
2.2, we instead model the SN distribution, by taking into account the event
rate and completeness of the observations performed.

Despite the analysis is performed in LDS, along the paper the redshift
notation is frequently used, since it is more familiar to the reader and more
similar to the one used when studying LSS. The consistency of the analysis is
preserved since the redshift is always computed from luminosity distance
assuming the Planck 2018 [18] Cosmology.

###  2.1 Gravitational Wave mergers

In our previous work [33], we explicitly derived the number distribution and
bias of DBH and DNS mergers; here, we also consider analogously BHNS systems.
In doing so, we rely on the simulations discussed in [5, 6] to get the
probability distribution of the mergers, as a function of their redshift $z$
and of the host galaxy stellar mass $M_{*}$ and star formation rate $SFR$. The
population-synthesis simulations and their coupling with the cosmological
simulation are the same as we described in [33]. We define the total number of
mergers in the simulations as:

$N_{m}^{SIM}(z)=\sum_{i}\sum_{j}\bigl{<}N^{SIM}_{m}(z)|M_{*}^{i},SFR^{j}\bigr{>}\
,$ (2.1)

where $m$ indicates DBH, BHNS or DNS, while $i$, $j$ indicate respectively the
$M_{*}$ and $SFR$ bins used in the simulation. Selection effects due to the ET
Signal-to-Noise Ratio (SNR) are taken into account inside the simulations; the
detection rate is then directly computed by converting the time intervals from
the source to the observer rest frame, through:

$N_{m}(z)\ =\ T^{OBS}_{m}\ \frac{N_{m}^{SIM}(z)}{T^{SIM}(z)}\ \frac{1}{1+z}\
.$ (2.2)

$T^{SIM}(z)$ represents the width of the redshift bins used for the simulation
snapshots $[z-\delta z,z+\delta z]$ expressed in Gyr, while $T^{OBS}_{m}$ is
the duration of the observational campaign, which we assume to be 3 years
long.

The observed number density per unit redshift and solid angle can be expressed
as:

$\frac{d^{2}N_{m}}{dzd\Omega}=N_{m}(z)\frac{c}{\ell\
H(z)}\biggl{(}\frac{D_{L}(z)}{\ell\ (1+z)^{2}}\biggr{)}^{2}\ ,$ (2.3)

where $\ell=25$Mpc is the overall size of the simulation box, $H(z)$ is the
Hubble factor and $c$ is the speed of light. Here, we adopt the same procedure
used in [33] and interpolate the source number density measured from
simulations with a skewed Gaussian; the resulting distributions are shown in
fig. 1.

The probability distribution
$\bigl{<}N^{SIM}_{m}(z)|M_{*}^{i},SFR^{j}\bigr{>}$ can be used also to compute
the bias of the mergers using an HOD based approach. We follow the two-step
procedure we discussed in detail in [33] and start by computing the bias of
host galaxies in each $M_{*}$ and $SFR$ bin as:

$b_{g}(z,M_{*},SFR)=\int_{M_{h}^{min,(M_{*},SFR)}}^{+\infty}dM_{h}\
n_{h}(z,M_{h})\
b_{h}(z,M_{h})\frac{\bigl{<}N_{g}(M_{*},SFR)|M_{h}\bigr{>}}{n_{g}(z,M_{*},SFR)}\
,$ (2.4)

where $n_{h}(z,M_{h})$, $b_{h}(z,M_{h})$ are the halo mass function and the
bias factor [46] respectively, while $M_{h}^{min,(M_{*},SFR)}$ represents the
minimum mass required for an halo to form a galaxy with a given stellar mass
$M_{*}$ and star formation rate $SFR$. Moreover,
$\bigl{<}N_{g}(M_{*},SFR)|M_{h}\bigr{>}$ is the galaxy Halo Occupation
Distribution (HOD) in the eagle simulation [44], defined as the number of
galaxies with stellar mass $M_{*}$ and star formation rate $SFR$, within a
dark matter halo of mass $M_{h}$. The galaxy mean number density,
$n_{g}(z,M_{*},SFR)$, is defined as:

$n_{g}(z,M_{*},SFR)=\int_{M_{h}^{min,(*,SFR)}}^{+\infty}dM_{h}\
n_{h}(z,M_{h})\ \bigl{<}N_{g}(M_{*},SFR)|M_{h}\bigr{>}\ .$ (2.5)

The merger bias can then be computed in the second step, through:

$b_{m}(z)=\int_{M_{*}^{min}}^{M_{*}^{max}}dM_{*}\
\int_{SFR^{min}}^{SFR^{max}}dSFR\ n_{g}(z,M_{*},SFR)\
b_{g}(z,M_{*},SFR)\frac{\bigl{<}N_{m}(z)|M_{*},SFR\bigr{>}}{n_{m}(z)}\ .$
(2.6)

In the previous equation, $n_{m}(z)$ represents the merger mean number
density, which is computed in analogy to eq. (2.5) through:

$n_{m}(z)=\int_{M_{*}^{min}}^{M_{*}^{max}}\int_{SFR^{min}}^{SFR^{max}}dM_{*}\
dSFR\ n_{g}(z,M_{*},SFR)\ \bigl{<}N_{m}^{SIM}(z)|M_{*},SFR\bigr{>}\ .$ (2.7)

Results obtained in [33] show that the bias for both DBH and DNS can be
modelled as:

$b_{m}(z)=A_{m}(z+B_{m})^{P_{m}}\ ,$ (2.8)

where the slope $P_{m}$ is well described by a linear behaviour, i.e.
$P_{m}=1$. $A_{m}$ and $B_{m}$ fiducial values are extracted from the
simulations, obtaining for the different kinds of mergers the following
values: $A_{DBH}=0.7$, $B_{DBH}=2.68$, $A_{DNS}=0.76$, $B_{DNS}=2.46$. Here,
we verify that eq. (2.8) well describes also the BHNS bias, once that the
fiducial values $A_{BHNS}=0.71$, $B_{BHNS}=2.63,P_{BHNS}=1$ are used. These
results are summarized in fig. 2.

###  2.2 Supernovae IA

To model the SN number distribution, we combine the SN rate provided in [17]
with the Status Quo completeness that the work by [23] describes for the Vera
Rubin Observatory. We assume a $T^{OBS}_{SN}=5$yr observational run. The rate
in the source rest frame is modelled as:

$r_{IA}=2.6\cdot 10^{-5}(1+z)^{2.5}\ h^{3}_{70}\ \text{SN
yr}^{-1}\text{Mpc}^{-3}\ ,$ (2.9)

where $h_{70}$ is the rescaling of the Hubble parameter determined by the
fiducial Cosmology. The conversion to the observer rest frame is performed
using the same factor as in eq. (2.2). The completeness of the survey, i.e.,
the fraction of events that are detected, is modelled as:

$\mathcal{C}=\gamma_{SN}\biggl{(}\frac{z}{\alpha_{SN}}\biggr{)}^{\beta_{SN}-1}e^{(-z/\alpha_{SN})^{\beta_{SN}}}\
,$ (2.10)

where $\alpha_{SN}$, $\beta_{SN}$, $\gamma_{SN}$ are chosen to reproduce the
completeness trend in [23]. A comoving box having size $\ell=1$Mpc is
considered to compute the observed number of SN events as:

$N_{SN}(z)=T^{OBS}_{SN}\mathcal{C}\ \ell^{3}\ r_{IA}\ .$ (2.11)

The distribution we obtain in this way is shown in fig. 1; the total number of
events observed is $\sim 10^{4.86}$, which is slightly more pessimistic than
what both [20] and [23] report.

To model the SN bias, we assume a constant value in the full observed redshift
range (i.e. $z\in[0,1]$). This follows the conservative prescription made in
[37], where SN are assumed to follow galaxy bias. The value we choose is
$b_{SN}=1.9$, which is the $b_{g}(z)$ value obtained by integrating eq. (2.4)
over $M_{*}$ and $SFR$ and averaging it between $z=0$ and $z=1$. In order to
use the same model as in eq. (2.8), the parameters $A_{SN}=1.9$, $B_{SN}=0$
and $P_{SN}=0$ are defined. Fig. 2 shows the bias we obtain.

Figure 1: $d^{2}N/dzd\Omega$ for the sources considered in the analysis: SN
observed by the Vera Rubin Observatory in the Status Quo configuration (green
dashed, see sec. 2.2), DNS (red, see sec. 2.1) and DBH observed by ET (blue,
see sec. 2.1). $z_{max}$ per each source depends on the detector horizon (see
tab. 2). BHNS observed by ET are not showed for clarity; their model is in
between DNS and DBH. Figure 2: Fiducial bias models considered in the
analysis. The color legend is the same as in fig. 1.

##  3 Observational effects in Luminosity Distance Space

When we study galaxy clustering using observational survey catalogs, we need
to include corrections to account for redshift space distortions and
relativistic effects, which alter the observed density field with respect to
the intrinsic one (see e.g. [27, 52, 14, 11]). When working in Luminosity
Distance Space, the situation is of course analogous. However, the mapping
between the intrinsic density field and the observed one is different, with
respect to the redshift space case. Large scale structures induce a
perturbation $\epsilon$ on luminosity distance, defined at leading order as
follows (cfr [28, 41]):

$\frac{\bar{D}_{L}+\delta D_{L}}{\bar{D}_{L}}=1+\epsilon=1+2v-\kappa\ ,$ (3.1)

where $\bar{D}_{L}$ is the background, unperturbed luminosity distance,
$\delta D_{L}$ is the induced perturbation, $v$ is the radial velocity of the
source along the LoS and $\kappa$ is the lensing convergence. For a given
comoving distance $\chi$ and LoS $\hat{n}$, the latter is defined as:

$\kappa=\int_{0}^{\chi}d\chi^{\prime}\frac{(\chi-\chi^{\prime})\chi^{\prime}}{\chi}\nabla^{2}_{\hat{n}}\
\psi[\chi^{\prime},\hat{n}\chi^{\prime}]\ ,$ (3.2)

where $\nabla^{2}_{\hat{n}}\ \psi[\chi^{\prime},\hat{n}\chi^{\prime}]$ are the
angular components of the covariant derivative of the Weyl potential on a
sphere having radius $\chi^{\prime}$.

In our previous work [33], we considered only the effects of LDSD coming from
peculiar velocities $v$, thus neglecting $\kappa$. This was motivated by the
fact that - as discussed for example in [53] \- peculiar velocities were
expected to be the dominant source of distortion. After our previous analysis
was completed, the lensing contribution was explicitly computed in [39],
showing that, while indeed subdominant, such contribution is not negligible on
large scales.

Other effects in LDS such as gravitational potential perturbations and
relativistic effects, are included neither in eq. (3.1) nor in what follows.
Such effects are still not well modelled in the literature but, despite being
subdominant [28], their analysis will be very important when considering full
sky surveys and large scales, which we do not include in this work. For all
these reasons, while beyond the scope of the present study, we plan on
exploring these issues in a future, dedicated analysis.

###  3.1 Lensing

Our analysis of lensing contributions to LDSD starts from the results of [39],
which slightly differ from e.g. those reported in [28]. In particular the
author of [39] obtains that the lensing contribution to LDSD includes both the
standard term $-2\kappa$, due to the lensing distortion of the angular
position, and an extra-term that depends on the LoS derivative of the lensing
convergence with respect to the comoving distance. Both these terms, together
with the peculiar velocity contribution, alter the observed density
fluctuations $\delta_{LDS}$ with respect to the intrinsic ones $\delta$.
Naming $\delta_{LDS}^{v}$ and $\delta_{LDS}^{\kappa}$ the peculiar velocity
and lensing distortions on $\delta$ (i.e. propagating the terms in eq. (3.1)
through the number conservation law as in the standard Kaiser approach [31]),
the observed density fluctuations can be described as:

$\displaystyle\delta_{LDS}$
$\displaystyle=\delta+\delta_{LDS}^{v}+\delta_{LDS}^{\kappa}$ (3.3)
$\displaystyle=\delta-2\alpha
v-\frac{f_{D_{L}}}{aH}\frac{dv}{d\chi}+\biggl{(}-2+\alpha+\frac{f_{D_{L}}}{2aH}\frac{d}{d\chi}\biggr{)}\kappa\
,$

where $a$ is the scale factor and $H$ represents the Hubble parameter
(redshift dependencies are omitted for clarity). The parameters $f_{D_{L}}$
and $\alpha$ are defined as:

$\displaystyle f_{D_{L}}=$ $\displaystyle\
\frac{2\bar{D}_{L}}{1+z}\biggl{(}\frac{\partial\bar{D}_{L}}{\partial
z}\biggr{)}^{-1}\ ,$ (3.4) $\displaystyle\alpha=$ $\displaystyle\
\biggl{[}3+\frac{d\ln n_{i}}{d\ln\chi}-\frac{f_{D_{L}}}{2}+\frac{\chi
f_{D_{L}}}{2aH}\frac{\partial(aH)}{\partial\eta}{\color[rgb]{0,.5,.5}\
}\biggr{]}\frac{1}{1+aH\chi}\ ,$ (3.5)

where $\eta$ is the conformal time and $n_{i}$ the comoving number density of
the source ($i=$ DBH, DNS, BHNS or SN). The latter is computed by dividing
$N_{m}$ from eq. (2.2) – for mergers – or $N_{SN}$ from eq. (2.11) – for SN –
by the comoving volume. The obtained distribution is interpolated with a
skewed gaussian and analytically differentiated in order to get $d\ln
n_{i}/d\ln\chi$. Eq. (3.3) is the same as the one originally derived in [39],
where $\gamma={f_{D_{L}}/2}$ is used. The term $f_{D_{L}}$ displays an
interesting behaviour, as discussed in [53, 33]: it approaches zero at low $z$
and increases at higher redshift, crossing the value $f_{D_{L}}=1$ at $z\simeq
1.7$. This means that at high redshifts ($z\gtrsim 2$), LDSD are larger than
redshift space distortions, whereas at lower redshift, LDSD turn out to be
very small. Since, as showed in sec. 2.2, the SN distribution is relevant only
for $z\lesssim 1$, this implies that SN observations are almost unaffected by
LDSD. Fig. 3 shows the results we obtained for $\alpha$, which slightly differ
from those in [39], because of the different comoving number densities used in
our work.

Figure 3: $\alpha$ computed using eq. (3.5). Each curve assumes a different
$n_{i}$ evolution in $\chi$. The distributions for SN (green dashed), DNS
(red) and DBH (blue) are computed as described in the text, while for non
evolving sources (black dotted) $d\ln n_{i}/d\ln\chi=0$. BHNS are not showed
for clarity; their $\alpha$ is intermediate between DBH and DNS. In the case
$z=0$, all the distributions consider $\alpha=3$: this directly follows from
the number conservation of sources between LDS and real space (see [39] for
the complete derivation).

In [39], it is shown that the lensing contribution to the overdensity (see eq.
(3.3)) can be written as:

$\delta_{LDS}^{\kappa}=\int_{0}^{\infty}d\chi^{\prime}w(\chi,\chi^{\prime})\nabla_{\hat{n}}^{2}\
\psi[\eta(\chi^{\prime}),\hat{n}\chi^{\prime}]\ ,$ (3.6)

where the lensing kernel $w(\chi,\chi^{\prime})$ is computed using eq. (3.2),
its derivative and eq. (3.3):

$w(\chi,\chi^{\prime})=\frac{\chi^{\prime}}{\chi}\biggl{[}(-2+\alpha)(\chi-\chi^{\prime})+\frac{f_{D_{L}}}{2aH}\frac{\chi^{\prime}}{\chi}\biggr{]}\
.$ (3.7)

This expression differs from that of the “standard” redshift space kernel,
which is
$w(\chi,\chi^{\prime})=(-2+5s_{*})(\chi-\chi^{\prime})\chi^{\prime}/\chi$
(cfr. [14]). We note in particular that in this formula the $f_{DL}$ and
$\alpha$ terms defined in eq. (3.4) and (3.5) introduce the dependence on the
distance from the observer, on the lensing convergence derivative along the
LoS and on the source mean number density $n_{i}$. Moreover, the magnification
bias, $s_{*}$, is not included. In the case of GW surveys, this was firstly
computed in [42] as:

$s_{*}=-\frac{d\log\bigl{[}d^{2}N\bigl{(}z,\sqrt{\bigl{<}\rho^{2}\bigr{>}}>\rho_{lim})/dzd\Omega\bigr{]}}{d\sqrt{\langle\rho^{2}\rangle}}\Biggl{|}_{\rho_{lim}}\
$ (3.8)

where $\rho$ is the signal-to-noise ratio (SNR),
$\sqrt{\langle\rho^{2}\rangle}$ represents its averaged value and
$\rho_{lim}=8$ is the detection threshold. This term is neglected in [39]: we
decided to keep this formulation at this level, but we plan to include a more
refined model for $s_{*}$ in a forthcoming paper, where the full analysis of
the simulated catalog and of LDSD will be performed. In the case of SN, the
magnification bias is commonly neglected in literature.

Figure 4: Overdensity (black), LDS distortions due to peculiar velocities
(purple) and lensing (green) contribution to the $C_{l}$ when Limber is not
considered. Dashed lines are used for negative contributions. The left plot
refers to low redshift $z\in[0.6,0.8]$, while the right one to high redshift
$z\in[2.6,2.8]$. All the curves refer to the DBH distribution, without
including the smoothing effect from eq. (4.7). Figure 5: Comparison between
peculiar velocities distortions in LDS (purple) and redshift space (pink) and
between lensing in LDS (green) and redshift space (light green). On the left,
the DBH distribution at low redshift $z\in[0.6,0.8]$ is considered, while on
the right the same at $z\in[2.6,2.8]$ is used, without including the smoothing
effect described in eq. (4.7). Dashed lines are used for negative
contributions. The lower plots show the residuals relative to the absolute
values of quantities in LDS.

In [39], the lensing correction is included in the computation of the two
point correlation function, assuming the Limber approximation [34]. In the
current analysis, we implement the new lensing kernel from eq. (3.7) inside
CAMB555https://github.com/cmbant/CAMB [14] – as we did in our previous work
[33] for LDSD from peculiar velocities. Since CAMB assumes the Limber
approximation only on very small scales (i.e. $l>100$ in our case), our
implementation properly accounts for full-sky effects in LDS. Fig. 4 shows the
resulting contribution of peculiar velocity and lensing to LDSD, with respect
to the intrinsic overdensity, in the case of low and high distances. These
results are consistent with the ones from [39]. The shift in the lensing term
when moving from redshift to LDS is bigger than the shift in the peculiar
velocity case, as fig. 5 shows; however, the lensing term still remains
subdominant in LDS.

##  4 Methodology

###  4.1 Multitracer analysis

The multitracer technique was firstly developed in [45] (see also [36, 26, 4])
to overcome the cosmic variance problem. This is relevant at large scales,
where the stochastic distribution of dark matter overdensities can be measured
in only a few realizations: since different tracers map the same underlying
distribution, by comparing them it is possible to measure the bias ratio
without modelling the dark matter field itself (cfr. [50]). Therefore,
combining and cross-correlating data from different surveys provide better
parameter estimates and forecasts.

Our forecast analysis expands and improves over our previous work [33] both by
including the lensing perturbation described in sec. 3.1 and by applying the
multitracer technique. The multitracer analysis in our case brings together
the GW survey performed by ET/ET$\times$3 and VRO observations of SN. As a
further case, we also check the results for the combination of different
merger events, e.g. DBH and DNS. Despite being performed with the same
detector (i.e. ET/ET$\times$3), GW from different mergers can be considered as
observed by independent surveys since their signal can be easily
distinguished. For this reason, we refer as well to combinations of different
kinds of mergers with the name of "multitracer analysis".

In order to test the stability of our results, we compute the Fisher matrix in
two different, yet equivalent, ways. First of all, we consider the Fisher
matrix for the overdensity field, which is Gaussian distributed with zero
average:

$F_{\alpha\beta}=\sum_{l}\frac{2l+1}{2}f_{sky}\text{Tr}\bigl{[}\partial_{\alpha}\mathbf{C}_{l}\
\Gamma_{l}^{-1}\partial_{\beta}\mathbf{C}_{l}\ \Gamma_{l}^{-1}\bigr{]}\ ,$
(4.1)

where $f_{sky}$ is the observed sky fraction and
$\partial_{\alpha,\beta}\mathbf{C}_{l}$ is the matrix of the derivatives with
respect to the parameters of interest (see sec. 4.3 for details) of the auto-
and cross- power spectra. $\Gamma_{l}=\mathbf{C}_{l}+\mathbf{N}_{l}$ is
instead the field covariance matrix, including noise contributions. We only
assume shot noise contributions, uncorrelated between different bins and
tracers, i.e.
$N_{ij}^{s_{a},s_{b}}=\delta^{K}_{s_{a},s_{b}}\delta^{K}_{ij}[N_{s_{a}}(z_{i})]^{-1}$
(being $s_{a,b}=$ SN or DNS, BHNS, DBH) . We construct $\mathbf{C}_{l}$ (and
analogously $\Gamma_{l}$) as a block matrix combining all the tracers and
distance bins.666For example, in the case of $2$ tracers and $2$ bins
$\mathbf{C}_{l}$ is defined as:
$\mathbf{C}_{l}=\left(\begin{array}[]{cc|cc}C_{00}^{s_{1},s_{1}}&C_{00}^{s_{1},s_{2}}&C_{01}^{s_{1},s_{1}}&C_{01}^{s_{1},s_{2}}\\\
C_{00}^{s_{2},s_{1}}&C_{00}^{s_{2},s_{2}}&C_{01}^{s_{2},s_{1}}&C_{01}^{s_{2},s_{2}}\\\
\hline\cr
C_{10}^{s_{1},s_{1}}&C_{10}^{s_{1},s_{2}}&C_{11}^{s_{1},s_{1}}&C_{11}^{s_{1},s_{2}}\\\
C_{10}^{s_{2},s_{1}}&C_{10}^{s_{2},s_{2}}&C_{11}^{s_{2},s_{1}}&C_{11}^{s_{2},s_{2}}\\\
\end{array}\right)_{l}$ (4.2) We then take an "estimator perspective", i.e.,
we consider as our actual observables the estimates $\hat{C}$ of the various
auto- and cross- power spectra, rather than the overdensity field. Power
spectrum estimates are then characterized by a Wishart sampling distribution.
Under a Gaussian approximation, one can write the corresponding Fisher matrix
as [13, 29]:

$F_{\alpha\beta}=\sum_{l}F_{\alpha\beta,l}=\sum_{l}(2l+1)f_{sky}\
\bigl{[}(\partial_{\alpha}\mathsf{C}_{l})^{T}\ \Gamma_{l}^{{}^{\prime}-1}\
\partial_{\beta}\mathsf{C}_{l}\bigl{]}\ ,$ (4.3)

$\mathsf{C}_{l}=\bigl{[}\hat{C}_{00}^{s_{1},s_{2}}\
\hat{C}_{00}^{s_{1},s_{2}}\ ...\bigr{]}_{l}$ being the vector containing all
the auto- and cross- power spectra and $\Gamma_{l}^{{}^{\prime}}$ being the
power spectrum covariance matrix, defined as:

$\displaystyle\Gamma_{l}^{\prime}$
$\displaystyle=\text{cov}\bigl{[}\hat{C}_{ij}^{s_{a},s_{b}}\hat{C}_{km}^{s_{c},s_{d}}\bigr{]}_{l}$
(4.4)
$\displaystyle=\bigl{[}(C_{ij}^{s_{a},s_{b}}+N_{ij}^{s_{a},s_{b}})(C_{km}^{s_{c},s_{d}}+N_{km}^{s_{c},s_{d}})+(C_{ik}^{s_{a},s_{c}}+N_{ik}^{s_{a},s_{c}})(C_{jm}^{s_{c},s_{d}}+N_{jm}^{s_{c},s_{d}})\bigr{]}_{l}\
.$

The number of spectra involved in the computation scales as
$n_{i}^{2}n_{b}^{2}$, where $n_{i}$, $n_{b}$ are the number of sources and the
number of bins used in the analysis, respectively.777Considering for example
$1$ tracer in $2$ bins, the terms $F_{\alpha\beta,l}$ in eq. (4.3) would be
rewritten as:  $\displaystyle F_{\alpha\beta,l}$
$\displaystyle=\sum_{l}(2l+1)f_{sky}\left[\begin{array}[]{cccc}\hat{C}_{00}^{s_{1},s_{1}}&\hat{C}_{01}^{s_{1},s_{1}}&\hat{C}_{10}^{s_{1},s_{1}}&\hat{C}_{11}^{s_{1},s_{1}}\end{array}\right]\times$
(4.5) $\displaystyle\times$
$\displaystyle\left[\begin{array}[]{cccc}\text{cov}[\hat{C}_{00}^{s_{1},s_{1}}\hat{C}_{00}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{00}^{s_{1},s_{1}}\hat{C}_{01}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{00}^{s_{1},s_{1}}\hat{C}_{10}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{00}^{s_{1},s_{1}}\hat{C}_{11}^{s_{1},s_{1}}]\\\
\text{cov}[\hat{C}_{01}^{s_{1},s_{1}}\hat{C}_{00}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{01}^{s_{1},s_{1}}\hat{C}_{01}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{01}^{s_{1},s_{1}}\hat{C}_{10}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{01}^{s_{1},s_{1}}\hat{C}_{11}^{s_{1},s_{1}}]\\\
\text{cov}[\hat{C}_{10}^{s_{1},s_{1}}\hat{C}_{00}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{10}^{s_{1},s_{1}}\hat{C}_{01}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{10}^{s_{1},s_{1}}\hat{C}_{10}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{10}^{s_{1},s_{1}}\hat{C}_{11}^{s_{1},s_{1}}]\\\
\text{cov}[\hat{C}_{11}^{s_{1},s_{1}}\hat{C}_{00}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{11}^{s_{1},s_{1}}\hat{C}_{01}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{11}^{s_{1},s_{1}}\hat{C}_{10}^{s_{1},s_{1}}]&\text{cov}[\hat{C}_{11}^{s_{1},s_{1}}\hat{C}_{11}^{s_{1},s_{1}}]\\\
\end{array}\right]\times\left[\begin{array}[]{c}\hat{C}_{00}^{s_{1},s_{1}}\\\
\hat{C}_{01}^{s_{1},s_{1}}\\\ \hat{C}_{10}^{s_{1},s_{1}}\\\
\hat{C}_{11}^{s_{1},s_{1}}\end{array}\right].$

We checked that changing $F$ from eq. (4.1) to eq. (4.3) produces fully
consistent results. All the results described in sec. 5.1 and 5.2 are obtained
using $F$ from eq. (4.1).

###  4.2 Survey specifications and analysis framework

Auto- and cross- angular power spectra used to get results in sec. 5 are
tomographically computed for the different tracer combinations reported in
tab. 1. For each survey we assume the specifications described in tab. 2: for
SN, we rely on the VRO survey, while for mergers we check both the single ET
scenario and the ET$\times$3 scenario that we introduced in [33]. The only
difference between the two resides in the sky localization uncertainty, which
in the ET$\times$3 case allows us to probe smaller scales; it is interesting
to explore whether similar sky localization improvements could possibly be
obtained using statistical techniques or cross-correlations with galaxies.

Analysis type | Survey | Tracer
---|---|---
Single tracer | VRO | SN
ET | DNS, BHNS, DBH
|  ET$\times$3 | DNS, BHNS, DBH
Multitracer |  ET$\times$3 | DNS + DBH
ET $+$ VRO | DBH + SN
|  ET$\times$3 $+$ VRO | DBH + SN
Table 1: Single and multitracer combinations considered in the analysis. As sec. 5.1 describes, VRO is the Vera Rubin Observatory, ET stands for Einstein Telescope, and ET$\times$3 for a network network made by 3 ET detectors. For each tracer the survey specifications are reported in tab. 2. BHNS are not included in the multitracer analysis, as sec. 5.2 explains. Survey | Source |  Area $[$deg${}^{2}]$ | $\Delta D_{L}/D_{L}$ |  $\Delta\Omega$ $[$deg${}^{2}]$ | $z_{max}$ | $T^{OBS}$ | $N^{TOT}$
---|---|---|---|---|---|---|---
VRO | SN | $18000$ | $0.15$ | $\sim 0$ | $1$ | 5yr | $10^{4.86}$
ET | DNS | Full sky | $0.3$ | $100$ | $2$ | 3yr | $10^{4.14}$
BHNS | $0.3$ | $3$ | $10^{4.37}$
DBH | $0.1$ | $5$ | $10^{4.79}$
ET$\times$3 | DNS | Full sky | $0.3$ | $10$ | $2$ | 3yr | $10^{4.14}$
BHNS | $0.3$ | $10$ | $3$ | $10^{4.37}$
DBH | $0.1$ | $3$ | $5$ | $10^{4.79}$
Table 2: Survey specifications (specs): observed sky area, luminosity distance
and sky localization uncertainties ($\Delta D_{L}/D_{L}$, $\Delta\Omega$,
respectively), detector horizon $z_{max}$. VRO specs are in accordance with
[20]; the sky localization is extremely precise for SN, therefore we assume
$\Delta\Omega\sim 0$. ET specs are compatible with [54] (full description
depending on frequency and distance for DNS and BHNS) and [49] (full analysis
for DBH), while for ET$\times$3 we assumed the same specs as in [33]; the only
difference between ET and ET$\times$3 relies in $\Delta\Omega$. $N^{TOT}$ is
computed as described in sec. 2 assuming $T^{OBS}$.

At each redshift, we consider only modes in the linear regime. More
specifically we choose a cut-off:

$l_{max}(z_{i},z_{j})=k_{nl}^{0}\min\bigl{[}(1+z_{i,j})^{2/(2+n_{s})}\chi(z_{i,j})\bigr{]}\
,$ (4.6)

where $k_{nl}^{0}=0.1\ h\text{Mpc}^{-1}$ is the non-linear cut-off scale at
$z=0$. We finally account for the sky localization uncertainty $\Delta\Omega$
by considering for each couple of tracers $s_{a},\ s_{b}$, a Gaussian beam
profile with $\sigma_{s_{a,b}}=\sqrt{\Delta\Omega_{s_{a,b}}/8\ln{2}}$, with
$\Delta\Omega_{s_{a,b}}$ from tab. 2 and converted to sterradians:

$B^{s_{a}s_{b}}_{l}=\exp{\biggl{[}{-\frac{l(l+1)}{2}\sigma_{s_{1}}^{2}}\biggr{]}}\exp{\biggl{[}{-\frac{l(l+1)}{2}\sigma_{s_{2}}^{2}}\biggr{]}}\
.$ (4.7)

In the case of SN, since the sky localization is extremely precise compared to
GW observations, we assume $\Delta\Omega_{SN}\sim 0$. Fig. 6 shows the
smoothing eq. (4.7) provides on the different surveys.

Figure 6: Beam window defined through eq. (4.7) for SN (green dashed) and
sources observed through ET (red for DNS, blue and overlapping for DBH, the
same is used for BHNS). The ET$\times$3 scenario is described by the red
dotted (DNS) and blue dotted (DBH) lines. Each line is computed considering
$l_{max}$ as the one associated through eq. (4.6) to $z_{max}$ reached for the
source (see tab. 2).

###  4.3 Parameters of interest

In our Fisher analysis, we consider the following set of parameters (i.e.
$\alpha,\beta$ in eq. (4.3)):

  * $\bullet$

Cosmological parameters: $[H_{0}$, $\Omega_{c}h^{2}$, $w_{0}$, $w_{a}$,
$\Omega_{b}h^{2}$, $n_{s}$, $A_{s}]$. Fiducial values are the ones from Planck
2018 [18].

  * $\bullet$

A pair of bias related parameters $[A_{i},P_{i}]$ for each source $i=$ SN,
DNS, BHNS, DBH (see eq. (2.8)). $A_{i}$ and $P_{i}$ represent respectively the
amplitude and the slope of the bias. The parameter $B_{i}$ in eq. (2.8), which
describes the bias local value, is fixed to the value extracted from the
simulations. The reason for this choice is that, using as fiducial values the
ones described in sec. 2 and including both $A_{i}$, $B_{i}$ in our analysis
would make the constraints on the two fully degenerate. As described in sec.
5.3, in the case of DBH we also check the effect of including a Gaussian prior
on the slope parameter $P_{DBH}$ (owing to the fact that results from
simulations strongly indicate a linear dependence of the $b_{i}(z)$ on
redshift).

Derivatives are computed numerically for cosmological parameters, analytically
for bias.

##  5 Results

As tab. 1 shows, we analyse different surveys, both considering one tracer
(sec. 5.1) and combining many of them to apply the multitracer technique (sec.
5.2). The specifications assumed for each kind of survey are described in tab.
2. Marginalized errors are obtained for cosmological and bias parameters using
the Fisher matrix formalism described in sec. 4; uniform priors are assumed
for all the parameters, except in sec. (5.3) where the effects of a prior on
the bias slope are analysed. We checked that including Planck 2018 [18] priors
on cosmological parameters does not produce a significant improvement in the
final constraints in the multitracer analysis with ET$\times$3, whereas in a
single tracer GW analysis, the Planck prior would simply dominate over the
signal.

###  5.1 Single tracer analysis

Tab. 3 and tab. 4 show the forecasted marginalized $1-\sigma$ errors for
cosmological and bias parameters, respectively.

BHNS mergers do not produce any significant improvement over DNS constraints,
despite the slightly larger number of these events (cfr. tab. 2) at slightly
higher redshifts ($2<z<3$). For this reason, we do not explicitly report BHNS
forecasts in tab. 3 and 4 and we do not include them in the multitracer
analysis in the next section.

| $H_{0}$ | $\Omega_{c}h^{2}$ | $w_{0}$ | $w_{a}$
---|---|---|---|---
SN | $7.864$ | $1.168\cdot 10^{-2}$ | $1.070\cdot 10^{-1}$ | $9.828\cdot 10^{-1}$
DNS ET$\times$3 | $76.67$ | $1.470\cdot 10^{-2}$ | $1.798$ | $12.23$
DBH ET$\times$x3 | $16.28$ | $2.446\cdot 10^{-2}$ | $2.552\cdot 10^{-1}$ | $1.622$

| $\Omega_{b}h^{2}$ | $n_{s}$ | $A_{s}$
---|---|---|---
SN | $1.468\cdot 10^{-3}$ | $8.924\cdot 10^{-2}$ | $3.315\cdot 10^{-10}$
DNS ET$\times$3 | $2.651\cdot 10^{-2}$ | $3.613\cdot 10^{-1}$ | $2.947\cdot 10^{-9}$
DBH ET$\times$3 | $5.134\cdot 10^{-3}$ | $1.783\cdot 10^{-1}$ | $5.976\cdot 10^{-10}$

Table 3: Marginalized $1-\sigma$ errors for cosmological parameters, single
tracer analysis.

As for cosmological parameters, the low-precision sky localization assumed for
ET prevents us from having good constraints on cosmological parameters for
both DNS and DBH. As we also discussed in [33], considering ET$\times$3 leads
to significant improvements, even though the overall constraining power
remains relatively small, compared to other cosmological probes. Forecasts for
SN are much better instead, since the VRO allows us to probe smaller (linear)
scales in the analysis, compared to GW. Fig. 8 compares confidence ellipses
for DBH, SN and for the multitracer analysis (see sec. 5.2) in constraining
the cosmological parameters (we marginalized over $H_{0}$, for which our
analysis have quite small constraining power).

We now turn our attention to bias parameters ($A_{i}$, $P_{i}$), and show in
fig. 7 the uncertainties which they produce on the estimates of the bias
coefficients, $b_{i}(z)$. Since the bias for each source is computed as
$b_{i}(z)=A_{i}(z+B_{i})^{P_{i}}$ (where $B_{i}$ is held fixed at 2.68 for
DBH, 2.46 for DNS and 0 for SN, as described in sec. 2), $\sigma_{b_{i}}$ is
computed by error propagation:

$\sigma_{b_{i}}^{2}=\begin{bmatrix}\dfrac{\partial b_{i}}{\partial
A_{i}}&\dfrac{\partial b_{i}}{\partial
P_{i}}\end{bmatrix}\begin{bmatrix}\sigma_{A_{i}}^{2}&\text{cov}(A_{i},P_{i})\\\
\text{cov}(A_{i},P_{i})&\sigma_{P_{i}}^{2}\end{bmatrix}\begin{bmatrix}\dfrac{\partial
b_{i}}{\partial A_{i}}\\\ \dfrac{\partial b_{i}}{\partial P_{i}}\end{bmatrix}\
.$ (5.1)

The values of $\sigma_{A_{i}}$ and $\sigma_{P_{i}}$ are described in tab. 4.
The single tracer analysis provides very good constraints for SN bias, since
its errors relative to the fiducial is $\sigma_{b_{SN}(z)}/b_{SN}\lesssim 7\%$
over all the redshifts considered. As for GW events instead, both ET and
ET$\times$3 forecasts are a bit worse with respect to our previous work [33]:
this depends on the degeneracy between the amplitude and slope parameters;
this is further studied in sec. 5.3, where ways to lift such degeneracy are
also discussed.

| $A_{DBH}$ | $P_{DBH}$ | $A_{DNS}$ | $P_{DNS}$ | $A_{SN}$ | $P_{SN}$
---|---|---|---|---|---|---
SN | $-$ | $-$ | $-$ | $-$ | $0.156$ | $0.020$
DNS | $-$ | $-$ | $12.80$ | $14.65$ | $-$ | $-$
DNS ET$\times$3 | $-$ | $-$ | $4.572$ | $5.014$ | $-$ | $-$
DBH | $1.446$ | $1.607$ | $-$ | $-$ | $-$ | $-$
DBH ET$\times$3 | $0.382$ | $0.404$ | $-$ | $-$ | $-$ | $-$
Table 4: Marginalized $1-\sigma$ errors for bias parameters computed, single
tracer analysis. Figure 7: Bias marginalized $1-\sigma$ errors at low $z$ for
SN (green dashed, upper panel), DNS (red, central panel) and DBH (blue, lower
panel) single tracer analysis. The bold line represent the fiducial models
described in sec. 2, while the shaded areas are included between
$b_{i}\pm\sigma_{b_{i}}$, with $\sigma_{b_{i}}$ from eq. (5.1). For DNS and
DBH both the ET (lighter) and ET$\times$3 (darker) constraints are showed.

###  5.2 Multitracer analysis

In our multitracer analysis we consider different combinations of tracers and
detectors, to test their impact on the final forecasts. The use of two or more
tracers makes it possible to cover a larger redshift range (with the density
distributions of different sources peaking at different $z$), add cross-power
spectrum information, break degeneracies and lower the effect of cosmic
variance, through the different bias dependence on redshift [45].

Tab. 5 and 6 show the forecasted $1-\sigma$ marginalized errors obtained
respectively for cosmological and bias parameters. We decided not to show the
analysis for DNS$+$BHNS$+$ DBH$+$SN since the DNS and BHNS distributions do
not introduce noticeable improvements.

Figure 8: Confidence $1-\sigma$ ellipses for cosmological parameters computed
for DBH (ET$\times$3, blue), SN (green) and DBH(ET$\times$3)$+$SN (black).
$H_{0}$, for which DBH$+$SN have quite poor constraining power, has been
marginalized.

| $H_{0}$ | $\Omega_{c}h^{2}$ | $w_{0}$ | $w_{a}$
---|---|---|---|---
DNS $+$ DBH (ET$\times$3) | $7.472$ | $1.987\cdot 10^{-2}$ | $1.938\cdot 10^{-1}$ | $1.035$
DBH (ET) $+$ SN | $6.386$ | $8.378\cdot 10^{-3}$ | $9.135\cdot 10^{-2}$ | $1.158\cdot 10^{-1}$
DBH (ET$\times$3) $+$ SN | $3.012$ | $6.983\cdot 10^{-3}$ | $6.701\cdot 10^{-2}$ | $2.532\cdot 10^{-1}$

| $\Omega_{b}h^{2}$ | $n_{s}$ | $A_{s}$
---|---|---|---
DNS $+$ DBH (ET$\times$3) | $3.772\cdot 10^{-3}$ | $1.603\cdot 10^{-1}$ | $3.405\cdot 10^{-10}$
DBH (ET) $+$ SN | $2.108\cdot 10^{-3}$ | $6.716\cdot 10^{-2}$ | $1.023\cdot 10^{-10}$
DBH (ET$\times$3) $+$ SN | $5.675\cdot 10^{-4}$ | $5.554\cdot 10^{-2}$ | $5.842\cdot 10^{-11}$

Table 5: Marginalized $1-\sigma$ errors for cosmological parameters,
multitracer analysis.

On one hand, the DNS$+$DBH(ET$\times$3) scenario only slightly improves the
forecasts with respect to sec. 5.1: this happens because DNS and DBH bias
models are close one to the other (cfr fig. 2) and therefore their combination
can not disentangle the power spectra dependence on the different parameters.
On the other hand, when including SN, forecasts largely improve. The SN
contribution is dominant with respect to the DBH one when considering a single
ET detector, making a joint analysis in this case not very meaningful:
forecasts for DBH(ET)$+$SN are in fact very close to those obtained in the
single-tracer, SN case. If we account for ET$\times$3, instead, the parameter
sensitivity achievable with either DBH or SN alone is more comparable.
Therefore, the multitracer analysis allows us in this case to significantly
improve the final constraints, over the single-tracer scenario. As a
benchmark, we can compare results from the multitracer DBH(ET$\times$3)$+$SN
case with those obtainable for galaxy clustering in a future Euclid-like
survey [10]. Our multitracer forecasts are comparable with the Euclid-like
ones for most of the parameters, the main exceptions being $H_{0}$ – for which
GW$+$SN have quite poor constraining power – and dark energy parameters
$w_{0}$ and $w_{a}$ – for which we obtain somewhat tighter constraints in the
GW $+$ SN analysis, thanks to the small scales reached by SN combined with the
large volumes probed by DBH.

If we now turn our attention to bias parameter constraints, we see again
significant improvements coming from the multitracer approach, considering
both ET and ET$\times$3\. Tab. 6 shows that the constraining power on $A_{i}$
and $P_{i}$ (hence on $b_{i}(z)$) improves considerably with respect to the
single tracer scenario, especially when we consider DBH$+$SN analysis (the
same would happen for DNS$+$SN). In particular, in the DBH(ET)$+$SN scenario,
$\sigma_{b_{DBH}(z)}/b_{DBH}\lesssim 35\%$ over all the distribution, while
for DBH(ET$\times$3)$+$SN $\sigma_{b_{DBH}(z)}/b_{DBH}\lesssim 15\%$. Also in
the case of bias parameters, since DBH and DNS bias fiducial models are very
similar, the DNS+DBH combination provides forecasts for $A_{DBH},\ P_{DBH},\
A_{DNS},\ P_{DNS}$ that are similar to those obtained with the single tracer
analysis in tab. 4.

| $A_{DBH}$ | $P_{DBH}$ | $A_{DNS}$ | $P_{DNS}$ | $A_{SN}$ | $P_{SN}$
---|---|---|---|---|---|---
DNS $+$ DBH (ET$\times$3) | $0.303$ | $0.322$ | $0.379$ | $0.456$ | $-$ | $-$
DBH (ET) $+$ SN | $0.308$ | $0.384$ | $-$ | $-$ | $0.085$ | $0.018$
DBH (ET$\times$3) $+$ SN | $0.126$ | $0.145$ | $-$ | $-$ | $0.058$ | $0.013$
Table 6: Marginalized $1-\sigma$ errors for bias parameters, multitracer
analysis. Figure 9: Bias error forecasts for SN and DBH in the DBH(ET)$+$SN
and DBH(ET$\times$3)$+$SN cases; the legend is the same of fig. 8.

###  5.3 Effects of a prior on bias parameters

The analysis of merger bias is particularly interesting to understand
clustering properties of binary mergers and disentangle different formation
scenarios. Since the ET catalogs will be dominated by DBH observations, we are
particularly interested in understanding the detectability of their bias. As
described in sec. 5.1 and 5.2, the bias error $\sigma_{b_{DBH}}$ depends on
both the uncertainties on the amplitude $A_{DBH}$ and on the slope $P_{DBH}$
parameters. The two, however, are degenerate one with respect to the other and
the degeneracy worsens the capability of this kind of analysis in detecting
the bias.

This can be checked by comparing single tracer bias forecasts described in
sec. 5.1 with results from our previous paper [33]. In our previous work, the
use of different bias parameters in each $D_{L}$ bin allowed us to detect bias
at low redshift. Using the $[A_{DBH},P_{DBH}]$ parameterization both in the
single ET and the ET$\times$3 makes $b_{DBH}(z)$ however always undetectable
in the single-tracer part of the current analysis. This issue can actually be
effectively solved by using the multitracer technique, as sec. 5.2 shows.
Moreover, we note that in sec. 5.1 and 5.2 we assume uninformative priors on
both $A_{DBH}$ and $P_{DBH}$. However, in our analysis in sec. 2 the slope
parameter turned out to be well described by a linear trend: we obtained this
result through an HOD analysis of hydrodynamical simulations [5, 6], but also
other works in the literature, despite following different approaches, found
the same behaviour, at least up to $z\sim 3.5$ (see e.g. [43, 38]). For this
reason, it seems reasonable to increase the level of "reliability" of the
fiducial value $P_{DBH}=1$ by associating a Gaussian prior to it.

Tab. 7 shows the effect of including different $P_{DBH}$ priors in the single
and full multitracer cases. All the measurements refer to the single ET or
ET$\times$3 scenarios. As expected, the more information the survey already
provides, the less impact the prior inclusion has in the analysis. Fig. 10
shows how the error propagation through eq. (5.1) changes depending on the
$P_{DBH}$ prior in the DBH and DBH+SN cases, both considering a single ET
detector. Assuming $\sigma_{P_{DBH}}^{prior}/b_{DBH}=50\%$ in the prior allows
us to achieve $\sigma_{b_{DBH}}/b_{DBH}\lesssim 50\%$ in the single-tracer DBH
case, while in the DBH(ET)$+$SN case we obtain
$\sigma_{b_{DBH}}/b_{DBH}\lesssim 25\%$, over the entire redshift range. These
constraints obviously further improve considering ET$\times$3, as tab. 7
shows.

$\sigma_{P_{DBH}}^{prior}/P_{DBH}$ | tracer | ET | ET$\times$3
---|---|---|---
$A_{DBH}$ | $P_{DBH}$ | $A_{DBH}$ | $P_{DBH}$
$100\%$ | DBH | $1.446$ | $1.607$ | $0.382$ | $0.404$
DBH $+$ SN | $0.308$ | $0.384$ | $0.126$ | $0.145$
$80\%$ | DBH | $0.663$ | $0.716$ | $0.341$ | $0.361$
DBH $+$ SN | $0.278$ | $0.346$ | $0.124$ | $0.142$
$50\%$ | DBH | $0.461$ | $0.477$ | $0.299$ | $0.314$
DBH $+$ SN | $0.245$ | $0.305$ | $0.121$ | $0.139$
$20\%$ | DBH | $0.249$ | $0.198$ | $0.176$ | $0.179$
DBH $+$ SN | $0.146$ | $0.177$ | $0.103$ | $0.117$
Table 7: Effects of different $P_{DBH}$ Gaussian priors on the marginalized
errors of bias parameters. The first line sums up results from tab. 4 and 6,
where uninformative, uniform prior is considered.

Figure 10: DBH bias error obtained in the DBH, single ET case (upper plot) and
in the DBH(ET)$+$SN case (lower plot) assuming the different priors on
$P_{DBH}$ described in tab. 7. The darker area is related to the larger prior.

##  6 Conclusions

In this work, we studied the power spectrum of different cosmic tracers in
Luminosity Distance Space. Our goal was that of studying the constraining
power of these observables on cosmological and bias parameters, by means of a
Fisher matrix approach. We considered future GW and SN observations, such as
those which will be made by Einstein Telescope (either in a single detector
configuration or in a network with $3$ detectors) and the Vera Rubin
Observatory.

In a previous analysis [33] we performed a similar forecast, but we considered
only GW mergers, in a single tracer approach. Consistently with our previous
analysis, we showed that the lack of sky localization in the single ET
scenario prevents us from observing scales small enough to provide good
constraints on cosmological parameters, while ET$\times$3 can do much better.
Using a multitracer analysis in which Supernovae are combined with GW
observations from ET$\times$3 leads to further significant improvements.888In
the single ET scenario, SN have instead a completely dominant constraining
power over GW, making the multitracer approach not so useful. The main source
of such improvements can be traced to the breaking of significant degeneracies
between different bias coefficients and cosmological parameters. This happens
because the information encoded in the SN and DBH datasets are complementary:
SN probe lower redshifts, while DBH probe large volumes and very high
redshifts. In the majority of cases, DNS mergers give a less important
contribution to this type of analysis, due to the lower resolution in $D_{L}$
of their dataset and to the fact that DNS do not probe the high redshift
region, where only DBH can be found. Moreover, the degeneracy between DNS and
DBH bias makes the combination DNS$+$DBH only slightly more predictive,
compared to the single-tracer analysis for either of them.

Our main results for cosmological parameters are shown in tab. 5. If, as a
benchmark test, we compare the DBH(ET$\times$3)$+$SN multitracer results to
the forecasts for a future Euclid-like survey for galaxy clustering, reported
in [10], we find similar constraining power in for most parameters, with the
exception of $H_{0}$ – which is constrained poorly by DBH(ET$\times$3)$+$SN –
and of dark energy parameters, $w_{0}$ and $w_{a}$ – which are instead
constrained better by DBH(ET$\times$3)$+$SN. These results might seem a little
surprising at first sight, because the expected number of galaxies in
forthcoming surveys, such as Euclid, is expected to be significantly larger
than the number of GW mergers or SN, even considering third generation GW
experiments (as we are doing here with ET). However, this is compensated by
the much larger volumes covered in GW merger surveys, especially when
considering DBH events, combined with the increase in the number of tracers
when we combine DBH and SN. Independently of the relative constraining power
of galaxy and GW merger clustering, the most interesting aspect is the
complementarity between different datasets. Moreover, the analysis of merger
and SN clustering in LDS presents specific peculiarities, compared to the
redshift space analysis of galaxy clustering, which are further discussed
below and in sec. 3.

In a cosmological parameter analysis, bias parameters are just treated as
nuisance parameters and marginalized out. GW merger bias coefficients,
however, are interesting per se, since they can provide insights into the
nature of the physical processes underlying the mergers themselves. For this
reason, a study of the detectability of $b_{i}(z)$, with $i=$ DBH, DNS, BHNS,
after marginalizing out the cosmological parameters, is also interesting [43].
Our fiducial model for merger bias is based on the outputs of hydrodynamical
N-body simulations [5, 6], processed via an HOD approach (see sec. 2.1 and
[33]). The main results of the bias forecast analysis have been summarized in
fig. 10: in the case of GW sources, when considering the single tracer
analysis the degeneracy between the amplitude $A_{i}$ and slope $P_{i}$
parameters makes the bias detection challenging. We showed however that in the
case of DBH both the multitracer technique (sec. 4) and the inclusion of a
Gaussian prior on the slope (sec. 5.3) make the bias well detectable even at
high redshift. This improves the results of our previous work [33], where bias
turned out to be detectable only at low redshift.

Besides the inclusion of SN and the combination of different tracers, another
significant new aspect of our analysis is the inclusion of lensing
contributions in the treatment of distortions in Luminosity Distance Space. In
our previous analysis, [33] we neglected lensing terms in LDSD because,
according to [53], they were expected to be subdominant with respect to to
peculiar velocity contributions. In the meantime, an explicit calculation of
lensing terms in the LDSD expression at first order was performed in [39],
where a numerical evaluation of the resulting GW merger power spectra was
obtained in Limber approximation and it was shown that the lensing
contribution was indeed smaller, but not negligible with respect to the
peculiar velocity one, especially on large scales. In this work, we modified
the publicly available CAMB code [14] to implement the LDSD formula of [39],
including lensing terms and without resorting to any flat-sky approximation.
The effects of these additional terms were then propagated to the final
parameter forecasts, showing that their impact is indeed very small for this
specific analysis. A full treatment of all GR effects on very large scales, in
luminosity distance space, is still missing in the literature and will be
pursued in a forthcoming publication [9]. When these will be available, it
will be possible to study GW surveys properly even on the largest scales: by
doing so, it will be interesting to test the forecasting power of ET in
measuring the non-Gaussianity parameter $f_{NL}$, via the merger power
spectrum at low-$k$. This would be in fact another suitable target for this
kind of surveys: since only large scales are required, on one hand, full sky
observations are needed, while, on the other, it is not necessary to have a
good sky localization. We are currently developing studies on the constraining
power of ET on $f_{NL}$ and we will provide our results in a future, dedicated
work.

## Acknowledgments

The authors thank D. Bertacca, A. Ricciardone and A. Raccanelli for useful
discussions. We thank A. Ghosh for the internal review provided for the
LIGO/Virgo Collaboration. ML was supported by the project "Combining Cosmic
Microwave Background and Large Scale Structure data: an Integrated Approach
for Addressing Fundamental Questions in Cosmology", funded by the MIUR
Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017 - grant
2017YJYZAH. ML, NB and SM acknowledge partial financial support by ASI Grant
No. 2016-24-H.0. MM, YB and NG acknowledge financial support from the European
Research Council for the ERC Consolidator grant DEMOBLACK, under contract no.
770017. MCA and MM acknowledge financial support from the Austrian National
Science Foundation through FWF stand-alone grant P31154-N27. DK acknowledge
financial support from the South African Radio Astronomy Observatory (SARAO)
and the National Research Foundation (Grant No. 75415).

## References

  * Abbott et al. [2021a] R. Abbott et al.  GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, 2021a.
  * Abbott et al. [2021b] R. Abbott et al.  GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run.  _Physical Review X_ , 11(2), 2021b.  doi: 10.1103/PhysRevX.11.021053.
  * Abbott et al. [2021c] R. Abbott et al.  Observation of gravitational waves from two neutron star–black hole coalescences.  _The Astrophysical Journal Letters_ , 915(1), 2021c.  doi: 10.3847/2041-8213/ac082e.
  * Abramo et al. [2016] L. Abramo et al.  Fourier analysis of multitracer cosmological surveys.  _Mon. Not. Roy. Astron. Soc._ , 455(4), 2016.
  * Artale et al. [2018] M.C. Artale et al.  The impact of assembly bias on the halo occupation in hydrodynamical simulations.  _Monthly Notices of the Royal Astronomical Society_ , 480(3):3978–3992, 2018.
  * Artale et al. [2019] M.C. Artale et al.  Mass and star formation rate of the host galaxies of compact binary mergers across cosmic time.  _Monthly Notices of the Royal Astronomical Society_ , 491(3):3419–3434, 2019.
  * Belgacem et al. [2018a] E. Belgacem et al.  Modified gravitational-wave propagation and standard sirens.  _Physical Review D_ , 98(2), 2018a.
  * Belgacem et al. [2018b] E. Belgacem et al.  Gravitational-wave luminosity distance in modified gravity theories.  _Physical Review D_ , 97(10), 2018b.
  * Bertacca et al. [in prep.] D. Bertacca et al.  , in prep.
  * Blanchard et al. [2020] A. Blanchard et al.  Euclid preparation.  _Astronomy $\&$ Astrophysics_, 642, 2020.  doi: 10.1051/0004-6361/202038071.
  * Breton et al. [2018] M. A. Breton et al.  Imprints of relativistic effects on the asymmetry of the halo cross-correlation function: from linear to non-linear scales.  _Monthly Notices of the Royal Astronomical Society_ , 483(2), 2018.  doi: 10.1093/mnras/sty3206.
  * Calore et al. [2020] F. Calore et al.  Cross-correlating galaxy catalogs and gravitational waves: A tomographic approach.  _Phys. Rev. Research_ , 2, 2020.  doi: 10.1103/PhysRevResearch.2.023314.
  * Carron [2013] J. Carron.  On the assumption of gaussianity for cosmological two-point statistics and parameter dependent covariance matrices.  _Astronomy $\&$ Astrophysics_, 551, 2013.  ISSN 1432-0746.  doi: 10.1051/0004-6361/201220538.
  * Challinor and Lewis [2011] A. Challinor and A. Lewis.  Linear power spectrum of observed source number counts.  _Physical Review D_ , 84(4), 2011.  ISSN 1550-2368.
  * Collaboration [2020a] LIGO Scientific Collaboration.  Prospects for observing and localizing gravitational-wave transients with advanced ligo, advanced virgo and kagra.  _Living Reviews in Relativity_ , 23(1), 2020a.
  * Collaboration and Collaboration [2016] LIGO Scientific Collaboration and Virgo Collaboration.  Observation of gravitational waves from a binary black hole merger.  _Phys. Rev. Lett._ , 116:061102, 2016.
  * Collaboration [2009] LSST Science Collaboration.  Lsst science book, version 2.0, 2009.
  * Collaboration [2020b] Planck Collaboration.  Planck 2018 results - vi. cosmological parameters.  _Astronomy and Astrophysics_ , 641, 2020b.  ISSN 1432-0746.
  * Cremonese et al. [2021] P. Cremonese, J. M. Ezquiaga, and V. Salzano.  Breaking the mass-sheet degeneracy with gravitational wave interference in lensed events.  _Physical Review D_ , 104(2), 2021.  doi: 10.1103/physrevd.104.023503.
  * Cullan et al. [2017] H. Cullan et al.  Measuring the growth rate of structure with type IA supernovae from LSST.  _The Astrophysical Journal_ , 847(2), sep 2017.
  * Dalang et al. [2020] C. Dalang et al.  Horndeski gravity and standard sirens.  _Physical Review D_ , 102(4), 2020.
  * Diaz and Mukherjee [2021] C. C. Diaz and S. Mukherjee.  Mapping the cosmic expansion history from LIGO-Virgo-KAGRA in synergy with DESI and SPHEREx.  _arXiv:2107.12787_ , 7 2021.
  * Garcia et al. [2020] K. Garcia et al.  On the amount of peculiar velocity field information in supernovae from lsst and beyond.  _Physics of the Dark Universe_ , 29:100519, 2020.
  * Garoffolo et al. [2020] A. Garoffolo et al.  Gravitational waves and geometrical optics in scalar-tensor theories.  _Journal of Cosmology and Astroparticle Physics_ , 2020(11), 2020.
  * Garoffolo et al. [2021] A. Garoffolo et al.  Detecting dark energy fluctuations with gravitational waves.  _Physical Review D_ , 103(8), 2021.
  * Hamaus et a.l. [2011] N. Hamaus et a.l.  Optimal constraints on local primordial non-gaussianity from the two-point statistics of large-scale structure.  _Phys. Rev. D_ , 84, 2011.
  * Hamilton [1998] A. J. S. Hamilton.  Linear redshift distortions: A review.  _The Evolving Universe_ , 1998.  doi: 10.1007/978-94-011-4960-0_17.
  * Hui and Greene [2006] L. Hui and P.B. Greene.  Correlated fluctuations in luminosity distance and the importance of peculiar motion in supernova surveys.  _Physical Review D_ , 73(12), 2006.
  * Ilić et al. [2021] S. Ilić et al.  Euclid preparation: Xv. forecasting cosmological constraints for the Euclid and cmb joint analysis, 2021.
  * Jensen [2004] J.W. Jensen.  Supernovae light curves: An argument for a new distance modulus, 2004\.
  * Kaiser [1987] N. Kaiser.  Clustering in real space and in redshift space.  _Monthly Notices of the Royal Astronomical Society_ , 227(1):1–21, 07 1987.  ISSN 0035-8711.  doi: 10.1093/mnras/227.1.1.
  * Khokhlov et al. [1993] A. Khokhlov, E. Mueller, and P. Hoeflich.  Light curves of type IA supernova models with different explosion mechanisms.  _Astronomy and Astrophysics_ , 270:223–248, 1993.
  * Libanore et al. [2021] S. Libanore et al.  Gravitational wave mergers as tracers of large scale structures.  _Journal of Cosmology and Astroparticle Physics_ , 2021(02):035–035, 2021.
  * Limber [1953] D. N. Limber.  The analysis of counts of the extragalactic nebulae in terms of a fluctuating density field.  _The Astrophysical Journal_ , 117, 1953.
  * Maggiore et al. [2020] M. Maggiore et al.  Science case for the Einstein telescope.  _ournal of Cosmology and Astroparticle Physics_ , 2020(3), 2020.  doi: 10.1088/1475-7516/2020/03/050.
  * McDonald and Seljak [2009] P. McDonald and U. Seljak.  How to evade the sample variance limit on measurements of redshift-space distortions.  _Journal of Cosmology and Astroparticle Physics_ , 2009(10), 2009.
  * Mukherjee and Wandelt [2018] S. Mukherjee and B. D. Wandelt.  Beyond the classical distance-redshift test: cross-correlating redshift-free standard candles and sirens with redshift surveys.  _arXiv:1808.06615_ , 8 2018.
  * Mukherjee et al. [2021] S. Mukherjee et al.  Accurate precision cosmology with redshift unknown gravitational wave sources.  _Physical Review D_ , 103(4), 2021.
  * Namikawa [2021] T. Namikawa.  Analyzing clustering of astrophysical gravitational-wave sources: luminosity-distance space distortions.  _Journal of Cosmology and Astroparticle Physics_ , 2021(01), 2021.  ISSN 1475-7516.
  * Namikawa et al. [2016] T. Namikawa et al.  Anisotropies of gravitational-wave standard sirens as a new cosmological probe without redshift information.  _Physical Review Letters_ , 116(12), 2016.
  * Sasaki [1987] M. Sasaki.  The magnitude-redshift relation in a perturbed friedmann universe.  _Monthly Notices of the Royal Astronomical Society_ , 228, 1987.
  * Scelfo et al. [2018] G. Scelfo et al.  GW×LSS: chasing the progenitors of merging binary black holes.  _Journal of Cosmology and Astroparticle Physics_ , 2018(09), 2018.
  * Scelfo et al. [2020] G. Scelfo et al.  Exploring galaxies-gravitational waves cross-correlations as an astrophysical probe.  _Journal of Cosmology and Astroparticle Physics_ , 2020(10), 2020.
  * Schaye et al. [2015] J. Schaye et al.  The EAGLE project: simulating the evolution and assembly of galaxies and their environments.  _Monthly Notices of the Royal Astronomical Society_ , 446(1), 2015.  doi: 10.1093/mnras/stu2058.
  * Seljak [2009] U. Seljak.  Extracting primordial non-gaussianity without cosmic variance.  _Physical Review Letters_ , 102(2), 2009.
  * Tinker et al. [2008] J. Tinker et al.  Toward a halo mass function for precision cosmology: The limits of universality.  _The Astrophysical Journal_ , 688(2):709–728, 2008.
  * Veitch et al. [2015] J. Veitch et al.  Parameter estimation for compact binaries with ground-based gravitational-wave observations using the lalinference software library.  _Physical Review D_ , 91(4), 2015.  ISSN 1550-2368.
  * Vijaykumar et al. [2020] A. Vijaykumar et al.  Probing the large scale structure using gravitational-wave observations of binary black holes, 2020.
  * Vitale and Evans [2017] S. Vitale and M. Evans.  Parameter estimation for binary black holes with networks of third-generation gravitational-wave detectors.  _Phys. Rev. D_ , 95, 2017.
  * Witzemann et al. [2019] A. Witzemann et al.  Simulated multitracer analyses with hi intensity mapping.  _Monthly Notices of the Royal Astronomical Society_ , 485(4), 2019.
  * Yang et al. [2019] T. Yang et al.  Constraints on the cosmic distance duality relation with simulated data of gravitational waves from the einstein telescope, 2019.
  * Yoo [2009] J. Yoo.  Complete treatment of galaxy two-point statistics: Gravitational lensing effects and redshift-space distortions.  _Physical Review D_ , 79(2), 2009.  doi: 10.1103/physrevd.79.023517.
  * Zhang [2018] P. Zhang.  The large scale structure in the 3d luminosity-distance space and its cosmological applications, 2018.
  * Zhao and Wen [2018] W. Zhao and L. Wen.  Localization accuracy of compact binary coalescences detected by the third-generation gravitational-wave detectors and implication for cosmology.  _Physical Review D_ , 97(6), 2018.