© The Authors 2025

1. Introduction

The current best practice for testing the cosmological model of the Universe using the large-scale structure (LSS) are two-point functions. These measurements typically involve the coherent weak lensing distortion of galaxy images and the correlation of galaxy positions. All of them have been successfully measured and used for cosmological inference by the Kilo-Degree Survey (KiDS, see e.g. Hildebrandt et al. 2017; Asgari et al. 2021b; Heymans et al. 2021)¹, the Dark Energy Survey (DES, see e.g. Abbott et al. 2018, 2022; Amon et al. 2022)², and the Subaru Hyper Suprime-Cam (HSC, see e.g. Sugiyama et al. 2022; Hamana et al. 2020; Dalal et al. 2023)³. These experiments have laid the groundwork for upcoming stage 4 surveys such as Euclid⁴, The Legacy Survey of Space and Time (LSST)⁵, and Roman⁶. Leveraging galaxy positions and shapes enables the measurement of three unique two-point functions, shape-shape (cosmic shear, CS), shape-position (galaxy-galaxy-lensing, GGL), and position-position (galaxy clustering, GC). Aggregating these into a data vector establishes the standard 3×2-point analysis. The main ingredient in this kind of analysis is the likelihood, which is commonly assumed to be Gaussian, so that it can be completely characterised by the data vector, model prediction, and covariance.

There are various approaches to computing or estimating the covariance matrix. The most direct method involves utilising the data itself, generating pseudo-realisations through subsampling (Norberg et al. 2009; Friedrich et al. 2016; Mohammad & Percival 2022; Trusov et al. 2024). Alternatively, simulation suites can be used to estimate the covariance if enough realisations are available and feasible (Hartlap et al. 2007; Harnois-Déraps et al. 2012, 2018; Dodelson & Schneider 2013; Heymans et al. 2013; Taylor et al. 2013; Taylor & Joachimi 2014; Joachimi 2017; Sellentin & Heavens 2018). Lastly, theoretical modelling is usually the most efficient and arguably the most fundamental way to estimate the covariance matrix (Takada & Jain 2004; Eifler et al. 2009; Krause & Eifler 2017; Reischke et al. 2017; Joachimi et al. 2021; Friedrich et al. 2021). Each approach has its own advantages and disadvantages. Simulations tend to be relatively accurate, with systematic effects such as masking or variable survey depth easily accounted for. However, obtaining an unbiased covariance estimate necessitates numerous realisations, with at least as many simulations required as there are effective degrees of freedom in the data for the covariance matrix to be invertible. Current surveys typically involve around 10²−10³ data points, while upcoming stage 4 surveys will boast 10⁴ data points for standard two-point statistics. Consequently, significant data compression (e.g. Heavens et al. 2000, 2017; Asgari & Schneider 2015) or running numerous computationally intensive simulations is necessary. Conversely, theoretical covariances are comparatively inexpensive to compute, although modelling the aforementioned systematic effects is difficult and a simulation-based covariance matrix estimation might be more accurate. While subsampling techniques offer a potential quick way to self-consistently estimate the covariance directly from the data, they tend to be biased on large scales where the subsamples are not independent due to cosmic variance.

Due to the symmetries of the Friedmann-Lemaître-Robertson-Walker metric, theoretical computations are typically conducted in harmonic (Fourier) space. Rotational symmetry results in the decoupling of different angular multipoles, ℓ, and allow considerable data compression. However, observations of a masked sky compromise the orthogonality of the spherical harmonic basis, leading to mode coupling that necessitates careful modelling, known as pseudo-C_ℓ (see e.g. Hikage et al. 2011, 2019; Becker et al. 2016; Alonso et al. 2019; Nicola et al. 2021; Loureiro et al. 2022; von Wietersheim-Kramsta et al. 2025). Unlike in harmonic space, standard real space correlation functions do not face these limitations to the first order in signal modelling, they encompass w(θ), γ_t(θ), and ξ_±(θ) for GC, GGL, and CS respectively⁷ However, photometric galaxy surveys introduce a broad radial window function, mixing different physical scales and complicating efforts to disentangle them and mitigate astrophysical effects on small scales. Moreover, correlation functions integrate over the angular power spectra with an oscillatory but broad integration kernel, resulting in limited control over the physical scales influencing the measurements. Leakage from very low multipoles can alter the Gaussianity of the likelihood at large angular scales (Sellentin et al. 2018; Louca & Sellentin 2020; Oehl & Tröster 2024), while significant contributions from highly non-linear and theoretically uncertain scales occur at small angular separations. Importantly, the cosmological lensing signal does not produce B-modes to first order, as it represents the gradient of the lensing potential. Consequently, the B-mode signal serves as a consistency check for systematic effects. However, real space correlation functions mix E- and B-modes, complicating the interpretation of these tests (see e.g. Schneider & Kilbinger 2007).

The aforementioned problems of real space correlation functions are partially alleviated by considering different derived summary statistics. In weak lensing surveys and in particular KiDS, two alternatives to real space correlation functions used are band powers (Schneider et al. 2002) and COSEBIs (Schneider et al. 2010). While the former approximately separates E and B modes, the latter yields a complete separation by definition. The demonstration of consistent parameter constraints over a set of summary statistics adds significant confirmation to the robustness of the pipelines used, the systematics involved and the overall fidelity of the cosmological inference. Although many of these methods have been integrated within the community for a 3×2-point analysis and are publicly available for model predictions, a unified framework for the covariance matrix in the cosmological community for all kinds of summary statistics has been absent.

While some of the previously mentioned tools for calculating the covariance matrix are publicly available (see e.g. Krause & Eifler 2017; Fang et al. 2020), they possess limited flexibility or lack comprehensive functionality when it comes to exchanging tracer fields, summary statistics, weight functions, bias, or halo occupation distribution prescriptions. Furthermore, it is not always straightforward to input files of pre-calculated quantities into these codes. Therefore, this paper presents the analytical covariance matrices alongside a Python code, the OneCovariance code, encompassing all functionalities used in the standard analysis of KiDS, ranging from CS (Asgari et al. 2021b) alone to a full 3×2-point analysis (Joachimi et al. 2021; Heymans et al. 2021), incorporating the stellar mass function (Dvornik et al. 2023), and employing all three types of summary statistics (i.e. real space correlation functions, bandpowers, and COSEBIs/Ψ-stats⁸) used in KiDS for CS, galaxy clustering, and galaxy–galaxy lensing as well as arbitrary cross-correlations between them. Moreover, it remains flexible enough to accommodate very general inputs and therefore can be used by any survey or collaboration to obtain quick and easy covariance matrices for photometric LSS observables.

The goal of this paper is threefold: (i) to furnish a comprehensive and instructive document encompassing all components required for theoretical computations of LSS covariances in photometric galaxy surveys, consolidating everything into a unified code, which we call OneCovariance, accessible to the cosmological community for generating covariance matrices effortlessly for any projected LSS parameter; (ii) to introduce the mocks used for the validation of the covariance matrix in the final CS analysis of the Kilo-Degree Survey (KiDS-Legacy) featuring, among other things, a realistic footprint and variable number density; and (iii) to employ the OneCovariance code in different scenarios. Our focus here, in addition to CS in KiDS-Legacy, includes the covariance matrix of clustering redshifts (van den Busch et al. 2020) and the cross-correlation between different summary statistics (Asgari et al. 2021b). It should be highlighted that the OneCovariance code, while showcased for KiDS, is completely general and can be used for all different types of analyses via simple file inputs to acquire fast covariance estimates.

Since the CS catalogue of KiDS-Legacy is not final at this stage, we only performed simulations and comparisons for KiDS-Legacy-like settings. However, we do not anticipate any major changes that would influence the covariance validation. We refer to data products that are final as of the writing of this paper as KiDS-Legacy, and to simulations or data products with nearly final settings as KiDS-Legacy-like. Furthermore, we generally refer to KiDS-Legacy as the entirety of the final analysis of the last data release of KiDS.

The manuscript is organised as follows. In Section 2 we discuss some general properties of covariance matrices in cosmology, its definition, the different contributions, and the flat sky approximation. In Section 3 we discuss the prescription for non-linear structure formation. Here we focus on a halo model-based description of structure formation, and we discuss all its ingredients. The projection along the line of sight and the corresponding covariance in harmonic space is discussed in Section 4. We then discuss the projection to observables in Section 5 and describe the motivation and features of the OneCovariance code in Section 6. The subsequent sections feature applications of the code to specific examples. We introduce a KiDS-Legacy-like CS sample in Section 7 and describe the general characteristics of the covariance matrix as well as the influence of its modelling choices on parameter inference. Then, in Section 8 we present the KiDS-Legacy-like mocks, incorporating variable depth, and we demonstrate the overall agreement between the analytical and numerical covariances. The results section concludes with Section 9, which showcases two additional applications of the covariance matrix utilised in KiDS-Legacy: covariances between different summary statistics for consistency tests and the covariance for clustering redshift calibration. To facilitate swift access, we highlight key plots at the outset of each results section. Finally, in Section 10 and Section 11 we examine the current status and capabilities of the code, along with future directions and our conclusions.

2. Covariance matrices in photometric galaxy surveys

In this section, we outline the general theoretical modelling involved in constructing the covariance matrix. We separate this discussion from the precise prescription of non-linear structure formation, which is described in more detail in Section 3. This separation allows us to distinguish observational definitions from the theoretical framework as clearly as possible.

2.1. General considerations

Cosmological analyses of the LSS hinge on estimating a two-point function between any pair of observables (tracers), denoted as t_i[δ], which are functionals of the underlying matter field, δ(x,χ). Here, x represents a spatial co-moving vector (i.e. a 3-dimensional vector in a spatial hyper-surface), and χ denotes the co-moving distance. These two-point functions can be aggregated into a vector $\text{O}$ with components

$${\mathcal{O}}_{\alpha }≔\langle t_i[\delta ]t_j[\delta ]\rangle \hspace{0.5em}\alpha \in \{\text{all unique combinations of}(t_i,t_j)\}.$$(1)

Here, the angular brackets, 〈·〉, denote an ensemble average, akin to a spatial average (for more details on this ergodicity argument see e.g. Desjacques et al. 2021). The number of unique combinations depends on the tracer used; if the underlying tracer is the same and i simply labels a subset such as tomographic bins, i≤j, otherwise, all combinations can exist. The components of the covariance matrix are then defined as

$$(\text{C}{)}_{\alpha \beta }\equiv \mathrm{Cov}({\text{O}}_{\alpha },{\text{O}}_{\beta })≔\langle ({\text{O}}_{\alpha }-\langle {\text{O}}_{\alpha }\rangle )({\text{O}}_{\beta }-\langle {\text{O}}_{\beta }\rangle )\rangle .$$(2)

Typically, we are provided with the statistical properties of the density field δ from theoretical calculations or numerical simulations, rather than those of its tracer. If the mapping from δ to t_i is linear, all the statistical properties of δ extend accordingly to t_i. However, this is not always the case, such as in a perturbative bias expansion or the reduced shear beyond linear order. For simplicity, we assume 〈t_i[δ]〉=0, resulting in the following breakdown of Equation (2) using Equation (1) and Wick's theorem:

$$C_{\mathit{ijmn}}=\langle t_i{t}_m\rangle \langle t_j{t}_n\rangle +\langle t_i{t}_n\rangle \langle t_j{t}_m\rangle +\langle t_i{t}_j{t}_m{t}_n{\rangle }_{\mathrm{c}}.$$(3)

Here, we have omitted the argument of the different tracers, [δ], for notational simplicity. The first two terms are known as Gaussian terms, deriving from products of the original observable, while the subscript ‘c’ in the last term denotes the connected part of the correlator, i.e. the term originating from a genuine four-point function and not the product of two-point function. This quantity vanishes for Gaussian fields. For surveys with incomplete sky coverage, the connected part of the covariance can be decomposed into a part arising from modes within the survey, termed the non-Gaussian (nG) component, and a term generated by modes larger than the survey, known as the super-sample covariance (SSC, Takada & Hu 2013):

$$\langle t_i{t}_j{t}_m{t}_n{\rangle }_{\mathrm{c}}=\langle t_i{t}_j{t}_m{t}_n{\rangle }_{\mathrm{nG}}+\langle t_i{t}_j{t}_m{t}_n{\rangle }_{\mathrm{SSC}}.$$(4)

It should be noted that both these terms are entirely arising from the connected term. Particularly, the SSC term stems from the cosmic variance's reaction to unobservable modes, manifested as the squeezed limit of a connected four-point function (Linke et al. 2024).

Real observables invariably serve as discrete tracers of the underlying smooth field, hence they are subject to stochastic noise due to the finite number of tracers, such as galaxies. Hence, one can schematically express

$${\hat{t}}_i=t_i+n_i,$$(5)

where n_i denotes some stochastic noise component, this component emerges only when the same object correlates with itself and is thus also referred to as the shot or shape noise component for galaxy clustering or CS, respectively. With this understanding, assuming that the signal does not correlate with the noise, the Gaussian terms, i.e., the first two terms in Equation (3), each decompose into three components,

$$\langle {\hat{t}}_i{\hat{t}}_m\rangle \langle {\hat{t}}_j{\hat{t}}_n\rangle \sim \underset{\text{sample variance}}{\underbrace{\langle t_i{t}_m\rangle \langle t_j{t}_n\rangle }}+\underset{\text{mixed}}{\underbrace{\langle t_i{t}_m\rangle \langle n_j{n}_n\rangle +\langle t_j{t}_n\rangle \langle n_i{n}_m\rangle }}+\underset{\mathrm{shot}/\mathrm{shape}\hspace{0.5em}\mathrm{noise}}{\underbrace{\langle n_j{n}_n\rangle \langle n_i{n}_m\rangle }},$$(6)

that are distinguished by their scaling with the signal. We note that the sample variance term is often called cosmic variance in the literature.

2.2. Angular polyspectra and Limber approximation

Observables of (photometric) galaxy surveys are projected quantities of three-dimensional fields. Consider such a three-dimensional field A(x,f_Kχ(z))∈{t_i} at position x and redshift z. Its projected version, $a(\hat{n})$, is

$$a(\hat{n})={\displaystyle \int }\mathrm{d}\chi \hspace{0.5em}{W}_A(\chi )A(f_{\mathrm{K}}(\chi )\hat{n},\chi ),$$(7)

where W_A(χ) is some weighting function whose exact form is probe specific and will be discussed in Section 4. $\hat{n}$ is the unit vector on the sphere and f_K(χ) is the co-moving angular diameter distance, such that $x=f_{\mathrm{K}}(\chi )\hat{n}$ and χ the co-moving radial distance. In Appendix B the definitions for a spherical harmonic decomposition of the projected field, a, are given. Leveraging those relations, the two-point function defines the angular power spectrum, ${\mathcal{P}}_{a_1{a}_2}({\ell }_1)$, of the two-dimensional fields a₁ and a₂ as follows:

$$\begin{array}{cc}\hfill \langle a_{1,{\ell }_1{m}_1}^{}{a}_{2,{\ell }_2{m}_2}^{}\rangle ≔& \hspace{0.5em}{\mathcal{P}}_{a_1{a}_2}({\ell }_1){\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}{\delta }_{m_1{m}_2}^{\mathrm{K}}\hfill \\ \hfill =& \hspace{0.5em}\frac{2}{\pi }{\displaystyle \int }\mathrm{d}{\chi }_1{W}_{a_1}({\chi }_1){\displaystyle \int }\mathrm{d}{\chi }_2{W}_{a_2}({\chi }_2){\displaystyle \int }{k}^2\mathrm{d}k\hspace{0.5em}{P}_{A_1{A}_2}(k,{\chi }_1,{\chi }_2)j_{\ell }({\mathit{kf}}_{\mathrm{K}}({\chi }_1))j_{\ell }({\mathit{kf}}_{\mathrm{K}}({\chi }_2))\hspace{0.5em}{\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}{\delta }_{m_1{m}_2}^{\mathrm{K}}.\hfill \end{array}$$(8)

Here the dependence on χ₁ and χ₂ of the three-dimensional power spectrum, $P_{A_1{A}_2}(k,{\chi }_1,{\chi }_2)$, was made explicit. ${\delta }_{\mathit{ij}}^{\mathrm{K}}$ is the Kronecker delta. The more general equation for higher-order correlators is provided in Appendix C. Throughout this paper the geometric mean is used to approximate the unequal time correlator such that for the power spectrum

$$P_{A_1{A}_2}(k,{\chi }_1,{\chi }_2)\approx {\left[P_{A_1{A}_2}(k,{\chi }_1)P_{A_1{A}_2}(k,{\chi }_2)\right]}^{1/2},$$(9)

which is an excellent approximation for almost all practical purposes (Castro et al. 2005; Kitching & Heavens 2017; Kilbinger et al. 2017; Chisari & Pontzen 2019; de la Bella et al. 2021).

In the flat sky limit (in which case the transform to spherical harmonics essentially becomes a two-dimensional Fourier transform, see Appendix D) angular polyspectra can be approximated using the Limber approximation (Limber 1953; Loverde & Afshordi 2008; Leonard et al. 2023):

$$\langle a_{{\ell }_1{m}_1}\dots a_{{\ell }_n{m}_n}\rangle \sim {\displaystyle \int }\frac{\mathrm{d}\chi }{f_{\mathrm{K}}(\chi {)}^{2n-2}}{W}_{A_1}(\chi )\dots W_{A_n}(\chi )P_{A_1\dots A_n}\left(\frac{{\ell }_1+0.5}{f_{\mathrm{K}}(\chi )},\dots ,\frac{{\ell }_n+0.5}{f_{\mathrm{K}}(\chi )},\chi \right).$$(10)

Here geometrical factors occurring in Equation (C.7) are disregarded due to the flat sky approximation. For n>3 this expression assumes that the angular average over possible configurations has already been carried out on the flat sky. The use of the Limber approximation introduces negligible biases in cosmological parameter inference when considering CS only (ss e.g. Joachimi et al. 2021; Leonard et al. 2023). Therefore, unless stated otherwise, we utilise Equation (10) for projected fields.

3. Non-linear evolution: Halo model

Due to the non-linearity of gravity and the complex interaction between cold dark matter and baryons, understanding structure formation in cosmology presents a formidable challenge. Various approaches have been developed to tackle this issue. These include standard perturbation theory (SPT, for an exhaustive review see Bernardeau et al. 2002), halo model approaches (HM, see e.g. Seljak 2000; Cooray & Sheth 2002; Asgari et al. 2023), effective field theory (EFT, Carrasco et al. 2012), or kinetic field theory (KFT, e.g. Bartelmann et al. 2021).

In addition to modelling the non-linear evolution of the smooth dark matter density field, there are additional complications: (i) The total matter power spectrum encapsulates the clustering of both cold dark matter and baryons. Processes such as star formation, supernovae, active galactic nuclei (AGN), and cooling can displace baryons relative to dark matter, leading to a spatial distribution of total matter that differs from that of cold dark matter alone. This phenomenon is collectively known as baryonic feedback (see e.g. Semboloni et al. 2013; Somerville & Davé 2015; Vogelsberger et al. 2020; Chisari et al. 2019b; Nicola et al. 2022). (ii) The tracers discussed in Section 2 are biased representations of the underlying matter distribution. This bias is most prominent in galaxy surveys (referred to as galaxy bias, see e.g. Desjacques et al. 2018, for a review), but also affects CS measurements through intrinsic alignment (IA) effects (see e.g. Schaefer 2009; Joachimi et al. 2015; Vlah et al. 2020; Fortuna et al. 2020, for reviews), which can be treated similarly to galaxy bias using an effective field theory (EFT) approach.

While perturbative approaches such as EFT are technically more rigorous than non-perturbative (but phenomenological) halo model calculations, they do not allow progressing into the highly non-linear regime, which is crucial for signals with low signal-to-noise ratios such as CS. In other words, the significance of the CS measurement is not large enough to fit all the EFT nuisance parameters and all cosmological information would be lost. Hence, incorporating prior knowledge regarding the statistical properties of the matter distribution via a halo model is key if one wants to access highly non-linear scales. Although KFT offers, in theory, an alternative route, it still has to be shown to provide good fits for models including baryonic feedback and higher-order statistics necessary for the covariance. Since the halo model has been shown to be flexible enough to describe different cosmological observations into the non-linear regime⁹ (see Mead et al. 2020, and references therein), we stick to a halo model approach in this paper, following the methodology of KiDS-1000 (Joachimi et al. 2021). However, it is worth noting that the code presented here is sufficiently flexible to accommodate any external power spectra, power spectrum response, and/or tri-spectra for subsequent calculations. The halo model assumes that the cosmic density field is entirely composed of dark matter halos, whose distribution is governed by a mass function. By assuming a density profile and a halo bias, we can calculate the statistical properties of the matter field. Populating each halo with galaxies using a halo occupation distribution (HOD) enables a versatile model that accurately predicts the statistical properties of CS, GC, GGL for the cross-correlation between the two galaxy samples.

3.1. Halo model ingredients: Galaxy bias and halo occupation distribution

Here, we briefly summarise the components involved in the halo model calculations of matter and galaxy polyspectra. Galaxy bias is addressed through an HOD prescription (van den Bosch et al. 2013; Dvornik et al. 2018, 2023; Lacasa et al. 2014; Lacasa 2018; Reischke et al. 2020). To achieve this, one employs a biasing function, b(M), depending on halo mass M, defined via the average number density of galaxies 〈N|M〉 in a halo of mass M

$$b(M)≔\frac{{\overline{\rho }}_{\mathrm{m}}(z)}{{\overline{n}}_{\mathrm{g}}(z)}\frac{\langle N|M\rangle }{M},$$(11)

where ${\overline{n}}_{\mathrm{g}}$ is the average number density of galaxies and ${\overline{\rho }}_{\mathrm{m}}(z)$ is the mean matter density in the Universe at redshift z. Average quantities in the occupation distribution are defined as

$$\langle N|M\rangle ≔{\displaystyle \sum _{N=0}^{\infty }}\mathit{NP}(N|M),$$(12)

where P(N|M) is the probability of a halo to be occupied with N galaxies. Moments of the bias at a redshift z, or any halo-related quantity, A(M,z), are calculated via the differential number density of halos of mass M: the halo mass function n_h(M,z)

$$\langle A\rangle (z)≔{\displaystyle \int }\mathrm{d}M\hspace{0.5em}{n}_{\mathrm{h}}(M,z)A(M,z).$$(13)

The halo model assumes that all matter in the Universe is bound in halos. We assume those halos to be spherically symmetric and to have an average density profile ρ_h(r|M)=Mu_h(r|M), such that u_h is normalised and can be used in Equation (13) without introducing additional normalisation factors. The galaxy populations can be split into central, ‘c’, and satellite, ‘s’, galaxies with different halo profiles u_c/s and halo occupations 〈N_c/s|M〉. If there is a central galaxy in a halo on average, 〈N_c|M〉=1 the distribution of additional galaxies, i.e. the satellites, is assumed to be Poisonnian. This allows us to calculate the expected number of n-tuples, 〈N⁽ⁿ⁾|M〉 in a halo, as required for the 1, …, 4-halo terms

