# KiDS-Legacy: Covariance validation and the unified OneCovariance framework for projected large-scale structure observables

Robert Reischke, Sandra Unruh, Marika Asgari, Andrej Dvornik, Hendrik

Hildebrandt, Benjamin Joachimi, Lucas Porth, Maximilian von

Wietersheim-Kramsta, Jan Luca van Den Busch, Benjamin Stölzner, et al.

► To cite this version:

Robert Reischke, Sandra Unruh, Marika Asgari, Andrej Dvornik, Hendrik Hildebrandt, et al. KiDS-Legacy: Covariance validation and the unified OneCovariance framework for projected large-scale structure observables. Astronomy & Astrophysics - A&A, 2025, 699, pp.A124. 10.1051/0004-6361/202452592. hal-04757910

HAL Id: hal-04757910

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

Submitted on 16 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# KiDS-Legacy: Covariance validation and the unified ONECOVARIANCE framework for projected large-scale structure observables

Robert Reischke<sup>1,2,*</sup>, Sandra Unruh<sup>1</sup>, Marika Asgari<sup>3,7,8</sup>①, Andrej Dvornik<sup>2</sup>①, Hendrik Hildebrandt<sup>2</sup>①,
Benjamin Joachimi<sup>4</sup>①, Lucas Porth<sup>1</sup>①, Maximilian von Wietersheim-Kramsta<sup>4,5,6</sup>, Jan Luca van den Busch<sup>2</sup>①,
Benjamin Stölzner<sup>2</sup>, Angus H. Wright<sup>2</sup>①, Ziang Yan<sup>2</sup>, Maciej Bilicki<sup>9</sup>①, Pierre Burger<sup>10,1</sup>①, Nora Elisa Chisari<sup>11,15</sup>①,
Joachim Harnois-Déraps<sup>3</sup>①, Christos Georgiou<sup>11</sup>, Catherine Heymans<sup>2,12</sup>, Priyanka Jalan<sup>9</sup>①, Shahab Joudaki<sup>13,14</sup>,
Konrad Kuijken<sup>15</sup>①, Shun-Sheng Li<sup>15,16</sup>①, Laila Linke<sup>17</sup>①, Constance Mahony<sup>18,2</sup>①, Davide Sciotti<sup>19,20</sup>①,
Tilman Tröster<sup>21</sup>①, and Mijin Yoon<sup>15</sup>

*(Affiliations can be found after the references)*

Received 12 October 2024 / Accepted 16 March 2025

## ABSTRACT

We introduce ONECOVARIANCE, an open-source software designed to accurately compute covariance matrices for an arbitrary set of two-point summary statistics across a variety of large-scale structure tracers. Utilising the halo model, we estimated the statistical properties of matter and biased tracer fields, incorporating all Gaussian, non-Gaussian, and super-sample covariance terms. The flexible configuration permits user-specific parameters, such as the complexity of survey geometry, the halo occupation distribution employed to define each galaxy sample, or the form of the real-space and/or Fourier space statistics to be analysed. We illustrate the capabilities of ONECOVARIANCE within the context of a cosmic shear analysis of the final data release of the Kilo-Degree Survey (KiDS-Legacy). Upon comparing our estimated covariance with measurements from mock data and calculations from independent software, we ascertain that ONECOVARIANCE achieves accuracy at the per cent level. When assessing the impact of ignoring complex survey geometry in the cosmic shear covariance computation, we discover misestimations at approximately the 10% level for cosmic variance terms. Nonetheless, these discrepancies do not significantly affect the KiDS-Legacy recovery of cosmological parameters. We derive the cross-covariance between real-space correlation functions, bandpowers, and COSEBIs, facilitating future consistency tests among these three cosmic shear statistics. Additionally, we calculate the covariance matrix of photometric-spectroscopic galaxy clustering measurements, validating the jackknife covariance estimates for calibrating KiDS-Legacy redshift distributions. The ONECOVARIANCE can be found on GitHub (https://github.com/rreischke/OneCovariance) together with comprehensive documentation and examples.

**Key words.** cosmological parameters – cosmology: observations – cosmology: theory – large-scale structure of Universe

## 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)<sup>1</sup>, the Dark Energy Survey (DES, see e.g. Abbott et al. 2018; Amon et al. 2022; Abbott et al. 2022)<sup>2</sup>, and the Subaru Hyper Suprime-Cam (HSC, see e.g. Sugiyama et al. 2022; Hamana et al. 2020; Dalal et al. 2023)<sup>3</sup>. These experiments have laid the groundwork for upcoming stage 4 surveys such as *Euclid*<sup>4</sup>, The Legacy Survey of Space and Time (LSST)<sup>5</sup>, and *Roman*<sup>6</sup>. 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

<footnote>* Corresponding author: rreischke@astro.uni-bonn.de, reischke@posteo.net</footnote>

<footnote n="1">1 https://kids.strw.leidenuniv.nl/</footnote>

<footnote n="2">2 https://www.darkenergysurvey.org/</footnote>

<footnote n="3">3 https://hsc.mtk.nao.ac.jp/ssp/</footnote>

<footnote n="4">4 https://www.euclid-ec.org/</footnote>

<footnote n="5">5 https://www.lsst.org/</footnote>

<footnote n="6">6 https://roman.gsfc.nasa.gov/</footnote>
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^2-10^3$ data points, while upcoming stage 4 surveys will boast $10^4$ data points for standard two-point statistics. Consequently, significant data compression (e.g. Heavens et al. 2000; Asgari & Schneider 2015; Heavens et al. 2017) 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; Becker et al. 2016; Hikage et al. 2019; 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(\theta)$, $\gamma_t(\theta)$, and $\xi_{\pm}(\theta)$ for GC, GGL, and CS respectively<fnref n="7" /> 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<fnref n="8" />) 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.

<footnote n="7">⁷ See Equations (59), (58) and (56) for their respective definitions.</footnote>

<footnote n="8">⁸ Ψ-stats are the COSEBIs-equivalent for GGL and GC equivalent.</footnote>
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[\delta]$, which are functionals of the underlying matter field, $\delta(\mathbf{x}, \chi)$. Here, $\mathbf{x}$ represents a spatial co-moving vector (i.e. a 3-dimensional vector in a spatial hyper-surface), and $\chi$ denotes the co-moving distance. These two-point functions can be aggregated into a vector $\mathcal{O}$ with components

$$O_{\alpha} := \langle t_i[\delta] t_j[\delta] \rangle \quad \alpha \in \{\text{all unique combinations of } (t_i, t_j)\}. \qquad (1)$$

Here, the angular brackets, $\langle\cdot\rangle$, 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 \le j$, otherwise, all combinations can exist. The components of the covariance matrix are then defined as

$$(\mathbf{C})_{\alpha\beta} \equiv \mathrm{Cov}(O_{\alpha}, O_{\beta}) := \langle (O_{\alpha} - \langle O_{\alpha} \rangle)(O_{\beta} - \langle O_{\beta} \rangle) \rangle. \qquad (2)$$

Typically, we are provided with the statistical properties of the density field $\delta$ from theoretical calculations or numerical simulations, rather than those of its tracer. If the mapping from $\delta$ to $t_i$ is linear, all the statistical properties of $\delta$ 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 $\langle t_i[\delta] \rangle = 0$, resulting in the following breakdown of Equation (2) using Equation (1) and Wick's theorem:

$$C_{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_c. \qquad (3)$$

Here, we have omitted the argument of the different tracers, $[\delta]$, 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_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}}. \qquad (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, \qquad (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 \underbrace{\langle t_i t_m \rangle \langle t_j t_n \rangle}_{\text{sample variance}} + \underbrace{\langle t_i t_m \rangle \langle n_j n_n \rangle}_{\text{mixed}} + \underbrace{\langle t_j t_n \rangle \langle n_i n_m \rangle}_{\text{shot/shape noise}} + \underbrace{\langle n_j n_n \rangle \langle n_i n_m \rangle}_{\text{shot/shape noise}}, \qquad (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(\mathbf{x}, f_K\chi(z)) \in \{t_i\}$ at position $\mathbf{x}$ and redshift $z$. Its projected version, $a(\hat{\mathbf{n}})$, is

$$a(\hat{\mathbf{n}}) = \int d\chi W_A(\chi) A(f_K(\chi)\hat{\mathbf{n}}, \chi), \qquad (7)$$

where $W_A(\chi)$ is some weighting function whose exact form is probe specific and will be discussed in Section 4. $\hat{\mathbf{n}}$ is the unit vector on the sphere and $f_K(\chi)$ is the co-moving angular diameter distance, such that $\mathbf{x} = f_K(\chi)\hat{\mathbf{n}}$ and $\chi$ 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, $P_{a_1a_2}(\ell_1)$, of the two-dimensional fields $a_1$ and $a_2$ as follows:

$$\begin{aligned} \langle a_{1,\ell_1 m_1} a_{2,\ell_2 m_2} \rangle &:= P_{a_1 a_2}(\ell_1) \delta_{\ell_1 \ell_2}^{K} \delta_{m_1 m_2}^{K} \\ &= \frac{2}{\pi} \int d\chi_1 W_{a_1}(\chi_1) \int d\chi_2 W_{a_2}(\chi_2) \int k^2 dk P_{A_1 A_2}(k, \chi_1, \chi_2) j_\ell(k f_K(\chi_1)) j_\ell(k f_K(\chi_2)) \delta_{\ell_1 \ell_2}^{K} \delta_{m_1 m_2}^{K}. \end{aligned} \qquad (8)$$

Here the dependence on $\chi_1$ and $\chi_2$ of the three-dimensional power spectrum, $P_{A_1A_2}(k, \chi_1, \chi_2)$, was made explicit. $\delta_{ij}^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_1A_2}(k, \chi_1, \chi_2) \approx [P_{A_1A_2}(k, \chi_1) P_{A_1A_2}(k, \chi_2)]^{1/2}, \qquad (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 \int \frac{d\chi}{f_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_K(\chi)}, \dots, \frac{\ell_n + 0.5}{f_K(\chi)}, \chi \right). \qquad (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<fnref n="9" /> (see Mead et al. 2020, and references therein), we stick to a halo model approach in

<footnote n="9">⁹ It should be noted that the vanilla halo model requires some more attention to accurately fit the transition region between the one-halo and two-halo term. Furthermore, for higher-order statistics these tweaks become more involved (see Takahashi et al. 2020 for the bispectrum). However, for the trispectrum contribution to the covariance, the accuracy of the halo model is, at least currently, sufficient.</footnote>
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; Lacasa et al. 2014; Lacasa 2018; Reischke et al. 2020; Dvornik et al. 2023). To achieve this, one employs a biasing function, $b(M)$, depending on halo mass $M$, defined via the average number density of galaxies $\langle N|M \rangle$ in a halo of mass $M$

$$b(M) := \frac{\bar{\rho}_m(z)}{\bar{n}_g(z)} \frac{\langle N|M \rangle}{M}, \qquad (11)$$

where $\bar{n}_g$ is the average number density of galaxies and $\bar{\rho}_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 := \sum_{N=0}^{\infty} N P(N|M), \qquad (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) := \int dM \, n_h(M, z) A(M, z). \qquad (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 $\rho_h(r|M) = M u_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 $\langle N_{c/s} | M \rangle$. If there is a central galaxy in a halo on average, $\langle N_c | M \rangle = 1$ the distribution of additional galaxies, i.e. the satellites, is assumed to be Poissonian. This allows us to calculate the expected number of $n$-tuples, $\langle N^{(n)} | M \rangle$ in a halo, as required for the 1, ..., 4-halo terms

$$\begin{align}
\langle N|M \rangle &= \langle N_c|M \rangle + \langle N_s|M \rangle & n=1, \nonumber \\
\langle N(N-1)|M \rangle &= \langle N_c|M \rangle (\langle N_s|M \rangle^2 + 2\langle N_s|M \rangle) & n=2, \nonumber \\
\langle N(N-1)(N-2)|M \rangle &= \langle N_c|M \rangle (\langle N_s|M \rangle^3 + 3\langle N_s|M \rangle) & n=3, \nonumber \\
\dots, & & \nonumber
\end{align} \qquad (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 $\bar{\rho}_m$, i.e. $\langle N_m | M \rangle = 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_1, \dots, X_\mu$:

$$I_{\mu, X_{[\mu]}}^{\alpha}(k_1, \dots, k_{\mu}, z) := \frac{1}{\bar{n}_{X_1}(z) \dots \bar{n}_{X_{\mu}}(z)} \int dM \, n_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). \qquad (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):

$$\bar{n}_X(z) = \int dM \, n_h(M, z) \langle N_X | M \rangle. \qquad (16)$$

We note again that $\bar{n}_m(z) = \bar{\rho}_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

$$ \Phi(M_{\star}|M) = \Phi_{c}(M_{\star}|M) + \Phi_{s}(M_{\star}|M). \qquad (17) $$

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

$$ \begin{aligned} \Phi_c(M_{\star}|M) &= \frac{1}{\sqrt{2\pi}\ln(10)\sigma_c M_{\star}} \exp\left[-\frac{\log(M_{\star}/M_c^*)^2}{2\sigma_c^2}\right], \\ \Phi_s(M_{\star}|M) &= \frac{\phi_s^*}{M_s^*} \left(\frac{M_{\star}}{M_s^*}\right)^{\alpha_s} \exp\left[-\left(\frac{M_{\star}}{M_s^*}\right)^2\right], \end{aligned} \qquad (18) $$

for centrals and satellites respectively<fnref n="10" />. Here $\sigma_c$ describes the scatter of stellar mass given a halo mass while $\alpha_s$ describes the power law scaling of the satellite abundance. It should be noticed that if $\alpha_s < 0$ the expression for $\Phi_s(M_{\star}|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 $\Phi_c(M_{\star}|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{aligned} M_c^* &= M_0 \frac{(M/M_1)^{\gamma_1}}{(1 + [M/M_1]^{\gamma_1 - \gamma_2})}, \\ M_s^* &= 0.56 M_c^*(M), \\ \log \phi_s^* &= b_0 + b_1 \log \left( \frac{M}{10^{13} M_{\odot} h^{-1}} \right). \end{aligned} \qquad (19) $$

Integrating the CSMF over the HMF yields the galaxy SMF:

$$ \Phi_X(M_{\star}) = \int dM \, \Phi_X(M_{\star}|M) n_h(M, z). \qquad (20) $$

Likewise, the occupation number in a stellar mass bin [$M_{\star,1}$, $M_{\star,2}$] is given by

$$ \langle N_X | M \rangle = \int_{M_{\star,1}}^{M_{\star,2}} dM_{\star} \Phi_X(M_{\star} | M), \qquad (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}^{1h}(k, z) + P_{A_1 A_2}^{2h}(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_{\text{lin}}(k, z), \qquad (22) $$

with the 1- and 2-halo term as well as the linear matter power spectrum $P_{\text{lin}}$. The one-halo term will tend to be a constant, as $k \to 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}^{1h}(k, z) \to P_{A_1 A_2}^{1h}(k, z) \text{erf}(k/k_{\text{damp}}), \qquad (23) $$

with a fiducial damping scale of $k_{\text{damp}} = 0.1$ h/Mpc. Likewise, the halo model trispectrum is split into 1-4-halo terms

$$ T_{X_1 \dots X_4}(k_1, k_2, k_3, k_4, z) = T_{A_1 \dots A_4}^{1h} + T_{A_1 \dots A_4}^{2h} + T_{A_1 \dots A_4}^{3h} + T_{A_1 \dots A_4}^{4h}, \qquad (24) $$

where it is understood that the right-hand side depends as well on $(k_1, k_2, k_3, k_4, z)$. The individual terms are given by

$$ T_{A_1 \dots A_4}^{1h} = I_{4,A_1 \dots A_4}^{0}(k_1, k_2, k_3, k_4), \qquad (25) $$

$$ T_{A_1 \dots A_4}^{2h} = P_{\text{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 \text{ perm.} + P_{\text{lin}}(k_1) I_{2,A_1 A_2}^{1}(k_1, k_2) I_{2,A_3 A_4}^{1}(k_3, k_4) + 2 \text{ perm.}, \qquad (26) $$

$$ T_{A_1 \dots A_4}^{3h} = B_{\text{tree}}(k_1, k_2, k_3) 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_{\text{lin}}(k_1) P_{\text{lin}}(k_2) \\ + 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 \text{ perm.}, \qquad (27) $$

$$ T_{A_1 \dots A_4}^{4h} = \prod_{\mu} I_{1,A_1 \mu}^{1}(k_{\mu}) T_{\text{tree}}(k_1, k_2, k_3, k_4). \qquad (28) $$

<footnote n="10">¹⁰ We refer to the natural logarithm as ln(x) and to the logarithm to base ten as log(x).</footnote>
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_{\text{tree}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = 2F_2^s(\mathbf{k}_1, \mathbf{k}_2)P_{\text{lin}}(k_1)P_{\text{lin}}(k_2) + 2 \text{perm.}, \qquad (29)$$

$$T_{\text{tree}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) = 4 \left[ F_2^s(\mathbf{k}_{12}, -\mathbf{k}_1) F_2^s(\mathbf{k}_{12}, \mathbf{k}_3) P_{\text{lin}}(k_1) P_{\text{lin}}(k_{12}) P_{\text{lin}}(k_3) + 12 \text{perm.} \right] \\ + 6 \left[ P_{\text{lin}}(k_1) P_{\text{lin}}(k_2) P_{\text{lin}}(k_3) F_3^s(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) + 4 \text{perm.} \right], \qquad (30)$$

where we omitted the explicit dependence on redshift. For the symmetric perturbation theory kernels, $F_i^s$, we refer to Bernardeau et al. (2002). The trispectrum relates to the nG term in Equation (4) for $\mathbf{k}_1 = -\mathbf{k}_2$ and $\mathbf{k}_2 = -\mathbf{k}_4$. 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}^{\text{SSC}}(k_1, k_2, z) = \frac{\partial P_{A_1 A_2}}{\partial \delta_{\text{bg}}}(k_1, z) \frac{\partial P_{A_3 A_4}}{\partial \delta_{\text{bg}}}(k_2, z) \sigma_{\text{bg}, (A_1 A_2) (A_3 A_4)}^2(z). \qquad (31)$$

We note that $k_1$ corresponds to the correlated pair $A_1, A_2$ and $k_2$ to $A_3, A_4$. The responses are given as

$$\frac{\partial P_{A_1 A_2}}{\partial \delta_{\text{bg}}}(k, z) = \left[ \frac{68}{21} - \frac{1}{3} \frac{\text{d} \log k^3 P_{A_1 A_2}(k, z)}{\text{d} \log k} \right] I_{1, A_1}^1(k, z) I_{1, A_2}^1(k, z) P_{\text{lin}}(k, z) \\ + I_{2, A_1 A_2}^1(k, z) - ([b_{A_1} + b_{A_2}] P_{A_1 A_2}(k, z))_{\text{if } A_1, A_2 \neq m}. \qquad (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_{\text{bg}, (A_1 A_2) (A_3 A_4)}^2(z) = \frac{1}{A_{(A_1 A_2)} A_{(A_3 A_4)}} \sum_{\ell} P_{\text{lin}}(\ell/\chi) \sum_{m} a_{\ell m}^{*(A_1 A_2)} a_{\ell m}^{(A_3 A_4)}. \qquad (33)$$

The bracket notation $(A_1 A_2)$ denotes the footprint of the survey over which the summary statistic between the tracers $A_1$ and $A_2$ has been evaluated. As expected $\sigma_{\text{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_{\text{bg}, (A_1 A_2) (A_3 A_4)}^2(z) = \chi(z)^2 \int \frac{dk}{(2\pi)^2} P_{\text{lin}}(k, z) \left[ \frac{2J_1(k\chi(z)\theta_s)}{k\chi(z)\theta_s} \right]^2, \qquad (34)$$

with a cylindrical Bessel function, $J_1$, and the survey size $\theta_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_{g_i}(\chi) = n_1^{(i)}(\chi), \qquad (35)$$

$$W_{m_i}(\chi) = \frac{3 f_K(\chi) \Omega_m}{2 \chi_H^2 a} \int_{\chi}^{\chi_H} d\chi' \frac{f_K(\chi' - \chi)}{f_K(\chi')} n_s^{(i)}(\chi) - A_{\text{IA}} \left[ \frac{1+z(\chi)}{1+z_{\text{pivot}}} \right]^{\eta_{\text{IA}}} \frac{C_1 \rho_{\text{cr}} \Omega_m}{D_+ [a(\chi)]} n_s^{(i)}(\chi). \qquad (36)$$

Here $\chi_H$ is the co-moving distance to the horizon, $n_1^{(i)}(\chi)$ and $n_s^{(i)}(\chi)$ are the normalised redshift distributions of the lens- and source sample respectively. $D_+(a)$ is the linear growth factor and $C_1\rho_{cr} \approx 0.0134$. The second term in the $W_{m_i}$ weight function corresponds to the non-linear linear alignment model (NLA, Hirata & Seljak 2004; Bridle & King 2007) with alignment amplitude $A_{\text{IA}}$ and redshift dependence $\eta_{\text{IA}}$ at a pivot redshift $z_{\text{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_1$ 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}, \qquad (37)$$
with $\langle a_{i,\ell m} n_{j,\ell'm'} \rangle = 0$ and noise statistic $\langle n_{i,\ell m} n_{j,\ell'm'} \rangle = N_i \delta_{ij}^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 $N_i$ determines the overall noise level depending on the observable under consideration:

$$ N_i = \begin{cases} \frac{\sigma_{\epsilon_1, i}^2}{\bar{n}_i} & \text{if } i \in \text{source} \\ \frac{1}{\bar{n}_i} & \text{if } i \in \text{lens} \end{cases}. \qquad (38) $$

Here $\bar{n}_i$ is the average number density of the tracers and $\sigma_{\epsilon_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)} \sum_{m=-\ell}^{\ell} \hat{a}_{\ell m}^{*} \hat{a}_{\ell m}^{2}, \qquad (39) $$

such that

$$ \langle \hat{C}_{a_1 a_2}(\ell) \rangle = \mathcal{P}_{a_1 a_2}(\ell) + N_{a_1} \delta_{a_1 a_2}^{K} \equiv C_{a_1 a_2}(\ell). \qquad (40) $$

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

$$ \hat{C}_{a_1 a_2}(\ell) := \frac{1}{A_s N_{\text{modes}}(\ell)} \sum_{\ell' \in \ell_{\text{ring}}} \hat{a}_{\ell}^{1} \hat{a}_{-\ell}^{2}, \qquad (41) $$

where we defined an annulus in Fourier space with volume $2\pi\ell\Delta\ell$, so that the number of available modes over the survey area $A_s$ is

$$ N_{\text{modes}}(\ell) = \frac{2\pi\ell\Delta\ell}{(2\pi)^2/A_s} = 2\ell\Delta\ell f_{\text{sky}}, \qquad (42) $$

where $\Delta\ell \ll \ell$ was assumed. The covariance matrix between two estimators is then, again by noting that $\ell_1$ corresponds to the pair $a_1, a_2$ and $\ell_2$ corresponds to the pair $a_3, a_4$:

$$ \begin{aligned} \text{Cov}[\hat{C}_{a_1 a_2}(\ell_1) \hat{C}_{a_3 a_4}(\ell_2)] &= \frac{1}{A_{(12)}N_{\text{modes}}(\ell_1)} \frac{1}{A_{(34)}N_{\text{modes}}(\ell_2)} \sum_{\tilde{\ell}_1 \in \ell_{1,\text{ring}}} \sum_{\tilde{\ell}_2 \in \ell_{2,\text{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) \\ &= \frac{1}{N_{\text{modes}}(\ell_1)} \frac{1}{N_{\text{modes}}(\ell_2)} \sum_{\tilde{\ell}_1 \in \ell_{1,\text{ring}}} (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)) \delta_{\ell_1 \ell_2}^K \\ &\quad + \frac{1}{\max(A_{(12)}A_{(34)}) N_{\text{modes}}(\ell_1) N_{\text{modes}}(\ell_2)} \sum_{\tilde{\ell}_1 \in \ell_{1,\text{ring}}} \sum_{\tilde{\ell}_2 \in \ell_{2,\text{ring}}} T^{(a_1 a_2)(a_3 a_4)}(\ell_1, -\ell_1, \ell_2, -\ell_2). \end{aligned} \qquad (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_1 a_2)(a_3 a_4)}(\ell_1, -\ell_1, \ell_2, -\ell_2)$ 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 $\ell_{\text{ring}}$, the Gaussian term in the flat sky approximation is given by the commonly used expression

$$ \text{Cov}_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}^K}{\max(A_{(12)}A_{(34)})2\ell_1\Delta\ell_1} [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)]. \qquad (44) $$

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

$$ \text{Cov}_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}^K}{\max(A_{(12)}A_{(34)})(2\ell_1 + 1)} [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)], \qquad (45) $$

where we introduced the sky fraction $\max(A_{(12)}A_{(34)})/(4\pi)$ 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

$$ \text{Cov}_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}^K}{\max(A_{(12)}A_{(34)})(2\ell_1 + 1)\Delta\ell_1} [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)], \qquad (46) $$

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

$$ \text{Cov}_{\text{nG}}[\hat{C}_{a_1 a_2}(\ell_1) \hat{C}_{a_3 a_4}(\ell_2)] = \frac{1}{\max(A_{(12)}A_{(34)})} \int_{\tilde{\ell}_1 \in \ell_1} \frac{d^2 \tilde{\ell}_1}{A(\ell_1)} \int_{\tilde{\ell}_2 \in \ell_2} \frac{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), \qquad (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), X} = \int \frac{d\chi}{f_K(\chi)^6} W_{a_1}(\chi) \dots W_{a_4}(\chi) \bar{T}_{(A_1 A_2)(A_3 A_4)}^{X} \left( \frac{\ell_1 + 0.5}{f_K(\chi)}, \frac{\ell_2 + 0.5}{f_K(\chi)}, \chi \right), \quad (48)$$

where

$$\bar{T}_{(A_1 A_2)(A_3 A_4)}^{X}(k_1, k_2, z(\chi)) = \begin{cases} T_{(A_1 A_2)(A_1 A_2)}^{\text{SSC}}(k_1, k_2, z(\chi)) & \text{if } X = \text{SSC}, \\ \int_0^{\pi} \frac{d\phi_{\ell}}{\pi} T_{A_1 \dots A_4}(\mathbf{k}_1, -\mathbf{k}_1, \mathbf{k}_2, -\mathbf{k}_2, z(\chi)) & \text{if } X = \text{nG}, \end{cases} \quad (49)$$

where $\phi_{\ell}$ is the angle between the two wave vectors $\mathbf{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}_G[\hat{C}_{a_1 a_2}(\ell_1) \hat{C}_{a_3 a_4}(\ell_2)] + \frac{1}{\max(A_{(12)}A_{(34)})} T_{\ell_1 \ell_2}^{(a_1 a_2)(a_3 a_4), \text{nG}} + T_{\ell_1 \ell_2}^{(a_1 a_2)(a_3 a_4), \text{SSC}} \quad (50)$$

at each multipole $\ell$. For an $\ell$-band averaged version, we use Equations (44) and (47). This is only important if a pure $\ell$-space covariance is required, since the averaging over scales is taken care of when mapping the theoretical $\ell$-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(\ell)$, 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 $\ell$ modes become intermingled, and the partial sky $C_\ell$, often referred to as pseudo-$C_\ell$, are derived from the full-sky $C_\ell$ 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_\ell$ 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(\chi)$ 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 $\mathcal{O}$ between two tracers $a_1$ and $a_2$ as

$$\mathcal{O}(L) = \int \frac{\ell d\ell}{2\pi} \mathbf{W}_L(\ell) \mathrm{vec}(\mathbf{C}_{a_1 a_2})(\ell). \quad (51)$$

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

$$\mathbf{C}_{a_1 a_2}(\ell_1) \delta_{\ell_1 \ell_2}^{\mathbf{K}} \delta_{m_1 m_2}^{\mathbf{K}} = \langle \mathbf{a}_{1, \ell_1 m_1} \mathbf{a}_{2, \ell_2 m_2}^{\dagger} \rangle, \quad (52)$$

on the full sky for a homogeneous and isotropic random field. The CS submatrix of $\mathbf{W}_L(\ell)$ 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(\ell)$ of all probes (for a standard 3 × 2 analysis) in a vector by using the notation of Equation 35 and 36

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

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

$$[\mathrm{vec}(\mathcal{O})](L) := (\mathcal{O}_{gg}, \mathcal{O}_{gm}, \mathcal{O}_{mm,EE}, \mathcal{O}_{mm,BB})^{\mathrm{T}}(L). \quad (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 $\mathbf{R}_L(\theta)$

$$ \mathrm{vec}(\mathbf{O})(L) = \int \theta d\theta \mathbf{R}_L(\theta) \begin{pmatrix} w(\theta) \\ \gamma_t(\theta) \\ \xi_+(\theta) \\ \xi_-(\theta) \end{pmatrix}, \qquad (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_+^{(ij)}(\theta) = \int \frac{\ell d\ell}{2\pi} [\mathcal{P}_{\epsilon_i \epsilon_j, \mathrm{E}}(\ell) + \mathcal{P}_{\epsilon_i \epsilon_j, \mathrm{B}}(\ell)] J_0(\ell\theta), \qquad (56) $$

$$ \xi_-^{(ij)}(\theta) = \int \frac{\ell d\ell}{2\pi} [\mathcal{P}_{\epsilon_i \epsilon_j, \mathrm{E}}(\ell) - \mathcal{P}_{\epsilon_i \epsilon_j, \mathrm{B}}(\ell)] J_4(\ell\theta), \qquad (57) $$

where $J_\mu$ is a cylindrical Bessel function of the first kind of order $\mu$. 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_t^{(ij)}(\theta) = \int \frac{\ell d\ell}{2\pi} \mathcal{P}_{n_i n_j}(\ell) J_2(\ell\theta). \qquad (58) $$

The galaxy correlation function can be calculated as

$$ w^{(ij)}(\theta) = \int \frac{\ell d\ell}{2\pi} \mathcal{P}_{n_i n_j}(\ell) J_0(\ell\theta). \qquad (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_m^{a_1}$ for source bin $a_1$. Considering the multiplicative shear bias uncertainty as an additive contribution to the covariance matrix (Joachimi et al. 2021), so that

$$ \mathrm{Cov}_{\mathrm{mult}}[O_{a_1 a_2}(L_1), O_{a_3 a_4}(L_2)] = O_{a_1 a_2}(L_1) O_{a_3 a_4}(L_2) [\sigma_m^{a_1} \sigma_m^{a_3} + \sigma_m^{a_1} \sigma_m^{a_4} + \sigma_m^{a_2} \sigma_m^{a_3} + \sigma_m^{a_2} \sigma_m^{a_4}]. \qquad (60) $$

This approximation holds for $\sigma_m^{a_1} \ll 1$, for $\langle m \rangle = 0$ after correction and fully correlated m-bias corrections across tomographic bins. We note that only source samples have a $\sigma_m^{a_1}$ by definition. The assumption $\langle m \rangle = 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 $\bar{\theta}_i$ with boundaries $[\theta_{i,i}, \theta_{i,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}(\ell\bar{\theta}_i) := \frac{2}{\theta_{u,i}^2 - \theta_{l,i}^2} \int_{\theta_{l,i}}^{\theta_{u,i}} d\theta' \theta' \mathbf{J}_{\mu}(\ell\theta') = \frac{2}{(\theta_{u,i}^2 - \theta_{l,i}^2) \ell^2} \times \begin{cases} [xJ_1(x)]_{\ell\theta_{l,i}}^{\ell\theta_{u,i}} & \mu = 0, \\ [-xJ_1(x) - 2J_0(x)]_{\ell\theta_{l,i}}^{\ell\theta_{u,i}} & \mu = 2, \\ [(x - \frac{8}{x})J_1(x) - 8J_2(x)]_{\ell\theta_{l,i}}^{\ell\theta_{u,i}} & \mu = 4. \end{cases} \qquad (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):

$$ \{w, \langle \gamma_t \rangle, \xi_+, \xi_-\} \leftrightarrow \{\Xi_0, \Xi_2, \Xi_0, \Xi_4\}. $$

We note that $w$ and $\xi_+$ 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 $\ell$-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,sva}} [\Xi_{\mu}^{(ij)}(\bar{\theta}_p); \Xi_{\nu}^{(mn)}(\bar{\theta}_q)] = \frac{1}{2\pi \max(A_{(ij)}A_{(mn)})} \int d\ell \ell \mathcal{K}_{\mu}(\ell\bar{\theta}_p) \mathcal{K}_{\nu}(\ell\bar{\theta}_q) [\mathcal{P}_{im}(\ell)\mathcal{P}_{jn}(\ell) + \mathcal{P}_{in}(\ell)\mathcal{P}_{jm}(\ell)], \quad (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_{\text{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 $\bar{\theta}_p$:

$$ N_{\text{pair, ideal}}^{(ij)}(\bar{\theta}_p) = \pi(\theta_{u,p}^2 - \theta_{l,p}^2) \bar{n}_i \bar{n}_j A_{(ij)}. \quad (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 $\propto \theta^2$ 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_{\text{pair}}^{(ij)}(\bar{\theta}_p)$. Altogether,

$$ \mathrm{Cov}_{\mathrm{G,sn}} [\Xi_{\mu}^{(ij)}(\bar{\theta}_p); \Xi_{\nu}^{(mn)}(\bar{\theta}_q)] = \delta_{\mu\nu}^{K} \delta_{\bar{\theta}_p \bar{\theta}_q}^{K} (\delta_{im}^{K} \delta_{jn}^{K} + \delta_{in}^{K} \delta_{jm}^{K}) \frac{\mathcal{T}_{(ij)(mn)}^{\text{sn}}}{N_{\text{pair}}^{(ij)}(\bar{\theta}_p)}, \quad (64) $$

where we defined the noise levels

$$ \mathcal{T}_{(ij)(mn)}^{\text{sn}} := \begin{cases} 2\sigma_{\epsilon_1, i}^2, & \sigma_{\epsilon_1, j}^2 \quad \text{if } i, j, m, n \in \text{source,} \\ \sigma_{\epsilon_1, i}^2, & \sigma_{\epsilon_1, j}^2 \quad \text{if both } (ij) \text{ and } (mn) \text{ contain a lens and source sample each,} \\ 1, & \sigma_{\epsilon_1, j}^2 \quad \text{if } i, j, m, n \in \text{lens.} \end{cases} \quad (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_{ij}^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 $\sim 1/N^2$, 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^2$ 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,mix}} [\Xi_{\mu}^{(ij)}(\bar{\theta}_p); \Xi_{\nu}^{(mn)}(\bar{\theta}_q)] = \delta_{jn} \frac{\mathcal{T}_j^{\text{mix}}}{2\pi n_{\text{eff}}^{(j)} \max(A_{(ij)}A_{(mn)})} \int d\ell \ell \mathcal{K}_{\mu}(\ell\bar{\theta}_p) \mathcal{K}_{\nu}(\ell\bar{\theta}_q) \mathcal{P}_{im}(\ell) + 3 \text{perm.}, \quad (66) $$

The noise level of the mixed term is defined as

$$ \mathcal{T}_j^{\text{mix}} := \begin{cases} \sigma_{\epsilon_1, j}^2 & \text{if } j \in \text{source,} \\ 1 & \text{if } j \in \text{lens.} \end{cases} \quad (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/SSC}} [\Xi_{\mu}^{(ij)}(\bar{\theta}_p); \Xi_{\nu}^{(mn)}(\bar{\theta}_q)] = \frac{1}{4\pi^2} \int d\ell_1 \ell_1 \mathcal{K}_{\mu}(\ell_1 \bar{\theta}_p) \int d\ell_2 \ell_2 \mathcal{K}_{\nu}(\ell_2 \bar{\theta}_q) \\ \times \begin{cases} \frac{1}{\max(A_{(ij)}A_{(mn)})} \int_0^{\pi} \frac{d\phi_c}{\pi} T^{(ij)(mn)}(\ell_1, -\ell_1, \ell_2, -\ell_2) & \text{for nG} \\ T_{\text{SSC}}^{(ij)(mn)}(\ell_1, \ell_2) & \text{for SSC} \end{cases}. \quad (68) $$
### 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,sva}} [C_{\mu}^{(ij)}(L_1), C_{\nu}^{(mn)}(L_2)] = \frac{2\pi}{N_{L_1} N_{L_2} \max(A_{(ij)}A_{(mn)})} \int d\ell \ell' \mathcal{W}_{\mu}^{L_1}(\ell) \mathcal{W}_{\nu}^{L_2}(\ell) [\mathcal{P}_{im}(\ell)\mathcal{P}_{jn}(\ell) + \mathcal{P}_{in}(\ell)\mathcal{P}_{jm}(\ell)] \quad (69) $$

$$ \mathrm{Cov}_{\mathrm{G,mix}} [C_{\mu}^{(ij)}(L_1), C_{\nu}^{(mn)}(L_2)] = \delta_{jn} \frac{2\pi \mathcal{T}_j^{\mathrm{mix}}}{N_{L_1} N_{L_2} n_{\mathrm{eff}}^{(j)} \max(A_{(ij)}A_{(mn)})} \int d\ell \ell' \mathcal{W}_{\mu}^{L_1}(\ell) \mathcal{W}_{\nu}^{L_2}(\ell) \mathcal{P}_{im}(\ell) + 3 \, \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{cases} W_{EE}^{L}(\ell)/2 & \text{if } \mu = \epsilon\epsilon E \\ W_{EE}^{L}(\ell) & \text{if } \mu = \text{nn} \\ W_{nE}^{L}(\ell) & \text{if } \mu = n\epsilon \\ W_{EB}^{L}(\ell) & \text{if } \mu = \epsilon\epsilon B \end{cases}. \quad (70) $$

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

$$ \mathrm{Cov}_{\mathrm{G,sn}} [C_{\mu}^{(ij)}(L_1), C_{\nu}^{(mn)}(L_2)] = \frac{\pi^2 \delta_{\mu\nu}^{K}}{N_{L_1} N_{L_2}} (\delta_{im}^{K} \delta_{jn}^{K} + \delta_{in}^{K} \delta_{jm}^{K}) \int \frac{\theta^2 d\theta}{n_{\mathrm{pair}}^{(ij)}(\theta)} \quad (71) $$

$$ \times \begin{cases} 2\sigma_{\epsilon_1, i}^2 \sigma_{\epsilon_1, j}^2 (g_+^{L_1}(\theta)g_+^{L_2}(\theta) + g_-^{L_1}(\theta)g_-^{L_2}(\theta)) & \text{if } \mu = \epsilon\epsilon \\ 4\sigma_{\epsilon_1, i}^2 h^{L_1}(\theta)h^{L_2}(\theta) & \text{if } \mu = n\epsilon \\ g_+^{L_1}(\theta)g_+^{L_2}(\theta) & \text{if } \mu = \text{nn} \end{cases} $$

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

$$ N_{\mathrm{pair}}^{(ij)} = \int_{\theta_1}^{\theta_2} d\theta' n_{\mathrm{pair}}^{(ij)}(\theta') \quad (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}} [C_{\mu}^{(ij)}(L_1), C_{\nu}^{(mn)}(L_2)] = \frac{1}{N_{L_1} N_{L_2} \max(A_{(ij)}A_{(mn)})} \int d\ell_1 \ell_1 \mathcal{W}_{\mu}^{L_1}(\ell_1) \int d\ell_2 \ell_2 \mathcal{W}_{\nu}^{L_2}(\ell_2) \\ \times \int_0^{\pi} \frac{d\phi_{\ell}}{\pi} T^{(ij)(mn)}(\ell_1, -\ell_1, \ell_2, -\ell_2), \quad (73) $$

$$ \mathrm{Cov}_{\mathrm{SSC}} [C_{\mu}^{(ij)}(L_1), C_{\nu}^{(mn)}(L_2)] = \frac{1}{N_{L_1} N_{L_2}} \int d\ell_1 \ell_1 \mathcal{W}_{\mu}^{L_1}(\ell_1) \int d\ell_2 \ell_2 \mathcal{W}_{\nu}^{L_2}(\ell_2) T_{\mathrm{SSC}}^{(ij)(mn)}(\ell_1, \ell_2). \quad (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,sva}} [E_a^{(ij)}, E_b^{(mn)}] = \frac{1}{2\pi \max(A_{(ij)}A_{(mn)})} \int d\ell \ell' W_a(\ell) W_b(\ell) [\mathcal{P}_{im}(\ell)\mathcal{P}_{jn}(\ell) + \mathcal{P}_{in}(\ell)\mathcal{P}_{jm}(\ell)] \quad (75) $$

$$ \mathrm{Cov}_{\mathrm{G,mix}} [E_a^{(ij)}, E_b^{(mn)}] = \delta_{jn} \frac{\mathcal{T}_j^{\mathrm{mix}}}{2\pi n_{\mathrm{eff}}^{(j)} \max(A_{(ij)}A_{(mn)})} \int d\ell \ell' W_a(\ell) W_b(\ell) \mathcal{P}_{im}(\ell) + 3 \, \mathrm{perm.} \quad (76) $$

$$ \mathrm{Cov}_{\mathrm{NG}} [E_a^{(ij)}, E_b^{(mn)}] = \frac{1}{2\pi \max(A_{(ij)}A_{(mn)})} \int d\ell_1 \ell_1 W_a(\ell_1) \int d\ell_2 \ell_2 W_b(\ell_2) \int_0^{\pi} \frac{d\phi_{\ell}}{\pi} T^{(ij)(mn)}(\ell_1, -\ell_1, \ell_2, -\ell_2), \quad (77) $$

$$ \mathrm{Cov}_{\mathrm{SSC}} [E_a^{(ij)}, E_b^{(mn)}] = \frac{1}{2\pi} \int d\ell_1 \ell_1 W_a(\ell_1) \int d\ell_2 \ell_2 W_b(\ell_2) T_{\mathrm{SSC}}^{(ij)(mn)}(\ell_1, \ell_2) \quad (78) $$

$$ \mathrm{Cov}_{\mathrm{G,sn}} [E_a^{(ij)}, E_b^{(mn)}] = \frac{\sigma_{\epsilon_1, i}^2 \sigma_{\epsilon_1, j}^2}{2} (\delta_{im}^{K} \delta_{jn}^{K} + \delta_{in}^{K} \delta_{jm}^{K}) \int \frac{\theta^2 d\theta}{n_{\mathrm{pair}}^{(ij)}(\theta)} [T_{+a}(\theta)T_{+b}(\theta) + T_{-a}(\theta)T_{-b}(\theta)]. \quad (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{\Phi}^{(i)}(M_{\star}^{\mu}) := \frac{1}{\Delta M_{\star}^{\mu}} \frac{1}{V_{\max}^{(i)}(M_{\star}^{\mu})} \sum_{a=1}^{N_{\text{tot}}} \Theta_{\text{H}}(z_{\text{photo},a}, z_{\text{photo}}^{(i)}) \Theta_{\text{H}}(M_{\star,a}, M_{\star}^{\mu}). \quad (80) $$

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

$$ \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), \hat{\Phi}^{(j)}(M_{\star}^{\nu})] = \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), \hat{\Phi}^{(j)}(M_{\star}^{\nu})]_{\text{sn}} + \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), \hat{\Phi}^{(j)}(M_{\star}^{\nu})]_{\text{SSC}}, \quad (81) $$

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

