# Impact of weak lensing on bright standard siren analyses

Charlie T. Mpetha1<EMAIL_ADDRESS>Giuseppe Congedo1 Andy Taylor1 Martin A.
Hendry2 1Institute for Astronomy, School of Physics and Astronomy, University
of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, United
Kingdom
2SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow G12
8QQ, United Kingdom

###### Abstract

Gravitational waves from binary mergers at cosmological distances will
experience weak lensing by large scale structure. This causes a
(de-)magnification, $\mu$, of the wave amplitude, and a completely degenerate
modification to the inferred luminosity distance $d_{L}$. The customary method
to address this is to increase the uncertainty on $d_{L}$ according to the
dispersion of the magnification distribution at the source redshift,
$\sigma_{\mu}(z)$. But this term is dependent on the cosmological parameters
that are being constrained by gravitational wave “standard sirens”, such as
the Hubble parameter $H_{0}$, and the matter density fraction $\Omega_{\rm
m}$. The dispersion $\sigma_{\mu}(z)$ is also sensitive to the resolution of
the simulation used for its calculation. Tension in the measured value of
$H_{0}$ from independent datasets, and the present use of weak lensing fitting
functions calibrated using outdated cosmological simulations, suggest
$\sigma_{\mu}(z)$ could be underestimated. This motivates an investigation
into the consequences of mischaracterising $\sigma_{\mu}(z)$. We consider two
classes of standard siren, supermassive black hole binary and binary neutron
star mergers. Underestimating $H_{0}$ and $\Omega_{\rm m}$ when calculating
$\sigma_{\mu}(z)$ increases the probability of finding a residual lensing bias
on these parameters greater than $1\sigma$ by $1.5-3$ times. Underestimating
$\sigma_{\mu}(z)$ by using low resolution$/$small sky-area simulations can
also significantly increase the probability of biased results, the probability
of a $1\sigma\,(2\sigma)$ bias in $H_{0}$ and $\Omega_{\rm m}$ found from
binary neutron star mergers is $54\%\,(19\%)$ in this case. For neutron star
mergers, the mean bias on $H_{0}$ caused by magnification selection effects is
$\Delta H_{0}=-0.1\,$km s${}^{-1}\,$Mpc-1. The spread around this mean
bias—determined by assumptions on $\sigma_{\mu}(z)$—is $\Delta H_{0}=\pm
0.25\,$km s${}^{-1}\,$Mpc-1, comparable to the forecasted uncertainty. These
effects do not impact merging neutron stars’ utility for addressing the
$H_{0}$ tension, but left uncorrected they limit their use for precision
cosmology. For supermassive black hole binaries, the spread of possible biases
on $H_{0}$ is significant, $\Delta H_{0}=\pm 5\,{\rm km}\,{\rm s}^{-1}\,{\rm
Mpc}^{-1}$, but $\mathcal{O}(200)$ observations are needed to reduce the
variance below the bias. To achieve accurate sub-percent level precision on
cosmological parameters using standard sirens, first much improved knowledge
on the form of the magnification distribution and its dependence on cosmology
is needed.

††preprint: APS/123-QED

##  I Introduction

Waves propagating through the Universe are magnified by the gravitational
potential of large scale structure. This weak lensing effect on light,
observed through its effect on galaxy shapes, has been intensely studied as a
source of cosmological information [1, 2, 3, 4, 5]. Similarly to light,
gravitational wave weak lensing (GW-WL) modifies the observed wave amplitude
$h$, to $h^{\prime}=h\sqrt{\mu}$, where a prime denotes a lensed quantity [6,
7, 8, 9, 10]. The magnification $\mu$ depends on the matter along the line of
sight where $\mu=1$ corresponds to no effect, $\mu<1$ a net underdensity of
matter, and $\mu>1$ a net overdensity.

The magnification distribution over the sky, $p(\mu,z)$, broadens towards
higher redshifts, and is peaked at values below one as there are more voids
than clusters. The tail of $p(\mu,z)$ is weighted towards values greater than
one since clustered matter can have a larger effect on the magnification than
voids. Gravitational waves from low redshift merging binaries, such as those
observed using the LIGO-Virgo-Kagra (LVK) network [11, 12, 13, 14], will not
experience impactful weak lensing, though they may still be strongly lensed
[15, 16]. The space-based GW detector LISA [17] expects to observe between a
few and a few tens of super massive black hole binary (SMBHB) mergers in the
mid-2030’s over its nominal 5-year survey time with a duty cycle of $80\%$
[18, 19, 20]. The effect of weak lensing on these observations will be
significant, even for individual sources due to small uncertainties on their
measured distances. A future network of $3^{\rm rd}$ generation ground-based
detectors (3G), such as the Einstein Telescope [21] and Cosmic Explorer [22]
is also anticipated to begin operations in the mid-2030’s. These detectors
will observe binary neutron stars (BNSs) at $z\\!\sim\\!2$ [23]. Weak lensing
will not be as important for individual BNSs due to their larger distance
uncertainties compared to those of SMBHBs. On the other hand, with the 3G
network expected to detect many more BNSs than the numbers of SMBHBs observed
by LISA, the effect of instrumental scatter on the cosmological parameters
will be reduced. The bias due to lensing is likely to be even more
significant.

For a single GW observation, the weak lensing effect can not be estimated from
GW data alone. Past works have studied how successfully the GW can be
“delensed” using external data to estimate $\mu$ for each observed source [24,
25, 26]. Ref. [26] found that sufficient delensing using galaxy surveys for
all SMBHB mergers observed with LISA will be extremely challenging. For both
SMBHBs and BNSs, we can not expect to remove or even meaningfully reduce the
weak lensing effect, unless further advances are made in delensing methods.

Measuring the distance and redshift of a binary merger is the core of using
these sources for cosmology, as “standard sirens” [27]. Bright standard sirens
(BSS) are the ideal case where the GW is observed alongside an electromagnetic
counterpart. Another method for using standard sirens for cosmology utilises
the weak lensing of a large number of GWs as information instead of noise,
constraining structure as well as geometry [28, 29, 30, 31, 32, 33, 34]. But
this is only viable with sufficient source numbers, likely not until at least
the 2040s. In the meantime, many works have demonstrated how weak lensing can
bias the distance redshift relation for the case of supernovae standard
candles and GW standard sirens [35, 36, 37, 38, 39, 40, 41, 42].

The luminosity distance $d_{L}$ is an observable quantity with gravitational
waves from binary mergers, without needing any external calibration from the
cosmic distance ladder. At redshift $z$,

$d_{L}(z)=(1+z)\int_{0}^{z}\frac{c\,dz^{\prime}}{H(z^{\prime})}\,,$ (1)

where the Hubble parameter is

$\displaystyle\left(\frac{H(z)}{H_{0}}\right)^{2}$ $\displaystyle=\Omega_{\rm
m}(1+z)^{3}+\Omega_{K}(1+z)^{2}$ $\displaystyle\quad+\Omega_{\rm
DE}(1+z)^{3(1+w_{0}+w_{a})}e^{\frac{-3w_{a}z}{1+z}}\,.$ (2)

Here $H_{0}=100\,h\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$ is the present
expansion rate, $\Omega_{\rm m}$, $\Omega_{K}$ and $\Omega_{\rm DE}$ are the
matter, curvature and dark energy density fractions respectively, and $w_{0}$
and $w_{a}$ are equation of state parameters for dark energy [43, 44].

The luminosity distance is inversely proportional to the observed wave
amplitude. Therefore the observed, lensed distance is

$d_{L}^{\prime}=\frac{d_{L}}{\sqrt{\mu}}\,.$ (3)

The total distance uncertainty is

$\sigma_{d_{L}}^{2}=\sigma_{\rm GW}^{2}+\sigma_{\rm
WL}^{2}+\left(\frac{\partial d_{L}}{\partial
z}\right)^{2}\left(\sigma_{z}^{2}+\sigma_{v_{\rm pec}}^{2}\right)\,,$ (4)

where $\partial d_{L}/\partial z$ converts redshift to distance uncertainties.
$\sigma_{\rm GW}$ and $\sigma_{v_{\rm pec}}$ are the uncertainty from the
detector and source peculiar velocity respectively. The weak lensing term is

$\frac{\sigma_{\rm
WL}(z)}{d_{L}(z)}=\frac{1}{2}\frac{\sigma_{\mu}(z)}{\langle\mu(z)\rangle}=\frac{1}{2}\sigma_{\mu}(z)\,,$
(5)

$\sigma_{\mu}(z)$ is the uncertainty in the magnification of the gravitational
wave. It has been assumed that the mean of the magnification distribution
$\langle\mu\rangle=1$. For real observations this is in general not the case
due to detector selection effects as explored in Ref. [45], henceforth CT21.
The consequence is $\langle\mu\rangle\\!\sim\\!1.02$ at $z=2$.

Weak lensing could lead to a net over- or under-estimation of luminosity
distances at different (true) redshifts, so that the inferred cosmological
parameters in Eq. (2) are prone to biases [42]. A mitigation might be to
include a weak lensing uncertainty according to Eq. (5), using an assumed
functional form for $\sigma_{\mu}(z)$ [46, 47, 48, 45, 49, 50]. The possible
lensing bias is washed out by increasing distance errors. A caveat of this
added weak lensing uncertainty is that the form of the magnification function
is not well known, is dependent on the uncertain values of cosmological
parameters, and is affected by poorly understood physics influencing the
growth of structure on small scales. There is tension in the value of $H_{0}$
[51] (recently exacerbated by JWST observations [52]), and present weak
lensing fitting functions have been created with cosmological simulations that
could be missing large and small-scale power important in estimating
$\sigma_{\mu}(z)$ [53] (henceforth T11). Both of these facts motivate an
investigation into how under- or overestimating the weak lensing uncertainty
can impact a cosmological analysis using GWs.

We will evaluate the residual lensing bias, the lensing bias remaining even
after including a weak lensing uncertainty term. Ref. [42] found the residual
lensing bias for one idealised source population, demonstrating for the first
time that the lensing bias can be greater than the statistical uncertainty for
BNS mergers. Our goals differ from that work in that we are broadly assessing
the impact of mischaracterising the weak lensing uncertainty term for
realistic observations from future space and ground-based GW detectors.

We find that underestimating the weak lensing uncertainty by using incorrect
values for cosmological parameters leads to a factor $1.5$ $(3)$ increase in
the probability of obtaining results biased by more than $1\sigma$ for BNSs
(SMBHBs). We also demonstrate the importance of the resolution of simulations
used to determine $\sigma_{\mu}(z)$, pointing out that missing power due to
simulations not including very small and large scale magnification power could
mean up to a factor of $10$ greater probability of obtaining results biased by
more than $1\sigma$. For constraints from BNSs, we find the bias on $H_{0}$
caused by lensing is too small to affect their application to the $H_{0}$
tension. But the bias is at a similar level to the variance and so limits
their use for precision cosmology. For SMBHBs the lensing bias is significant
with a range $\sim\\!5\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$, but will only
become relevant with $\mathcal{O}(200)$ observations. If GWs from merging
binaries are to be used to their full potential by the next generation of
detectors, much improved knowledge of $\sigma_{\mu}(z)$ is needed.

Section II describes the method used to evaluate the residual lensing bias. In
Section III we introduce the binary merger catalogues used in our analysis.
Section IV details the generation of magnification distributions and weak
lensing uncertainties, and Section V gives our results. We conclude in Section
VI.

##  II Method

Our goal is to evaluate the probability of finding cosmological constraints
that are biased due to lensing, in particular the residual lensing bias caused
by underestimating the weak lensing uncertainty. The following steps outline
how we evaluate this residual lensing bias.

  1. 1.

Set a true cosmology by choosing values for cosmological parameters in Eq.
(2).

  2. 2.

Choose a population: supermassive black hole binaries (SMBHBs) with a
particular seed model observed with LISA [20] (see Section III.1), or binary
neutron stars (BNSs) observed with third generation ground based detectors
(see Section III.2). Steps $3-7$ are repeated over many catalogue
realisations: independent sets of observations drawn from the population
model.

  3. 3.

Choose a catalogue realisation. For each merger in the catalogue, randomly
draw a magnification $\mu$ as described in Section IV from the redshift and
cosmology dependent magnification distribution. Applying the magnification to
the measured merger distance and its uncertainty $\sigma_{\rm GW}$ (a more
magnified source has a higher SNR and smaller uncertainty, see Appendix A)
gives a lensed source, to be compared to the unlensed case. The unlensed case
has $\mu=1$ for all sources. Only sources with a SNR greater than a chosen
threshold after applying the magnification go into the catalogue.

  4. 4.