$$\begin{array}{cc}\hfill \langle N|M\rangle & =\langle N_{\mathrm{c}}|M\rangle +\langle N_{\mathrm{s}}|M\rangle \hspace{0.5em}\hspace{0.5em}\hspace{0.5em}n=1,\hfill \\ \hfill \langle N(N-1)|M\rangle & =\langle N_{\mathrm{c}}|M\rangle \left(\langle N_{\mathrm{s}}|M{\rangle }^2+2\langle N_{\mathrm{s}}|M\rangle \right)n=2,\hfill \\ \hfill \langle N(N-1)(N-2)|M\rangle & =\langle N_{\mathrm{c}}|M\rangle \left(\langle N_{\mathrm{s}}|M{\rangle }^3+3\langle N_{\mathrm{s}}|M\rangle \right)n=3,\hfill \\ \hfill & \dots ,\hfill \end{array}$$(14)

where the number of tuples for each individual type e.g. ‘s’ and ‘c’ can be read off. In the case of matter, ‘m’, the HOD is set to the halo mass M for each point and the normalisation is taken care of by ${\overline{\rho }}_{\mathrm{m}}$, i.e. 〈N_m|M〉=M. It should be noted that sub-Poissonian behaviour was found in Dvornik et al. (2023) for the satellite population. This could, in principle, be accounted for by including an additional free parameter on the two-point level. For the connected non-Gaussian term, this might lead to more complications. However, we do not expect that the latter is dominating and would change cosmological results significantly. It is instructive to introduce the following integrals (Cooray & Hu 2001; Takada & Hu 2013; Lacasa 2018) for the species X₁,…,X_μ:

$$I_{\mu ,X_{\{\mu \}}}^{\alpha }(k_1,\dots ,k_{\mu },z)≔\frac{1}{{\overline{n}}_{X_1}(z)\dots {\overline{n}}_{X_{\mu }}(z)}\int \mathrm{d}M\hspace{0.5em}{n}_{\mathrm{h}}(M,z)\langle N_{\{X\}}^{(\mu )}|M\rangle b_{\mathrm{h}}^{\alpha }(M,z){\tilde{u}}_{X_1}(k_1|M,z)\dots {\tilde{u}}_{X_{\mu }}(k_{\mu }|M,z).$$(15)

Here $b_{\mathrm{h}}^{\alpha }$ is the halo bias, for which we use the fitting function from Tinker & Wetzel (2010). $\tilde{u}(k|M)$ is the Fourier transform of the normalised halo profile. In Appendix E we provide $\tilde{u}(k|M)$ for an NFW profile. However, they are in general arbitrary. Lastly, the average number density n_X(z) can be calculated using Equation (13):

$${\overline{n}}_X(z)=\int \mathrm{d}M\hspace{0.5em}{n}_{\mathrm{h}}(M,z)\langle N_X|M\rangle .$$(16)

We note again that ${\overline{n}}_{\mathrm{m}}(z)={\overline{\rho }}_{\mathrm{m}}(z)$ with this definition.

3.2. Conditional stellar mass function

It remains to specify the occupation statistics of a galaxy with halo mass M. As the fiducial model, we chose the conditional stellar mass function (CSMF, Yang et al. 2009; Cacciato et al. 2009, 2013; van Uitert et al. 2016; Dvornik et al. 2018, 2023), specifying the average distribution of galaxies of stellar mass M_★ given they occupy a halo with mass M. It is also split into centrals and satellites so that the total CSMF is

$$\mathrm{\Phi }(M_{★}|M)={\mathrm{\Phi }}_{\mathrm{c}}(M_{★}|M)+{\mathrm{\Phi }}_{\mathrm{s}}(M_{★}|M).$$(17)

The two contributions are modelled as a log-normal and a modified Schechter function

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\mathrm{c}}(M_{★}|M)& =\frac{1}{\sqrt{2\pi }ln(10){\sigma }_{\mathrm{c}}{M}_{★}}exp\left[-\frac{log(M_{★}/M_{\mathrm{c}}^{*}{)}^2}{2{\sigma }_{\mathrm{c}}^2}\right],\hfill \\ \hfill {\mathrm{\Phi }}_{\mathrm{s}}(M_{★}|M)& =\frac{{\varphi }_{\mathrm{s}}^{*}}{M_{\mathrm{s}}^{*}}{\left(\frac{M_{★}}{M_{\mathrm{s}}^{*}}\right)}^{{\alpha }_{\mathrm{s}}}exp\left[-{\left(\frac{M_{★}}{M_{\mathrm{s}}^{*}}\right)}^2\right],\hfill \end{array}$$(18)

for centrals and satellites respectively¹⁰. Here σ_c describes the scatter of stellar mass given a halo mass while α_s describes the power law scaling of the satellite abundance. It should be noticed that if α_s<0 the expression for Φ_s(M_★|M) would formally diverge for low masses. However, the existence of satellites in the HOD are always conditioned on the existence of a central galaxy, thus ensuring convergence at low masses due to the exponential cut-off of Φ_c(M_★|M). All free parameters in this model can, in principle, be arbitrary functions of the halo mass, M. For this fiducial case, we follow Yang et al. (2009), Dvornik et al. (2018, 2023) and use the following parametrisation:

$$\begin{array}{cc}\hfill M_{\mathrm{c}}^{*}& =\hspace{0.5em}{M}_0\frac{(M/M_1{)}^{{\gamma }_1}}{(1+[M/M_1]{)}^{{\gamma }_1-{\gamma }_2}},\hfill \\ \hfill M_{\mathrm{s}}^{*}& =\hspace{0.5em}0.56M_{\mathrm{c}}^{*}(M),\hfill \\ \hfill log{\varphi }_{\mathrm{s}}^{*}& =\hspace{0.5em}{b}_0+b_{1}log\left(\frac{M}{{10}^{13}{M}_{\odot }{h}^{-1}}\right).\hfill \end{array}$$(19)

Integrating the CSMF over the HMF yields the galaxy SMF:

$${\mathrm{\Phi }}_X(M_{★})=\int \mathrm{d}M\hspace{0.5em}{\mathrm{\Phi }}_X(M_{★}|M)n_{\mathrm{h}}(M,z).$$(20)

Likewise, the occupation number in a stellar mass bin [M_★,1,M_★,2] is given by

$$\langle N_X|M\rangle ={\displaystyle {\int }_{M_{★,1}}^{M_{★,2}}}\mathrm{d}{M}_{★}{\mathrm{\Phi }}_X(M_{★}|M),$$(21)

completing the HOD prescription. It should be noted that the covariance code presented here is modular and can support other HOD prescriptions.

3.3. Halo model polyspectra

With the ingredients defined in the previous section, one can calculate the power spectrum between species A₁ and A₂, ignoring non-linear bias terms (Mead & Verde 2021), as

$$P_{A_1{A}_2}(k,z)=P_{A_1{A}_2}^{1\mathrm{h}}(k,z)+P_{A_1{A}_2}^{2\mathrm{h}}(k,z)\equiv I_{2,A_1,A_2}^0(k,k,z)+I_{1,A_1}^1(k,z)I_{1,A_2}^1(k,z)P_{\mathrm{lin}}(k,z),$$(22)

with the 1- and 2-halo term as well as the linear matter power spectrum P_lin. The one-halo term will tend to be a constant, as k→0 due to the infinite support of the halo profile. Due to this, it will overcome the two-halo term on large scales. To remove this unphysical behaviour, we introduce a large-scale damping for the 1h term in Equation (22)

$$P_{A_1{A}_2}^{1\mathrm{h}}(k,z)\hspace{0.5em}\to \hspace{0.5em}{P}_{A_1{A}_2}^{1\mathrm{h}}(k,z)\hspace{0.5em}\mathrm{erf}(k/k_{\mathrm{damp}}),$$(23)

with a fiducial damping scale of k_damp=0.1 h/Mpc. Likewise, the halo model trispectrum is split into 1-4-halo terms

$$\begin{array}{cc}\hfill T_{X_1\dots A_4}(k_1,k_2,k_3,k_4,z)=& \hspace{0.5em}{T}_{A_1\dots A_4}^{1\mathrm{h}}+T_{A_1\dots A_4}^{2\mathrm{h}}+T_{A_1\dots A_4}^{3\mathrm{h}}+T_{A_1\dots A_4}^{4\mathrm{h}}\hspace{0.5em},\hfill \end{array}$$(24)

where it is understood that the right-hand side depends as well on (k₁,k₂,k₃,k₄,z). The individual terms are given by

$$T_{A_1\dots A_4}^{1\mathrm{h}}=\hspace{0.5em}{I}_{4,A_1\dots A_4}^0(k_1,k_2,k_3,k_4),$$(25)

$$T_{A_1\dots A_4}^{2\mathrm{h}}=P_{\mathrm{lin}}(k_1)I_{3,A_2{A}_3{A}_4}^1(k_2,k_3,k_4)I_{1,A_1}^1(k_1)+3\hspace{0.5em}\mathrm{perm}.+P_{\mathrm{lin}}(k_{12})I_{2,A_1{A}_2}^1(k_1,k_2)I_{2,A_3{A}_4}^1(k_3,k_4)+2\hspace{0.5em}\mathrm{perm}.,$$(26)

$$T_{A_1\dots A_4}^{3\mathrm{h}}=\hspace{0.5em}{B}_{\mathrm{tree}}(k_1,k_2,k_{34})I_{2,A_3{A}_4}^1(k_3,k_4)I_{1,A_1}^1(k_1)I_{1,A_2}^1(k_2)+P_{\mathrm{lin}}(k_1)P_{\mathrm{lin}}(k_2)$$(27)

$$\begin{array}{cc}\hfill & +I_{2,A_3{A}_4}^1(k_3,k_4)I_{1,A_1}^1(k_1)I_{1,A_2}^1(k_2)+5\hspace{0.5em}\mathrm{perm}.,\hfill \\ \hfill T_{A_1\dots A_4}^{4\mathrm{h}}=& \hspace{0.5em}\prod _{\mu }{I}_{1,A_{\mu }}^1(k_{\mu })T_{\mathrm{tree}}(k_1,k_2,k_3,k_4).\hfill \end{array}$$(28)

We neglect the explicit dependence on redshift here and refer to Lacasa (2018) for an in-depth discussion. The perturbation theory bispectrum and trispectrum is given by (Fry 1984)

$$B_{\mathrm{tree}}(k_1,k_2,k_3)=\hspace{0.5em}2F_2^{\mathrm{s}}(k_1,k_2)P_{\mathrm{lin}}(k_1)P_{\mathrm{lin}}(k_2)+2\hspace{0.5em}\mathrm{perm}.,$$(29)

$$T_{\mathrm{tree}}(k_1,k_2,k_3,k_4)=\hspace{0.5em}4\left[F_2^{\mathrm{s}}(k_{12},-k_1)F_2^{\mathrm{s}}(k_{12},k_3)P_{\mathrm{lin}}(k_1)P_{\mathrm{lin}}(k_{12})P_{\mathrm{lin}}(k_3)+12\hspace{0.5em}\mathrm{perm}.\right]$$(30)

$$+\hspace{0.5em}6\left[P_{\mathrm{lin}}(k_1)P_{\mathrm{lin}}(k_2)P_{\mathrm{lin}}(k_3)F_3^{\mathrm{s}}(k_1,k_2,k_3)+4\hspace{0.5em}\mathrm{perm}.\right],$$

where we omitted the explicit dependence on redshift. For the symmetric perturbation theory kernels, $F_i^{\mathrm{s}}$, we refer to Bernardeau et al. (2002). The trispectrum relates to the nG term in Equation (4) for k₁=−k₂ and k₂=−k₄. Lastly, we specify the SSC (Takada & Hu 2013; Krause & Eifler 2017; Barreira et al. 2018), i.e. the second term in Equation (4):

$$T_{A_1\dots A_4}^{\mathrm{SSC}}(k_1,k_2,z)=\frac{\partial P_{A_1{A}_2}}{\partial {\delta }_{\mathrm{bg}}}(k_1,z)\frac{\partial P_{A_3{A}_4}}{\partial {\delta }_{\mathrm{bg}}}(k_2,z){\sigma }_{\mathrm{bg},(A_1{A}_2)(A_3{A}_4)}^2(z).$$(31)

We note that k₁ corresponds to the correlated pair A₁,A₂ and k₂ to A₃,A₄. The responses are given as

$$\begin{array}{cc}\hfill \frac{\partial P_{A_1{A}_2}}{\partial {\delta }_{\mathrm{bg}}}(k,z)=& \hspace{0.5em}\left[\frac{68}{21}-\frac{1}{3}\frac{\mathrm{d}log{k}^3{P}_{A_1{A}_2}(k,z)}{\mathrm{d}logk}\right]I_{1,A_1}^1(k,z)I_{1,A_2}^1(k,z)P_{\mathrm{lin}}(k,z)\hfill \\ \hfill & +I_{2,A_1{A}_2}^1(k,k,z)-{\left([b_{A_1}+b_{A_2}]P_{A_1{A}_2}(k,z)\right)}_{\mathrm{if}\hspace{0.5em}{A}_1,A_2\ne \mathrm{m}}.\hfill \end{array}$$(32)

We note that in the last term the contribution of $b_{A_i}$ vanishes if the mean of A_i is not constructed in the survey window, such as for CS. The responses are defined with respect to fluctuations in the background, ‘bg’, matter field, which can be calculated directly from the survey mask and its spherical harmonic decomposition:

$${\sigma }_{\mathrm{bg},(A_1{A}_2)(A_3{A}_4)}^2(z)=\frac{1}{A_{(A_1{A}_2)}{A}_{(A_3{A}_4)}}{\displaystyle \sum _{\ell }}{P}_{\mathrm{lin}}(\ell /\chi ){\displaystyle \sum _m}{a}_{\ell m}^{*(A_1{A}_2)}{a}_{\ell m}^{(A_3{A}_4)}.$$(33)

The bracket notation (A₁A₂) denotes the footprint of the survey over which the summary statistic between the tracers A₁ and A₂ has been evaluated. As expected ${\sigma }_{\mathrm{bg}}^2\to 0$ as $A_{(A_1{A}_2)}\to \infty$ due to statistical homogeneity and isotropy. If the survey footprint is not available via a mask file, we assume a circular footprint and follow Li et al. (2014)

$${\sigma }_{\mathrm{bg},(A_1{A}_2)(A_3{A}_4)}^2(z)=\chi (z{)}^2\int \frac{k\mathrm{d}k}{(2\pi {)}^2}{P}_{\mathrm{lin}}(k,z){\left[\frac{2{\mathrm{J}}_1(k\chi (z){\theta }_{\mathrm{s}})}{k\chi (z){\theta }_{\mathrm{s}}}\right]}^2,$$(34)

with a cylindrical Bessel function, J₁, and the survey size ${\theta }_{\mathrm{s}}=\sqrt{\max (A_{(A_1{A}_2)}{A}_{(A_3{A}_4)})/\pi }$, such that the larger survey area is used.

4. Harmonic space covariance

The components outlined in Section 3 can now be translated to harmonic space through the Limber projection, as given in Equation (10). To achieve this, the line-of-sight weights, $W_{a_i}(\chi )$, still have to be specified. The OneCovariance code typically employs two types of tracers, a_i, which can represent either the CS measured from a source galaxy sample, denoted m_i or the positions of a lens galaxy sample g_i, respectively. In the literature, different terminologies exist. With a ‘source’ we refer to a galaxy whose lensed ellipticity will be used. Conversely, a ‘lens’ is a galaxy whose position will be used. To summarise, GC would be the ‘lens’-‘lens’ auto-correlation, CS the ‘source’-‘source’ auto-correlation and GGL the ‘source’-‘lens’ cross-correlation. The index i labels the tomographic bin of the corresponding sample. With this, the weight functions assume the following form

$$W_{{\mathrm{g}}_i}(\chi )=\hspace{0.5em}{n}_{\mathrm{l}}^{(i)}(\chi ),$$(35)

$$W_{{\mathrm{m}}_i}(\chi )=\frac{3f_{\mathrm{K}}(\chi ){\mathrm{\Omega }}_{\mathrm{m}}}{2{\chi }_{\mathrm{H}}^{2}a}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}}\mathrm{d}{\chi }^{\prime }\hspace{0.5em}\frac{f_{\mathrm{K}}({\chi }^{\prime }-\chi )}{f_{\mathrm{K}}({\chi }^{\prime })}{n}_{\mathrm{s}}^{(i)}(\chi )-A_{\mathrm{IA}}{\left[\frac{1+z(\chi )}{1+z_{\mathrm{pivot}}}\right]}^{{\eta }_{\mathrm{IA}}}\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D_{+}[a(\chi )]}{n}_{\mathrm{s}}^{(i)}(\chi ).$$(36)

Here χ_H is the co-moving distance to the horizon, $n_{\mathrm{l}}^{(i)}(\chi )$ and $n_{\mathrm{s}}^{(i)}(\chi )$ are the normalised redshift distributions of the lens- and source sample respectively. D₊(a) is the linear growth factor and C₁ρ_cr≈0.0134. The second term in the $W_{{\mathrm{m}}_i}$ weight function corresponds to the non-linear linear alignment model (NLA, Hirata & Seljak 2004; Bridle & King 2007) with alignment amplitude A_IA and redshift dependence η_IA at a pivot redshift z_pivot. There are two remarks in order here: (i) The NLA model is the fiducial choice of the OneCovariance code. However, the user can include any alignment model (at least in the dominating Gaussian part) by providing the corresponding C_ℓ as input. Therefore, including, for example, the Tidal Alignment-Tidal Torquing (TATT, Blazek et al. 2019) is straightforwardly done (see Section 6). (ii) For the fiducial NLA model, any lensing tracer automatically assumes a linear response to the non-linearly evolved density field. Consequently, the alignment is also modelled in the SSC and NG terms of the covariance matrix and not just in the Gaussian part.

In addition to the cosmological signal ($a_{\ell m}^i$), observed modes (${\hat{a}}_{\ell m}^i$) contain a noise realisation ($n_{\ell m}^i$) since the continuous field is sampled by a finite number of tracers (galaxies),

$${\hat{a}}_{i,\ell m}=a_{i,\ell m}+n_{i,\ell m},$$(37)

with $\langle a_{i,\ell m}{n}_{j,{\ell }^{\prime }{m}^{\prime }}\rangle =0$ and noise statistic $\langle n_{i,\ell m}{n}_{j,{\ell }^{\prime }{m}^{\prime }}\rangle ={\mathcal{N}}_i^{}{\delta }_{\mathit{ij}}^{\mathrm{K}}$, as the noise in different maps is uncorrelated, and the noise is scale independent. Strictly speaking, this relation only holds for full sky coverage where all the modes are independent. Therefore, the noise component will be treated in real space directly to avoid this complication. We come back to this issue in Section 5. The factor ${\mathcal{N}}_i$ determines the overall noise level depending on the observable under consideration:

$${\mathcal{N}}_i=\{\begin{array}{cc}\hfill \frac{{\sigma }_{{ϵ}_1,\hspace{0.5em}i}^2}{{\overline{n}}_i}\hfill & \mathrm{if}\hspace{0.5em}i\in \hspace{0.5em}\mathrm{source}\hfill \\ \hfill \frac{1}{{\overline{n}}_i}\hfill & \mathrm{if}\hspace{0.5em}i\in \hspace{0.5em}\mathrm{lens}\hfill \end{array}.$$(38)

Here ${\overline{n}}_i$ is the average number density of the tracers and ${\sigma }_{{ϵ}_1}^2$ is the single component ellipticity dispersion of the sources. The standard idealised angular power spectrum estimator at each multipole is defined as

$${\hat{C}}_{a_1{a}_2}(\ell )≔\frac{1}{(2\ell +1)}{\displaystyle \sum _{m=-\ell }^{\ell }}{\hat{a}}_{\ell m}^{1*}{\hat{a}}_{\ell m}^2,$$(39)

such that

$$\langle {\hat{C}}_{a_1{a}_2}(\ell )\rangle ={\mathcal{P}}_{a_1{a}_2}(\ell )+{\mathcal{N}}_{a_1}{\delta }_{a_1{a}_2}^{\mathrm{K}}\equiv C_{a_1{a}_2}(\ell )\hspace{0.5em}.$$(40)

In the flat sky approximation, the same estimator assumes the form

$${\hat{C}}_{a_1{a}_2}(\ell )≔\frac{1}{A_{\mathrm{s}}{N}_{\mathrm{modes}}(\ell )}{\displaystyle \sum _{\ell \in {\ell }_{\mathrm{ring}}}}{\hat{a}}_{\ell }^1{\hat{a}}_{-\ell }^2,$$(41)

where we defined an annulus in Fourier space with volume 2πℓΔℓ, so that the number of available modes over the survey area A_s is

$$N_{\mathrm{modes}}(\ell )=\frac{2\pi \ell \mathrm{\Delta }\ell }{(2\pi {)}^2/A_{\mathrm{s}}}=2\ell \mathrm{\Delta }\ell f_{\mathrm{sky}},$$(42)

where Δℓ≪ℓ was assumed. The covariance matrix between two estimators is then, again by noting that ℓ₁ corresponds to the pair a₁,a₂ and ℓ₂ corresponds to the pair a₃,a₄:

$$\begin{array}{cc}\hfill \mathrm{Cov}[{\hat{C}}_{a_1{a}_2}({\ell }_1){\hat{C}}_{a_3{a}_4}({\ell }_2)]=& \hspace{0.5em}\frac{1}{A_{(12)}{N}_{\mathrm{modes}}({\ell }_1)}\frac{1}{A_{(34)}{N}_{\mathrm{modes}}({\ell }_2)}\sum _{{\tilde{\ell }}_1\in {\ell }_{1,\mathrm{ring}}}\sum _{{\tilde{\ell }}_2\in {\ell }_{2,\mathrm{ring}}}\left(\langle {\hat{a}}_{{\ell }_1}^1{\hat{a}}_{-{\ell }_1}^2{\hat{a}}_{{\ell }_2}^3{\hat{a}}_{-{\ell }_2}^4\rangle -\langle {\hat{a}}_{{\ell }_1}^1{\hat{a}}_{-{\ell }_1}^2\rangle \langle {\hat{a}}_{{\ell }_2}^3{\hat{a}}_{-{\ell }_2}^4\rangle \right)\hfill \\ \hfill =& \hspace{0.5em}\frac{1}{N_{\mathrm{modes}}({\ell }_1)}\frac{1}{N_{\mathrm{modes}}({\ell }_2)}\sum _{{\tilde{\ell }}_1\in {\ell }_{1,\mathrm{ring}}}\left(C_{a_1{a}_3}({\ell }_1)C_{a_2{a}_4}({\ell }_2)+C_{a_1{a}_4}({\ell }_1)C_{a_2{a}_3}({\ell }_2)\right){\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}\hfill \\ \hfill & +\hspace{0.5em}\frac{1}{\max (A_{(12)}{A}_{(34)})N_{\mathrm{modes}}({\ell }_1)N_{\mathrm{modes}}({\ell }_2)}\sum _{{\tilde{\ell }}_1\in {\ell }_{1,\mathrm{ring}}}\sum _{{\tilde{\ell }}_2\in {\ell }_{2,\mathrm{ring}}}{T}^{(a_1{a}_2)(a_3{a}_4)}({\ell }_1,-{\ell }_1,{\ell }_2,-{\ell }_2).\hfill \end{array}$$(43)

Here A_ij,s is the survey area over which the angular power spectrum of the two fields a_i and a_j is estimated. Furthermore, T^{(a₁a₂)(a₃a₄)}(ℓ₁,−ℓ₁,ℓ₂,−ℓ₂) is the trispectrum, see Equation (49) below. We note that the area cancels out in the Gaussian term due to the application of two correlators (compare Equation (D.3)).

Assuming that the angular power spectra do not vary significantly over the bandwidth of the ring ℓ_ring, the Gaussian term in the flat sky approximation is given by the commonly used expression

$${\mathrm{Cov}}_{\mathrm{G}}[{\hat{C}}_{a_1{a}_2}({\ell }_1){\hat{C}}_{a_3{a}_4}({\ell }_2)]\approx \frac{4\pi {\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}}{\max (A_{(12)}{A}_{(34)})2{\ell }_1\mathrm{\Delta }{\ell }_1}\left[C_{a_1{a}_3}({\ell }_1)C_{a_2{a}_4}({\ell }_2)+C_{a_1{a}_4}({\ell }_1)C_{a_2{a}_3}({\ell }_2)\right].$$(44)

This can be compared to the full sky version of the Gaussian term which, at each multipole, is given by

$${\mathrm{Cov}}_{\mathrm{G}}[{\hat{C}}_{a_1{a}_2}({\ell }_1){\hat{C}}_{a_3{a}_4}({\ell }_2)]=\frac{4\pi {\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}}{\max (A_{(12)}{A}_{(34)})(2{\ell }_1+1)}\left[C_{a_1{a}_3}({\ell }_1)C_{a_2{a}_4}({\ell }_2)+C_{a_1{a}_4}({\ell }_1)C_{a_2{a}_3}({\ell }_2)\right],$$(45)

where we introduced the sky fraction max(A₍₁₂₎A₍₃₄₎)/(4π) to mimic incomplete sky coverage (van Uitert et al. 2018). Averaging over multipole bands and assuming again that the spectra do not vary significantly across the bands yields

$${\mathrm{Cov}}_{\mathrm{G}}[{\hat{C}}_{a_1{a}_2}({\ell }_1){\hat{C}}_{a_3{a}_4}({\ell }_2)]\approx \frac{4\pi {\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}}{\max (A_{(12)}{A}_{(34)})(2{\ell }_1+1)\mathrm{\Delta }{\ell }_1}\left[C_{a_1{a}_3}({\ell }_1)C_{a_2{a}_4}({\ell }_2)+C_{a_1{a}_4}({\ell }_1)C_{a_2{a}_3}({\ell }_2)\right],$$(46)

which is the same expression as Equation (44) for ℓ≫1, as expected. Equation (45) amounts to the splitting discussed in Equation (6). We now turn back to the connected term in Equation (43):

$${\mathrm{Cov}}_{\mathrm{nG}}[{\hat{C}}_{a_1{a}_2}({\ell }_1){\hat{C}}_{a_3{a}_4}({\ell }_2)]=\frac{1}{\max (A_{(12)}{A}_{(34)})}{\displaystyle {\int }_{{\tilde{\ell }}_1\in {\ell }_1}}\frac{{\mathrm{d}}^2{\tilde{\ell }}_1}{A({\ell }_1)}{\displaystyle {\int }_{{\tilde{\ell }}_2\in {\ell }_2}}\frac{{\mathrm{d}}^2{\tilde{\ell }}_2}{A({\ell }_2)}{T}^{(a_1{a}_2)(a_3{a}_4)}({\tilde{\ell }}_1,-{\tilde{\ell }}_1,{\tilde{\ell }}_2,-{\tilde{\ell }}_2),$$(47)

where we moved to the continuous version of the two-dimensional Fourier transformation. All trispectra are projected with the Limber projection and the angular average is carried out already in k-space:

$$T_{{\ell }_1{\ell }_2}^{(a_1{a}_2)(a_3{a}_4),\hspace{0.5em}X}=\int \frac{\mathrm{d}\chi }{f_{\mathrm{K}}(\chi {)}^6}{W}_{a_1}(\chi )\dots W_{a_4}(\chi ){\overline{T}}_{(A_1{A}_2)(A_3{A}_4)}^X\left(\frac{{\ell }_1+0.5}{f_{\mathrm{K}}(\chi )},\frac{{\ell }_2+0.5}{f_{\mathrm{K}}(\chi )},\chi \right),$$(48)

where

$${\overline{T}}_{(A_1{A}_2)(A_3{A}_4)}^X(k_1,k_2,z(\chi ))=\{\begin{array}{cc}{T}_{(A_1{A}_2)(A_1{A}_2)}^{\mathrm{SSC}}(k_1,k_2,z(\chi ))\hfill & \text{if}\hspace{0.5em}X=\mathrm{SSC},\hfill \\ {\int }_0^{\pi }\frac{\mathrm{d}{\varphi }_{\ell }}{\pi }{T}_{A_1\dots A_4}(k_1,-k_1,k_2,-k_2,z(\chi ))\hfill & \text{if}\hspace{0.5em}X=\mathrm{nG},\hfill \end{array}$$(49)

where ϕ_ℓ is the angle between the two wave vectors k_1,2 in the plane perpendicular to the line of sight. Altogether, the covariance in harmonic space is given by

$$\mathrm{Cov}[{\hat{C}}_{a_1{a}_2}({\ell }_1){\hat{C}}_{a_3{a}_4}({\ell }_2)]={\mathrm{Cov}}_{\mathrm{G}}\left[{\hat{C}}_{a_1{a}_2}({\ell }_1){\hat{C}}_{a_3{a}_4}({\ell }_2)\right]+\frac{1}{\max (A_{(12)}{A}_{(34)})}{T}_{{\ell }_1{\ell }_2}^{(a_1{a}_2)(a_3{a}_4),\hspace{0.5em}\mathrm{nG}}+T_{{\ell }_1{\ell }_2}^{(a_1{a}_2)(a_3{a}_4),\hspace{0.5em}\mathrm{SSC}}\hspace{0.5em}$$(50)

at each multipole ℓ. For an ℓ-band averaged version, we use Equations (44) and (47). This is only important if a pure ℓ-space covariance is required, since the averaging over scales is taken care of when mapping the theoretical ℓ-space covariance to the observables in the next section.

5. From harmonic space to observables

Whilst theoretical modelling is most straightforward in harmonic space, the theoretical power spectra, denoted as C(ℓ), are technically defined solely on the full sky. As real surveys typically cover only a portion of the sky, the advantage of harmonic space diminishes. This is because different ℓ modes become intermingled, and the partial sky C_ℓ, often referred to as pseudo-C_ℓ, are derived from the full-sky C_ℓ via convolution with the mode mixing matrix (see e.g. Hivon et al. 2002; Kogut et al. 2003; Reinecke & Seljebotn 2013; Elsner et al. 2017; Alonso et al. 2019).

Deconvolution may not always be feasible for intricate sky masks, rendering the C_ℓ indirectly observable and serving primarily as an intermediary product. Another challenge arises from the projection, as outlined in Equation (7), which combines various physical scales if the window function W_A(χ) possesses broad support in the co-moving distance. Additionally, the shear field's spin-2 nature permits its separation into E and B modes.

To mitigate these difficulties, a variety of summary statistics are available, including the probes of the classical 3×2-point analysis: GC, CS, and GGL. Accordingly, the OneCovariance code can accommodate diverse inputs for three-dimensional power spectra, angular power spectra, or weight functions for the projection, rendering the discussion applicable to all general tracers or summary statistics of two-point functions in the LSS. In a broader context, any two-point summary statistic represents a linear transformation of the underlying two-point statistic in harmonic space, as nonlinear transformations would involve contributions from higher-order cumulants. Therefore, we write the summary statistic $\text{O}$ between two tracers a₁ and a₂ as

$$\text{O}(L)=\int \frac{\ell \mathrm{d}\ell }{2\pi }\hspace{0.5em}{\text{W}}_L(\ell )\mathrm{vec}({\text{C}}_{a_1{a}_2})(\ell ).$$(51)

Here L can be a discrete label or a continuous variable for some Fourier filter W_L(ℓ). Importantly, $\mathrm{vec}({\text{C}}_{a_1{a}_2})$, concatenates all unique ${\text{C}}_{a_1{a}_2}$ into a vector which is mapped to the new summary statistic $\text{O}$ via the aforementioned linear map. Throughout this section, we omit the formal integration limits from zero to infinity, but all integrals over ℓ are understood to run over this range if not specified differently. It should be noted that the internal structure of ${\text{C}}_{a_1{a}_2}$ depends on the underlying fields a₁ and a₂. If the spin of a₁ and a₂ is non-zero, such as for CS, they both consist of two complex components so that C is a real 2×2 matrix given by

$${\text{C}}_{a_1{a}_2}({\ell }_1){\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}{\delta }_{m_1{m}_2}^{\mathrm{K}}=\langle a_{1,{\ell }_1{m}_1}^{\phantom{†}}{a}_{2,{\ell }_2{m}_2}^{†}\rangle ,$$(52)

on the full sky for a homogeneous and isotropic random field. The CS submatrix of W_L(ℓ) is a symmetric and real 2×2 matrix with two independent contributions: the curl-free E-mode signal and the divergence-free B-mode signal, where we assume that there is no measurable contribution to the EB spectrum. We thus bundle the combined C(ℓ) of all probes (for a standard 3×2 analysis) in a vector by using the notation of Equation 35 and 36

$$[\mathrm{vec}(C)](\ell )≔(C_{\mathrm{gg}},C_{\mathrm{gm}},C_{\mathrm{mm},\mathrm{EE}},C_{\mathrm{mm},\mathrm{BB}}{)}^{\mathrm{T}}(\ell ),$$(53)

where tomographic bin indices have been omitted for clarity. The resulting summary statistics are bundled in the same fashion:

$$[\mathrm{vec}(\text{O})](L)≔({\mathcal{O}}_{\mathrm{gg}},{\mathcal{O}}_{\mathrm{gm}},{\mathcal{O}}_{\mathrm{mm},\mathrm{EE}},{\mathcal{O}}_{\mathrm{mm},\mathrm{BB}}{)}^{\mathrm{T}}(L).$$(54)

In this manner, we can express the transformation from Fourier space to the summary statistic using Equation (51). Real space correlation functions decouple the shot or shape noise contribution from different scales, L. Furthermore, the shot or shape noise contribution can be precisely estimated from the data itself, regardless of complex masks. To leverage those properties the shot or shape noise levels defined in Equation (38), which only apply to full-sky data, are improved by also providing a relation between real space correlation (see Section 5.2) and the observables. Therefore, we also require the following correlation for some real space filter function R_L(θ)

$$\mathrm{vec}(\text{O})(L)=\int \theta \mathrm{d}\theta \hspace{0.5em}{\text{R}}_L(\theta )\left(\begin{array}{c}\hfill w(\theta )\hfill \\ \hfill {\gamma }_t(\theta )\hfill \\ \hfill {\xi }_{+}(\theta )\hfill \\ \hfill {\xi }_{-}(\theta )\hfill \end{array}\right),$$(55)

with the real space correlation functions defined via Hankel transformations of the angular power spectrum, Equation (8). For CS, this amounts to the two correlation functions (Bartelmann & Schneider 2001)

$${\xi }_{+}^{(\mathit{ij})}(\theta )=\int \frac{\ell \mathrm{d}\ell }{2\pi }\hspace{0.5em}\left[{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{E}}(\ell )+{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{B}}(\ell )\right]\hspace{0.5em}{\mathrm{J}}_0(\ell \theta ),$$(56)

$${\xi }_{-}^{(\mathit{ij})}(\theta )=\int \frac{\ell \mathrm{d}\ell }{2\pi }\hspace{0.5em}\left[{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{E}}(\ell )-{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{B}}(\ell )\right]\hspace{0.5em}{\mathrm{J}}_4(\ell \theta ),$$(57)

where J_μ is a cylindrical Bessel function of the first kind of order μ. The angular power spectrum ${\mathcal{P}}_{a_1{a}_2}(\ell )$ was decomposed into an E and B-mode component. Recall that we explicitly use the noise-free version of the angular power spectra here as defined in the estimator in Equation (40). For GGL one measures the tangential shear correlation function

$${\gamma }_{\mathrm{t}}^{(\mathit{ij})}(\theta )=\int \frac{\ell \mathrm{d}\ell }{2\pi }{\mathcal{P}}_{{\mathrm{n}}_i{ϵ}_j}(\ell )\hspace{0.5em}{\mathrm{J}}_2(\ell \theta ).$$(58)

The galaxy correlation function can be calculated as

$$w^{(\mathit{ij})}(\theta )=\int \frac{\ell \mathrm{d}\ell }{2\pi }\hspace{0.5em}{\mathcal{P}}_{{\mathrm{n}}_i{\mathrm{n}}_j}(\ell )\hspace{0.5em}{\mathrm{J}}_0(\ell \theta ).$$(59)

We note that these are all the continuous versions and large multipole approximations of the discrete transformations. We continue to review the definitions of the three most commonly used summary statistics (real space correlation functions, band powers and COSEBIs) and the corresponding covariances. A more detailed discussion of the general setting can be found in Appendix I.

5.1. Multiplicative shear bias

Shear measurements for source samples are calibrated on image simulations (see e.g. Kannawadi et al. 2019; Li et al. 2023). Residual biases are captured in a multiplicative and additive shear bias, which needs to be propagated into the cosmological inference. We assume that there are no residual spatial patterns in the additive contribution (e.g. from point-spread function leakage) and hence no correlation with the cosmological signal. Thus, the only source of error which has to be propagated in the covariance is the multiplicative m-correction. The residual error on the multiplicative shear bias after calibration is labelled ${\sigma }_{\mathrm{m}}^{a_1}$ for source bin a₁. Considering the multiplicative shear bias uncertainty as an additive contribution to the covariance matrix (Joachimi et al. 2021), so that

$${\mathrm{Cov}}_{\mathrm{mult}}\left[{\mathcal{O}}_{a_1{a}_2}(L_1),{\mathcal{O}}_{a_3{a}_4}(L_2)\right]={\mathcal{O}}_{a_1{a}_2}(L_1){\mathcal{O}}_{a_3{a}_4}(L_2)\left[{\sigma }_{\mathrm{m}}^{a_1}{\sigma }_{\mathrm{m}}^{a_3}+{\sigma }_{\mathrm{m}}^{a_1}{\sigma }_{\mathrm{m}}^{a_4}+{\sigma }_{\mathrm{m}}^{a_2}{\sigma }_{\mathrm{m}}^{a_3}+{\sigma }_{\mathrm{m}}^{a_2}{\sigma }_{\mathrm{m}}^{a_4}\right].$$(60)

This approximation holds for ${\sigma }_{\mathrm{m}}^{a_1}\ll 1$, for 〈m〉=0〉 after correction and fully correlated m-bias corrections across tomographic bins. We note that only source samples have a ${\sigma }_{\mathrm{m}}^{a_1}$. by definition. The assumption 〈m〉=0〉 is not crucial for the covariance, as this can be addressed by directly shifting the data vector while keeping the covariance unchanged (Li et al. 2023).

5.2. Real space correlation functions

When propagating the covariance matrix from the angular power spectra to the real space correlation functions in Equations (56), (57), (58), and (59) we take into account that the correlation functions are measured over finite angular separation bins centred at ${\overline{\theta }}_i$ with boundaries [θ_l,i,θ_u,i], the bin average needs to be carried out explicitly. To capture the different weight functions, we use the notation from Joachimi et al. (2021)

$${\mathcal{K}}_{\mu }\left(\ell {\overline{\theta }}_i\right):=\frac{2}{{\theta }_{\mathrm{u},i}^2-{\theta }_{\mathrm{l},i}^2}{\displaystyle {\int }_{{\theta }_{\mathrm{l},i}}^{{\theta }_{\mathrm{u},i}}}\mathrm{d}\theta \prime \theta \prime {\mathrm{J}}_{\mu }(\ell \theta \prime )=\frac{2}{\left({\theta }_{\mathrm{u},i}^2-{\theta }_{\mathrm{l},i}^2\right){\ell }^2}\times \{\begin{array}{cc}{\left[x{\mathrm{J}}_1(x)\right]}_{\ell {\theta }_{\mathrm{l},i}}^{\ell {\theta }_{\mathrm{u},i}}\hfill & \mu =0,\hfill \\ {\left[-x{\mathrm{J}}_1(x)-2{\mathrm{J}}_0(x)\right]}_{\ell {\theta }_{\mathrm{l},i}}^{\ell {\theta }_{\mathrm{u},i}}\hfill & \mu =2,\hfill \\ {\left[\left(x-\frac{8}{x}\right){\mathrm{J}}_1(x)-8{\mathrm{J}}_2(x)\right]}_{\ell {\theta }_{\mathrm{l},i}}^{\ell {\theta }_{\mathrm{u},i}}\hfill & \mu =4.\hfill \end{array}\hspace{0.5em},$$(61)

ignoring the weights of the pair counts (Asgari et al. 2019). The set of possible correlation functions is denoted accordingly in a compact form (see the notation used in Joachimi et al. 2021):

$$\left\{w,\langle {\gamma }_{\mathrm{t}}\rangle ,{\xi }_{+},{\xi }_{-}\right\}\leftrightarrow \left\{{\mathrm{\Xi }}_0,{\mathrm{\Xi }}_2,{\mathrm{\Xi }}_0,{\mathrm{\Xi }}_4\right\}.$$

We note that w and ξ₊ have the same weight function. As discussed, the Gaussian covariance is split into three contributions as in Equation (6). For the pure shot/shape noise (sn) term, this procedure is a numerical necessity, as an ℓ-independent integrand will lead to a δ_D-contribution to the covariance. Since we reconsider the idealised term in Appendix F for real space correlation functions, this splitting is also necessary here as the ‘mix’ term is the most important while also still being a somewhat uncertain term in the analytical covariance modelling. Numerically, it is, however, usually advantageous to combine the sample variance (‘sva’) and mixed terms to speed up the convergence of the integrals.

Using Equation (51), the definitions of the weights above and Equation (50), the ‘sva’ term is given by

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sva}}\left[{\mathrm{\Xi }}_{\mu }^{(\mathit{ij})}\left({\overline{\theta }}_p\right);{\mathrm{\Xi }}_{\nu }^{(\mathit{mn})}\left({\overline{\theta }}_q\right)\right]=\frac{1}{2\pi \hspace{0.5em}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}\ell \ell \hspace{0.5em}{\mathcal{K}}_{\mu }\left(\ell {\overline{\theta }}_p\right){\mathcal{K}}_{\nu }\left(\ell {\overline{\theta }}_q\right)\left[{\mathcal{P}}_{\mathit{im}}(\ell ){\mathcal{P}}_{\mathit{jn}}(\ell )+{\mathcal{P}}_{\mathit{in}}(\ell ){\mathcal{P}}_{\mathit{jm}}(\ell )\right],$$(62)

where the tomographic indices i,j,m,n label either a source or a lens sample. We measure all areas from a binary healpix mask down to a certain angular scale so that masked stars and other features only reduce the effective number density and not the survey area. This assumption is valid as long as this dilution scale is smaller than the smallest scale over which the cosmological signal is measured. Typically, one assumes N_side=4096 for the healpix (Górski et al. 2005) evaluation, corresponding to a pixel size of just below one arcmin.

In an ideal survey, the pure noise contribution is just the product of the two individual noise contributions, Equation (38), to the measurement of the two-point statistic amounting to a term proportional to the number of (non-unique) pairs in bin ${\overline{\theta }}_p$:

$$N_{\mathrm{pair},\hspace{0.5em}\mathrm{ideal}}^{(\mathit{ij})}\left({\overline{\theta }}_p\right)=\pi \left({\theta }_{\mathrm{u},p}^2-{\theta }_{\mathrm{l},p}^2\right){\overline{n}}_i{\overline{n}}_j{A}_{(\mathit{ij})}.$$(63)

This equation is accurate in the absence of survey boundaries. For finite areas, however, the number of pairs as a function of angular scale differs from the expected ∝θ² scaling. To obtain a more accurate prescription of the pure shot noise contribution, the number of pairs is directly measured at the catalogue level, providing $N_{\mathrm{pair}}^{(\mathit{ij})}\left({\overline{\theta }}_p\right)$. Altogether,

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sn}}\left[{\mathrm{\Xi }}_{\mu }^{(\mathit{ij})}\left({\overline{\theta }}_p\right);{\mathrm{\Xi }}_{\nu }^{(\mathit{mn})}\left({\overline{\theta }}_q\right)\right]={\delta }_{\mu \nu }^{\mathrm{K}}{\delta }_{{\overline{\theta }}_p{\overline{\theta }}_q}^{\mathrm{K}}\left({\delta }_{\mathit{im}}^{\mathrm{K}}{\delta }_{\mathit{jn}}^{\mathrm{K}}+{\delta }_{\mathit{in}}^{\mathrm{K}}{\delta }_{\mathit{jm}}^{\mathrm{K}}\right)\frac{{\mathcal{T}}_{(\mathit{ij})(\mathit{mn})}^{\mathrm{sn}}}{N_{\mathrm{pair}}^{(\mathit{ij})}\left({\overline{\theta }}_p\right)},$$(64)

where we defined the noise levels

$${\mathcal{T}}_{(\mathit{ij})(\mathit{mn})}^{\mathrm{sn}}≔\{\begin{array}{cc}2{\sigma }_{{ϵ}_1,\hspace{0.5em}i}^2{\sigma }_{{ϵ}_1,\hspace{0.5em}j}^2\hfill & \text{if}\hspace{0.5em}i,j,m,n\in \text{source},\hfill \\ {\sigma }_{{ϵ}_1,\hspace{0.5em}i}^2\hfill & \text{if both}\hspace{0.5em}(\mathit{ij})\hspace{0.5em}\text{and}\hspace{0.5em}(\mathit{mn})\hspace{0.5em}\text{contain a lens and source sample each},\hfill \\ 1\hfill & \text{if}\hspace{0.5em}i,j,m,n\in \text{lens}.\hfill \end{array}$$(65)

We recall that ‘source’ refers to the galaxy shape and ‘lens’ to the galaxy position being used. Furthermore, we note that the Latin indices carry a hidden label “lens” or “source” so that ${\delta }_{\mathit{ij}}^{\mathrm{K}}=1$ only when both indices come from the same set. It should be noted that, while we use in the equations here the non-unique pairs, for any correlation measurement, however, only the unique pairs are relevant. This is accounted for in the covariance calculations. Since this often leads to confusion, let us use tomographic clustering with equi-populated bins as an example. The noise level will always be ∼1/N², with N being the number of galaxies in each bin. For auto-correlations, however, the Kronecker symbols in the brackets in Equation (64) make sure that only unique pairs are counted, effectively increases the covariance by a factor of two. In contrast, for cross-correlations, the N² is actually the number of unique pairs.

The mixed term is calculated in the same fashion as the sample variance contribution, Equation (62):

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{mix}}\left[{\mathrm{\Xi }}_{\mu }^{(\mathit{ij})}\left({\overline{\theta }}_p\right);{\mathrm{\Xi }}_{\nu }^{(\mathit{mn})}\left({\overline{\theta }}_q\right)\right]={\delta }_{\mathit{jn}}\frac{{\mathcal{T}}_j^{\mathrm{mix}}}{2\pi n_{\mathrm{eff}}^{(j)}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}\ell \ell \hspace{0.5em}{\mathcal{K}}_{\mu }\left(\ell {\overline{\theta }}_p\right){\mathcal{K}}_{\nu }\left(\ell {\overline{\theta }}_q\right){\mathcal{P}}_{\mathit{im}}(\ell )+3\hspace{0.5em}\mathrm{perm}.,$$(66)

The noise level of the mixed term is defined as

