Joint cosmological inference of standard sirens and gravitational wave weak lensing

TL;DR

This paper demonstrates a novel joint inference method combining standard sirens and gravitational wave weak lensing to improve cosmological parameter constraints without external priors, especially in a future post-LISA scenario.

Contribution

It introduces the first combined analysis of standard sirens and gravitational wave weak lensing from the same dataset, enhancing constraints on cosmological parameters.

Findings

Achieves 2-3% precision on dark matter and clustering parameters

Constrains Hubble constant to better than 1%

Provides independent cosmological constraints that could resolve existing tensions

Abstract

We present the first joint inference of standard sirens and gravitational wave weak lensing by filtering of the same dataset. We imagine a post-LISA scenario emerging around the late 2030s when LISA will have accumulated a number of detections at high redshift; LIGO-VIRGO will have finished observing at low redshift, and Einstein Telescope will have started making new observations out to redshifts possibly overlapping with LISA. Euclid and other cosmological probes will have provided constraints at percent level by then, but mostly exhausted their ability to improve any further. We derive forecasts assuming $\sim1\,\text{deg}^{-2}$ detected sources, in conjunction with a spectroscopic follow-up (e.g. by Euclid, DESI, or ATHENA). Thanks to the statistical power of standard sirens as a geometry probe -- lifting key degeneracies in the gravitational wave weak lensing -- and no external priors assumed, the constraints on dark matter and its clustering, namely $\Omega_m$ and $\sigma_8$, could be achieved to $2\%$ and $3\%$. The Hubble constant could be constrained to better than $1\%$ in all cases; the dark energy density, $\Omega_\Lambda$, to $2\%$, and curvature, $\Omega_K$, to $0.02$; the amplitude and spectral tilt of the scalar fluctuations, $\ln(10^{10}A_s)$ and $n_s$, to $2\%$ and $7\%$. As a completely independent cosmological probe, with less calibration requirements, the joint inference of standard sirens and gravitational wave weak lensing might help solve the tensions currently observed between other cosmological probes, such as CMB, galaxy lensing, and Type Ia SNs, and distinguish between residual systematics and new physics.

Tables

Table 1. Table 1: Forecasted 1D marginal constraints for flat Λ Λ \Lambda CDM (left columns) and Λ Λ \Lambda CDM with curvature (right columns). Shown are fractional percent errors for each of the relevant cosmological parameters that are considered in this paper. Standard sirens, d L subscript 𝑑 𝐿 d_{L} , bring in the best constraint on geometry ( h ℎ h and Ω i subscript Ω 𝑖 \Omega_{i} ), but do not say anything about clustering. However, when combined with gravitational wave weak lensing, 𝒞 ℓ κ subscript superscript 𝒞 𝜅 ℓ \mathcal{C}^{\kappa}_{\ell} , clustering ( σ 8 subscript 𝜎 8 \sigma_{8} , A s subscript 𝐴 𝑠 A_{s} , and n s subscript 𝑛 𝑠 n_{s} ) can also be improved to percent level.

	$h$		${\Omega }_m$		${\sigma }_8$		$\ln \left({10}^{10}{A}_s\right)$		$n_s$		${\Omega }_{\Lambda }$	${\Omega }_K$
$d_L$	$0.21$	$1.1$	$1.2$	$2.4$	-	-	-	-	-	-	$2.0$	$0.020$
${𝒞}_{\ell }^{\kappa }$	$77$	$390$	$61$	$44$	$27$	$23$	$22$	$290$	$34$	$130$	$26$	$0.16$
$d_L$ + ${𝒞}_{\ell }^{\kappa }$	$0.21$	$1.1$	$1.1$	$2.4$	$2.7$	$2.7$	$1.9$	$2.2$	$6.7$	$6.7$	$2.0$	$0.019$

Equations

$$d_L(z)=(1+z)f_K(\chi ),$$

$$f_K(\chi )=\frac{1}{\sqrt{-K}}\sinh \left[\sqrt{-K}\chi \right].$$