For each merger generate an observed distance and redshift by drawing from a
Gaussian distribution with a width equal to the instrumental uncertainty and
mean equal to the true distance$/$redshift. Also add scatter to the redshift
caused by peculiar velocities, motion of the binaries not due to Hubble flow.
This is done for both the lensed and unlensed sources, with each having the
same random draw to ensure we are only comparing the effects of lensing
($d_{L}^{\prime{\rm obs}}=d_{L}^{\rm obs}/\sqrt{\mu}$).

  5. 5.

In both the lensed and unlensed cases include a weak lensing uncertainty
$\sigma_{\rm WL}(z)$, as in Eq. (4), to the distance. The assumptions going
into generating $\sigma_{\rm WL}(z)$ are the focus of this investigation.

  6. 6.

Given the observed values of $d_{L}$ and $z$ for each merger in the catalogue,
perform Monte Carlo Markov Chain sampling of $d_{L}-z$ in Eq. (1) to
reconstruct cosmological parameters in the lensed and unlensed cases. In this
work we use the emcee package [54].

  7. 7.

The absolute bias is the difference in the best fit recovered parameter value
in the lensed and unlensed case. The bias in units of uncertainty is found
from the parameter differences method of the tensiometer
package111https://github.com/mraveri/tensiometer [55, 56]. This uses
information from the full $h-\Omega_{\rm m}$ 2D posterior to quantify the size
of the bias, and does not assume Gaussianity in the parameter space.

  8. 8.

Aggregate the biases for each catalogue realisation of a population. The
resulting distribution of biases from all realisations characterises the
residual lensing bias that can be expected for that population, given the true
cosmology and assumed weak lensing uncertainty.

  9. 9.

Repeat steps $3-8$ a further $\sim\\!300$ times, and combine the distributions
from each repetition. There is possible variation in observations due to the
random drawing of magnifications and random scatter caused by distance and
redshift uncertainties. By doing this, we have ensured we sample over both
statistical variation in the catalogues, and statistical variation in the
observations.

##  III Binary Merger catalogues

Our goal is to evaluate the residual lensing bias on cosmological parameters
constrained with GWs. To this end, realistic catalogues of observations by the
next generation of gravitational wave detectors are needed. For the following
statistical tests $100/500$ catalogue realisations are used for BSS observed
with LISA$/$3G. Section III.1 details the subset of SMBHB catalogues from Ref.
[20] used in this work, while the BNS catalogues we generated are described in
Section III.2. An example of one of the catalogues for each population,
including instrumental distance$/$redshift uncertainties and magnification and
instrumental scatter, can be seen in Fig. 1. The number of observations will
depend on many factors; the performance of future GW and EM detectors, the
fraction of the sky area covered by electromagnetic detectors for the
multimessenger follow-up, the merger rate, and properties of the EM
counterpart. For this reason we remain agnostic of rates and integration
times, and instead use a broad range of source numbers in each catalogue. We
choose not to consider dark standard sirens in this analysis, for example 3G
observed binary black holes, or extreme-mass-ratio inspirals observed with
LISA. The large redshift uncertainty of these sources will dominate over the
weak lensing uncertainty for most cases. While Ref. [57] demonstrate that
using several hundred 3G-observed ${\rm SNR}>300$ binary black holes could
achieve comparable precision on $H_{0}$ to bright standard sirens, most of
these will be too low redshift for important magnification effects.

All population catalogues have been generated assuming a Planck-like cosmology
[58]. In this work we evaluate the residual lensing bias in different true
cosmologies, and while the assumed cosmology would impact the source
distributions and uncertainties, we note that this effect is subdominant to
details on the seed model and host galaxy evolution.

Figure 1: Representative catalogues for each binary merger population used in
this work. Every gravitational wave has been observed alongside an
electromagnetic counterpart. Sources have instrumental uncertainty and scatter
in distance and redshift, and scatter caused by weak lensing and peculiar
velocity. Weak lensing causes large scatter towards high redshift. The popIII
light-seed SMBHB model typically predicts fewer, lower redshift mergers than
the Q3 heavy-seed models.

###  III.1 LISA: Supermassive black hole binaries

We use a subset of the catalogues presented in Refs [20, 59], refer there for
more detail on the catalogue generation, and to Ref. [59] for cosmological
forecasts using these catalogues. As is common in the literature, uncertainty
on the seed model for supermassive black holes motivates including three
possibilities; black holes are formed from high-redshift high-mass metal poor
“popIII” stars [60], the “Q3” heavy seed model where massive black holes are
formed from the collapse of proto-galactic disks [61, 62] with a time delay
between the formation time of a binary and its merger (Q3d), and the Q3 model
with no delay time between formation and merger (Q3nod). Only mergers with an
electromagnetic counterpart are considered. Other methods of estimating the
redshift of a single source often rely on galaxy catalogues which will not be
possible for high-redshift SMBHBs. In these catalogues there are three
possibilities for the redshift uncertainty depending on the magnitude and
redshift of the host galaxy observed with the ELT222https://elt.eso.org/. If
the host galaxy magnitude $m_{\rm gal,ELT}<27.2$ and $0.5<z_{\rm gal}\leq 5$
then $\sigma_{z}=5\times 10^{-4}$. If $27.2<m_{\rm gal,ELT}<31.3$ and
$0.5<z_{\rm gal}\leq 5$ ($z_{\rm gal}>5$) then $\sigma_{z}=0.1$
($\sigma_{z}=0.25$). See Ref. [20] for more detail. We also include a peculiar
velocity drawn from a Gaussian centred at zero with a width of $500\,$km s-1.
Sources (from all seed models) range from $z=0.8$ to $z=9.2$, with distance
uncertainties from $0.02\%$ to $30\%$ having a median value of $0.8\%/0.4\%$
for the popIII$/$Q3 seed models. The average number of sources (median
redshift) is 7 (2.4), 15 (3.2) and 20 (2.9) for popIII, Q3d and Q3nod
respectively.

###  III.2 3G: Binary neutron stars

For $3^{\rm rd}$ generation ground-based detectors (3G), the most likely BSS
are binary neutron star (BNS) mergers. It is not expected that many black
hole-neutron star mergers will be observed alongside an EM counterpart [63],
so this population is not included. BNSs are generated using the assumptions
in Table 1 of Ref. [64], which involve uniform sampling over sky location and
inclination, uniform sampling of spins aligned with the orbital angular
momentum in the range $[-0.05,0.05]$, and uniform sampling of mass in the
range $[1,2.5]\,M_{\odot}$. The adopted redshift distribution of BSS is
described in Ref. [32]. The Fisher-matrix code GWFAST [65] is used to find
measurement uncertainties of each source in the catalogue. We assume a network
of Einstein Telescope and two Cosmic Explorers (ET+2CE) for the main analysis,
considering other configurations in Appendix B.

$\mathcal{O}(10^{5})$ BNSs will be observed every year by a 3G detector
network [66, 64]. Predictions on the corresponding number of BSS vary
significantly, from 15 yr-1 to 3500 yr-1 [67, 23, 68, 69, 70, 71, 72]
depending on assumptions on the GW and EM detectors, merger rates and
counterparts. To reflect this our 500 generated catalogues have source numbers
sampled uniformly between $100$ and $1000$. The minimum is chosen to be $100$
as this could be realised even with the lowest rate estimates after 10 years
of runtime.

We impose a GW detector network SNR threshold of $\rho_{\rm lim}=12$ [73]. Due
to the requirement of an EM counterpart, these observations come with
selection effects beyond the requirements for a detection of the gravitational
wave alone. We add a cut on sky localisation of $20\,$deg2 [74, 75]. We also
limit the distance uncertainty to $\sigma_{d_{L}}/d_{L}<30\%$, around the
upper limit of the $z-\sigma_{d_{L}}/d_{L}$ fitting function in Ref. [76] for
BSS observations with ET. As pointed out in Ref. [10], knowledge of the sky
location and inclination of a merger from its EM counterpart can improve
constraints on GW parameters, quoting an order of magnitude improvement in
$\sigma_{d_{L}}/d_{L}$. More recent studies use information on the inclination
angle $\iota$ from EM observations to break the $d_{L}-\iota$ degeneracy and
improve $\sigma_{d_{L}}/d_{L}$ by several factors [77, 78, 79]. This was
performed for GW170817 using the radio counterpart, improving the $H_{0}$
constraint by a factor of two [80]. In the context of LISA, Ref. [81] show how
from sky localisation the wave amplitude and inclination uncertainties can
both be improved by a factor of $\sim\\!2$ when $\iota\nsim 90^{o}$. This
factor increases to $\sim\\!10$ if there is also some prior knowledge from EM
observations that constrains the inclination. Therefore we make the
simplifying assumption that, for all sources with $\iota\leq 60^{o}$ and
$\iota\geq 120^{o}$ remaining after making the SNR and sky localisation cuts,
the distance uncertainty is reduced by a factor of two using localisation
information from the EM counterpart. We note that there could also be a
constraint on the inclination angle for a fraction of observations, but do not
explicitly include this effect.

The final catalogues contain sources with redshifts from $0.018$ to $2$ and
instrumental $d_{L}$ uncertainties ranging from $0.008\%$ to $30\%$, with a
median value of $4\%$. Each merger has a peculiar velocity drawn from a
Gaussian centred at zero with width $500\,{\rm km}\,{\rm s}^{-1}$. Finally, we
assume spectroscopic redshift uncertainties of $\sigma_{z}=5\times
10^{-4}\,(1+z)$, based on electromagnetic follow up from present and
forthcoming spectroscopic redshift surveys e.g DESI [82], ELT, and MegaMapper
[83]. The impact of this uncertainty is negligible.

##  IV Magnification distributions

The following section details the redshift and cosmology dependent
magnification distribution, and how it is used to generate a weak lensing
fitting function. We will show how the width of the magnification distribution
depends strongly on the assumed cosmological parameters and modelling
assumptions. This lack of knowledge of the true form motivates an
investigation into the effect on a cosmological analysis of mischaracterising
the weak lensing uncertainty. The various weak lensing fitting functions used
in this work are summarised in Fig. 2.

Figure 2: The weak lensing uncertainty fitting functions used in this work.
Three have been created numerically in camb using different cosmologies:
Planck 2018 [58] and the Pantheon+ (P+) $\Lambda$CDM and $w$CDM maximum
posterior values [84]. Supernovae only constrain geometry, so Planck 2018
values for $n_{s}$ and $A_{s}$ are used. Also shown are results from CT21
using numerical simulations in T11, specifically for popIII / Q3 MBHBs
observed with LISA, including selection effects.

###  IV.1 Background

The magnification is given by

$\mu=\frac{1}{(1-\kappa)^{2}-\gamma^{2}}\,.$ (6)

The convergence $\kappa$ is the integrated surface mass density along the line
of sight, and $\gamma$ is the shear field. We are interested in the variance
of the magnification at a particular redshift $\sigma_{\mu}^{2}(z)$, as this
provides $\sigma_{\rm WL}(z)$ in Eq. (4). For the rest of this work, unless
otherwise stated, we will drop the small shear corrections which have little
impact on gravitational waves [85]. In that case, at a given redshift, the
variance in the convergence field is directly related to the variance in the
magnification [50],

$\sigma^{2}_{\mu}(z\,|\,\bm{\theta})\approx
4\sigma^{2}_{\kappa}(z\,|\,\bm{\theta})\,.$ (7)

The convergence field is cosmology dependent. The cosmological parameters
influencing the convergence field are $\bm{\theta}=\\{h,\Omega_{\rm
m},\Omega_{\rm b},\Sigma m_{\nu},n_{s},A_{s},\Omega_{K},w_{0},w_{a}\\}$.
$\Omega_{\rm b}$ is the baryon density fraction, $\Sigma m_{\nu}$ the sum of
neutrino masses, $A_{s}$ is the amplitude of scalar perturbations and $n_{s}$
the spectral tilt. Now, using Eq. (5), the variance in the convergence field
can be related to the weak lensing uncertainty we assume for a gravitational
wave source,

$\frac{\sigma_{\rm
WL}(z\,|\,\bm{\theta})}{d_{L}(z\,|\,\bm{\theta})}\approx\sigma_{\kappa}(z\,|\,\bm{\theta})\,.$
(8)

Knowing the dispersion, a magnification probability distribution,
$p(\mu,z\,|\,\bm{\theta})$, can be constructed. Poorly understood small-scale
physics significantly impacts the high-magnification tail of
$p(\mu,z\,|\,\bm{\theta})$. To avoid introducing modelling uncertainty, we
restrict our analysis to the weak lensing regime, with $\mu\leq 1.75$.
Therefore a simple shifted lognormal model as presented in Ref. [50] is
sufficient, but for comparison we also use $p(\mu,z)$ from the $N_{\rm
side}=4096$ and $N_{\rm side}=16384$ simulations described in Ref.
[86]333http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allsky_raytracing/
(henceforth T17).

