# Tomographic weak lensing bispectrum: a thorough analysis towards the next generation of galaxy surveys

Matteo Rizzato, Karim Benabed, Francis Bernardeau, Fabien Lacasa

► To cite this version:

Matteo Rizzato, Karim Benabed, Francis Bernardeau, Fabien Lacasa. Tomographic weak lensing bispectrum: a thorough analysis towards the next generation of galaxy surveys. Monthly Notices of the Royal Astronomical Society, 2019, 490, pp.4688-4714. 10.1093/mnras/stz2862. hal-01975051

HAL Id: hal-01975051

https://hal.science/hal-01975051v1

Submitted on 6 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Tomographic weak lensing bispectrum: a thorough analysis towards the next generation of galaxy surveys

Matteo Rizzato,$^{1,2*}$ Karim Benabed,$^{1,2}$ Francis Bernardeau$^{1,2,3}$ and Fabien Lacasa$^4$

Accepted 2019 October 6. Received 2019 October 2; in original form 2019 February 22

## ABSTRACT

We address key points for an efficient implementation of likelihood codes for modern weak lensing large-scale structure surveys. Specifically, we focus on the joint weak lensing convergence power spectrum–bispectrum probe and we tackle the numerical challenges required by a realistic analysis. Under the assumption of (multivariate) Gaussian likelihoods, we have developed a high performance code that allows highly parallelized prediction of the binned tomographic observables and of their joint non-Gaussian covariance matrix accounting for terms up to the six-point correlation function and supersample effects. This performance allows us to qualitatively address several interesting scientific questions. We find that the bispectrum provides an improvement in terms of signal-to-noise ratio (S/N) of about 10 per cent on top of the power spectrum, making it a non-negligible source of information for future surveys. Furthermore, we are capable to test the impact of theoretical uncertainties in the halo model used to build our observables; with presently allowed variations we conclude that the impact is negligible on the S/N. Finally, we consider data compression possibilities to optimize future analyses of the weak lensing bispectrum. We find that, ignoring systematics, five equipopulated redshift bins are enough to recover the information content of a Euclid-like survey, with negligible improvement when increasing to 10 bins. We also explore principal component analysis and dependence on the triangle shapes as ways to reduce the numerical complexity of the problem.

**Key words:** gravitational lensing: weak–methods: analytical–methods: statistical–cosmology: theory–large-scale structure of Universe.

## 1 INTRODUCTION

In the coming decades, large-scale structure (LSS) surveys will probe the late time Universe with an unprecedented precision, allowing us to constrain dark energy models. Euclid, the Dark Energy Survey (DES), and the Large Synoptic Survey Telescope (LSST) are well known examples of observational programs aiming at this scientific goal (The Dark Energy Survey Collaboration 2005; LSST Science Collaboration et al. 2009; Laureijs et al. 2011). One of their key probes will be the detection of the weak lensing signal from distant sources. Unlike galaxy clustering measurements, weak lensing probes the total amount of matter in the Universe (Bartelmann & Schneider 2001; Kilbinger 2015) giving us unbiased information on its distribution and evolution over time. In preparation for future missions, we need to understand the performance of this probe in terms of cosmological parameter error forecasts, crucially for deviations away from the ΛCDM scenario.

A standard way to extract cosmological information from the matter field, is through the computation of its 2-point correlation function and associated errors. So far, weak lensing analyses have employed this basic approach, the field of interest being the projection of the matter field along the line-of-sight. At low redshift, the non-linear evolution of the LSS of the Universe skews the distribution of the matter field and the 2-point correlation function is not sufficient for a complete description of the statistical properties of the field. In order to recover the cosmological information not accessible via standard approaches, different techniques have been advocated (Biswas, Alizadeh & Wandelt 2010; Kratochvil, Haiman & May 2010; Carron & Szapudi 2013; Pisani et al. 2015; Barthelemy et al. 2019). In this work, we exploit the

<footnote>$^1$Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote>$^2$Sorbonne Universités, Institut Lagrange de Paris (ILP), 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote>$^3$Institut de Physique Théorique, Université Paris-Saclay, CEA, CNRS, UMR3681, F-91191 Gif-sur-Yvette, France</footnote>

<footnote>$^4$Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ernest Ansermet, CH-1211 Geneva, Switzerland</footnote>

<footnote>* E-mail: matteo.rizzato@iap.fr</footnote>
higher order correlation functions of the weak lensing field to catch part of the missing cosmological information. We rely on the assumption that their joint statistical distribution is a multivariate Gaussian.

The constraining power of this approach can be addressed either via Fisher forecast (Tegmark 1997), DALI forecast (Sellentin, Quartin & Amendola 2014; Sellentin 2015), or Markov Chain Monte Carlo (MCMC) analyses (Christensen et al. 2001; Lewis & Bridle 2002; Dunkley et al. 2005; Akeret et al. 2013). The core part of these approaches lies in the computation (and inversion) of covariance matrices to quantify the correlations between different binned modes of the field. As for weak lensing, modern galaxy surveys will be capable to detect the position of the sources in different tomographic redshift bins probing the time evolution of the matter field. The complexity of our covariance matrices is enhanced by the correlations between observables sourced by galaxies in different bins: with the final set-up considered in our computations, we compute covariance matrices with $\gtrsim 10^4 \times 10^4$ elements. Future data analyses based on the joint weak lensing power spectrum–bispectrum probe requires the manipulation of these matrices. Optimization via the reduction of their dimensionality is an extremely important path to explore (Heavens et al. 2017). Bearing this idea in mind, we present possible summary statistics to be applied to both the power spectrum and the bispectrum and we explore possible paths to further simplify the complexity of future forecasts. In order to make these analyses feasible (at least for a small cluster), an original algorithm had to be developed: the code will shortly be made available publicly, together with its documentation.

The weak lensing signal is sensitive to very small scales. For a Euclid-like survey, Kitching & Taylor (2011) proved that the power spectrum has to be accurately known to 1 per cent down to $k \approx 50\,h\,\text{Mpc}^{-1}$ to saturate the dark energy figure of merit. This study is in agreement with the previous one from Doré, Lu & Pen (2009). Similarly, Huterer & Takada (2005) argued that we typically need the power spectrum within a few per cent accuracy at the $k \approx 10\,h\,\text{Mpc}^{-1}$ in order for the uncertainty on the power spectrum itself not to degrade cosmological constraints in DES- and LSST-like survey. In Eifler (2011), Hearin, Zentner & Ma (2012) a similar conclusion is obtained. At these scales linear theory for the evolution of the matter perturbations cannot be trusted to build our observables. Perturbative approaches (e.g. perturbation theory up to 2 loops) starts deviating by more than 1 per cent at scales $k \gtrsim 0.1\,h\,\text{Mpc}^{-1}$ at $z=0$ (Taruya et al. 2012). In the above regime cosmological N-body simulations are usually employed to study the non-linear gravitational evolution and eventually used to tune phenomenological models or fitting formulae for the power spectrum. For example, Peacock & Dodds (1996) provided a fitting formula for the power spectrum based on a scaling ansatz presented in Hamilton et al. (1991). Later, Smith et al. (2003) proposed a new model of the power spectrum, the so-called halofit model, which is based on a well known halo model of structure formation (e.g. Ma & Fry 2000; Seljak 2000; Cooray & Hu 2001). Mead et al. (2015) presented an optimized variant of the halo model, designed to produce accurate matter power spectra well into the non-linear regime for a wide range of cosmological models, including baryonic feedback. The halo model, along with its variants, has nowadays become a standard in the computation of covariance matrices within galaxy survey pipelines (Hildebrandt et al. 2017; Krause & Eifler 2017). All the approaches mentioned above are meant to provide a good fit for the power spectrum alone. However, they give a poor description of the true underlying physics (Reimberg et al. 2018). In particular, they cannot be used to compute higher order statistics, which are key ingredients in our analysis. Here, we rather build our observables in the standard halo model framework. Recently, Kayo, Takada & Jain (2013) proved that the halo model has a precision of 20 per cent up to scales $k \approx 10\,h\,\text{Mpc}^{-1}$ at $z=0$ in terms of weak lensing convergence power spectrum and bispectrum and of the joint covariance. Though not sufficient for actual data analysis, this precision is sufficient for our purpose of an S/N evaluation. While we perform several approximations to optimize the computational resources at our disposal, the mentioned model accuracy is actually the real bottleneck of our implementation. Even if promising, the halo model is not the final theory for dark matter clustering: its precision might start degrading before reaching the desired scales, assumptions have to be used to build the observables and different effects can be added or neglected according to their importance for the forecast. Thanks to the high performance of the tool developed within the present project, we test the robustness of our forecast against one of the main uncertainty of the model. We underline that the main goal of this work is the validation and the test of the algorithm developed for the computation and compression of the joint weak lensing convergence power spectrum-bispectrum covariance matrix, being a challenging numerical task by itself. We leave sophisticated parameter forecasts (which anyway require the evaluation of the covariance) for future works. Instead we use here a signal-to-noise ratio (S/N) analysis as suggested by Tegmark, Taylor & Heavens (1997), which has become a standard tool in the literature to gauge the information content of weak lensing observables (Rimes & Hamilton 2005; Sato et al. 2009; Takada & Jain 2009; Kayo et al. 2013).

This paper is organized as follows. In Section 2, we introduce our notation for the theory of weak lensing. In Section 3 we provide the expressions for the observables used in this work (power spectrum and bispectrum of weak lensing convergence) and for their covariance matrix. We also list further parameters related to the modelling of the latter. The expressions provided in this section are meant to be generic: no use of any specific physical model is required. In Section 4 we give a quick review of the halo model and we build the 3D spectra in this frame: at this point, all the assumptions and the analytical expressions required for our implementation are in place. Sections 5 and 6 contain the core part of this paper. After introducing the joint S/N analysis for the power spectrum and the bispectrum, we derive our results and analyse the possibility for future data compression via the principal component analysis (PCA) of our covariance matrices. Finally, in Section 7 we summarize our main findings and the approximations underlying our analyses.

The tools we developed are extremely flexible and can be applied to different experimental setup. In this work, we derive our results according to the specificities of the ESA Euclid mission (Laureijs et al. 2011), listed in Appendix A. We build our observables in flat-sky exploiting the Limber approximation when moving to the statistics of the fields projected along the line of sight. The good performance of this approximation has been well investigated at the level of power spectra for weak lensing in Kilbinger et al. (2017): we assume that this result is reliable also at the level of higher order correlation functions given that we are interested in the information content at small
scales. Verifying these approximations, however, is an interesting path to follow in future works. We assume a six-parameter spatially flat ΛCDM with $\Omega_m = 0.238$, $\Omega_b = 0.042$, $\Omega_\Lambda = 0.762$, $h = 0.732$, $\sigma_8 = 0.76$, $n_s = 0.958$, and we work in units $c = 1$. Lengths and masses are respectively measured in Mpc h⁻¹ and Mpc h⁻¹.

## 2 WEAK LENSING

The images of far sources are distorted due to the underlying matter distribution. The matter fractional density field $\delta$ locally perturbs the background metric through the Weyl scalar potential. Given this perturbed metric, we can solve the geodesic equation for a photon emitted in a given angular direction $\theta$ at a specific time, and evaluate its angular deviation from an unperturbed trajectory at every moment. This deviation is directly related to the first derivative of the post-Newtonian gravitational potentials. The magnification of the sources in a redshift bin $b_i = [z_i, z_{i+1}]$ induced by the cosmological weak lensing in a given angular direction can be estimated via the tomographic convergence field $\kappa_{(i)}(\theta)$, which is in particular related to the Laplacian of the potential. The full calculation proves that we can express this field as the line-of-sight integration of the matter density contrast, properly convoluted with a geometrical kernel which accounts for the position of the sources (Bartelmann & Schneider 2001; Kilbinger 2015). Precisely, for sources selected in a redshift bin $b_i$, the convergence reads

$$ \kappa_{(i)}(\theta) = \int_{0}^{\chi(z_{i+1})} d\chi W_{(i)}(\chi) \delta[\theta, \chi], \qquad (1) $$

where $\chi(z)$ is the comoving distance from the observer and $H(z)$ is the time-dependent Hubble factor. The function $W_{(i)}(\chi)$ is the aforementioned lensing kernel associated to the projected distribution of the sources placed within the $i$th redshift bin. We express this last quantity as the convolution of a cosmology (only)-dependent function $S(z_s, z)$ and a function $\mathcal{F}(z_s, i)$ which is dependent on the properties of the survey

$$ W_{(i)}(\chi(z)) = \int_{z}^{\infty} dz_s \mathcal{F}(z_s, i) S(z_s, z). \qquad (2) $$

In equation (2), $z$ is the redshift that enters the line-of-sight integration in equation (1) and $z_s$ is one of the sources whose contribution to the signal is considered. The function $S(z_s, z)$ is also known as lensing efficiency

$$ S(z_s, z) = \frac{3}{2} \Omega_{m,0} H_0^2 (1+z) \chi(z) \frac{\chi(z_s) - \chi(z)}{\chi(z_s)}, \qquad (3) $$

while, $\mathcal{F}(z_s, i)$ is simply the expected projected number density of the sources in the $i$th bin $n_i(z_s)$ and it is vanishing outside the considered bin. In a more general approach, this last function should also account for the photometric errors in the detection of the source positions. Given the complexity of the computation we perform, we omit this contribution assuming we are capable of measuring the realistic position of the sources without any errors.

## 3 CORRELATION FUNCTIONS AND COVARIANCE ESTIMATION

### 3.1 Higher order correlation functions for weak lensing

Convergence maps exhibit large-scale correlations that reflect those of the matter field. One of the main results of this work consists in the analysis of the level of correlation between power spectra and bispectra for the tomographic weak lensing convergence over different scales and spatial configurations. In Section 3.2 and in Appendix B, we show that we need the correlation functions of the convergence field up to the six-point one. Indeed, the non-linear gravitational evolution couples different modes of the matter field leading to non-vanishing higher order correlation functions. In Fourier space, we define the connected part of the $n$th order correlation function of the tomographic convergence field (1)

$$ \langle \kappa_{\ell_1}^{(i_1)} \dots \kappa_{\ell_n}^{(i_n)} \rangle_c = (2\pi)^2 P_{i_1, \dots, i_n}(\ell_1, \dots, \ell_n) \delta_D(\ell_1 + \dots + \ell_n), \qquad (4) $$

where $\kappa_{\ell}^{(i)} = \int d\theta \kappa_{(i)}(\theta) e^{-i\ell\theta}$, is the Fourier transform of the convergence field (1) and $\delta_D$ is the Dirac delta. We call the quantity $P_{i_1, \dots, i_n}(\ell_1, \dots, \ell_n)$ polyspectrum of order $n$. In equation (4) we described the field on a flat sky, i.e. we approximate the full spherical harmonics decomposition of the real field with a simple 2D Fourier transform. Such approximation is valid at percent level for $\ell > 100$ (Kilbinger et al. 2017). By replacing the convergence field definition into equation (4), we can derive the expression for the flat-sky polyspectra. In general, the redshift integration (line-of-sight integration) appearing in equation (1) would naturally translate into a complex $n$-dimensional one. To simplify this calculation, we make use of the Limber approximation (Loverde & Afshordi 2008): we assume that the 3D matter polyspectra have a weak dependence on the momenta component corresponding to the line-of-sight direction. Consequently, the projection collapses into a simple 1D redshift integration and we can relate the angular multipoles $\ell$ to the 3D momenta $k$ via the well known Limber relation $k(\ell, z) \approx \{\ell/\chi(z), 0\}$. Finally, the general $n$-order tomographic convergence polyspectrum $P_{i_1, \dots, i_n}(\ell_1, \dots, \ell_n)$ relates to the same
order matter one $P(k_1, \dots, k_n)$ via

$$ 
\begin{aligned}
P_{i_1 \dots i_n}^{(n)}(\ell_1, \dots, \ell_n) &= \int_0^\infty d\chi \chi^{2-2n} \left[ \prod_{i=1}^{i_n} W_{(i)}(\chi) \right] P(k(\ell_1, \chi), \dots, k(\ell_n, \chi)) \\
&= \int_0^\infty d\chi \mathcal{T}(i_1, \dots, i_n; \chi) P(k(\ell_1, \chi), \dots, k(\ell_n, \chi)).
\end{aligned} 
\quad (5) 
$$

For consistency with the literature, we call the 2-, the 3-, and the 4-order polyspectrum, respectively, power spectrum, bispectrum, and trispectrum

$$ P_{ij}(\ell) \equiv P_{ij}^{(2)}(\ell), \quad (6) $$

$$ B_{ijk}(\ell_1, \ell_2, \ell_3) \equiv P_{ijk}^{(3)}(\ell_1, \ell_2, \ell_3), \quad (7) $$

$$ T_{ijkl}(\ell_1, \ell_2, \ell_3, \ell_4) \equiv P_{ijkl}^{(4)}(\ell_1, \ell_2, \ell_3, \ell_4). \quad (8) $$