$$\chi (z)=c{\int }_0^z\frac{\text{d}{z}^{\prime }}{H(z^{\prime })},$$

$${\sigma }_z^2=\left({\sigma }_{\text{sp}}^2+\frac{{\sigma }_v^2}{c^2}\right)(1+z{)}^2,$$

$$\mu =\frac{1}{(1-\kappa {)}^2-|\gamma {|}^2}\sim 1+2\kappa ,$$

$$d_L^{\prime }=(1-\kappa )d_L,$$

$${\mathcal{C}}_{\ell }^{\kappa }={\int }_0^{{\chi }_d}\text{d}\chi W^2(\chi ){\mathcal{P}}_{\delta }\left(k=\frac{\ell +1/2}{f_K(\chi )},\chi \right),$$

$$W(\chi )=\frac{3}{2}{\left(\frac{H_0}{c}\right)}^2\frac{{\Omega }_m}{a(\chi )}{\int }_{\chi }^{{\chi }_d}\text{d}{\chi }^{\prime }n({\chi }^{\prime })\frac{f_K({\chi }^{\prime }-\chi )}{f_K({\chi }^{\prime })},$$

$${\Sigma }_{\ell \ell }^{\kappa }=\frac{2}{f_{\text{sky}}(2\ell +1)\delta \ell }{\left({\mathcal{C}}_{\ell }^{\kappa }+\frac{{\sigma }_{\kappa }^2}{\bar{n}}\right)}^2,$$

$$\begin{array}{rl}p({\mathcal{C}}_{\ell }^{\kappa },\kappa ,d_L|\theta )& =p({\mathcal{C}}_{\ell }^{\kappa },\kappa |d_L,\theta )p(d_L|\theta )\\ & =p({\mathcal{C}}_{\ell }^{\kappa }|\kappa ,d_L,\theta )p(\kappa |d_L,\theta )p(d_L|\theta ).\end{array}$$

Full text

I Introduction

The gravitational wave (GW) distance-redshift relation can be exploited to derive independent cosmological constraints on $H_{0}$, as demonstrated by the LIGO-VIRGO (LV) collaboration with the first multi-messenger event Abbott et al. (2017). Moreover, GW measurements can now be combined with Planck Aghanim et al. (2018a) to further improve constraints on other parameters, such as neutrino mass, curvature, and dark energy equation of state Di Valentino and Melchiorri (2018). In future, we will have collected a plethora of GW events, attained from the combination of different source catalogues, namely by LV, Einstein Telescope (ET) Punturo et al. (2010), and LISA Amaro-Seoane et al. (2017), reaching median redshifts of $\sim$2. Precious cosmological information will then be extracted by having access to luminosity distance, redshift, and sky position data.

In constrast to Type Ia supernovae (SN) analysis which is subjected to external calibration (the cosmic distance ladder) Riess et al. (2016, 2018), the distance to the GW source can instead be measured directly and accurately with GW data alone Schutz (1986). Owing to the degeneracy between mass and redshift, the corresponding redshift must however be derived from electromagnetic counterparts through multi-messenger events. The weak lensing (WL) from large-scale structure (LSS) is also known to affect the distance-redshift relation Holz and Linder (2005). The farther the sources are, the bigger the effect will be, hence the bigger the scatter of the luminosity distance being randomly magnified/demagnified by over/under matter densities. This percent level effect is comparable with the GW measurement error and has traditionally been seen as a source of potential systematic that could in principle be corrected for Hilbert et al. (2011); Hirata et al. (2010), or more recently as an actual cosmological probe on its own for future, more dedicated, GW experiments Camera and Nishizawa (2013); Cutler and Holz (2009).