To include the impact of selection effects, we must modify the magnification
probability distribution. Following Refs [87, 45],

$\mathcal{P}(\mu,z\,|\,\bm{\theta})=\mathcal{C}\,p(\mu,z\,|\,\bm{\theta})\frac{d\mathcal{N}(\mu,z)}{dz}\,,$
(9)

where

$\frac{d\mathcal{N}(\mu,z)}{dz}=\int_{\rho_{\rm
lim}/\sqrt{\mu}}^{\infty}N(z,\rho)\,d\rho\,.$ (10)

The SNR is given by $\rho$, $N(z,\rho)$ is the number density of sources as a
function of redshift and SNR, $\rho_{\rm lim}$ is the SNR threshold and
$\mathcal{C}$ is a normalisation constant. This modified magnification
distribution is drawn from randomly to find the magnification of sources. Fig.
16 of Ref. [20] shows that BSS SMBHBs have SNRs far above any possible SNR
threshold, hence have no selection effects. This is not the case for BSS BNSs.

The dispersion of the magnification distribution as a function of redshift can
be found using observations, simulations, or numerical calculation using a
halo model. These possibilities are explored in the following sections.

###  IV.2 Numerical calculation

Ref. [50] provide a tool to generate the dispersion in a redshift and
cosmology dependent convergence field, $\sigma_{\kappa}(z\,|\,\bm{\theta})$,
for varying cosmological parameters using the Boltzmann code
camb444https://camb.info/. We made minor modifications to their notebook;
using the HMcode-2020 nonlinear model [88], setting $k_{\rm
max}=100\,h\,\rm{Mpc}^{-1}$ as this is where HMcode switches from calling the
linear power spectrum to using a linear extrapolation, and including neutrinos
and dark energy to the fitting. The main uncertainty in this approach is the
modelling of non-linear scales. HMcode is known to not be reliable for $z>2$
[88], and calculating the convergence power spectrum requires a power-law
extrapolation of the Weyl potential at high $k-$modes which is not well
motivated. This calculation involves integrating the convergence power
spectrum over the multipole $\ell$ to find the total variance. This integral
convergences for $\ell_{\rm max}=10^{7}$ [50], and while we can not expect to
accurately predict the power spectrum at such high multipoles, omitting small
scales could lead to underpredicting the variance. Thus, knowledge of much
smaller scales than presently available is needed for accurate modelling of
$\sigma_{\kappa}(z\,|\,\bm{\theta})$.

Using this modified notebook, $\sigma_{\kappa}(z\,|\,\bm{\theta})$ is found
for different choices of cosmological parameters: Planck 2018 base
$\Lambda$CDM constraints from TT,TE,EE+lowE+lensing+BAO [58] (henceforth
labelled ‘Planck’), and the best fit parameters from the Pantheon+ analysis
assuming either $\Lambda$CDM or $w$CDM [84]. These choices highlight the
current tension in constraints of $H_{0}$. Supernovae only constrain geometry
so $n_{s}$ and $A_{s}$ are fixed to their Planck values when generating these
fitting functions. The weak lensing uncertainty for various cosmologies is
plotted in Fig. 2. The value of $h$, $\Omega_{\rm m}$, $w_{0}$ and $w_{a}$
used for each cosmology is given in Table 1.

Table 1:  We will evaluate the residual lensing bias when one of these three cosmologies is used to generate a weak lensing fitting function, while the truth is given by another. Cosmology | $h$ | $\Omega_{\rm m}$ | $w_{0}$ | $w_{a}$
---|---|---|---|---
Planck | $0.677$ | $0.311$ | $-1$ | $0$
Pantheon+ $\Lambda$CDM | $0.736$ | $0.344$ | $-1$ | $0$
Pantheon+ $w$CDM | $0.733$ | $0.403$ | $-0.93$ | $-0.1$

###  IV.3 Simulations

The present industry standard weak lensing fitting function [89, 47, 46, 48,
45, 49] was created in 2010 [89] using the magnification distribution of the
1998 work Ref. [90]. It has since been modified in CT21 using the ray-tracing
simulations presented in T11. The use of these simulations gives cause for
concern; they have limited resolution over a very small patch of the sky and
therefore could under-predict the width of the magnification distribution,
they are dark matter only simulations neglecting baryonic physics, and they
were created using outdated values of cosmological parameters ($h=0.705$,
$\Omega_{\rm m}=0.274$). Ref. [91] demonstrate how the inclusion of baryons in
hydrodynamic simulations modifies the high magnification tail ($\mu>3$), past
the weak lensing regime. This suggests that dark matter only simulations are
sufficient for evaluating the effects of weak lensing. But more recent works
using the IllustrisTNG [92] and MilleniumTNG [93] simulation suites suggest
that baryons, and in particular feedback effects, have a larger impact,
causing a narrower $p(\mu,z)$ as matter is redistributed from high to low
density regions. Furthermore, the magnification distribution is highly
dependent on cosmological parameters. Lower values of $h$ and $\Omega_{\rm m}$
predict a smaller $\sigma_{\kappa}(z\,|\,\bm{\theta})$. More recent
constraints on $\Omega_{\rm m}$ favour larger values [58, 84] than used in
T11, and there is the possibility $h$ favours the near-Universe measurement
[94, 84]. Therefore current weak lensing fitting functions in the literature
could be underestimating the dispersion of $p(\mu,z\,|\,\bm{\theta})$.

The function calibrated using T11, as presented in CT21, where the authors
consider lensing selection effects for the popIII and Q3 SMBHB seed models
separately, is plotted as a green solid and green dotted line in Fig. 2. The
fitting function, whose form was originally found in Ref. [89], is given by

$\frac{\sigma_{\rm
WL}}{d_{L}}=\frac{C}{2}\left(\frac{(1-(1+z)^{-\beta})}{\beta}\right)^{\alpha}\,,$
(11)

where $C=0.061$, $\beta=0.264$, $\alpha=1.89$ for the popIII seed model and
$C=0.096$, $\beta=0.62$, $\alpha=2.36$ for the Q3 model. Note the factor of
$1/2$, which is missing in Ref. [47]. This is because the original expression
is for the uncertainty in $\sigma({\rm
ln}\,d_{L}^{2})=2\sigma(d_{L})/d_{L}$555Incidentally, including this factor of
$1/2$ would improve results in numerous forecasting papers, e.g. Refs [47, 30,
48, 68, 95, 32]. [45, 57].

The difference between these curves and those generated using camb is
significant. This is due to the simulation’s resolution. As previously
mentioned, for the numerical calculation there is an integral over the
convergence power spectrum from $\ell_{\rm min}=1$ to $\ell_{\rm max}=10^{7}$.
The ray-tracing simulations of T11 use a $1\times 1\,$deg2 box with a grid
resolution of $3\,$kpc$/h$. The small box will fail to include super-clusters
or huge voids, both of which will broaden the magnification distribution. This
sets a high $\ell_{\rm min}$, while $\ell_{\rm max}$ increases with redshift
but is smaller than $10^{7}$. Hence the variance found from these simulations
will be smaller. Setting $\ell_{\rm min}$ and $\ell_{\rm max}$ in the
numerical calculation based on the simulation parameters leads to comparable
predictions for $\sigma_{\rm WL}(z)$ to Ref. [89] and CT21.

One could also consider using the ray-tracing simulations of T17 to find the
dispersion in the magnification distribution over finer redshift slices. They
are full-sky, providing access to lower multipoles, but have significantly
lower resolution in the high $\ell$ regime, and are available up to a maximum
redshift of $5.3$, compared to $z_{\rm max}=20$ in T11 (though there are only
7 redshift bins from $z=1$ to $z=20$). To demonstrate the importance of
simulation resolution, the magnification distribution from simulations of T17
for the same snapshot at $z=1.0334$, but with a different $N_{\rm side}$, can
be seen in Fig. 3. Both the convergence and shear fields are used to find
$p(\mu)$ from these simulations. There is a marked difference in the width of
$p(\mu)$, highlighting the importance of resolution (or equivalently,
$\ell_{\rm max}$) in estimating the weak lensing uncertainty. Also shown is
$p(\mu)$ from T11 at $z=1$, which shows close similarity to the $N_{\rm
side}=16384$ case. This suggests that even finer resolution full-sky
simulations could give a broader distribution than found in T11. Corresponding
weak lensing uncertainties found using the method of Ref. [89] can be seen in
Fig. 4. The weak lensing uncertainty found from the full-sky T17 $N_{\rm
side}=16384$ simulations follows the fitting function of Eq. (11) with
$C=0.076$, $\beta=0.344$ and $\alpha=2.25$ with no selection effects, and
$C=0.092$, $\beta=0.164$, $\alpha=2.09$ considering selection effects from a
ET+2CE network. The finer redshift slices and full-sky nature make this
fitting function more appropriate for the redshift regime where it is
applicable, $z\leq 5.3$. We suggest this fitting function should be used for
the lensing uncertainty on BNS mergers. However higher resolution full-sky
simulations will likely lead to a wider predicted magnification distribution.

Figure 3: The magnification probability distribution at $z=1.0334$ found from
the full-sky ray-tracing simulations of T17. Increasing the $N_{\rm side}$ of
the simulation from $4096$ to $16384$ leads to a broader distribution, and
therefore a larger assumed weak lensing uncertainty. Also overlaid is the
$z=1$ pdf from the T11 $1\times 1\,$deg2 ray-tracing simulations. Figure 4: A
comparison of weak lensing uncertainties from full-sky ray tracing simulations
with different resolutions (T17), high resolution ray-tracing simulations over
a $1\times 1\,$deg2 patch of sky (T11), and from the analytic pdf in Ref. [90]
used in Ref. [89] (H10). No selection effects have been included.

Weak lensing data products from other cosmological simulations exist [96, 97,
98, 99, 100, 101] and it would be interesting to see how choice of simulation
affects the derived weak lensing uncertainty. To our knowledge, the only other
example of a weak lensing uncertainty predicted from simulations is shown in
Ref. [46]. The authors use the ELUCID simulations [101] to estimate
$\sigma_{\rm WL}(z)$ for observations of binary black holes at $z\leq 1$. For
$z>1$ they also use results from T11, hence find general agreement with the
curves in CT21.

###  IV.4 Observations

Ideally, $p(\mu,z)$ would be inferred from observations. Through the weak
lensing of galaxies we measure reduced shear $g=\gamma/(1-\kappa)$, which is
either approximated as shear, or $\gamma$ is recovered through an iterative
scheme [102]. The convergence is then reconstructed from the shear. This is a
very noisy process, and measuring enough galaxies for accurate reconstruction
requires deep observations. Measuring $p(\mu,z)$ directly over the full sky in
fine redshift slices is not achievable in the near future. The best we can do
is calibrate simulations to lensing statistics from forthcoming large scale
galaxy surveys, such as the Vera C. Rubin Observatory [103] or Euclid [104].

##  V Results

###  V.1 Evaluating the residual lensing bias

For each combination of true and assumed cosmology, the steps in Section II
are performed to gain a large distribution of possible lensing biases for each
population. Using this distribution, we can characterise the impact of weak
lensing through the mean bias, and the dispersion of biases. Positively
magnified sources are more likely to be observed, this is a weak lensing
selection effect. When there are no selection effects $\langle\mu\rangle=1$
and the mean bias$\,\sim 0$. But with low SNR sources that can be shifted
inside the SNR limit through a weak lensing magnification,
$\langle\mu\rangle\neq 1$, and there is a bias in the recovered cosmological
parameters. This depends primarily on the source redshift distribution and the
magnification distribution. For a chosen population, the mean residual lensing
bias found by averaging over all statistical samples of a single catalogue
realisation shows little variation from realisation to realisation. Even
though the realisations have a range of source numbers, this will only impact
the parameter uncertainties and not the mean lensing bias. For the BNS
catalogues considered in this work, we report the mean absolute bias on $h$
and $\Omega_{\rm m}$ due to selection effects in Table 2. Results are shown
for the three cosmologies used to generate the numerical estimate of
$\sigma_{\rm WL}(z)$, and when magnifications are drawn directly from the
$N_{\rm side}=4096$ and $N_{\rm side}=16384$ simulations of T17. In the
$N_{\rm side}=4096$ case the mean bias is smaller as the limited resolution
leads to a smaller $\ell_{\rm max}$ than in the $N_{\rm side}=16384$ and
numerical cases. BSS SMBHBs need a high SNR to accurately localise and
associate with an EM counterpart so there are no magnification selection
effects on these sources, and their mean bias fluctuates around zero across
each catalogue.