$${\mathcal{T}}_j^{\mathrm{mix}}≔\{\begin{array}{cc}{\sigma }_{{ϵ}_1,j}^2\hfill & \text{if}\hspace{0.5em}j\in \text{source},\hfill \\ 1\hfill & \text{if}\hspace{0.5em}j\in \text{lens}.\hfill \end{array}$$(67)

We note that this is an idealised setting for a rescaled full sky survey without a complicated mask. In Appendix F we revisit this term using the triplet counts of a catalogue to estimate the effect on the covariance and cosmological inference.

The non-Gaussian covariance and SSC are

$${\mathrm{Cov}}_{\mathrm{NG}/\mathrm{SSC}}\left[{\mathrm{\Xi }}_{\mu }^{(\mathit{ij})}\left({\overline{\theta }}_p\right);{\mathrm{\Xi }}_{\nu }^{(\mathit{mn})}\left({\overline{\theta }}_q\right)\right]=\hspace{0.5em}\frac{1}{4{\pi }^2}{\displaystyle \int }\mathrm{d}{\ell }_1{\ell }_1\hspace{0.5em}{\mathcal{K}}_{\mu }\left({\ell }_1{\overline{\theta }}_p\right){\displaystyle \int }\mathrm{d}{\ell }_2{\ell }_2\hspace{0.5em}{\mathcal{K}}_{\nu }\left({\ell }_2{\overline{\theta }}_q\right)$$(68)

$$\times \{\begin{array}{cc}\frac{1}{\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}{\int }_0^{\pi }\frac{\mathrm{d}{\varphi }_{\ell }}{\pi }\hspace{0.5em}{T}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,-{\ell }_1,{\ell }_2,-{\ell }_2)\hfill & \mathrm{for}\hspace{0.5em}\mathrm{nG}\hfill \\ T_{\mathrm{SSC}}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,{\ell }_2)\hfill & \mathrm{for}\hspace{0.5em}\mathrm{SSC}\hfill \end{array}\hspace{0.5em}.$$

5.3. Bandpowers

We define the bandpower signal their weights and further ingredients in Appendix G and focus on the covariance here. In contrast to the real space approach in Joachimi et al. (2021) we calculate the bandpower covariance from Fourier space directly using Equations (G.13), (G.14), (G.15), and (G.16). This is done to resemble the analytical calculations of the predicted signal and for numerical stability and speed. In complete analogy to the real space correlation functions, we write

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sva}}\left[{\mathcal{C}}_{\mu }^{(\mathit{ij})}(L_1),{\mathcal{C}}_{\nu }^{(\mathit{mn})}(L_2)\right]=\frac{2\pi }{{\mathcal{N}}_{L_1}{\mathcal{N}}_{L_2}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}\ell \ell \hspace{0.5em}{\mathcal{W}}_{\mu }^{L_1}(\ell ){\mathcal{W}}_{\nu }^{L_2}(\ell )[{\mathcal{P}}_{\mathit{im}}(\ell ){\mathcal{P}}_{\mathit{jn}}(\ell )+{\mathcal{P}}_{\mathit{in}}(\ell ){\mathcal{P}}_{\mathit{jm}}(\ell )]$$(69)

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{mix}}\left[{\mathcal{C}}_{\mu }^{(\mathit{ij})}(L_1),{\mathcal{C}}_{\nu }^{(\mathit{mn})}(L_2)\right]={\delta }_{\mathit{jn}}\frac{2\pi \hspace{0.5em}{\mathcal{T}}_j^{\mathrm{mix}}}{{\mathcal{N}}_{L_1}{\mathcal{N}}_{L_2}{n}_{\mathrm{eff}}^{(j)}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}\ell \ell \hspace{0.5em}{\mathcal{W}}_{\mu }^{L_1}(\ell ){\mathcal{W}}_{\nu }^{L_2}(\ell ){\mathcal{P}}_{\mathit{im}}(\ell )+3\hspace{0.5em}\mathrm{perm}.,$$

assuming that the cosmological B-mode signal vanishes. The following shorthand notation for the weights was used:

$${\mathcal{W}}_{\mu }^L(\ell )=\{\begin{array}{cc}{W}_{\mathit{EE}}^L(\ell )/2\hfill & \text{if}\mu =ϵϵ\text{E}\hfill \\ W_{\mathit{EE}}^L(\ell )\hfill & \text{if}\mu =\hspace{0.5em}\text{nn}\hfill \\ W_{\mathit{nE}}^L(\ell )\hfill & \text{if}\mu =\mathrm{n}ϵ\hfill \\ W_{\mathit{EB}}^L(\ell )\hfill & \text{if}\mu =ϵϵ\text{B}\hfill \end{array}.$$(70)

To properly incorporate the pair counts the pure shot-noise contributions are calculated from real space:

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sn}}\left[{\mathcal{C}}_{\mu }^{(\mathit{ij})}(L_1),{\mathcal{C}}_{\nu }^{(\mathit{mn})}(L_2)\right]=\hspace{0.5em}\frac{{\pi }^2{\delta }_{\mu \nu }^{\mathrm{K}}}{{\mathcal{N}}_{L_1}{\mathcal{N}}_{L_2}}\left({\delta }_{\mathit{im}}^{\mathrm{K}}{\delta }_{\mathit{jn}}^{\mathrm{K}}+{\delta }_{\mathit{in}}^{\mathrm{K}}{\delta }_{\mathit{jm}}^{\mathrm{K}}\right)\int \frac{{\theta }^2\mathrm{d}\theta }{n_{\mathrm{pair}}^{(\mathit{ij})}(\theta )}$$(71)

$$\times \hspace{0.5em}\{\begin{array}{cc}2{\sigma }_{{ϵ}_1,\hspace{0.5em}i}^2{\sigma }_{{ϵ}_1,\hspace{0.5em}j}^2\left(g_{+}^{L_1}(\theta )g_{+}^{L_2}(\theta )+g_{-}^{L_1}(\theta )g_{-}^{L_2}(\theta )\right)\hfill & \hspace{0.5em}\text{if}\hspace{0.5em}\mu =ϵϵ\hfill \\ 4{\sigma }_{{ϵ}_1,\hspace{0.5em}i}^2{h}^{L_1}(\theta )h^{L_2}(\theta )\hfill & \hspace{0.5em}\text{if}\hspace{0.5em}\mu =\mathrm{n}ϵ\hfill \\ g_{+}^{L_1}(\theta )g_{+}^{L_2}(\theta )\hfill & \hspace{0.5em}\text{if}\hspace{0.5em}\mu =\mathrm{nn}\hfill \\ \hfill & \hfill \end{array}.$$

The differential pair counts, $n_{\mathrm{pair}}^{(\mathit{ij})}(\theta )$, are defined such that the number of pairs, $N_{\mathrm{pair}}^{(\mathit{ij})}$, in an angular bin [θ_l,θ_u] is

$$N_{\mathrm{pair}}^{(\mathit{ij})}={\displaystyle {\int }_{{\theta }_{\mathrm{l}}}^{{\theta }_{\mathrm{u}}}}\mathrm{d}{\theta }^{\prime }\hspace{0.5em}{n}_{\mathrm{pair}}^{(\mathit{ij})}({\theta }^{\prime })$$(72)

directly from the, possibly weighted, pair counts of the catalogue. Lastly, the nG and SSC terms can be calculated as follows:

$${\mathrm{Cov}}_{\mathrm{NG}}\left[{\mathcal{C}}_{\mu }^{(\mathit{ij})}(L_1),{\mathcal{C}}_{\nu }^{(\mathit{mn})}(L_2)\right]=\frac{1}{{\mathcal{N}}_{L_1}{\mathcal{N}}_{L_2}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}{\ell }_1{\ell }_1{\mathcal{W}}_{\mu }^{L_1}({\ell }_1)\int \mathrm{d}{\ell }_2{\ell }_2\hspace{0.5em}{\mathcal{W}}_{\nu }^{L_2}({\ell }_2)$$(73)

$$\times \hspace{0.5em}{\displaystyle {\int }_0^{\pi }}\frac{\mathrm{d}{\varphi }_{\ell }}{\pi }\hspace{0.5em}{T}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,-{\ell }_1,{\ell }_2,-{\ell }_2),$$

$${\mathrm{Cov}}_{\mathrm{SSC}}\left[{\mathcal{C}}_{\mu }^{(\mathit{ij})}(L_1),{\mathcal{C}}_{\nu }^{(\mathit{mn})}(L_2)\right]=\frac{1}{{\mathcal{N}}_{L_1}{\mathcal{N}}_{L_2}}\int \mathrm{d}{\ell }_1{\ell }_1{\mathcal{W}}_{\mu }^{L_1}({\ell }_1)\int \mathrm{d}{\ell }_2{\ell }_2\hspace{0.5em}{\mathcal{W}}_{\nu }^{L_2}({\ell }_2)T_{\mathrm{SSC}}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,{\ell }_2).$$(74)

5.4. COSEBIs and Ψ statistics

Due to the similar structure of the COSEBIs we repeat the same calculation as for the band power covariance, i.e. projecting all terms from harmonic space except for the pure shot noise term. With the definition of COSEBIs, see Appendix H, this yields the following:

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sva}}\left[E_a^{(\mathit{ij})},E_b^{(\mathit{mn})}\right]=\hspace{0.5em}\frac{1}{2\pi \hspace{0.5em}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}\ell \ell \hspace{0.5em}{W}_a(\ell )W_b(\ell )\left[{\mathcal{P}}_{\mathit{im}}(\ell ){\mathcal{P}}_{\mathit{jn}}(\ell )+{\mathcal{P}}_{\mathit{in}}(\ell ){\mathcal{P}}_{\mathit{jm}}(\ell )\right]$$(75)

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{mix}}\left[E_a^{(\mathit{ij})},E_b^{(\mathit{mn})}\right]=\hspace{0.5em}{\delta }_{\mathit{jn}}\frac{\hspace{0.5em}{\mathcal{T}}_j^{\mathrm{mix}}}{2\pi \hspace{0.5em}{n}_{\mathrm{eff}}^{(j)}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}\ell \ell \hspace{0.5em}{W}_a(\ell )W_b(\ell ){\mathcal{P}}_{\mathit{im}}(\ell )+3\hspace{0.5em}\mathrm{perm}.$$(76)

$${\mathrm{Cov}}_{\mathrm{NG}}\left[E_a^{(\mathit{ij})},E_b^{(\mathit{mn})}\right]=\hspace{0.5em}\frac{1}{2\pi \hspace{0.5em}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\int \mathrm{d}{\ell }_1{\ell }_1{W}_a({\ell }_1)\int \mathrm{d}{\ell }_2{\ell }_2\hspace{0.5em}{W}_b({\ell }_2){\displaystyle {\int }_0^{\pi }}\frac{\mathrm{d}{\varphi }_{\ell }}{\pi }\hspace{0.5em}{T}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,-{\ell }_1,{\ell }_2,-{\ell }_2),$$(77)

$${\mathrm{Cov}}_{\mathrm{SSC}}\left[E_a^{(\mathit{ij})},E_b^{(\mathit{mn})}\right]=\hspace{0.5em}\frac{1}{2\pi \hspace{0.5em}}\int \mathrm{d}{\ell }_1{\ell }_1{W}_a({\ell }_1)\int \mathrm{d}{\ell }_2{\ell }_2\hspace{0.5em}{W}_b({\ell }_2)T_{\mathrm{SSC}}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,{\ell }_2)$$(78)

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sn}}\left[E_a^{(\mathit{ij})},E_b^{(\mathit{mn})}\right]=\hspace{0.5em}\frac{{\sigma }_{{ϵ}_1,\hspace{0.5em}i}^2{\sigma }_{{ϵ}_1,\hspace{0.5em}j}^2}{2}\left({\delta }_{\mathit{im}}^{\mathrm{K}}{\delta }_{\mathit{jn}}^{\mathrm{K}}+{\delta }_{\mathit{in}}^{\mathrm{K}}{\delta }_{\mathit{jm}}^{\mathrm{K}}\right)\int \frac{{\theta }^2\mathrm{d}\theta }{n_{\mathrm{pair}}^{(\mathit{ij})}(\theta )}\left[T_{+a}(\theta )T_{+b}(\theta )+T_{-a}(\theta )T_{-b}(\theta )\right].$$(79)

The expressions of the covariance of Ψ statistics is in full analogy to Equations (75), (76), (77), (78) (79).

5.5. Stellar mass function

To complement the halo model, we also implement the SMF covariance, Equation (20), which was already used in Dvornik et al. (2023). For a flux-limited sample, we consider the standard V_max(M_★)-estimator, where V_max is the maximum volume out to which a galaxy with a given mass can be observed given the limiting magnitude of the survey. The estimator for the SMF is then (see e.g. Smith 2012):

$${\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu })≔\frac{1}{\mathrm{\Delta }{M}_{★}^{\mu }}\frac{1}{V_{\max }^{(i)}(M_{★}^{\mu })}{\displaystyle \sum _{a=1}^{N_{\mathrm{tot}}}}{\mathrm{\Theta }}_{\mathrm{H}}(z_{\mathrm{photo},\mathrm{a}},z_{\mathrm{photo}}^{(i)}){\mathrm{\Theta }}_{\mathrm{H}}(M_{★,a},M_{★}^{\mu }).$$(80)

Here we assumed that the mapping from observed to real stellar mass is very tight and well approximated by a δ distribution. If the relation is less well known, this would amount to an additional integration. The indices i,μ,a label a possible splitting in tomographic bins, the stellar mass bin and the individual galaxy in each bin respectively. Θ_H is the Heaviside function. The auto-correlation consists of a noise term and an SSC contribution (Takada & Bridle 2007; Smith 2012),

$$\mathrm{Cov}\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\hat{\mathrm{\Phi }}}^{(j)}(M_{★}^{\nu })\right]=\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\hat{\mathrm{\Phi }}}^{(j)}(M_{★}^{\nu })\right]}_{\mathrm{sn}}+\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\hat{\mathrm{\Phi }}}^{(j)}(M_{★}^{\nu })\right]}_{\mathrm{SSC}},$$(81)

with the halo occupation variance being neglected as it is subdominant (Smith 2012) and

$$\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\hat{\mathrm{\Phi }}}^{(j)}(M_{★}^{\nu })\right]}_{\mathrm{sn}}=\hspace{0.5em}{\delta }_{\mathit{ij}}^{\mathrm{K}}{\delta }_{\mu \nu }^{\mathrm{K}}\frac{{\mathrm{\Phi }}^{(i)}(M_{★}^{\mu })}{\mathrm{\Delta }{M}_{★}^{\mu }{V}_{\max ,\mu }^{(i)}}$$(82)

$$\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\hat{\mathrm{\Phi }}}^{(j)}(M_{★}^{\nu })\right]}_{\mathrm{SSC}}=\hspace{0.5em}\frac{A^2{f}^{(i)}{f}^{(j)}}{V_{\max ,\mu }^{(i)}{V}_{\max ,\nu }^{(j)}}\int \mathrm{d}\chi \hspace{0.5em}\frac{p^{(i)}(\chi )p^{(j)}(\chi )}{p_{\mathrm{tot}}^2(\chi )}{f}_{\mathrm{k}}^2(\chi ){\sigma }_{\mathrm{bg},A}^2(\chi ){\tilde{\mathrm{\Phi }}}_{\mu }(\chi ){\tilde{\mathrm{\Phi }}}_{\nu }(\chi ),$$(83)

where f⁽ⁱ⁾ is the fraction of galaxies in tomographic bin i, A is the survey area, p⁽ⁱ⁾ the redshift distribution of the i-th tomographic bin and σ_bg,A is given by Equation (33) with the same footprint. Lastly, the quantity ${\tilde{\mathrm{\Phi }}}_{\mu }(\chi )$ is defined as

$${\tilde{\mathrm{\Phi }}}_{\mu }(\chi )≔\int \mathrm{d}M\hspace{0.5em}{n}_{\mathrm{h}}(M,z(\chi ))\mathrm{\Phi }(M_{★}^{\mu }|M)b_{\mathrm{h}}(M,z(\chi )).$$(84)

The SMF is also correlated with the LSS and their cross-variance includes a sample variance and SSC term:

$$\mathrm{Cov}\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\mathcal{O}}_{a_1{a}_2}(L)\right]=\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\mathcal{O}}_{a_1{a}_2}(L)\right]}_{\mathrm{sva}}+\hspace{0.5em}\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\mathcal{O}}_{a_1{a}_2}(L)\right]}_{\mathrm{SSC}}.$$(85)

This is in analogy to Equation (51) with the two contributions reading

$$\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\mathcal{O}}_{a_1{a}_2}(L)\right]}_{\mathrm{sva}}=f^{(i)}\int \mathrm{d}\ell \hspace{0.5em}{W}_L(\ell )\int \mathrm{d}\chi \hspace{0.5em}\frac{p^{(i)}(\chi )}{p_{\mathrm{tot}}(\chi )}\frac{W_{A_1}(\chi )W_{A_2}(\chi )}{{\chi }^2}{B}_{\mathrm{cmm},\mu }\left(\frac{\ell +1/2}{\chi },z(\chi )\right),$$(86)

with the bispectrum of counts and matter for a collapsed triangle configuration given by

$$B_{\mathrm{cmm},\mu }\left(k,z\right)=\hspace{0.5em}\int \mathrm{d}M\hspace{0.5em}{n}_{\mathrm{h}}(M,z)\mathrm{\Phi }(M_{★}^{\mu }|M){\left(\frac{M}{{\overline{\rho }}_{\mathrm{m}}}\right)}^2{\tilde{u}}^2(k|M,z)$$(87)

$$+2P_{\mathrm{lin}}(k,z)\int \mathrm{d}M\hspace{0.5em}{n}_{\mathrm{h}}(M,z)\mathrm{\Phi }(M_{★}^{\mu }|M)\frac{M}{{\overline{\rho }}_{\mathrm{m}}}{b}_{\mathrm{h}}(M,z)\tilde{u}(k|M,z)\int \mathrm{d}M\hspace{0.5em}{n}_{\mathrm{h}}(M,z)b_{\mathrm{h}}(M,z)\tilde{u}(k|M,z).$$

Finally, the SSC term is given by

$$\mathrm{Cov}{\left[{\hat{\mathrm{\Phi }}}^{(i)}(M_{★}^{\mu }),{\mathcal{O}}_{a_1{a}_2}(L)\right]}_{\mathrm{SSC}}=A{f}^{(i)}\int \mathrm{d}\ell \hspace{0.5em}{W}_L(\ell )\int \mathrm{d}\chi \hspace{0.5em}\frac{p^{(i)}(\chi )}{p_{\mathrm{tot}}(\chi )}\frac{W_{A_1}(\chi )W_{A_2}(\chi )}{{\chi }^2}\frac{\partial P_{\mathrm{mm}}}{\partial {\delta }_{\mathrm{bg}}}\left(\frac{\ell +1/2}{\chi },z(\chi )\right){\sigma }_{\mathrm{bg},A}^2(\chi ){\tilde{\mathrm{\Phi }}}_{\mu }(\chi ),$$(88)

with the power spectrum responses from Equation (32).

6. The OneCovariance code

This section aims to provide a brief rationale for the initial development of the code. For an overview of the code's general structure, please see Appendix A. To enhance accessibility, a flowchart illustrating the typical workflow of the OneCovariance code¹¹ is presented in Figure A.1.

Several outstanding public codes are available for computing covariance matrices in harmonic space, notably packages like ccl (Chisari et al. 2019a) and its derived harmonic space covariance code TJPCov¹², CosmoLike (Krause & Eifler 2017), or pySSC (Lacasa & Grain 2019), offering comprehensive tools for constructing idealised harmonic space covariance matrices. For real space correlation functions, the excellent cosmocov (Fang et al. 2020) utilises fast logarithmic Fourier transforms to compute real space covariance without the flat sky approximation adopted in this paper, rendering it extremely efficient for this purpose.

However, these tools are either focused on harmonic space or tailored to a specific observable or setup. Consequently, adapting the code to use different summary statistics, observables, or external inputs can be cumbersome. Integrating theory power spectra, trispectra, or power spectrum responses from files, or employing different weighting schemes such as the Bernardeau-Nishimichi-Taruya transformation (Bernardeau et al. 2014) for lensing efficiency, is not straightforward. Similarly, projecting an existing harmonic space covariance to a summary statistic presents challenges.

It is this need for input flexibility and user-friendliness that drove the development of the OneCovariance code. With this objective in mind, the code was designed to offer three key features:

• (i)

  Easy to use: OneCovariance requires standard Python packages and includes its own conda environment to ensure stability. Running the code involves executing a single Python script, while code inputs are specified in a .ini (configuration) file read by the code using the configparser framework. The configuration file's design closely resembles that of class (Lesgourgues 2011; Blas et al. 2011) or CosmoSIS (Zuntz et al. 2015). Sample configuration files, including a comprehensive config.ini file with detailed parameter explanations, are provided.

• (ii)

  Adaptability: The OneCovariance code incorporates a default halo model and HOD-based prescription for biased tracer statistical properties. It communicates with camb (Lewis et al. 2000; Lewis & Bridle 2002; Howlett et al. 2012) for the matter power spectrum, providing both linear and nonlinear corrections. While these ingredients are modular and easily exchangeable, the code accepts almost all critical quantities as input files, enabling flexibility in various scenarios.

  • (a)

    Given an alternative HOD or a mass-concentration relation, various options are available. One can implement them in the hod- and halomodel-class, or encode them and save them into a file to pass to the OneCovariance code. Alternatively, one can calculate the 3-dimensional power spectra, for instance via camb, class, hmcode and baccoemu (Aricò et al. 2021; Angulo et al. 2021) to name a few, to directly provide them to the code. It should be noted that the code itself communicates directly with camb and therefore provides all power spectra implemented there natively.

  • (b)

    It is straightforward to provide harmonic space covariances, i.e., C_ℓ, along with the SSC and NG contribution, for any tracer to the OneCovariance code and project them to an observable.

  • (c)

    Complete freedom is granted in choosing the summary statistics for each tracer. While four hard-coded cases are included: bandpowers, COSEBIs, real space correlation functions, and harmonic space covariance, any summary statistic can be passed as a file to the code, as long as it represents a linear transformation of the C_ℓ. For example, measuring clustering with a real space correlation function, GGL with bandpowers and CS with COSEBIs.

  • (d)

    Consistency checks between different summary statistics are supported. Two summary statistics can be provided to any given tracer, enabling, for instance, the analysis of CS with both COSEBIs or bandpowers while accounting for their cross-covariance.

  • (e)

    Galaxy bias can be determined either by the default HOD prescription or by supplying a file containing the galaxy bias as a function of redshift for each lens bin considered, facilitating the incorporation of numerous (linear) bias models. For non-linear bias, one can pass the resulting power spectrum directly, as discussed in (b).

  This list only scratches the surface of the calculations achievable with the OneCovariance code. A comprehensive array of examples demonstrating the code's functionality can be found on the readthedocs webpage¹³.

• (iii)

  Legacy: In the KiDS collaboration, the steadfast aim has consistently been to utilise diverse sets of summary statistics and analytical methodologies to deduce cosmological parameters reliably, thereby attaining resilient constraints (e.g. Asgari et al. 2021b). Through the provision of a publicly accessible code, we empower the wider scientific community to conduct such analyses, not only with existing KiDS data but also with forthcoming datasets. Furthermore, our endeavour guarantees the applicability of the methodologies honed over years in the KiDS collaboration to future analyses.

The code has been validated against the previously used covariance matrix code (Joachimi et al. 2021) for all statistics (Asgari et al. 2021b) and tracers (Dvornik et al. 2018, 2023), finding per cent agreement. Furthermore, we compared the harmonic space lensing spectra against ccl for further cross-checks and found that they match below one per cent. We show and discuss the results of these tests in Section 7.3.

7. The KiDS-Legacy covariance