Just as the cosmic microwave background (CMB) and the lensing of the CMB jointly provide additional constraining power for cosmogical inference, as demonstrated by the Planck satellite Aghanim et al. (2018b), in this paper we propose a joint analysis of standard sirens and weak lensing by filtering of the same GW data, such that the weak lensing inference effectively becomes conditional to the standard sirens data. We derive the first joint cosmological forecasts for data as it will be observed post the space-based detector LISA in the late 2030s, in conjunction with other detectors (e.g. LV and ET), and redshifts from spectroscopic follow-up surveys such as Euclid Laureijs et al. (2011), DESI Aghamousa et al. (2016), or ATHENA Nandra et al. (2014). We show how the method will provide competitive constraints, which will be alternative and complementary to galaxy surveys and CMB experiments. By that time the same experiments will have exhausted their ability to improve their cosmological constraints significantly below the percent level as it becomes more and more apparent in the era of systematics-dominated cosmology.

Throughout the rest of the paper, we assume the 2018 Planck cosmology, $H_{0}=100\,h\,\mathrm{km\,s^{-1}Mpc^{-1}}$, $h=0.673$, $h^{2}\Omega_{m}=0.143$, $h^{2}\Omega_{\Lambda}=0.310$, $\Omega_{K}=0$, $\ln(10^{10}A_{s})=3.05$, and $n_{s}=0.965$ Aghanim et al. (2018a).

II Standard sirens and weak lensing

Interferometric detectors are sensitive to the derivative of the light frequency shift induced by incoming GWs Congedo (2015); Congedo et al. (2013). Owing to amplitude and frequency modulations induced by the motion of the detector around the Sun (LISA) or Earth (LV and ET), the amplitude of the wave is particularly well measured. As the observed GW strain is inversely proportional to the luminosity distance Congedo (2017); Maggiore (2009), the luminosity distance is determined very accurately.

The standard sirens inference is based on relating the GW luminosity distance, $d_{L}$, to the optical redshift, $z$, as follows

[TABLE]

where the comoving angular distance is given by

[TABLE]

In the above equation, $K$ is the curvature and the comoving distance is

[TABLE]

where $H(z)=H_{0}[\Omega_{m}(1+z)^{3}+\Omega_{K}(1+z)^{2}+\Omega_{\Lambda}]^{1/2}$ and $\Omega_{K}=-(c/H_{0})^{2}K$.

With a nominal sky resolution of $\sim 1\,\mathrm{deg^{2}}$, it would, in principle, be possible to determine the redshift for the detected source, either from an electromagnetic counterpart, or an ensemble average over a population of galaxies selected in the same sky bin. Potentially, the redshift could also be inferred directly from GW data (see e.g. Ref. Bonvin et al. (2017) and references therein). Given that the redshift selection function of GW observations will be required to have good overlap with optical follow-up surveys, it is reasonable to assume in our analysis that (i) the electromagnetic counterpart can indeed be identified in the provided position error box; (ii) the redshift can be determined by a spectroscopic follow-up survey, which is feasible given the relative small number of GW sources involved in this analysis compared to galaxy surveys.

We adopt the following scaling relation for the expected total redshift error,

[TABLE]

which accounts for both the spectroscopic error $\sigma_{\text{sp}}=0.001$ as per the Euclid/DESI nominal requirements, and the peculiar velocity dispersion, being $\sigma_{v}\sim 300\,\mathrm{km/s}$ a reasonable estimate 111We note however that the peculiar motion can affect the total redshift error budget only for very shallow observations, whereas it becomes totally negligible at, e.g., the LISA depth..

Current forecasts for LISA Tamanini (2017); Tamanini et al. (2016) are based only on the $d_{L}$ versus $z$ relation with reasonable expectations of measurement errors similarly to the above. However, the forecasts for models beyond flat $\Lambda$CDM are usually derived by setting key parameters (e.g. $h$) to nominal values, instead of correctly marginalising over their probability distribution, and the constraints on the dark sector are generally not very strong.

As noted in Ref. Cutler and Holz (2009), the WL fluctuations can be directly inferred from GW data. The WL magnification factor is defined as

[TABLE]