Table 2:  The magnification selection effect on gravitational waves causes a bias $\overline{b}$ on recovered cosmological parameters. This depends on the source redshift distribution: results for the BNS catalogues used in this work are given here. The bias has been averaged over $300$ statistical realisations of $500$ catalogues. The bias also depends on the width of magnification distribution, which is affected by cosmological parameters and simulation resolution. $p(\mu)$ Cosmology | $\overline{b}(h)$ | $\overline{b}(\Omega_{\rm m})$
---|---|---
Planck Log Norm | $-0.0011$ | $0.0074$
Pantheon+ $\Lambda$CDM Log Norm | $-0.0015$ | $0.0098$
Pantheon+ $w$CDM Log Norm | $-0.0020$ | $0.0153$
T17 ($N_{\rm side}=4096$) | $-0.0007$ | $0.0044$
T17 ($N_{\rm side}=16384$) | $-0.0013$ | $0.0083$

Next, we go on to consider the impact of the assumption on $\sigma_{\rm
WL}(z)$. This will primarily affect the width of the bias distribution. From
Table 2 it can already be seen how mischaracterising the true cosmology, or
using low resolution simulations, when correcting for magnification effects
could lead to a residual lensing bias. For the SMBHB catalogues, the high
redshift sources lead to large magnifications and large absolute biases. But
the small source numbers mean large parameter uncertainties as the bias is
below the level of the variance. For BNSs the opposite is true: they have both
smaller magnifications and uncertainties. The absolute bias is small, but its
impact is larger. We quote results in two ways, first the probability of the
bias on $h$ ($\Omega_{\rm m}$) exceeding a value of $\pm 0.0025$ ($\pm 0.015$)
for BNSs and $\pm 0.02$ $(\pm 0.05$) for SMBHBs. This indicates the spread of
biases across catalogues. Then the probability of there being a discordance in
the parameter space greater than a certain number of standard deviations. This
takes into account the parameter uncertainties to show the impact of the
residual lensing bias on parameter constraints.

####  V.1.1 No added weak lensing uncertainty

First we will evaluate the impact of not including a weak lensing uncertainty
term in Eq. (4). If this term is not added, the probability of a lensing bias
is given in Table 3, where the true cosmology is Planck. It can be seen that
the probability of obtaining highly significant biases from LISA sources is
very large, due to the small uncertainties on the high redshift, high SNR
SMBHBs and the lack of a weak lensing uncertainty term. Even for the ground-
based observed BNSs the probability is high, demonstrating that weak lensing
is an important consideration for these sources.

Table 3:  Mean probability of finding a lensing bias when a lensing uncertainty term on the distance is not included. Both the probability of the bias on $h$ and $\Omega_{\rm m}$ exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D $h-\Omega_{\rm m}$ posterior, are shown. The true cosmology used is Planck 2018 [58]. For BNSs $x_{h}=0.0025$, $x_{\Omega_{\rm m}}=0.015$. For SMBHBs $x_{h}=0.02$, $x_{\Omega_{\rm m}}=0.05$.  Population | $P(b(h)\geq\pm x_{h})$ | $P(b(\Omega_{\rm m})\geq\pm x_{\Omega_{\rm m}})$ | $P(b\geq 1\sigma[3\sigma])$
---|---|---|---
popIII | $0.65$ | $0.72$ |  $0.90$ $[0.69]$
Q3 | $0.78$ | $0.77$ |  $0.98$ $[0.79]$
Q3nod | $0.56$ | $0.62$ |  $0.98$ $[0.81]$
BNS | 0.62 | 0.55 |  $0.79$ $[0.44]$

This scenario highlights the necessity of including $\sigma_{\rm WL}(z)$. Any
cosmological analysis using well measured binary mergers beyond
$z\\!\sim\\!0.5$ that does not could be significantly biased.

####  V.1.2 Correct uncertainty implemented

Now we test the situation where $\sigma_{\rm WL}(z)$ is generated using the
same distribution that magnifications are drawn from. There is perfect
knowledge of the true $p(\mu,z\,|\,\bm{\theta})$. This does not mean the
magnification can be removed from the GW, for this an external measurement of
$\mu$ is needed for each source. Magnifications are drawn from
$p(\mu,z\,|\,\bm{\theta})$ generated using Planck parameters and the added
uncertainty is the Planck curve in Fig. 2. Results are presented in Table 4.

In all populations, the probability of a large parameter shift remains high.
The absolute bias due to lensing is primarily dependent on the shape of the
magnification distribution, which is the same as it was in Table 3.

In the case of SMBHBs, for a bias $b\geq 1\sigma$ the probability is $\leq
0.05$, and for $b\geq 2\sigma$ the probability is $\leq 0.01$. The lensing
bias has been effectively removed by increasing uncertainties according to the
dispersion of $p(\mu,z\,|\,\bm{\theta})$.

For BNSs we are more likely to observe sources with a positive magnification.
There is an average negative$/$positive shift on $h/\Omega_{\rm m}$, also
shown in Ref. [42]. Due to this selection effect the probability of a residual
bias $b\geq 1\sigma\,(b\geq 2\sigma)$ is $0.18\,(0.01)$. The lensing selection
effect could in principle be corrected for using the method of CT21, outlined
in their Eq. (10). Following this we confirmed that, if source numbers are
large enough, binning them by redshift, finding $\langle d_{L}^{2}(z)\rangle$
in each bin and multiplying by $\langle\mu(z)\rangle$ for that bin removed the
lensing bias caused by selection effects. The key caveats of this method are
needing sufficient source numbers to obtain a reliable $\langle
d_{L}(z)\rangle$, and crucially, accurate knowledge of
$p(\mu,z\,|\,\bm{\theta})$. As will be demonstrated in the next section,
mischaracterising the magnification distribution could exacerbate the lensing
bias instead of mitigating it.

Table 4:  Mean probability of finding a lensing bias. Both the probability of the bias on $h$ and $\Omega_{\rm m}$ exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D $h-\Omega_{\rm m}$ posterior, are shown. The weak lensing uncertainty used is based on perfect knowledge of the shape of the magnification distribution, created using Planck 2018 parameters. For BNSs $x_{h}=0.0025$, $x_{\Omega_{\rm m}}=0.015$. For SMBHBs $x_{h}=0.02$, $x_{\Omega_{\rm m}}=0.05$. Population | $P(b(h)\geq\pm x_{h})$ | $P(b(\Omega_{\rm m})\geq\pm x_{\Omega_{\rm m}})$ | $P(b\geq 1\sigma[2\sigma])$
---|---|---|---
popIII | $0.45$ | $0.47$ | $0.03\,[0.00]$
Q3 | $0.54$ | $0.53$ | $0.03\,[0.00]$
Q3nod | $0.32$ | $0.40$ | $0.05\,[0.01]$
BNS | $0.21$ | $0.14$ | $0.18\,[0.01]$

####  V.1.3 Imperfect correction

Figure 5: The Lensed and Unlensed contours are constraints from the same
binary neutron star catalogue with $641$ sources. The true cosmology is the
the best fit Pantheon+ $\Lambda$CDM parameters. In both cases, the distance
uncertainty is increased according to a weak lensing uncertainty fitting
function generated using Planck 2018 parameters: this is the assumed
cosmology. The green “Unlensed” contour is the constraint on $h$ and
$\Omega_{\rm m}$ when all sources experience magnification $\mu=1$. The purple
“Lensed” contour is when each source experiences a magnification drawn from a
redshift dependent magnification distribution created using the true
cosmology. The parameters of this model give a broader magnification
distribution than the weak lensing uncertainty fitting function predicts,
hence the lensing bias is not completely mitigated; there is a total bias of
$2.0\sigma$, or an absolute bias of $b(h)=-0.007$ and $b(\Omega_{\rm
m})=0.031$.

In reality, we do not know the correct form of $p(\mu,z\,|\,\bm{\theta})$.
First we will consider the case when the uncertainty is underestimated, and
then when it is overestimated. If the weak lensing uncertainty is
underestimated, then the lensing bias demonstrated in Table 3 will not be
mitigated, as it has mostly been in Section V.1.2. There are two main reasons
the weak lensing uncertainty could be underestimated:

  * A.

Using too small values of $H_{0}$ or $\Omega_{\rm m}$ to predict $\sigma_{\rm
WL}(z)$.

  * B.

Incorrect modelling of the shape of the magnification distribution.

We investigate case A by assuming the true cosmology used to generate
$p(\mu,z\,|\,\bm{\theta})$ is either the Pantheon+ $\Lambda$CDM or Pantheon+
$w$CDM best fit parameters, and the weak lensing uncertainty is given by the
Planck curve in Fig. 2. Fig. 5 shows an example of the residual lensing bias
from ET+2CE observations of BNSs with an electromagnetic counterpart.
Probabilities of obtaining a lensing bias are presented in Table 5. A key
result is the importance of accurate modelling of $p(\mu,z\,|\,\bm{\theta})$
even for large numbers of BSS BNSs, where the lensing bias is greater than
that caused by selection effects alone. We also tested the bias when BNSs are
used to probe extensions to the base $\Lambda$CDM, including $k\Lambda$CDM
(curvature) and $w_{0}$CDM (dark energy) in Appendix C. For LISA SMBHBs,
despite the large probability of a significant parameter shift, the broad
parameter uncertainties make the probability of a large discordance caused by
lensing small.

Table 5:  Mean probability of finding a lensing bias. Both the probability of the bias on $h$ and $\Omega_{\rm m}$ exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D $h-\Omega_{\rm m}$ posterior, are shown. The weak lensing uncertainty is a fitting function based on the Planck 2018 cosmological parameters. The column headers give the true cosmology. For BNSs $x_{h}=0.0025$, $x_{\Omega_{\rm m}}=0.015$. For SMBHBs $x_{h}=0.02$, $x_{\Omega_{\rm m}}=0.05$. Population | Pantheon+ $\Lambda$CDM | Pantheon+ $w$CDM
---|---|---
$P(b(h)\geq\pm x_{h})$ | $P(b(\Omega_{\rm m})\geq\pm x_{\Omega_{\rm m}})$ | $P(b\geq 1\sigma[2\sigma])$ | $P(b(h)\geq\pm x_{h})$ | $P(b(\Omega_{\rm m})\geq\pm x_{\Omega_{\rm m}})$ | $P(b\geq 1\sigma[2\sigma])$
popIII | $0.55$ | $0.56$ | $0.07\,[0.02]$ | $0.62$ | $0.65$ | $0.15\,[0.04]$
Q3 | $0.66$ | $0.64$ | $0.10\,[0.02]$ | $0.72$ | $0.73$ | $0.18\,[0.05]$
Q3nod | $0.44$ | $0.53$ | $0.13\,[0.02]$ | $0.49$ | $0.64$ | $0.21\,[0.05]$
BNS | $0.32$ | $0.24$ | $0.29\,[0.02]$ | $0.50$ | $0.53$ | $0.49\,[0.10]$
Table 6:  Mean probability of finding a lensing bias. Both the probability of the bias on $h$ and $\Omega_{\rm m}$ exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D $h-\Omega_{\rm m}$ posterior, are shown. The true cosmology is given by the Planck parameters and the weak lensing uncertainty is given by CT21. For BNSs $x_{h}=0.0025$, $x_{\Omega_{\rm m}}=0.015$. For SMBHBs $x_{h}=0.02$, $x_{\Omega_{\rm m}}=0.05$. Population | Planck & CT21
---|---
$P(b(h)\geq\pm x_{h})$ | $P(b(\Omega_{\rm m})\geq\pm x_{\Omega_{\rm m}})$ | $P(b\geq 1\sigma[2\sigma])$
popIII | $0.53$ | $0.55$ | $0.32\,[0.09]$
Q3 | $0.64$ | $0.64$ | $0.38\,[0.11]$
Q3nod | $0.37$ | $0.44$ | $0.46\,[0.15]$
BNS | $0.32$ | $0.25$ | $0.54\,[0.19]$

To assess the impact of underestimating cosmological parameters, probabilities
in Table 5 are compared with those in Table 4, which describe when the weak
lensing uncertainty is correctly estimated. First consider if the true
cosmology is Pantheon+ $\Lambda$CDM, and assumed is Planck. For all
populations, the probability of the absolute bias on $h$ and $\Omega_{\rm m}$
of exceeding a certain value has increased. This is due to the wider
magnification distribution leading to more significant magnification effects.
For SMBHBs, the probability of obtaining results biased by more than $1\sigma$
is increased by a factor $\sim\\!2.5-3$, but still remains low due to the
large parameter uncertainties caused by few observations. For BNSs, the
probability of results biased by more than $1\sigma$ has increased by a factor
of $1.5$, and is close to $1/3$. If the true cosmology is in fact Pantheon+
$w$CDM, probabilities increase by an extra factor of $\sim\\!2$.