In this section, we provide a concise overview of a KiDS-Legacy-like CS sample and outline the error modelling strategies adopted for the associated CS legacy analysis. We elucidate how these choices impact the inference of the structure growth parameter, $S_8≔{\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, as illustrated in Figure 1. Moreover, we underscore the significance of the various terms comprising the covariance matrix, as demonstrated in Figure 2. Finally, we conduct a comparative analysis between the OneCovariance code and selected existing codes, as depicted in Figure 3.

7.1. KiDS-Legacy

The final data release of KiDS (de Jong et al. 2013) and VIKING (VISTA Kilo-degree INfrared Galaxy, Edge et al. 2013) is described in detail in Wright et al. (2024) and is referred to as the fifth data release (dr5). Covering a survey area of 1347 deg², it encompasses 9 photometric bands spanning from optical to near-infrared, with a 5σ limiting magnitude of 24.8 in the r-band. The footprint of the main survey includes a 4deg² overlap with existing deep spectroscopic surveys. This data is complemented by an additional 23 deg² KiDS- and VIKING-like imaging over additional deep spectroscopic fields, yielding a total of about 126 000 sources with both spectroscopy and photometry, crucial for robust redshift calibration (Wright et al. 2025a) across the dr5 footprint. Relative to the fourth data release (Kuijken et al. 2019), dr5 represents a roughly 30 per cent increase in area coverage and around 0.4 magnitude deepening in the i-band.

Significantly, the photometric redshift calibration has been refined, as outlined in Wright et al. (2025a). Notably, a new matching algorithm enables the generation of highly realistic mock catalogues for both photometric and spectroscopic sources. Redshift calibration employs two techniques: colour-based self-organising maps (SOM) and a clustering redshift-based approach, leveraging extensively the spectroscopic samples from dr5. The SOM method yields residual shifts in the mean redshift of each tomographic bin 〈δ_z〉≤0.01, serving as a conservative error floor on the priors for the mean redshifts of each bin. These enhancements, combined with the dr5 data, facilitate the calibration of an additional tomographic bin with redshift up to z≈2, resulting in a total of six bins compared to five in the previous CS analysis in KiDS-1000. The KiDS-Legacy is then a subsample of the whole galaxy catalogue with certain quality requirements, for example for the shape measurements, as well as selection criteria such as masking or blended sources. This leads to a final catalogue of around 4×10⁷ galaxies for CS over an effective area of roughly 1000 deg².

7.2. KiDS-Legacy covariance

In this part, the covariance modelling choices in a KiDS-Legacy-like analysis as it will be carried out in Wright et al. (2025b) and Stölzner et al. (2025) are discussed. While the cosmology remains fixed to the fiducial values presented in Table 1 for all plots, it is important to note that for the actual cosmological analysis, an iterative approach is adopted for covariance calculation. Initially, a fiducial covariance matrix is utilised, followed by maximising of the posterior and iterative updates to the covariance. This iterative process, typically converges after one iteration, as shown in van Uitert et al. (2018), resulting in final contours that exhibit negligible changes, as demonstrated in the subsequent results.

In Figure 1 we show how the inference of S₈ is affected by different modelling choices of the covariance for real space correlation functions and bandpowers (COSEBIs almost pick the same scales as bandpowers, so we do not show them here explicitly). The grey band corresponds to the 68 per cent interval together with the posterior maximum of the marginal distribution for the fiducial choice for the covariance modelling whose settings we are summarising in the following. We model the Gaussian term with an idealised mixed term (see Figure 5), but shape noise is estimated from the catalogue. The covariance always includes a multiplicative shear bias uncertainty, since it is not marginalised over in the inference process via sampling. We include the non-Gaussian covariance and the SSC term with a non-binary mask with the realistic survey footprint to calculate its variance. We use the standard NLA model in the covariance as discussed in Section 4, but stress again that the particular choice does not matter and that the code is flexible enough to deal with any alignment model on the Gaussian covariance level. Lastly, the matter power spectrum is modelled using hmcode2020 (Mead et al. 2020) with feedback controlled via log₁₀T_AGN. Each blue dot with error bars corresponds to an analysis where a specific attribute in the covariance modelling is varied.

Figure 1 shows the influence of individual modelling choices on the marginal constraints on S₈ Nonetheless, it is important to acknowledge altering multiple fiducial assumptions about the covariance, e.g. ignoring the realistic pair counts, changing the IA model and removing the SSC contribution, could lead to a bias in the final value of S₈. However, ignoring the non-Gaussian and SSC contribution on its own has very little effect on the constraints on S₈. This is in particular reassuring for the non-Gaussian term, which is only accurate to 20 per cent (Joachimi et al. 2021). Especially the addition of the sixth tomographic bin in KiDS-Legacy compared to KiDS-1000 will make the nG contribution less pronounced as more linear scales with high signals are added to the analysis. Furthermore, due to the increase in effective area, the survey response is also smaller than in KiDS-1000.

Interestingly, we find that the use of the idealised pair counts (i.e. not using the pair counts from the catalogue but instead the analytic formula) has little effect on the S₈ inference. This might be surprising at first, since using the analytic pair counts will change diagonal elements of the covariance matrix by up to 20 per cent (Troxel et al. 2018). However, we can explain this by the fact that the signal-to-noise ratio of the measurement in this case is almost unchanged. The reason for this is that changing the variance will also change the amount of correlation between the different data points. So if the analytic pair counts over (under) estimate the number of pairs, the variance in the fiducial setting will be larger (smaller), which makes the individual data point less significant for the inference. However, this effect also reduces (increases) correlations between the different data points, thus increasing the information content. Those two effects seem to balance each other such that the S₈ constraints do not change. Interestingly, we also find that using a circular mask instead of the real footprint has the same effect as ignoring the SSC term altogether. IA seems to have a larger effect on real space correlation functions than on bandpowers which could be indicative of the fact that this is driven by contributions from large multipoles to ξ_±. At the same time, however, we find that the case with no feedback changes both real space correlation functions and bandpowers in the same way. In general, we find the consistent trend that contributions decreasing the covariance will bias the constraint on S₈ to larger values, as expected by signal-to-noise considerations. It should also be highlighted that, while the marginal posterior maximum can shift, the upper limit of the 68 per cent interval is very robust. Here we find the largest effect from the new mixed term (compare Figure 5 and Appendix F), since it strongly changes the off-diagonal elements of the covariance matrix while leaving the diagonals untouched. We note that we do not include the new mixed term in the COSEBIs or bandpower covariance matrix at this stage, since they are mapped directly from harmonic space and not from real space to avoid discretisation effects or the computation of a very finely sampled real space covariance matrix.

All those considerations are of course very specific to the KiDS-Legacy survey. They show, however, that our inference of S₈ and in particular the upper limit of the corresponding error is very robust when changing the covariance modelling. While the impact of the mixed term seems very daunting in the light of next-generation surveys it should be highlighted that KiDS is a ground-based, single-visit survey, thus being much more inhomogeneous than Euclid or LSST, which are either in space or scan the same patch of the sky multiple times respectively. It is thus expected that the idealised mixed term is more accurate for these surveys.

In Figure 2, we scrutinise the different contributions of the covariance matrix for the KiDS-Legacy analysis, focusing on the COSEBIs covariance. Since we have shown ξ_± and bandpowers so far, we discuss this using the COSEBIs covariance. Specifically, we plot the Pearson correlation coefficient. The structure of each of the six subplots is the following: the COSEBI mode, n, varies in each small square from n=1,…,5 while the tomographic bin combination i≤j, from one square to the next. We also do not show the covariance of the B_n modes, as they are only given by the shape noise contribution. Lastly, we only include cosmological contributions and do not consider the multiplicative shear bias uncertainty. The six large subsquares of Figure 2 show the six terms of the covariance matrix discussed in Section 2.1. The covariance is dominated completely by the Gaussian term, with the SSC and nG contribution only making up to around 35 and 15 per cent each on the off-diagonals. The Gaussian contribution itself is dominated by the sample variance terms on the off-diagonals, while the shape noise and the mixed term dominate the diagonals, with relative importance increasing towards tomographic bins at higher redshift.

7.3. Comparison with selected existing codes

The analytical covariance utilised in KiDS-1000 has undergone extensive testing against mock data, as detailed in Joachimi et al. (2021). While some of these tests are reiterated in the following sections using updated simulations, we present a selection of results to validate the OneCovariance code.

Given that the OneCovariance code interfaces with camb for linear and non-linear matter power spectra, and other codes like ccl employ a modified version of our halo model trispectrum implementation, we do not validate these quantities directly. The first step of our validation is therefore to compare the CS angular power spectra, or in other words, the line-of-sight projection. On the left side of Figure 3 we demonstrate the relative accuracy of the OneCovariance code compared to ccl. A KiDS-Legacy-like setup with six source bins was used and the non-linear matter power spectrum has been calculated using halofit (Takahashi et al. 2012). The tomographic bin combination index is colour-coded, that is labelling each bin combination, i≤j with a unique number starting at one. We find that the agreement is well below one per cent, except for the very first tomographic bin. This difference originates from slightly different extrapolations used for the non-linear matter power spectrum (as can be seen by the steep rise at large multiples) as well as the kind of interpolation used for the source redshift distributions. The small visible spikes are of order 10⁻³ and can be traced back to different accuracy settings when CAMB is called via ccl or the OneCovariance as well as slightly different interpolation schemes. However, we find excellent agreement with ccl.

Another critical validation step involves comparing our code with the covariance code used in the KiDS-1000 analysis. The right panel of Figure 3 displays the percentage difference between the OneCovariance code and the previously employed code for the diagonal elements of the ξ_± covariance matrix, excluding shape noise. We conducted this comparison across nine logarithmically spaced angular bins ranging from 0.5 to 300 arcmin. Once again, we observe excellent agreement, with differences within a few per cent. Slightly larger disparities are noted at small angular scales, primarily attributed to differences in integration accuracy and cut-off limits for the integrals, particularly for the non-Gaussian (nG) and super-sample covariance (SSC) contributions. However, these discrepancies, while reaching up to 15 per cent, are not significant as small scales are shot-noise dominated, thereby minimally impacting previous analyses. Similar validation procedures were repeated for bandpowers and COSEBIs, yielding agreement below a per cent with previous studies (Joachimi et al. 2021; Asgari et al. 2021b). Notably, due to the more compact Fourier filters, the agreement for these cases is even better than for real space statistics. Additionally, we verified that we obtain identical results when passing Fourier filter functions externally to the code for different statistics (for how to pass external information, see Appendix A).

8. Validation against simulations

Testing the CS signal and noise modelling for the final release of the KiDS data, dr5, necessitates a substantial number of realisations of the data to accurately characterise the covariance. Achieving this requires forward simulations that closely resemble real data, incorporating factors such as the survey mask and variable depth. In this section, we detail the construction of the forward simulations employed to validate the error modelling in KiDS-Legacy. Section 8.1 provides an overview of the log-normal mocks and Section 8.2 briefly describes how the mocks are populated with galaxies and how the shape catalogues are created. Lastly, in Section 8.3 we compare the simulations against the analytical covariance.

Central to our assessment are key comparison figures that directly juxtapose our most realistic mock with the OneCovariance code. As all summary statistics can be derived from the real space correlation function, we only measure those on the mocks. Figures 4, 5, and 6 show the comparison for the signal, the correlation coefficient and diagonal covariance matrix elements respectively.

8.1. GLASS simulations

The forward simulations for covariance validation in this study are founded on the KiDS-SBI¹⁴ suite of simulations detailed in von Wietersheim-Kramsta et al. (2025). At the core of the pipeline lie log-normal mocks of the 3-dimensional matter distribution in radial shells, which are subsequently projected along the line of sight to create shear maps. This process is facilitated within the GLASS¹⁵ (Generator for Large-Scale Structure) framework (Tessore et al. 2023). GLASS generates log-normal realisations to match the two-point statistics of the given input power spectrum. While higher-order statistics are not precisely recovered due to the inherently non-log-normal nature of the cosmic density field, this nonlinear transformation captures a significant portion of the trispectrum (Friedrich et al. 2021; Hall & Taylor 2022) relevant for the covariance matrix.

Within this framework, 4224 simulations are created for two distinct survey configurations, akin to the one selected in Joachimi et al. (2021), a discussion of which is provided in further detail in Section 8.2. The cosmological parameters remain fixed at the values specified in Table 1. The 3-dimensional power spectrum of the input matter is computed using camb (Lewis et al. 2000; Lewis & Bridle 2002; Howlett et al. 2012), with non-linear corrections computed from hmcode2020 (Mead et al. 2020). After computing the angular power spectra of the matter distribution in 22 shells of roughly equal thickness spanning the KiDS-Legacy redshift range, GLASS generates corresponding log-normal matter fields in each shell. The shell thickness ranges between 100 and 200 h⁻¹Mpc, ensuring a reasonably accurate description of the matter field as log-normal (Hall & Taylor 2022; Piras et al. 2023), while also being sufficiently finely sampled to accurately capture the matter overdensity along the line of sight. The computationally intensive integration for the C_ℓ in the matter shells is executed using Levin's method (Levin 1996; Zieser & Merkel 2016; Leonard et al. 2023). Subsequently, the concentric shells are integrated along the line of sight, weighted by the lensing efficiency function (refer to Tessore et al. 2023, for detailed information), resulting in a CS field realisation covering the entire sky. It is noteworthy that intrinsic alignments are not included in the covariance matrix validation. However, this omission does not compromise the validity of our tests, given that intrinsic alignment exerts minimal influence on the covariance matrix (see Figure 1) even if removing it completely and generally shares the same functional form as the lensing signal, thus being equally well accounted for.

8.2. KiDS-Legacy-like mocks

Given a realisation of the density and shear fields, galaxies are sampled within those, including systematic effects such as the survey footprint and the effect of variable depth on the redshift distribution and the shape noise with the help of SALMO (Speedy Acquisition for Lensing and Matter Observables, Joachimi et al. 2021).

To create a realistic number density of galaxies, we utilise organised randoms (ORs, introduced in Johnston et al. 2021). ORs employ self-organising-maps (SOMs, Kohonen 1982) to generate random samples that follow the spatial variability induced by a high-dimensional systematic space. SOMs are unsupervised neural network algorithms which are designed to map a high-dimensional data space onto a two-dimensional map while preserving topological features of the input data space. We define a systematic data vector including atmospheric seeing, sky background light, stellar number density, dust extinction and variations in the point-spread function. The SOM provides a map of Tiaogeng (TG, see Yan et al. 2025, for a detailed explanation) weights, w_TG, as a function of pixel on the sky. The TG weights quantify the anisotropic selection function on the galaxy sample observed within KiDS as a function of the aforementioned observational systematics in a manner that is independent of the cosmological signal in the data. We found the TG weights to be a good indicator of the local seeing, atmospheric transparency and signal-to-noise, so it can be considered a good estimator of variable depth. In particular, the TG weights are uncorrelated with the r-band magnitude measurements of galaxies and their photometric redshift, while being highly correlated with the magnitude limit in the r-band which KiDS uses for shape measurements, the degree of background noise, and the PSF variations. Hence, w_TG can robustly predict the local observational depth independently of the galaxy position and shape measurements. Next, we partition each of the six tomographic bins into ten equi-populated bins in w_TG. For each ${\mathit{w}}_{\mathrm{TG}}^i$ we independently recompute the effective number density, n_eff, and the total ellipticity dispersion σ_ϵ, from the galaxy population. It is possible to linearly interpolate n_eff and σ_ϵ as a function of w_TG allowing the evaluations of both quantities at every pixel in the KiDS-Legacy-like footprint from the single spatial w_TG map. We use a footprint of approximately 1000 deg² with an N_side=1024.

The redshift calibration itself uses the following photometric bin boundaries: [0.1, 0.42, 0.58, 0.71, 0.9, 1.14, 2), adding a sixth tomographic bin in KiDS-Legacy compared to KiDS-1000 (see Section 7.1) By applying this binning scheme to the mock catalogues, we obtain the redshift distributions discussed above from the SOM, mapping the observed photometric redshifts to spectroscopic redshifts based on a reference spectroscopic sample (Wright et al. 2020; Hildebrandt et al. 2021).

Having set up the forward simulations, the two-point correlation functions ξ_± are measured on each realisation using TreeCorr (Jarvis et al. 2004), following the same implementation described in Giblin et al. (2021), Joachimi et al. (2021). Initially, ξ_± are measured in the north and south separately using 300 logarithmically spaced bins in θ∈[0.5, 300] arcmin. The estimated correlation functions, ${\hat{\xi }}_{\pm }$ are then combined as a weighted mean,

$${\hat{\xi }}_{\pm }^{(\mathit{ij})}(\theta )=\frac{N_{\mathrm{pairs},\mathrm{N}}^{(\mathit{ij})}{\hat{\xi }}_{\pm ,\mathrm{N}}^{(\mathit{ij})}(\theta )+N_{\mathrm{pairs},\mathrm{S}}^{(\mathit{ij})}{\hat{\xi }}_{\pm ,\mathrm{S}}^{(\mathit{ij})}(\theta )}{N_{\mathrm{pairs},\mathrm{N}}^{(\mathit{ij})}+N_{\mathrm{pairs},\mathrm{S}}^{(\mathit{ij})}},$$(89)

with the weighted galaxy pairs and measured correlation functions in tomographic bins i and j in the north (N) and south (S) fields respectively. Lastly, the finely binned correlation functions are then re-binned into the final nine logarithmically spaced θ bins covering the same range.

In Joachimi et al. (2021), three kinds of mocks, Buceros, Cygnus and Egretta, with increasing complexity were distinguished. Buceros is an idealised survey with a simple rectangular footprint and a homogeneous number density which serves as a baseline check of the covariance modelling, which was already found to agree within a few per cent with the analytical covariance (Joachimi et al. 2021). Given the agreement of the OneCovariance with the KiDS-1000 covariance matrix, see Figure 3, this mock as such is not considered. In case of Cygnus, we mask the simulations according to a realistic KiDS-Legacy-like mask. However, the survey is kept homogenous in n_eff and σ_ϵ. In particular, the intrinsic ellipticities are sampled from a Gaussian distribution with fixed σ_ϵ and zero mean for each tomographic bin. In the most complex mock, Egretta, the homogeneity and isotropy assumption in n(z), n_eff and σ_ϵ are dropped, and they are sampled from realisations in the w_TG maps. Therefore, Egretta is a realistic representation of the real KiDS-Legacy-like catalogue, including all systematics discussed before. It will therefore be used as the benchmark for the analytical covariance matrix. For more details on the simulations, we refer the reader to Joachimi et al. (2021).

8.3. Comparison with mocks

We limit the discussion to the most realistic mock, Egretta, since a general comparison between the different mock settings was already done in Joachimi et al. (2021) which the OneCovariance code was validated against. Having said that, we want to briefly summarise those results here again. For the Buceros mock with an idealised footprint and homogeneous survey properties, we find agreement and the percentage with slightly larger difference at large separations. Due to the more complex survey mask in Cygnus, those difference become more pronounced, leading to ten per cent differences at 300 arcmin separation. At smaller scales, however, the difference is still well below five per cent.

Before delving into the covariance comparison with Egretta we show the relative difference between the signal measured in the simulations for ξ_± and calculated from theory via the OneCovariance code is shown in Figure 4. Dashed lines represent ξ₋ and solid lines ξ₊ while the colour bar indicates the tomographic bin index combination. The four panels show different settings of the theory calculation. In particular, the resolution of the simulation is partially mimicked by taking the pixel window into account which dampens the power on scales, ℓ≤N_side and we explicitly average over angular bins in which the correlation function is measured. The grey band shows a 5 per cent difference between the simulation and the theory predictions. In general, we see good agreement between the theoretical predictions and the simulation, within 5 per cent. Incorporating the damping of small-scale power into the theoretical predictions slightly enhances the agreement on small angular scales. However, because of the broad filter function in Fourier space of ξ_±, the small angular bins necessitate information from relatively small scales and are thus not accurately represented in the simulation. This effect is particularly noticeable for ξ₋ which requires integration up to ℓ>10⁴ to reach convergence for θ<arcmin, although it is not evident for ξ₊ at the scales displayed in this figure for clarity. For ξ₊, this effect can be observed at scales below a few arcminutes. Nevertheless, for the covariance comparison, this limitation should not present a significant issue, as the covariance will be predominantly influenced by shape noise on small scales. Additionally, we observe that bin-averaging introduces differences of a few per cent across all angular bins, with no significant variation as a function of θ, consistent with expectations due to the logarithmic binning scheme.

On the left side of Figure 5, we present the covariance for the setup outlined in Section 8.2. The lower triangle displays the analytic prediction, while the upper triangle shows the covariance estimated from the mock data. Visually, there is excellent agreement between the simulation and theory for the correlation coefficient. The elements are ordered in such a way that each square corresponds to a unique combination of tomographic bins, arranged in increasing order such that i≤j. Overall, it is evident that the covariance is dominated by shape noise. However, for the tomographic bins at higher redshifts, as discussed in more detail in Section 7.2, other contributions become more pronounced.

The right panel of Figure 5 illustrates the impact of this modelling choice for the mixed term on the covariance matrix. Specifically, it shows the relative difference between these two cases for the Gaussian covariance matrix, i.e. considering an idealised mixed-term with homogeneous number density and no masking effects relative to using the triplet counts directly from the catalogue (compare Equation (6) and Appendix F respectively). Generally, we observe that the effect is most pronounced at large or small angular separations, while intermediate separations are less affected. This behaviour is expected because large scales are influenced by the extent of the survey, while very small scales can be affected by masking effects within the survey area. Additionally, the effect on the diagonal is less noticeable due to the contribution of shape noise.

For ξ₊, the effect of the updated mixed term reaches a maximum of 10 per cent. However, for ξ₋ and its cross-correlation with ξ₊, the effect can be significant for individual elements of the covariance matrix, exceeding 50 per cent. Nonetheless, since the signal strength for ξ₋ is relatively small, this effect is mitigated by the shape noise (its impact on inference is discussed in Section 7.2). It is worth noting that although there are individual points showing a substantial effect in the cross-covariance between ξ₊ and ξ₋, this is attributed to numerical noise and does not influence the results since the corresponding covariance matrix elements are very small.

We anticipate this effect to be less pronounced for larger and more homogeneous surveys than KiDS-Legacy. Surveys like Euclid and LSST are expected to have fewer spatial variations in the number density of sources. This statement, however, bears the caveat that the other components of the covariance matrix will become more significant for these surveys as they are deeper, resulting in a stronger CS signal compared to the shape noise. Therefore, a concluding analysis of the impact of the mixed term in a stage 4 survey will be the subject of future work.

Finally, Figure 6 illustrates the relative difference in the standard deviation, i.e. the square root of the diagonal elements of the covariance matrix, between Egretta and the OneCovariance code. There is excellent agreement at small angular scales due to the contribution of shape noise, as the pair counts are matched to those measured in the Egretta mocks. It is worth noting that we address the most significant effect of variable depth, which generally changes the shape-noise in KiDS-Legacy compared to a uniform galaxy distribution. This effect is driven by the distribution of w_TG, which skews above the mean, resulting in most pixels in the field being under-dense compared to the mean galaxy density. Generally, the agreement between the mocks and the analytical covariance is of a slightly better level than in previous analysis (Joachimi et al. 2021). Due to the increased sensitivity of KiDS-Legacy compared to KiDS-1000 the question arises whether this agreement is good enough. As discussed above, the main effect is the survey geometry which we address in Appendix F reconsideration of the mix term and via the pair counts in Equation (64). In the ‘sva’ term, however, the survey geometry is only accounted for via the f_sky approximation. To asses the accuracy we analysed mock data vectors with the fiducial KiDS-Legacy covariance (compare to the grey band in Figure 1), including the new mixed term, and directly with Egretta. We found that using the Egretta covariance gives a shift in χ² similar to the inclusion of the survey geometry in the mix term. From this we conclude that the f_sky approximation in the ‘sva’ term is still accurate enough for KiDS-Legacy. If the sensitivity with upcoming surveys is further increased, this approximation might need to be reconsidered. This point is discussed a bit more in the conclusions.

9. OneCovariance applications

To demonstrate the versatility of the OneCovariance code, we employ it to address two additional key aspects of the KiDS-Legacy analysis: consistency tests among different summary statistics in Section 9.1 and clustering redshifts in Section 9.2. These aspects are integral to the studies presented in Stölzner et al. (2025) and Wright et al. (2025b), respectively. In Figure 7, we present the correlation coefficient between the three main summary statistics utilised in KiDS-Legacy. For the clustering redshift covariance, our primary findings are depicted in Figures 8 and 10. These illustrations offer a comparison with mocks and validate the fundamental structure of the jackknife covariance employed for the redshift calibration.

9.1. Consistency tests

As an initial showcase of the OneCovariance code's capabilities for KiDS-Legacy CS, we present the covariance matrix between various summary statistics. This calculation, facilitated by the ARBsummary class of the code (see Appendix A), serves as a foundational component for the consistency checks (or cross-validation) that will be carried out in Stölzner et al. (2025), a process that hinges not only on calculating the covariance of each summary statistic but also on their cross-covariance. Given that the different summary statistics probe nearly identical physical scales, the resulting covariance matrix exhibits significant cross-correlation and can be close to singular.

As a specific example and to demonstrate the versatility of the OneCovariance code, we undertake the covariance matrix calculation for a KiDS-1000 setup. This pairwise computation encompasses our principal CS summary statistics. It is, however, available for arbitrary summary statistics of all tracers supported by the code.

The resulting correlation coefficient is depicted in Figure 7, showing large off-diagonal elements between different summary statistics originating both from sample variance terms as well as shape noise. Moreover, the pronounced correlations highlight that the diverse summary statistics employed in the KiDS-Legacy analysis are almost completely degenerate.

The OneCovariance code can straightforwardly accommodate any combination of summary statistics envisioned by the user since it propagates the relation between summary statistics and the harmonic covariance through simple integral transformations (compare Equation (51)). Consequently, users can seamlessly integrate different summary statistics for various tracers or explore the impact of different scale cuts for any summary statistic, as done for COSEBIs in Dark Energy Survey and Kilo-Degree Survey Collaboration (2023).

9.2. Clustering redshifts

As the final application of the OneCovariance code in this study, we investigate the covariance matrix of clustering redshift estimation (see e.g. Lima et al. 2008; Bonnett et al. 2016; van den Busch et al. 2020; Hildebrandt et al. 2021) which utilises a spectroscopic reference sample (s) and a target sample with noisy redshift estimates (p). In Hildebrandt et al. (2021) a covariance using jackknife resampling was used. So far, however, it was never compared to analytical expectation. In this section we are filling this gap. By measuring the correlation function, Equation (59), between these two samples at fixed physical scale r at different redshifts, the redshift distribution n(z) can be recovered up to an irrelevant multiplicative constant (as the redshift distribution will be normalised in the end) and bias evolution of the target sample,

$$n_{\mathrm{p}}(z)=\frac{{\mathit{w}}_{\mathrm{sp}}(\theta (r),z)}{\mathrm{\Delta }z\sqrt{{\mathit{w}}_{\mathrm{ss}}(\theta (r),z)}},$$(90)

where Δz is the redshift width of each cross-correlation measurement. By labelling tomographic bin indices of the target sample with Greek letters and reference sample redshift bins with Latin letters, the covariance for this estimator at first order is

$$\begin{array}{cc}\hfill \mathrm{Cov}\left[n_{\alpha }(z_i),n_{\beta }(z_j)\right]& \equiv \mathrm{Cov}\left[\frac{{\mathit{w}}_{i\alpha }}{\mathrm{\Delta }{z}_i\sqrt{{\mathit{w}}_{\mathit{ii}}}}\frac{{\mathit{w}}_{j\beta }}{\mathrm{\Delta }{z}_j\sqrt{{\mathit{w}}_{\mathit{jj}}}}\right]\hfill \\ \hfill & \approx \frac{1}{\mathrm{\Delta }{z}_i\mathrm{\Delta }{z}_j\sqrt{{\mathit{w}}_{\mathit{ii}}{\mathit{w}}_{\mathit{jj}}}}\left[C_{i\alpha j\beta }-\frac{{\mathit{w}}_{i\alpha }}{2{\mathit{w}}_{\mathit{ii}}}{C}_{\mathit{iij}\beta }-\frac{{\mathit{w}}_{j\beta }}{2{\mathit{w}}_{\mathit{jj}}}{C}_{i\alpha \mathit{jj}}+\frac{{\mathit{w}}_{i\alpha }{\mathit{w}}_{j\beta }}{4{\mathit{w}}_{\mathit{ii}}{\mathit{w}}_{\mathit{jj}}}{C}_{\mathit{iijj}}\right].\hfill \end{array}$$(91)

Here C_ijkm denotes the covariance between clustering measurement w_ij and w_km at fixed scale r. The last three terms are only non-zero if i=j, with some small off-diagonal elements from the full non-Limber expression. Only the first term can have important off-diagonal elements, in particular where the redshift of the reference sample i approaches the support of the target sample β (and the same for j and α).

The performance evaluation of the analytic covariance for clustering redshifts is conducted using the synthetic data outlined in Hildebrandt et al. (2021) and van den Busch et al. (2020). This dataset comprises mock galaxy samples closely resembling the KiDS+VIKING-450 (KV450) CS sample. The basis of these mock samples is the MICE simulation (Fosalba et al. 2015). The cosmological model adopted in these simulations adheres to a ΛCDM model with parameters: Ω_m=0.25, Ω_Λ=1−Ω_m, Ω_b=0.044, σ₈=0.8, and h=0.7. These simulations are accompanied by the corresponding galaxy catalogues (Crocce et al. 2015), serving as the foundation for the publicly accessible pipeline for generating KiDS mock samples¹⁶. Due to the nature of the estimator, the covariance retains a residual dependence on the bias evolution of the target sample. Davis et al. (2018) characterised this bias using the parametrisation ℬ_α(z)=(1+z)^α, we use the values from van den Busch et al. (2020), Hildebrandt et al. (2021) for the residual bias evaluation in the covariance modelling.

Figure 8 illustrates the resulting redshift distributions for the five tomographic bins of KiDS1000, represented by solid colour-coded lines on the left side, denoting the target sample. The redshift distribution of the reference sample is depicted as a bar chart. On the right, the Pearson correlation coefficient of the corresponding covariance matrix is displayed for scales r∈[0.1,1] h⁻¹ Mpc. In the upper triangle, we show the analytic prediction of Equation (91), while the lower triangle uses the mock catalogue and jackknife re-sampling over the individual subsamples of the simulation. The general structure is very similar, in particular, one can see that the previously mentioned non-vanishing off-diagonals at i≠j are in the correct position. Moreover, the diagonal of each block matrix exhibits comparable structures, providing valuable qualitative validation of the jackknife covariance outlined in Hildebrandt et al. (2021).

With the narrow redshift bins employed in cross-correlation measurements with spectroscopic surveys, the applicability range of the Limber approximation, as denoted by Equation (10), is significantly constrained. As discussed in Loverde & Afshordi (2008), the next-to-leading order term in the expansion becomes subdominant for ℓ>χ/Δχ, where Δχ is the width of a Gaussian redshift bin at co-moving distance χ. For the present situation, this requires the non-Limber calculation up to ℓ∼10³. In the left plot of Figure 9, we depict the angular power spectrum of spectroscopic sample auto-correlations. Dashed lines represent the complete integral, as described in Equation (8), while solid lines utilise the Limber approximation. All lines are colour-coded with respect to the spectroscopic redshift bins (as illustrated in the bar plot of Figure 8). It is evident from the figure that the non-Limber projection significantly impacts the integrand of ${\mathit{w}}_{{\mathrm{s}}_i{\mathrm{s}}_j}$. On the right side of Figure 9, we illustrate the impact of the non-Limber approximation on the Gaussian error of the covariance, revealing approximately 50 per cent effects for all tomographic bins.

Figure 10 presents, on the left side, the contributions to the covariance as a function of redshift, colour-coded according to the tomographic bin index. It is evident that low redshifts are significantly influenced by the non-Gaussian covariance (as indicated by dashed lines), while the shot noise contribution remains negligible for all redshifts. On the right side of Figure 10, we juxtapose the analytical covariance with simulations (indicated by symbols with thin lines). We observe a good agreement for the lower tomographic bins, with the increase attributed to the non-Gaussian contribution clearly visible. However, this increase is less pronounced in the fourth and fifth tomographic bins for the analytical covariance, whereas it remains evident in the mock data. This difference is likely due to uncertainties in the non-linear modelling of the trispectrum, especially as we measure signals deep within the non-linear regime, k>1h/Mpc.

It is worth noting that while the non-Limber effect is easily detectable for the Gaussian covariance, its influence on the non-Gaussian term presents a more challenging task and has thus far only been explored in Lee & Dvorkin (2020) for trispectra using specific kernels from perturbation theory and the galaxy bias expansion. Their analysis suggests significant effects on the trispectrum. However, a similar analysis for general trispectra, such as those from the halo model, remains unexplored and is beyond the scope of this section for several reasons. Firstly, while the OneCovariance can compute the non-Limber projection of power spectra, its primary utility lies in photometric surveys with broader window functions, where Limber's approximation is more readily applicable. Secondly, spectroscopic samples, necessitating non-Limber calculations, are typically conducted on quasi-linear scales, unlike the one presented here. On those scales, the non-Gaussian contribution is generally less significant. Lastly, this application serves as a comprehensive test and sanity check for clustering redshift covariance, where jackknife resampling is readily available. Furthermore, realistic samples are likely to be more influenced by noise, as galaxy density tends to be substantially lower, at least by a factor of ten.

10. Present and future of the OneCovariance code

In this section, we offer a summary of the capabilities of the code, delineate the modelling choices incorporated, underscore its flexibility, and contemplate potential avenues for future enhancements.

10.1. Current status of OneCovariance code

We have provided a comprehensive review of the covariance matrix modelling in KiDS so far in Sections 3 and 4. Employing a halo model approach, we have consolidated all relevant equations used in various analyses, encompassing CS (Hildebrandt et al. 2020; Asgari et al. 2021b), GGL, GC, and the CSMF (Heymans et al. 2021; Dvornik et al. 2018, 2023). This effort culminates in establishing the legacy of analytic covariance modelling within KiDS. While our primary focus lies on CS for the direct validation of the code in this study, we have implemented all pertinent quantities into the new OneCovariance code from scratch and cross-validated them against previous results. Moreover, we have reviewed and incorporated all summary statistics utilised within KiDS in Section 5, encompassing real space correlation functions, COSEBIs, and bandpowers, each applicable across all three tracers. Through the development of a unified code for covariance matrices, we have streamlined the calculation process, now directly linking harmonic space statistics to observables via integral transformations. This approach differs from that of KiDS-1000, where the bandpowers covariance was derived from the ξ_± covariance. This modification brings the analytic covariance for bandpowers closer to the theoretical modules employed in KiDS, enhancing the consistency and efficiency of our analyses.

The OneCovariance offers users the capability to effortlessly download and compute simulation-validated covariance matrices for different traces across various statistical representations including harmonic space, real space, COSEBIs, and bandpowers, all under the flat sky approximation. While we focused on the comparison using CS in this paper, the code has been validated against all previous covariance codes employed in KiDS also for the other tracers. Minimal user input is required beyond specifying the redshift distributions, as the code operates out of the box. Moreover, the code's flexibility allows for extensive customisation. Users can readily adjust inputs, such as modifying the HOD prescription or directly altering galaxy power spectra or bias parameters. Additionally, it offers the flexibility to alter line-of-sight weight functions through input, facilitating the definition of entirely new tracers. For instance, by defining a new line-of-sight weight function and adjusting the HOD prescription, the code can compute the covariance matrix for the Cosmic Infrared Background or any other intensity mapping quantity. Alternatively, users can provide harmonic-space quantities directly, allowing the code to perform the necessary projections to observables. Furthermore, the code includes scripts to compute the weight functions used for various summary statistics in KiDS-Legacy.

The overarching philosophy of the OneCovariance code is centred on the replaceability of all essential quantities by file inputs. If a file input is provided, it supersedes the standard configuration of the code, utilising the file input instead. This approach permeates the entirety of the code, from fundamental quantities required for theory calculations, such as line-of-sight weights and bias prescriptions, to higher-level quantities including power spectra, their responses, and trispectra. The final elements of the covariance are effectively arrays for each tracer combination with the general shape

$$({\mathrm{spatial}}_{t_1,t_2},\hspace{0.5em}{\mathrm{spatial}}_{t_3,t_4},\hspace{0.5em}{\mathrm{mass}}_{t_1,t_2}\hspace{0.5em},{\mathrm{mass}}_{t_3,t_4}\hspace{0.5em},{\mathrm{tomo}}_{i,t_1}\hspace{0.5em},{\mathrm{tomo}}_{j,t_2}\hspace{0.5em},{\mathrm{tomo}}_{i,t_3}\hspace{0.5em},{\mathrm{tomo}}_{j,t_4})\hspace{0.5em}.$$

Here, ${\mathrm{spatial}}_{t_1,t_2}$ labels the spatial index of the tracer pair t₁,t₂ over which its two-point function is measured. This can be for example the θ bin, the band power mode, the order of the COSEBIs, or any other label of an arbitrary summary statistic passed to the code. ${\mathrm{mass}}_{t_1,t_2}$ is the corresponding stellar mass bin (this could be also modified to cluster number counts, for example). For CS, this always has a length of one. Lastly, ${\mathrm{tomo}}_{i,t_1}$ and ${\mathrm{tomo}}_{j,t_2}$ is the tomographic bin combination for the pair (t₁,t₂). With this structured approach, the OneCovariance code possesses the versatility required to handle various scenarios in cosmology, allowing users to specify inputs via file paths in a configuration file. This means that users can specify the necessary inputs, thereby solving the often cumbersome task of ordering covariance matrix elements and performing the projection to observables. In essence, the OneCovariance code provides a comprehensive solution for covariance estimation whenever a HOD prescription-based halo model is employed for the signal, offering a valuable resource for the cosmological community.

10.2. Future developments of the code

The future development roadmap for the OneCovariance code includes both incremental improvements and significant expansions to its functionality. Here is an outline of the planned enhancements:

• i)

  Currently, the code employs the flat sky approximation for defining observables. In the next version, full sky transformations will be implemented, especially for the Gaussian terms, as they dominate the covariance on scales where curved sky effects become significant.