where $\kappa$ and $\gamma$ are the convergence and shear fields owing to the lensing by the LSS, and the approximation holds in the weak regime, $\kappa,|\gamma|\ll 1$. The magnification factor describes how solid angles, hence angular diameter distances, are transformed under lensing and is proportional, to leading order, to the convergence field. Likewise, the luminosity distance is transformed as follows

[TABLE]

where primed denotes a lensed quantity. Note that $d_{L}^{\prime}$ is the only observed quantity, as the detected GW strain is $h\propto 1/d_{L}^{\prime}$. In effect, the GW luminosity is $L\propto\dot{h}^{2}=f^{2}/{d_{L}^{\prime}}^{2}$, where $f$ is the frequency of the wave. However, in general, $d_{L}$ is also affected by, e.g. cosmological and gravitational redshift Yoo and Scaccabarozzi (2016). Similarily, the phase of the GW signal suffers from analogous deviations Bonvin et al. (2017). These effects are not included in our analysis. They may well be treated as systematic errors or even, potentially, as source of extra information Bonvin et al. (2006).

As the rms of the $\kappa$ field is $\sim 0.05$ at $z\sim 2$, this induces a typical lensing error on $d_{L}$ of $5\%$ (incidentally, this is comparable with or bigger than the GW measurement error on $d_{L}$ that LISA will be able to achieve, as discussed later on in this section). We adopt, however, the more realistic fitting formula of Refs Tamanini et al. (2016); Hirata et al. (2010), which predicts the lensing fluctuation on $d_{L}$, as a function of redshift.

We can now derive an estimate of the GW-WL convergence by adopting a fiducial model for the unlensed $d_{L}$, and inverting Eq. (6). In fact, assuming a typical measurement error on $d_{L}$, we would easily get a point estimate of $\kappa$ of significance $S/N>1$ with a single GW detection, as opposed to optical WL where an ensemble average of $\sim 10^{3}$ identical sources or more are required to get the same statistical power. The consequent reduction in the effective survey sample is therefore $\sim 10^{3}$.

Although noisy and conditional to the measured $d_{L}$, the GW estimate of $\kappa$ would be sufficient to calculate a 2D power spectrum. We define the convergence power spectrum as follows Kilbinger (2015)

[TABLE]

where $\mathcal{P}_{\delta}(k,\chi)$ is the 3D power spectrum of matter density fluctuations, $k$ is the Fourier mode, $\ell$ is the spherical harmonic multipole, $\chi_{d}=\chi(z_{d})$, and $z_{d}$ is the survey depth. The weight function defining the lensing efficiency is given by

[TABLE]

where $a(\chi)$ is the scale factor, $n(\chi)=n(z)|d\chi/dz|^{-1}$ is the source redshift distribution, and $\bar{n}=\int n(z)\text{d}z$ is the survey mean number density per steradian.

The lowest multipole accessible with a WL analysis of GW observations is limited by cosmic variance, which is determined by the sky coverage of the optical survey, say $\ell_{\min}\sim 2$. The highest multipole is instead limited by shot noise, which is given by the GW angular resolution. Assuming a realistic error box of $\sim 1\,\mathrm{deg^{-2}}$, that brings us to $\ell_{\max}\sim 100$, which is in contrast to other WL forecasts that assume $\ell_{\max}\sim 1,000$ or more. Therefore a WL analysis of GW data would only probe the large scales, and not be sensitive to the fully non-linear scales typical of galaxy WL surveys. In general, the total GW-WL noise variance is given by

[TABLE]

where $\delta\ell$ is the multipole resolution, $f_{\text{sky}}\sim 0.4$ is the sky coverage that would be attained by an optical follow-up by a galaxy survey like Euclid, $\sigma_{\kappa}$ is the total rms error on our estimate of convergence. The first term is the cosmic variance and the second term is the shot noise. Because the shot noise variance is $\propto\sigma_{\kappa}^{4}/\bar{n}^{2}$ and $\bar{n}$ is expected be small, GW-WL will in general be shot noise dominated. In this case, the total signal-to-noise ($S/N$) will be $\propto\bar{n}$.