$$ \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), \hat{\Phi}^{(j)}(M_{\star}^{\nu})]_{\text{sn}} = \delta_{ij}^{K} \delta_{\mu\nu}^{K} \frac{\Phi^{(i)}(M_{\star}^{\mu})}{\Delta M_{\star}^{\mu} V_{\max,\mu}^{(i)}} \quad (82) $$

$$ \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), \hat{\Phi}^{(j)}(M_{\star}^{\nu})]_{\text{SSC}} = \frac{A^2 f^{(i)} f^{(j)}}{V_{\max,\mu}^{(i)} V_{\max,\nu}^{(j)}} \int d\chi \frac{p^{(i)}(\chi) p^{(j)}(\chi)}{p_{\text{tot}}^2(\chi)} f_{\chi}^2(\chi) \sigma_{\text{bg},A}^2(\chi) \tilde{\Phi}_{\mu}(\chi) \tilde{\Phi}_{\nu}(\chi), \quad (83) $$

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

$$ \tilde{\Phi}_{\mu}(\chi) := \int dM \, n_h(M, z(\chi)) \Phi(M_{\star}^{\mu} | M) b_h(M, z(\chi)). \quad (84) $$

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

$$ \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), O_{a_1 a_2}(L)] = \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), O_{a_1 a_2}(L)]_{\text{sva}} + \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), O_{a_1 a_2}(L)]_{\text{SSC}}. \quad (85) $$

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