• ii)

  The code already supports calculating the non-Limber covariance for Gaussian contributions. The next step involves extending this capability to cover connected terms as well, enhancing the accuracy of covariance estimates.

• iii)

  Incorporating survey window effects into the harmonic space covariance and propagating them to observables, particularly for sample variance, will be a crucial addition to enhancing the accuracy of covariance estimates. This could be done by streamlining the code with NaMaster (Alonso et al. 2019).

• iv)

  A longer list of pre-implemented models will be added to reduce the number of required input quantities, including alternative HOD and bias prescriptions, streamlining the user experience and expanding the code's usability.

• v)

  Interfacing CCL (Chisari et al. 2019a) will be implemented, enabling the use of pre-defined quantities within CCL and leveraging its extensive functionalities for enhanced accuracy and efficiency in covariance calculations.

11. Conclusion

In this study, we have provided a comprehensive account of the covariance methodology utilised in and up to the latest data release of KiDS. This encompasses the validation of the CS covariance matrix utilised in Wright et al. (2025b) and Stölzner et al. (2025), as well as its application to clustering redshifts as employed in Wright et al. (2025a). Accompanying this paper is the OneCovariance code, a unified covariance tool tailored for photometric LSS surveys. We have conducted a thorough validation of the OneCovariance code for analysing KiDS-Legacy data, with a specific focus on CS. Our validation efforts have yielded several noteworthy conclusions. We have demonstrated that the theoretical predictions generated by the OneCovariance code align closely with those obtained from external libraries, such as ccl, at the per cent level. This underscores the robustness and accuracy of the code's theoretical predictions. Furthermore, comparisons with the independent implementation of the covariance matrix calculations from KiDS-1000, as detailed in Joachimi et al. (2021), being typically at the per cent level, reinforces the accuracy of the OneCovariance code. Any larger discrepancies identified were traced back to differences in accuracy settings or extrapolation schemes used for integrations when transitioning from harmonic space to observables. Importantly, we have verified that these discrepancies would not have significantly impacted any previous analyses conducted on the KiDS data. This ensures the consistency and reliability of previous findings derived from KiDS data analyses.

Our validation efforts provide strong evidence supporting the efficacy and accuracy of the OneCovariance code for CS analyses within the KiDS framework. This instils confidence in the code's utility for future cosmological investigations, ensuring robust and reliable covariance estimation in the analysis of LSS data. Next, we established a simulation pipeline to further validate the covariance matrix. We used log-normal realisations of the density field with both a realistic mask and homogeneous sources (Cygnus), as well as variations in source density (Egretta). On these mocks, we measured the real space correlation functions of the shear, ξ_±, and found excellent agreement between the theoretical signal and the signal measured in the mocks. This agreement was observed across nine logarithmically spaced angular bins ranging from 0.5 to 300 arcmin. By leveraging the exact pair counts from the Egretta mocks, we also found similar agreement between the simulation and the analytic covariance as in KiDS-1000. Particularly, we observed agreement well within ten per cent on most scales, with only the largest angular scales exhibiting more significant deviations due to survey boundary effects. Importantly, these differences do not significantly impact the final constraints derived from the analysis.

A significant deviation from KiDS-1000 was the re-evaluation of the mixed term, namely the cross-correlation between shape noise and two-point statistics. This term is influenced by both the response of the signal to the survey geometry and the number of pairs used to estimate the signal itself, owing to the complex survey mask of KiDS. We found that assuming a homogeneous survey misestimates the Gaussian term by roughly 10 per cent in the case of ξ₊ for off-diagonal elements, with this effect potentially being both positive and negative. For ξ₋ and its cross-correlation with ξ₊, the deviation can be substantial on off-diagonal elements. This is attributed to the broad support of the Fourier filter for ξ₋ and its sensitivity to many scales, thus being affected by the inhomogeneity of the KiDS-Legacy sample and the intricate survey footprint. However, despite the potentially large effects on individual covariance matrix elements, the low signal of ξ₋ and high shape noise mitigate these discrepancies. For upcoming surveys, the situation is less clear. On the one hand, shape noise effects will become less important and the surveys will become more homogeneous as well as larger, reducing the importance of edge effects and variable depth. On the other hand, however, due to the increasing importance of sample variance and the increased overall sensitivity, the usage of an ideal mixed term might still lead to more significant effects on parameter inference. In future work, we plan an investigation of this issue along with the inclusion of the survey geometry in the ‘sva’ term. This can be done by using the pseudo-C_ℓ covariance directly to incorporate the effect of the mask.

We examined the influence of various modelling choices for the covariance on the final inference of S₈, demonstrating their robustness. These choices included the baryonic feedback prescription, the Super Sample Covariance (SSC) and non-Gaussian (nG) terms, the updated mixed term, and the intrinsic alignment model, among others. The upper bound of the 68 per cent interval of S₈ remained largely unaffected by these modelling choices. Furthermore, we compared the impact of different terms on the covariance matrix and found that the Gaussian term overwhelmingly dominated on almost all scales and tomographic bins. The non-Gaussian and SSC contributions were of minor importance, with the SSC term slightly more significant across tomographic bins. One contributing factor to the diminished impact of these terms in KiDS-Legacy is the inclusion of the sixth tomographic bin, which enhances sensitivity to more linear scales, especially given the substantial CS signal in the sixth bin.

Furthermore, we demonstrated that the OneCovariance code can serve as a valuable tool for conducting consistency tests between different summary statistics or different surveys. That is, we calculate the covariance matrix for the same tracer (CS in this case), but two different summary statistics. We confirmed that the covariance matrices required for the KiDS-Legacy consistency tests (Stölzner et al. 2025) are invertible and well-behaved for pairs of summary statistics. This functionality holds significant importance when performing consistency tests between large-scale and small-scale analyses of CS data, particularly in addressing tensions such as the S₈ tension. By enabling such tests, the OneCovariance code enhances the capability to probe and validate cosmological models across different scales and datasets.

In our final application of the OneCovariance code, we turn our attention to the estimation of the covariance matrix for clustering redshifts, a critical component utilised in Wright et al. (2025b) to validate the jackknife covariance employed for n(z) calibration. Our approach involved applying the code to a KiDS-1000 redshift distribution, slated for calibration against an idealised reference sample derived from the MICE2 simulations, approximately ten times denser than the actual dataset. Our analysis revealed encouraging results; there was a generally good agreement observed for the correlation coefficient between the analytical and jackknife covariance derived directly from the mocks. However, the narrow redshift shells over which the signal is averaged at a given physical scale posed a challenge, rendering the Limber expression inapplicable over a broad range of scales. Consequently, the covariance experienced deviations of up to 50 per cent. To address this issue, we incorporated the full non-Limber expression in the Gaussian covariance, while assuming the Limber expression for the non-Gaussian component. However, we recognise that a comprehensive treatment necessitates a full non-Limber modelling of the non-Gaussian covariance, a task beyond the scope of this paper. Nevertheless, our analytical covariance matrix effectively captured the salient features of the jackknife covariance, thereby serving as a valuable cross-check for the clustering redshift calibration process.

In summary, our study confirms the robustness of the KiDS-Legacy covariance in delivering accurate results for the S₈ parameter. By providing the cosmological community with a versatile tool capable of calculating covariance matrices for a wide array of photometric(-like) LSS observables in both current and future surveys, we hope to facilitate further advancements in cosmological research.

Acknowledgments