We assume a nominal $\sim 1\%$ error on luminosity distance from GW detection. To put this figure into context, the LV multi-messenger event had an error of $\sim 15\%$ Abbott et al. (2017), whereas this is expected be around $1\%$ for stellar binary black holes Kyutoku and Seto (2017) and $5\%$ for extreme mass ratio inspirals Babak et al. (2017), both in the LISA band. Supermassive binary black holes will likely be observed in a greater number and with a much higher $S/N$, hence resulting in a similar or even better precision. We include spectroscopic redshift errors and peculiar velocities as they also contribute to the total error budget, but are in general subdominant. We adopt an $n(z)$ selection function as customary in galaxy surveys (see e.g. Ref. Laureijs et al. (2011)), but with a peak at $z\sim 2$, which agrees with predictions for e.g. LISA Tamanini (2017). We therefore estimate that a conservative number density $\bar{n}\sim 1\,\mathrm{deg^{-2}}$ would be enough to allow the detection of the GW-WL signal at a $S/N\sim 20$ significance level.

III Joint Fisher forecast

We adopt a Bayesian approach to combine the standard sirens and WL inference. We wish to derive the cosmological posterior probability given joint measurements of $d_{L}$, $\kappa$, and $\mathcal{C}^{\kappa}_{\ell}$, call this $p(\theta|\mathcal{C}^{\kappa}_{\ell},\kappa,d_{L})$, where $\theta$ is the set of cosmological parameters. Thanks to the Bayes theorem, this probability is $\propto p(\mathcal{C}^{\kappa}_{\ell},\kappa,d_{L}|\theta)\,p(\theta)$, where the former is the joint likelihood and the latter is the cosmological prior. By applying the conditional probability twice, the joint likelihood becomes

[TABLE]

In the first line, $p(\mathcal{C}^{\kappa}_{\ell},\kappa|d_{L},\theta)$ is the joint likelihood of power spectrum and convergence field conditional to standard sirens; $p(d_{L}|\theta)$ is the standard sirens likelihood. Upon further applying the conditional probability in the second line, $p(\mathcal{C}^{\kappa}_{\ell}|\kappa,d_{L},\theta)$ becomes the power spectrum likelihood conditional on both standard sirens and convergence field; finally, $p(\kappa|d_{L},\theta)$ is the convergence field likelihood conditional on standard sirens. The equation shows how a joint likelihood for power spectrum, convergence field and standard sirens can be derived. For simplicity, we will not, however, model $p(\kappa|d_{L},\theta)$ in our analysis, and defer that to future work. We just note that, as we are effectively restricted to the linear regime ($\ell_{\max}\sim 100$), a measure of the variance of the convergence map as constructed above can be used to constrain the amplitude of the clustering, through $\sigma_{8}$ (see e.g. Ref. Patton et al. (2017)). Moreover, $\mathcal{C}^{\kappa}_{\ell}$ and $\kappa$ (via $d_{L}$) are in general correlated. It turns out that this happens at the level of a 3-point statistic (or bi-spectrum), which may be significant only at highly non-linear scales. Once again, because we work in the linear regime, it is safe to ignore this correlation and assume that individual Fisher matrices can be safely summed up together.

For the purpose of this investigation, we define two baseline cosmologies: (a) the concordance flat $\Lambda$CDM model with free parameters $\theta=\{h,h^{2}\Omega_{m},\ln(10^{10}A_{s}),n_{s}\}$ and constraint $\Omega_{K}=0$; (b) the extended $\Lambda$CDM model with curvature, i.e. $\theta=\{h,h^{2}\Omega_{m},h^{2}\Omega_{K},\ln(10^{10}A_{s}),n_{s}\}$. In either cases, $h^{2}\Omega_{\Lambda}$ and $\sigma_{8}$, the amplitude of the linear matter power spectrum at $8\,\mathrm{h^{-1}Mpc}$, are derived parameters. In particular, $\sigma_{8}$ correlates strongly with $A_{s}$ and $n_{s}$, the amplitude and spectral tilt of the primordial scalar fluctuations. We note that standard sirens are a geometry probe, in that their parameter space is restricted to $\theta=\{h,h^{2}\Omega_{m},h^{2}\Omega_{K}\}$ and therefore this does not bring any information about the clustering.