We underline that the assumption of an isotropic and homogeneous Universe (Cosmological Principle) allows us to reduce the actual dependencies of the polyspectra. The power spectrum is expressed as function of the module $\ell$ of the momentum $\ell$ (we define $\ell \equiv |\ell|$) and the bispectrum has a dependence on just three degrees of freedom, e.g. the edges of the associated triangular configuration (Marcori & Pereira 2017). The power spectrum measurements for an actual survey are also affected by intrinsic shape noise due to the finite number of sources and the intrinsic variability of galaxy shapes. Assuming that the orientation of intrinsic galaxy shapes is random and the shapes of different galaxies are uncorrelated, this component is Gaussian<fnref n="1" /> and we account for this effect in the following way

$$ P_{ij}^{s}(\ell) \equiv P_{ij}(\ell) + \frac{\sigma_{\epsilon}^{2}}{\bar{n}_{(i)}} \delta_{ij}^{K}. \quad (9) $$

In equation (9), $\bar{n}_{(i)}$ is the expected projected number of sources per unit of solid angle within the $i$th redshift bin and $\delta_{ij}^K$ is the Kronecker delta for the two indices $i, j$. In the following numerical calculations, we use the value $\sigma_\epsilon = 0.3$ which is representative of the expected Euclid sample (Amendola et al. 2013) and we ignore weak lensing systematics like intrinsic alignments (Schneider & Bridle 2010; Joachimi et al. 2013; Kiessling et al. 2015; Troxel & Ishak 2015; Blazek et al. 2017; Schmitz et al. 2018), photometric redshift errors (Ma, Hu & Huterer 2006), blending (Hartlap et al. 2011), and point spread function (Jarvis & Jain 2008; Chang et al. 2013; Hamana et al. 2013).

## 3.2 Covariance matrix for the observables

### 3.2.1 Structure of the covariance matrix

In this work, we want to address the cosmological information content of the joint tomographic convergence power spectrum–bispectrum probe. We accomplish this task via a signal-to-noise ratio analysis (described in Section 5) for which a data vector of estimators and the corresponding covariance matrix is required. For our computation, we chose the binned estimators $\hat{P}_{ij}^W(\ell^b)$ (21) and $\hat{B}_{ijk}^W(\ell_1^b, \ell_2^b, \ell_3^b)$ (B1) for the true underlying tomographic power spectrum (6) and bispectrum (7) respectively. We bin the magnitude $\ell$ of the angular momenta $\ell$ in equally spaced bins in log $\ell$. At the level of notation, we define the bin $\ell^b$ as the symmetric interval of width $\Delta\ell^b$ around the central value $\ell$. Therefore, a multipole $\ell'$ belongs to the bin $\ell^b$ when its magnitude $\ell' \in [\ell - \Delta\ell^b/2, \ell + \Delta\ell^b/2]$. We refer to Appendix A for a detailed description of the chosen binning. Moving to the covariance of the estimators, we split it in the following way (Cooray & Hu 2001; Takada & Bridle 2007; Sato et al. 2009; Takada & Jain 2009; Kayo et al. 2013; Kayo & Takada 2013)