To test case B, we compare the fitting functions of CT21 with curves created
using camb. $\sigma_{\rm WL}(z\,|\,\bm{\theta})$ generated using the numerical
method of Section IV.2 and the cosmological parameters in T11 and T17 are very
similar to the Planck curve in Fig. 2 (orange dash-dot line). So using Planck
as the true cosmology to draw magnifications, and the fitting functions of
CT21 (green solid and dotted lines) as the lensing uncertainty acts as a model
comparison. The probability of obtaining a residual lensing bias in this case
is presented in Table 6 for each population. Once again, as the true cosmology
is the same as in Section V.1.2, we can compare bias probabilities with Table
4. In this case the uncertainty is too small to wash out the bias. For SMBHBs
the probability of results biased by $\geq 1\sigma$ increases by a factor of
$\sim\\!10$, and there is now a significant probability of results biased by
$\geq 2\sigma$. For BNSs, there is now greater than $50\%$ probability of
finding results biased by more than $1\sigma$ and for $b\geq 2\sigma$ the
probability has jumped from $1\%$ to $19\%$. Due to the impact of selection
effects on BNSs and their stronger constraining power, combining the SMBHB and
BNS catalogues has a small impact on the BNS results.

If the assumption on $\sigma_{\rm WL}(z\,|\,\bm{\theta})$ in fact
overestimates the dispersion in $p(\mu,z\,|\,\bm{\theta})$, then the
uncertainty on the cosmological parameters is excessively diluted. When the
true cosmology is Planck and the assumed cosmology is Pantheon+ $\Lambda$CDM
($w$CDM), the increased uncertainty at low redshift causes a degradation of
parameter constraints of $10\%\,(15\%)/8\%\,(12\%)$ for SMBHBs$/$BNSs.

###  V.2 Reducing the probability of a lensing bias

####  V.2.1 Setting requirements

In Section V.1 we calculated the residual lensing bias in several test cases,
demonstrating the bias can exceed the variance significantly if either
cosmological parameters are underestimated, or too-low resolution simulations
are underestimating the weak lensing uncertainty. We now set requirements on
the required precision’s to which the input cosmology should be known, and
discuss prospects for using state-of-the-art simulations to prevent the
residual lensing bias becoming a limiting factor for bright standard siren
analyses.

The dependence of the residual lensing bias on $H_{0}$ and $\Omega_{\rm m}$ to
the biased input parameters can be easily tested in our framework. For a fixed
true cosmology, the average lensing bias is constant as it only depends on the
magnification and source distributions. When we vary the input cosmology used
to calculate the weak lensing uncertainty, we are just scaling the average
uncertainty on recovered cosmological parameters (smaller assumed $\sigma_{\rm
WL}(z\,|\,\bm{\theta})$ leads to smaller parameter uncertainties). And so the
scaling of the average bias in units of $\sigma$ will simply depend on the
scaling of $\sigma$, the parameter uncertainties. The reason we have cast the
change of bias in this way is to remain agnostic of the true source
numbers$/$distribution and the true cosmology.

In Fig. 6 we show the change in the lensing bias (the colour of the heatmap)
as a function of the input bias on cosmological parameters (the $x-$ and
$y-$axes). The bias change is computed as a fractional change relative to the
lensing bias when the input bias is zero, this is the $(0,0)$ point on these
plots. When the input bias is zero, this can be considered a best-case
scenario. If the input bias is negative, the lensing bias is increased and the
bias change is greater than one. If the input bias is positive, then the bias
change is less than one, but the uncertainties on the recovered $H_{0}$ and
$\Omega_{\rm m}$ have been increased, weakening their constraints. We
confirmed that the results in Fig. 6 are general for different source numbers
and input cosmologies. The left plot shows that if the input $\Omega_{\rm m}$
is known to a precision of $0.02$ or better, the increase in bias is kept at
or below the $5\%$ level. The $H_{0}$ offset is less impactful; weak lensing
is less sensitive to the Hubble parameter [105], but $H_{0}$ should be known
in combination with $\Omega_{\rm m}$ to $\sim\\!2\,{\rm km}\,{\rm
s}^{-1}\,{\rm Mpc}^{-1}$. In the right plot we also consider the amplitude of
the linear matter power spectrum in spheres of radius $8\,{\rm
Mpc}^{-1}\,h^{-1}$, $\sigma_{8}$. We find that, if $\Omega_{\rm m}$ is to be
known to $0.02$ or better, then $\sigma_{8}$ should be known to a precision of
$0.07$ to keep the increase in bias below the $5\%$ level. The required
precision’s for $\sigma_{8}$ and $\Omega_{\rm m}$ are well within parameter
uncertainties from e.g. Planck, and are consistent with differences in these
parameters observed from different probes (e.g Refs [84, 5, 58]). While the
required precision on $H_{0}$ is also well within current constraints, it does
not span the $H_{0}$ tension.

It is reassuring that the precision to which we need to know $H_{0}$ and
$\Omega_{\rm m}$ is much larger than the expected constraints from standard
sirens on these parameters. But we only have one Universe and one set of
observations. As demonstrated in the previous section, the probability for
large lensing biases can still be high, even when we are close to being
correct to the input cosmology.

Figure 6: Heatmap displaying the lensing bias change due to incorrect
cosmological parameters when generating the weak lensing uncertainty. The $x-$
and $y-$axes show the difference between the parameter used to generate the
weak lensing uncertainty, and the true value of that parameter. The colour
shows the resulting change in the average residual lensing bias, relative to
the bias when the input cosmology is correct, $(0,0)$ on these plots. Also
overlaid is a contour showing the boundary where the change in the bias is at
the $5\%$ level.

Testing the required simulation resolution is more challenging, as we can not
know what the true magnification distribution should be in this case. We use
the full-sky simulations of Ref [86], using the expected weak lensing
uncertainty from each $N_{\rm side}$. We also downsample the $N_{\rm
side}=4096$ simulations to $N_{\rm side}=2048$ and $N_{\rm side}=1024$. We
find a linear scaling between the average bias and ${\rm log}_{2}(N_{\rm
side})$ (with and without the downsampled simulations), which is largely the
same, regardless of the true cosmology. This is due to a constant factor
difference between these weak lensing fitting functions at $z\geq 1$. For
every power of 2 increase in $N_{\rm side}$, the average bias drops by
$\sim\\!0.1\sigma$ where $\sigma$ is the recovered parameter uncertainty.

To fully investigate the effect of the simulation resolution on the
magnification distribution, and relate this to expected magnifications of real
sources, baryonic effects at small scales are required. The AbacusSummit full
simulation suite includes simulations with baryons, and weak lensing maps of
these are forthcoming [96]. Weak lensing maps for full hydro-dynamical
simulations with baryons have been calculated in Ref. [93], but for a single
resolution of $N_{\rm side}=12288$. In that work, the authors find that
baryons reduce the width of the magnification distribution. If this effect
could be probed down to smaller scales / greater resolution, this may indicate
a convergence of the width of the magnification distribution with resolution,
which would provide a best-guess weak lensing uncertainty.

A similar result could also be obtained using baryonic feedback effects in the
halo model used to numerically calculate $\sigma_{\rm WL}(z\,|\,\bm{\theta})$.
A feedback model, HMcode2020_feedback, is already included in camb [88], but
this model was not used in this work due to a divergence in the variance of
the power spectrum at small scales.

####  V.2.2 Other approaches

We consider two methods to address imperfect knowledge on the weak lensing
uncertainty. The first assumes we are confident in the modelling of
$\sigma_{\rm WL}(z)$ and know how it varies with cosmology, but are ignorant
of the true values of cosmological parameters. This could be achieved either
through improved modelling of small scales in numerical codes such as camb, or
a thorough investigation on how the width of the magnification distribution
varies across high resolution full-sky simulations using different assumed
cosmologies, such as is available with AbacusSummit. The generation of
$\sigma_{\rm WL}(z\,|\,\bm{\theta})$ is added as an ingredient to the MCMC
analysis with each new sampling of $h$, $\Omega_{\rm m}$ and any other
parameters being constrained. This does not guarantee results are unbiased,
but removes the sensitivity to an assumed cosmology.

Alternatively, we are not confident in our modelling of $\sigma_{\rm WL}(z)$,
and we do not know how it varies with cosmology. Here, we consider introducing
an extra parameter into the fitting which gives an amplitude modification to
the variance, similar to that given in Ref. [106].

$\sigma_{\rm WL}^{\prime}(z)=A\,\sigma_{\rm WL}(z)\,.$ (12)

A realistic prior on $A$ is needed. The fitting will try to increase the
variance as much as it can, as this will minimise $({\rm data}-{\rm
model})^{2}/{\rm variance}$. The prior on $A$ could be based on the difference
in predictions for the weak lensing uncertainty across simulations or
cosmological parameters. This is a simple choice, and neglects that we are
also uncertain on the variation of uncertainty with redshift. A second shape
changing parameter could be introduced, but this would also require a well
justified prior. This method attempts to find an acceptable bias-variance
trade-off, increasing uncertainties within a defined range as you fit the data
to a model. It also allows you to jointly constrain cosmology with the form of
the weak lensing uncertainty.

If the goal is reliability of parameter constraints, the simplest approach may
be to use a large enough weak lensing uncertainty to encompass all possible
curves, increasing parameter uncertainties by $\sim\\!15\%$ but strongly
reducing the probability of obtaining biased results.

##  VI Conclusions

Gravitational waves from binary mergers will be weakly lensed. To address the
unknown magnification of each source, a weak lensing uncertainty term must be
added to the observations. We aim to investigate the impact of
mischaracterising this term on a cosmological analysis.

Our key conclusions are:

  * •

Any cosmological analysis using binary mergers at $z\gtrsim 0.5$ not including
a weak lensing term has a high probability of being biased by more than
$1\sigma$.

  * •

Present weak lensing uncertainty fitting functions could be drastically
underestimating the width of the redshift dependent magnification
distribution. This is caused by the use of $1\times 1\,$deg2 simulations using
too small cosmological parameters. We provide a new form for the weak lensing
fitting formula in Eq. (11) using the full-sky ray-tracing simulations of T17
in Section IV.3, noting this reduces the probability of a residual lensing
bias for $z\leq 5.3$ bright standard sirens, such as merging binary neutron
stars.

  * •

If cosmological parameters are underestimated when generating $\sigma_{\rm
WL}(z)$, the probability of parameter constraints biased by more than
$1\sigma$ increases by a factor $1.5-3$ compared to the case when these
parameters are not underestimated.

  * •

If the magnification distribution is in fact much broader than predicted by
small sky-area ray-tracing simulations, but they are used to estimate the weak
lensing uncertainty, bias probabilities increase by a factor $3$ ($10$) for
BNSs (SMBHBs) compared to the case when the correct $p(\mu)$ is used. The
probability of a $1\sigma\,(2\sigma)$ bias in $H_{0}$ and $\Omega_{\rm m}$
found from BNS mergers is $54\%\,(19\%)$.

  * •

Considering the small absolute value of the bias for BNSs, with a spread of
$\Delta b(H_{0})\sim 0.25\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$,
magnification selection effects or incorrect assumptions on the true form of
$p(\mu\,|\,z,\bm{\theta})$ are unlikely to impact the efficacy of BNSs on
addressing the $H_{0}$ tension. However should more BSS SMBHBs be detected
than anticipated, the absolute bias having a spread $\Delta b(H_{0})\sim
5\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$ would be significant. An order of
magnitude more SMBHBs than considered in these catalogues (from $\sim\\!20$ to
$\sim\\!200$) would be needed for the bias to become a serious consideration.

  * •

Utilising the full potential of bright standard sirens requires accurate
modelling of $p(\mu\,|\,z,\bm{\theta})$ using $H_{0}$, $\Omega_{\rm m}$ and
$\sigma_{8}$ known to a precision better than $2\,$km s${}^{-1}\,$Mpc-1,
$0.02$ and $0.07$ respectively.

  * •

The most promising approach to estimate $\sigma_{\rm WL}(z)$ is using high-
resolution weak lensing maps, including baryons, over a range of resolutions.
A convergence of the width of $p(\mu,z)$ with resolution caused by baryonic
feedback effects would provide a reliable estimate for the weak lensing
uncertainty. This investigation will be possible in the near-future [96].