We wrote Python Foundation (2018) code to derive our forecasts. We calculate the luminosity distance and its derivatives with respect to cosmology semi-analytically, whereas we get the power spectrum from CLASS Blas et al. (2011), whose numerical derivatives are robustly estimated with the numdifftools package Brodtkorb (2017). We assume that the likelihoods are Gaussian in their data, and all the measurement errors are Gaussian and uncorrelated. In particular, we assume the WL noise covariance in Eq. (9), and nominal GW errors of $1\%$ on $d_{L}$, plus redshift errors, peculiar velocities [see Eq. (4) and text thereafter], and lensing errors for standard sirens. We then calculate the individual Fisher matrices for standard sirens and WL, and the Fisher matrix for the corresponding joint analysis as the sum of the two. These Fisher matrices are derived for the physical parameters, $h^{2}\Omega_{i}$, as customary for standard cosmological measurements, and then mapped to density parameters, $\Omega_{i}$. To do that, we first draw samples from the initial Fisher matrices in the physical parameter space, and then map those samples to the density space via their non-linear transformation. In doing so, the following flat priors are assumed: $0\leq h\leq 1$, $0\leq\Omega_{m}\leq 1+\Omega_{K}$, $|\Omega_{K}|\leq 0.3$, $2\leq\ln(10^{10}A_{s})\leq 4$, and $0.5\leq n_{s}\leq 1.5$. This general procedure has the advantage that it can reproduce the typical degeneracies that are seen between, e.g., $\Omega_{m}$ and $\sigma_{8}$ when the noise is large. At the same time this procedure does not affect the final Fisher matrix estimation in any case, and can also be applied to any n-dimensional cosmologies. It does however require a Monte Carlo simulation over multiple cosmologies, hence calling the power spectrum calculation multiple times to compute derived parameters such as $\sigma_{8}$.

IV Results and discussion

We constructed 2D marginal contour plots of the relevant parameters, whose samples were drawn with a Monte Carlo simulation, following the procedure described in the previous section. In what follows, shown are one and two sigma contour plots for individual and combined analyses for the curved cosmology. For the sake of completeness, results are also summarised in Table 1 where 1D marginal errors are reported for both cosmologies and individual probes ($d_{L}$, $\mathcal{C}^{\kappa}_{\ell}$, and jointly).

We start from the $\Omega_{m}$-$h$ contour plot [see Fig. 1, panel (a)]. We note that because of shot noise at such a small number density, the WL constraints are generally poorer than standard sirens. Also, because of noise and the non-linear mapping from physical parameters to density parameters, the WL contour emerges as slightly shifted (i.e. biased) at $\sim 0.5\,\sigma$ level from the nominal cosmology. This is a noise bias effect that is absolutely expected here given the low $S/N$ of the WL observable. We checked that this is automatically reduced with bigger number densities, and therefore higher $S/N$. By quantifying this error, one could also think about correcting it in the first place. In the shot noise limit, the information on both parameters is primarily driven by standard sirens. The joint error on $h$ is $0.21\%$ for flat cosmology, degraded to $1.1\%$ with curvature. Similarly, the error on $\Omega_{m}$ is $1.1\%$, degraded to $2.4\%$.

The second contour plot, $\Omega_{m}$-$\Omega_{\Lambda}$ [see Fig. 1, panel (b)] is the primary constraint on dark energy. Here we find a situation that is very similar to $\Omega_{m}$-$h$ – as a geometry probe standard sirens drive most of the information, and the constraint is much better than that of WL. We predict a joint error on $\Omega_{\Lambda}$ of $2\%$ from $26\%$ with only WL. Our constraint on $\Omega_{K}$ is $\sim 0.02$.