$$ \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), O_{a_1 a_2}(L)]_{\text{sva}} = f^{(i)} \int d\ell \, W_L(\ell) \int d\chi \frac{p^{(i)}(\chi) W_{A_1}(\chi) W_{A_2}(\chi)}{p_{\text{tot}}(\chi) \chi^2} B_{\text{cmm},\mu} \left( \frac{\ell + 1/2}{\chi}, z(\chi) \right), \quad (86) $$

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

$$ B_{\text{cmm},\mu}(k, z) = \int dM \, n_h(M, z) \Phi(M_{\star}^{\mu} | M) \left(\frac{M}{\bar{\rho}_{\text{m}}}\right)^2 \bar{u}^2(k|M, z) + 2P_{\text{lin}}(k, z) \int dM \, n_h(M, z) \Phi(M_{\star}^{\mu} | M) \frac{M}{\bar{\rho}_{\text{m}}} b_h(M, z) \bar{u}(k|M, z) \int dM \, n_h(M, z) b_h(M, z) \bar{u}(k|M, z). \quad (87) $$

Finally, the SSC term is given by

$$ \text{Cov}[\hat{\Phi}^{(i)}(M_{\star}^{\mu}), O_{a_1 a_2}(L)]_{\text{SSC}} = Af^{(i)} \int d\ell \, W_L(\ell) \int d\chi \frac{p^{(i)}(\chi) W_{A_1}(\chi) W_{A_2}(\chi)}{p_{\text{tot}}(\chi) \chi^2} \frac{\partial P_{\text{mm}}}{\partial \delta_{\text{bg}}} \left( \frac{\ell + 1/2}{\chi}, z(\chi) \right) \sigma_{\text{bg},A}^2(\chi) \tilde{\Phi}_{\mu}(\chi), \quad (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<fnref n="11" /> 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<fnref n="12" />, COSMOLIKE (Krause & Eifler 2017), or PYSSC

<footnote n="11">¹¹ Available on https://github.com/rreischke/OneCovariance</footnote>

<footnote n="12">¹² https://github.com/LSSTDESC/TJPCov</footnote>
(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_\ell$, 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_\ell$. 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<fnref n="13" />.

(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{\Omega_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 4 deg² overlap with existing deep spectroscopic surveys. This data is complemented by

<footnote n="13">¹³ https://onecovariance.readthedocs.io/en/latest/index.html</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Effect of different choices for the covariance modelling for the inference of $S_8$ using the KiDS-Legacy data. Shown is the marginal maximum posterior and the corresponding 68% intervals (one $\sigma$) interval normalised to the fiducial settings assuming realistic pair counts and including feedback parametrised by $T_{\text{AGN}}$ via HMCODE2020 (Mead et al. 2020), the SSC and nG contributions, an idealised mixed term, a realistic survey mask, and the NLA model. Each blue data point replaces one of those assumptions at a time. Left: real space correlation functions. Right: bandpowers.</figcaption></figure>

an additional 23 deg$^2$ 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 ($\delta z$) $\le 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 \approx 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 \times 10^7$ galaxies for CS over an effective area of roughly 1000 deg$^2$.

## 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_8$ 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_{10} T_{\text{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_8$. 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_8$. However, ignoring the non-Gaussian and SSC contribution on its own has very little effect on the constraints on $S_8$. 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.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** KiDS-Legacy-like covariance matrix for COSEBIs. The diagonals are in the following order $i \le j$ with the tomographic bin indices $i$ and $j$ and in each little sub-block the order of the COSEBIs runs over $n = 1, \dots, 5$. Furthermore, we omit the B-mode signal here, since it is pure shape noise in the case of COSEBIs. The six different panels show the relative contribution of each term to the total covariance. In particular, we show the super-sample covariance (SSC), the non-Gaussian (nG), and the Gaussian (G) contribution in the upper three panels (compare Equations (2), and (3)). The lower three panels show the three components of the Gaussian contribution, the sample variance (sva), the mixed term (mix), and the shape or shot noise term (sn), as described in Equation (6).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Cosmic shear setup with the six KiDS-Legacy tomographic bins. The colour bar shows the tomographic bin index combination, *I* (i.e. for tomographic bins *i*, *j* the index is *I* = $\sum_{\alpha}^{i} \sum_{\beta=\alpha}^{j}$). Left: comparison between the ONECOVARIANCE code and CCL. The relative difference between the angular power spectrum calculation is shown as a percentage. Matter power spectra have been modelled using the Takahashi et al. (2012) version of the halo model. Right: relative difference in per cent between the diagonal elements of the covariance matrix for $\xi_{\pm}$ calculated with the KiDS-1000 covariance code (Joachimi et al. 2021) and the ONECOVARIANCE code. The covariance shown here does not contain shape noise.</figcaption></figure>
<table><caption>**Table 1.** Fiducial choice of cosmological parameters within ΛCDM used for the forward simulations to test the signal and noise modelling for KiDS-Legacy-like. The given cosmology is equivalent to $S_8 \equiv \sigma_8[(\Omega_b + \Omega_c)/0.3]^{1/2} = 0.82$. The parameters are set in accordance with the constraints from KiDS-1000 (Asgari et al. 2021b), while $\sigma_8$ is set to be the mid-value between KiDS-1000 and Planck (Planck Collaboration VI 2020).</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Hubble constant</td><td><i>h</i></td><td>0.670</td></tr><tr><td>Cold dark matter density</td><td><i>Ω</i><sub>c</sub></td><td>0.276</td></tr><tr><td>Baryonic matter density</td><td><i>Ω</i><sub>b</sub></td><td>0.050</td></tr><tr><td>Curvature density</td><td><i>Ω</i><sub>k</sub></td><td>0</td></tr><tr><td>Spectral index</td><td><i>n</i><sub>s</sub></td><td>0.960</td></tr><tr><td>Matter fluctuation amplitude at 8 Mpc <i>h</i><sup>-1</sup></td><td><i>σ</i><sub>8</sub></td><td>0.786</td></tr><tr><td>Baryon feedback amplitude</td><td><i>A</i><sub>bary</sub></td><td>3.1</td></tr><tr><td>Sum of the neutrino masses</td><td>Σ <i>m</i><sub>ν</sub></td><td>0.06 eV</td></tr></tbody></table>

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_8$ 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_8$ 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 $\xi_{\pm}$. 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_8$ 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_8$ 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 $\xi_{\pm}$ 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, \dots, 5$ while the tomographic bin combination $i \le 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 \le 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
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Relative difference between the signal measured in the 4224 *Egretta* mocks (realistic mask and depth variations) with the signal prediction of ONECOVARIANCE code. The colour bar indicates the different unique tomographic bin combinations, the same as in Figure 3. The dashed lines show $\xi^-$ and solid lines $\xi^+$. The grey band indicates a five per cent relative difference. The different plots show varying settings in the ONECOVARIANCE code. In particular, we distinguish whether the averaging over the $\theta$ bin (see Equation (61)) is carried out and if the pixel window due to the finite resolution of the HEALPIX map is taken into account, i.e. damping power on small scales ($N_{side}$ correction).</figcaption></figure>

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^{-3}$ 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 $\xi_{\pm}$ 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<fnref n="14" /> 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<fnref n="15" /> (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.

<footnote n="14">¹⁴ https://github.com/mwiet/kids_sbi</footnote>

<footnote n="15">¹⁵ https://glass.readthedocs.io/stable/</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** KiDS-Legacy covariance matrix for real space correlation functions. Each square represents a unique combination of tomographic bins and contains the covariance for that combination in nine angular bins between 0.5 and 300 arcmin. The diagonals are in the order $i \le j$ with the tomographic bin indices *i* and *j*. **Left:** correlation coefficient of the full covariance matrix without multiplicative shear bias uncertainty. In the lower triangle, the analytic covariance is shown, while the upper triangle shows the covariance as measured in the *Egretta* mocks (realistic mask and depth variations). **Right:** relative impact of the new mixed term (see Appendix F) on the Gaussian part covariance matrix compared to an idealised mixed term with a homogeneous pair and triplet counts. The latter were directly estimated from the KiDS-Legacy-like catalogue.</figcaption></figure>

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^{-1}$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_\ell$ 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 random (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 $w_{TG}^t$ we independently recompute the effective number density,
<figure><img src="image_7.png" /><figcaption>Fig. 6. Relative difference between the standard deviation of the *Egretta* mocks (realistic mask and depth variations) and that from the analytic covariance matrix from the ONECOVARIANCE code for the same setting as Figure 5. The dashed lines show $\xi_-$ and solid lines depict $\xi_+$. The tomographic bin combination is shown in the top left corner of each subplot.</figcaption></figure>

$n_{\text{eff}}$, and the total ellipticity dispersion $\sigma_{\epsilon}$, from the galaxy population. It is possible to linearly interpolate $n_{\text{eff}}$ and $\sigma_{\epsilon}$ as a function of $w_{\text{TG}}$ allowing the evaluations of both quantities at every pixel in the KiDS-Legacy-like footprint from the single spatial $w_{\text{TG}}$ map. We use a footprint of approximately 1000 deg$^2$ with an $N_{\text{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 $\xi_{\pm}$ 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, $\xi_{\pm}$ are measured in the north and south separately using 300 logarithmically spaced bins in $\theta \in [0.5, 300]$ arcmin. The estimated correlation functions, $\hat{\xi}_{\pm}$ are then combined as a weighted mean,

$$ \hat{\xi}_{\pm}^{(ij)}(\theta) = \frac{N_{\text{pairs, N}}^{(ij)} \hat{\xi}_{\pm, N}^{(ij)}(\theta) + N_{\text{pairs, S}}^{(ij)} \hat{\xi}_{\pm, S}^{(ij)}(\theta)}{N_{\text{pairs, N}}^{(ij)} + N_{\text{pairs, S}}^{(ij)}}, \qquad (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 $\theta$ 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_{\text{eff}}$ and $\sigma_{\epsilon}$. In particular, the intrinsic ellipticities are sampled from a Gaussian distribution with fixed $\sigma_{\epsilon}$ and zero mean for each tomographic bin. In the most complex mock, *Egretta*, the homogeneity and isotropy assumption in $n(z)$, $n_{\text{eff}}$ and $\sigma_{\epsilon}$ are dropped, and they are sampled from realisations in the $w_{\text{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 $\xi_{\pm}$ and calculated from theory via the ONECOVARIANCE code is shown in Figure 4. Dashed lines represent $\xi_{-}$ and solid lines $\xi_{+}$ 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, $\ell \le N_{\text{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 $\xi_{\pm}$, 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 $\xi_{-}$ which requires integration up to $\ell > 10^4$ to reach convergence for $\theta < \text{arcmin}$, although it is not evident for $\xi_{+}$ at the scales displayed in this figure for clarity. For $\xi_{+}$, 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 $\theta$, 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 \le 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 $\xi_{+}$, the effect of the updated mixed term reaches a maximum of 10 per cent. However, for $\xi_{-}$ and its cross-correlation with $\xi_{+}$, the effect can be significant for individual elements of the covariance matrix, exceeding 50 per cent. Nonetheless, since the signal strength for $\xi_{-}$ 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 $\xi_{+}$ and $\xi_{-}$, 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_{\text{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-
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Showcase of the consistency tests possible with the ONECOVARIANCE code. This is a cosmic shear setting as in KiDS-1000 (Asgari et al. 2021b) (i.e. five tomographic bins, eight band power bins, five modes for COSEBIs, and nine bins for $\xi_{\pm}$. The code was used to create the full covariance matrix for any pair of those three statistics. Since these are all summary statistics for the same tracer and almost over the same physical scales, the matrix is close to being singular. The full matrix shown here has a single negative eigenvalue originating from numerical noise. However, the subcovariance matrices between each pair of summary statistics are still positive-definite. We do not show the B modes for COSEBIs and bandpowers since they mainly consist of shape noise (although they are still highly correlated with the shape noise of the other statistics).</figcaption></figure>

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_{\text{sky}}$ approximation. To assess 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 $\chi^2$ similar to the inclusion of the survey geometry in the mix term. From this we conclude that the $f_{\text{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
<figure><img src="image_9.png" /><figcaption>Fig. 8. Left: redshift distribution of the target sample for five tomographic bins of KiDS-1000 as solid lines. The grey bars indicate the reference sample. Right: Pearson correlation coefficient with the lower triangle showing the predictions of the analytical prescription, Equation (91), while the upper triangle shows the results obtained by jackknife resampling from Hildebrandt et al. (2021) (compare their Figure 3).</figcaption></figure>

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_p(z) = \frac{w_{sp}(\theta(r), z)}{\Delta z \sqrt{w_{ss}(\theta(r), z)}}, \qquad (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{aligned} \mathrm{Cov}[n_{\alpha}(z_i), n_{\beta}(z_j)] &\equiv \mathrm{Cov}\left[\frac{w_{i\alpha}}{\Delta z_i \sqrt{w_{ii}}} \frac{w_{j\beta}}{\Delta z_j \sqrt{w_{jj}}}\right] \\ &\approx \frac{1}{\Delta z_i \Delta z_j \sqrt{w_{ii} w_{jj}}} \left[C_{i\alpha j\beta} - \frac{w_{i\alpha}}{2w_{ii}} C_{ii j\beta} - \frac{w_{j\beta}}{2w_{jj}} C_{i\alpha jj} + \frac{w_{i\alpha} w_{j\beta}}{4w_{ii}w_{jj}} C_{ii jj}\right]. \end{aligned} \qquad (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
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Left: auto-correlation angular power spectra of the reference sample as a function of spectroscopic redshift (colour bar). The solid lines use the Limber approximation, Equation (10), while the dashed lines make use of the full expression, Equation (8). Right: fractional difference as a percentage of the Gaussian covariance term when using Limber vs non-Limber.</figcaption></figure>

cosmological model adopted in these simulations adheres to a $\Lambda$CDM model with parameters: $\Omega_m = 0.25$, $\Omega_\Lambda = 1 - \Omega_m$, $\Omega_b = 0.044$, $\sigma_8 = 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<fnref n="16" />. 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 $\mathcal{B}_\alpha(z) = (1+z)^\alpha$, 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 \in [0.1, 1] h^{-1} \text{ 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 \neq 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 $\ell > \chi/\Delta\chi$, where $\Delta\chi$ is the width of a Gaussian redshift bin at co-moving distance $\chi$. For the present situation, this requires the non-Limber calculation up to $\ell \sim 10^3$. 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 $w_{s,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/\text{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.

<footnote n="16">¹⁶ Accessible at https://github.com/KiDS-WL/MICE2_mock</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Left: contribution to the standard deviation of the clustering redshift measurements. The solid lines correspond to the total variance, the dashed lines to the non-Gaussian contribution, and the dotted lines to the shot noise. For all tests, we ignored the SSC term since it is suppressed for galaxy clustering. Right: comparison of the analytical standard deviation (solid lines) with the simulations (symbols).</figcaption></figure>

# 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 $\xi_{\pm}$ 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

$$ (\text{spatial}_{t_1,t_2}, \text{spatial}_{t_3,t_4}, \text{mass}_{t_1,t_2}, \text{mass}_{t_3,t_4}, \text{tomo}_{i,t_1}, \text{tomo}_{j,t_2}, \text{tomo}_{i,t_3}, \text{tomo}_{j,t_4}). $$

Here, $\text{spatial}_{t_1,t_2}$ labels the spatial index of the tracer pair $t_1, t_2$ over which its two-point function is measured. This can be for example the $\theta$ bin, the band power mode, the order of the COSEBIs, or any other label of an arbitrary summary statistic passed to the code. $\text{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, $\text{tomo}_{i,t_1}$ and $\text{tomo}_{j,t_2}$ is the tomographic bin combination for the pair $(t_1, t_2)$. 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, $\xi_+$, 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 $\xi_+$ for off-diagonal elements, with this effect potentially being both positive and negative. For $\xi_-$ and its cross-correlation with $\xi_+$, the deviation can be substantial on off-diagonal elements. This is attributed to the broad support of the Fourier filter for $\xi_-$ 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 $\xi_-$ 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_l$ 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_8$, 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_8$ 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_8$ 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_8$ 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.

**Acknowledgements.** 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. 23EXPR0SS1) 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

* Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526

* Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520

* Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration 2019, MNRAS, 484, 4127

* Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514

* Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, MNRAS, 507, 5869

* Aricò, G., Angulo, R. E., Contreras, S., et al. 2021, MNRAS, 506, 4070

* Asgari, M., & Schneider, P. 2015, A&A, 578, A50

* Asgari, M., Heymans, C., Hildebrandt, H., et al. 2019, A&A, 624, A134

* Asgari, M., Friswell, I., Yoon, M., et al. 2021a, MNRAS, 501, 3003

* Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021b, A&A, 645, A104

* Asgari, M., Mead, A. J., & Heymans, C. 2023, Open J. Astrophys., 6, 39

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

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

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

* Bartelmann, M., Dombrowski, J., Konrad, S., et al. 2021, SciPost Phys., 10, 153

* Becker, M. R., Troxel, M. A., MacCrann, N., et al. 2016, Phys. Rev. D, 94, 022002

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

* Bernardeau, F., Nishimichi, T., & Taruya, A. 2014, MNRAS, 445, 1526

* Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 07, 034

* Blazek, J. A., MacCrann, N., Troxel, M. A., & Fang, X. 2019, Phys. Rev. D, 100, 103506

* Bonnett, C., Troxel, M. A., Hartley, W., et al. 2016, Phys. Rev. D, 94, 042005
Bridle, S., & King, L. 2007, New J. Phys., 9, 444

Buddendiek, A., Schneider, P., Hildebrandt, H., et al. 2016, MNRAS, 456, 3886

Cacciato, M., van den Bosch, F. C., More, S., et al. 2009, MNRAS, 394, 929

Cacciato, M., van den Bosch, F. C., More, S., Mo, H., & Yang, X. 2013, MNRAS, 430, 767

Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, JHEP, 2012, 82

Castro, P. G., Heavens, A. F., & Kitching, T. D. 2005, Phys. Rev. D, 72, 023516

Chisari, N. E., & Pontzen, A. 2019, Phys. Rev. D, 100, 023543

Chisari, N. E., Alonso, D., Krause, E., et al. 2019a, ApJS, 242, 2

Chisari, N. E., Mead, A. J., Joudaki, S., et al. 2019b, Open J. Astrophys., 2, 4

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

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

Crocce, M., Castander, F. J., Gaztañaga, E., Fosalba, P., & Carretero, J. 2015, MNRAS, 453, 1513

Dalal, R., Li, X., Nicola, A., et al. 2023, Phys. Rev. D, 108, 123519

Dark Energy Survey and Kilo-Degree Survey Collaboration (Abbott, T. M. C., et al.) 2023, Open J. Astrophys., 6, 36

Davis, C., Rozo, E., Roodman, A., et al. 2018, MNRAS, 477, 2196

de Jong, J. T. A., Verdoes Kleijn, G. A., Kuijken, K. H., & Valentijn, E. A. 2013, Exp. Astron., 35, 25

de la Bella, L. F., Tessore, N., & Bridle, S. 2021, J. Cosmol. Astropart. Phys., 2021, 001

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

Desjacques, V., Ginat, Y. B., & Reischke, R. 2021, MNRAS, 504, 5612

Dodelson, S., & Schneider, M. D. 2013, Phys. Rev. D, 88, 063537

Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64

Dvornik, A., Hoekstra, H., Kuijken, K., et al. 2018, MNRAS, 479, 1240

Dvornik, A., Heymans, C., Asgari, M., et al. 2023, A&A, 675, A189

Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32

Eifler, T., Schneider, P., & Hartlap, J. 2009, A&A, 502, 721

Elsner, F., Leistedt, B., & Peiris, H. V. 2017, MNRAS, 465, 1847

Fang, X., Eifler, T., & Krause, E. 2020, MNRAS, 497, 2699

Fortuna, M. C., Hoekstra, H., Joachimi, B., et al. 2020, MNRAS, 501, 2983

Fosalba, P., Crocce, M., Gaztañaga, E., & Castander, F. J. 2015, MNRAS, 448, 2987

Friedrich, O., Seitz, S., Eifler, T. F., & Gruen, D. 2016, MNRAS, 456, 2662

Friedrich, O., Andrade-Oliveira, F., Camacho, H., et al. 2021, MNRAS, 508, 3125

Fry, J. N. 1984, ApJ, 279, 499

Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Hall, A., & Taylor, A. 2022, Phys. Rev. D, 105, 123527

Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2020, PASJ, 72, 16

Harnois-Déraps, J., Vafaee, S., & Van Waerbeke, L. 2012, MNRAS, 426, 1262

Harnois-Déraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399

Heavens, A. F., Jimenez, R., & Lahav, O. 2000, MNRAS, 317, 965

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

Heymans, C., Grocutt, E., Heavens, A., et al. 2013, MNRAS, 432, 2433

Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140

Hikage, C., Takada, M., Hamana, T., & Spergel, D. 2011, MNRAS, 412, 65

Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43

Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454

Hildebrandt, H., Köhlinger, F., Busch, J. L. V. D., et al. 2020, A&A, 633, A69

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Hirata, C. M., & Seljak, U. 2004, Phys. Rev. D, 70, 063526

Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2

Howlett, C., Lewis, A., Hall, A., & Challinor, A. 2012, J. Cosmol. Astropart. Phys., 2012, 027

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338

Joachimi, B. 2017, MNRAS, 466, L83

Joachimi, B., Cacciato, M., Kitching, T. D., et al. 2015, Space Sci. Rev., 193, 1

Joachimi, B., Lin, C. A., Asgari, M., et al. 2021, A&A, 646, A129

Johnston, H., Wright, A. H., Joachimi, B., et al. 2021, A&A, 648, A98

Kannawadi, A., Hoekstra, H., Miller, L., et al. 2019, A&A, 624, A92

Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126

Kitching, T. D., & Heavens, A. F. 2017, Phys. Rev. D, 95, 063522

Kogut, A., Spergel, D. N., Barnes, C., et al. 2003, ApJS, 148, 161

Kohonen, T. 1982, Biol. Cybern., 43, 59

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

Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2

Lacasa, F. 2018, A&A, 615, A1

Lacasa, F., & Grain, J. 2019, A&A, 624, A61

Lacasa, F., Pénin, A., & Aghanim, N. 2014, MNRAS, 439, 123

Lange, J. U. 2023, MNRAS, 525, 3181

Lee, H., & Dvorkin, C. 2020, JCAP, 2020, 044

Leonard, C. D., Ferreira, T., Fang, X., et al. 2023, Open J. Astrophys., 6, 8

Lesgourgues, J. 2011, ArXiv e-prints [arXiv:1104.2932]

Levin, D. 1996, J. Comput. Appl. Math., 67, 95

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

Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473

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

Li, S.-S., Kuijken, K., Hoekstra, H., et al. 2023, A&A, 670, A100

Lima, M., Cunha, C. E., Oyaizu, H., et al. 2008, MNRAS, 390, 118
Limber, D. N. 1953, ApJ, 117, 134

Linke, L., Burger, P. A., Heydenreich, S., Porth, L., & Schneider, P. 2024, A&A, 681, A33

Louca, A. J., & Sellentin, E. 2020, Open J. Astrophys., 3, 11

Loureiro, A., Whittaker, L., Spurio Mancini, A., et al. 2022, A&A, 665, A56

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

Mead, A. J., & Verde, L. 2021, MNRAS, 503, 3095

Mead, A. J., Tröster, T., Heymans, C., Van Waerbeke, L., & McCarthy, I. G. 2020, A&A, 641, A130

Mohammad, F. G., & Percival, W. J. 2022, MNRAS, 514, 1289

Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astron. Comput., 3, 23

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nicola, A., García-García, C., Alonso, D., et al. 2021, J. Cosmol. Astropart. Phys., 2021, 067

Nicola, A., Villaescusa-Navarro, F., Spergel, D. N., et al. 2022, JCAP, 2022, 046

Norberg, P., Baugh, C. M., Gaztañaga, E., & Croton, D. J. 2009, MNRAS, 396, 19

Oehl, V., & Tröster, T. 2024, ArXiv e-prints [arXiv:2407.08718]

Piras, D., Joachimi, B., & Villaescusa-Navarro, F. 2023, MNRAS, 520, 668

Planck Collaboration VI. 2020, A&A, 641, A6

Porth, L., Heydenreich, S., Burger, P., Linke, L., & Schneider, P. 2024, A&A, 689, A227

Reinecke, M., & Seljebotn, D. S. 2013, A&A, 554, A112

Reischke, R., Kiessling, A., & Schäfer, B. M. 2017, MNRAS, 465, 4016

Reischke, R., Desjacques, V., & Zaroubi, S. 2020, MNRAS, 491, 1079

Schaefer, B. M. 2009, Int. J. Mod. Phys. D, 18, 173

Schneider, P., & Kilbinger, M. 2007, A&A, 462, 841

Schneider, P., van Waerbeke, L., Kilbinger, M., & Mellier, Y. 2002, A&A, 396, 1

Schneider, P., Eifler, T., & Krause, E. 2010, A&A, 520, A116

Seljak, U. 2000, MNRAS, 318, 203

Sellentin, E., & Heavens, A. F. 2018, MNRAS, 473, 2355

Sellentin, E., Heymans, C., & Harnois-Déraps, J. 2018, MNRAS, 477, 4879

Semboloni, E., Hoekstra, H., & Schaye, J. 2013, MNRAS, 434, 148

Slepian, Z., & Eisenstein, D. J. 2015, MNRAS, 454, 4142

Smith, R. E. 2012, MNRAS, 426, 531

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Spurio Mancini, A., Taylor, P. L., Reischke, R., et al. 2018, Phys. Rev. D, 98, 103507

Spurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771

Stölzner, B., Wright, A. H., Asgari, M., et al. 2025, A&A, submitted [arXiv:2503.19442]

Sugiyama, S., Takada, M., Miyatake, H., et al. 2022, Phys. Rev. D, 105, 123537

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

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

Takada, M., & Jain, B. 2004, MNRAS, 348, 897

Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152

Takahashi, R., Nishimichi, T., Namikawa, T., et al. 2020, ApJ, 895, 113

Taylor, A., & Joachimi, B. 2014, MNRAS, 442, 2728

Taylor, A., Joachimi, B., & Kitching, T. 2013, MNRAS, 432, 1928

Tessore, N., Loureiro, A., Joachimi, B., von Wietersheim-Kramsta, M., & Jeffrey, N. 2023, Open J. Astrophys., 6, 11

Tinker, J. L., & Wetzel, A. R. 2010, AJ, 719, 88

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528

Trusov, S., Zarrouk, P., Cole, S., et al. 2024, MNRAS, 527, 9048

van den Bosch, F. C., More, S., Cacciato, M., Mo, H., & Yang, X. 2013, MNRAS, 430, 725

van den Busch, J. L., Hildebrandt, H., Wright, A. H., et al. 2020, A&A, 642, A200

van Uitert, E., Cacciato, M., Hoekstra, H., et al. 2016, MNRAS, 459, 3251

van Uitert, E., Joachimi, B., Joudaki, S., et al. 2018, MNRAS, 476, 4662

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Vlah, Z., Chisari, N. E., & Schmidt, F. 2020, J. Cosmol. Astropart. Phys., 2020, 025

Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nat Rev Phys, 2, 42

von Wietersheim-Kramsta, M., Lin, K., Tessore, N., et al. 2025, A&A, 694, A223

Wright, A. H., Hildebrandt, H., Busch, J. L. V. D., & Heymans, C. 2020, A&A, 637, A100

Wright, A. H., Kuijken, K., Hildebrandt, H., et al. 2024, A&A, 686, A170

Wright, A. H., Hildebrandt, H., van den Busch, J. L., et al. 2025a, A&A, submitted [arXiv:2503.19440]

Wright, A. H., Stölzner, B., Asgari, M., et al. 2025b, A&A, submitted [arXiv:2503.19441]

Yan, Z., Wright, A. H., Elisa Chisari, N., et al. 2025, A&A, 694, A259

Yang, X., Mo, H. J., & van den Bosch, F. C. 2009, AJ, 695, 900

Zieser, B., & Merkel, P. M. 2016, MNRAS, 459, 1586

Zuntz, J., Paterno, M., Jennings, E., et al. 2015, Astron. Comput., 12, 45

<footnote n="1">¹ Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany</footnote>

<footnote n="2">² Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum, Germany</footnote>

<footnote n="3">³ School of Mathematics, Statistics and Physics, Newcastle University, Herschel Building, NE1 7RU, Newcastle-upon-Tyne, UK</footnote>

<footnote n="4">⁴ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="5">⁵ Department of Physics, Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="6">⁶ Department of Physics, Centre for Extragalactic Astronomy, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="7">⁷ E. A. Milne Centre, University of Hull, Cottingham Road, Hull, HU6 7RX, UK</footnote>

<footnote n="8">⁸ Centre of Excellence for Data Science, AI, and Modelling (DAIM), University of Hull, Cottingham Road, Hull, HU6 7RX, UK</footnote>

<footnote n="9">⁹ Center for Theoretical Physics, Polish Academy of Sciences, al. Lotników 32/46, 02-668 Warsaw, Poland</footnote>

<footnote n="10">¹⁰ Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada</footnote>
¹¹ Institute for Theoretical Physics, Utrecht University, Princetonplein 5, 3584CC Utrecht, The Netherlands

¹² Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, UK

¹³ Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Av. Complutense 40, E-28040 Madrid, Spain

¹⁴ Institute of Cosmology & Gravitation, Dennis Sciama Building, University of Portsmouth, Portsmouth, PO1 3FX, United Kingdom

¹⁵ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands

¹⁶ Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

¹⁷ Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria

¹⁸ Donostia International Physics Center, Manuel Lardizabal Ibilbidea, 4, 20018 Donostia, Gipuzkoa, Spain

¹⁹ INAF – Osservatorio Astronomico di Roma, Via Frascati 33, 00040 Monteporzio Catone (Roma), Italy

²⁰ INFN – Sezione di Roma, Piazzale Aldo Moro, 2 – c/o Dipartimento di Fisica, Edificio G. Marconi, I-00185 Roma, Italy

²¹ Institute for Particle Physics and Astrophysics, ETH Zürich Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland
<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** General flowchart of the ONECOVARIANCE code. The green boxes indicate files that are fed to the code (visit the READTHEDOCS webpage, https://onecovariance.readthedocs.io/en/latest/index.html, for a more detailed discussion). The dashed lines with arrows indicate that files or functionalities are included, but not inherited. The solid arrows instead indicate inheritance. Each blue box indicates a PYTHON module in the form of a class with the corresponding name. The red box is a C++ class wrapped into PYTHON with PYBIND to carry out some of the heavy lifting in the code. Section numbers indicate where the corresponding equations and description of the content of each module can be found in the paper.</figcaption></figure>

# 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_1$ and $k_2$ 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 $\theta$, 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, THETA SPACE, and COSEBI classes all work in the same way. They use the harmonic space expressions and map them into the desired L, $\theta$ 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<fnref n="17" /> 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}) = \sum_{\ell,m} a_{\ell m} Y_{\ell m}(\hat{n}). \qquad (B.1)$$

<footnote n="17">¹⁷ Ten because it includes B modes for lensing. Instead, for pure $C_\ell$, there are only six.</footnote>
The spherical harmonics, $Y_{\ell m}(\hat{\mathbf{n}})$, satisfy the orthogonality relation

$$ \int d\Omega_n Y_{\ell_1 m_1}(\hat{\mathbf{n}}) Y_{\ell_2 m_2}^*(\hat{\mathbf{n}}) = \delta_{\ell_1 \ell_2} \delta_{m_1 m_2}^K, \qquad (B.2) $$

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

$$ A(\mathbf{x}) = \int_{\mathbf{k}} A(\mathbf{k}) e^{-i\mathbf{k}\cdot\mathbf{x}} := \int \frac{d^3\mathbf{k}}{(2\pi)^3} A(\mathbf{k}) e^{-i\mathbf{k}\cdot\mathbf{x}}. \qquad (B.3) $$

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

$$ \langle A_1(\mathbf{k}_1, \chi_1) \dots A_1(\mathbf{k}_n, \chi_n) \rangle = (2\pi)^3 P_{A_1 \dots A_n}(\mathbf{k}_1, \chi_1, \dots, \mathbf{k}_n, \chi_n) \delta_D^{(3)}(\mathbf{k}_1 \dots n), \qquad (B.4) $$

where $\mathbf{k}_{1 \dots n} := \mathbf{k}_1 + \dots + \mathbf{k}_n$. Using the Rayleigh expansion, the exponential can be written as

$$ e^{i\mathbf{k}\cdot\mathbf{x}} = 4\pi \sum_{\ell, m} i^\ell j_\ell(kx) Y_{\ell m}(\hat{\mathbf{k}}) Y_{\ell m}^*(\hat{\mathbf{n}}), \qquad (B.5) $$

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

## Appendix C: Polyspectra on the sphere

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

$$ \langle a_{1\ell_1 m_1} \dots a_{n\ell_n m_n} \rangle = (i)^{\ell_1 n} (2\pi)^3 \prod_{\alpha=1}^{n} \left( 4\pi \int d\chi_\alpha W_{A_\alpha}(\chi_\alpha) \int \frac{k_\alpha^2}{(2\pi)^3} \frac{dk_\alpha}{d\Omega_{k_\alpha}} j_{\ell_\alpha}(k_\alpha \chi_\alpha) Y_{\ell_\alpha m_\alpha}^*(\hat{\mathbf{k}}_\alpha) \right) \times P_{A_1 \dots A_n}(\mathbf{k}_1, \dots, \mathbf{k}_n) \delta_D^{(3)}(\mathbf{k}_{1 \dots n}), \qquad (C.1) $$

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

$$ \langle a_{1\ell_1 m_1} \dots a_{n\ell_n m_n} \rangle = (i)^{\ell_1 n} (2\pi)^3 \prod_{\alpha=1}^{n} \left( (4\pi)^2 \int d\chi_\alpha W_{A_\alpha}(\chi_\alpha) \int \frac{k_\alpha^2}{(2\pi)^3} \frac{dk_\alpha}{d\Omega_{k_\alpha}} j_{\ell_\alpha}(k_\alpha \chi_\alpha) j_{\ell_\alpha}(k_\alpha x) Y_{\ell_\alpha m_\alpha}^*(\hat{\mathbf{k}}_\alpha) \right) \times P_{A_1 \dots A_n}(\mathbf{k}_1, \dots, \mathbf{k}_n) \prod_{\beta=1}^{n} \sum_{\ell_\beta \ell_\beta'} \left( (i)^{\ell_\beta} \mathcal{H}_{\ell_1 \dots \ell_n}^{\ell_1 \dots \ell_n} \int x^2 dx \, j_{\ell_\beta}(k_\beta x) Y_{\ell_\beta \ell_\beta'}(\hat{\mathbf{k}}_\beta) \right). \qquad (C.2) $$

Here we defined the integral over spherical harmonics,

$$ \mathcal{H}_{\ell_1 \dots \ell_n}^{m_1 \dots m_n} \equiv \int d\Omega_x \prod_{\alpha=1}^{n} Y_{\ell_\alpha m_\alpha}(\hat{\mathbf{x}}), \qquad (C.3) $$

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

$$ Y_{\ell_1 m_1}(\hat{\mathbf{x}}) Y_{\ell_2 m_2}(\hat{\mathbf{x}}) = \sum_{\ell m} c_{\ell m}(\ell_1, \ell_2, m_1, m_2) Y_{\ell m}(\hat{\mathbf{x}}), \qquad (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} \qquad (C.5) $$

and the Gaunt integral $\mathcal{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 d\Omega_x Y_{\ell_1 m_1}(\hat{\mathbf{x}}) Y_{\ell_2 m_2}(\hat{\mathbf{x}}) Y_{\ell_3 m_3}(\hat{\mathbf{x}}) = \sqrt{\frac{(2\ell_1 + 1)(2\ell_2 + 1)(2\ell_3 + 1)}{4\pi}} \begin{pmatrix} \ell_1 & \ell_2 & \ell_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_1 & \ell_2 & \ell_3 \\ m_1 & m_2 & m_3 \end{pmatrix}. \qquad (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}. \qquad (C.7) $$

For $n=2$, 3 the functional form of the angular polyspectrum is completely fixed by rotational symmetry since the integral over $\Omega_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_s} \sum_{\ell} \tilde{a}_{\ell} e^{-i\ell\theta}, \qquad (D.1)$$

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

$$a(\theta) = \int d^2\ell \tilde{a}_{\ell} e^{-i\ell\theta}. \qquad (D.2)$$

The corresponding angular power spectra are defined via

$$\langle a_{1\ell_1} a_{2\ell_2} \rangle = A_s \delta^{\mathrm{K}}_{\ell_1 \ell_2} \mathcal{P}_{a_1 a_2}(\ell_1), \quad \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) \qquad (D.3)$$

for the discrete and continuous cases respectively, where $\delta^\mathrm{D}(x - y)$ is the Dirac $\delta$-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(kr_s) [\mathrm{Ci}(k(1+c)r_s) - \mathrm{Ci}(kr_s)] - \frac{\sin(ckr_s)}{kr_s(1+c)} + \frac{\sin(kr_s)(\mathrm{Si}(kr_s(1+c)) - \mathrm{Si}(kr_s))}{\frac{1}{1+c} + \ln(1+c) - 1}, \qquad (E.1)$$

where $\mathrm{Si}(x)$ and $\mathrm{Ci}(x)$ are sine and cosine integrals. The scaling radius

$$r_s = \frac{r_{\mathrm{vir}}}{c} = \left( \frac{3M}{4\pi\Delta_V \bar{\rho}_m c^3} \right)^{1/3} \qquad (E.2)$$

defines the empirical concentration relation $c(M)$; here $\Delta_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}. \qquad (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 $\tilde{\xi}_+$ for a tomographic setup,

$$\mathrm{Cov}_{\mathrm{G,mix}}[\tilde{\xi}_+^{(a_2 a_2)}(\bar{\theta}_1); \tilde{\xi}_+^{(a_3 a_4)}(\bar{\theta}_2)] \equiv \frac{\sigma_{\tilde{\epsilon}_1}^2}{N_{\mathrm{pair}}^{(a_1 a_2)}(\bar{\theta}_1) N_{\mathrm{pair}}^{(a_3 a_4)}(\bar{\theta}_2)} \qquad (F.1)$$

$$\times \left[ \sum_{i,j,k}^{N_{a_1}, N_{a_2}, N_{a_3}} w_i^2 w_j w_k \xi_+^{(a_2 a_3)}(\theta_{jk}) \Delta_{\theta_1}(\theta_{ij}) \Delta_{\theta_2}(\theta_{ik}) \delta_{a_1 a_4}^{\mathcal{K}} + 3 \mathrm{perm.} \right], \qquad (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_i$th tomographic redshift bin and $\Delta_{\theta_1}(\theta_{ij})$ denotes the angular bin selection function which is unity iff $\theta_{ij} \equiv |\theta_i - \theta_j| \in \theta$.

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{aligned} M_{++}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2, \phi) &\equiv \sum_{i,j,k}^{N_{a_1}, N_{a_2}, N_{a_3}} w_i^2 w_j w_k \Delta_{\theta_1}(ij) \Delta_{\theta_2}(ik) \delta_{\phi_{jk}}^{\mathcal{K}} \\ &= \frac{1}{2\pi} \sum_{n=-\infty}^{\infty} M_{++,n}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2) e^{-in\phi}, \end{aligned} \qquad (F.3)$$
where $\varphi_{ji}$ denotes the polar angle of $\theta_{ij}$, $\phi \equiv \varphi_{ki} - \varphi_{ji}$ and the multipole components can be calculated as

$$ 
\begin{aligned}
M_{++,n}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2) &= \sum_{i}^{N_{a_1}} w_i^2 c_{i,n}^{(a_1 a_2)}(\bar{\theta}_1) (c_{i,n}^{(a_1 a_3)}(\bar{\theta}_2))^* ; \\
c_{i,n}^{(a_1 a_2)}(\bar{\theta}_1) &= \sum_{j=1}^{N_{a_2}} w_j e^{i n \varphi_{ji}} \Delta_{\theta_1}(\theta_{ij}) .
\end{aligned} 
\quad (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)}(\bar{\theta}_1, \bar{\theta}_2) \equiv \int_0^{2\pi} d\phi M_{++}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2, \phi) \xi_+^{(a_2 a_3)}(\bar{\theta}_3), \quad (F.5) $$

where we evaluate the shear correlation function at $\bar{\theta}_3 = \sqrt{\bar{\theta}_1^2 + \bar{\theta}_2^2 - 2\bar{\theta}_1\bar{\theta}_2\cos(\phi)}$. Plugging everything together we arrive at an approximation of Equation (F.1) in which $\xi_+$ is not evaluated for every triplet, but instead the mean distance across the bin is taken<fnref n="18" />,

$$ \mathrm{Cov}_{\mathrm{G,mix}}[\hat{\xi}_{+}^{(a_2 a_2)}(\bar{\theta}_1); \hat{\xi}_{+}^{(a_3 a_4)}(\bar{\theta}_2)] \approx \frac{\sigma_{\epsilon_1}^2}{N_{\mathrm{pair}}^{(a_1 a_2)}(\bar{\theta}_1) N_{\mathrm{pair}}^{(a_3 a_4)}(\bar{\theta}_2)} [T_{++}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2) \delta_{a_1 a_4}^{\mathcal{K}} + 3 \mathrm{perm.}], \quad (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,mix}}[\hat{\xi}_{-}^{(a_2 a_2)}(\bar{\theta}_1); \hat{\xi}_{-}^{(a_3 a_4)}(\bar{\theta}_2)] \approx \frac{\sigma_{\epsilon_1}^2}{N_{\mathrm{pair}}^{(a_1 a_2)}(\bar{\theta}_1) N_{\mathrm{pair}}^{(a_3 a_4)}(\bar{\theta}_2)} [T_{--}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2) \delta_{a_1 a_4}^{\mathcal{K}} + 3 \mathrm{perm.}], \quad (F.7) $$

where the $T_{--}$ are given as

$$ T_{--}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2) \equiv \int_0^{2\pi} d\phi M_{--}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2, \phi) \xi_+^{(a_2 a_3)}(\bar{\theta}_3); $$

$$ M_{--}^{(a_1 a_2 a_3)}(\bar{\theta}_1, \bar{\theta}_2, \phi) \equiv \frac{1}{2} \sum_{i=1}^{N_{a_1}} w_i^2 \{ c_{i,4+n}^{(a_1 a_2)}(\bar{\theta}_1) [c_{i,4+n}^{(a_1 a_3)}(\bar{\theta}_2)]^* + c_{i,n-4}^{(a_1 a_2)}(\bar{\theta}_1) [c_{i,n-4}^{(a_1 a_3)}(\bar{\theta}_2)]^* \}. \quad (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

$$ C_{(a_1 a_2)}(L) := \frac{1}{N_L} \int \ell d\ell S_L(\ell) \mathcal{P}_{a_1 a_2}(\ell), \quad (G.1) $$

where $S_L$ is the band power response function and $N_L$ is the normalisation

$$ N_L = \int \frac{d\ell}{\ell} S_L(\ell), \quad (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$, $\gamma_t$ and $\xi_\pm$. Thus, by using Equation (G.1), one finds (Schneider et al. 2002)

$$ C_{n_i n_j}(L) = \frac{\pi}{N_L} \int \theta d\theta T(\theta) w^{(ij)}(\theta) g_+^L(\theta), \quad (G.3) $$

$$ C_{n_i \epsilon_j}(L) = \frac{2\pi}{N_L} \int \theta d\theta T(\theta) \gamma_i^{(ij)}(\theta) h^L(\theta), \quad (G.4) $$

$$ C_{\epsilon_i \epsilon_j \mathrm{E/B}}(L) = \frac{\pi}{N_L} \int \theta d\theta T(\theta) [\xi_+^{(ij)}(\theta) g_+^L(\theta) \pm \xi_-^{(ij)}(\theta) g_-^L(\theta)], \quad (G.5) $$

<footnote n="18">¹⁸ For ensuring a higher accuracy one could also divide the $\theta_i$ into sub-intervals which are then put into the corresponding $(\bar{\theta}_1, \bar{\theta}_2)$ combination.</footnote>
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) = \int \ell \, d\ell S_{L}(\ell) J_{0/4}(\ell\theta), \qquad (G.6)$$

$$h^{L}(\theta) = \int \ell \, d\ell S_{L}(\ell) J_{2}(\ell\theta). \qquad (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{cases} \cos^2\left(\frac{\pi[x-(x_{lo}+\Delta_x/2)]}{2\Delta_x}\right) & \text{if } x_{lo} - \frac{\Delta_x}{2} \le x < x_{lo} + \frac{\Delta_x}{2} \\ 1 & \text{if } x_{lo} + \frac{\Delta_x}{2} \le x < x_{up} - \frac{\Delta_x}{2} \\ \cos^2\left(\frac{\pi[x-(x_{up}-\Delta_x/2)]}{2\Delta_x}\right) & \text{if } x_{up} - \frac{\Delta_x}{2} \le x < x_{up} + \frac{\Delta_x}{2} \\ 0 & \text{else} \end{cases}, \qquad (G.8)$$

with $x = \log \theta$ and $\Delta x$ the logarithmic bandwidth. The apodisation is such that the lower and upper limit is centred on $x_{lo} = \log \theta_{lo}$ and $x_{up} = \log \theta_{up}$, respectively.

Here we chose the band power response to be a top hat between $\ell_{lo,l}$ and $\ell_{up,l}$ leading to a normalisation $N_l = \ln(\ell_{up,l}/\ell_{lo,l})$ and a closed form for $g_{\pm}^L$ and $h^L$:

$$g_{+}^{L}(\theta) = \frac{1}{\theta^2} \left[ \theta \ell_{up,L} J_{1}(\theta \ell_{up,L}) - \theta \ell_{lo,L} J_{1}(\theta \ell_{lo,L}) \right], \qquad (G.9)$$

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

$$h^{L} = -\frac{1}{\theta^2} \left[ \theta \ell_{up,L} J_{1}(\theta \ell_{up,L}) - \theta \ell_{lo,L} J_{1}(\theta \ell_{lo,L}) + 2J_{0}(\theta \ell_{up,L}) - 2J_{0}(\theta \ell_{lo,L}) \right], \qquad (G.11)$$

where

$$\mathcal{G}_{-}(x) = \left(x - \frac{8}{x}\right) J_{1}(x) - 8J_{2}(x). \qquad (G.12)$$

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

$$C_{n,n_j}(L) = \frac{1}{N_L} \int \ell \, d\ell W_{EE}^L(\ell) \mathcal{P}_{n,n_j}(\ell), \qquad (G.13)$$

$$C_{n_i \epsilon_j}(L) = \frac{1}{N_L} \int \ell \, d\ell W_{nE}^L(\ell) \mathcal{P}_{n_i \epsilon_j}(\ell), \qquad (G.14)$$

$$C_{\epsilon_i \epsilon_j E}(L) = \frac{1}{2N_L} \int \ell \, d\ell \left[ W_{EE}^L(\ell) \mathcal{P}_{\epsilon_i \epsilon_j, E}(\ell) + W_{EB}^L(\ell) \mathcal{P}_{\epsilon_i \epsilon_j, B}(\ell) \right], \qquad (G.15)$$

$$C_{\epsilon_i \epsilon_j B}(L) = \frac{1}{2N_L} \int \ell \, d\ell \left[ W_{BE}^L(\ell) \mathcal{P}_{\epsilon_i \epsilon_j, E}(\ell) + W_{BB}^L(\ell) \mathcal{P}_{\epsilon_i \epsilon_j, B}(\ell) \right], \qquad (G.16)$$

with kernels given by

$$W_{EE}^L(\ell) = W_{BB}^L(\ell) = \int \theta \, d\theta T(\theta) \left[ J_0(\ell\theta) g_+^L(\theta) + J_4(\ell\theta) g_-^L(\theta) \right], \qquad (G.17)$$

$$W_{EB}^L(\ell) = W_{BE}^L(\ell) = \int \theta \, d\theta T(\theta) \left[ J_0(\ell\theta) g_+^L(\theta) - J_4(\ell\theta) g_-^L(\theta) \right], \qquad (G.18)$$

$$W_{nE}^L(\ell) = \int \theta \, d\theta T(\theta) J_2(\ell\theta) h^L(\theta). \qquad (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^{(ij)} = \frac{1}{2} \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta [T_{+n}(\theta) \xi_{+}^{(ij)}(\theta) + T_{-n}(\theta) \xi_{-}^{(ij)}(\theta)], \qquad (H.1)$$

$$B_n^{(ij)} = \frac{1}{2} \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta [T_{+n}(\theta) \xi_{+}^{(ij)}(\theta) - T_{-n}(\theta) \xi_{-}^{(ij)}(\theta)]. \qquad (H.2)$$
The filter functions $T_{\pm}(\theta)$ are defined for a given angular range, $\theta \in [\theta_{\min}, \theta_{\max}]$. Two families of COSEBIs have been used in the past, linear-COSEBIs and log-COSEBIs, defining whether the oscillations of $T_{\pm}(\theta)$ 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^{(ij)} = \int_0^\infty \frac{\ell \, d\ell}{2\pi} \mathcal{P}_{\epsilon_i \epsilon_j, \mathbf{E}}(\ell) W_n(\ell), \qquad (H.3)$$

$$B_n^{(ij)} = \int_0^\infty \frac{\ell \, d\ell}{2\pi} \mathcal{P}_{\epsilon_i \epsilon_j, \mathbf{B}}(\ell) W_n(\ell), \qquad (H.4)$$

where the weight functions, $W_n(\ell)$, are Hankel transforms of $T_{\pm}(\theta)$

$$W_n(\ell) = \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta \, T_{+n}(\theta) J_0(\ell\theta) = \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta \, T_{-n}(\theta) J_4(\ell\theta). \qquad (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 $[\theta_{\min}, \theta_{\max}]$. Similarly to the COSEBIs, they are defined as

$$\Psi_m^{n_1 n_2} = \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta \, U_m(\theta) w^{(ij)}(\theta) \qquad (H.6)$$

$$\Psi_m^{n_1 \epsilon_j} = \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta \, Q_m(\theta) \gamma_t^{(ij)}(\theta), \qquad (H.7)$$

so that $U_m(\theta)$ is compensated and orthogonal and $Q_m(\theta)$ given by the expression

$$Q_m(\theta) = \frac{2}{\theta^2} \int_0^\theta \theta' d\theta' U_m(\theta') - U_m(\theta). \qquad (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^{\Psi}(\ell) = \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta \, U_m(\theta) J_0(\ell\theta) = \int_{\theta_{\min}}^{\theta_{\max}} \theta \, d\theta \, Q_m(\theta) J_2(\ell\theta). \qquad (H.9)$$

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

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

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

Here $\pi$ is a top-hat filter centred around $L$ with bandwidth $\Delta L$. This includes the averaging over $\theta$ bins in the weights. For bandpowers one finds

$$\mathbf{W}_L^{\text{BP}}(\ell) = \frac{\pi}{N_L} \begin{pmatrix} 2W_{\text{EE}}^L(\ell) & 0 & 0 & 0 \\ 0 & 2W_{\text{nE}}^L(\ell) & 0 & 0 \\ 0 & 0 & W_{\text{EE}}^L(\ell) & W_{\text{EB}}^L(\ell) \\ 0 & 0 & W_{\text{EB}}^L(\ell) & W_{\text{BB}}^L(\ell) \end{pmatrix} \quad \text{and} \quad \mathbf{R}_L^{\text{BP}}(\theta) = \frac{\pi}{N_L} T(\theta) \begin{pmatrix} g_+^L(\theta) & 0 & 0 & 0 \\ 0 & 2h^L(\theta) & 0 & 0 \\ 0 & 0 & g_+^L(\theta) & g_-^L(\theta) \\ 0 & 0 & g_+^L(\theta) & -g_-^L(\theta) \end{pmatrix}, \qquad (I.2)$$

and for COSEBIs:

$$\mathbf{W}_L^{\text{COSEBI}}(\ell) = \begin{pmatrix} W_L^{\Psi}(\ell) & 0 & 0 & 0 \\ 0 & W_L^{\Psi}(\ell) & 0 & 0 \\ 0 & 0 & W_L(\ell) & 0 \\ 0 & 0 & 0 & W_L(\ell) \end{pmatrix} \quad \text{and} \quad \mathbf{R}_L^{\text{COSEBI}}(\theta) = \frac{1}{2} \begin{pmatrix} 2U_L(\theta) & 0 & 0 & 0 \\ 0 & 2Q_L(\theta) & 0 & 0 \\ 0 & 0 & T_{+L}(\theta) & T_{-L}(\theta) \\ 0 & 0 & T_{+L}(\theta) & -T_{-L}(\theta) \end{pmatrix}, \qquad (I.3)$$

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

$$\text{Cov}_{\text{G,sva}}[\mathcal{O}_{p_1 L_1}^{ij} \mathcal{O}_{p_2 L_2}^{mn}] = \frac{1}{2\pi \max(A_{(ij)} A_{(mn)})} \sum_{q_1, q_2} \int \ell d\ell (\mathbf{W}_{L_1}(\ell))_{p_1 q_1} (\mathbf{W}_{L_2}(\ell))_{p_2 q_2} \qquad (I.4)$$

$$\times [\mathcal{P}^{(im)}(\ell)\mathcal{P}^{(jn)}(\ell) + \mathcal{P}^{(in)}(\ell)\mathcal{P}^{(jm)}(\ell)]_{q_1 q_2} \qquad (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}}[\mathcal{O}_{p_1,L_1}^{ij} \mathcal{O}_{p_2,L_2}^{mn}] = \sum_q N_q^{ij} \int \frac{\theta^2 d\theta}{n_q^{ij}(\theta)} (\mathbf{R}_{L_1}(\theta))_{p_1 q} (\mathbf{R}_{L_2}(\theta))_{p_2 q} [\delta_{im}^K \delta_{jn}^K + \delta_{in}^K \delta_{jm}^K], \quad (1.6) $$

where $n_q^{ij}(\theta)$ is the differential pair-count density and $N_q^{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}}[\mathcal{O}_{p_1,L_1}^{ij} \mathcal{O}_{p_2,L_2}^{mn}] = \frac{1}{4\pi^2 \max(A_{(ij)}A_{(mn)})} \sum_{q_1,q_2} \int d\ell_1 \ell_1 \int d\ell_2 \ell_2 (\mathbf{W}_{L_1}(\ell_1))_{p_1 q_1} (\mathbf{W}_{L_2}(\ell_2))_{p_2 q_2} \\ \times \int_0^{\pi} \frac{d\phi_{\ell}}{\pi} T_{q_1 q_2}^{(ij)(mn)}(\ell_1, -\ell_1, \ell_2, -\ell_2), \quad (1.7) $$

$$ \mathrm{Cov}_{\mathrm{SSC},\mathrm{sn}}[\mathcal{O}_{p_1,L_1}^{ij} \mathcal{O}_{p_2,L_2}^{mn}] = \frac{1}{4\pi^2} \int d\ell_1 \ell_1 \int d\ell_2 \ell_2 (\mathbf{W}_{L_1}(\ell_1))_{p_1 q_1} (\mathbf{W}_{L_2}(\ell_2))_{p_2 q_2} (T_{\mathrm{SSC}}^{(ij)(mn)}(\ell_1, \ell_2))_{q_1 q_2}. \quad (1.8) $$