The authors would like to thank an anonymous referee for valuable comments. RR, AD, JLvdB, HH and CM are supported by a European Research Council Consolidator Grant (No. 770935). SU, BS and ZY acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. HH is supported by a DFG Heisenberg grant (Hi 1495/5-1), the DFG Collaborative Research Center SFB1491 and the DLR project 50QE2305. BJ acknowledges support from the ERC-selected UKRI Frontier Research Grant EP/Y03015X/1 and by STFC Consolidated Grant ST/V000780/1. MA is supported by the UK Science and Technology Facilities Council (STFC) under grant number ST/Y002652/1 and the Royal Society under grant numbers RGSR2222268 and ICAR1231094. MvWK acknowledges the support by the Science and Technology Facilities Council. AHW is supported by the Deutsches Zentrum für Luft- und Raumfahrt (DLR), made possible by the Bundesministerium für Wirtschaft und Klimaschutz, and acknowledges funding from the German Science Foundation DFG, via the Collaborative Research Center SFB1491 ‘Cosmic Interacting Matters – From Source to Signal’. MB and PJ are supported by the Polish National Science Center through grant no. 2020/38/E/ST9/00395. MB is supported by the Polish National Science Center through grants no. 2018/30/E/ST9/00698, 2018/31/G/ST9/03388 and 2020/39/B/ST9/03494. PB acknowledges financial support from the Canadian Space Agency (Grant No. 23EXPROSS1) and the Waterloo Centre for Astrophysics. JHD acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1). CG acknowledges support from the project ‘A rising tide: Galaxy intrinsic alignments as a new probe of cosmology and galaxy evolution’ (with project number 1393VI.Vidi.203.011) of the Talent programme Vidi which is (partly) financed by the Dutch Research Council (NWO). SJ acknowledges the Dennis Sciama Fellowship at the University of Portsmouth and the Ramón y Cajal Fellowship from the Spanish Ministry of Science. KK acknowledges support from the Royal Society and Imperial College. LL is supported by the Austrian Science Fund (FWF) [ESP 357-N]. CM acknowledges support under the grant number PID2021-128338NB-I00 from the Spanish Ministry of Science. TT acknowledges funding from the Swiss National Science Foundation under the Ambizione project PZ00P2_193352 MY acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 101053992). Software: The figures in this work were created with matplotlib (Hunter 2007), making use of the NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Cosmosis (Zuntz et al. 2015), Nautilus (Lange 2023) and CosmoPower (Spurio Mancini et al. 2022) software packages. Author contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (RR, SU) followed by two alphabetical groups. The first alphabetical group includes those who are key contributors to both the scientific analysis and the data products. The second group covers those who have either made a significant contribution to the data products, or to the scientific analysis.

---

References

Appendix A: Code structure

The general code structure is depicted in the flowchart of Figure A.1. Dashed black arrows provide input, while solid black arrows indicate inheritance. Each blue square is a class object in the code, while green boxes indicate input files and their description. Here we give a brief explanation of each of the different modules and inputs to provide an idea of the philosophy of the code.

• Input, FileInput, and configuration files

  In the code, users specify all parameters through the configuration file, encompassing cosmological, HOD, bias, precision, and other parameters, as well as settings for considered summary statistics or survey specifications. Certain inputs, such as the redshift distribution of sources and lenses, are handled through input files, which are mandatory. While some settings are mandatory, others have default values. The Input and FileInput classes play a crucial role in managing user input. They communicate with the user to ensure that all mandatory information is provided and that settings are compatible. Upon reading all information from the configuration file via the Input class, variables are stored in dictionaries and passed to the FileInput class. Subsequently, the FileInput class reads in all specified files and stores them in a specified dictionary. These dictionaries containing all necessary information are then passed through the code. It is worth noting that any input-related messages or warnings include factors of h in the units, where applicable.

• Setup

  The Setup class sets the cosmology via the astropy library (Astropy Collaboration 2022) and contains a lot of auxiliary functions which are called later in the code, such as the mode calculation of the survey footprint or the binning of the number of pairs. Furthermore, it contains a lot of the functionality to perform consistency checks of the input parameters, for example whether the k range is large enough to support a required multipole range.

• HOD and HaloModel

  In the HOD class, all the components of the HOD and CSMF (Section 3.2) are defined. This class is then instantiated in the HaloModel class, where all the halo model quantities are implemented (see Section 3.1). For the halo mass function calculations, the hmf package (Murray et al. 2013) is used. All quantities are calculated on a grid in k at a single redshift.

• PolySpectra

  Uses the functionality of the HaloModel class together with the perturbation theory kernels to calculate all power spectra, their responses and trispectra in Fourier space (Section 3.3) for all observables considered at a single redshift. Up to this point, all calculations are completely independent of any survey specifications. Quantities in this class live on a k-grid (k₁ and k₂ in the case of trispectra) and are eventually divided into stellar mass bins if required.

• CovELLSpace

  Carries out the projections from three-dimensional Fourier space to harmonic space as described in Section 4. The code inherits the functionality from the PolySpectra class, updates all the quantities that depend on redshift, and creates splines for the line-of-sight integration. It is here where the survey specifications enter the stage via the redshift distribution, survey mask, for example. The CovELLSpace is also the first class that can produce a covariance matrix. All covariance matrix elements are stored in numpy arrays.

• Integrators

  The red box in the middle of Figure A.1 is an external C++ code which has been wrapped into python for the numerical integration of the highly oscillatory integrals. For the Bessel functions occurring in some of the Fourier weights, we use either cquad from the Gnu Scientific Library (GSL) or Levin's method (Levin 1996; Zieser & Merkel 2016; Spurio Mancini et al. 2018; Leonard et al. 2023). The latter is in particular used for large θ, i.e. when many oscillations of the Bessel function are crossed. In particular, we use cquad only when there are less than 100 extrema of the Bessel functions in the integral domain.

  For most of the Fourier weights, however, Levin's method is not applicable as it requires an analytic relation between the oscillatory part of the integral and its derivative. This is only given via the recurrence relations of the Bessel functions, but not for the Fourier filter functions. Therefore, the integration domain is simply split into different intervals, each covering roughly a specified (in the configuration file) number of extrema of the Filter functions. These integrals are then again evaluated using cquad. The integrator is parallelised over all combinations of tomographic bins and can therefore benefit greatly from a large number of cores. Convergence of the integrals is reached if the final result does not change by a user-specified relative amount.

• CovBandPowers, CovThetaSpace, and CovCOSEBI

  The BandPowers, ThetaSpace, and COSEBI classes all work in the same way. They use the harmonic space expressions and map them into the desired L, θ or n-range respectively (Section 5.3, Section 5.2 and Section 5.4) as specified in the configuration file with all Fourier and real space filter functions computed internally. The main method in this class is always the calc_covbandpowers method (and similarly for the other classes). They provide lists of 10 arrays containing all unique combinations of probes for a standard 3×2 analysis¹⁷ with the different covariance matrix components.

• ARBsummary

  This class works in principle in the same way as the previous classes for the covariance calculation, however, it does not calculate any filters (see Equations 51 and 55) for the summary statistics internally but requires them to be passed as file inputs via the configuration file. It is therefore completely agnostic to the kind of summary statistic which is passed. In particular, you can combine different summary statistics and also check the consistency between summary statistics of different probes via this class by just passing two kinds of Fourier and real space filters for each summary statistic.

• Output

  The last class of the code just takes the list of covariance matrix entries computed by one of the previous four classes and does three things. It writes this long list of all components into a file, depending on the setting they can be split into ‘sva’, ‘mix’, ‘sn’, ‘SSC’, ‘nG’ and ‘cov’ (the sum of all) or just ‘gauss’, ‘SSC’, ‘nG’ and ‘cov’. The order of the elements can be specified so that either the spatial index is the fastest or the slowest varying index. Secondly, the long list is brought into a matrix format, which is also written into a file. Here the order of the indices is always the probe, the tomographic bin index, and the spatial index labelled from the slowest to the fastest. Lastly, this matrix is also used to create a plot of the correlation coefficient.

The code is executed by just running the main script python covariance.py your_config.ini. It will then go through the flowchart depending on the settings you have chosen and communicate via the terminal until it concludes with the output.

Appendix B: Spherical harmonic decomposition

A two-dimensional field $a(\hat{n})$ can be expanded in a spherical harmonic basis:

$$a(\hat{n})={\displaystyle \sum _{\ell ,m}}{a}_{\ell m}^{}{}_{}\mathit{Y}{}_{\ell m}^{}(\hat{n}).$$(B.1)

The spherical harmonics, ${}_{}\mathit{Y}{}_{\ell m}^{}(\hat{n})$, satisfy the orthogonality relation

$$\int \mathrm{d}{\mathrm{\Omega }}_n{}_{}\mathit{Y}{}_{{\ell }_1{m}_1}^{}(\hat{n}){}_{}\mathit{Y}{}_{{\ell }_2{m}_2}^{*}(\hat{n})={\delta }_{{\ell }_1{\ell }_2}^{\mathrm{K}}{\delta }_{m_1{m}_2}^{\mathrm{K}},$$(B.2)

with the solid angle element dΩ_n associated with $\hat{n}$. Furthermore, we use the following Fourier space convention:

$$A(x)={\displaystyle {\int }_k}A(k){\mathrm{e}}^{-\mathrm{i}k·x}≔\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}A(k){\mathrm{e}}^{-\mathrm{i}k·x}.$$(B.3)

With this definition, the Fourier components of statistically homogeneous fields have the following n-point statistics

$$\langle A_1(k_1,{\chi }_1)\dots A_1(k_n,{\chi }_n)\rangle =(2\pi {)}^3{P}_{A_1\dots A_n}(k_1,{\chi }_1,\dots ,k_n,{\chi }_1,\dots ,{\chi }_n){\delta }_{\mathrm{D}}^{(3)}(k_{1\dots n}),$$(B.4)

where k_1…n≔k₁+…+k_n. Using the Rayleigh expansion, the exponential can be written as

$${\mathrm{e}}^{\mathrm{i}k·x}=4\pi {\displaystyle \sum _{\ell ,m}}{\mathrm{i}}^{\ell }{j}_{\ell }^{}(\mathit{kx}){}_{}\mathit{Y}{}_{\ell m}^{}(\hat{k}){}_{}\mathit{Y}{}_{\ell m}^{*}(\hat{n}),$$(B.5)

where x=|x| is the absolute value and $\hat{n}$ the angular part of the vector x, likewise for k.

Appendix C: Polyspectra on the sphere

Expanding the expression using Equation (B.1), Equation (B.4) and Equation (B.5) yields

$$\begin{array}{cc}\hfill \langle a_{1_{{\ell }_1{m}_1}}\dots a_{n_{{\ell }_n{m}_n}}\rangle =& \hspace{0.5em}(\mathrm{i}{)}^{{\ell }_{1·n}}\hspace{0.5em}(2\pi {)}^3\prod _{\alpha =1}^n\left(4\pi \int \mathrm{d}{\chi }_{\alpha }\hspace{0.5em}{W}_{A_{\alpha }}({\chi }_{\alpha })\int \frac{k_{\alpha }^2\mathrm{d}{k}_{\alpha }}{(2\pi {)}^3}\mathrm{d}{\mathrm{\Omega }}_{k_{\alpha }}{j}_{{\ell }_{\alpha }}(k_{\alpha }{\chi }_{\alpha })Y_{{\ell }_{\alpha }{m}_{\alpha }}^{*}({\hat{k}}_{\alpha })\right)\hspace{0.5em}\hfill \\ \hfill & \times P_{A_1\dots A_n}(k_1,\dots ,k_n){\delta }_{\mathrm{D}}^{(3)}(k_{1\dots n}),\hfill \end{array}$$(C.1)

where we carried out the angular integrals $\mathrm{d}{\mathrm{\Omega }}_{n_i}$ and used the orthonormality of the spherical harmonics. Angular integrations over spherical harmonics are trivially evaluated casting Equation (C.1) into a form respecting rotational symmetries. For Equation (8) this is for example realised by the diagonality in ℓ and m. To make further progress, we expand the delta distribution into plane waves as well:

$$\begin{array}{cc}\hfill \langle a_{1_{{\ell }_1{m}_1}}\dots a_{n_{{\ell }_n{m}_n}}\rangle =& \hspace{0.5em}(\mathrm{i}{)}^{{\ell }_{1·n}}(2\pi {)}^3\prod _{\alpha =1}^n\left((4\pi {)}^2{\displaystyle \int }\mathrm{d}{\chi }_{\alpha }\hspace{0.5em}{W}_{A_{\alpha }}({\chi }_{\alpha }){\displaystyle \int }\frac{k_{\alpha }^2\mathrm{d}{k}_{\alpha }}{(2\pi {)}^3}\mathrm{d}{\mathrm{\Omega }}_{k_{\alpha }}{j}_{{\ell }_{\alpha }}(k_{\alpha }{\chi }_{\alpha })j_{{\ell }_{\alpha }}(k_{\alpha }x)Y_{{\ell }_{\alpha }{m}_{\alpha }}^{*}({\hat{k}}_{\alpha })\right)\hspace{0.5em}\hfill \\ \hfill & \times P_{A_1\dots A_n}(k_1,\dots ,k_n)\hspace{0.5em}\prod _{\beta =1}^n\sum _{{\tilde{\ell }}_{\beta }{\tilde{m}}_{\beta }}(\mathrm{i}{)}^{{\tilde{\ell }}_{\beta }}{\mathcal{H}}_{{\ell }_1\dots {\tilde{\ell }}_n}^{{\tilde{m}}_1\dots {\tilde{m}}_n}{\displaystyle \int }{x}^2\mathrm{d}x\hspace{0.5em}{j}_{{\tilde{\ell }}_{\beta }}(k_{\beta }x)Y_{{\tilde{\ell }}_{\beta }{\tilde{m}}_{\beta }}({\hat{k}}_{\beta }).\hfill \end{array}$$(C.2)

Here we defined the integral over spherical harmonics,

$${\mathcal{H}}_{{\ell }_1\dots {\ell }_n}^{m_1\dots m_n}\equiv \int \mathrm{d}{\mathrm{\Omega }}_x{\displaystyle \prod _{\alpha =1}^n}{Y}_{{\ell }_{\alpha }{m}_{\alpha }}(\hat{x}),$$(C.3)

which can be decomposed into integrals over two and three spherical harmonics using

$$Y_{{\ell }_1{m}_1}(\hat{x})Y_{{\ell }_2{m}_2}(\hat{x})={\displaystyle \sum _{\ell m}}{c}_{\ell m}({\ell }_1,{\ell }_2,m_1,m_2)Y_{\ell m}(\hat{x}),$$(C.4)

with the Klebsch-Gordon coefficients

$$c_{\ell m}({\ell }_1,{\ell }_2,m_1,m_2)=(-1{)}^m{\mathcal{G}}_{{\ell }_1{\ell }_2\ell }^{m_1{m}_2-m}$$(C.5)

and the Gaunt integral $G_{{\ell }_1{\ell }_2{\ell }_3}^{m_1{m}_2{m}_3}$:

$${\mathcal{G}}_{{\ell }_1{\ell }_2{\ell }_3}^{m_1{m}_2{m}_3}=\int \mathrm{d}{\mathrm{\Omega }}_x{Y}_{{\ell }_1{m}_1}(\hat{x})Y_{{\ell }_2{m}_2}(\hat{x})Y_{{\ell }_3{m}_3}(\hat{x})=\sqrt{\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{4\pi }}\left(\begin{array}{ccc}\hfill {\ell }_1\hfill & \hfill {\ell }_2\hfill & \hfill {\ell }_3\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{ccc}\hfill {\ell }_1\hfill & \hfill {\ell }_2\hfill & \hfill {\ell }_3\hfill \\ \hfill m_1\hfill & \hfill m_2\hfill & \hfill m_3\hfill \end{array}\right).$$(C.6)

The angular polyspectrum, ${\mathcal{P}}_{{\ell }_1\dots {\ell }_n}^{a_1\dots a_n}$, is then defined via

$$\langle a_{1_{{\ell }_1{m}_1}}\dots a_{n_{{\ell }_n{m}_n}}\rangle ={\mathcal{H}}_{{\ell }_1\dots {\ell }_n}^{m_1\dots m_n}{\mathcal{P}}_{{\ell }_1\dots {\ell }_n}^{a_1\dots a_n}.$$(C.7)

For n=2, 3 the functional form of the angular polyspectrum is completely fixed by rotational symmetry since the integral over Ω_k can always be evaluated.

Appendix D: Flat sky approximation

For computational simplicity, one often relies not on a spherical harmonic decomposition but rather a flat two-dimensional discrete Fourier transform by replacing the directional vector $\hat{n}$ by a unit vector in the flat sky $\theta \perp {\hat{e}}_z$. In this flat sky, one can decompose the projected field a into Fourier modes $\tilde{a}$

$$a(\theta )=\frac{1}{A_{\mathrm{s}}}{\displaystyle \sum _{\ell }}{\tilde{a}}_{\ell }\hspace{0.5em}{\mathrm{e}}^{-\mathrm{i}\ell ·\theta },$$(D.1)

where A_s is the area of the flat sky patch. The sum runs over all ${\ell }^{\mathrm{T}}=2\pi /\sqrt{A_{\mathrm{s}}}(i_x,i_y)$, where $i_x,i_y\in \mathbb{N}$. For A_s→∞ we can replace the discrete by a continuous Fourier transformation, such that

$$a(\theta )=\int {\mathrm{d}}^2\ell {\tilde{a}}_{\ell }\hspace{0.5em}{\mathrm{e}}^{-\mathrm{i}\ell ·\theta }.$$(D.2)

The corresponding angular power spectra are defined via

$$\langle a_{1{\ell }_1}{a}_{2{\ell }_2}\rangle =A_{\mathrm{s}}{\delta }_{{\ell }_1-{\ell }_2}^{\mathrm{K}}{\mathcal{P}}_{a_1{a}_2}({\ell }_1),\langle a_{1{\ell }_1}{a}_{2{\ell }_2}\rangle =(2\pi {)}^2{\delta }^{\mathrm{D}}({\ell }_1-{\ell }_2){\mathcal{P}}_{a_1{a}_2}({\ell }_1)$$(D.3)

for the discrete and continuous cases respectively, where δ^D(x−y) is the Dirac δ-distribution. Higher-order correlators follow analogously.

Appendix E: Halo profiles

For an NFW profile (Navarro et al. 1997) the Fourier transform assumes the form

$$\tilde{u}(k|M,z)=cos({\mathit{kr}}_{\mathrm{s}})\left[\mathrm{Ci}(k(1+c)r_{\mathrm{s}})-\mathrm{Ci}({\mathit{kr}}_{\mathrm{s}})\right]-\frac{sin({\mathit{ckr}}_{\mathrm{s}})}{{\mathit{kr}}_{\mathrm{s}}(1+c)}+\frac{sin({\mathit{kr}}_{\mathrm{s}})\left(\mathrm{Si}({\mathit{kr}}_{\mathrm{s}}(1+c))-\mathrm{Si}({\mathit{kr}}_{\mathrm{s}})\right)}{\frac{1}{1+c}+ln(1+c)-1},$$(E.1)

where Si(x) and Ci(x) are sine and cosine integrals. The scaling radius

$$r_{\mathrm{s}}=\frac{r_{\mathrm{vir}}}{c}={\left(\frac{3M}{4\pi {\mathrm{\Delta }}_{\mathrm{V}}{\overline{\rho }}_{\mathrm{m}}{c}^3}\right)}^{1/3}$$(E.2)

defines the empirical concentration relation c(M); here Δ_V is the virial overdensity. It should be noted that Equation (15) can, in principle, include an integral over the concentration. However, assuming a mean relation between c and M, allows for analytic integration. For the baseline functional form, we assume the relation provided in Duffy et al. (2008):

$$c(M,z)=10.14{\left[\frac{M}{2\times {10}^{12}{h}^{-1}{M}_{\odot }}\right]}^{-0.81}(1+z{)}^{-1.01}.$$(E.3)

Appendix F: Mixed term reconsidered

It should be noted that the equations in Sect 5.2 for the mixed term Equation (66) and the non-Gaussian term assume the tracers to be uniformly distributed. This will not be fulfilled when estimating the statistics in a realistic survey with a nontrivial angular mask and depth variations across the survey footprint. Schneider et al. (2002) have derived the corresponding equations for the shear correlation functions from which we reproduce the mixed term of ${\hat{\xi }}_{+}$ for a tomographic setup,

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{mix}}[{\hat{\xi }}_{+}^{(a_2{a}_2)}({\overline{\theta }}_1);{\hat{\xi }}_{+}^{(a_3{a}_4)}({\overline{\theta }}_2)]\equiv \hspace{0.5em}\frac{{\sigma }_{{ϵ}_1}^2}{N_{\mathrm{pair}}^{(a_1{a}_2)}({\overline{\theta }}_1)\hspace{0.5em}\hspace{0.5em}{N}_{\mathrm{pair}}^{(a_3{a}_4)}({\overline{\theta }}_2)}\_$$(F.1)

$$\times \hspace{0.5em}[{\displaystyle \sum _{i,j,k}^{N_{a_1},N_{a_2},N_{a_3}}}{w}_i^2{w}_j{w}_k\hspace{0.5em}{\xi }_{+}^{(a_2{a}_3)}({\theta }_{\mathit{jk}})\hspace{0.5em}{\mathrm{\Delta }}_{{\theta }_1}({\theta }_{\mathit{ij}}){\mathrm{\Delta }}_{{\theta }_2}({\theta }_{\mathit{ik}})\hspace{0.5em}{\delta }_{a_1{a}_4}^{\mathcal{K}}+\hspace{0.5em}3\hspace{0.5em}\mathrm{perm}.]\hspace{0.5em},$$(F.2)

where the indices i,j,k stand for the individual galaxies having weights w_i, the summation runs over all $N_{a_i}$ galaxies residing in the a_ith tomographic redshift bin and ${\mathrm{\Delta }}_{{\theta }_1}({\theta }_{\mathit{ij}})$ denotes the angular bin selection function which is unity iff θ_ij≡|θ_i−θ_j|∈θ.

For an efficient evaluation of the triplet sums we proceed along the lines of Slepian & Eisenstein (2015) and decompose the three-point correlator in its multipoles:

$$\begin{array}{cc}\hfill M_{++}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2,\varphi )& \equiv \sum _{i,j,k}^{N_{a_1},N_{a_2},N_{a_3}}{\mathit{w}}_i^2{\mathit{w}}_j{w}_k\hspace{0.5em}{\mathrm{\Delta }}_{{\theta }_1}(\mathit{ij})\hspace{0.5em}{\mathrm{\Delta }}_{{\theta }_2}(\mathit{ik})\hspace{0.5em}{\delta }_{{\phi }_{\mathit{kj}}\varphi }^{\mathcal{K}}\hfill \\ \hfill & =\frac{1}{2\pi }\sum _{n=-\infty }^{\infty }{M}_{++,n}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2)\hspace{0.5em}{\mathrm{e}}^{-\mathrm{i}n\varphi }\hspace{0.5em},\hfill \end{array}$$(F.3)

where φ_ji denotes the polar angle of θ_ij, ϕ≡φ_ki−φ_ji and the multipole components can be calculated as

$$\begin{array}{cc}\hfill M_{++,n}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2)& =\sum _i^{N_{a_1}}{\mathit{w}}_i^2\hspace{0.5em}{c}_{i,n}^{(a_1{a}_2)}({\overline{\theta }}_1)\hspace{0.5em}{\left(c_{i,n}^{(a_1{a}_3)}({\overline{\theta }}_2)\right)}^{*}\hspace{0.5em}\hspace{0.5em};\hfill \\ \hfill c_{i,n}^{(a_1{a}_2)}({\overline{\theta }}_1)& \equiv \sum _{j=1}^{N_{a_2}}{\mathit{w}}_j\hspace{0.5em}{\mathrm{e}}^{\mathit{in}{\phi }_{\mathit{ji}}}\hspace{0.5em}{\mathrm{\Delta }}_{{\theta }_1}({\theta }_{\mathit{ij}})\hspace{0.5em}.\hfill \end{array}$$(F.4)

To further speed up the computation we make use of the combined estimator method introduced in Porth et al. (2024) that amends the equations presented here with grid-based methods for which the $c_{i,n}^{(a_1{a}_2)}$ can be computed via FFTs.

Once the multipoles are computed, we can perform the numerical integral to obtain

$$T_{++}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2)\equiv {\displaystyle {{\displaystyle \int }}_0^{2\pi }}\mathrm{d}\varphi \hspace{0.5em}{M}_{++}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2,\varphi )\hspace{0.5em}{\xi }_{+}^{(a_2{a}_3)}({\overline{\theta }}_3)\hspace{0.5em},$$(F.5)