A similar, independent analysis was performed in Ref. [42]. There the true
cosmology is that used in the AbacusSummit simulations [96] with $h=0.6736$,
$\Omega_{\rm m}=0.3138$, similar to the Planck values used in this work, and
the weak lensing uncertainty is the fitting function of Ref. [47]. The
magnification distributions from AbacusSummit, as seen in Fig. 1 of Ref. [42],
are broader than those from the simulations of T11 from which the WL
uncertainty is derived, contributing to their derived bias. Ref. [42] include
larger lensing magnifications while we limit to $\mu\leq 1.75$ to avoid
modelling uncertainties. Therefore we could be under-representing the
probability of a bias.

With enough data, there is the possibility of constraining the parameters of
the weak lensing fitting function in Eq. (11) jointly with cosmology. We
briefly explored this in Eq. (12), where the amplitude of $\sigma_{\rm WL}(z)$
is treated as a free parameter. This is a hybrid approach, using the
information from observations to avoid fixing $\sigma_{\rm WL}(z)$ based only
on simulations, and was originally envisioned in Ref. [30]. Also, as pointed
out in Ref. [28], the relative abundance of magnified and de-magnified sources
reveals the non-Gaussianity of $p(\mu\,|\,z)$. We leave as future work a
detailed investigation into the prospects of jointly constraining cosmology
and the variance of the magnification distribution.

An interesting consideration is how these results impact delensing prospects.
Present delensing studies assume a form for the magnification distribution,
which if incorrect could lead to a biased delensing procedure. Future
delensing prospects, and the impact of the magnification distribution on
delensing, will be explored in a future work. Other future work includes
assessing the impact of weak lensing on numerous high redshift dark standard
sirens, and investigating the effect of the residual lensing bias for analyses
focusing on high redshift SMBHBs to study, for example, modified gravity [68,
95, 107].

###### Acknowledgements.

We thank Alberto Mangliagli for use of their SMBHB catalogues. CTM thanks
Ryuichi Takahashi, Jack Elvin-Poole and Pierre Burger for helpful discussion.
CTM is supported by a Science and Technology Facilities Council (STFC)
Studentship Grant. AT is supported by a STFC Consolidate Grant. M. H.
acknowledges the support of the Science and Technology Facilities Council
(Grant Ref ST/V005634/1). The python packages numpy, scipy, matplotlib,
getdist, emcee, gwfast and tensiometer have been used in this work. For the
purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

##  Appendix A Lensed distance uncertainty

Consider two cases, a source at $d_{L}$ with no magnification ($\mu=1$), and
the same source detected by the same detector but with a magnification $\mu$.

From Fisher formalism,

$\displaystyle\sigma_{d_{L}}$
$\displaystyle=\frac{1}{\sqrt{\Gamma_{d_{L}}}}\,,$ (13)
$\displaystyle\Gamma_{d_{L}}$
$\displaystyle=\left(\frac{\partial\tilde{h}(f)}{\partial
d_{L}}\mid\frac{\partial\tilde{h}(f)}{\partial d_{L}}\right)\,,$ (14)
$\displaystyle\Gamma_{d_{L}}$ $\displaystyle=4\int_{f_{\rm min}}^{f_{\rm
max}}\frac{1}{S_{h}(f)}\frac{\partial\tilde{h}^{*}(f)}{\partial
d_{L}}\frac{\partial\tilde{h}(f)}{\partial d_{L}}\,df\,.$ (15)

For simplicity assume

$\tilde{h}(f)=\frac{A}{d_{L}}M_{z}^{5/6}f^{-7/6}e^{i\bm{\Phi}(f)}\,.$ (16)

Then

$\displaystyle\Gamma_{d_{L}}$ $\displaystyle=\int_{f_{\rm min}}^{f_{\rm
max}}\frac{1}{S_{h}(f)}\frac{A^{2}}{d_{L}^{4}}M_{z}^{5/3}f^{-7/6}\,df\,,$ (17)
$\displaystyle\sigma_{d_{L}}$
$\displaystyle=\frac{d_{L}^{2}}{A}\left(M_{z}^{5/3}\int_{f_{\rm min}}^{f_{\rm
max}}\frac{1}{S_{h}(f)}f^{-7/6}\,df\right)^{-1/2}\,.$ (18)

Hence, the uncertainty on $d_{L}$

$\sigma_{d_{L}}\propto d_{L}^{2}\,.$ (19)

Therefore, using

$d_{L}^{\prime}=\frac{d_{L}}{\sqrt{\mu}}\,,$ (20)

where ′ corresponds to lensed quantities,

$\displaystyle\frac{\sigma_{d_{L}}^{\prime}}{(d_{L}^{\prime})^{2}}$
$\displaystyle=\frac{\sigma_{d_{L}}}{d_{L}^{2}}\,,$ (21)
$\displaystyle\sigma_{d_{L}}^{\prime}$
$\displaystyle=\frac{\sigma_{d_{L}}\left(d_{L}/\sqrt{\mu}\right)^{2}}{d_{L}^{2}}\,,$
(22) $\displaystyle\sigma_{d_{L}}^{\prime}$
$\displaystyle=\frac{\sigma_{d_{L}}}{\mu}\,.$ (23)

This relation has been verified by passing a catalogue, both before and after
lensing, to the GW Fisher parameter uncertainty code GWFish [108].

##  Appendix B Other ground-based detector networks

###  B.1 ET and CE

If there is only one CE like detector in the network, the sky localisation is
degraded. We use a single catalogue of 500 sources to test this case, finding
a factor 2 increase in the probability of exceeding $1\sigma$ when the Planck
curve in Fig. 2 is used for $\sigma_{\rm WL}(z)$ and the true cosmology is
given by the Pantheon+ $\Lambda$CDM constraints. There is a factor 5 increase
in bias probabilities when CT21 is used for $\sigma_{\rm WL}(z)$ and Planck is
the true cosmology. The residual lensing bias is less impactful due to the
degraded measurement uncertainties, but still significant.

###  B.2 ET and Voyager

Here we consider a combination of Einstein Telescope and a next generation LVK
detector [12], to aid sky localisation. Due to the EM counterpart selection
effects, such as a sky localisation $<20\,$deg2 and a distance uncertainty
$<30\%$, the redshift of BSS BNS from this detector combination are $\leq
0.5$. Selection effects now cause a small bias of $b(h)=-0.0001$ and
$b(\Omega_{\rm m})=0.0013$.

If we relax the assumption on the sky localisation to $<200\,$deg2, the bias
due to selection effects increases by a factor of 5 due to the inclusion of
higher redshift sources. But because the uncertainty on $h$ and $\Omega_{\rm
m}$ when being jointly constrained by BSS observed in this GW detector network
is large, the probability of the residual lensing bias exceeding $1\sigma$ is
$0.04$ even when no weak lensing uncertainty term is included. The measurement
uncertainty dominates.

We can conclude that the impact of residual lensing biases are most important
for high-redshift binaries observed by LISA, or BNS mergers observed through a
ground-based network of two or more 3G detectors.

##  Appendix C Residual lensing bias in extensions to $\Lambda$CDM

###  C.1 Curvature

We consider the scenario where the true cosmology is the Pantheon+
$\Lambda$CDM best fit values, and the weak lensing uncertainty is derived
using Planck 2018 parameters. The source population are BNSs observed by the
3G detector network. We find the mean parameter biases caused by magnification
selection effects are $\overline{b}(h)=0.0004$, $\overline{b}(\Omega_{\rm
m})=0.0161$, $\overline{b}(\Omega_{k})=-0.03$. We find the probabilities
$P(b(h)\geq\pm 0.0025)=0.86$, $P(b(\Omega_{\rm m})\geq\pm 0.015)=0.60$,
$P(b(\Omega_{k})\geq\pm 0.1)=0.52$. It can be seen that, for $\Omega_{k}$
particularly, the probability of a residual lensing bias that is both large in
terms of the parameter, and large in terms of the parameter uncertainties, is
high. The probability of a discordance in the data set $\geq 1\sigma\,(\geq
2\sigma)$ is $0.19\,(0.02)$

Figure 7: Similar to Fig. 5 but for the $k\Lambda$CDM cosmological model,
where we are also constraining curvature. In this case lensing leads to a
$2.7\sigma$ bias in this parameter space, and a $>2\sigma$ preference for
negative curvature.

###  C.2 Dark Energy

We use the same scenario as in the case of curvature—Pantheon+ $\Lambda$CDM
for the true cosmology and Planck to generate the assumed weak lensing
uncertainty. We find the mean parameter biases caused by magnification
selection effects are $\overline{b}(h)=0.0007$, $\overline{b}(\Omega_{\rm
m})=0.0236$, $\overline{b}(w_{0})=-0.0574$. We find the probabilities
$P(b(h)\geq\pm 0.0025)=0.72$, $P(b(\Omega_{\rm m})\geq\pm 0.015)=0.77$,
$P(b(w_{0})\geq\pm 0.1)=0.54$. Once again, the parameter extension to
$\Lambda$CDM causes significant biases. The probability of a discordance in
the data set $\geq 1\sigma\,(\geq 2\sigma)$ is $0.54\,(0.01)$

Figure 8: Similar to Fig. 5 but for the $w_{0}$CDM cosmological model, where
we are also constraining a constant dark energy equation of state. In this
case lensing leads to a $1.2\sigma$ bias in this parameter space, and a
$>2\sigma$ preference for a phantom dark energy.