$$ \mathrm{Cov}[\hat{P}_{ij}^W(\ell^b), \hat{P}_{i'j'}^W(\ell'^b)] = \mathrm{Cov}[\dots]_{\mathrm{Gauss}} + \mathrm{Cov}[\dots]_{\mathrm{NGins}} + \mathrm{Cov}[\dots]_{\mathrm{NGssc}}, \quad (10) $$

$$ \mathrm{Cov}[\hat{B}_{ijk}^W(\ell_1^b, \ell_2^b, \ell_3^b), \hat{B}_{i'j'k'}^W(\ell_1'^b, \ell_2'^b, \ell_3'^b)] = \mathrm{Cov}[\dots]_{\mathrm{Gauss}} + \mathrm{Cov}[\dots]_{\mathrm{NGins}} + \mathrm{Cov}[\dots]_{\mathrm{NGssc}}, \quad (11) $$

$$ \mathrm{Cov}[\hat{P}_{ij}^W(\ell^b), \hat{B}_{i'j'k'}^W(\ell_1'^b, \ell_2'^b, \ell_3'^b)] = \mathrm{Cov}[\dots]_{\mathrm{NGins}} + \mathrm{Cov}[\dots]_{\mathrm{NGssc}}. \quad (12) $$

In equations (10) and (12), we label with the subscript Gauss the covariance terms containing only two-point statistics, which are non-vanishing only for correlations within the same $\ell$-bin. The other covariance terms arise due to the non-Gaussian statistics of the convergence field and correlate modes in different $\ell$-bins and the different probes, i.e. power spectrum and bispectrum. We distinguish two classes of terms, respectively, labelled via the subscripts NGins and NGssc. The former is sourced by correlations between observed intrasurvey modes, while the latter is sourced by correlations between observed modes and background supersurvey modes and it is known in the literature as supersample covariance (SSC) (Hamilton, Rimes & Scoccimarro 2006; Takada & Jain 2009; Takada & Hu 2013; Barreira, Krause & Schmidt 2018; Chan, Moradinezhad Dizgah & Noreña 2018; Barreira 2019). For notation purpose, in Appendix B we will further split the

<footnote n="1">¹A more realistic case would have a binomial shot noise, leading to a non-vanishing bispectrum component. We leave this consideration to future works.</footnote>
NGins-like terms in different components according to the order of the correlations sourcing them. Also, we will generically dub as NG the total non-Gaussian covariance, regardless of the type of correlations involved. In the following section we will briefly review the main logical and mathematical steps needed for deriving the covariance matrix for the power spectrum (10), in all its components. As not to break the flow of the paper, we defer to Appendix B a similar analysis for the bispectrum and the power spectrum-bispectrum cross-covariance.

### 3.2.2 Power spectrum covariance matrix: Gauss and NGins terms

We give here few details about the computation of the power spectrum covariance while deferring the discussion about the bispectrum and power spectrum-bispectrum cross-covariance to Appendix B. We start by defining the binned estimator for tomographic power spectrum

$$ \hat{P}_{ij}(\ell^b) = \frac{1}{\Omega_{\text{sky}} N(\ell^b)} \sum_{\ell'} \kappa_{\ell'}^{(i)} \kappa_{-\ell'}^{(j)} \Delta_{\ell b}^{(2)}(\ell'). \quad (13) $$

The sum in equation (13) runs over discrete modes which are integer multiples of the fundamental frequency of our survey $\ell_f \equiv 2\pi/\Theta_{\text{sky}}$, $\Theta_{\text{sky}}$ being the survey footprint angular size. In particular, the finite real-space domain of our observation implies that we cannot measure angular modes $\ell \le \ell_f$ and also $\ell_f$ provides a minimum resolution for the module of the multipoles we can access

$$ \ell \sim \ell_{n_x}^{n_y} = \ell_f \cdot \{n_x, n_y\}, \quad n_x, n_y \in \mathbb{N}^x. \quad (14) $$

The estimator (13) differs from the one introduced in equations (10) and (12). In particular the estimator (13) does not account for the impact of the survey mask function, whose effect is considered when introducing the estimator $\hat{P}_{ij}^W(\ell^b)$ (21). The link between the two will be clearer in the following. We choose this presentation for didactic purposes. The selection function $\Delta_{\ell b}^{(2)}(\ell') = 1$ when the magnitude $\ell'$ of the mode $\ell'$ falls into the required bin $\ell^b$. Also, the quantity $N(\ell^b) \approx 2\ell \Delta \ell^b f_{\text{sky}}$ (in the limit $\ell \gg \ell_f$) gives the number of vector pairs $\hat{l}, -\hat{l}$ whose magnitude $\hat{l}$ is within the bin $\ell^b$, each pair being discriminated by a deviation in the module of the vectors of a unit of the survey fundamental mode $\ell_f$ (Joachimi, Schneider & Eifler 2008; Joachimi, Shi & Schneider 2009; Kayo et al. 2013). In equation (13) we have introduced the survey angular coverage $\Omega_{\text{sky}}$, which is related to the linear angular size of the survey $\Theta_{\text{sky}}$ via

$$ \Omega_{\text{sky}} = 2\pi (1 - \cos \Theta_{\text{sky}}). \quad (15) $$

Often, it is useful to refer to the angular coverage in terms of the fraction $f_{\text{sky}}$ of the sky which is observed by the survey. In our case, it is defined as

$$ f_{\text{sky}} = \frac{\Omega_{\text{sky}}}{4\pi}. \quad (16) $$

The estimator (13) is proved unbiased (Takada & Bridle 2007). Also, assuming that the true power spectrum varies slowly within the bin width, it can be approximated by the power spectrum (6) itself as calculated at the central value of the corresponding bin, i.e.

$$ \langle \hat{P}_{ij}(\ell^b) \rangle \approx P_{ij}(\ell). \quad (17) $$

We can derive the covariance of the estimator (13) by applying its standard definition (Cooray & Hu 2001; Takada & Bridle 2007; Sato et al. 2009; Takada & Jain 2009; Kayo et al. 2013)

$$ \begin{aligned} \text{Cov} \left[ \hat{P}_{ij}(\ell^b), \hat{P}_{i'j'}(\ell'^b) \right] &= \langle \hat{P}_{ij}(\ell^b) \hat{P}_{i'j'}(\ell'^b) \rangle - \langle \hat{P}_{ij}(\ell^b) \rangle \langle \hat{P}_{i'j'}(\ell'^b) \rangle \\ &= \frac{1}{\Omega_{\text{sky}} N(\ell^b)} \frac{1}{\Omega_{\text{sky}} N(\ell'^b)} \sum_{\bar{\ell}, \ell'} \langle \kappa_{\bar{\ell}}^{(i)} \kappa_{-\bar{\ell}}^{(j)} \kappa_{\ell'}^{(i')} \kappa_{-\ell'}^{(j')} \rangle \Delta_{\ell b}^{(2)}(\bar{\ell}) \Delta_{\ell' b}^{(2)}(\ell') - P_{ij}(\ell) P_{i'j'}(\ell'). \end{aligned} \quad (18) $$

In particular, within the covariance, we account for the shot noise term introduced in equation (9). The Gaussian component can be easily derived via the Wick's theorem by decomposing the 4-point correlator in equation (18) into products of power spectra. A detailed calculation leads to

$$ \text{Cov} \left[ \hat{P}_{ij}(\ell^b), \hat{P}_{i'j'}(\ell'^b) \right]_{\text{Gauss}} = \frac{\delta_{\ell\ell'}^K}{N(\ell^b)} \left[ P_{ii'}^s(\ell) P_{jj'}^s(\ell) + P_{ij'}^s(\ell) P_{ji'}^s(\ell) \right]. \quad (19) $$

The connected component leads instead to the NGins term in equation (10)

$$ \text{Cov} \left[ \hat{P}_{ij}(\ell^b), \hat{P}_{i'j'}(\ell'^b) \right]_{\text{NGins}} = \frac{1}{N(\ell^b) N(\ell'^b) \Omega_{\text{sky}}} \sum_{\bar{\ell}, \ell'} \left[ T_{ij i' j'} (\bar{\ell}, -\bar{\ell}, \ell', \ell') \right] \Delta_{\ell b}^{(2)}(\bar{\ell}) \Delta_{\ell' b}^{(2)}(\ell') \approx \frac{1}{\Omega_{\text{sky}}} T_{ij i' j'}(\ell, -\ell, \ell', -\ell'). \quad (20) $$

In equation (20) the exact covariance evaluation would require an average of the trispectrum over the two bins $\ell^b, \ell'^b$. However, as we will explain in Section 4.1, we work in a regime for which the 1-halo term is a good approximation to the trispectrum. This component does not depend on the angles between the wavevectors. On top of that, we assume that the trispectrum does not vary significantly within the bins, dropping the average over the modules as well.
### 3.2.3 Power spectrum covariance matrix: NGssc term

The supersample covariance term, i.e. the second non-Gaussian component in equation (10), is rooted in the intrinsic nature of our observations. As a matter of fact, the presence (as it is the case in all observations) of a finite survey mask function $W(\theta)$ ($\tilde{W}(\ell)$ in Fourier space) of angular size $\Theta_{\text{sky}}$ induces correlations with modes which are of the order of the fundamental length of its domain in Fourier space, i.e. the fundamental frequency $\ell_{f}$. This component cannot be derived from the estimator (13) since it does not account for couplings with modes $\ell \le \ell_{f}$. We can introduce a more sophisticated one accounting for the impact of the window function $W(\theta)$ on the observed field

$$ \hat{P}_{ij}^{W}(\ell^{b}) \equiv \frac{1}{\Omega_{\text{sky}} N(\ell^{b})} \sum_{\ell'} \kappa^{W(i)}_{\ell'} \kappa^{W(j)}_{-\ell'} \Delta_{\ell b}^{(2)}(\ell'), \qquad (21) $$

where $\kappa^W$ is the Fourier transform of the observed convergence field: the convolution of the underlying convergence field $\kappa$ with the window function of the survey $\tilde{W}$

$$ \kappa_{(i)}^{W}(\theta) = W(\theta) \kappa_{(i)}(\theta), \qquad \kappa_{\ell}^{W(i)} = \int \frac{d^2\ell'}{(2\pi)^2} \tilde{W}(\ell') \kappa_{\ell-\ell'}^{(i)}. \qquad (22) $$

At the observed scales $\ell \gg \ell_{f}$, a detailed computation proves that the estimator $\hat{P}_{ij}^{W}(\ell^{b})$ is unbiased (Takada & Hu 2013; Barreira et al. 2018) so that

$$ \langle \hat{P}_{ij}^{W}(\ell^{b}) \rangle \approx P_{ij}(\ell), \qquad (23) $$

where once again we approximated the value of the estimator with the power spectrum itself on the central value of the corresponding bin (under the assumption of negligible variations within the bin). At the level of covariance we have (Takada & Hu 2013; Barreira et al. 2018; Barreira 2019)

$$ \text{Cov} \left[ \hat{P}_{ij}^{W}(\ell^{b}), \hat{P}_{i'j'}^{W}(\ell'^{b}) \right]_{\text{Gauss}} \approx \text{Cov} \left[ \hat{P}_{ij}(\ell), \hat{P}_{i'j'}(\ell'^{b}) \right]_{\text{Gauss}}, \qquad (24) $$

$$ \text{Cov} \left[ \hat{P}_{ij}^{W}(\ell^{b}), \hat{P}_{i'j'}^{W}(\ell'^{b}) \right]_{\text{NG}} \approx \frac{1}{N(\ell^{b}) N(\ell'^{b}) \Omega_{\text{sky}}} \sum_{\ell, \ell'} \int \frac{d^2 \ell''}{(2\pi)^2} |W(\ell'')|^2 T_{ij i' j'}(\ell, -\ell + \ell'', \ell', -\ell' - \ell'') \Delta_{\ell b}^{(2)}(\bar{\ell}) \Delta_{\ell' b}^{(2)}(\bar{\ell}') \qquad (25) $$

always in the limit $\ell, \ell' \gg \ell_{f}$. In equation (25), the subscript NG refer to the general non-Gaussian contribution to the covariance regardless of the type of correlations. Since we are working in the Limber approximation, the trispectrum appearing in equation (25) can be derived via the projection of the 3D matter one

$$ T_{ij i' j'}(\ell, -\ell + \ell'', \ell', -\ell' - \ell'') = \int_0^\infty d\chi \, \chi^{-6} \, \mathcal{T}(i, j, i', j'; \chi) \, T\left[\mathbf{k}(\ell, z), -\mathbf{k}(\ell, z) + \mathbf{k}(\ell'', z), \mathbf{k}(\ell', z), -\mathbf{k}(\ell', z) - \mathbf{k}(\ell'', z)\right]. \qquad (26) $$

The 3D momenta are evaluated through the Limber relation introduced around equation (5). At every redshift, the matter trispectrum can be evaluated via the consistency relations introduced by Takada & Hu (2013)

$$ T(k_1, -k_1 + p, k_2, -k_2 - p) \approx T(k_1, -k_1, k_2, -k_2) + \frac{\partial P(k_1|\delta_b)}{\partial \delta_b} \frac{\partial P(k_2|\delta_b)}{\partial \delta_b} P^{\text{lin}}(p). \qquad (27) $$

The quantities $\partial P/\partial \delta_b$ are the power spectrum responses to a change in the background matter density induced by a long background mode $\delta_b$

$$ \delta_b \equiv \frac{1}{V} \int_V d^3 \mathbf{x} \, \delta_D(\mathbf{x}). \qquad (28) $$

V being the volume accessible by the survey. Matter responses, such as those required in equation (27), can be either measured from simulations or computed within a given theoretical framework via the so-called separate Universe ansatz (Sirko 2005; Baldauf et al. 2011; Sherwin & Zaldarriaga 2012; Li, Hu & Takada 2014; Wagner et al. 2015b; Baldauf et al. 2016). We employ the second approach and we briefly go through their derivation as from the halo model in Section 4.4. After replacing equation (27) within equation (25), the first term in equation (27) easily leads to the expression for the intrasurvey covariance NGins in equation (20). The second term instead involves couplings with supersample modes $p \chi \le \ell_{f}$, commonly assumed to be in the linear regime

$$ \begin{aligned} \text{Cov} \left[ \hat{P}_{ij}^{W}(\ell^{b}), \hat{P}_{i'j'}^{W}(\ell'^{b}) \right]_{\text{NGssc}} & \approx \frac{1}{N(\ell^{b}) N(\ell'^{b}) \Omega_{\text{sky}}} \int_{0}^{\infty} d\chi \, \chi^{-6} \, \mathcal{T}(i, j, i', j'; \chi) \\ & \times \sum_{\ell, \ell'} \frac{\partial P(\mathbf{k}(\ell, \chi)|\delta_b)}{\partial \delta_b} \frac{\partial P(\mathbf{k}(\ell', \chi)|\delta_b)}{\partial \delta_b} \Delta_{\ell b}^{(2)}(\bar{\ell}) \Delta_{\ell' b}^{(2)}(\bar{\ell}') \int \frac{d^2 \ell''}{(2\pi)^2} |\tilde{W}(\ell'')|^2 P^{\text{lin}}(\mathbf{k}(\ell'', \chi)) \\ & \approx \frac{1}{\Omega_{\text{sky}}} \int_{0}^{\infty} d\chi \, \chi^{-6} \, \mathcal{T}(i, j, i', j'; \chi) \frac{\partial P(\mathbf{k}(\ell, \chi)|\delta_b)}{\partial \delta_b} \frac{\partial P(\mathbf{k}(\ell', \chi)|\delta_b)}{\partial \delta_b} \sigma_{W}^{2}(\chi). \end{aligned} \qquad (29) $$
We assumed slowly varying responses within the bins. The quantity $\sigma_W^2(\chi)$ is the time-dependent variance of the linearly evolved matter field within the observed disc-like volume at comoving distance $\chi$ from the observer

$$ \sigma_W^2(\chi) = \int \frac{d^2\ell}{(2\pi)^2} |\tilde{W}(\ell)|^2 P^{\text{lin}}(k(\ell, \chi)). \quad (30) $$

We will give a detailed expression of this term in Appendix B.

# 4 THE HALO MODEL

## 4.1 Halo model performance

The evaluation of equations (10)–(12) requires the computation of the polyspectra (4) up to the sixth order one plus the responses of the observables to the background modes in the case of the supersample terms. The halo model (Cooray & Sheth 2002) provides an effective and physically motivated ansatz to compute the different polyspectra we need. As anticipated in Section 1, this theoretical framework has already been used in the context of weak lensing power spectrum and bispectrum error estimation and we expect it to be accurate at 20 per cent up to $\ell \approx 7000$ (see Kayo et al. 2013 for a test without source tomography) with a minor dependence on the details of the implementation (see for example Fosalba, Pan & Szapudi 2005 for a discussion at the level of two- and three-point correlation function). Even if the above stated precision may not be enough for the final target of future galaxy surveys, however the halo model represents the best available choice in the literature to model the matter clustering up to the six-point correlation function from an analytical prospective. Therefore, we consider the above performance satisfactory and in particular good enough for the S/N analysis performed in Section 5. Furthermore, the halo model, along with its variants, has nowadays become a standard in the computation of covariance matrices within galaxy survey pipelines (Hildebrandt et al. 2017; Krause & Eifler 2017). This description of matter clustering relies on the fact, supported by numerical simulations, that we can model the statistical properties of the matter field via haloes of dark matter of different masses, redshifts, and positions. For the practical halo model implementation, we used the Sheth and Tormen mass function (Sheth & Tormen 1999), the NFW halo profile (Navarro, Frenk & White 1996), and the concentration–mass relation $c_v-m$ from Bullock et al. (2001). In the following we investigate the impact of the uncertainties of the parameters in the $c_v-m$ relation on our covariances. Our computational framework is flexible enough to allow for the estimation of the impact of other assumptions of the model, such as the parametrization of the mass function and the halo profile. This will be left for subsequent works. In the halo model the n-point correlation function is described as the sum of terms accounting for all the possible distributions of the points within the haloes: from the one-halo term where all the points lie within a single halo, dominant at the smallest scales, up to the n-halo term where every point lies in a different halo, dominant at the largest scales. For the data vector, we use all the n-halo terms, both for the power spectrum (one- and two-halo) and bispectrum (one-, two-, and three-halo). We show their analytical expressions in Section 4.2. We restrict the polyspectra beyond the second order one (power spectrum) within the covariance to the one-halo term only, the expressions being provided in Section 4.3 instead. In the simpler case of the bispectrum only covariance, the two- and three-halo terms contribute mainly by increasing the correlations between squeezed triangular configurations (i.e. $B(\ell_i, \ell_j, \ell_k)$ with $\ell_k \sim \ell_j \gg \ell_i$), and change the overall signal-to-noise ratio by a few per cent. Those squeezed configurations have a low information content. Beside, we verified in the joint covariance that they are strongly correlated with power spectrum modes, and thus does not contribute significantly in our final signal-to-noise figure. For this reason, we simplify our computations and assume that the one-halo approximation is sufficient for these analyses.

## 4.2 The halo model: matter polyspectra for the data vector

### 4.2.1 Halo bias

In the halo model framework the matter polyspectra are computed as a sum over all the possible point allocations within one or more haloes. While the value of the matter field is related to the halo properties via the halo density profile, the statistical properties of the contrast matter density field do not reflect directly those of the number density of haloes: matter haloes are indeed biased tracers of the underlying matter density field. If we define the halo density contrast for haloes of mass $m$ as (Mo & White 1996; Mo, Jing & White 1997)

$$ \delta_h(\mathbf{x}, t, m) = \frac{n_h(\mathbf{x}, t, m)}{\bar{n}_h(t, m)} - 1, \quad (31) $$

where $\bar{n}_h(t, m)$ is the mean comoving number density for haloes of mass $m$ and $n_h(\mathbf{x}, t, m)$ is the actual one at position $\mathbf{x}$, then $\delta_h(\mathbf{x}, t, m)$ is not equal to $\delta(\mathbf{x}, t)$ and it can be expressed as an expansion in powers of the latter (see Desjacques, Jeong & Schmidt 2018 for a thorough review on the subject). In order to describe the power spectrum and the bispectrum at leading order in perturbation theory (PT) (see Bernardeau et al. 2002 for a thorough review), we need terms up the second power in the field $\delta$. In particular, two local biases $b_1(m, t)$, $b_2(m, t)$ (Fry & Gaztanaga 1993) and a non-local bias term $b_{s2}(m, t)$ (Baldauf et al. 2012; Chan, Scoccimarro & Sheth 2012) are required

$$ \delta_h(\mathbf{k}, t, m) = b_1(m, t) \delta(\mathbf{k}, t) + \frac{b_2(m, t)}{2} \int \frac{d^3\mathbf{q}}{(2\pi)^3} \delta(\mathbf{q}, t) \delta(\mathbf{k} - \mathbf{q}, t) + \frac{b_{s2}(m, t)}{2} \int \frac{d^3\mathbf{q}}{(2\pi)^3} \delta(\mathbf{q}, t) \delta(\mathbf{k} - \mathbf{q}, t) S_2(\mathbf{q}, \mathbf{k} - \mathbf{q}), \quad (32) $$
with

$$S_2(\mathbf{k}_1, \mathbf{k}_2) = \frac{(\mathbf{k}_1 \cdot \mathbf{k}_2)^2}{k_1^2 k_2^2} - \frac{1}{3}, \quad b_{s2}(m, t) = -\frac{2}{7}(b_1(m, t) - 1). \qquad (33)$$

The local biases $b_1(m, t)$ and $b_2(m, t)$ can instead be predicted via the peak-background split approach from the halo mass function $f_m(m, t)$ (Kaiser 1984; Bardeen et al. 1986; Cole & Kaiser 1989; Mo & White 1996). In the following expressions we will omit the time dependence to simplify the notation.

### 4.2.2 Halo model matter power spectrum

Let us start from the power spectrum. It is well approximated at every redshift (a dependence that we will omit in the following) as the sum of the one-halo term and of the two-halo term which, respectively, captures the contribution given by the two points being in the same halo and in two separate ones. In Fourier space, the one-halo term is simply the product of the two Fourier-transformed profiles $u(m, c_v, k)$ (we recall we use the NFW profile), convoluted with the halo mass function

$$P^{\text{1h}}(k) = \int_{m}^{m_{\text{Max}}} dm \left(\frac{m}{\rho_{\text{com}}}\right)^2 f_m(m) \int d c_v p(c_v, m) u^2(m, c_v, k). \qquad (34)$$

In equation (34) $\rho_{\text{com}}$ is the comoving background matter density and $m^{\text{Max}}$, $m^{\text{Min}}$ specify the integration range for the mass integral, which has to be evaluated numerically. We will comment further on how to choose these values in Section 4.5. We also accounted for the uncertainty on the halo concentration parameter via the convolution with the probability density distribution $p(c_v, m)$. It gives the probability that a virialized halo of mass $m$ has a concentration parameter $c_v$. As anticipated, we employ the results from Bullock et al. (2001) who found that this distribution can be well approximated by a lognormal

$$p(c_v, m) = \frac{1}{c_v \sqrt{2\pi \sigma_{\text{ln}c_v}^2}} \exp \left[ - \frac{(\ln c_v - \ln \bar{c}_v(m))^2}{2\sigma_{\text{ln}c_v}^2} \right], \qquad (35)$$

where $\bar{c}_v(m)$ is the median concentration parameter for every redshift and mass and $\sigma_{\text{ln}c_v} = 0.18$, independent from the redshift. The two-halo term describes the correlation between two points hosted in two different haloes of mass (e.g.) $m_1$ and $m_2$. Then it depends on the halo-halo correlation function. Qualitatively we write

$$\langle \delta_h(\mathbf{k}_1, m_1) \delta_h(\mathbf{k}_2, m_2) \rangle \approx b_1(m_1) b_1(m_2) \langle \delta(\mathbf{k}_1) \delta(\mathbf{k}_2) \rangle \approx b_1(m_1) b_1(m_2) P^{\text{lin}}(k_1) \delta_D(\mathbf{k}_1 + \mathbf{k}_2). \qquad (36)$$

In equation (36) we stopped the bias expansion at the linear level, i.e. $\sim(b_1(m_2)\delta)$, since it is a common assumption to consider the scales here involved (beyond the virialization halo radius) in the linear regime (Cooray & Hu 2001). For this reason, we neglect quadratic corrections and the power spectra are the linear ones. A more quantitative evaluation leads to (Cooray & Hu 2001)

$$P^{\text{2h}}(k) = \left[ \int_{m}^{m_{\text{Max}}} dm \, b_1(m) \frac{m}{\rho_{\text{com}}} f_m(m) \left( \int d c_v p(c_v, m) u(m, c_v, k) \right) \right]^2 P^{\text{lin}}(k) \qquad (37)$$

and the total matter power spectrum is simply the sum of the two terms derived above

$$P^{\text{HM}}(k) = P^{\text{1h}}(k) + P^{\text{2h}}(k). \qquad (38)$$

To simplify the expression for the general matter polyspectrum, we can introduce the following quantity

$$I_{\mu}^{\beta}(k_1, \dots, k_{\mu}) = \int_{m}^{m_{\text{Max}}} dm \, b_{\beta}(m) \left(\frac{m}{\rho_{\text{com}}}\right)^{\mu} f_{m}(m) \left(\int d c_{v} p(c_{v}, m) \left[\prod_{i=1}^{\mu} u(m, c_{v}, k_{i})\right]\right), \qquad (39)$$

where $b_0 \equiv 1$. Then, the matter power spectrum (38) can be written in a more synthetic way as

$$P^{\text{HM}}(k) = I_2^0(k, k) + [I_1^1(k)]^2 P^{\text{lin}}(k). \qquad (40)$$

### 4.2.3 Halo model matter bispectrum

We can write the 3D matter bispectrum as the sum of three multihalo terms corresponding to the following cases (Cooray & Hu 2001):

* $B^{\text{1h}}$, the one-halo term: all the points are within the same halo,

* $B^{\text{2h}}$, the two-halo term: two out of three points are in the same halo while the third is not,

* $B^{\text{3h}}$, the three-halo term: each point is hosted in a different halo.
In analogy with the calculations performed at the level of power spectrum, the one-halo term for the bispectrum is simply the convolution of the halo mass function with the third power of the halo profile (in Fourier space)

$$B^{\text{1h}}(k_1, k_2, k_3) = \int_{m_{\text{Min}}}^{m_{\text{Max}}} dm \left(\frac{m}{\rho_{\text{com}}}\right)^3 f_m(m) \left(\int d\mathbf{c}_v p(\mathbf{c}_v, m) \left[\prod_{i=1}^3 u(m, \mathbf{c}_v, k_i)\right]\right) = I_3^0(k_1, k_2, k_3). \quad (41)$$

The two-halo term depends qualitatively on the matter density field as in equation (36) where a quantitative derivation leads to

$$B^{\text{2h}}(k_1, k_2, k_3) = I_1^1(k_1) I_2^1(k_2, k_3) P^{\text{lin.}}(k_1) + I_1^1(k_3) I_2^1(k_1, k_2) P^{\text{lin.}}(k_3) + I_1^1(k_2) I_2^1(k_3, k_1) P^{\text{lin.}}(k_2). \quad (42)$$

Finally, the three-halo term has the following dependence

$$\langle \delta_h(\mathbf{k}_1, m_1) \delta_h(\mathbf{k}_2, m_2) \delta_h(\mathbf{k}_3, m_3) \rangle \approx \left\langle \left( b_1(m_1) \delta(\mathbf{k}_1) + \frac{b_2(m_1)}{2} \delta(\mathbf{q}) \cdot \delta(\mathbf{k}_1 - \mathbf{q}) + \frac{b_{s_2}(m_1)}{2} \delta(\mathbf{q}) \cdot_s \delta(\mathbf{k}_1 - \mathbf{q}) \right) \right. \\ \left. \times \left( b_1(m_2) \delta(\mathbf{k}_1) + \frac{b_2(m_2)}{2} \delta(\mathbf{q}) \cdot \delta(\mathbf{k}_2 - \mathbf{q}) + \frac{b_{s_2}(m_2)}{2} \delta(\mathbf{q}) \cdot_s \delta(\mathbf{k}_2 - \mathbf{q}) \right) \right. \\ \left. \times \left( b_1(m_3) \delta(\mathbf{k}_1) + \frac{b_2(m_3)}{2} \delta(\mathbf{q}) \cdot \delta(\mathbf{k}_3 - \mathbf{q}) + \frac{b_{s_2}(m_3)}{2} \delta(\mathbf{q}) \cdot_s \delta(\mathbf{k}_3 - \mathbf{q}) \right) \right\rangle, \quad (43)$$

where we wrote in a synthetic way the convolutions (32) via the operators $\cdot$ and $\cdot_s$ and $m_{1, 2, 3}$ are the masses of the three haloes hosting the points used for computing the correlation. According to the Perturbation Theory (PT) paradigm (see Bernardeau et al. 2002 for a thorough review), we can expand the field $\delta(\mathbf{k})$ at different orders

$$\delta(\mathbf{k}) = \sum_{p} \delta^{(p)}(\mathbf{k}), \quad (44)$$

each term of the expansion being proportional to the pth power of the linearly evolved initial overdensity (Goroff et al. 1986; Makino, Sasaki & Suto 1992; Jain & Bertschinger 1996). The leading non-vanishing term in equation (43) includes at least one mode at second order in PT $\delta^{(2)}$. Overall, the three-halo term at leading order is proportional to the fourth power of the initial linearly evolved contrast density field, i.e.

$$\langle \delta_h(\mathbf{k}_1, m_1) \delta_h(\mathbf{k}_2, m_2) \delta_h(\mathbf{k}_3, m_3) \rangle|_{4\text{th}} \approx b_1(m_1) b_1(m_2) b_1(m_3) \langle \delta^{(1)}(\mathbf{k}_1) \delta^{(1)}(\mathbf{k}_2) \delta^{(2)}(\mathbf{k}_3) \rangle + \text{cycles over } \{\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3\}. \quad (45)$$

The correlation (qualitatively) $\langle \delta^{(1)} \delta^{(1)} \delta^{(2)} \rangle$ can be related to the tree-level PT matter bispectrum $B^{\text{PT}}(k_1, k_2, k_3)$. By looking at equation (43), we can see that at the same order more components are present, specifically those including $b_2 \delta^2$ and $b_{s_2} \delta^2$ with $\delta$ in the linear regime, i.e. $\delta = \delta^{(1)}$. These correlations of linear modes can be reduced via the Wick theorem in products of two-point correlation functions as

$$\begin{aligned} \langle \delta_h(\mathbf{k}_1, m_1) \delta_h(\mathbf{k}_2, m_2) \delta_h(\mathbf{k}_3, m_3) \rangle|_{4\text{th}} & \approx b_1(m_1) b_1(m_2) b_{2,(s_2)}(m_3) \langle \delta^{(1)}(\mathbf{k}_1) \delta^{(1)}(\mathbf{k}_2) \delta^{(1)}(\mathbf{q}) \delta^{(1)}(\mathbf{k}_3 - \mathbf{q}) \rangle + \text{cycles over } \{\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3\} \\ & = b_1(m_1) b_1(m_2) \frac{b_{2,(s_2)}(m_3)}{2} \left[ \langle \delta^{(1)}(\mathbf{k}_1) \delta^{(1)}(\mathbf{q}) \rangle \langle \delta^{(1)}(\mathbf{k}_2) \delta^{(1)}(\mathbf{k}_3 - \mathbf{q}) \rangle \right. \\ & \qquad \left. + \langle \delta^{(1)}(\mathbf{k}_1) \delta^{(1)}(\mathbf{k}_3 - \mathbf{q}) \rangle \langle \delta^{(1)}(\mathbf{k}_2) \delta^{(1)}(\mathbf{q}) \rangle \right] + \text{cycles over } \{\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3\} \\ & \approx b_1(m_1) b_1(m_2) b_{2,(s_2)}(m_3) P^{\text{lin.}}(k_1) P^{\text{lin.}}(k_2) + \text{cycles over } \{\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3\}, \end{aligned} \quad (46)$$

where we did not consider the contributions forcing $\mathbf{k}_i = \mathbf{k}_j$, ($i, j = 1, 2, 3$) since they are associated to degenerate triangular configurations in the original bispectrum. A detailed evaluation of the mass integration over the halo profiles leads to

$$\begin{aligned} B^{\text{3h}}(k_1, k_2, k_3) = & \prod_{i=1}^{3} I_i^1(k_i) B^{\text{PT}}(k_1, k_2, k_3) + \\ & + I_1^1(k_1) I_1^1(k_2) I_1^2(k_3) P^{\text{lin.}}(k_1) P^{\text{lin.}}(k_2) + 2 \text{ terms from cycles over } \{k_1, k_2, k_3\} + \\ & + \frac{4}{7} \left[ I_1^1(k_1) I_1^1(k_2) S_2(\mathbf{k}_1, \mathbf{k}_2) P^{\text{lin.}}(k_1) P^{\text{lin.}}(k_2) (1 - I_1^1(k_3)) + 2 \text{ terms from cycles over } \{\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3\} \right]. \end{aligned} \quad (47)$$

and we can recognize the tree-level PT power spectra (linear power spectra) and bispectrum induced by the correlations obtained in equation (46) and in equation (45), respectively. The total bispectrum is then the sum of the terms derived above

$$B^{\text{HM}}(k_1, k_2, k_3) = B^{\text{1h}}(k_1, k_2, k_3) + B^{\text{2h}}(k_1, k_2, k_3) + B^{\text{3h}}(k_1, k_2, k_3). \quad (48)$$

## 4.3 The halo model: matter polyspectra for the covariance matrix

As anticipated in Section 4.1, we approximate the polyspectra within the covariance matrix (beyond the two-point correlation function) with the respective one-halo terms. We give here their general expression for an n-order polyspectrum. While in real space they would require an
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left-hand panel: Autotomographic polyspectra of different orders (from 2 to 6) evaluated on equilateral configurations. We are considering autocorrelations for the first and the tenth tomographic bin, in a Euclid-like photometry (Appendix A). Left: polyspectra employed in the computation of the covariance, i.e. power spectrum including the 1- + 2-halo term and higher order polyspectra approximated to the one-halo component. In the legend we are employing the notation $P_i^{(n)}(\ell) \equiv P_{i_1 \dots i_n}(\ell_1, \dots, \ell_n)$ with $i_1 = \dots = i_n$ and $\ell = \ell_i = \dots = \ell_n$, the second member of the equivalence being introduced in equation (4). Right-hand panel: polyspectra employed in the computation of the data vector: power spectrum including the 1-+2-halo term (we also depict the separate behaviour for each of them) and bispectrum including 1-+2- + 3-halo term (we also depict the separate behaviour for each of them). In this panel we switched to the standard notation for the bispectrum, i.e. $B = P^{(3)}$.</figcaption></figure>

$n$-dimensional convolution of halo profiles, the Fourier analysis makes it much simpler turning it into a 1D (mass) integration of the product of $n$ Fourier-transformed halo profiles, shortly

$$P^{1h}(k_1, \dots, k_n) = I_n^0(k_1, \dots, k_n). \qquad (49)$$

Also, the dependence on solely the magnitude of the vectors derives from the fact that chosen single halo profile (NFW profile) is spherical in real space. As a summary, in Fig. 1 we show the behaviour of the polyspectra required for our analyses, both at the level of covariance (left-hand panel) and at the level of data vector (right-hand panel).

## 4.4 The halo model: matter polyspectrum responses

The halo model also provides a powerful recipe for the computation of the responses in equations (29), (B10), and (B11) required for the forward modelling of the supersample covariance. We can obtain them by taking the derivative of the power spectrum (38) and of the bispectrum (48) with respect to the long mode $\delta_b$ (28). We assume that the single halo profile is not affected by the long mode while it impacts the distribution of haloes on larger scales via the bias and the mass function. Following Chan et al. (2018), we consider linear responses in $\delta_b$ and up to the linear bias (the second order ones were proved to be subleading in Chan et al. 2018 and we neglect for consistency the tidal component of the response being also quadratic in $\delta$). We can formally express the impact of the long mode with the help of equation (39)

$$\frac{\partial I_{\mu}^{\beta}(k_1, \dots, k_{\mu} | \delta_b)}{\partial \delta_b} = \int_{m_{\text{Min}}}^{m_{\text{Max}}} dm \left(\frac{m}{\rho_{\text{com}}}\right)^{\mu} \frac{\partial}{\partial \delta_b} [b_{\beta}(m) f_{\text{m}}(m)] \left( \int dc_v p(c_v, m) \left[ \prod_{i=1}^{\mu} u(m, c_v, k_i) \right] \right) = I_{\mu+1}^{\beta}(k_1, \dots, k_{\mu}), \qquad (50)$$

where the derivative appearing within the expression above can be computed via the peak-background split approach as (Press & Schechter 1974; Cole & Kaiser 1989; Mo & White 1996; Sheth & Tormen 1999; Desjacques et al. 2018)

$$b_{\beta}(m) = \frac{1}{f_{\text{m}}(m)} \frac{\partial f_{\text{m}}(m)}{\partial \delta_b}. \qquad (51)$$
Therefore, the halo model power spectrum linear response is related to the one of the linear matter power spectrum as

$$ \left. \frac{d P^{\mathrm{HM}}(k | \delta_b)}{d \delta_b} \right|_{\delta_b=0} \approx \left[ I_1^\dagger(k) \right]^2 \left. \frac{d P^{\mathrm{lin.}}(k | \delta_b)}{d \delta_b} \right|_{\delta_b=0} + I_2^\dagger(k). \qquad (52) $$

In equation (52) we neglected the term proportional to $I_2^2$ begin smaller than those already included (Chiang et al. 2014; Wagner et al. 2015a; Chan et al. 2018). Under the same approximations, the halo model bispectrum linear response can be expressed as

$$ \begin{aligned} \left. \frac{d B^{\mathrm{HM}}(k_1, k_2, k_3 | \delta_b)}{d \delta_b} \right|_{\delta_b=0} & \approx I_1^\dagger(k_1) I_1^\dagger(k_2) I_1^\dagger(k_3) \left. \frac{d B^{\mathrm{PT}}(k_1, k_2, k_3 | \delta_b)}{d \delta_b} \right|_{\delta_b=0} + \\ & + \left[ I_1^\dagger(k_1) I_2^\dagger(k_2, k_3) P^{\mathrm{lin.}}(k) + I_1^\dagger(k_1) I_2^\dagger(k_2, k_3) \left. \frac{d P^{\mathrm{lin.}}(k_1 | \delta_b)}{d \delta_b} \right|_{\delta_b=0} \right] + 2 \text{ cycles} + \\ & + I_3^\dagger(k_1, k_2, k_3). \end{aligned} \qquad (53) $$

Equations (52) and (53) require the evaluation at tree-level in perturbation theory of the linear response for the power spectrum and for the bispectrum. The derivation is quite technical and we refer to Chan et al. (2018) for a detailed step-by-step explanation

$$ \left. \frac{d P^{\mathrm{lin.}}(k | \delta_b)}{d \delta_b} \right|_{\delta_b=0} \approx \frac{47}{21} P^{\mathrm{lin.}}(k) - \frac{1}{3} \frac{d P^{\mathrm{lin.}}(k)}{d \ln k}, \qquad (54) $$

$$ \left. \frac{d B^{\mathrm{PT}}(k_1, k_2, k_3 | \delta_b)}{d \delta_b} \right|_{\delta_b=0} \approx \frac{433}{126} B^{\mathrm{PT}}(k_1, k_2, k_3) + \frac{5}{126} B_{G_2}(k_1, k_2, k_3) - \frac{1}{3} \sum_{i=1}^{3} \frac{d B^{\mathrm{PT}}(k_1, k_2, k_3)}{d \ln k_i}. \qquad (55) $$

By replacing the responses (54) and (55) into equations (52) and (53), we obtain the desired expressions.

## 4.5 Halo model numerical implementation

We would like to conclude this section by discussing more on a particular point: the mass numerical integration within equation (39). The main consequence of having a numerical integration over a finite domain is the exclusion of haloes whose masses fall outside the range ($m^{\mathrm{Min}}$, $m^{\mathrm{Max}}$). This cut-off is intrinsically present in all known halo model implementations given that they are calibrated against simulations over a finite mass range. However, we do not expect this feature to have any impacts: very light haloes, extending over extremely small scales, should not contribute to our cosmological observables while heavy haloes are exponentially suppressed due to the shape of the mass function. The sensitivity to the mass cut-off is in particular problematic for the lower bound. We can easily see this point by looking at the integrand $\mathcal{I}_\beta(m)$ in the following consistency relations (Cooray & Sheth 2002; Takada & Jain 2003)

$$ \int_{m^{\mathrm{Min}}}^{m^{\mathrm{Max}}} \frac{m}{\rho_{\mathrm{com}}} dm \, b_{\beta}(m) \, f_{\mathrm{m}}(m) \equiv \int_{m^{\mathrm{Min}}}^{m^{\mathrm{Max}}} \mathcal{I}_{\beta}(m) \, dm = \begin{cases} 1 & \text{if } \beta \le 1, \\ 0 & \text{if } \beta \ge 2, \end{cases} \qquad (56) $$

where the case for $\beta = 0$ is actually the consistency relation for the mass function. In physical terms, the relations equation (56) requires that the mass of the Universe is entirely enclosed in haloes ($\beta = 0$) of mass $m \in (m^{\mathrm{Min}}, m^{\mathrm{Max}})$ and that the overall distribution of haloes is not biased compared to the total matter distribution ($\beta > 0$), at each order in the bias expansion.

In the top left-hand panel of Fig. 2, we plot the integrand $\mathcal{I}_\beta(m)$ defined in equation (56) for $\beta = 0, 1, 2$ as function of the variable

$$ v(m, z) \equiv \left( \delta_{\mathrm{sp.}}^{c}(z) / \sigma(m) \right)^{2}. \qquad (57) $$

$\delta_{\mathrm{sp.}}^c$ and $\sigma^2(m)$ are, respectively, the barrier for the spherical collapse (Press & Schechter 1974) and the variance of the linear matter field smoothed over spheres enclosing the mass $m$. We see that the integrals we are trying to evaluate are slowly convergent for $v \to 0$ ($m \to 0$ implies $\sigma(m) \to \infty$). Therefore we would expect every halo model implementation to be very sensitive to extremely light haloes, which contradicts the physical intuition described above and consequently undermines halo models themselves, not being testable at these small masses. The problem comes from the extrapolation of the model beyond its regime of validity assuming that arbitrarily small mass haloes are present in our simulation accounting for the whole matter content and biases. Clearly equation (56) would not be satisfied since we are excluding these haloes from the integration. However, once a minimum mass is set, lighter haloes are excluded while the total matter budget may be filled by non-virialized dust. This contribution cannot be caught by the mass function $f_m(m)$, by definition. In order to fulfil the normalization constraints (56) we simply assume that the non-virialized matter content provides an effective contribution via the following regularization of the mass function and of the biases (Schmidt 2016)

$$ f_{\mathrm{m}}(m) \to f_{\mathrm{m}}(m) + \alpha_0 \delta_{\mathrm{D}}(m - m^{\mathrm{Min}}), \qquad (58) $$

$$ b_{\beta}(m) \to \begin{cases} b_{\beta}(m) & \text{if } m > m^{\mathrm{Min}}, \\ \alpha_{\beta} & \text{if } m = m^{\mathrm{Min}}. \end{cases} \qquad (59) $$

The parameters $\alpha_\beta$ ($\alpha_0$) then are fixed in order to satisfy the consistency relations (56). We do not consider any corrections related to the upper bound $m^{\mathrm{Max}}$ because heavy haloes are strongly suppressed by the mass function itself, as we can see from the top left-hand panel in Fig. 2. Our numerical integrals converge for $m^{\mathrm{Max}} \approx 10^{16} M_{\odot} h^{-1}$. In the computation of the weak lensing observables, we pushed our
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top-left: integrate $\mathcal{I}_\beta$ in equation (39) plotted against $v$ (see the main text for definition) for $\beta = 0, 1, 2$. Bottom-left: inverse halo comoving radius as a function of the virialized mass at different redshift. Right: impact on the polyspectra, from the order 2 to 6, from excluding light haloes from the mass integrations required in their computation. The minimum mass considered in the different cases is given in the legend (and included as further dependence on the spectra) and we are not differentiating polyspectra of different orders since interested in the overall behaviour. We are considering autotomographic polyspectra for the first tomographic bin evaluated on equilateral configuration $P_{1,...,1}(\ell, ..., \ell)$ (on the y-axis label we are omitting the dependence on the multipoles) and we are including a dependence on the lower integration bound $m^{\text{Min}}$ employed. We investigate the sensitivness of the polyspectra to the low-mass cut-off by computing the percentage fractional difference with a conservative case $\bar{P}_{1,...,1}(\ell, ..., \ell)$, derived by integrating over haloes of mass $m \ge 10^2 M_{\odot} h^{-1}$.</figcaption></figure>

line-of-sight integration up to $k_{\text{Max}} = 50 h \text{ Mpc}^{-1}$ (Kitching & Taylor 2011). Therefore, we do not expect haloes whose comoving radius is smaller than this scale to be significant for our analyses. In the bottom left-hand panel in Fig. 2 we show the (inverse) comoving halo virialization radius as function of their virialized mass. The radius corresponding to the above scale of $k_{\text{Max}} = 50 h \text{ Mpc}^{-1}$ (indicated by the horizontal blue dotted line) encloses a mass $m^{\text{kMax}} \approx 10^{10} M_{\odot} h^{-1}$ (with a small dependence on the redshift). Therefore, if our set-up (once applied the regularization (58)) is consistent, our polyspectra should not be sensitive to haloes lighter than $m^{\text{kMax}}$. This is indeed what we prove with the right-hand panel in Fig. 2. Here, we look at the impact on the polyspectra, from the order 2 to 6, from excluding light haloes from the mass integrations required in their computation. The minimum mass considered in the different cases is given in the legend and we are not differentiating polyspectra of different orders since interested in the overall behaviour. Also, we are presenting this analysis with single points in order to have a cleaner figure. We are considering autotomographic polyspectra for the first tomographic bin (being the ones more sensitive to light haloes) evaluated on equilateral configuration (for simplicity), i.e. $P_{1,...,1}(\ell, ..., \ell)$ (we are omitting the dependence on the multipoles on the y-axis label) according to the notation given in equation (4). We investigate the sensitivness of the polyspectra to the low-mass cut-off by comparing them with a conservative case $\bar{P}_{1,...,1}(\ell, ..., \ell)$ which is derived by integrating over haloes of mass $m \ge 10^2 M_{\odot} h^{-1}$. Our implementation is solid: we observe a first deviation by more than 1 per cent when integrating over haloes of mass $m \gtrsim 10^{10} M_{\odot}/h = m^{\text{kMax}}$. We chose for the mass integrals $m^{\text{Min}} = 10^9 M_{\odot} h^{-1}$ to avoid numerical inaccuracies with a negligible extra computational price. This is extremely interesting also because this resolution in mass can be easily achieved by present smoothed-particle hydrodynamics simulations (Genel et al. 2014; Hopkins et al. 2014; Khandai et al. 2015) allowing precise tests of the statistical properties of the halo model at an accuracy which is sufficient for the joint power spectrum–bispectrum analysis. In our specific case, we employed a halo model implementation based on the work of Bullock et al. (2001). In their work, the distribution of haloes was tested at the level of mass function with a simulation covering a mass range $10^{11} M_{\odot}/h - 10^{14} M_{\odot} h^{-1}$, from redshift $z = 40$ to the present.

# 5 ANALYSIS OF THE INFORMATION CONTENT IN THE WEAK LENSING OBSERVABLES

The major result of the work presented in this paper is the computation of remarkably large covariance matrices involved in the exploitation of the joint tomographic weak lensing convergence power spectrum–bispectrum probe and the design of an algorithm which will be a key element for future galaxy surveys. In order to deeply understand the actual benefit in terms of cosmological parameter estimation, Fisher
forecast (Tegmark 1997), DALI forecast (Sellentin et al. 2014; Sellentin 2015), or MCMC analyses are clearly advised: we leave these ambitious steps for future works. We underline that the covariance matrix computation is a key step in these directions and this work already finalized this calculation. Still, we would like to have an idea of the information content achievable via the joint probe and have insights on the level of correlations between the observables. Following Kayo & Takada (2013), we define the information content of an observable as the expected inverse variance of its amplitude A from a set of measured values x, assuming a fixed shape. In other terms, the second power of the signal-to-noise ratio (S/N) is the inverse of the unique element of the Fisher information matrix in such a single parameter experiment. Specifically, we can write

$$ \left(\frac{S}{N}\right)^2 = \frac{\partial(\mathbf{x} - A\mathbf{D})^t}{\partial A} \cdot \mathbf{Cov}^{-1} \cdot \frac{\partial(\mathbf{x} - A\mathbf{D})}{\partial A} = \mathbf{D}^t \cdot \mathbf{Cov}^{-1} \cdot \mathbf{D}, \quad (60) $$

where $A\mathbf{D} = \langle \mathbf{x} \rangle$ and $\mathbf{Cov}_{ij} = \langle (x_i - \langle x_i \rangle)(x_j - \langle x_j \rangle) \rangle$.

For a joint analysis of the binned tomographic power spectrum and bispectrum, we modify equation (60) to account for correlations up to a maximum binned measured angular multipole $\ell_{max}$. The S/N ratio as cumulative function of the maximum multipole included can be written as

$$ \left(\frac{S}{N}\right)^2_{P+B} = \sum_{\ell(i), \ell(j) < \ell_{\max}}^{i,j} D_i [C^D]_{ij}^{-1} D_j, \quad (61) $$

where we need (1) a vector of observables and (2) their inverse covariance matrices

$$ \mathbf{D} = \{\mathbf{P}, \mathbf{B}\}, \quad (62) $$

$$ \mathbf{C}^D = \begin{pmatrix} \mathbf{C}^{PP} & \mathbf{C}^{PB} \\ \mathbf{C}^{PB} & \mathbf{C}^{BB} \end{pmatrix}. \quad (63) $$

The matrices $\mathbf{C}^{PP}$, $\mathbf{C}^{BB}$, and $\mathbf{C}^{PB}$, are respectively the covariance of the power spectrum, bispectrum, and the cross-covariance between the two as defined in equations (10)–(12). The spectra are ordered within the vector $\mathbf{D}$ as suggested in Kayo & Takada (2013). Specifically, the power spectra evaluated on the $i^{th}$ $\ell$-bin are placed for increasing value of $i$. For each multipole, the tomographic indexes $(i, j)$ are ordered such that $i \le j$, $j$ being the faster varying index through the vector. The bispectra, evaluated over the binned configuration $(\ell_i, \ell_j, \ell_k)$, are placed in the order which satisfies $\ell_i \le \ell_j \le \ell_k$ where the index $j$ is the fastest and the index $k$ is the slowest varying one while moving along the vector. For each triangular configuration, the bispectra associated to the tomographic bins $(i, j, k)$ are ordered such that $k$ is the fastest index and $i$ the slowest index. While not imposing any constraint on these index for scalene triangles, we exploit symmetries at the level of triangular configurations to neglect some tomographic combinations which might eventually lead to double counting the information. We refer for more details to Kayo & Takada (2013). In both cases, tomographic spectra for the same binned $\ell$-configuration are contiguous. When interested in accessing the cosmological information in the single observable, it is just matter of assuming $\mathbf{D} = \mathbf{P}$ (resp. $\mathbf{B}$) and $\mathbf{C}^D = \mathbf{C}^{PP}$ (resp. $\mathbf{C}^{BB}$). Given that one of the main problems for future weak lensing surveys is the modelling of the non-linear scales, we exploit the S/N (1) to understand how well the overall parameter space is constrained by our observations up to a given angular scale $\ell_{max}$, (2) to test how much the uncertainties of our theoretical model can degrade the information on the cosmological parameters and (3) to understand the impact of different approximations to the covariance matrices. On top of that, motivated by the effort in simplifying the complexity of the covariance matrices, we use this quantifier to (4) verify the robustness of simpler statistics when compared with the most complete analysis. We would like to remind that the formalism presented so far is general and provide the theoretical background for forecasting the information content of the observables from any galaxy survey. The results showed in the following were specifically produced implementing the specificities for a Euclid-like survey described in Appendix A.

## 5.1 Power spectrum signal-to-noise ratio

In Fig. 3 we start our analysis from the simplest case: the power spectrum. For this first analysis we are keeping the redshift binning for the sources as in a Euclid-like tomography, i.e. 10 equi-populated bins (see Appendix A). We can see that, starting from a Gaussian covariance not contaminated by shot noise (dashed violet line), the main impact on the maximum information content is produced by accounting for the shot noise which degrades the S/N on all the scales. This is expected being a scale-independent contributions to the errors. We then compare the S/N as reconstructed from different approximations to the covariance, with and without the shot noise. We follow the notation introduced in equation (10) labelling as G, NGins, and NGssc, respectively, the Gaussian, the intrasurvey, and the supersurvey component of the covariance. In both cases, and with respect to the Gaussian approximation, the biggest loss of information is induced by the supersample covariance NGssc. Focusing on the analysis including the shot noise (the more realistic case), the NGssc induces a 40 per cent reduction on the maximum S/N. The in-survey component NGins (alone) leads to a degradation of about 30 per cent. Finally, the combined effect of the two leads to a loss of about 45 per cent.

In Fig. 4 we investigate at the level of power spectrum a first way to compress our data by looking at the behaviour of the maximum information content depending on the redshift binning of the sources. From the original 10 Euclid-like redshift bins, we test a first possible compression by re-binning the sources into respectively five (redshift bin edges: 0.001, 0.560, 0.789, 1.019, 1.324, 2.500), two (redshift bin
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left-hand panel: S/N analysis for the power spectrum of the weak lensing convergence field when considering different contributions to the covariance matrix and the full power from a 10 bin Euclid-like tomography (as specified by the superscripts on the y-axis labels). To label them, we follow the notation given in equation (10). For every case, we plot (1) the signal assuming our observations not being contaminated by shot noise (dashed lines) and (2) the signal when the shot noise is considered (solid lines). Top right: S/N analysis when the shot noise is considered, for different approximations to the covariance matrix. Bottom-right: fractional differences between the S/N (shot noise included) for a purely Gaussian covariance and the S/N as computed from different approximations to the covariance.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** S/N analysis for the power spectrum of weak lensing convergence when considering different approximations to the covariance matrix. For them, we follow the notation introduced in equation (10). Also, we allow for different binning in redshift of the expected sources. The number of bins indicated on the title to each panel are constructed as equipopulated, starting from the survey specificities listed in Appendix A. We give the edges for the different binning choices in the main text. As for the angular multipole, we keep our usual binning in $\log \ell$ (14 regularly spaced bins from $\ell = 10$ to $\ell = 5000$) and all the survey related parameters as chosen in Appendix A.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** S/N analysis for the bispectrum alone of weak lensing convergence when considering different approximations to the covariance matrix. We are here inspecting different effects. As in Fig. 4, we are looking (1) at the impact of different approximations to the covariance matrix and (2) at the behaviour of the information content as function of the binning of the sources in redshift. For the first point, we label the different components of the covariance matrix according to the notation introduced in equation (11). On top of that, we study the impact of approximating the data vector **B** by neglecting different multihalo configurations. In this analysis, we employ the same binning in *l* and in redshift as for the respective panels in Fig. 4. The shot noise is included in all the forecasts.</figcaption></figure>

edges: 0.001, 0.900, 2.500), and one (redshift bin edges: 0.001, 2.500) broader equipopulated intervals. In the case of the power spectrum analysis our data vector has dimension 770, 210, 42, and 14 when considering 10, 5, 2, or 1 tomographic bins, respectively. The edges of the bins involved are given in the caption to Fig. 4. For all the possible approximations to the covariance, we see that the maximum information content decreases (as expected) while moving to smaller numbers of bins. Also, the way the cumulative S/N is affected by the binning is independent on which covariance approximation is considered. For the most realistic case where no approximations on the covariance are assumed, the maximum S/N in the three situations proposed in Fig. 3 (from left to the right), is reduced by roughly 1 per cent, 5 per cent, and 25 per cent when compared to the 10 bin case, whose maximum value is represented by red lines within the three panels on the right. Therefore, this analysis shows that a forecast based on a five bin tomography does not spoil our knowledge on the cosmological parameters for more than 1 per cent while simplifying our computation a lot. Indeed, the vector of power spectra in the five and ten bin cases is made of 210 and 770 elements, respectively, so that the five bins case gives a ~3.7 reduction of the data vector size and a ~13.4 reduction of the covariance matrix size.

## 5.2 Bispectrum signal-to-noise ratio

In Fig. 5 we show for the bispectrum alone a similar analysis as the one proposed in Fig. 4 for the power spectrum. On top of that, we are also considering different approximations to the bispectrum in the data vector **B**, namely excluding different multihalo contributions to its description. By comparing the black and yellow solid lines (no approximations made on the data vector), we can inspect instead the impact of different approximations to the covariance matrix. Also, we investigate the possibility to compress the data vector via re-binning of the sources.

Starting from the possibility of data compression, we can see that a forecast based on five equipopulated redshift bins still allows us to recover the full information content for a Euclid-like survey. In the case of the bispectrum, the data vector **B** has dimension 72 280, 9790, 776, and 130 when considering 10, 5, 2, or 1 tomographic bins, respectively. The five bin tomography then allows us for a reduction of the size of the data vector and of the covariance matrix by a factor ~7.4 and ~54.8, respectively. The edges of the bins involved are given in the caption to Fig. 4.

Secondly, we analyse the impact on the maximum S/N of different approximations both at the level of covariance (while not having any approximations on the data vector **B** – solid lines) and at the level of data vector **B** (while keeping the same covariance – lines of the same colour). Independently of the tomographic analysis, neglecting the three-halo term leads to a very small effect on the information content at all scales (~1 per cent). The analysis is instead more sensitive to the two-halo term. The supersample covariance clearly leads to a major degradation of the information content in the bispectrum. By comparing the analyses when all the multihalo configurations are included in
<figure><img src="image_8.png" /><figcaption>**Figure 6.** S/N analysis for the weak lensing convergence power spectrum–bispectrum joint probe. We compare the analyses when including (black solid line) and excluding (yellow solid line) the NGssc component. The G and NGins components are kept in both cases. We are including all the required multihalo configurations in the bispectrum data vector **B**. As in Fig. 4, we are also focusing on the behaviour of the information content as function of the binning of the sources in redshift. We employ the same bins in *ℓ* and the same survey properties as in the three respective panels in Figs 4 and 5. In all the three panels we also depict the cumulative S/N for the power spectrum probe, when including all the components to the associated covariance (blue solid line). The green solid line represent the hypothetical case of a pure Gaussian weak lensing convergence field contaminated with shot noise.</figcaption></figure>

the bispectra (solid lines), we can see that the supersample covariance reduces the information content by ~20 per cent for the case of one bin tomography, ~15 per cent for the two and five bins cases and of ~13 per cent for the ten bins tomography.<fnref n="2" />

## 5.3 Joint signal-to-noise ratio

In Fig. 6, we finally move to the joint study of the probes. A joint analysis allows us to improve the forecast by accessing the cosmological information that was lost in mode couplings. Motivated by our previous findings, we restrict the joint analysis to five tomographic bins. In the case of the joint power spectrum–bispectrum S/N our data vector has dimension 10 000, 818, and 144 when considering 5, 2, or 1 tomographic bins, respectively. Comparing the joint S/N, NGssc included (black solid line), with the power spectrum alone case, no approximations assumed (blue solid line), we find that the maximum information content increases by ~10 per cent for all the considered tomographic analysis. This confirms the need for the inclusion of the bispectrum analysis for future weak lensing analyses. Similar to the previous paragraph, we address the effect of the supersample covariance: when included, the maximum achievable information content is reduced (compared to the yellow solid line) by about 30 per cent in the one bin tomographic case and by about 25 per cent in the two and five bin analyses. As a final remark on this joint analysis, we compare the information content in all the three tomographic cases with the hypothetical case of a Gaussian field contaminated with shot noise (dark green solid lines). If the convergence field were Gaussian, this line would represent a perfect reconstruction of the cosmological information in the field. It can be observed that, in the most informative case with five tomographic bins, we recover about 60 per cent of this ideal cosmological information.

## 5.4 Uncertainty of the theoretical model: scatter of the halo concentration parameter

We investigate the robustness of our S/N-based forecasts against the uncertainties of the halo profile properties. In particular, we want to test the impact of the convolution of the haloes with the function $p(c_v, m, z)$ (35) which gives the probability that a halo of mass *m* has a concentration parameter *c*<sub>v</sub>. Specifically it allows us to account for the scatter of haloes, as observed in simulations, around the mean relation $c_v - m$ when fitted with an NFW profile (Jing 2000; Bullock et al. 2001). In this section we make the redshift dependence explicit again. We compare the observables as built from two models. In the first model (model A), we evaluate $I_{\mu}^{\beta}$ (39) from the full integration over $c_v$. In the

<footnote n="2">²The results presented in this section represent a novelty with respect to the first version of the article. In the first version, due to a numerical error in the computation of the bispectrum supersample covariance, we erroneously detected no effect on the signal-to-noise ratio from these extra correlations. After correcting for it, we find that the supersample covariance cannot be removed from the analysis.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Fractional differences between observables (data vector and covariance) as derived from the models A and B (defined in the main text). From the top to the bottom: (1) vector of power spectra **P**, (2) vector of bispectra **B**, (3) (left) power spectrum covariance matrix, (3) (right) bispectrum covariance matrix, (4) power spectrum and joint power spectrum–bispectrum cumulative S/N. For this analysis we keep our usual binning in log *ℓ* (14 regularly spaced bins from *ℓ* = 10 to *ℓ* = 5000) and all the survey related parameter as chosen in Appendix A. Specifically, we use a five bin tomography for the sources.</figcaption></figure>

second model (model B) we simply assume $p(c_v, m, z) = \delta_D(c_v - \bar{c}_v(m, z))$. This last approximation is the one mostly used in the literature. The median value $\bar{c}_v(m, z)$ is given by the following fitting formula (Bullock et al. 2001)

$$ \bar{c}_v(m, z) = K \frac{a(z)}{a_c}, \qquad (64) $$

$$ m_*(a_c) = F m, \quad \nu(m_*, a_c) = 1, \qquad (65) $$

where $a_c$ is the epoch at which the typical collapsing mass $m_*$ (defined by $\nu = 1$) equals a fixed fraction of the halo mass $m$ at the same epoch. The best-fitting values for the free parameters $F, K$ are respectively 0.01 and 4.0 for a ΛCDM Cosmology. In Fig. 7 (starting from the top), we can see that the vector of power spectra **P** differs at maximum by 1 per cent between the two models (first panel). The data vector **B** is more affected with a peak at 5 per cent (second panel) for the smallest scale. We recall from the introduction to this section that the observables in our data vector are ordered for increasing value of the multipoles used within the configurations. The jagged profile instead corresponds to different tomographic dispositions of the sources per fixed spatial configuration. In the third row we compare the power spectrum (left-hand panel) and the bispectrum (right-hand panel) covariance matrix as derived from the two models. They differ at maximum by 4 per cent and 25 per cent, respectively. However, there is no impact at all at the level of reconstructed S/N: we can see in the bottom panel in Fig. 7 that the fractional differences between the two models is well within the numerical precision of our pipeline (~1 per cent).

For completeness, we refer to Cooray & Hu (2001) for a similar analysis. In this work, the impact of marginalizing over the concentration parameter was performed at the level of 3D, 1-halo power spectrum and trispectrum. Specifically they found an impact of ~5 per cent and ~20 per cent, respectively, for a distribution $p(c_v, m, z)$ of variance $\sigma_{\ln c_v} = 0.2$.
## 6 PRINCIPAL COMPONENT ANALYSIS

The bispectrum only brings a relatively small improvement of information, mainly due to the effect of the supersample covariance. Thus, it seems desirable to simplify our problem and restrict our data vector to the most informative modes. The PCA is a simple way to exhibit which observables, or linear combinations of them, are the most informative. It can be achieved via an eigenvalue decomposition of the $n \times n$ covariance $\mathbf{C}$ on to an orthonormal basis of vectors $\mathbf{v}_{i=1,...,n}$

$$C_{ij} = \sum_{a} S_{ai} S_{aj} \lambda_{a}. \qquad (66)$$

In the above equation, $S_{ij} \equiv \mathbf{v}_{a,j}$ and $\lambda_a$ is the $a$th eigenvalue associated to the mode $\mathbf{v}_a$. Also, we are omitting the superscript D: from now on we apply this method only to the joint covariance implying $\mathbf{C} \equiv \mathbf{C}^D$. The decomposition (66) allows us to find linear combinations of the original observables in the data vector $\{D_i\}_{i=1,...,n}$ which define a new set of uncorrelated observables $\{\hat{D}_i\}_{i=1,...,n}$

$$\hat{D}_i = \sum_j S_{ij} D_j. \qquad (67)$$

At this point, the signal-to-noise ratio for each of them is simply the ratio between the signal and its own error

$$\left(\frac{S}{N}\right)_a^2 = \frac{\hat{D}_a^2}{\lambda_a}. \qquad (68)$$

In the PCA formalism, the eigenvalue $\lambda_a$ is indeed the variance associated to $\hat{D}_a$. From a much more interesting perspective, we can analyse the information content in the different modes $\mathbf{v}_a$ of the covariance as cumulative function of the new observables

$$\left(\frac{S}{N}\right)_{c,a_{\max}}^2 \equiv \sum_{a=1}^{a_{\max}} \left(\frac{S}{N}\right)_a^2, \qquad (69)$$

up to a mode $\mathbf{v}_{a_{\max}}$. In the following, we analyse the cosmological information both as quantified from a single mode (equation 68) and as cumulative function of a set of eigenmodes (equation 69). In this kind of analyses, the smallest eigenvalues (associated to the most important eigenmodes in terms of information reconstruction) are potentially affected by numerical errors. Due to the wide dynamics of the observables used for this work, the covariance matrices have large condition numbers: the eigenvalues span over a range of about 30 orders of magnitude and the smallest ones can be affected by numerical errors if a too naive eigenmode decomposition algorithm is used. For this reason we rely on a specific high performance routine able to search for the eigenvalues (and associated eigenmodes) in a large dynamical range.<fnref n="3" />

We present the main results of the PCA analyses in Figs 8 and 9 where the eigenmodes $\mathbf{v}_a$ are, respectively, ordered for increasing value of their variance $\lambda_a$ and decreasing information content $S/N_a$. The analyses have been applied on the full joint covariance, including supersample contributions and all the multihalo configurations for the bispectrum vector. In the first row of both figures, we display the reconstruction of the signal-to-noise ratio as cumulative function of the eigenmodes included in the covariance. In the second row of both figures instead, we show the information content per eigenmode $S/N_a$. We refer to the corresponding captions for more details. From a parallel study of Figs 8 and 9 we can investigate the possibility of further reducing the dimensionality of our analysis while preserving the maximum information content. These figures ought to be analysed in the light of Fig. 10 where we show the ratio between the (original) observables within the data vector $\mathbf{D}$ and their standard deviations $\sigma_i^D \equiv \sqrt{C_{ii}^D}$ (for a two bin tomography, for simplicity). In particular, Fig. 10 would represent the information content of our observables if they were not correlated. Even though it does not provide a reliable insight on the strength of our measurements in the regime we are exploring, it can still help in understanding what are the configurations expected to be more (less) informative. Also, it helps in giving a physical understanding to the results of the PCA analyses.

First of all, by looking at the bottom row in Fig. 8, we can identify a fraction of the eigenmodes which have a very poor information content $S/N_a$. On the top row, they correspond to a red plateau, i.e. eigenmodes mainly associated to bispectra (in the original data vector $\mathbf{D}$) in the linear/mildly non-linear regime via the linear combination (67) (see caption for more details). This feature is present for all the tomographic redistributions of the sources. This result can be easily explained by looking at Fig. 10. As a matter of fact, the bispectra sourcing these modes have a very low signal-to-noise ratio since at these scales the matter field has a very weak deviation from the Gaussian statistics.

Going back to the top-row panels of Fig. 8, the eigenmodes located in blue regions are characterized by the largest variance, being on the right end of all the panels. However, they are also the most important in terms of recovering the total information content, carrying the largest signal-to-noise ratio per mode (bottom row). In particular, they are at the left end of the panels in Fig. 9. Following the procedure described in detail within the caption of Fig. 8, we find that these modes are mainly linear combinations of power spectra. Then, we can easily understand the location of these modes within both Figs 8 and 9. Power spectra have an absolute standard deviation much larger than the bispectra ($\sim 10^{-9}-10^{-13}$ and $\sim 10^{-12}-10^{-20}$, respectively), which motivates their location on the right end of Fig. 8 (we recall that, approximately, $C_{ii}^{PP} \sim P^2$). On the other hand, power spectra have a much larger signal-to-noise ratio when compared to all the other observables in the vector $\mathbf{D}$. By comparing the first and the second row in Fig. 10 we clearly see that power spectra can be measured with a much smaller statistical uncertainty when compared to the bispectra.

<footnote n="3">³Specifically we made use of Intel MKL Extended Eigensolver dfeast_syev and implemented a search of the desired eigenvalues for each order of magnitude.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** PCA analysis for the joint five tomographic bin covariance matrix with no approximations on the covariance nor on the bispectrum data vector. Top: $S/N_{c,a}$ (69) as cumulative function of the eigenmodes $v_a$. The eigenmodes $v_a$ are ordered for increasing value of the associated variance $\lambda_a$. The different colours on the background label modes which corresponds to observables $\hat{D}_a$ which are mostly combinations of power spectra (blue), bispectra in the linear/mildly non-linear regime [$\ell < 400$] (red) and bispectra in the non-linear regime [$\ell > 400$] (green). We perform this classification for a given observable $\hat{D}_a$ by looking at which are the most representative observables $D_i$ contributing via the matrix $S_{ai}$. Considering only observables $D_i$ for which the corresponding matrix element $|S_{ai}| \ge 0.05$, the mode **a** is then classified according to which of the above classes represents more than the 90 per cent of them. If there is not a specific preference, we assign the **grey** colour. Bottom: S/N ratio per eigenmode. The eigenmodes $v_a$ are ordered for increasing value of the associated variance $\lambda_a$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The analysis here proposed is the same as in Fig. 8. However, the eigenmodes are now ordered for decreasing value of the associated $S/N_a$ and we are reporting our results in terms of the fraction of the total number of eigenmodes used. This is just motivated by graphical reasons. Top: S/N as cumulative function of the eigenmodes (69). The colour code is the same as in Fig. 8. Bottom: $S/N_a$ per eigenmode (67). The vertical lines in the different panels indicate the fraction of eigenmodes required for recovering the 99 per cent of the maximum achievable S/N.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Classical signal-to-noise ratio analysis for the original data vector **D**. We display the ratio between the value of an observable, distinguishing between power spectra and bispectra, and the associated standard deviation $\sigma_i^{P(B)} \equiv \sqrt{C_{ii}^{PP(BB)}}$. We assumed a two bin tomography. For every configuration in Fourier space, the different tomographic distributions of the sources are indicated by the same background colour. *Top*: binned tomographic power spectrum. *Bottom*: binned tomographic bispectrum. In the latter case, different colours on the background correspond to triangular configurations sharing the largest edge. Therefore, according to the ordering of the data vector, within each of these bends we move from squeezed to equilateral configurations.</figcaption></figure>

Finally, the modes located in the green bands are mainly sourced by bispectra evaluated on modes deep into the non-linear regime (see caption to Fig. 8). Compared to the modes associated to the bispectra in the linear/mildly non-linear regime (red), they have a lower variance, as they are located at the left end of the plots in Fig. 8 while carrying more information: they represent a transition between the blue and the red modes in the top panels of Fig. 9 and they are crucial for improving the information carried by the first modes associated to the power spectra (blue). Once again, we can understand these dynamics with the help of Fig. 10, where the bispectrum at non-linear scales has a higher signal-to-noise ratio compared to the one at larger scales. On the other hand, the smaller variance of these configurations is due to the fact that the bispectrum signal is much weaker the more we measure it in the non-linear regime, as we can see in Fig. 1 (we recall that, approximately, $C^{BB} \sim \rho^3 + B^2$).

The vertical lines in Fig. 9 indicate the fraction of eigenmodes required to recover 99 per cent of the full information content, once the PCA modes have been reordered by decreasing value of the associated signal-to-noise ratio. The result is remarkable: the higher is the number of tomographic bins, the higher is the compression efficiency. Specifically, just $\lesssim 40$ per cent, $\lesssim 20$ per cent, and $\lesssim 10$ per cent of the modes are required, respectively, for the one, two, and five bin analyses.

In Figs 11 and 12 we visualize the absolute value of the elements within the projection matrix $S_{ai}$ (66) for different choices of $a$ (left and central panel) and the total matrix itself (right-hand panel). The index $a$ indicate a specific mode $v_a$ of the original covariance matrix. In the right-hand panels, the modes are ordered for decreasing value of the information content $S/N_a$. The index $i$ runs instead over the different elements of the data vector **D**. Therefore, we want here to visualize how much each of the original observables contributes to the newly defined ones $\hat{D}$. The colours on the background of the left and central panels have to be read as follows. Different tiny bends (alternatively white and coloured) refer to the same configurations in Fourier space, then spanning over the different tomographic contributions to it. We recall that in our data vector **D** these contributions are contiguous. In the region of the x-axis corresponding to the sub-vector **B** (on the right of the vertical violet line), different colours of the bends indicate Fourier configurations sharing the same larger multipole. In particular, while moving from left to the right within each of these macrobends, we are actually spanning bispectrum configurations from squeezed triangles to equilateral ones. In Fig. 11, we first focus on the one tomographic bin joint covariance matrix. In the left-hand panel we depict the values $|S_{ai}|$ as function of the index $i$ for the most informative modes $v_a$ recovering 90 per cent of the total information content. We can see that they are mainly sourced by the power spectra in the non-linear regime. In the central panel we show the values $|S_{ai}|$ as function of the index $i$ for the modes $v_a$ required to add an extra 5 per cent on top of the previous ones. They are mainly linear combinations of bispectra in the non-linear regime. Furthermore, the projection matrix shown in the rightmost panel is extremely sparse, peaking on just few configurations: overall, just
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Left-hand panel: values of the projection matrix $|S_{ai}|$ for different choices of $a$. For this plot, we picked the most informative modes $v_a$ which allow us to recover the 90 per cent of the total information content. We depict our results for the one tomographic bin joint covariance matrix. The vertical blue lines separate the vectors $\mathbf{P}$ (on the left) and $\mathbf{B}$ (on the right) within the vector of observables, labelled by $i$. Centre: same analysis as in the left-hand panel. Here we picked the most informative modes $v_a$ that allows us to recover an extra 5 per cent of the total information on top of the modes displayed in the left-hand panel. Right-hand panel: representation of the full projection matrix $|S_{ai}|$, $a$ being ordered for decreasing value of the information content per mode $S/N_a$. The red horizontal line separate power spectra and bispectra in the data vector.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Left-Centre: same analysis as in the left-central panels in Fig. 11 but for the five tomographic bin covariance matrix. The tiny different bends on the background include different tomographic contributions for the same Fourier configuration. Different colours on the background correspond to triangular configuration sharing the largest multipole. Therefore, according to the ordering of the data vector $\mathbf{D}$, within each of these bends we move from squeezed to equilateral configurations (more details within the main text). Right-hand panel: sparsity pattern for the full projection matrix $|S_{ai}|$ with $a$ being ordered for decreasing value of the information content per mode $S/N_a$. The points have been selected whenever $|S_{ai}| \ge 0.05$.</figcaption></figure>

53 observables (36 per cent) are involved in recovering 95 per cent of the information content.<fnref n="4" /> This opens a possibility for a second level of data compression beyond the one identified at the previous step at the level of modes. In Fig. 12 we perform the same analysis for the five bin tomography, reaching similar conclusions. In particular the tomography does not affect the possibility of dimensional reduction: we can discard 2093 observables from the original data vector (i.e. 80 per cent of it) and still recover 95 per cent of the information<fnref n="4" />.

In Fig. 13, we show the profiles $|S_{ai}|$ zooming on to smaller fraction of the data vector $\mathbf{D}$ (top row) and on to single triangular configurations (bottom row). We are here considering the five bin tomography covariance. The first panel on the left in the top row focus on the region of x-axis associated to the vector $\mathbf{P}$. These zoomed-in panels are useful in order to understand the distribution of the cosmological information among different tomographic bins, per fixed spatial configuration. While it is not possible to identify dominant tomographic configurations at the level of the power spectra (top-left most panel), at the bispectrum level the tomographic configurations that are most informative are, in general, those involving high redshift sources. We can clearly see this feature in the bottom row of Fig. 13 where we focus on single spatial configurations. The jagged profile mirrors the ordering of the tomographic distributions of the sources within the data vector. We recall that for the bispectrum, the observables associated to the tomographic bins $(i, j, k)$ are ordered such that $k$ is the fastest index and $i$ the slowest index. Then we can see a peak in the profile whenever, for a fixed value of $k > i, j, i$ and $j$ move closer to $k$. The maximum information for a given spatial configuration happens to be when all the tomographic indices have reached the maximum value corresponding to the furthest sources.

Finally, in Fig. 14, we show the profile of the projection matrix $|S_{ai}|$ for a single mode $a$, picked among the most informative ones

<footnote n="4">⁴We consider configurations $D_i$ with an associated weight $|S_{ai}| < 0.05$ as negligible. The same threshold has been applied to reproduce the sparsity pattern for the five tomographic bin covariance matrix in Fig. 12.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Same analysis as in the left and central panel in Fig. 12. In particular we are considering a five tomographic bin joint covariance matrix. In the bottom panels, we are zooming on just one triangular configuration per panel. The colour map used in these plots is the same as in the left and central panels in Fig. 12, to facilitate the comparison.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Profile of the projection matrix $|S_{\text{proj}}|$ for a specific choice of mode $v_a$ picked from the most informative ones accounting for 95 per cent of the total information content. The colour code is the same as in Fig. 12 in order to facilitate the comparison.</figcaption></figure>

recovering 95 per cent of the total information content in a five bin tomography. In particular we picked a mode mainly dominated by bispectra to inspect the possibility of reducing the computational burden for this part of the data vector. Unfortunately, we can see that there are no dominant configurations and bispectra evaluated on several different spatial configurations are required to access the information content of the mode analysed. Also, for a specific triangular configuration (single background bend), it is not possible to identify a dominant contribution from a specific distribution of the sources.

# 7 DISCUSSIONS AND CONCLUSIONS

As part of this final section, we summarize the approximations adopted throughout our analysis along with their range of validity. As we will explain, the precision of the halo model in describing the matter clustering can represent a bottleneck for the accuracy of the whole forecast. Therefore, while lacking a better semi-analytical model for the gravitational collapse in the non-linear regime, it is important to remember that the halo model itself is an approximation. To start with, in deriving the statistical properties of the projected fields we assumed the Limber and the flat-sky approximations. Kilbinger et al. (2017) thoroughly explored these approximations at the level of two-point statistics and compared their predictions against a full-sky approach for the cosmic shear field with a CFHTLenS-like galaxy distribution (Kilbinger et al. 2013). For the joint approximation under exam, they found it to be accurate to better than 10 per cent for $\ell > 3$, converging slowly to the true projection with per cent level precision at $\ell > 100$. Given that the accuracy of the halo model predicted polyspectra is in general much lower, we consider the above performance satisfactory for our work. Second, while computing the covariance matrix for the power spectrum and bispectrum we assumed that the trispectra, five- and six-order polyspectra in equations (20), (B9), and (B8) (respectively) are slowly varying within the considered $\ell$-bins. To our knowledge, no tests in the literature were performed on this matter for the observables and the scales
of interests. However, given the dynamics of the spectra in Fig. 1 (left-hand panel) we advise further tests of this assumption, especially for high *ℓ*. On the other hand, the high number of configurations required for our analyses makes the full integrations over the bin width unfeasible given the computational resources at our disposal. Still at the level of covariance matrix we approximated all the polyspectra with their one-halo component (starting from the three-point one). Once again, we are not aware of any studies on the impact of neglecting these higher halo terms. Given our discussion at the beginning of Section 4.1, we do not expect them to be important for cosmological analyses even though they may have a role on regularizing the covariance matrix reducing numerical errors at the inversion. Finally, in the computation of the supersample covariance (B14), we are again assuming the flat-sky and Limber approximations for both the intra- and supersurvey modes. At the power spectrum level, Barreira et al. (2018) analysed the impact of these approximations against a full spherical analysis of the supersurvey modes. For Euclid- and LSST-like tomographies and at the power spectrum covariance level, they report that for surveys covering less than ~5 per cent of the sky the two results agree to better than 1 per cent. However, when moving to the expected realistic coverage of $f_{sky} \approx 0.3-0.4$ the use of the flat-sky expressions results in an underestimation of the SSC contribution of about 10 per cent. Again, this performance has to be considered satisfactory for this work given that the precision bottleneck is mainly given by the accuracy of the matter clustering model.

For completeness, we also want to compare the analyses performed in this paper with similar works in the literature. The impact of the correlations between observed and supersample modes has already been addressed in the literature. Specifically, Kayo & Takada (2013) performed a preliminary study on the information content of the weak lensing bispectrum, which was further developed in Kayo et al. (2013) where three tomographic bins where considered for cosmological parameter forecast. However, they included the supersample correlations at very small scales. These terms are known in the literature as halo sample variance and they contribute to the correlations when all the points in a given configuration are inside the same halo. Formally, the halo sample variance is just a part of the components we obtained via the response approach used for our implementation. In particular they correspond to the last terms in equations (52) and (53) which are the responses of the one-halo term to the long mode $\delta_b$. In a different study, Barreira (2019) accounts for the supersample covariance in the PT response formalism by making different approximations. In particular, he considers only correlations between bispectra on squeezed configurations. This approximation results in a negligible impact on the signal-to-noise ratio when adding the supersample covariance. Our study shows that this is not the case when all the triangular configurations are included in the analysis.

Finally, we summarize the results of our work. We have presented the first full joint analysis of the information content for the convergence weak lensing power spectrum and bispectrum for a Euclid-like survey. We modelled the observables with the halo model and developed a high performance code capable to output fast and precise covariance matrices for the binned tomographic spectra. The covariance matrices calculated for this work have been used for signal-to-noise ratio analyses, but can benefit any forecast based on a Gaussian likelihood.

At the power spectrum level and including shot noise, we were capable to show that the supersample covariance is the main source of error, leading to a reduction of 40 per cent of the maximum achievable signal-to-noise ratio compared to the Gaussian case. The non-Gaussian cross-correlations between in-survey modes account for a loss of information of about 30 per cent instead when compared to the Gaussian case. The combined effect of these two sources of error leads to a loss of about 45 per cent on the signal-to-noise ratio. An important result of our analysis is the possibility to recover the cosmological information content of a Euclid-like survey by using five equipopulated tomographic redshift bins, instead of 10. This result is in particular insensitive to the angular scale and to the components included in the error budget of the observables and does not account for systematics like PSF, photo-z, blending, or intrinsic alignments. The same phenomenology was found at the bispectrum level. On the bispectrum signal-to-noise ratio the supersample covariance has an impact of about ~13 per cent (10 bin tomography) when the observables have been estimated via all the multihalo configurations. We tested the impact of these configurations in the modelling of the vector of bispectra and we found that neglecting the three-halo terms is a good approximation for our analysis, both when including or excluding the supersample covariance.

Motivated by our previous findings, we performed a joint power spectrum–bispectrum analysis on five equipopulated tomographic redshift bins. We proved that this combined approach can improve the information content by ~10 per cent with respect to the power spectrum alone. The supersample covariance of the bispectrum cannot be ignored and reduces the maximum information achievable of about ~25 per cent. In these analyses the two- and three-halo terms have been taken into account for the bispectra in the vector, but we restricted the computation of the covariance to the one-halo terms which dominates the covariance for the most relevant configurations. These are the most important results of this work.

We found that the halo modelling uncertainty due to the scatter of the concentration parameter in the simulations does not affect the joint signal-to-noise ratio by more than 1 per cent.

In this work, a preliminary study on the possibility of further compressing our covariance matrices while preserving the cosmological information is considered. In particular, by performing a principal component analysis on our covariance matrices, we found that a very small fraction of the eigenmodes ($\lesssim 10$ per cent for five bin tomography) carries most of the information, and that not all the configurations equally contribute to the full information content. Indeed, only 20 per cent of the data vector is of importance for the linear combinations that form this 10 per cent of eigenmodes, and thus a large fraction of the vector of the observables (~80 per cent for five bin tomography) is not significant to reconstruct the signal-to-noise ratio at the different scales considered in this work.

Starting from the results presented in this paper, several possible research paths are open. The natural and most important step forward is a full cosmological parameter forecast in order to translate the analyses on the signal-to-noise ratio into actual error bars on the parameters. Secondly, the expressions used to model the supersample correlations are based on the flat-sky and the Limber approximations. However, they might not be precise for the volumes accessible with future galaxy surveys and further tests are advised. Finally, the actual implementation
of data compression techniques qualitatively identified in the last part of this work would lead to great benefit for future data analyses. In particular, with regards to the numerical implementation of future likelihoods, the simplification of covariance matrices for joint two- and three-point statistics should be considered as a high priority task to face in order to be ready for upcoming data.

## ACKNOWLEDGEMENTS

We would like to express our very great appreciation to Alexandre Barreira, whose crucial suggestions were of great help in identifying an important weakness in a first version of the manuscript. We would also like to thank Sandrine Codis, Silvia Galli, Benjamin Joachimi, and Susan Pyne for fruitful discussions and interesting comments on our results. We convey our gratitude to Doogesh Kodi Ramanah for constructive suggestions to improve the manuscript. This work has made use of the Horizon Cluster hosted by Institut d'Astrophysique de Paris. We thank Stephane Rouberol for maintenance and running of this computing cluster. We made use of the public code HMcode developed by Mead et al. (2015) and available at https://github.com/alexander-mead/hmcode. We acknowledge financial support from the Institut Lagrange de Paris Laboratory of Excellence (ILP LABEX) (under reference ANR-10-LABX-63) which is financed by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02. In particular, MR acknowledges financial support from the Centre National d'Etudes Spatiales (CNES) fellowship program. FL acknowledges support by the Swiss National Science Foundation.

## REFERENCES

Akeret J., Seehars S., Amara A., Refregier A., Csillaghy A., 2013, CosmoHammer: Cosmological parameter estimation with the MCMC Hammer, record ascl:1303.003

Amendola L. et al., 2013, Living Rev. Relativ., 16, 6

Baldauf T., Seljak U., Senatore L., Zaldarriaga M., 2011, J. Cosmol. Astropart. Phys., 2011, 031

Baldauf T., Seljak U., Desjacques V., McDonald P., 2012, Phys. Rev. D, 86, 083540

Baldauf T., Seljak U., Senatore L., Zaldarriaga M., 2016, J. Cosmol. Astropart. Phys., 2016, 007

Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15

Barreira A., 2019, J. Cosmol. Astropart. Phys., 2019, 008

Barreira A., Krause E., Schmidt F., 2018, J. Cosmol. Astropart. Phys., 6, 015

Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291

Barthelemy A., Codis S., Uhlemann C., Bernardeau F., Gavazzi R., 2019, preprint (arXiv:1909.02615)

Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Phys. Rep., 367, 1

Biswas R., Alizadeh E., Wandelt B. D., 2010, Phys. Rev. D, 82, 023002

Blazek J., MacCrann N., Troxel M. A., Fang X., 2017, preprint (arXiv:1708.09247)

Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov A. V., Klypin A. A., Primack J. R., Dekel A., 2001, MNRAS, 321, 559

Carron J., Szapudi I., 2013, MNRAS, 434, 2961

Chang C. et al., 2013, MNRAS, 428, 2695

Chan K. C., Scoccimarro R., Sheth R. K., 2012, Phys. Rev. D, 85, 083509

Chan K. C., Moradinezhad Dizgah A., Noreña J., 2018, Phys. Rev. D, 97, 043532

Chiang C.-T., Wagner C., Schmidt F., Komatsu E., 2014, J. Cosmol. Astropart. Phys., 2014, 048

Christensen N., Meyer R., Knox L., Luey B., 2001, Class. Quantum Gravity, 18, 2677

Cole S., Kaiser N., 1989, MNRAS, 237, 1127

Cooray A., Hu W., 2001, ApJ, 554, 56

Cooray A., Sheth R., 2002, Phys. Rep., 372, 1

Desjacques V., Jeong D., Schmidt F., 2018, Phys. Rep., 733, 1

Doré O., Lu T., Pen U.-L., 2009, preprint (arXiv:0905.0501)

Dunkley J., Bucher M., Ferreira P. G., Moodley K., Skordis C., 2005, MNRAS, 356, 925

Eifler T., 2011, MNRAS, 418, 536

Fosalba P., Pan J., Szapudi I., 2005, ApJ, 632, 29

Fry J. N., Gaztanaga E., 1993, ApJ, 413, 447

Genel S. et al., 2014, MNRAS, 445, 175

Goroff M. H., Grinstein a., Rey S. J., Wise M. B., 1986, ApJ, 311, 6

Hamana T., Miyazaki S., Okura Y., Okamura T., Futamase T., 2013, PASJ, 65, 104

Hamilton A. J. S., Kumar P., Lu E., Matthews A., 1991, ApJ, 374, L1

Hamilton A. J. S., Rimes C. D., Scoccimarro R., 2006, MNRAS, 371, 1188

Hartlap J., Hilbert S., Schneider P., Hildebrandt H., 2011, A&A, 528, A51

Hearin A. P., Zentner A. R., Ma Z., 2012, J. Cosmol. Astropart. Phys., 4, 034

Heavens A. F., Sellentin E., de Mijolla D., Vianello A., 2017, MNRAS, 472, 4244

Hildebrandt H. et al., 2017, MNRAS, 465, 1454

Hopkins P. F., Kereš D., Ónorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Huterer D., Takada M., 2005, Astropart. Phys., 23, 369

Jain B., Bertschinger E., 1996, ApJ, 456, 43

Jarvis M., Jain B., 2008, J. Cosmol. Astropart. Phys., 2008, 003

Jing Y. P., 2000, ApJ, 535, 30

Joachimi B., Schneider P., Eifler T., 2008, A&A, 477, 43

Joachimi B., Shi X., Schneider P., 2009, A&A, 508, 1193
Joachimi B., Semboloni E., Hilbert S., Bett P. E., Hartlap J., Hoekstra H., Schneider P., 2013, MNRAS, 436, 819

Kaiser N., 1984, ApJ, 284, L9

Kayo I., Takada M., 2013, preprint (arXiv:1306.4684)

Kayo I., Takada M., Jain B., 2013, MNRAS, 429, 344

Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf C., Liu M.-S., 2015, MNRAS, 450, 1349

Kiessling A. et al., 2015, Space Sci. Rev., 193, 67

Kilbinger M., 2015, Rep. Prog. Phys., 78, 086901

Kilbinger M. et al., 2013, MNRAS, 430, 2200

Kilbinger M. et al., 2017, MNRAS, 472, 2126

Kitching T. D., Taylor A. N., 2011, MNRAS, 416, 1717

Kratochvil J. M., Haiman Z., May M., 2010, Phys. Rev. D, 81, 043519

Krause E., Eifler T., 2017, MNRAS, 470, 2100

Lacasa F., Rosenfeld R., 2016, J. Cosmol. Astropart. Phys., 8, 005

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511

Li Y., Hu W., Takada M., 2014, Phys. Rev. D, 89, 083519

Lima M., Hu W., 2007, Phys. Rev. D, 76, 123013

Loverde M., Afshordi N., 2008, Phys. Rev. D, 78, 123506

LSST Science Collaboration et al., 2009, preprint (arXiv:0912.0201)

Ma Z., Hu W., Huterer D., 2006, ApJ, 636, 21

Makino N., Sasaki M., Suto Y., 1992, Phys. Rev. D, 46, 585

Ma C.-P., Fry J. N., 2000, ApJ, 543, 503

Marcori O. H., Pereira T. S., 2017, J. Cosmol. Astropart. Phys., 2017, 032

Mead A. J., Peacock J. A., Heymans C., Joudaki S., Heavens A. F., 2015, MNRAS, 454, 1958

Mo H. J., White S. D. M., 1996, MNRAS, 282, 347

Mo H. J., Jing Y. P., White S. D. M., 1997, MNRAS, 284, 189

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Peacock J. A., Dodds S. J., 1996, MNRAS, 280, L19

Pisani A., Sutter P. M., Hamaus N., Alizadeh E., Biswas R., Wandelt B. D., Hirata C. M., 2015, Phys. Rev. D, 92, 083531

Press W. H., Schechter P., 1974, ApJ, 187, 425

Reimberg P., Bernardeau F., Nishimichi T., Rizzato M., 2018, preprint (arXiv:1811.02976)

Rimes C. D., Hamilton A. J. S., 2005, MNRAS, 360, L82

Sato M., Hamana T., Takahashi R., Takada M., Yoshida N., Matsubara T., Sugiyama N., 2009, ApJ, 701, 945

Schmidt F., 2016, Phys. Rev. D, 93, 063512

Schmitz D. M., Hirata C. M., Blazek J., Krause E., 2018, J. Cosmol. Astropart. Phys., 2018, 030

Schneider M. D., Bridle S., 2010, MNRAS, 402, 2127

Seljak U., 2000, MNRAS, 318, 203

Sellentin E., 2015, MNRAS, 453, 893

Sellentin E., Quartin M., Amendola L., 2014, MNRAS, 441, 1831

Sherwin B. D., Zaldarriaga M., 2012, Phys. Rev. D, 85, 103523

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Sirko E., 2005, ApJ, 634, 728

Smith R. E. et al., 2003, MNRAS, 341, 1311

Takada M., Bridle S., 2007, New J. Phys., 9, 446

Takada M., Hu W., 2013, Phys. Rev. D, 87, 123504

Takada M., Jain B., 2003, MNRAS, 340, 580

Takada M., Jain B., 2009, MNRAS, 395, 2065

Taruya A., Bernardeau F., Nishimichi T., Codis S., 2012, Phys. Rev. D, 86, 103528

Tegmark M., 1997, Phys. Rev. Lett., 79, 3806

Tegmark M., Taylor A. N., Heavens A. F., 1997, ApJ, 480, 22

The Dark Energy Survey Collaboration, 2005, preprint (arXiv:astro-ph/0510346)

Troxel M. A., Ishak M., 2015, Phys. Rep., 558, 1

Wagner C., Schmidt F., Chiang C.-T., Komatsu E., 2015a, J. Cosmol. Astropart. Phys., 8, 042

Wagner C., Schmidt F., Chiang C.-T., Komatsu E., 2015b, MNRAS, 448, L11

## APPENDIX A: A EUCLID-LIKE SURVEY

While analytical results in this work apply to any weak lensing survey, we illustrated them with specific application to the Euclid mission. To this end, we make use of the requirements presented in Laureijs et al. (2011). Specifically, the angular multipole range we investigate is [10, 5000], and we chose to bin it in 14 regularly spaced intervals in log $ℓ$, then with extrema (approximately): 10, 16, 24, 38, 59, 92, 143, 224, 348, 543, 847, 1320, 2058, 3208, 5000. For the sky coverage, we use $Θ_{sky} = 1.29$ rad, $Ω_{sky} = 15.000 \text{ deg}^2 \approx 4.57$ sterad, $f_{sky} = 0.36$. For the photometric properties of the survey, we use a total comoving number of observed sources of $n_{tot} = 30 \text{ gal arcmin}^{-2}$ from $z_{min} = 0.001$ up to $z_{max} = 2.500$. The distribution of the sources is $n(z) \propto (z/z_o)^2 \exp[-(z/z_o)^{1.5}]$ where $z_o = 0.9/\sqrt{2}$, and it is normalized over the observed range. The sources are then split in (up to) 10 equi-populated redshift bins with extrema: 0.001, 0.418, 0.560, 0.678, 0.789, 0.900, 1.019, 1.155, 1.324, 1.576, 2.500.
APPENDIX B: JOINT COVARIANCE MATRIX FOR THE WEAK LENSING CONVERGENCE POWER SPECTRUM AND BISPECTRUM

To complement the discussion of Section 3, we give here in detail the Gaussian and non-Gaussian error contributions to the tomographic bispectrum signal, including the cross-covariance with the power spectrum (Kayo & Takada 2013). While in Section 3.2 we followed a didactic approach by first introducing the estimator $\hat{B}_{ijk}^W(\ell_1^b, \ell_2^b, \ell_3^b)$ (13) and in a second moment the masked estimator $\hat{B}_{ij}^W(\ell^b)$ (21), here we introduce the bispectrum covariance directly with the general approach of the masked fields. Therefore, borrowing the notation from Section 3.2, we introduce the following binned bispectrum estimator

$$ \hat{B}_{ijk}^W (\ell_1^b, \ell_2^b, \ell_3^b) \equiv \frac{1}{\Omega_{\text{sky}} N_{\text{tri.}} (\ell_1^b, \ell_2^b, \ell_3^b)} \int \frac{d^2 \ell_1''}{(2\pi)^2} \int \frac{d^2 \ell_2''}{(2\pi)^2} \int \frac{d^2 \ell_3''}{(2\pi)^2} \bar{W}(\ell_1'') \bar{W}(\ell_2'') \bar{W}(\ell_3'') \\ \times \sum_{\ell'_1, \ell'_2, \ell'_3} \kappa_{\ell'_1}^{(i)} \kappa_{\ell'_2}^{(j)} \kappa_{\ell'_3}^{(k)} \Delta_{\ell_1^b, \ell_2^b, \ell_3^b}^{(3)} (\ell'_1, \ell'_2, \ell'_3). \qquad (B1) $$

The newly introduced selection function $\Delta_{\ell_1^b, \ell_2^b, \ell_3^b}^{(3)}(\ell'_1, \ell'_2, \ell'_3)$ forces on the modes $\ell'_1, \ell'_2, \ell'_3$ the constraint $\ell'_i \in \ell_i^b$ for $i=1, 2, 3$ simultaneously and further requires that $\ell'_1 + \ell'_2 + \ell'_3 = 0$. $N_{\text{tri.}}(\ell_1^b, \ell_2^b, \ell_3^b)$ normalizes the sum over the total number of independent triangular configurations built with the modes available within the combined bins (Joachimi et al. 2008, 2009; Kayo et al. 2013).

$$ N_{\text{tri.}}(\ell_1^b, \ell_2^b, \ell_3^b) = \sum_{\ell'_1, \ell'_2, \ell'_3} \Delta_{\ell_1^b, \ell_2^b, \ell_3^b}^{(3)} (\ell'_1, \ell'_2, \ell'_3) \approx \frac{2\Omega_{\text{sky}} \ell_1 \ell_2 \ell_3 \Delta_{\ell_1}^{\ell b} \Delta_{\ell_2}^{\ell b} \Delta_{\ell_3}^{\ell b}}{\sqrt{2 \ell_1^2 \ell_2^2 \ell_3^2 - \ell_1^4 - \ell_2^4 - \ell_3^4}}, \qquad (B2) $$

where the approximation assumed is $\Delta \ell_i^b \gg \ell_f$. The calculation starts from the basic definition of covariance of two estimators $O_1, O_2$

$$ \text{Cov}[O_1, O_2] = \langle O_1 O_2 \rangle - \langle O_1 \rangle \langle O_2 \rangle. \qquad (B3) $$

By replacing the power spectrum and the bispectrum estimators (21)–(B1) (resp.) within the covariance definition, we can derive the full bispectrum covariance and power spectrum–bispectrum cross-covariance. Similarly to equation (18), we reduce the correlations between the windowed instances of the convergence field into their connected components. In particular, the Cov[BB] (schematically) requires the decomposition of a six-point correlation function into the sum of its irreducible components: (2-×2-×2-), (3-×3-), (2-×4-), and six-order polyspectra (qualitatively). We label them in the following with the further subscript Gauss, BB, PT, and 6P, respectively. The Cov[PB] requires the decomposition of a five-point correlation function into the sum of its five-point connected component and into a series of power spectrum–bispectrum products. Analogously, we will label them with the further subscript 5P and PB, respectively. At the level of matter power spectrum–bispectrum covariance, Barreira (2019) proved that the window function does only impact the six- and five-order polyspectra related terms, respectively, in Cov[BB] and Cov[PB]. Since we are working under the flat-sky approximation, the conclusion drawn from the work of Barreira (2019) naturally extends to the weak lensing field (by moving from 3D Fourier integrals to 2D ones). A detailed calculation leads to the following expressions for the intrasurvey part (NGins) of the bispectrum covariance (Kayo et al. 2013; Kayo & Takada 2013)

$$ \text{Cov}[\hat{B}_{ijk}^W(\ell_1^b, \ell_2^b, \ell_3^b), \hat{B}_{i'j'k'}^W(\ell_1'^b, \ell_2'^b, \ell_3'^b)]_{\text{Gauss}} \qquad (B4) $$

$$ = \frac{\Omega_{\text{sky}}}{N_{\text{tri.}}(\ell_1^b, \ell_2^b, \ell_3^b)} \left[ P_{ii'}^{\text{s.}}(\ell_1) \delta_{\ell_1 \ell_1'}^{\text{K}} \left\{ P_{jj'}(\ell_2) P_{kk'}^{\text{s.}}(\ell_3) \delta_{\ell_2 \ell_2'}^{\text{K}} \delta_{\ell_3 \ell_3'}^{\text{K}} + P_{jk'}^{\text{s.}}(\ell_2) P_{kj'}^{\text{s.}}(\ell_3) \delta_{\ell_2 \ell_3'}^{\text{K}} \delta_{\ell_3 \ell_2'}^{\text{K}} \right\} \right. \\ \left. +2 \text{ terms obtained from perm. of } (i' \leftrightarrow j', \ell_1' \leftrightarrow \ell_2') \right. \\ \left. +2 \text{ terms obtained from perm. of } (i' \leftrightarrow k', \ell_1' \leftrightarrow \ell_3') \right], \qquad (B5) $$

$$ \text{Cov}[\hat{B}_{ijk}^W(\ell_1^b, \ell_2^b, \ell_3^b), \hat{B}_{i'j'k'}^W(\ell_1'^b, \ell_2'^b, \ell_3'^b)]_{\text{NGins, BB+PT}} \\ = \frac{2\pi}{\Omega_{\text{sky}}} \left[ \frac{1}{\ell_1 \Delta \ell_1^b} \left( B_{i'jk} (\ell_1', \ell_2, \ell_3) B_{ij'k'} (\ell_1, \ell_2', \ell_3') \delta_{\ell_1 \ell_1'}^{\text{K}} + B_{j'jk} (\ell_2', \ell_2, \ell_3) B_{i'ik'} (\ell_1', \ell_1, \ell_3') \delta_{\ell_1 \ell_2'}^{\text{K}} + B_{k'jk} (\ell_3', \ell_2, \ell_3) B_{i'j'i} (\ell_1', \ell_2', \ell_1) \delta_{\ell_1 \ell_3'}^{\text{K}} \right) \right. \\ \left. +3 \text{ terms obtained from perm. of } (i \leftrightarrow j, \ell_1 \leftrightarrow \ell_2) \right. \\ \left. +3 \text{ terms obtained from perm. of } (i \leftrightarrow k, \ell_1 \leftrightarrow \ell_3) \right] \\ + \frac{2\pi}{\Omega_{\text{sky}}} \left[ \frac{1}{\ell_1 \Delta \ell_1^b} \left( P_{ii'}^{\text{s.}}(\ell_1) T_{jkj'k'} (\ell_2, \ell_3, \ell_2', \ell_3') \delta_{\ell_1 \ell_1'}^{\text{K}} + P_{ij'}^{\text{s.}}(\ell_1) T_{jki'k'} (\ell_2, \ell_3, \ell_1', \ell_3') \delta_{\ell_1 \ell_2'}^{\text{K}} + P_{ik'}^{\text{s.}}(\ell_1) T_{jki'j'} (\ell_2, \ell_3, \ell_1', \ell_2') \delta_{\ell_1 \ell_3'}^{\text{K}} \right) \right. \\ \left. +3 \text{ terms obtained from perm. of } (i \leftrightarrow j, \ell_1 \leftrightarrow \ell_2) \right. \\ \left. +3 \text{ terms obtained from perm. of } (i \leftrightarrow k, \ell_1 \leftrightarrow \ell_3) \right], \qquad (B6) $$
while the covariances between the binned tomographic power spectrum and bispectrum are

$$
\begin{align}
\mathrm{Cov} \left[ \hat{P}_{ij}^{\mathrm{W}} (\ell^b), \hat{B}_{i'j'k'}^{\mathrm{W}} (\ell_1^b, \ell_2^b, \ell_3^b) \right]_{\mathrm{NGins,PB}} &= \nonumber \\
&= \frac{2\pi}{\Omega_{\mathrm{sky}}} \left[ \frac{1}{\ell_1 \Delta \ell_1^b} (P_{i'j}^{\mathrm{s}}(\ell) B_{ij'k'}(\ell, \ell_2, \ell_3) \delta_{\ell\ell_1}^{\mathrm{K}} + P_{ii'}^{\mathrm{s}}(\ell_1) B_{jj'k'}(\ell, \ell_2, \ell_3) \delta_{\ell\ell_1}^{\mathrm{K}}) + \right. \nonumber \\
&\qquad \left. + 2 \text{ terms obtained from perm. of } (i' \leftrightarrow j', \ell_1 \leftrightarrow \ell_2) + \right. \nonumber \\
&\qquad \left. + 2 \text{ terms obtained from perm. of } (i' \leftrightarrow k', \ell_1 \leftrightarrow \ell_3) \right]. \tag{B7}
\end{align}
$$

For the moment, we did not include the components related to the six- and five-order polyspectra. As anticipated, these terms are indeed affected by the window function of our survey (see Chan et al. (2018), Barreira (2019) for a 3D analysis) as

$$
\begin{equation}
\mathrm{Cov} \left[ \hat{B}_{ijk}^{\mathrm{W}} (\ell_1^b, \ell_2^b, \ell_3^b), \hat{B}_{i'j'k'}^{\mathrm{W}} (\ell_1'^b, \ell_2'^b, \ell_3'^b) \right]_{\mathrm{NG,P6}} = \nonumber \\
\qquad \sum_{\bar{\ell}_1, \bar{\ell}_2, \bar{\ell}_3} \sum_{\ell_1', \ell_2', \ell_3'} \int \frac{d^2 \ell''}{(2\pi)^2} |W(\ell'')|^2 P_{ijklj'k'} (\bar{\ell}_1, \bar{\ell}_2, \bar{\ell}_3 + \ell'', \bar{\ell}_1', \bar{\ell}_2', \bar{\ell}_3' - \ell'') \frac{\Delta_{\ell_1^b, \ell_2^b, \ell_3^b}^{(3)} (\bar{\ell}_1, \bar{\ell}_2, \bar{\ell}_3 + \ell'') \Delta_{\ell_1'^b, \ell_2'^b, \ell_3'^b}^{(3)} (\bar{\ell}_1', \bar{\ell}_2', \bar{\ell}_3' - \ell'')}{\Omega_{\mathrm{sky}}^2 N_{\mathrm{tri.}} (\ell_1^b, \ell_2^b, \ell_3^b) N_{\mathrm{tri.}} (\ell_1'^b, \ell_2'^b, \ell_3'^b)} \tag{B8}
\end{equation}
$$

$$
\mathrm{Cov} \left[ \hat{P}_{ij}^{\mathrm{W}} (\ell^b), \hat{B}_{i'j'k'}^{\mathrm{W}} (\ell_1^b, \ell_2^b, \ell_3^b) \right]_{\mathrm{NG,P5}} = \sum_{\ell'} \sum_{\ell_1' \ell_2' \ell_3'} \int \frac{d^2 \ell''}{(2\pi)^2} |W(\ell'')|^2 P_{ijklj'k'} (\bar{\ell}_1, \bar{\ell}_2, \bar{\ell}_3 + \ell'', \bar{\ell}', -\bar{\ell}' - \ell'') \frac{\Delta_{\ell^b}^{(2)}(\ell) \Delta_{\ell_1^b, \ell_2^b, \ell_3^b}^{(3)}(\ell_1, \ell_2, \ell_3)}{\Omega_{\mathrm{sky}}^2 N_{\mathrm{tri.}}(\ell^b) N_{\mathrm{tri.}}(\ell_1^b, \ell_2^b, \ell_3^b)}. \tag{B9}
$$

Similar to the power spectrum case, the Limber approximation allows us to determine the convergence polyspectra via a line-of-sight integration of the matter density ones, at every redshift. The consistency relations allow the following decomposition (Chan et al. 2018)

$$
P(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3 + \mathbf{p}, \mathbf{k}'_1, \mathbf{k}'_2, \mathbf{k}'_3 - \mathbf{p}) \approx P(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}'_1, \mathbf{k}'_2, \mathbf{k}'_3) + \frac{\partial B(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3 | \delta_b)}{\partial \delta_b} \frac{\partial B(\mathbf{k}'_1, \mathbf{k}'_2, \mathbf{k}'_3 | \delta_b)}{\partial \delta_b} P^{\text{lin}}(p), \quad (B10)
$$

$$
P(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3 + \mathbf{p}, \mathbf{k}', -\mathbf{k}' - \mathbf{p}) \approx P(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}', -\mathbf{k}') + \frac{\partial B(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3 | \delta_b)}{\partial \delta_b} \frac{\partial P(\mathbf{k}' | \delta_b)}{\partial \delta_b} P^{\text{lin}}(p). \quad (B11)
$$

Therefore, we can complement equations (B6) and (B7) with the standard intrasurvey covariance derived from the first terms

$$
\mathrm{Cov} \left[ \hat{B}_{ijk}^{\mathrm{W}} (\ell_1, \ell_2, \ell_3), \hat{B}_{i'j'k'}^{\mathrm{W}} (\ell_1', \ell_2', \ell_3') \right]_{\mathrm{NGins},6\mathrm{P}} \approx \frac{1}{\Omega_{\mathrm{sky}}} P_{ijklj'k'} (\ell_1, \ell_2, \ell_3, \ell_1', \ell_2', \ell_3'), \quad (\mathrm{B12})
$$

$$
\mathrm{Cov} \left[ P_{ij}(\ell), B_{i'j'k'}(\ell_1, \ell_2, \ell_3) \right]_{\mathrm{NGins},5\mathrm{P}} \approx \frac{1}{\Omega_{\mathrm{sky}}} P_{ijklj'k'}(\ell, -\ell, \ell_1, \ell_2, \ell_3). \quad (\mathrm{B13})
$$

The supersurvey ones instead come from the projected responses and we can nicely introduce a much more general expression for them

$$
\begin{align}
\mathrm{Cov} \left[ \hat{P}_{i_1 \dots i_n}^{\mathrm{W}} (\ell_1, \dots, \ell_n), \hat{P}_{i'_1 \dots i'_n}^{\mathrm{W}} (\ell'_1, \dots, \ell'_n) \right]_{\mathrm{NGscc}} &= \nonumber \\
&= \int_0^\infty \mathrm{d}\chi \ \mathcal{T}(i_1, \dots, i_n; \chi) \ \mathcal{T}(i'_1, \dots, i'_n; \chi) \frac{\partial P(\mathbf{k}(\ell_1, \chi), \dots, \mathbf{k}(\ell_n, \chi) | \delta_b)}{\partial \delta_b} \frac{\partial P(\mathbf{k}(\ell'_1, \chi), \dots, \mathbf{k}(\ell'_n, \chi) | \delta_b)}{\partial \delta_b} \sigma_W^2(\chi), \tag{B14}
\end{align}
$$

where the projection functions $\mathcal{T}$ are defined in equation (5). Also in equations (B12) and (B13), we assumed the six- and the five-point correlation functions to be well approximated by the 1-halo term. Therefore, we kept their dependence on just the magnitudes of the angular multipoles involved. The expression for the supersample term (B14) is consistent with the more general result in Lacasa & Rosenfeld (2016).

To estimate the variance $\sigma_W^2$ of the long mode at every redshift, we use a cylindrical mask in real space for the single observed patch, the Fourier expansion of which can be derived as done in Lima & Hu (2007)

$$
W_{\mathrm{cyl}}(\mathbf{k}, \delta_{\chi}, \hat{\chi}) = 2 \exp(i k_{\parallel} \hat{\chi}) j_0 \left( \frac{1}{2} k_{\parallel} \delta_{\chi} \right) \frac{J_1(k_{\perp} \hat{\chi} \Theta_{\mathrm{sky}})}{k_{\perp} \hat{\chi} \Theta_{\mathrm{sky}}}. \quad (\mathrm{B15})
$$

In particular equation (B15) is the Fourier transform of the selection function for a comoving cylindrical volume of depth $\delta_\chi$, centred in $\hat{\chi}$ and derived under the assumption of a slowly varying $\chi$ and Hubble factor $H(\chi)$ within $\delta_\chi$. Compared to the expression in Lima & Hu (2007) we are omitting the photometric error related component. The special functions $j_0$ and $J_1$ are respectively the 0th order spherical Bessel function and the 1st order Bessel function of the first kind. Since we are interested in computing the variance of the matter field at a specific redshift, we compute the variance $\sigma_W^2(\chi)$ over a disc-like volume in the limit $\delta_\chi \to 0$, $\hat{\chi} \to \chi$ of equation (B15)

$$
\sigma_W^2(\chi) = 4 \int \frac{d^2 \epsilon_{\perp}}{(2\pi)^2} \left[ \frac{J_1(\epsilon_{\perp} \chi \Theta_{\text{sky}})}{\epsilon_{\perp} \chi \Theta_{\text{sky}}} \right]^2 P^{\text{lin.}}(\epsilon_{\perp}, \chi). \quad (\text{B16})
$$

This paper has been typeset from a TeX/LA TeX file prepared by the author.