where we evaluate the shear correlation function at ${\overline{\theta }}_3=\sqrt{{\overline{\theta }}_1^2+{\overline{\theta }}_2^2-2{\overline{\theta }}_1{\overline{\theta }}_2\cos (\varphi )}$. Plugging everything together we arrive at an approximation of Equation (F.1) in which ξ₊ is not evaluated for every triplet, but instead the mean distance across the bin is taken¹⁸,

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{mix}}\left[{\hat{\xi }}_{+}^{(a_2{a}_2)}({\overline{\theta }}_1);{\hat{\xi }}_{+}^{(a_3{a}_4)}({\overline{\theta }}_2)\right]\approx \frac{{\sigma }_{{ϵ}_1}^2}{N_{\mathrm{pair}}^{(a_1{a}_2)}({\overline{\theta }}_1)\hspace{0.5em}\hspace{0.5em}{N}_{\mathrm{pair}}^{(a_3{a}_4)}({\overline{\theta }}_2)}\hspace{0.5em}\hspace{0.5em}\left[T_{++}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2)\hspace{0.5em}{\delta }_{a_1{a}_4}^{\mathcal{K}}+\hspace{0.5em}3\hspace{0.5em}\mathrm{perm}.\right]\hspace{0.5em}.$$(F.6)

Performing similar computations we can also work out a multipole-based decomposition for the mixed covariance of ${\hat{\xi }}_{-}$:

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{mix}}\left[{\hat{\xi }}_{-}^{(a_2{a}_2)}({\overline{\theta }}_1);{\hat{\xi }}_{-}^{(a_3{a}_4)}({\overline{\theta }}_2)\right]\approx \frac{{\sigma }_{{ϵ}_1}^2}{N_{\mathrm{pair}}^{(a_1{a}_2)}({\overline{\theta }}_1)\hspace{0.5em}\hspace{0.5em}{N}_{\mathrm{pair}}^{(a_3{a}_4)}({\overline{\theta }}_2)}\hspace{0.5em}\hspace{0.5em}\left[T_{--}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2)\hspace{0.5em}{\delta }_{a_1{a}_4}^{\mathcal{K}}+\hspace{0.5em}3\hspace{0.5em}\mathrm{perm}.\right]\hspace{0.5em},$$(F.7)

where the T₋₋ are given as

$$\begin{array}{cc}\hfill T_{--}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2)& \equiv {{\displaystyle \int }}_0^{2\pi }\mathrm{d}\varphi \hspace{0.5em}{M}_{--}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2,\varphi )\hspace{0.5em}{\xi }_{+}^{(a_2{a}_3)}({\overline{\theta }}_3)\hspace{0.5em};\hfill \\ \hfill M_{--}^{(a_1{a}_2{a}_3)}({\overline{\theta }}_1,{\overline{\theta }}_2,\varphi )& \equiv \frac{1}{2}\sum _{i=1}^{N_{a_1}}{\mathit{w}}_i^2\left\{c_{i,4+n}^{(a_1,a_2)}({\overline{\theta }}_1){\left[c_{i,4+n}^{(a_1,a_3)}({\overline{\theta }}_2)\right]}^{*}+c_{i,n-4}^{(a_1,a_2)}({\overline{\theta }}_1){\left[c_{i,n-4}^{(a_1,a_3)}({\overline{\theta }}_2)\right]}^{*}\right\}\hspace{0.5em}.\hfill \end{array}$$(F.8)

Appendix G: Bandpowers

Band powers are angular averages over the angular power spectrum, Equation (8), and can as such make more use of the easier calculations in harmonic space, while including effects of survey mask and finite coverage of scales more accessible. Generally, they are defined as

$${\mathcal{C}}_{(a_1{a}_2)}(L)≔\frac{1}{{\mathcal{N}}_L}\int \ell \mathrm{d}\ell \hspace{0.5em}{S}_L(\ell )\hspace{0.5em}{\mathcal{P}}_{a_1{a}_2}(\ell ),$$(G.1)

where S_L is the band power response function and ${\mathcal{N}}_L$ is the normalisation

$${\mathcal{N}}_L=\int \frac{\mathrm{d}\ell }{\ell }{S}_L(\ell ),$$(G.2)

such that the band powers trace the angular power ${\ell }^2{\mathcal{P}}_{\ell }^{a_1{a}_2}$ at the log-centered bins in L. Due to the orthogonality of the Bessel functions, Equation (56) and Equation (57) can be inverted to express the band powers in terms of the real space correlation functions w, γ_t and ξ_±. Thus, by using Equation (G.1), one finds (Schneider et al. 2002)

$${\mathcal{C}}_{{\mathrm{n}}_i{\mathrm{n}}_j}(L)=\hspace{0.5em}\frac{\pi }{{\mathcal{N}}_L}\hspace{0.5em}\int \theta \mathrm{d}\theta \hspace{0.5em}T(\theta )w^{(\mathit{ij})}(\theta )g_{+}^L(\theta ),$$(G.3)

$${\mathcal{C}}_{{\mathrm{n}}_i{ϵ}_j}(L)=\hspace{0.5em}\frac{2\pi }{{\mathcal{N}}_L}\hspace{0.5em}\int \theta \mathrm{d}\theta \hspace{0.5em}T(\theta ){\gamma }_{\mathrm{t}}^{(\mathit{ij})}(\theta )h^L(\theta ),$$(G.4)

$${\mathcal{C}}_{{ϵ}_i{ϵ}_j\mathrm{E}/\mathrm{B}}(L)=\hspace{0.5em}\frac{\pi }{{\mathcal{N}}_L}\hspace{0.5em}\int \theta \mathrm{d}\theta \hspace{0.5em}T(\theta )\left[{\xi }_{+}^{(\mathit{ij})}(\theta )\hspace{0.5em}{g}_{+}^L(\theta )\pm {\xi }_{-}^{(\mathit{ij})}(\theta )\hspace{0.5em}{g}_{-}^L(\theta )\right],$$(G.5)

for GC, GGL, and CS respectively. Where the E-mode corresponds to the sum and the B-mode to the difference. With ‘+’ corresponding to J₀, the kernels are given by

$$g_{\pm }^L(\theta )=\hspace{0.5em}\int \ell \hspace{0.5em}\mathrm{d}\ell S_L(\ell )\hspace{0.5em}{\mathrm{J}}_{0/4}(\ell \theta ),$$(G.6)

$$h^L(\theta )=\hspace{0.5em}\int \ell \hspace{0.5em}\mathrm{d}\ell S_L(\ell )\hspace{0.5em}{\mathrm{J}}_2(\ell \theta ).$$(G.7)

If the real-space correlation functions are only accessible over a finite range of angular separations, one apodises the kernels with a Hann window (see Joachimi et al. 2021):

$$T(\theta )=\{\begin{array}{cc}\hspace{0.5em}{cos}^2\left(\frac{\pi [x-(x_{\mathrm{lo}}+{\mathrm{\Delta }}_x/2)]}{2{\mathrm{\Delta }}_x}\right)\hfill & \text{if}{x}_{\mathrm{lo}}-\frac{{\mathrm{\Delta }}_x}{2}\le x<x_{\mathrm{lo}}+\frac{{\mathrm{\Delta }}_x}{2}\hfill \\ \hspace{0.5em}1\hfill & \text{if}{x}_{\mathrm{lo}}+\frac{{\mathrm{\Delta }}_x}{2}\le x<x_{\mathrm{up}}-\frac{{\mathrm{\Delta }}_x}{2}\hfill \\ \hspace{0.5em}{cos}^2\left(\frac{\pi [x-(x_{\mathrm{up}}-{\mathrm{\Delta }}_x/2)]}{2{\mathrm{\Delta }}_x}\right)\hfill & \text{if}{x}_{\mathrm{up}}-\frac{{\mathrm{\Delta }}_x}{2}\le x<x_{\mathrm{up}}+\frac{{\mathrm{\Delta }}_x}{2}\hfill \\ \hspace{0.5em}0\hfill & \text{else}\hfill \end{array},$$(G.8)

with x=logθ and Δx the logarithmic bandwidth. The apodisation is such that the lower and upper limit is centred on x_lo=logθ_lo and x_up=logθ_up, respectively.

Here we chose the band power response to be a top hat between ℓ_lo,l and ℓ_up,l leading to a normalisation ${\mathcal{N}}_l=ln({\ell }_{\mathrm{up},i}/{\ell }_{\mathrm{lo},i})$ and a closed form for $g_{\pm }^L$ and h^L:

$$g_{+}^L(\theta )=\frac{1}{{\theta }^2}\left[\theta {\ell }_{\mathrm{up},L}\hspace{0.5em}{\mathrm{J}}_1(\theta {\ell }_{\mathrm{up},L})-\theta {\ell }_{\mathrm{lo},L}\hspace{0.5em}{\mathrm{J}}_1(\theta {\ell }_{\mathrm{lo},L})\right],$$(G.9)

$$g_{-}^L(\theta )=\frac{1}{{\theta }^2}\left[{\mathcal{G}}_{-}(\theta {\ell }_{\mathrm{up},L})-{\mathcal{G}}_{-}(\theta {\ell }_{\mathrm{lo},L})\right],$$(G.10)

$$h^L=-\frac{1}{{\theta }^2}\left[\theta {\ell }_{\mathrm{up},L}{\mathrm{J}}_1(\theta {\ell }_{\mathrm{up},L})-\theta {\ell }_{\mathrm{lo},L}{\mathrm{J}}_1(\theta {\ell }_{\mathrm{lo},L})+2{\mathrm{J}}_0(\theta {\ell }_{\mathrm{up},L})-2{\mathrm{J}}_0(\theta {\ell }_{\mathrm{lo},L})\right],$$(G.11)

where

$${\mathcal{G}}_{-}(x)=\left(x-\frac{8}{x}\right){\mathrm{J}}_1(x)-8{\mathrm{J}}_2(x).$$(G.12)

Linking the band powers directly to the angular power spectra yields the final expressions

$${\mathcal{C}}_{{\mathrm{n}}_i{\mathrm{n}}_j}(L)=\frac{1}{{\mathcal{N}}_L}\int \ell \hspace{0.5em}\mathrm{d}\ell W_{\mathrm{EE}}^L(\ell )\hspace{0.5em}{\mathcal{P}}_{{\mathrm{n}}_i{\mathrm{n}}_j}(\ell ),$$(G.13)

$${\mathcal{C}}_{{\mathrm{n}}_i{ϵ}_j}(L)=\frac{1}{{\mathcal{N}}_L}\int \ell \hspace{0.5em}\mathrm{d}\ell W_{\mathrm{nE}}^L(\ell )\hspace{0.5em}\mathcal{P}{\mathrm{n}}_i{ϵ}_j(\ell ),$$(G.14)

$${\mathcal{C}}_{{ϵ}_i{ϵ}_j\mathrm{E}}(L)=\frac{1}{{2\mathcal{N}}_L}\int \ell \hspace{0.5em}\mathrm{d}\ell \left[W_{\mathrm{EE}}^L(\ell )\hspace{0.5em}{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{E}}(\ell )+W_{\mathrm{EB}}^L(\ell )\hspace{0.5em}{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{B}}(\ell )\right],$$(G.15)

$${\mathcal{C}}_{{ϵ}_i{ϵ}_j\mathrm{B}}(L)=\frac{1}{2{\mathcal{N}}_L}\int \ell \hspace{0.5em}\mathrm{d}\ell \left[W_{\mathrm{BE}}^L(\ell )\hspace{0.5em}{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{E}}(\ell )+W_{\mathrm{BB}}^L(\ell )\hspace{0.5em}{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{B}}(\ell )\right],$$(G.16)

with kernels given by

$$W_{\mathrm{EE}}^L(\ell )=W_{\mathrm{BB}}^L(\ell )=\int \theta \hspace{0.5em}\mathrm{d}\theta T(\theta )\left[{\mathrm{J}}_0(\ell \theta )\hspace{0.5em}{g}_{+}^L(\theta )+{\mathrm{J}}_4(\ell \theta )\hspace{0.5em}{g}_{-}^L(\theta )\right],$$(G.17)

$$W_{\mathrm{EB}}^L(\ell )=W_{\mathrm{BE}}^L(\ell )=\int \theta \hspace{0.5em}\mathrm{d}\theta T(\theta )\left[{\mathrm{J}}_0(\ell \theta )\hspace{0.5em}{g}_{+}^L(\theta )-{\mathrm{J}}_4(\ell \theta )\hspace{0.5em}{g}_{-}^L(\theta )\right],$$(G.18)

$$W_{\mathrm{nE}}^L(\ell )=\int \theta \hspace{0.5em}\mathrm{d}\theta T(\theta )\hspace{0.5em}{\mathrm{J}}_2(\ell \theta )\hspace{0.5em}{h}^L(\theta ).$$(G.19)

Appendix H: COSEBIs and Ψ statistics

The complete orthogonal sets of E/B-integrals (COSEBIs, Schneider et al. 2010) are summary statistics for CS, avoiding leakage between E- and B-modes on a finite range of angular scales. These are discrete values and can be directly linked to the two-point correlation functions:

$$E_n^{(\mathit{ij})}=\frac{1}{2}{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}[T_{+n}(\theta ){\xi }_{+}^{(\mathit{ij})}(\theta )+T_{-n}(\theta ){\xi }_{-}^{(\mathit{ij})}(\theta )],$$(H.1)

$$B_n^{(\mathit{ij})}=\frac{1}{2}{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}[T_{+n}(\theta ){\xi }_{+}^{(\mathit{ij})}(\theta )-T_{-n}(\theta ){\xi }_{-}^{(\mathit{ij})}(\theta )].$$(H.2)

The filter functions T_±n(θ) are defined for a given angular range, θ∈[θ_min, θ_max]. Two families of COSEBIs have been used in the past, linear-COSEBIs and log-COSEBIs, defining whether the oscillations of T_±(θ) are linearly or logarithmically spaced respectively. The COSEBIs are labelled by n natural numbers, with filters having n+1 roots in their range of support. Here, we employ log-COSEBIs, but the code can use any filter function. The theoretical model prediction for COSEBIs is given by

$$E_n^{(\mathit{ij})}={\int }_0^{\infty }\frac{\ell \hspace{0.5em}\mathrm{d}\ell }{2\pi }{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{E}}(\ell )W_n(\ell ),$$(H.3)

$$B_n^{(\mathit{ij})}={\int }_0^{\infty }\frac{\ell \hspace{0.5em}\mathrm{d}\ell }{2\pi }{\mathcal{P}}_{{ϵ}_i{ϵ}_j,\mathrm{B}}(\ell )W_n(\ell ),$$(H.4)

where the weight functions, W_n(ℓ), are Hankel transforms of T_±(θ)

$$W_n(\ell )={\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}{T}_{+n}(\theta ){\mathrm{J}}_0(\ell \theta )={\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}{T}_{-n}(\theta ){\mathrm{J}}_4(\ell \theta ).$$(H.5)

A very similar set of observables can be constructed for GGL and GC (Buddendiek et al. 2016; Asgari et al. 2021a) to avoid including physical scales outside [θ_min,θ_max]. Similarly to the COSEBIs, they are defined as

$${\mathrm{\Psi }}_m^{{\mathrm{n}}_i{\mathrm{n}}_j}=\hspace{0.5em}{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}{U}_m(\theta )w^{(\mathit{ij})}(\theta )$$(H.6)

$${\mathrm{\Psi }}_m^{{\mathrm{n}}_i{ϵ}_j}=\hspace{0.5em}{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}{Q}_m(\theta ){\gamma }_{\mathrm{t}}^{(\mathit{ij})}(\theta ),$$(H.7)

so that U_m(θ) is compensated and orthogonal and Q_m(θ) given by the expression

$$Q_m(\theta )=\frac{2}{{\theta }^2}{\int }_0^{\theta }{\theta }^{\prime }\mathrm{d}{\theta }^{\prime }{U}_m({\theta }^{\prime })-U_m(\theta ).$$(H.8)

Their respective Fourier counterparts are calculated in the same manner as the COSEBIs E-mode, Equation (H.3), with a Fourier weight:

$$W_m^{\mathrm{\Psi }}(\ell )={\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}{U}_m(\theta ){\mathrm{J}}_0(\ell \theta )={\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\theta \hspace{0.5em}\mathrm{d}\theta \hspace{0.5em}{Q}_m(\theta ){\mathrm{J}}_2(\ell \theta ).$$(H.9)

Appendix I: Definition of Fourier and real space filters for summary statistics

As an example, we define the linear maps W_L(ℓ) and R_L(θ) for the three summary statistics considered in Section 5. For real space correlation functions, one finds the following:

$${\text{W}}_L^{2\mathrm{pcf}}(\ell )=\left(\begin{array}{cccc}\hfill {\mathcal{K}}_0(\ell L)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathcal{K}}_2(\ell L)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathcal{K}}_0(\ell L)\hfill & \hfill {\mathcal{K}}_0(\ell L)\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathcal{K}}_4(\ell L)\hfill & \hfill -{\mathcal{K}}_4(\ell L)\hfill \end{array}\right)\text{and}{\text{R}}_L^{2\mathrm{pcf}}(\theta )=1\frac{1}{\theta }\pi \left(\frac{\theta -L}{\mathrm{\Delta }L}\right).$$(I.1)

Here π is a top-hat filter centred around L with bandwidth ΔL. This includes the averaging over θ bins in the weights. For bandpowers one finds

$${\text{W}}_L^{\mathrm{BP}}(\ell )=\frac{\pi }{{\mathcal{N}}_L}\left(\begin{array}{cccc}\hfill 2W_{\mathrm{EE}}^L(\ell )\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 2W_{\mathrm{nE}}^L(\ell )\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill W_{\mathrm{EE}}^L(\ell )\hfill & \hfill W_{\mathrm{EB}}^L(\ell )\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill W_{\mathrm{EB}}^L(\ell )\hfill & \hfill W_{\mathrm{BB}}^L(\ell )\hfill \end{array}\right)\text{and}{\text{R}}_L^{\mathrm{BP}}(\theta )=\frac{\pi }{{\mathcal{N}}_L}T(\theta )\left(\begin{array}{cccc}\hfill g_{+}^L(\theta )\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 2h^L(\theta )\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill g_{+}^L(\theta )\hfill & \hfill g_{-}^L(\theta )\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill g_{+}^L(\theta )\hfill & \hfill -g_{-}^L(\theta )\hfill \end{array}\right),$$(I.2)

and for COSEBIs:

$${\text{W}}_L^{\mathrm{COSEBI}}(\ell )=\left(\begin{array}{cccc}\hfill W_L^{\mathrm{\Psi }}(\ell )\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill W_L^{\mathrm{\Psi }}(\ell )\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill W_L(\ell )\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill W_L(\ell )\hfill \end{array}\right)\text{and}{\text{R}}_L^{\mathrm{COSEBI}}(\theta )=\frac{1}{2}\left(\begin{array}{cccc}\hfill 2U_L(\theta )\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 2Q_L(\theta )\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill T_{+L}(\theta )\hfill & \hfill T_{-L}(\theta )\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill T_{+L}(\theta )\hfill & \hfill -T_{-L}(\theta )\hfill \end{array}\right),$$(I.3)

We now label components of vec(C)(ℓ) with Latin indices and similar for W_L(ℓ) and the summary statistic $\mathrm{vec}(\text{O})(L)$. Using Equation (51) and Equation (45), Gaussian covariances which are not pure shot noise terms can be expressed as

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sva}}\left[{\mathcal{O}}_{p_1,L_1}^{\mathit{ij}}{\mathcal{O}}_{p_2,L_2}^{\mathit{mn}}\right]=\hspace{0.5em}\frac{1}{2\pi \hspace{0.5em}\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\sum _{q_1,q_2}\int \ell \mathrm{d}\ell {\left({\text{W}}_{L_1}(\ell )\right)}_{p_1{q}_1}{\left({\text{W}}_{L_2}(\ell )\right)}_{p_2{q}_2}$$(I.4)

$$\times {\left[{\mathcal{P}}^{(\mathit{im})}(\ell ){\mathcal{P}}^{(\mathit{jn})}(\ell )+{\mathcal{P}}^{(\mathit{in})}(\ell ){\mathcal{P}}^{(\mathit{jm})}(\ell )\right]}_{q_1{q}_2}$$(I.5)

where it is understood that the summary statistic i is the one identified with the tracer combination a₁,a₂ and the same of j and a₃,a₄. The shot noise term can then be calculated as

$${\mathrm{Cov}}_{\mathrm{G},\mathrm{sn}}\left[{\mathcal{O}}_{p_1,L_1}^{\mathit{ij}}{\mathcal{O}}_{p_2,L_2}^{\mathit{mn}}\right]=\sum _q{N}_q^{\mathit{ij}}\int \frac{{\theta }^2\mathrm{d}\theta }{n_q^{\mathit{ij}}(\theta )}{\left({\text{R}}_{L_1}(\theta )\right)}_{p_{1}q}{\left({\text{R}}_{L_2}(\theta )\right)}_{p_{2}q}\left[{\delta }_{\mathit{im}}^{\mathrm{K}}{\delta }_{\mathit{jn}}^{\mathrm{K}}+{\delta }_{\mathit{in}}^{\mathrm{K}}{\delta }_{\mathit{jm}}^{\mathrm{K}}\right],$$(I.6)

where $n_q^{\mathit{ij}}(\theta )$ is the differential pair-count density and $N_q^{\mathit{ij}}$ the noise level in the q-th summary statistic and the i,j tomographic bin combination. Lastly, connected terms such as NG and SSC are given by

$${\mathrm{Cov}}_{\mathrm{NG},\mathrm{sn}}\left[{\mathcal{O}}_{p_1,L_1}^{\mathit{ij}}{\mathcal{O}}_{p_2,L_2}^{\mathit{mn}}\right]=\hspace{0.5em}\frac{1}{4{\pi }^2\max (A_{(\mathit{ij})}{A}_{(\mathit{mn})})}\sum _{q_1,q_2}\int \mathrm{d}{\ell }_1{\ell }_1\hspace{0.5em}\int \mathrm{d}{\ell }_2{\ell }_2\hspace{0.5em}{\left({\text{W}}_{L_1}({\ell }_1)\right)}_{p_1{q}_1}{\left({\text{W}}_{L_2}({\ell }_2)\right)}_{p_2{q}_2}$$(I.7)

$$\begin{array}{cc}\hfill & \times {\displaystyle {\int }_0^{\pi }}\frac{\mathrm{d}{\varphi }_{\ell }}{\pi }\hspace{0.5em}{T}_{q_1{q}_2}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,-{\ell }_1,{\ell }_2,-{\ell }_2),\hfill \\ \hfill {\mathrm{Cov}}_{\mathrm{SSC},\mathrm{sn}}\left[{\mathcal{O}}_{p_1,L_1}^{\mathit{ij}}{\mathcal{O}}_{p_2,L_2}^{\mathit{mn}}\right]=& \hspace{0.5em}\frac{1}{4{\pi }^2}\int \mathrm{d}{\ell }_1{\ell }_1\hspace{0.5em}\int \mathrm{d}{\ell }_2{\ell }_2\hspace{0.5em}{\left({\text{W}}_{L_1}({\ell }_1)\right)}_{p_1{q}_1}{\left({\text{W}}_{L_2}({\ell }_2)\right)}_{p_2{q}_2}{\left(T_{\mathrm{SSC}}^{(\mathit{ij})(\mathit{mn})}({\ell }_1,{\ell }_2)\right)}_{q_1{q}_2}.\hfill \end{array}$$(I.8)

All Tables

All Figures