## References

  * Bartelmann and Schneider [2001] M. Bartelmann and P. Schneider,  _Weak gravitational lensing_ , Phys. Rep. 340, 291–472 (2001).
  * Mandelbaum [2018] R. Mandelbaum,  _Weak lensing for precision cosmology_ , Annu. Rev. Astron. Astrophys. 56, 393–433 (2018).
  * Heymans  _et al._ [2021] C. Heymans, T. Tröster, M. Asgari, C. Blake, H. Hildebrandt, B. Joachimi, K. Kuijken, C.-A. Lin, _et al._ ,  _KiDS-1000 cosmology: Multi-probe weak gravitational lensing and spectroscopic galaxy clustering constraints_ , Astron. Astrophys. 646, A140 (2021).
  * Dalal  _et al._ [2023] R. Dalal, X. Li, A. Nicola, J. Zuntz, M. A. Strauss, S. Sugiyama, T. Zhang, M. M. Rau, _et al._ ,  _Hyper Suprime-Cam Year 3 results: Cosmology from cosmic shear power spectra_ , Phys. Rev. D 108, 123519 (2023).
  * Abbott  _et al._ [2022] T. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves, A. Amon, F. Andrade-Oliveira, J. Annis, _et al._ ,  _Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing_ , Phys. Rev. D 105, 023520 (2022).
  * Marković [1993] D. Marković,  _Possibility of determining cosmological parameters from measurements of gravitational waves emitted by coalescing, compact binaries_ , Phys. Rev. D 48, 4738–4756 (1993).
  * Wang  _et al._ [1996] Y. Wang, A. Stebbins, and E. L. Turner,  _Gravitational lensing of gravitational waves from merging neutron star binaries_ , Phys. Rev. Lett. 77, 2875–2878 (1996).
  * Takahashi and Nakamura [2003] R. Takahashi and T. Nakamura,  _Wave effects in the gravitational lensing of gravitational waves from chirping binaries_ , Astrophys. J. 595, 1039–1051 (2003).
  * Takahashi  _et al._ [2005] R. Takahashi, T. Suyama, and S. Michikoshi,  _Scattering of gravitational waves by the weak gravitational fields of lens objects_ , Astron. Astrophys. 438, L5–L8 (2005).
  * Holz and Hughes [2005] D. E. Holz and S. A. Hughes,  _Using gravitational-wave standard sirens_ , Astrophy. J. 629, 15–22 (2005).
  * The LIGO Scientific collaboration  _et al._ [2016] The LIGO Scientific collaboration, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, _et al._ ,  _Observation of Gravitational Waves from a Binary Black Hole Merger_ , Phys. Rev. Lett. 116, 061102 (2016).
  * The LIGO Scientific collaboration [2019] The LIGO Scientific collaboration, Gravitational wave astronomy with LIGO and similar detectors in the next decade (2019), arXiv:1904.03187 .
  * Aso  _et al._ [2013] Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto,  _Interferometer design of the KAGRA gravitational wave detector_ , Phys. Rev. D 88, 043007 (2013).
  * Unnikrishnan [2013] C. S. Unnikrishnan,  _IndIGO and LIGO-India: Scope and Plans for Gravitational Wave Research and Precision Metrology in India_ , Int. J. Mod. Phys. A 22, 1341010 (2013).
  * Abbott  _et al._ [2021] R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. X. Adhikari, _et al._ ,  _Search for lensing signatures in the gravitational-wave observations from the first half of LIGO–Virgo’s third observing run_ , Astrophys. J. 923, 14 (2021).
  * Liao  _et al._ [2017] K. Liao, X.-L. Fan, X. Ding, M. Biesiada, and Z.-H. Zhu,  _Precision cosmology from future lensed gravitational wave and electromagnetic signals_ , Nat. Commun. 8, 1148 (2017).
  * The eLISA Consortium  _et al._ [2013] The eLISA Consortium, P. A. Seoane, S. Aoudia, H. Audley, G. Auger, S. Babak, J. Baker, E. Barausse, S. Barke, _et al._ , The Gravitational Universe (2013), arXiv:1305.5720 .
  * Klein  _et al._ [2016] A. Klein, E. Barausse, A. Sesana, A. Petiteau, E. Berti, S. Babak, J. Gair, S. Aoudia, _et al._ ,  _Science with the space-based interferometer eLISA: Supermassive black hole binaries_ , Phys. Rev. D 93, 024003 (2016).
  * Bonetti  _et al._ [2019] M. Bonetti, A. Sesana, F. Haardt, E. Barausse, and M. Colpi,  _Post-Newtonian evolution of massive black hole triplets in galactic nuclei – IV. Implications for LISA_ , Mon. Not. R. Astron. Soc. 486, 4044–4060 (2019).
  * Mangiagli  _et al._ [2022] A. Mangiagli, C. Caprini, M. Volonteri, S. Marsat, S. Vergani, N. Tamanini, and H. Inchauspé ,  _Massive black hole binaries in LISA: Multimessenger prospects and electromagnetic counterparts_ , Phys. Rev. D 106, 103017 (2022).
  * Maggiore  _et al._ [2020] M. Maggiore, C. V. D. Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, M. Branchesi, _et al._ ,  _Science case for the Einstein Telescope_ , J. Cosmol. Astropart. 2020, 050–050 (2020).
  * Reitze  _et al._ [2019] D. Reitze, R. X. Adhikari, S. Ballmer, B. Barish, L. Barsotti, G. Billingsley, D. A. Brown, Y. Chen, _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, in  _Bull. Am. Astron. Soc._, Vol. 51 (2019) p. 35.
  * Belgacem  _et al._ [2019a] E. Belgacem, Y. Dirian, S. Foffa, E. J. Howell, M. Maggiore, and T. Regimbau,  _Cosmology and dark energy from joint gravitational wave-GRB observations_ , J. Cosmol. Astropart. 2019, 015–015 (2019a).
  * Hilbert  _et al._ [2011] S. Hilbert, J. R. Gair, and L. J. King,  _Reducing distance errors for standard candles and standard sirens with weak-lensing shear and flexion maps_ , Mon. Not. R. Astron. Soc. 412, 1023–1037 (2011).
  * Shapiro  _et al._ [2010] C. Shapiro, D. J. Bacon, M. Hendry, and B. Hoyle,  _Delensing gravitational wave standard sirens with shear and flexion maps_ , Mon. Not. R. Astron. Soc. 404, 858–866 (2010).
  * Wu  _et al._ [2023] Z.-F. Wu, L. W. L. Chan, M. Hendry, and O. A. Hannuksela,  _Reducing the impact of weak-lensing errors on gravitational-wave standard sirens_ , Mon. Not. R. Astron. Soc. 522, 4059–4077 (2023).
  * Schutz [1986] B. F. Schutz,  _Determining the Hubble constant from gravitational wave observations_ , Nature 323, 310–311 (1986).
  * Cutler and Holz [2009] C. Cutler and D. E. Holz,  _Ultrahigh precision cosmology from gravitational waves_ , Phys. Rev. D 80, 104009 (2009).
  * Camera and Nishizawa [2013] S. Camera and A. Nishizawa,  _Beyond Concordance Cosmology with Magnification of Gravitational-Wave Standard Sirens_ , Phys. Rev. Lett. 110, 151103 (2013).
  * Congedo and Taylor [2019] G. Congedo and A. Taylor,  _Joint cosmological inference of standard sirens and gravitational wave weak lensing_ , Phys. Rev. D 99, 083526 (2019).
  * Mukherjee  _et al._ [2020] S. Mukherjee, B. D. Wandelt, and J. Silk,  _Multimessenger tests of gravity with weakly lensed gravitational waves_ , Phys. Rev. D 101, 103509 (2020).
  * Mpetha  _et al._ [2023] C. T. Mpetha, G. Congedo, and A. Taylor,  _Future prospects on testing extensions to $\Lambda$CDM through the weak lensing of gravitational waves_, Phys. Rev. D 107, 103518 (2023).
  * Balaudo  _et al._ [2023] A. Balaudo, A. Garoffolo, M. Martinelli, S. Mukherjee, and A. Silvestri,  _Prospects of testing late-time cosmology with weak lensing of gravitational waves and galaxy surveys_ , J. Cosmol. Astropart. 2023, 050 (2023).
  * Scelfo  _et al._ [2023] G. Scelfo, M. Berti, A. Silvestri, and M. Viel,  _Testing gravity with gravitational waves × electromagnetic probes cross-correlations_ , J. Cosmol. Astropart. Phys. 2023 (2), 010.
  * Kaiser and Peacock [2015] N. Kaiser and J. A. Peacock,  _On the bias of the distance–redshift relation from gravitational lensing_ , Mon. Not. R. Astron. Soc. 455, 4518–4547 (2015).
  * Sarkar  _et al._ [2008] D. Sarkar, A. Amblard, D. E. Holz, and A. Cooray,  _Lensing and supernovae: Quantifying the bias on the dark energy equation of state_ , Astrophys. J. 678, 1–5 (2008).
  * Smith  _et al._ [2013] M. Smith, D. J. Bacon, R. C. Nichol, H. Campbell, C. Clarkson, R. Maartens, C. B. D'Andrea, B. A. Bassett, _et al._ ,  _The effect of weak lensing on distance estimates from Supernovae_ , Astrophys. J. 780, 24 (2013).
  * Shang and Haiman [2010] C. Shang and Z. Haiman,  _Cosmology with standard sirens: the importance of the shape of the lensing magnification distribution_ , Mon. Not. R. Astron. Soc. 411, 9–22 (2010).
  * Fleury  _et al._ [2017] P. Fleury, C. Clarkson, and R. Maartens,  _How does the cosmic large-scale structure bias the Hubble diagram?_ , J. Cosmol. Astropart. 2017, 062–062 (2017).
  * Adamek  _et al._ [2019] J. Adamek, C. Clarkson, L. Coates, R. Durrer, and M. Kunz,  _Bias and scatter in the Hubble diagram from cosmological large-scale structure_ , Phys. Rev. D 100, 021301 (2019).
  * Yoo  _et al._ [2007] C.-M. Yoo, K.-I. Nakao, H. Kozaki, and R. Takahashi,  _Lensing effects on gravitational waves in a clumpy universe: Effects of inhomogeneity on the distance-redshift relation_ , Astrophys. J. 655, 691–703 (2007).
  * Canevarolo and Chisari [2023] S. Canevarolo and N. E. Chisari, Lensing bias on cosmological parameters from bright standard sirens (2023), arXiv:2310.12764 .
  * Chevallier and Polarski [2001] M. Chevallier and D. Polarski,  _Accelerating Universes with scaling dark matter_ , Int. J. Mod. Phys. D 10, 213–223 (2001).
  * Linder [2003] E. V. Linder,  _Exploring the Expansion History of the Universe_ , Phys. Rev. Lett. 90, 091301 (2003).
  * Cusin and Tamanini [2021] G. Cusin and N. Tamanini,  _Characterization of lensing selection effects for LISA massive black hole binary mergers_ , Mon. Not. R. Astron. Soc. 504, 3610–3618 (2021).
  * Shan  _et al._ [2021] X. Shan, C. Wei, and B. Hu,  _Lensing magnification: gravitational waves from coalescing stellar-mass binary black holes_ , Mon. Not. R. Astron. Soc. 508, 1253–1261 (2021).
  * Tamanini  _et al._ [2016] N. Tamanini, C. Caprini, E. Barausse, A. Sesana, A. Klein, and A. Petiteau,  _Science with the space-based interferometer eLISA. III: probing the expansion of the universe using gravitational wave standard sirens_ , J. Cosmol. Astropart. 2016, 002 (2016).
  * Speri  _et al._ [2021] L. Speri, N. Tamanini, R. R. Caldwell, J. R. Gair, and B. Wang,  _Testing the quasar Hubble diagram with LISA standard sirens_ , Phys. Rev. D 103, 083526 (2021).
  * Dai  _et al._ [2017] L. Dai, T. Venumadhav, and K. Sigurdson,  _Effect of lensing magnification on the apparent distribution of black hole mergers_ , Phys. Rev. D 95, 044011 (2017).
  * Martinelli  _et al._ [2022] M. Martinelli, F. Scarcella, N. B. Hogg, B. J. Kavanagh, D. Gaggero, and P. Fleury,  _Dancing in the dark: detecting a population of distant primordial black holes_ , J. Cosmol. Astropart. 2022, 006 (2022).
  * Verde  _et al._ [2019] L. Verde, T. Treu, and A. G. Riess,  _Tensions between the early and late universe_ , Nat. Astron. 3, 891–895 (2019).
  * Riess  _et al._ [2023] A. G. Riess, G. S. Anand, W. Yuan, S. Casertano, A. Dolphin, L. M. Macri, L. Breuval, D. Scolnic, _et al._ ,  _Crowded No More: The Accuracy of the Hubble Constant Tested with High-resolution Observations of Cepheids by JWST_ , Astrophys. J. Lett. 956, L18 (2023).
  * Takahashi  _et al._ [2011] R. Takahashi, M. Oguri, M. Sato, and T. Hamana,  _Probability distribution functions of cosmological lensing: convergence, shear and magnification_ , Astrophys. J. 742, 15 (2011).
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman,  _emcee : The MCMC hammer_, Publ. Astron. Soc. Pac. 125, 306–312 (2013).
  * Raveri  _et al._ [2020] M. Raveri, G. Zacharegkas, and W. Hu,  _Quantifying concordance of correlated cosmological data sets_ , Phys. Rev. D 101, 103527 (2020).
  * Raveri and Doux [2021] M. Raveri and C. Doux,  _Non-gaussian estimates of tensions in cosmological parameters_ , Phys. Rev. D 104, 043504 (2021).
  * Muttoni  _et al._ [2023] N. Muttoni, D. Laghi, N. Tamanini, S. Marsat, and D. Izquierdo-Villalba,  _Dark siren cosmology with binary black holes in the era of third-generation gravitational wave detectors_ , Phys. Rev. D 108, 043543 (2023).
  * The Planck Collaboration  _et al._ [2020] The Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, _et al._ ,  _Planck 2018 results_ , Astron. Astrophys. 641, A6 (2020).
  * Mangiagli  _et al._ [2023] A. Mangiagli, C. Caprini, S. Marsat, L. Speri, R. R. Caldwell, and N. Tamanini, Massive black hole binaries in LISA: constraining cosmological parameters at high redshifts (2023), arXiv:2312.04632 .
  * Madau and Rees [2001] P. Madau and M. J. Rees,  _Massive black holes as Population III remnants_ , Astrophys. J. 551, L27–L30 (2001).
  * Begelman  _et al._ [2006] M. C. Begelman, M. Volonteri, and M. J. Rees,  _Formation of supermassive black holes by direct collapse in pre-galactic haloes_ , Mon. Not. R. Astron. Soc. 370, 289–298 (2006).
  * Volonteri  _et al._ [2007] M. Volonteri, G. Lodato, and P. Natarajan,  _The evolution of massive black hole seeds_ , Mon. Not. R. Astron. Soc. 383, 1079–1088 (2007).
  * Fragione [2021] G. Fragione,  _Black-hole–neutron-star mergers are unlikely multimessenger sources_ , Astrophys. J. Lett. 923, L2 (2021).
  * Iacovelli  _et al._ [2022a] F. Iacovelli, M. Mancarella, S. Foffa, and M. Maggiore,  _Forecasting the detection capabilities of third-generation gravitational-wave detectors using GWFAST_ , Astrophys. J. 941, 208 (2022a).
  * Iacovelli  _et al._ [2022b] F. Iacovelli, M. Mancarella, S. Foffa, and M. Maggiore,  _GWFAST: A fisher information matrix python code for third-generation gravitational-wave detectors_ , Astrophys. J., Suppl. Ser. 263, 2 (2022b).
  * Borhanian and Sathyaprakash [2022] S. Borhanian and B. S. Sathyaprakash, Listening to the universe with next generation ground-based gravitational-wave detectors (2022), arXiv:2202.11048 .
  * Stratta  _et al._ [2018] G. Stratta, R. Ciolfi, L. Amati, E. Bozzo, G. Ghirlanda, E. Maiorano, L. Nicastro, A. Rossi, _et al._ ,  _THESEUS: A key space mission concept for multi-messenger astrophysics_ , Adv. Space Res. 62, 662–682 (2018).
  * Yang [2021] T. Yang,  _Gravitational-wave detector networks: standard sirens on cosmology and modified gravity theory_ , J. Cosmol. Astropart. 2021, 044 (2021).
  * Ronchini  _et al._ [2022] S. Ronchini, M. Branchesi, G. Oganesyan, B. Banerjee, U. Dupletsa, G. Ghirlanda, J. Harms, M. Mapelli, _et al._ ,  _Perspectives for multimessenger astronomy with the next generation of gravitational-wave detectors and high-energy satellites_ , Astron. Astrophys. 665, A97 (2022).
  * Alfradique  _et al._ [2022] V. Alfradique, M. Quartin, L. Amendola, T. Castro, and A. Toubiana,  _The lure of sirens: joint distance and velocity measurements with third-generation detectors_ , Mon. Not. R. Astron. Soc. 517, 5449–5462 (2022).
  * Chen  _et al._ [2021a] H.-Y. Chen, P. S. Cowperthwaite, B. D. Metzger, and E. Berger,  _A program for multimessenger standard siren cosmology in the era of LIGO A+, Rubin Observatory, and beyond_ , Astrophys. J. 908, L4 (2021a).
  * Gardner  _et al._ [2023] J. W. Gardner, L. Sun, S. Borhanian, P. D. Lasky, E. Thrane, D. E. McClelland, and B. J. J. Slagmolen,  _Multimessenger astronomy with a Southern-hemisphere gravitational-wave observatory_ , Phys. Rev. D 108, 123026 (2023).
  * Sathyaprakash  _et al._ [2012] B. Sathyaprakash, M. Abernathy, F. Acernese, P. Ajith, B. Allen, P. Amaro-Seoane, N. Andersson, S. Aoudia, _et al._ ,  _Scientific objectives of Einstein Telescope_ , Class. Quantum Gravity 29, 124013 (2012).
  * Chen  _et al._ [2021b] H.-Y. Chen, P. S. Cowperthwaite, B. D. Metzger, and E. Berger,  _A program for multimessenger standard siren cosmology in the era of LIGO A $+$, Rubin Observatory, and beyond_, Astrophys. J. 908, L4 (2021b).
  * Andreoni  _et al._ [2022] I. Andreoni, R. Margutti, O. S. Salafia, B. Parazin, V. A. Villar, M. W. Coughlin, P. Yoachim, K. Mortensen, _et al._ ,  _Target-of-opportunity observations of gravitational-wave events with Vera C. Rubin Observatory_ , Astrophys. J., Suppl. Ser. 260, 18 (2022).
  * Zhao  _et al._ [2011] W. Zhao, C. V. D. Broeck, D. Baskaran, and T. G. F. Li,  _Determination of dark energy by the Einstein Telescope: Comparing with CMB, BAO, and SNIa observations_ , Phys. Rev. D 83, 023005 (2011).
  * Chen  _et al._ [2019] H.-Y. Chen, S. Vitale, and R. Narayan,  _Viewing angle of binary neutron star mergers_ , Phys. Rev. X 9, 031028 (2019).
  * Nakar and Piran [2021] E. Nakar and T. Piran,  _Afterglow constraints on the viewing angle of binary neutron star mergers and determination of the hubble constant_ , Astrophys. J. 909, 114 (2021).
  * de Souza and Sturani [2023] J. M. S. de Souza and R. Sturani,  _Luminosity distance uncertainties from gravitational wave detections of binary neutron stars by third generation observatories_ , Phys. Rev. D 108, 043027 (2023).
  * Hotokezaka  _et al._ [2019] K. Hotokezaka, E. Nakar, O. Gottlieb, S. Nissanke, K. Masuda, G. Hallinan, K. P. Mooley, and A. T. Deller,  _A Hubble constant measurement from superluminal motion of the jet in GW170817_ , Nat. Astron. 3, 940–944 (2019).
  * Shah  _et al._ [2013] S. Shah, G. Nelemans, and M. van der Sluys,  _Using electromagnetic observations to aid gravitational-wave parameter estimation of compact binaries observed with LISA_ , Astron. Astrophys. 553, A82 (2013).
  * Hahn  _et al._ [2023] C. Hahn, M. J. Wilson, O. Ruiz-Macias, S. Cole, D. H. Weinberg, J. Moustakas, A. Kremin, J. L. Tinker, _et al._ ,  _The desi bright galaxy survey: Final target selection, design, and validation_ , Astron. J. 165, 253 (2023).
  * Schlegel  _et al._ [2022] D. J. Schlegel, J. A. Kollmeier, G. Aldering, S. Bailey, C. Baltay, C. Bebek, S. BenZvi, R. Besuner, _et al._ , The MegaMapper: A Stage-5 Spectroscopic Instrument Concept for the Study of Inflation and Dark Energy (2022), arXiv:2209.04322 .
  * Brout  _et al._ [2022] D. Brout, D. Scolnic, B. Popovic, A. G. Riess, A. Carr, J. Zuntz, R. Kessler, T. M. Davis, _et al._ ,  _The Pantheon $+$ analysis: Cosmological constraints_, Astrophys. J. 938, 110 (2022).
  * Takahashi [2006] R. Takahashi,  _Amplitude and phase fluctuations for gravitational waves propagating through inhomogeneous mass distribution in the universe_ , Astrophys. J. 644, 80–85 (2006).
  * Takahashi  _et al._ [2017] R. Takahashi, T. Hamana, M. Shirasaki, T. Namikawa, T. Nishimichi, K. Osato, and K. Shiroyama,  _Full-sky gravitational lensing simulation for large-area galaxy surveys and cosmic microwave background experiments_ , Astrophys. J. 850, 24 (2017).
  * Cusin  _et al._ [2020] G. Cusin, R. Durrer, and I. Dvorkin, Strong and weak lensing of gravitational waves: a semi-analytical approach (2020), arXiv:1912.11916 .
  * Mead  _et al._ [2021] A. J. Mead, S. Brieden, T. Tröster, and C. Heymans,  _HMcode-2020: improved modelling of non-linear cosmological power spectra with baryonic feedback_ , Mon. Not. R. Astron. Soc. 502, 1401–1422 (2021).
  * Hirata  _et al._ [2010] C. M. Hirata, D. E. Holz, and C. Cutler,  _Reducing the weak lensing noise for the gravitational wave hubble diagram using the non-gaussianity of the magnification distribution_ , Phys. Rev. D 81, 124046 (2010).
  * Holz and Wald [1998] D. E. Holz and R. M. Wald,  _New method for determining cumulative gravitational lensing effects in inhomogeneous universes_ , Phys. Rev. D 58, 063501 (1998).
  * Castro  _et al._ [2018] T. Castro, M. Quartin, C. Giocoli, S. Borgani, and K. Dolag,  _The effect of baryons in the cosmological lensing PDFs_ , Mon. Not. R. Astron. Soc. 478, 1305–1325 (2018).
  * Osato  _et al._ [2021] K. Osato, J. Liu, and Z. Haiman,  _$\kappa$ TNG: effect of baryonic processes on weak lensing with IllustrisTNG simulations_, Mon. Not. R. Astron. Soc. 502, 5593–5602 (2021).
  * Ferlito  _et al._ [2023] F. Ferlito, V. Springel, C. T. Davies, C. Hernández-Aguayo, R. Pakmor, M. Barrera, S. D. M. White, A. M. Delgado, _et al._ ,  _The MillenniumTNG project: the impact of baryons and massive neutrinos on high-resolution weak gravitational lensing convergence maps_ , Mon. Not. R. Astron. Soc. 524, 5591–5606 (2023).
  * Wong  _et al._ [2019] K. C. Wong, S. H. Suyu, G. C.-F. Chen, C. E. Rusu, M. Millon, D. Sluse, V. Bonvin, C. D. Fassnacht, _et al._ ,  _H0LiCOW – XIII. A 2.4 per cent measurement of H0 from lensed quasars: 5.3 $\sigma$ tension between early- and late-Universe probes_, Mon. Not. R. Astron. Soc. 498, 1420–1439 (2019).
  * Belgacem  _et al._ [2019b] E. Belgacem, G. Calcagni, M. Crisostomi, C. Dalang, Y. Dirian, J. M. Ezquiaga, M. Fasiello, S. Foffa, _et al._ ,  _Testing modified gravity at cosmological distances with LISA standard sirens_ , J. Cosmol. Astropart. Phys. 2019 (07), 024–024.
  * Hadzhiyska  _et al._ [2023] B. Hadzhiyska, S. Yuan, C. Blake, D. J. Eisenstein, J. Aguilar, S. Ahlen, D. Brooks, T. Claybaugh, _et al._ ,  _Synthetic light-cone catalogues of modern redshift and weak lensing surveys waith ABACUSSUMMIT_ , Mon. Not. R. Astron. Soc. 525, 4367–4387 (2023).
  * Harnois-Déraps  _et al._ [2018] J. Harnois-Déraps, A. Amon, A. Choi, V. Demchenko, C. Heymans, A. Kannawadi, R. Nakajima, E. Sirks, _et al._ ,  _Cosmological simulations for combined-probe analyses: covariance and neighbour-exclusion bias_ , Mon. Not. R. Astron. Soc. 481, 1337–1367 (2018).
  * DeRose  _et al._ [2019] J. DeRose, R. H. Wechsler, M. R. Becker, M. T. Busha, E. S. Rykoff, N. MacCrann, B. Erickson, A. E. Evrard, _et al._ , The Buzzard Flock: Dark Energy Survey Synthetic Sky Catalogs (2019), arXiv:1901.02401 .
  * Korytov  _et al._ [2019] D. Korytov, A. Hearin, E. Kovacs, P. Larsen, E. Rangel, J. Hollowed, A. J. Benson, K. Heitmann, _et al._ ,  _CosmoDC2: A Synthetic Sky Catalog for Dark Energy Science with LSST_ , Astrophys. J., Suppl. Ser. 245, 26 (2019).
  * Kacprzak  _et al._ [2023] T. Kacprzak, J. Fluri, A. Schneider, A. Refregier, and J. Stadel,  _CosmoGridV1: a simulated $\Lambda$CDM theory prediction for map-level cosmological inference_, J. Cosmol. Astropart. Phys. 2023 (2), 050.
  * Wei  _et al._ [2018] C. Wei, G. Li, X. Kang, Y. Luo, Q. Xia, P. Wang, X. Yang, H. Wang, _et al._ ,  _Full-sky ray-tracing simulation of weak lensing using ELUCID simulations: Exploring galaxy intrinsic alignment and cosmic shear correlations_ , Astrophys. J. 853, 25 (2018).
  * Pires  _et al._ [2020] S. Pires, V. Vandenbussche, V. Kansal, R. Bender, L. Blot, D. Bonino, A. Boucaud, J. Brinchmann, _et al._ ,  _Euclid: Reconstruction of weak-lensing mass maps for non-gaussianity studies_ , Astron. Astrophys. 638, A141 (2020).
  * LSST Dark Energy Science Collaboration [2012] LSST Dark Energy Science Collaboration, Large Synoptic Survey Telescope: Dark Energy Science Collaboration (2012), arXiv:1211.0310 .
  * Laureijs  _et al._ [2011] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, _et al._ , Euclid Definition Study Report (2011), arXiv:1110.3193 .
  * Hall [2021] A. Hall,  _Cosmology from weak lensing alone and implications for the Hubble tension_ , Mon. Not. R. Astron. Soc. 505, 4935–4955 (2021).
  * Wagoner  _et al._ [2021] E. L. Wagoner, E. Rozo, X. Fang, M. Crocce, J. Elvin-Poole, and N. W. and,  _Linear systematics mitigation in galaxy clustering in the dark energy survey year 1 data_ , Mon. Not. R. Astron. Soc. 503, 4349–4362 (2021).
  * Corman  _et al._ [2022] M. Corman, A. Ghosh, C. Escamilla-Rivera, M. A. Hendry, S. Marsat, and N. Tamanini,  _Constraining cosmological extra dimensions with gravitational wave standard sirens: From theory to current and future multimessenger observations_ , Phys. Rev. D 105, 064061 (2022).
  * Dupletsa  _et al._ [2023] U. Dupletsa, J. Harms, B. Banerjee, M. Branchesi, B. Goncharov, A. Maselli, A. Oliveira, S. Ronchini, and J. Tissino,  _gwfish: A simulation software to evaluate parameter-estimation capabilities of gravitational-wave detector networks_ , Astron. Comput. 42, 100671 (2023).