Impact of weak lensing on bright standard siren analyses

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 degenerate modification to the inferred luminosity distance $d_L$. To address this the uncertainty on $d_L$ is increased according to the dispersion of the magnification distribution at the source redshift, $\sigma_\mu$. But this term is dependent on cosmological parameters that are being constrained by gravitational wave"standard sirens", such as the Hubble parameter $H_0$, and the matter density fraction $\Omega_m$. $\sigma_\mu$ 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 outdated cosmological simulations, suggest $\sigma_\mu$ could be underestimated. We consider two classes of standard siren, supermassive black hole binary and binary neutron star mergers. Underestimating $H_0$ and $\Omega_m$ when calculating $\sigma_\mu$ increases the probability of finding a residual lensing bias on these parameters greater than $1\sigma$ by 1.5-3 times. Underestimating $\sigma_\mu$ by using low resolution/small sky-area simulations can also significantly increase the probability of biased results. For neutron star mergers, the spread of possible biases is 0.25 km/s/Mpc, comparable to the forecasted uncertainty. Left uncorrected this effect limits the use of BNS mergers for precision cosmology. For supermassive black hole binaries, the spread of possible biases on $H_0$ is significant, 5 km/s/Mpc, but $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.


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 ′ = h √ µ, where a prime denotes a lensed quantity [6][7][8][9][10].The magnification µ depends on the matter along the line of sight where µ = 1 corresponds to no effect, µ < 1 a net underdensity of matter, and µ > 1 a net overdensity.The magnification distribution over the sky, p(µ, z), broadens towards higher redshifts, and is peaked at values below one as there are more voids than clusters.The tail of p(µ, 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 * c.mpetha@ed.ac.uk 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 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 ∼ 2 [23].Weak lensing will not be as important for individual BNSs due to their larger distance uncertainties compared to those of SMB-HBs.On the other hand, with the 3G network expected to detect many more BNSs than the numbers of SMB-HBs 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.
Past works have studied how successfully a GW can be "delensed" using external data to estimate µ 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.An alternative delensing prospect which uses GW data alone is through measuring wave-optics features in a single wave-form [27].However the expected probability for observing these features is low [28], and this approach is not able to directly probe the range of lens masses which contribute most to weak lensing magnification.
Measuring the distance and redshift of a binary merger is the core of using these sources for cosmology, as "standard sirens" [29].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 [30][31][32][33][34][35][36][37].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 [38][39][40][41][42][43][44][45][46][47][48][49][50].
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, where the Hubble parameter is + Ω DE (1 + z) 3(1+w0+wa) e −3wa z 1+z .
Here H 0 = 100 h km s −1 Mpc −1 is the present expansion rate, Ω m , Ω K and Ω 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 [51,52].
The luminosity distance is inversely proportional to the observed wave amplitude.Therefore the observed, lensed distance is The total distance uncertainty is where ∂d L /∂z converts redshift to distance uncertainties.σ GW and σ vpec are the uncertainty from the detector and source peculiar velocity respectively.The weak lensing term is σ µ (z) is the uncertainty in the magnification of the gravitational wave.It has been assumed that the mean of the magnification distribution ⟨µ⟩ = 1.For real observations this is in general not the case due to detector selection effects as explored in Ref. [53], henceforth CT21.The consequence is ⟨µ⟩ ∼ 1.02 at z = 2. Weak lensing could lead to a net over-or underestimation of luminosity distances at different (true) redshifts, so that the inferred cosmological parameters in Eq. ( 2) are prone to biases [50].A mitigation might be to include a weak lensing uncertainty according to Eq. ( 5), using an assumed functional form for σ µ (z) [53][54][55][56][57][58].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 [59] (recently exacerbated by JWST observations [60]), and present weak lensing fitting functions have been created with cosmological simulations that could be missing large and small-scale power important in estimating σ µ (z) [61] (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. [50] 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σ for BNSs (SMBHBs).We also demonstrate the importance of the resolution of simulations used to determine σ µ (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σ.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 ∼ 5 km s −1 Mpc −1 , but will only become relevant with 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 σ µ (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. Set a true cosmology by choosing values for cosmological parameters in Eq. ( 2).
2. Choose a population: supermassive black hole binaries (SMBHBs) with a particular seed model observed with LISA [20] (see Section III A), or binary neutron stars (BNSs) observed with third generation ground based detectors (see Section III B).
Steps 3 − 7 are repeated over many catalogue realisations: independent sets of observations drawn from the population model.
3. Choose a catalogue realisation.For each merger in the catalogue, randomly draw a magnification µ as described in Section IV from the redshift and cosmology dependent magnification distribution.Applying the magnification to the measured merger distance and its uncertainty σ 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 µ = 1 for all sources.Only sources with a SNR greater than a chosen threshold after applying the magnification go into the catalogue.
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 ′obs 5. In both the lensed and unlensed cases include a weak lensing uncertainty σ WL (z), as in Eq. ( 4), to the distance.The assumptions going into generating σ WL (z) are the focus of this investigation.

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 [62].
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 package1 [63,64].This uses information from the full h − Ω m 2D posterior to quantify the size of the bias, and does not assume Gaussianity in the parameter space.
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. Repeat steps 3 − 8 a further ∼ 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 A details the subset of SMBHB catalogues from Ref. [20] used in this work, while the BNS catalogues we generated are described in Section III B. 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. [65] demonstrate that using several hundred 3G-observed 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 [66].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.

A. LISA: Supermassive black hole binaries
We use a subset of the catalogues presented in Refs [20,67], refer there for more detail on the catalogue generation, and to Ref. [67] 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 [68], the "Q3" heavy seed model where massive black holes are formed from the collapse of proto-galactic disks [69,70] 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 coun-terpart 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 ELT2 .If the host galaxy magnitude m gal,ELT < 27.2 and 0.5 < z gal ≤ 5 then σ z = 5 × 10 −4 .If 27.2 < m gal,ELT < 31.3 and 0.5 < z gal ≤ 5 (z gal > 5) then σ z = 0.1 (σ 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.

B. 3G: Binary neutron stars
For 3 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 [71], so this population is not included.BNSs are generated using the assumptions in Table 1 of Ref. [72], 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 ⊙ .The adopted redshift distribution of BSS is described in Ref. [35].The Fisher-matrix code GWFAST [73] 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.
O(10 5 ) BNSs will be observed every year by a 3G detector network [72,74].Predictions on the corresponding number of BSS vary significantly, from 15 yr −1 to 3500 yr −1 [23,[75][76][77][78][79][80] 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 ρ lim = 12 [81].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 deg 2 [82,83].We also limit the distance uncertainty to σ d L /d L < 30%, around the upper limit of the z − σ d L /d L fitting function in Ref. [84] 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 σ d L /d L .More recent studies use information on the inclination angle ι from EM observations to break the d L − ι degeneracy and improve σ d L /d L by several factors [85][86][87].This was performed for GW170817 using the radio counterpart, improving the H 0 constraint by a factor of two [88].In the context of LISA, Ref. [89] show how from sky localisation the wave amplitude and inclination uncertainties can both be improved by a factor of ∼ 2 when ι ≁ 90 o .This factor increases to ∼ 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 ι ≤ 60 o and ι ≥ 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 km s −1 .Finally, we assume spectroscopic redshift uncertainties of σ z = 5 × 10 −4 (1 + z), based on electromagnetic follow up from present and forthcoming spectroscopic redshift surveys e.g DESI [90], ELT, and MegaMapper [91].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.

A. Background
The magnification is given by The convergence κ is the integrated surface mass density along the line of sight, and γ is the shear field.We are interested in the variance of the magnification at a particular redshift σ 2 µ (z), as this provides σ 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 [93].In that case, at a given redshift, the variance in the convergence field is directly related to the variance in the magnification [58], The convergence field is cosmology dependent.The cosmological parameters influencing the convergence field are θ = {h, Ω m , Ω b , Σm ν , n s , A s , Ω K , w 0 , w a }.Ω b is the baryon density fraction, Σm ν 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, Knowing the dispersion, a magnification probability distribution, p(µ, z | θ), can be constructed.Poorly understood small-scale physics significantly impacts the high-magnification tail of p(µ, z | θ).To avoid introducing modelling uncertainty, we restrict our analysis to the weak lensing regime, with µ ≤ 1.75 3 .Therefore a simple shifted lognormal model as presented in Ref. [58] is sufficient, but for comparison we also use p(µ, z) from the N side = 4096 and N side = 16384 simulations described in Ref. [94] 4 (henceforth T17).
To include the impact of selection effects, we must modify the magnification probability distribution.Following Refs [53,95], where The SNR is given by ρ, N (z, ρ) is the number density of sources as a function of redshift and SNR, ρ lim is the SNR threshold and 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.

B. Numerical calculation
Ref. [58] provide a tool to generate the dispersion in a redshift and cosmology dependent convergence field, σ κ (z | θ), for varying cosmological parameters using the Boltzmann code camb 5 .We made minor modifications to their notebook; using the HMcode-2020 nonlinear model [96], setting k max = 100 h 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 [96], 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 ℓ to find the total variance.This integral convergences for ℓ max = 10 7 [58], and while we can not expect to accurately predict the power spectrum at such Planck 2018 base ΛCDM constraints from TT,TE,EE+lowE+lensing+BAO [66] (henceforth labelled 'Planck'), and the best fit parameters from the Pantheon+ analysis assuming either ΛCDM or wCDM [92].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, Ω m , w 0 and w a used for each cosmology is given in Table I.

C. Simulations
The present industry standard weak lensing fitting function [53][54][55][56][57]97] was created in 2010 [97] using the magnification distribution of the 1998 work Ref. [98].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, Ω m = 0.274).Ref. [99] demonstrate how the inclusion of baryons in hydrodynamic simulations modifies the high magnification tail (µ > 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 [100] and MilleniumTNG [101] simulation suites suggest that baryons, and in particular feedback effects, have a larger impact, causing a narrower p(µ, 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 Ω m predict a smaller σ κ (z | θ).More recent constraints on Ω m favour larger values [66,92] than used in T11, and there is the possibility h favours the near-Universe measurement [92,102].Therefore current weak lensing fitting functions in the literature could be underestimating the dispersion of p(µ, z | θ).
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. [97], is given by where C = 0.061, β = 0.264, α = 1.89 for the popIII seed model and C = 0.096, β = 0.62, α = 2.36 for the Q3 model.Note the factor of 1/2, which is missing in Ref. [55].This is because the original expression is for the uncertainty in σ(ln 53,65].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 ℓ min = 1 to ℓ max = 10 7 .The ray-tracing simulations of T11 use a 1 × 1 deg 2 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 ℓ min , while ℓ max increases with redshift but is smaller than 10 7 .Hence the variance found from these simulations will be smaller.Setting ℓ min and ℓ max in the numerical calculation based on the simulation parameters leads to comparable predictions for σ WL (z) to Ref. [97] 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 fullsky, providing access to lower multipoles, but have significantly lower resolution in the high ℓ regime, and are available up to a maximum redshift of 5.3, compared to z max = 20 in T11 (though there are only 7 redshift bins from z = 1 to z = 20).They use a similar set of cosmological parameters to T11, with h = 0.7, Ω m = 0.279.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 side , can be seen in Fig. 3.Both the convergence and shear fields are used to find p(µ) from these simulations.There is a marked difference in the width of p(µ), highlighting the importance of resolution (or equivalently, ℓ max ) in estimating the weak lensing uncertainty.Also shown is p(µ) from T11 at z = 1, which shows close similarity to the N 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. [97] can be seen in Fig. 4. The weak lensing 6 Incidentally, including this factor of 1/2 would improve results in numerous forecasting papers, e.g.Refs [32,35,55,56,76,103]. uncertainty found from the full-sky T17 N side = 16384 simulations follows the fitting function of Eq. ( 11) with C = 0.076, β = 0.344 and α = 2.25 with no selection effects, and C = 0.092, β = 0.164, α = 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 ≤ 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.
Weak lensing data products from other cosmological simulations exist [104][105][106][107][108][109] 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. [54].The authors use the ELUCID simulations [109] to estimate σ WL (z) for observations of binary black holes at z ≤ 1.For z > 1 they also use results from T11, hence find general agreement with the curves in CT21.

D. Observations
Ideally, p(µ, z) would be inferred from observations.Through the weak lensing of galaxies we measure reduced shear g = γ/(1 − κ), which is either approxi- mated as shear, or γ is recovered through an iterative scheme [110].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(µ, 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 [111] or Euclid [112].
V. RESULTS

A. 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 ⟨µ⟩ = 1 and the mean bias ∼ 0. But with low SNR sources that can be shifted inside the SNR limit through a weak lensing magnification, ⟨µ⟩ ̸ = 1, and there is a bias in the recovered cosmological parameters.This depends primarily on the source redshift distribution and the magnification TABLE II.The magnification selection effect on gravitational waves causes a bias 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.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 Ω m due to selection effects in Table II.Results are shown for the three cosmologies used to generate the numerical estimate of σ WL (z), and when magnifications are drawn directly from the N side = 4096 and N side = 16384 simulations of T17.In the N side = 4096 case the mean bias is smaller as the limited resolution leads to a smaller ℓ max than in the N 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.
Next, we go on to consider the impact of the assumption on σ WL (z).This will primarily affect the width of the bias distribution.From Table II 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 (Ω m ) exceeding a value of ±0.0025 (±0.015) for BNSs and ±0.02 (±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.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 III, 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.
This scenario highlights the necessity of including σ WL (z).Any cosmological analysis using well measured binary mergers beyond z ∼ 0.5 that does not could be significantly biased.

Correct uncertainty implemented
Now we test the situation where σ WL (z) is generated using the same distribution that magnifications are drawn from.There is perfect knowledge of the true p(µ, z | θ).This does not mean the magnification can be removed from the GW, for this an external measurement of µ is needed for each source.Magnifications are drawn from p(µ, z | θ) generated using Planck parameters and the added uncertainty is the Planck curve in Fig. 2. Results are presented in Table IV.
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 III.
In the case of SMBHBs, for a bias b ≥ 1σ the probability is ≤ 0.05, and for b ≥ 2σ the probability is ≤ 0.01.The lensing bias has been effectively removed by increasing uncertainties according to the dispersion of p(µ, z | θ).
For BNSs we are more likely to observe sources with a positive magnification.There is an average negative/positive shift on h/Ω m , also shown in Ref. [50].Due to this selection effect the probability of a residual bias b ≥ 1σ (b ≥ 2σ) is 0.18 (0.01).The lensing se- lection 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 ⟨d 2 L (z)⟩ in each bin and multiplying by ⟨µ(z)⟩ 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 ⟨d L (z)⟩, and crucially, accurate knowledge of p(µ, z | θ).As will be demonstrated in the next section, mischaracterising the magnification distribution could exacerbate the lensing bias instead of mitigating it.

Imperfect correction
In reality, we do not know the correct form of p(µ, z | θ).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 III will not be mitigated, as it has mostly been in Section V A 2. There are two main reasons the weak lensing uncertainty could be underestimated: A. Using too small values of H 0 or Ω m to predict σ 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(µ, z | θ) is either the Pantheon+ ΛCDM or Pantheon+ wCDM 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 V.A key result is the importance of accurate modelling of p(µ, z | θ) 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 ΛCDM, including kΛCDM (curvature) and w 0 CDM (dark energy) in Appendix C. For LISA FIG. 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+ Λ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 Ωm when all sources experience magnification µ = 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σ, or an absolute bias of b(h) = −0.007and b(Ωm) = 0.031.
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.
To assess the impact of underestimating cosmological parameters, probabilities in Table V are compared with those in Table IV, which describe when the weak lensing uncertainty is correctly estimated.First consider if the true cosmology is Pantheon+ ΛCDM, and assumed is Planck.For all populations, the probability of the absolute bias on h and Ω 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σ is increased by a factor ∼ 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σ has increased by a factor of 1.5, and is close to 1/3.If the true cosmology is in fact Pantheon+ wCDM, probabilities increase by an extra factor of ∼ 2.
To test case B, we compare the fitting functions of CT21 with curves created using camb.σ WL (z | θ) generated using the numerical method of Section IV B and the cosmological parameters in T11 and T17 are very similar to the Planck curve in Fig. 2 (orange dash-dot line).This is because the increase in Ω m and decrease in h between the parameters used in the simulations and Planck are competing effects and there is a partial cancellation of biases.The difference in σ WL is ≲ 5%, 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 VI for each population.Once again, as the true cosmology is the same as in Section V A 2, we can compare bias probabilities with Table IV.In this case the uncertainty is too small to wash out the bias.For SMBHBs the probability of results biased by ≥ 1σ increases by a factor of ∼ 10, and there is now a significant probability of results biased by ≥ 2σ.For BNSs, there is now greater than 50% probability of finding results biased by more than 1σ and for b ≥ 2σ 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 σ WL (z | θ) in fact overestimates the dispersion in p(µ, z | θ), then the uncertainty on the cosmological parameters is excessively diluted.When the true cosmology is Planck and the assumed cosmology is Pantheon+ ΛCDM (wCDM), the increased uncertainty at low redshift causes a degradation of parameter constraints of 10% (15%)/8% (12%) for SMBHBs/BNSs.

Setting requirements
In Section V A 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 Ω 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 σ WL (z | θ) leads to smaller parameter uncertain- ties).And so the scaling of the average bias in units of σ will simply depend on the scaling of σ, 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 Ω 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 Ω 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 [113], but H 0 should be known in combination with Ω m to ∼ 2 km s −1 Mpc −1 .In the right plot we also consider the amplitude of the linear matter power spectrum in spheres of radius 8 Mpc −1 h −1 , σ 8 .We find that, if Ω m is to be known to 0.02 or better, then σ 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 σ 8 and Ω m are well within parameter uncertain-ties from e.g.Planck, and are consistent with differences in these parameters observed from different probes (e.g Refs [5,66,92]).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 Ω 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.
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 [94], using the expected weak lensing uncertainty from each N side .We also downsample the N side = 4096 simulations to N side = 2048 and N side = 1024.We find a linear scaling between the average bias and log 2 (N 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 ≥ 1.For every power of 2 increase in N side , the average bias drops by ∼ 0.1σ where σ 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 [104].Weak lensing maps for full hydro-dynamical simulations with baryons have been calculated in Ref. [101], but for a single resolution of N 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 σ WL (z | θ).A feedback model, HM-code2020 feedback, is already included in camb [96], but FIG. 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.
this model was not used in this work due to a divergence in the variance of the power spectrum at small scales.

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 σ 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 σ WL (z | θ) is added as an ingredient to the MCMC analysis with each new sampling of h, Ω 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 σ 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. [114].
A realistic prior on A is needed.The fitting will try to increase the variance as much as it can, as this will minimise (data − model) 2 /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 biasvariance 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 ∼ 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 ≳ 0.5 not including a weak lensing term has a high probability of being biased by more than 1σ.
• 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 × 1 deg 2 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 C, noting this reduces the probability of a residual lensing bias for z ≤ 5.3 bright standard sirens, such as merging binary neutron stars.
• If cosmological parameters are underestimated when generating σ WL (z), the probability of parameter constraints biased by more than 1σ 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 raytracing 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(µ) is used.The probability of a 1σ (2σ) bias in H 0 and Ω m found from BNS mergers is 54% (19%).
• Considering the small absolute value of the bias for BNSs, with a spread of ∆b(H 0 ) ∼ 0.25 km s −1 Mpc −1 , magnification selection effects or incorrect assumptions on the true form of p(µ | z, θ) 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 ∆b(H 0 ) ∼ 5 km s −1 Mpc −1 would be significant.An order of magnitude more SMBHBs than considered in these catalogues (from ∼ 20 to ∼ 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(µ | z, θ) using H 0 , Ω m and σ 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 σ WL (z) is using high-resolution weak lensing maps, including baryons, over a range of resolutions.A convergence of the width of p(µ, 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 [104].
A similar, independent analysis was performed in Ref. [50].There the true cosmology is that used in the AbacusSummit simulations [104] with h = 0.6736, Ω m = 0.3138, similar to the Planck values used in this work, and the weak lensing uncertainty is the fitting function of Ref. [55].The magnification distributions from AbacusSummit, as seen in Fig. 1 of Ref. [50], are broader than those from the simulations of T11 from which the WL uncertainty is derived, contributing to their derived bias.Ref. [50] include larger lensing magnifications while we limit to µ ≤ 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 σ WL (z) is treated as a free parameter.This is a hybrid approach, using the information from observations to avoid fixing σ WL (z) based only on simulations, and was originally envisioned in Ref. [32].Also, as pointed out in Ref. [30], the relative abundance of magnified and de-magnified sources reveals the non-Gaussianity of p(µ | 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 [76,103,115].

ACKNOWLEDGMENTS
We thank Alberto Mangliagli for use of their SMBHB catalogues.CTM thanks Ryuichi Consider two cases, a source at d L with no magnification (µ = 1), and the same source detected by the same detector but with a magnification µ.ΛCDM causes significant biases.The probability of a discordance in the data set ≥ 1σ (≥ 2σ) is 0.54 (0.01).An example of the lensing bias in this cosmological model can be seen in Fig. 8.

FIG. 3 .
FIG.3.The magnification probability distribution at z = 1.0334 found from the full-sky ray-tracing simulations of T17.Increasing the N 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 × 1 deg 2 ray-tracing simulations.

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

FIG. 8 .
FIG.8.Similar to Fig.5but for the w0CDM cosmological model, where we are also constraining a constant dark energy equation of state.In this case lensing leads to a 1.2σ bias in this parameter space, and a > 2σ preference for a phantom dark energy.

TABLE I .
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.

TABLE III .
[66] 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 Ωm exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D h − Ωm posterior, are shown.The true cosmology used is Planck 2018[66].For BNSs x h = 0.0025, xΩ m = 0.015.For SMBHBs x h = 0.02, xΩ m = 0.05.

TABLE IV .
Mean probability of finding a lensing bias.Both the probability of the bias on h and Ωm exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D h − Ω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Ω m = 0.015.For SMBHBs x h = 0.02, xΩ m = 0.05.

TABLE V .
Mean probability of finding a lensing bias.Both the probability of the bias on h and Ωm exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D h − Ω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Ω m = 0.015.For SMBHBs x h = 0.02, xΩ m = 0.05.≥ ±x h ) P (b(Ωm) ≥ ±xΩ m ) P (b ≥ 1σ[2σ]) P (b(h) ≥ ±x h ) P (b(Ωm) ≥ ±xΩ m ) P (b ≥ 1σ[2σ])

TABLE VI .
Mean probability of finding a lensing bias.Both the probability of the bias on h and Ωm exceeding a set value, and the probability of the bias exceeding a number of standard deviations in the 2D h − Ω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Ω m = 0.015.For SMBHBs x h = 0.02, xΩ m = 0.05.
Takahashi for helpful discussion and feedback, and Jack Elvin-Poole and Pierre Burger for useful discourse.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