The third and fourth contour plots, $\Omega_{m}$-$\sigma_{8}$ and $h$-$\sigma_{8}$ [see Fig. 1, panels (c) and (d)] – the main result of this paper – illustrate the benefit of combining standard sirens with GW-WL, which is sensitive to clustering. Here the constraint from WL alone would already be competitive with galaxy WL surveys to date, even without applying any informative prior on $h$ (apart hard bounds such as $0\leq h\leq 1$) as instead done in galaxy surveys (see e.g. Refs Abbott et al. (2018); Köhlinger et al. (2017)). By combining with standard sirens – a geometry probe with great statistical power on $h$ and $\Omega_{m}$ – the $\Omega_{m}$-$\sigma_{8}$ as well as $h$-$\sigma_{8}$ degeneracies are broken and the constraint on clustering is improved dramatically. Note that the WL constraint has been marginalised over a broad range on $h$, whereas the standard sirens one over a narrow range [compare with panels (a) and (b)]. This allows breaking the degeneracy and therefore significantly improving the constraint on $\sigma_{8}$. For instance, the error on $\sigma_{8}$ is $30\%$ from WL alone; this is reduced to $3\%$ for a joint analysis. Correlated with $\sigma_{8}$ and $\Omega_{m}$ are the $\ln(10^{10}A_{s})$ and $n_{s}$ parameters, which are respectively constrained to $2\%$ and $7\%$.

V Summary

This paper has set out the general framework with which a joint analysis of WL and standard sirens could be done in the foreseeable future, with realistic assumptions on the expected measurement errors. The benefits of this approach have been illustrated in the figures above. The improvement of the constraints on dark matter and dark energy, over a WL-only analysis, is evident and always at least an order of magnitude. The constraints from WL alone would already be competitive with ongoing surveys to date, however with a far smaller number of sources. The joint analysis with standard sirens would only bring further improvement owing to it adding more information and therefore lifting key degeneracies. For instance, one of the key findings of this paper is that, although standard sirens is significantly better than WL in constraining geometry parameters, the combination of the two observables breaks clustering parameter degeneracies – in this case, $\Omega_{m}$ and $\sigma_{8}$.

Our forecast is based on nominal realistic assumptions about measurement errors, namely luminosity distance, position error, optical redshift, and number density. We found that shot noise is the limiting factor in our WL analysis, thus not allowing easy access to non-linear scales. This could in principle be problematic because of the limited available information. On the other hand, the modelling is generally easier on larger scales, and therefore less prone to potential systematics.

A standard sirens analysis does not require distance ladder calibration – this dramatically reduces the need of external calibration. However, given the uncertainties related to the redshift determination, our joint analysis would probably suffer from systematics that are different from other probes. Although an accurate assessment of all the systematics and their impact to cosmology is not currently available to date, the benefit of our approach is obviously clear: the joint analysis of standard sirens and GW-WL might help solve the tensions in the $H_{0}$ and $\Omega_{m}$-$\sigma_{8}$ spaces between the various cosmological probes, such as CMB, galaxy WL, and Type Ia SNs, and distinguish between residual systematics and new physics. We conclude that it is not unrealistic to expect that our joint analysis will probably compete with (if not outperform with increasing number densities) cosmology experiments of the future.

Acknowledgements.

GC acknowledges the following people for fruitful discussion, interaction, and feedback at different stages of this research: S. Alam, M. Cataneo, J. Gair, A. Hall, C.-A. Lin, L. Lombriser, J. Peacock, J. H. Deraps, and all the WL lunch attendees (Edinburgh); E. N. Chisari (Oxford); A. Sesana (Birmingham). AT thanks the Royal Society for the support of a Wolfson Research Merit Award.