© The Authors 2025

1. Introduction

The photometric survey of the Euclid mission (Laureijs et al. 2011; Euclid Collaboration: Mellier et al. 2025) will infer cosmology using correlations between the observed angular positions of galaxies (angular galaxy clustering), their observed shapes (cosmic shear), and the cross-correlation between positions and shapes (galaxy–galaxy lensing). These two-point statistics are powerful probes of the late-time evolution of the Universe, both on their own and in a joint ‘3 × 2 pt’ analysis. As a result, two-point statistics have become the de facto standard observable for cosmological analysis in Stage III galaxy surveys such as the Kilo-Degree Survey (Heymans et al. 2021), the Dark Energy Survey (Abbott et al. 2022), and the Subaru Hyper Suprime-Cam Survey (More et al. 2023).

Angular correlations can be quantified and measured in a variety of ways. In real-space methods, correlations are measured in terms of real angular separation on the sky. Conversely, in harmonic-space methods, observations first undergo a spherical harmonic transform before two-point statistics are extracted. Examples of real-space methods include angular correlation functions (Peebles 1973; Schneider et al. 2002), COSEBIs (Schneider et al. 2010), and band powers (Schneider et al. 2002), while examples of harmonic-space methods include various flavours of angular power spectra. As we show below, there are exact mathematical relations to transform between real-space and harmonic-space observables. In practice, however, these transformations usually cannot be applied to measured data, so that real-space and harmonic-space methods are effectively slightly different probes of the same underlying information. For that reason, Euclid will deliver data products for all of the aforementioned methods. In what follows we describe the harmonic-space measurement; the real-space methods will be described separately (Euclid Collaboration: Kilbinger et al., in prep.).

Most current methodology for the measurement of angular power spectra for 3 × 2 pt cosmology comes from the analysis of the cosmic microwave background (CMB; e.g. Kamionkowski et al. 1997; Zaldarriaga & Seljak 1997; Wandelt et al. 2001; Hivon et al. 2002). However, the observables of the CMB are continuous temperature and polarisation fields, for which maps are created by carefully planned observations that are slightly oversampled with respect to the instrument’s beam size (Dupac & Tauber 2005). The same is not true for the observables in galaxy surveys such as Euclid: galaxy clustering is observed via the individual, discrete positions of galaxies, and cosmic shear probes the gravitational lensing fields through the ellipticities of galaxies at whatever positions they may be located.

To extract angular power spectra from a photometric galaxy survey, the typical approach is then to treat observations as if they were sampling continuous fields, much like the CMB

(Alonso et al. 2019; Nicola et al. 2021). For galaxy clustering, this requires an assumption that galaxies are discrete ‘tracers’ of an underlying galaxy density field. By making pixelated maps of observed galaxy number counts, the idea is to create a fair representation of this underlying field, up to a ‘shot noise’ contribution in each pixel. Similarly, for cosmic shear, the observed ellipticities of galaxies are considered tracers of the weak lensing signal. By averaging all the observed ellipticities in each pixel of a cosmic shear map, the assumption is that one recovers the underlying field, up to a ‘shape noise’ contribution due to the distribution of intrinsic galaxy shapes.

The approximation of having a smooth continuous map of an underlying field overlaid with noise starts to break down when observations are sparse with respect to the map resolution. For example, at the angular resolution required for Euclid’s ambitious science goals (Euclid Collaboration: Mellier et al. 2025), we expect about half of the observed pixels in the resulting maps to be empty. Our approach for Euclid is therefore to consider the angular power spectra of the discrete data itself, similar to the traditional analysis of spectroscopic galaxy catalogues (Heavens & Taylor 1995; Tadros et al. 1999; Percival et al. 2004). As we show here, this is possible without assuming that observations recover an underlying continuous field. In particular, the angular power spectra from discrete data points are essentially the spherical harmonics evaluated in these points, and can hence be calculated in practice. During implementation of the Euclid pipeline following this approach, the same idea had been independently developed in two other recent works, first by Baleato Lizancos & White (2024), and subsequently by Wolz et al. (2024).

In addition to practical computation, the discrete angular power spectra offer an additional advantage on the theoretical side: not having to assume the existence of intermediary maps with resolution-dependent noise greatly simplifies theoretical predictions for the measured spectra. Apart from a number of scientific assumptions, which we track and call out explicitly, this approach allows us to obtain an exact theory for the expectations of our measurements. For angular galaxy clustering, we find that the additive shot noise bias is not random, but is a known number, and we obtain an expected galaxy clustering signal that depends directly on the angular correlation function 𝔴(θ), as originally defined by Peebles (1973), instead of the two-point statistics of an ancillary galaxy density field. For cosmic shear, we obtain an expression for the additive shape noise bias that correctly treats the interplay between reduced shear and intrinsic galaxy shapes, as well as a novel method for removing the residual additive bias from the intrinsic variance of the cosmic shear field. In light of the stringent requirements on admissible biases in the Euclid data processing pipeline, these results allow us to validate our measurements to unprecedented levels of accuracy, which would otherwise be impossible to achieve, due to uncertainty in the expectation values.

Directly measuring angular power spectra from maps is fast, which makes it the de facto standard approach for obtaining spectra, despite the emergence of competing harmonic-space methods such as quadratic maximum likelihood (QML) estimators (Tegmark 1997; Tegmark & de Oliveira-Costa 2001; Maraio et al. 2023) or Bayesian hierarchical model (BHM) estimators (Alsing et al. 2016; Loureiro et al. 2023; Sellentin et al. 2023). Compared to discrete angular power spectra, the overall computational cost of map-based spectra is essentially a function of map resolution, and largely independent of the number of objects in the input catalogues. Map-based angular power spectra therefore remain an attractive computational option, particularly in the context of a large galaxy survey such as Euclid, where we eventually expect more than 1.5 billion galaxies to be observed.

For this reason, we investigate ways to obtain spectra from finite-resolution maps, while keeping the theoretical benefits of the discrete methodology. We can achieve this by using a formalism for spin-weighted spherical convolution with an exact convolution theorem. In principle, we are able to recover the discrete angular power spectra up to a resolution-dependent band limit, and hence apply the exact theoretical predictions to map-based spectra. To improve performance even further, we also show how this approach can be approximated using HEALPix maps (Górski et al. 2005), which do not have an exact convolution theorem, and require special handling of the additive bias.

A schema of our approach is shown in Fig. 1. The text is organised similarly. In Sect. 2 we review the theory of angular power spectra, and develop results on which we rely later. In Sect. 3 we compute the angular power spectra of discrete sets of observations. In Sect. 4 we obtain the expectations of the angular power spectra, for the cases in which the observations are generated by point processes or random fields. In Sect. 5 we show how to obtain the angular power spectra of discrete observations from the usual maps. In Sect. 6 we validate our results against simulations, and show that our methodology can be applied to Euclid’s first data release. We conclude with a brief discussion of our method in Sect. 7.

The methodology presented here, for both the discrete and map-based spectra, is implemented in a publicly available code called Heracles¹. This code is used for data processing in the 3 × 2 pt pipeline within the Euclid Science Ground Segment. However, it was designed from the ground up as a modular, adaptable, and user-friendly general-purpose utility that can be used for a multitude of probes and surveys.

2. Angular power spectra

In this section we state the key results and theorems regarding the two-point statistics of arbitrary spherical functions (i.e. functions on the sphere). The crucial point here will be that the concepts of angular power spectra and angular correlation functions are well-defined not only for the particular case of homogeneous random fields, but for any function on the sphere, such as an individual realisation of a random field.

In the following, we always deal with spin-weighted spherical functions, which sometimes have spin weight zero, and we follow the definitions of Boyle (2016). We parametrise the sphere using unit vectors, which we denote $\widehat{u}$ and $\widehat{u}\prime$. A spherical function f has spin weight s if the function value $f(\widehat{u})$ transforms under a rotation γ of the coordinate frame in $\widehat{u}$ as

$$\begin{array}{c}\hfill f(\widehat{u})\mapsto {\mathrm{e}}^{-\mathrm{i}s\gamma }f(\widehat{u}).\end{array}$$(1)

It follows that a spherical function with non-zero spin weight is necessarily complex-valued. Examples of a spin-weighted spherical functions are the global surface temperature on Earth (s = 0), wind speed and direction (s = 1), or the polarisation of the CMB (s = 2).

We generally only consider spherical functions f with spin weight s that have an expansion into spin-weighted spherical harmonics _sY_lm:

$$\begin{array}{c}\hfill f(\widehat{u})={\displaystyle \sum _{\mathit{lm}}}{f}_{\mathit{lm}}{}_s{Y}_{\mathit{lm}}(\widehat{u}).\end{array}$$(2)

Here and in the following, sums always extend over all admissible values l ≥ |s| and −l ≤ m ≤ l. The coefficients f_lm of the expansion are obtained by integration against the spherical harmonics,

$$\begin{array}{c}\hfill f_{\mathit{lm}}={\displaystyle \int f}(\widehat{u}){[}_s{Y}_{\mathit{lm}}(\widehat{u}){]}^{\ast }\mathrm{d}\widehat{u},\end{array}$$(3)

where the integral extends over the entire sphere, and an asterisk denotes complex conjugation. For s = 0 (i.e. no spin weight) the expansion is in the classical spherical harmonics Y_lm ≡ ₀Y_lm. In practice, we always have s = 0 or s = 2, but we treat s as an arbitrary integer spin weight as much as possible.

2.1. Two-point statistics

For any pair of spherical functions f and f′ with respective spin weights s and s′, where f = f′ is allowed, we can define the angular correlation C^ff′(θ) as the correlation of $f(\widehat{u})$ and $f\prime (\widehat{u}\prime )$ over all points $\widehat{u},\widehat{u}\prime$ on the sphere separated by the angle θ,

$$\begin{array}{c}\hfill C^{f{f}^{\prime }}(\theta )=\frac{1}{8{\pi }^2}{\displaystyle \int \underset{\widehat{u}·{\widehat{u}}^{\prime }=cos\theta }{\int }}{\mathrm{e}}^{\mathrm{i}s\alpha }{f}^{\ast }(\widehat{u})f^{\prime }({\widehat{u}}^{\prime }){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }^{\prime }}\mathrm{d}\widehat{u}\mathrm{d}{\widehat{u}}^{\prime }.\end{array}$$(4)

Here, if the spin weights s and s′ are non-zero, the angles α and α′, defined in Appendix C, describe a rotation of the respective coordinate frames in $\widehat{u}$ and $\widehat{u}\prime$ such that the resulting correlation is frame-independent. It is clear that the definition (4) of the angular correlation function does not require f or f′ to be random fields, or possess any kind of symmetry. If f and f′ are complex-valued, the two-point statistics are not fully characterised by Eq. (4) alone, but also by the correlations between f^* and f′, as described in Appendix D.

For any angular correlation function C^ff′, we can define an associated angular power spectrum C_l^ff′ as the coefficients of the expansion of C^ff′ into the Wigner d functions d_ss′^l (Edmonds 1960),

$$\begin{array}{c}\hfill C^{f{f}^{\prime }}(\theta )={\displaystyle \sum _l}\frac{2l+1}{4\pi }{C}_l^{f{f}^{\prime }}{d}_{s{s}^{\prime }}^l(\theta ).\end{array}$$(5)

As usual, the coefficients are obtained by projection against the basis functions,

$$\begin{array}{c}\hfill C_l^{f{f}^{\prime }}=2\pi {\displaystyle {\int }_0^{\pi }}{C}^{f{f}^{\prime }}(\theta )d_{s{s}^{\prime }}^l(\theta )sin(\theta )\mathrm{d}\theta .\end{array}$$(6)

The definition (5) of the angular power spectrum in terms of the angular correlation function makes once again no reference to random fields or symmetries.

To express the angular power spectrum C_l^ff′ directly in terms of the functions f and f′, we can replace C^ff′(θ) in the angular power spectrum (6) by its definition (4),

$$\begin{array}{c}\hfill C_l^{f{f}^{\prime }}=\frac{1}{4\pi }{\displaystyle \int \int f^{\ast }}(\widehat{u})f^{\prime }({\widehat{u}}^{\prime }){\mathrm{e}}^{\mathrm{i}s\alpha }{d}_{s{s}^{\prime }}^l(\theta ){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }^{\prime }}\mathrm{d}\widehat{u}\mathrm{d}{\widehat{u}}^{\prime },\end{array}$$(7)

where the angles α, θ, α′ still depend on $\widehat{u}$ and $\widehat{u}\prime$, but we now have two unrestricted integrals over the sphere. Inserting the spherical harmonic addition theorem

$$\begin{array}{c}\hfill {\displaystyle \sum _m}{}_s{Y}_{\mathit{lm}}(\widehat{u}){}_{s^{\prime }}{Y}_{\mathit{lm}}^{\ast }({\widehat{u}}^{\prime })=\frac{2l+1}{4\pi }{\mathrm{e}}^{\mathrm{i}s\alpha }{d}_{s{s}^{\prime }}^l(\theta ){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }^{\prime }}\end{array}$$(8)

into the integrand in definition (7), the two integrals decouple, and reduce to the spherical harmonic coefficients f_lm and $f{\prime }_{\mathit{lm}}$ given by definition (3). The angular power spectrum of spherical functions f and f′ is hence equivalently defined in terms of their spherical harmonic coefficients,

$$\begin{array}{c}\hfill C_l^{f{f}^{\prime }}=\frac{1}{2l+1}{\displaystyle \sum _m}{f}_{\mathit{lm}}^{\ast }{f}_{\mathit{lm}}^{\prime }.\end{array}$$(9)

This expression is sometimes called “the estimator of the angular power spectrum”, for reasons that are given below. However, expression (9) is in fact the true angular power spectrum of the particular spherical functions f and f′ (which, in cosmology, are a particular realisation from a random process), as C_l^ff′ contains exactly the same information as the angular correlation function (4). To fully describe the two-point statistics of complex-valued functions f and f′, we hence require both C_l^ff′ and the pseudo-spectrum C_l^f*f′ (see Appendix D).

2.2. Homogeneous random fields

The angular correlation function (4) is obtained by averaging a spherical function over all pairs of points at a given angular separation θ. There is an important class of fields where this averaging over direction does not remove information from the two-point statistics. These are the random fields that are invariant under rotations, which we call ‘homogeneous’ on the sphere².

Under a rotation R of the sphere the coefficients of the spherical harmonic expansion (2) transform as (Boyle 2016)

$$\begin{array}{c}\hfill f_{\mathit{lm}}\mapsto {\displaystyle \sum _{\mu }}{f}_{l\mu }{D}_{\mu m}^l(R),\end{array}$$(10)

where D_μm^l is the Wigner D function. The importance of homogeneous random fields on the sphere is closely related to this transformation: if f is a realisation of such a field, all of its statistical properties are by definition invariant under rotations, and both sides of transformation (10) have the same distribution. The distinction between the random field itself and its realisations is important here: the random field is invariant, but any given realisation that we may observe is a fixed non-random spherical function.

Homogeneity is a powerful tool: for example, consider the product $f_{\mathit{lm}}^{\ast }{f}_{l\prime m\prime }^{\prime }$ of modes from the spherical harmonic expansion (2) of a pair of functions f and f′. Under the rotation (10), the product transforms as

$$\begin{array}{c}\hfill f_{\mathit{lm}}^{\ast }{f}_{l^{\prime }{m}^{\prime }}^{\prime }\mapsto {\displaystyle \sum _{\mu {\mu }^{\prime }}}{f}_{l\mu }^{\ast }{f}_{l^{\prime }{\mu }^{\prime }}^{\prime }{D}_{\mu m}^{l\ast }(R)D_{{\mu }^{\prime }{m}^{\prime }}^{l^{\prime }}(R).\end{array}$$(11)

If f and f′ are realisations of jointly homogeneous random fields, both sides of transformation (11) must be equal in distribution. Taking the expectation over realisations, denoted ⟨  ⋅  ⟩, we find

$$\begin{array}{c}\hfill ⟨f_{\mathit{lm}}^{\ast }{f}_{l^{\prime }{m}^{\prime }}^{\prime }⟩={\displaystyle \sum _{\mu {\mu }^{\prime }}}⟨f_{l\mu }^{\ast }{f}_{l^{\prime }{\mu }^{\prime }}^{\prime }⟩D_{\mu m}^{l\ast }(R)D_{{\mu }^{\prime }{m}^{\prime }}^{l^{\prime }}(R).\end{array}$$(12)

Integrating out the rotation R on both sides using the orthogonality of the D functions (Edmonds 1960, Eq. 4.6.1), we recover expression (9), and thus obtain the well-known expectation

$$\begin{array}{c}\hfill ⟨f_{\mathit{lm}}^{\ast }{f}_{l^{\prime }{m}^{\prime }}^{\prime }⟩={\delta }_{l{l}^{\prime }}^{\mathrm{K}}{\delta }_{m{m}^{\prime }}^{\mathrm{K}}⟨C_l^{f{f}^{\prime }}⟩,\end{array}$$(13)

where δ^K is the Kronecker delta symbol. In other words, the modes of homogeneous random fields on the sphere are uncorrelated, unless their modes numbers coincide.

Note that the expectation (13) is sometimes used to define the angular power spectrum ⟨C_l⟩ of random fields, in which case the expression (9) is used as an estimator for ⟨C_l⟩. However, we prefer to think of the sum (9) as the actual, realised, observable angular power spectrum, and ⟨C_l⟩ as its expectation over realisations.

Having obtained the two-point expectation (13) in harmonic space, its equivalent $⟨f^{\ast }(\widehat{u})f\prime (\widehat{u}\prime )⟩$ in real-space can be obtained by computing the spherical harmonic expansion (2) of the product, substituting expectation (13), and using the complex conjugate of the spherical harmonic addition theorem (8):

$$\begin{array}{c}\hfill ⟨f^{\ast }(\widehat{u})f^{\prime }({\widehat{u}}^{\prime })⟩=\frac{2l+1}{4\pi }{\displaystyle \sum _l}⟨C_l^{f{f}^{\prime }}⟩{\mathrm{e}}^{-\mathrm{i}s\alpha }{d}_{s{s}^{\prime }}^l(\theta ){\mathrm{e}}^{\mathrm{i}{s}^{\prime }{\alpha }^{\prime }}.\end{array}$$(14)

Factoring out the exponentials, the remaining sum is precisely the expectation of the relation (5) between angular power spectrum and angular correlation function,

$$\begin{array}{c}\hfill ⟨C^{f{f}^{\prime }}(\theta )⟩=\frac{2l+1}{4\pi }{\displaystyle \sum _l}⟨C_l^{f{f}^{\prime }}⟩d_{s{s}^{\prime }}^l(\theta ),\end{array}$$(15)

and expectation (14) thus yields the expected two-point statistics in real space,

$$\begin{array}{c}\hfill ⟨f^{\ast }(\widehat{u})f^{\prime }({\widehat{u}}^{\prime })⟩={\mathrm{e}}^{-\mathrm{i}s\alpha }⟨C^{f{f}^{\prime }}(\theta )⟩{\mathrm{e}}^{\mathrm{i}{s}^{\prime }{\alpha }^{\prime }}.\end{array}$$(16)

Naturally, the inverse relations (6)–(15) holds in expectation as well,

$$\begin{array}{c}\hfill ⟨C_l^{f{f}^{\prime }}⟩=2\pi {\displaystyle {\int }_0^{\pi }}⟨C^{f{f}^{\prime }}(\theta )⟩d_{s{s}^{\prime }}^l(\theta )sin(\theta )\mathrm{d}\theta ,\end{array}$$(17)

and is consistent with expectations (16) and (14).

2.3. Mixing matrices

An important special case is a random field f that is the product of a homogeneous random field g and a non-stochastic spherical function w,

$$\begin{array}{c}\hfill f(\widehat{u})=g(\widehat{u})w(\widehat{u}).\end{array}$$(18)

We usually call w a weight function, but it is in fact arbitrary, and could in principle encode systematic effects, such as position-dependent multiplicative biases, including higher-order biases with non-zero spin weights (Kitching et al. 2021; Kitching & Deshpande 2022). The functions in the product (18) can each have an associated spin weight; if s, s₁, s₂ are the respective spin weights of f, g, w, it follows that s = s₁ + s₂ by the definition of the spin weight (1).

The angular correlation function (4) of f and a second such field f′ can be expressed in terms of g, g′ and w, w′,

$$\begin{array}{cc}\hfill C^{f{f}^{\prime }}(\theta )& =\frac{1}{8{\pi }^2}{\displaystyle \int \underset{\widehat{u}·{\widehat{u}}^{\prime }=cos\theta }{\int }}[{\mathrm{e}}^{\mathrm{i}{s}_1\alpha }{g}^{\ast }(\widehat{u})g^{\prime }({\widehat{u}}^{\prime }){\mathrm{e}}^{-\mathrm{i}{s}_1^{\prime }{\alpha }^{\prime }}]\hfill \\ \hfill & \times [{\mathrm{e}}^{\mathrm{i}{s}_2\alpha }{w}^{\ast }(\widehat{u})w^{\prime }({\widehat{u}}^{\prime }){\mathrm{e}}^{-\mathrm{i}{s}_2^{\prime }{\alpha }^{\prime }}]\mathrm{d}\widehat{u}\mathrm{d}{\widehat{u}}^{\prime }.\hfill \end{array}$$(19)

To compute the expectation of expression (19), we assume that g and g′ are independent of w and w′. The expectation can then be moved into the integral, and we recover the angular correlation function (16) of g and g′,

$$\begin{array}{c}\hfill {\mathrm{e}}^{\mathrm{i}{s}_1\alpha }⟨g^{\ast }(\widehat{u})g^{\prime }({\widehat{u}}^{\prime })⟩{\mathrm{e}}^{-\mathrm{i}{s}_1^{\prime }{\alpha }^{\prime }}=⟨C^{g{g}^{\prime }}(\theta )⟩.\end{array}$$(20)

We can factor ⟨C^gg′(θ)⟩ out of the integral, which reduces to the angular correlation function (4) of w and w′,

$$\begin{array}{c}\hfill ⟨C^{f{f}^{\prime }}(\theta )⟩=⟨C^{g{g}^{\prime }}(\theta )⟩C^{w{w}^{\prime }}(\theta ).\end{array}$$(21)

We thus find that the expected angular correlation of products of homogeneous random fields and weight functions is the product of their (expected) angular correlations.

Given expectation (21), we can also compute the expected angular power spectrum using relation (6),

$$\begin{array}{c}\hfill ⟨C_l^{f{f}^{\prime }}⟩=2\pi {\displaystyle {\int }_0^{\pi }}⟨C^{g{g}^{\prime }}(\theta )⟩C^{w{w}^{\prime }}(\theta )d_{s{s}^{\prime }}^l(\theta )sin(\theta )\mathrm{d}\theta .\end{array}$$(22)

We then expand the angular correlations ⟨C^gg′(θ)⟩ and C^ww′(θ) in the integral using relation (5). Since s = s₁ + s₂ and s′=s₁′+s₂′, we obtain Gaunt’s integral for the d functions (Edmonds 1960),

$$\begin{array}{cc}& \frac{1}{2}{\displaystyle {\int }_0^{\pi }}{d}_{s_1{s}_1^{\prime }}^{l_1}(\theta )d_{s_2{s}_2^{\prime }}^{l_2}(\theta )d_{s{s}^{\prime }}^l(\theta )sin(\theta )\mathrm{d}\theta \hfill \\ \hfill & ={(-1)}^{s-s^{\prime }}\left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1& s_2& -s\end{array}\right)\left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1^{\prime }& s_2^{\prime }& -s^{\prime }\end{array}\right),\hfill \end{array}$$(23)

where the right-hand side contains the Wigner 3j symbols. The result expresses the expected angular power spectrum ⟨C_l^ff′⟩ in terms of the angular power spectra ⟨C_l^gg′⟩ and C_l^ww′,

$$\begin{array}{cc}\hfill ⟨C_l^{f{f}^{\prime }}⟩=& {(-1)}^{s-s^{\prime }}{\displaystyle \sum _{l_1{l}_2}}\frac{(2l_1+1)(2l_2+1)}{4\pi }⟨C_{l_1}^{g{g}^{\prime }}⟩C_{l_2}^{w{w}^{\prime }}\hfill \\ & \times \left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1& s_2& -s\end{array}\right)\left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1^{\prime }& s_2^{\prime }& -s^{\prime }\end{array}\right).\hfill \end{array}$$(24)

There is hence a convolution theorem for angular power spectra and angular correlation functions, which more generally holds for expansions in Wigner d functions: the product of functions in the real-space expectation (21) corresponds to a convolution in the harmonic-space expectation (24).

In practice, we usually want to keep the weight functions fixed, and compute the expectation ⟨C_l^ff′⟩ as a function of the expected angular power spectrum ⟨C_l^gg′⟩ of the underlying random fields. In that case, the convolution (24) can be separated into a linear operator containing the sum over l₂,

$$\begin{array}{cc}\hfill M_{l{l}_1}^{w{w}^{\prime }}=& {(-1)}^{s-s^{\prime }}{\displaystyle \sum _{l_2}}\frac{(2l_1+1)(2l_2+1)}{4\pi }{C}_{l_2}^{w{w}^{\prime }}\hfill \\ & \times \left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1& s_2& -s\end{array}\right)\left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1^{\prime }& s_2^{\prime }& -s^{\prime }\end{array}\right),\hfill \end{array}$$(25)

which can subsequently be applied to any given ⟨C_l1^gg′⟩,

$$\begin{array}{c}\hfill ⟨C_l^{f{f}^{\prime }}⟩={\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{w{w}^{\prime }}⟨C_{l_1}^{g{g}^{\prime }}⟩.\end{array}$$(26)

We call the operator M^ww′ the mixing matrix of the weights w, w′ applied to the fields f, f′ and g, g′. This is slightly misleading, since expression (25) is merely a formal matrix with infinitely many rows and columns. However, in practice, it is always truncated to a finite size, and hence indeed a matrix. Note that the mixing matrix (25) not only depends on w and w′, but also on the full set of spin weights.

There is an important, non-trivial consequence of the above derivation: the mixing matrix only maps the expected angular power spectrum of a homogeneous random field to the expected angular power spectrum of its product with another function. The critical step occurs in the expectation (20), which only holds (i) in expectation and (ii) for homogeneous random fields. If either condition is not fulfilled, the mixing matrix formalism breaks down. In particular, it follows that mixing matrices for random fields cannot in general be multiplied: If the function w = w₁ w₂ is the product of spherical functions w₁ and w₂, then

$$\begin{array}{c}\hfill M^{w{w}^{\prime }}\ne M^{w_1{w}_1^{\prime }}{M}^{w_2{w}_2^{\prime }},\end{array}$$(27)

except for special cases. The reason is a lack of homogeneity: the random field w₂ g that yields the mixing matrix M^w₂w₂′ is no longer homogeneous, and the product w₁ (w₂ g) is hence not described by a second mixing matrix application. For example, consider the respective footprint of the northern and southern hemisphere. Individually, both footprints have the same angular correlation function, same angular power spectrum, and same non-vanishing mixing matrix. But since the product of the footprints is identically zero, so is their combined mixing matrix.

3. Discrete observations

Having reviewed the theory of angular power spectra, we now turn our attention towards creating the necessary spherical functions from sets of discrete observations. To this end, we consider two distinct types of observations:

• Points. The information lies in the distribution of the observed positions themselves, which have no further data attached.

• Fields. The information comes from the observed values of some underlying spherical function, which is observed in a discrete set of points.

Depending on which kind of data we wish to analyse, we must proceed in slightly different ways.

3.1. Points

We first consider the case where we observe a number of points ${\widehat{u}}_k$, k = 1, 2, …, on the sphere, as well as a set of weights w_k³. In the specific case of Euclid, this might be the observed angular positions of galaxies. We can represent the set of observed points as a sum of “point masses” using the Dirac delta function δ^D,

$$\begin{array}{c}\hfill n(\widehat{u})={\displaystyle \sum _k}{w}_k{\delta }^{\mathrm{D}}(\widehat{u}-{\widehat{u}}_k),\end{array}$$(28)

where the sum extends over the observed points. This turns the discrete observations into a function defined over the entire sphere. The spherical function n has spin weight s = 0 and is a true (weighted) number density, since the integral of the definition (28) over any given area of the sphere produces the contained (weighted) number of observed points.

The spherical harmonic expansion (2) of the observed number density n is readily obtained: inserting the function (28) into the definition (3) of the spherical harmonic coefficients, we can use the defining property of the delta function,

$$\begin{array}{c}\hfill n_{\mathit{lm}}={\displaystyle \sum _k}{w}_k{Y}_{\mathit{lm}}^{\ast }({\widehat{u}}_k).\end{array}$$(29)

The spherical harmonic coefficients of the number density n are hence simply the weighted, complex-conjugated values of the spherical harmonics in the observed points.

To compute the angular power spectrum (9) of n and a second set of points $\widehat{u}{\prime }_{k\prime }$ with weights w′_k′ and associated number density n′, where the two observed sets of points can be one and the same, it suffices to insert the sum (28) of delta functions for n and n′ into definition (7), set the spin weights to zero, and carry out the integration. The result is

$$\begin{array}{c}\hfill C_l^{n{n}^{\prime }}=\frac{1}{4\pi }{\displaystyle \sum _{k{k}^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }{P}_l(cos{\theta }_{k{k}^{\prime }}),\end{array}$$(30)

where P_l = d₀₀^l is the Legendre polynomial, and θ_kk′ is the angle between ${\widehat{u}}_k$ and $\widehat{u}{\prime }_{k\prime }$. This is the exact angular power spectrum for any two sets of points.

3.2. Fields

Next, we consider observations of a set of (complex) function values g_k, k = 1, 2, …, which are observed at points ${\widehat{u}}_k$ on the sphere, and given weights w_k. As in the case of the number density (28), we can construct a spherical function f from the discrete observations using the Dirac delta function δ^D,

$$\begin{array}{c}\hfill f(\widehat{u})={\displaystyle \sum _k}{g}_k{w}_k{\delta }^{\mathrm{D}}(\widehat{u}-{\widehat{u}}_k),\end{array}$$(31)

where the sum extends over all observed values. As before, we obtain a function which is defined over the entire sphere. The spin weight s of f is the sum of the respective spin weights s₁ and s₂ of g and w: if a rotation γ of the sphere in ${\widehat{u}}_k$ transforms g_k into e^−is₁γ g_k and w_k into e^−is₂γ w_k, the function value $f({\widehat{u}}_k)$ transforms into ${\mathrm{e}}^{-\mathrm{i}s\gamma }f({\widehat{u}}_k)$.

To compute the spherical harmonic expansion (2) of f, it once again suffices to insert the function (31) into the definition (3) of the spherical harmonic coefficient and use the defining property of the delta function,

$$\begin{array}{c}\hfill f_{\mathit{lm}}={\displaystyle \sum _k}{g}_k{w}_k{}_s{Y}_{\mathit{lm}}^{\ast }({\widehat{u}}_k).\end{array}$$(32)

The spherical harmonic coefficients are therefore the complex conjugate values of the spin-weighted spherical harmonics in the observed points, multiplied by the observed values and their weights.

To compute the angular power spectrum of f and a second, similarly defined function f′, where both functions can be one and the same, we proceed as above, inserting the sum (31) of delta functions for f and f′ into definition (7) and carrying out the integration. The resulting angular power spectrum for f and f′ is

$$\begin{array}{c}\hfill C_l^{f{f}^{\prime }}=\frac{1}{4\pi }{\displaystyle \sum _{k{k}^{\prime }}}{g}_k^{\ast }{g}_{k^{\prime }}^{\prime }{w}_k{w}_{k^{\prime }}^{\prime }{\mathrm{e}}^{\mathrm{i}s{\alpha }_{k{k}^{\prime }}}{d}_{s{s}^{\prime }}^l({\theta }_{k{k}^{\prime }}){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }_{k{k}^{\prime }}^{\prime }},\end{array}$$(33)

where the angles α_kk′, θ_kk′, α′_kk′ are defined for each pair of points ${\widehat{u}}_k,\widehat{u}{\prime }_{k\prime }$ as in definition (4). This is the exact angular power spectrum given two discrete sets of observed values on the sphere.

The discrete angular power spectrum (33) demonstrates the equivalence between harmonic and real space nicely: it is equivalent to the well-known real-space estimator (Schneider et al. 2002), transformed pair by pair to harmonic space using the transformation (6). In practice, however, the two measurements do contain slightly different information, since we cannot obtain them over all angular scales, which would be required to carry out the transformation mathematically.

3.3. Cross-correlations

Finally, we can consider the case where we wish to obtain the two-point statistics between discrete sets of measured points and measured function values. Following the preceding sections, we can construct spherical functions n and f′ using the sums (28) and (31) of delta functions, respectively. The angular power spectrum is once again obtained by inserting n and f′ into definition (7) and integrating out the delta functions,

$$\begin{array}{c}\hfill C_l^{n{f}^{\prime }}=\frac{1}{4\pi }{\displaystyle \sum _{k{k}^{\prime }}}{g}_{k^{\prime }}^{\prime }{w}_k{w}_{k^{\prime }}^{\prime }{d}_{0s^{\prime }}^l({\theta }_{k{k}^{\prime }}){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }_{k{k}^{\prime }}^{\prime }},\end{array}$$(34)

where s′ is the spin weight of f′, and the angles θ_kk′, α′_kk′ are defined as above. Naturally, the result (34) is merely the special case of the angular power spectrum (33) when setting f = n, and hence g_k ≡ 1 and s = 0.

4. Expectations

We are now able to compute expectations of the angular power spectra (30), (33), and (34) when the observations are random variates, such as the cosmological data observed by Euclid. We once again have to distinguish the cases where we observe points (e.g. galaxy positions) and fields (e.g. cosmic shear). In the first case, we have two-point statistics from observed points, which are generated by point processes on the sphere. In the second case, we have two-point statistics from observed function values, which are generated by random fields.

There is a subtle difference between point processes and random fields beyond the fact that we observe positions for one and function values for the other. It is encoded in what will be called Assumptions 1 and 6 below: to compute an expectation for point processes, we must allow the random positions to vary. This means that we require a priori information about the probability of observing a point anywhere on the sphere. For random fields, we are instead able to compute expectations conditional on the observed positions and weights.

4.1. Point processes, angular clustering

For observations generated by point processes, we compute the expectation of the angular power spectrum (30) for the observed number densities n and n′. To do so, the sum in expression (30) is split into separate sums over the set of true pairs of distinct points (denoted here by k ≢ k′, meaning ${\widehat{u}}_k$ and $\widehat{u}{\prime }_{k\prime }$ are not the same observed point) and over the set of degenerate pairs of identical points (denoted by k ≡ k′, meaning ${\widehat{u}}_k$ and $\widehat{u}{\prime }_{k\prime }$ are the same observed point),

$$\begin{array}{c}\hfill C_l^{n{n}^{\prime }}=\frac{1}{4\pi }\{{\displaystyle \sum _{k{k}^{\prime }}}+{\displaystyle \sum _{k\equiv k^{\prime }}}\}{w}_k{w}_{k^{\prime }}^{\prime }{P}_l(cos{\theta }_{k{k}^{\prime }}).\end{array}$$(35)

The second sum is sometimes empty, but not always, e.g. when computing an auto-correlation, where n and n′ describe the same observation. Since θ_kk′ = 0 for k ≡ k′, the second sum contains only P_l(1) = 1, and reduces to the total weight of degenerate pairs of points in n and n′, for which we define

$$\begin{array}{c}\hfill A^{n{n}^{\prime }}=\frac{1}{4\pi }{\displaystyle \sum _{k\equiv k^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }.\end{array}$$(36)

For an auto-correlation, A^nn′ is simply the total squared weight. Overall, we thus find that the angular power spectrum (30) can be written as

$$\begin{array}{c}\hfill C_l^{n{n}^{\prime }}=\frac{1}{4\pi }{\displaystyle \sum _{k{k}^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }{P}_l(cos{\theta }_{k{k}^{\prime }})+A^{n{n}^{\prime }},\end{array}$$(37)

where the remaining sum contains the two-point statistics from true pairs of distinct points. The term A^nn′ is an additive bias from degenerate pairs of identical points, which is often called the ‘noise bias’. However, even though A^nn′ is a stochastic quantity over realisations of the point processes, for any given realisation of points, the bias (36) is evidently a known number that we can compute exactly.

We thus subtract A^nn′ from both sides of expression (37) and compute the expectation of the bias-subtracted angular power spectrum,

$$\begin{array}{c}\hfill ⟨C_l^{n{n}^{\prime }}-A^{n{n}^{\prime }}⟩=\frac{1}{4\pi }\langle {\displaystyle \sum _{k{k}^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }{P}_l(cos{\theta }_{k{k}^{\prime }})\rangle .\end{array}$$(38)

Our goal is to express this expectation in terms of the intrinsic two-point statistics of the point process. The main difficulty lies in the fact that we may not have a complete sample of observations; for example, because we were only able to observe part of the sphere, as happens in any galaxy survey such as Euclid. In addition, there may be complicated observational effects at play, which result in some random points being missed even within the survey footprint. Any systematic removal of points from our sample affects the observed two-point statistics, and must hence be carefully taken into account.

To compute the expectation (38) with missing observations and systematic effects, we set w_k = 0 and w′_k′ = 0 for all unobserved points ${\widehat{u}}_k$ and $\widehat{u}{\prime }_{k\prime }$ in the (unknown) complete sample. We can then extend the sum in expression (38) to all points generated by the point process, both observed and unobserved, without changing its value,

$$\begin{array}{c}\hfill ⟨C_l^{n{n}^{\prime }}-A^{n{n}^{\prime }}⟩=\frac{1}{4\pi }\langle {\displaystyle \sum _{\mathrm{all}k{k}^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }{P}_l(cos{\theta }_{k{k}^{\prime }})\rangle .\end{array}$$(39)

Naturally, we have no knowledge about the unobserved points ${\widehat{u}}_k$ in the sum, but that will not be a problem for computing the expectation.

Since each realisation of the point process yields a different set of observed points, the weights w_k and w′_k′ in the expectation (39) are themselves random variables. We can use the law of total expectation to compute the expectation of w_k conditional on ${\widehat{u}}_k$, by making

Assumption 1. asmSPSHYP0 There exist functions v and v′ that describe the expected weight conditional on the observed position,

$$\begin{array}{c}\hfill ⟨w_k\mid {\widehat{u}}_k⟩=v({\widehat{u}}_k),\end{array}$$(40)

and similarly $⟨w{\prime }_{k\prime }\mid \widehat{u}{\prime }_{k\prime }⟩=v\prime (\widehat{u}{\prime }_{k\prime })$.

We call v and v′ the (weighted) ‘visibility’ of the respective observation; for unit weights, the value $v(\widehat{u})$ is a number between 0 and 1 that describes the a priori probability that a point sampled in a given position ${\widehat{u}}_k$ is observed. For general weights, the expectation is also taken over realisations of their values. In practice, estimating the visibility of a galaxy imaging survey is an open problem, and the subject of ongoing research (Johnston et al. 2021; Rodríguez-Monroy et al. 2022).

Using the visibility (40), the expectation (39) no longer depends on the exact set of observed points,

$$\begin{array}{c}\hfill ⟨C_l^{n{n}^{\prime }}-A^{n{n}^{\prime }}⟩=\frac{1}{4\pi }\langle {\displaystyle \sum _{\mathrm{all}k{k}^{\prime }}}v({\widehat{u}}_k)v^{\prime }({\widehat{u}}_{k^{\prime }}^{\prime })P_l(cos{\theta }_{k{k}^{\prime }})\rangle .\end{array}$$(41)

In fact, the expectation of the sum in expression (41) depends solely on the pairs of points ${\widehat{u}}_k,\widehat{u}{\prime }_{k\prime }$ in a given realisation. We can hence make

Assumption 2. asmSPSHYP1 All observed pairs of points have the same a priori distribution.

This is a weak assumption, since it is difficult to imagine how any specific pair of points in a realisation might be a priori distinguishable from the rest.

Under Assumption 2, all terms in the sum in expression (41) have the same expectation. If N and N′ are the respective total number of points for each point process, there are NN′ pairs of points⁴, and hence terms in the sum. Introducing functions $\overline{n}$ and $\overline{n}\prime$ with

$$\begin{array}{c}\hfill \overline{n}(\widehat{u})=\frac{N}{4\pi }v(\widehat{u}),\end{array}$$(42)

and similarly $\overline{n}\prime (\widehat{u}\prime )={(4\pi )}^{-1}N\prime v\prime (\widehat{u}\prime )$, the expectation (41) is

$$\begin{array}{c}\hfill ⟨C_l^{n{n}^{\prime }}-A^{n{n}^{\prime }}⟩=4\pi \langle \overline{n}(\widehat{u}){\overline{n}}^{\prime }({\widehat{u}}^{\prime })P_l(cos\theta )\rangle ,\end{array}$$(43)

where $\widehat{u},\widehat{u}\prime$ is a pair of random points, and θ is the angular separation between them. The functions $\overline{n}$ and $\overline{n}\prime$ can be understood as the position-dependent mean density of the observed points, taking the visibility into account. This has the conceptual advantage that we never have to define the exact sample of points to which N and $v(\widehat{u})$ refer, which would be difficult for Euclid with its complicated coverage from different ground-based surveys.

The remaining expectation on the right-hand side of expression (43) contains two random effects: one is the angular distribution of points, and the other is the random realisation of the mean densities. Here, we are only interested in the former, and we therefore make

Assumption 3. asmSPSHYP2 The expected angular power spectrum is conditional on the observed densities of points.

To see why the expectation over realisations with varying density is not very interesting, one can imagine a point process where the distribution of points is smoother or clumpier depending on the realised density. In that case, the expected two-point statistics over all densities can be arbitrarily different from the expectation conditional on the observed density, and we can extract essentially no information from our measurement. We hence want to compute an expectation that is, in some sense, close to our observation, except for the angular distribution of the points. This also agrees with intuition, since the (conditional) expectation of the observed density (28) over realisations of positions is then equal to the mean density (42),

$$\begin{array}{c}\hfill \langle n(\widehat{u})\rangle =\overline{n}(\widehat{u}).\end{array}$$(44)

However, our assumption comes with two important caveats: firstly, for galaxy clustering, the number of galaxies (as well as their weights, if given) will depend to some degree on the underlying realisation of the universe, and we are hence assuming that this correlation can be neglected. Secondly, in practice, we have no a priori knowledge about the mean density $\overline{n}$, and we must hence estimate it from the observations themselves. We check the impact of the latter point in Sect. 6.

Using Assumption 3, only the expectation over positions remains in expression (43), which is a double integral over the sphere,

$$\begin{array}{cc}& \langle \overline{n}(\widehat{u}){\overline{n}}^{\prime }(\widehat{u}\prime )P_l(cos\theta )\rangle \hfill \\ \hfill & ={\displaystyle \int \int }\overline{n}(\widehat{u}){\overline{n}}^{\prime }(\widehat{u}\prime )P_l(cos\theta )p(\widehat{u},\widehat{u}\prime )\mathrm{d}\widehat{u}\mathrm{d}\widehat{u}\prime ,\hfill \end{array}$$(45)

with $p(\widehat{u},\widehat{u}\prime )\mathrm{d}\widehat{u}\mathrm{d}\widehat{u}\prime$ the a priori probability of the point process to generate a pair of points in $\mathrm{d}\widehat{u}\mathrm{d}\widehat{u}\prime$. In the general case, this integral must be evaluated explicitly. But for the point processes in which we are interested here, we can make

Assumption 4. asmSPSHYP3 The point processes are homogeneous on the sphere, i.e. their distribution is unchanged under rotations of the sphere.

For galaxy clustering, this assumption is usually granted by the “cosmological principle”.

Under Assumption 4, the joint probability density $p(\widehat{u},\widehat{u}\prime )$ in the integral (45) depends only on the angular distance θ between the pair of points $\widehat{u}$ and $\widehat{u}\prime$. It can be written as (Landy & Szalay 1993; Peebles 1973)

$$\begin{array}{c}\hfill p(\widehat{u},{\widehat{u}}^{\prime })=\frac{1+\mathfrak{w}(\theta )}{{(4\pi )}^2},\end{array}$$(46)

where 𝔴 is the expected angular correlation function of density fluctuations in the observed point processes, which describes the clustering of points⁵.

Inserting the integral (45) and joint probability density (46) into expectation (43), we find that only $\overline{n}$ and $\overline{n}\prime$ depend explicitly on the positions $\widehat{u}$ and $\widehat{u}\prime$, while everything else depends on the angular separation θ alone,

$$\begin{array}{c}\hfill ⟨C_l^{n{n}^{\prime }}-A^{n{n}^{\prime }}⟩=\frac{1}{4\pi }{\displaystyle \int \int \overline{n}}(\widehat{u}){\overline{n}}^{\prime }({\widehat{u}}^{\prime })[1+\mathfrak{w}(\theta )]P_l(cos\theta )\mathrm{d}\widehat{u}\mathrm{d}{\widehat{u}}^{\prime }.\end{array}$$(47)

Writing the double integral over the sphere in terms of the angular separation θ recovers precisely the definition (4) of the angular correlation function $C^{\overline{n}\overline{n}\prime }$,

$$\begin{array}{c}\hfill ⟨C_l^{n{n}^{\prime }}-A^{n{n}^{\prime }}⟩=2\pi {\displaystyle {\int }_0^{\pi }{C}^{\overline{n}{\overline{n}}^{\prime }}(\theta )[1+\mathfrak{w}(\theta )]}{P}_l(cos\theta )sin(\theta )\mathrm{d}\theta .\end{array}$$(48)

Integrating the two terms of 1 + 𝔴(θ) separately, the former is the transformation (6) from $C^{\overline{n}\overline{n}\prime }(\theta )$ to $C_l^{\overline{n}\overline{n}\prime }$, while the latter is the convolution (22) of $C^{\overline{n}\overline{n}\prime }(\theta )$ and 𝔴(θ), which we can write in the form of a mixing matrix product (26). Overall, we can hence write the expectation (48) as

$$\begin{array}{c}\hfill ⟨C_l^{n{n}^{\prime }}-A^{n{n}^{\prime }}⟩=C_l^{\overline{n}{\overline{n}}^{\prime }}+{\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{\overline{n}{\overline{n}}^{\prime }}{\mathfrak{w}}_{l_1},\end{array}$$(49)

where $M_{l{l}_1}^{\overline{n}\overline{n}\prime }$ is the mixing matrix (25) due to the mean density functions $\overline{n}$ and $\overline{n}\prime$, and 𝔴_l is the angular power spectrum of the point processes, obtained from the intrinsic angular correlation function 𝔴 using relation (6).

We hence find that the expectation (49) contains the desired intrinsic two-point statistics of the point processes, in the form of 𝔴_l. However, the signal is doubly contaminated when the mean densities $\overline{n}$ and $\overline{n}\prime$ contain systematic variations, by both the angular power spectrum $C_l^{\overline{n}\overline{n}\prime }$ and by the associated mixing matrix $M_{l{l}_1}^{\overline{n}\overline{n}\prime }$. To remove these contaminations, we can directly manipulate expression (49) until it yields an estimator for 𝔴_l. While this approach is somewhat unusual in harmonic space, we show in Appendix B that it recovers well-known results from real space, such as the estimator of Landy & Szalay (1993).

In what follows, we focus instead on the more traditional approach for isolating the signal 𝔴_l in the expectation (49). We directly construct spherical functions δ and δ′ for the number density contrast of the observed points,

$$\begin{array}{c}\hfill \delta (\widehat{u})=\frac{n(\widehat{u})-\overline{n}(\widehat{u})}{{\overline{n}}_0},\end{array}$$(50)

and equivalently for $\delta \prime (\widehat{u}\prime )$, where ${\overline{n}}_0=N/(4\pi )$ denotes the total mean density over the sphere. Note that we divide here by a constant, and not by the function $\overline{n}$; for the alternative case (see Appendix B). The density contrast (50) is hence a linear combination of the spherical functions n and $\overline{n}$, and it follows that the angular power spectrum of δ and δ′ is

$$\begin{array}{c}\hfill C_l^{\delta {\delta }^{\prime }}=\frac{C_l^{n{n}^{\prime }}-C_l^{\overline{n}{n}^{\prime }}-C_l^{n{\overline{n}}^{\prime }}+C_l^{\overline{n}{\overline{n}}^{\prime }}}{{\overline{n}}_0{\overline{n}}_0^{\prime }}·\end{array}$$(51)

We therefore find that measuring the angular power spectrum of the number density contrast (50) yields a result that is equivalent to the partial-sky harmonic-space Landy–Szalay estimator (B.4). The expectation ⟨C_l^δδ′⟩ of the angular power spectrum (51) is readily computed using expressions (42), (44), and (49),

$$\begin{array}{c}\hfill ⟨C_l^{\delta {\delta }^{\prime }}-A^{\delta {\delta }^{\prime }}⟩={\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{v{v}^{\prime }}{\mathfrak{w}}_{l_1},\end{array}$$(52)

where M^vv′ is the mixing matrix for the visibilities v and v′, and $A^{\delta \delta \prime }={({\overline{n}}_0{\overline{n}}_{0\prime })}^{-1}{A}^{nn\prime }$ is the rescaled additive bias.

4.2. Random fields, cosmic shear

The second case of interest is where we observe values g_k which are the variates of an underlying random field, and use them to construct a spherical function f using definition (31). If there is a spherical function g such that the observations are the function values $g_k=g({\widehat{u}}_k)$ of g in the observed points, we can use the defining property of the delta function to factor $g(\widehat{u})$ out of the sum in definition (31),

$$\begin{array}{c}\hfill f(\widehat{u})=g(\widehat{u}){\displaystyle \sum _k}{w}_k{\delta }^{\mathrm{D}}(\widehat{u}-{\widehat{u}}_k).\end{array}$$(53)

For the remaining sum, we introduce a spherical function w, which we call the ‘weight function’ of the random field,

$$\begin{array}{c}\hfill w(\widehat{u})={\displaystyle \sum _k}{w}_k{\delta }^{\mathrm{D}}(\widehat{u}-{\widehat{u}}_k).\end{array}$$(54)

We can therefore write $f(\widehat{u})=g(\widehat{u})w(\widehat{u})$, and understand our constructed function f as the product of the function g under observation and a weight function w that encodes where and how well g has been observed. While the visibility (40) of a point process is, firstly, an expectation and, secondly, usually a relatively smooth function over the sphere, the weight function (54) of a random field consists of the given weights w_k in the observed positions ${\widehat{u}}_k$.

If the function g is the realisation of a random field, we want to use the mixing matrix formalism (26) to compute the expectation of the angular power spectrum (33) of f and a second such function f′ with $f\prime (\widehat{u}\prime )=g\prime (\widehat{u}\prime )w\prime (\widehat{u}\prime )$. To this end, we firstly require

Assumption 5. asmSPSHYP7 The functions g and g′ are realisations of jointly homogeneous random fields.

In the case of cosmic shear, this is once again a reasonable assumption by the cosmological principle. To apply the mixing matrix formalism, we further require

Assumption 6. asmSPSHYP4 The distribution of observed values g_k is conditional on the observed positions ${\widehat{u}}_k$ and weights w_k.

For cosmic shear, this assumption implies two approximations. Firstly, it ignores that the positions of galaxies are slightly correlated with their shears (source–lens clustering, Linke et al. 2024), since both positions and shears are ultimately connected to the large-scale structure of the universe. Secondly, the weights and values of shear observations are generally also slightly correlated, since more extreme galaxy shapes are harder to measure accurately, and thus given lower weights.

Under Assumptions 5 and 6, only the functions g and g′ are considered realisations of (homogeneous) random fields when computing the expectation of the angular power spectrum (33) for f = g w and f′=g′ w′, while w and w′ are considered fixed functions. We can hence use the mixing matrix formalism (26) to obtain the expected angular power spectrum of f and f′,

$$\begin{array}{c}\hfill ⟨C_l^{f{f}^{\prime }}⟩={\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{w{w}^{\prime }}⟨C_{l_1}^{g{g}^{\prime }}⟩,\end{array}$$(55)

where the mixing matrix is computed for the weight functions w and w′ of point masses following definition (54).

The situation is slightly more complicated if we observe the field g only indirectly via some intermediary observable. For cosmic shear, that observable is the galaxy ellipticity ϵ_k, which probes the cosmic shear field through the effect of weak gravitational lensing on the intrinsic galaxy shapes (e.g. Bartelmann & Schneider 2001),

$$\begin{array}{c}\hfill {ϵ}_k=\frac{{ϵ}_k^{\mathrm{i}}+g_k}{1+g_k^{\ast }{ϵ}_k^{\mathrm{i}}},\end{array}$$(56)

where ϵ_kⁱ is the intrinsic galaxy ellipticity that would have been observed without gravitational lensing. We say that the ellipticity ϵ_k traces the cosmic shear field g, because the conditional expectation of ϵ_k for a fixed value g_k and random orientations of the galaxy is (Seitz & Schneider 1997)

$$\begin{array}{c}\hfill \langle {ϵ}_k|g_k\rangle =g_k.\end{array}$$(57)

However, even though the observed ellipticity is an unbiased estimate of the cosmic shear field, the intrinsic variability of galaxy shapes leads to an increase in variance compared to the pure cosmic shear signal,

$$\begin{array}{c}\hfill \langle |{ϵ}_k{|}^2\rangle =\langle |g_k{|}^2\rangle +\langle |{ϵ}_k^{\mathrm{i}}{|}^2\frac{(1-|g_k{|}^2{)}^2}{1-|g_k{|}^2{|{ϵ}_k^{\mathrm{i}}|}^2}\rangle .\end{array}$$(58)

The second term in expectation (58) is an additional variance commonly called shape noise, and we see that the effect depends on both the variance of the intrinsic galaxy ellipticity and the one-point statistics of the cosmic shear field. In practice, there is a further contribution to shape noise due to the variance from imperfect shape measurement.

To understand the impact of noise on the expected angular power spectrum of a random field g, we make

Assumption 7. asmSPSHYP5 Observed values of the random field g have independent noise contributions.

Taken in isolation, this is not a good approximation for the shape noise of cosmic shear, since galaxies have intrinsic alignments (Joachimi et al. 2015; Kiessling et al. 2015; Kirk et al. 2015; Troxel & Ishak 2015). However, intrinsic alignments are generally absorbed into the theoretical prediction of the cosmic shear signal, so that our assumption is effectively a statement about our capability to model this effect.

Under Assumption 7, the expectation (55) of the angular power spectrum does not change its signal content, but picks up an additional variance term,

$$\begin{array}{c}\hfill ⟨C_l^{f{f}^{\prime }}⟩={\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{w{w}^{\prime }}⟨C_{l_1}^{g{g}^{\prime }}⟩+A^{f{f}^{\prime }},\end{array}$$(59)

where A^ff′ is the additive bias due to the noise variance σ_kk′² from degenerate pairs of identical objects (denoted as before by k ≡ k′)⁶,

$$\begin{array}{c}\hfill A^{f{f}^{\prime }}={\delta }_{s{s}^{\prime }}^{\mathrm{K}}\frac{1}{4\pi }{\displaystyle \sum _{k\equiv k^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }{\sigma }_{k{k}^{\prime }}^2,\end{array}$$(60)

where s and s′ are the spin weights of f and f′, respectively, as before. For random fields, the additive bias A^ff′ is therefore a true noise bias, in the sense that it is the expectation of a stochastic noise contribution, unlike the additive bias A^nn′ of the point process, which is a known number for each realisation.

For cosmic shear, we do not know, a priori, the additional variance σ_kk′² due to shape noise for each observed value g_k or g′_k′. In that situation, we can construct an estimate 𝒜^ff′ of the additive bias from the variance of the noisy observations (Nicola et al. 2021),

$$\begin{array}{c}\hfill {\mathcal{A}}^{f{f}^{\prime }}={\delta }_{s{s}^{\prime }}^{\mathrm{K}}\frac{1}{4\pi }{\displaystyle \sum _{k\equiv k^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }{ϵ}_k^{\ast }{ϵ}_{k^{\prime }}^{\prime }.\end{array}$$(61)

By expectation (58), this is a biased estimator for a non-vanishing A^ff, since it contains not only the variance due to shape noise, but the sum of intrinsic and noise variance,

$$\begin{array}{c}\hfill ⟨{\mathcal{A}}^{f{f}^{\prime }}⟩=A^{f{f}^{\prime }}+{\delta }_{s{s}^{\prime }}^{\mathrm{K}}\frac{1}{4\pi }{\displaystyle \sum _{k\equiv k^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }⟨C^{g{g}^{\prime }}(0)⟩,\end{array}$$(62)

where the expected zero-lag angular correlation ⟨C^gg′(0)⟩ is the intrinsic variance of the random fields g and g′. Subtracting 𝒜^ff′ from the measured angular power spectrum C_l^ff′ and taking the expectation using expressions (59) and (62), we obtain

$$\begin{array}{cc}\hfill ⟨C_l^{f{f}^{\prime }}-{\mathcal{A}}^{f{f}^{\prime }}⟩=& {\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{w{w}^{\prime }}⟨C_{l_1}^{g{g}^{\prime }}⟩\hfill \\ \hfill & -{\delta }_{s{s}^{\prime }}^{\mathrm{K}}\frac{1}{4\pi }{\displaystyle \sum _{k\equiv k^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }⟨C^{g{g}^{\prime }}(0)⟩.\hfill \end{array}$$(63)

Noting that the two-point statistics of g and g′ enter both terms on the right-hand side of the expectation, we use relation (5) and the properties of the Wigner d function to replace ⟨C^gg′(0)⟩ by a sum over the expected angular power spectrum,

$$\begin{array}{c}\hfill ⟨C^{g{g}^{\prime }}(0)⟩={\delta }_{s_1{s}_1^{\prime }}^{\mathrm{K}}{\displaystyle \sum _l}\frac{2l+1}{4\pi }⟨C_l^{g{g}^{\prime }}⟩,\end{array}$$(64)

with s₁ and s₁′ the respective spin weights of g and g′, as above. The expectation (63) is therefore equivalent to

$$\begin{array}{c}\hfill ⟨C_l^{f{f}^{\prime }}-{\mathcal{A}}^{f{f}^{\prime }}⟩={\displaystyle \sum _{l_1}}{\mathcal{M}}_{l{l}_1}^{w{w}^{\prime }}⟨C_{l_1}^{g{g}^{\prime }}⟩,\end{array}$$(65)

where we have introduced the reduced mixing matrix

$$\begin{array}{c}\hfill {\mathcal{M}}_{l{l}_1}^{w{w}^{\prime }}=M_{l{l}_1}^{w{w}^{\prime }}-{\delta }_{s{s}^{\prime }}^{\mathrm{K}}{\delta }_{s_1{s}_1^{\prime }}^{\mathrm{K}}\frac{2l_1+1}{4\pi }\frac{1}{4\pi }{\displaystyle \sum _{k\equiv k^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }.\end{array}$$(66)

In the expectation (65), the bias introduced by 𝒜^ff′ is thus completely absorbed into ℳ_ll1^ww′.

As it turns out, the reduced mixing matrix has a much simpler interpretation than the definition (66) suggests. Consider the angular power spectrum (33) for the pair of weight functions w and w′ with respective spin weights s₂ and s₂′. Following expression (35), we split C_l^ww′ into contributions from true pairs of distinct points (k ≢ k′) and degenerate pairs of identical points (k ≡ k′), so that we may define the known additive bias A^ww′ for the weight functions w and w′,

$$\begin{array}{c}\hfill A^{w{w}^{\prime }}={\delta }_{s_2{s}_2^{\prime }}^{\mathrm{K}}\frac{1}{4\pi }{\displaystyle \sum _{k\equiv k^{\prime }}}{w}_k{w}_{k^{\prime }}^{\prime }.\end{array}$$(67)

Since s = s₁ + s₂ and s′=s₁′+s₂, we can substitute A^ww′ for the sum in expression (66),

$$\begin{array}{c}\hfill {\mathcal{M}}_{l{l}_1}^{w{w}^{\prime }}=M_{l{l}_1}^{w{w}^{\prime }}-{\delta }_{s{s}^{\prime }}^{\mathrm{K}}{\delta }_{s_1{s}_1^{\prime }}^{\mathrm{K}}\frac{2l_1+1}{4\pi }{A}^{w{w}^{\prime }}.\end{array}$$(68)

Furthermore, we can substitute the Kronecker symbols by an identity for the Wigner 3j symbols,

$$\begin{array}{c}\hfill {\delta }_{s{s}^{\prime }}^{\mathrm{K}}{\delta }_{s_1{s}_1^{\prime }}^{\mathrm{K}}={\displaystyle \sum _{l_2\sigma }}(2l_2+1)\left(\begin{array}{ccc}l& l_1& l_2\\ -s& s_1& \sigma \end{array}\right)\left(\begin{array}{ccc}l& l_1& l_2\\ -s^{\prime }& s_1^{\prime }& \sigma \end{array}\right).\end{array}$$(69)

Using the fact that A^ww′ vanishes unless s₂ = s₂′, an equivalent way to write expression (68) is therefore

$$\begin{array}{cc}\hfill {\mathcal{M}}_{l{l}_1}^{w{w}^{\prime }}=& M_{l{l}_1}^{w{w}^{\prime }}-{\delta }_{s{s}^{\prime }}^{\mathrm{K}}{\displaystyle \sum _{l_2}}\frac{(2l_1+1)(2l_2+1)}{4\pi }{A}^{w{w}^{\prime }}\hfill \\ & \times \left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1& s_2& -s\end{array}\right)\left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1^{\prime }& s_2^{\prime }& -s^{\prime }\end{array}\right).\hfill \end{array}$$(70)

Comparing the result to the definition (25) of the mixing matrix, we indeed obtain a straightforward interpretation of the reduced mixing matrix,

$$\begin{array}{cc}\hfill {\mathcal{M}}_{l{l}_1}^{w{w}^{\prime }}=& {(-1)}^{s-s^{\prime }}{\displaystyle \sum _{l_2}}\frac{(2l_1+1)(2l_2+1)}{4\pi }(C_{l_2}^{w{w}^{\prime }}-{\delta }_{s{s}^{\prime }}^{\mathrm{K}}{A}^{w{w}^{\prime }})\hfill \\ & \times \left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1& s_2& -s\end{array}\right)\left(\begin{array}{ccc}{l}_1& l_2& l\\ s_1^{\prime }& s_2^{\prime }& -s^{\prime }\end{array}\right).\hfill \end{array}$$(71)

In other words, the reduced mixing matrix is the mixing matrix of the angular power spectrum C_l^ww′ with its additive bias A^ww′ subtracted.

To summarise, we obtain the following four key results. For noisy observations where the additive bias to the angular power spectrum is not known, which is the case for cosmic shear, we can construct the estimate (61) using the variance of the noisy observations. Subtracting the estimated additive bias from the measured angular power spectrum leads to a biased expectation (63) with respect to the mixing matrix formalism, since the estimate contains not only the additional variance due to noise, but also the intrinsic variance of the fields. However, we can return the expectation (65) to standard form by introducing a reduced mixing matrix, which implicitly removes the intrinsic variance of the random fields from the expected angular power spectrum. Finally, the reduced mixing matrix (71) is simply the mixing matrix with the additive bias of the weight functions, which is a known number, subtracted.

The nature of this correction becomes clear in real space. The unknown noise variance σ² is a delta-like contribution to the expected angular correlation function of the random fields,

$$\begin{array}{c}\hfill ⟨C^{g{g}^{\prime }}(\theta )⟩\mapsto ⟨C^{g{g}^{\prime }}(\theta )⟩+{\sigma }^2{\delta }^{\mathrm{D}}(cos\theta -cos0).\end{array}$$(72)

Subtracting the additive bias from the angular power spectrum is equivalent to subtracting the variance, which is the zero-lag correlation, from the angular correlation function. There is hence a correspondence

$$\begin{array}{c}\hfill C_l^{f{f}^{\prime }}-{\mathcal{A}}^{f{f}^{\prime }}⟺C^{f{f}^{\prime }}(\theta )-C^{f{f}^{\prime }}(0){\delta }^{\mathrm{D}}(cos\theta -cos0)\end{array}$$(73)

for the random fields, and

$$\begin{array}{c}\hfill C_l^{w{w}^{\prime }}-A^{w{w}^{\prime }}⟺C^{w{w}^{\prime }}(\theta )-C^{w{w}^{\prime }}(0){\delta }^{\mathrm{D}}(cos\theta -cos0)\end{array}$$(74)

for the weight functions. By expectation (21), the real-space equivalent of the reduced mixing matrix expectation (65) is hence

$$\begin{array}{cc}& \langle C^{f{f}^{\prime }}(\theta )-C^{f{f}^{\prime }}(0){\delta }^{\mathrm{D}}(cos\theta -cos0)\rangle \hfill \\ \hfill & =[⟨C^{g{g}^{\prime }}(\theta )⟨+{\sigma }^2{\delta }^{\mathrm{D}}(cos\theta -cos0)]\hfill \\ \hfill & \times [C^{w{w}^{\prime }}(\theta )-C^{w{w}^{\prime }}(0){\delta }^{\mathrm{D}}(cos\theta -cos0)],\hfill \end{array}$$(75)

where we can evaluate the right-hand side for all θ ≥ 0 without knowing the value of σ².

4.3. Cross-correlations, galaxy–galaxy lensing

The final case of interest is the cross-correlation of points ${\widehat{u}}_k$ generated by a point process, and observed values $g{\prime }_{k\prime }=g(\widehat{u}{\prime }_{k\prime })$ from the realisation g of a random field. The two observations define the spherical functions n and f′ as above.

For the expectation of the angular power spectrum (34) of n and f′, we again fundamentally rely on Assumption 6: the distribution of observed values g′_k′ is conditional on the observed points ${\widehat{u}}_k$ and weights w′_k′, which are held fixed. We assume that this remains true even when correlating positions and values from a single observation, in which case the observed positions are both random variates (within n) and fixed (within w′ and hence f′). For galaxy–galaxy lensing, the approximation performs worse than for cosmic shear; this is seen in Sect. 6. As in the case of intrinsic alignments, the assumption is therefore effectively a statement about our ability to model the effect of source–lens clustering in the theory part of the expectation.

To treat the point process in the expectation of the angular power spectrum (34), we proceed as before. We extend the sum over k to all points using the visibility (40), and replace w_k by $v({\widehat{u}}_k)$ under Assumption 1,

$$\begin{array}{c}\hfill ⟨C_l^{n{f}^{\prime }}⟩=\frac{1}{4\pi }\langle {\displaystyle \sum _{\mathrm{all}k}\sum _{k^{\prime }}}{g}_{k^{\prime }}^{\prime }v({\widehat{u}}_k)w_{k^{\prime }}^{\prime }{d}_{0s^{\prime }}^l({\theta }_{k{k}^{\prime }}){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }_{k{k}^{\prime }}^{\prime }}\rangle .\end{array}$$(76)

While Assumption 2 considers pairs of points, here we only have a single set, and hence make

Assumption 8. asmSPSHYP6 All random points in the cross-correlation have the same a priori distribution.

As in the case of pairs of points, this seems a weak assumption, since it is difficult to imagine how individual points might be a priori distinguishable from each other.

Under Assumption 8, the sum over k in expectation (76) reduces to N identically distributed terms. Using definition (42), we can replace the product of N and visibility v by the mean number density $\overline{n}$. Using the definition (54) of the weight function w′, we may also replace the remaining sum over k′ by an integral over $\widehat{u}\prime$,

$$\begin{array}{c}\hfill ⟨C_l^{n{f}^{\prime }}⟩={\displaystyle \int \langle }{g}^{\prime }({\widehat{u}}^{\prime })\overline{n}(\widehat{u})w^{\prime }({\widehat{u}}^{\prime })d_{0s^{\prime }}^l(\theta ){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }^{\prime }}\rangle \mathrm{d}{\widehat{u}}^{\prime },\end{array}$$(77)

where the angles θ and α′ now describe the relative orientation between the random point $\widehat{u}$ and $\widehat{u}\prime$.

Using Assumption 6, we can factor the weight $w\prime (\widehat{u}\prime )$ out of the integral in expectation (77). Furthermore, using Assumption 3, the expectation is conditional on the mean number density $\overline{n}$, and only the position $\widehat{u}$ in $\overline{n}(\widehat{u})$ is random. The remaining expectation in Eq. (77) therefore reduces to the random point $\widehat{u}$ and the random field g′. It can be computed in two steps using the law of total expectation. Using the angular correlation (16), the expectation of a homogeneous (by Assumption 5) random field g′ conditional on $\widehat{u}$ is

$$\begin{array}{c}\hfill ⟨g^{\prime }({\widehat{u}}^{\prime }){\mathrm{e}}^{-\mathrm{i}{s}^{\prime }{\alpha }^{\prime }}\mid \widehat{u}⟩=\gamma (\theta ),\end{array}$$(78)

where γ is the expected angular cross-correlation function⁷. For galaxy–galaxy lensing, the expected correlation is more commonly written in terms of a tangential component γ_t and cross-component γ_× as γ(θ) = γ_t(θ)+i γ_×(θ).

Combining expectations (78) and (77), it remains to compute the expectation over random positions $\widehat{u}$. Under Assumption 4, the point process is homogeneous, and hence

$$\begin{array}{c}\hfill ⟨C_l^{n{f}^{\prime }}⟩=\frac{1}{4\pi }{\displaystyle \int \int }\overline{n}(\widehat{u})w^{\prime }({\widehat{u}}^{\prime })\gamma (\theta )d_{0s^{\prime }}^l(\theta )\mathrm{d}\widehat{u}\mathrm{d}{\widehat{u}}^{\prime }.\end{array}$$(79)

As before, the double integral recovers the definition (4) of the angular correlation function for $\overline{n}$ and w′,

$$\begin{array}{c}\hfill ⟨C_l^{n{f}^{\prime }}⟩=2\pi {\displaystyle {\int }_0^{\pi }}{C}^{\overline{n}{w}^{\prime }}(\theta )\gamma (\theta )d_{0s^{\prime }}^l(\theta )sin(\theta )\mathrm{d}\theta ,\end{array}$$(80)

which in turn is the convolution (22) of $C^{\overline{n}w\prime }$ and γ that yields the mixing matrix (26),

$$\begin{array}{c}\hfill ⟨C_l^{n{f}^{\prime }}⟩={\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{\overline{n}{w}^{\prime }}{\gamma }_{l_1},\end{array}$$(81)

where γ_l is the angular power spectrum associated with the angular cross-correlation function γ.

Overall, we therefore obtain the intuitively clear result that the expected angular power spectrum is given by the intrinsic spectrum γ_l for point process and random field, modulated by a mixing matrix coming from the mean number density $\overline{n}$ (due to the point process) and weight function w′ (due to the random field).

In the case of two point processes, it was advantageous to correlate the density contrast δ instead of the number density n. Since the density contrast (50) is linear in n, the angular power spectrum of δ and f′ is

$$\begin{array}{c}\hfill C_l^{\delta f^{\prime }}=\frac{C_l^{n{f}^{\prime }}-C_l^{\overline{n}{f}^{\prime }}}{{\overline{n}}_0}·\end{array}$$(82)

The expectation ⟨C_l^δf′⟩ then follows immediately from the definition (42) of the mean density and expectation (81),

$$\begin{array}{c}\hfill ⟨C_l^{\delta f^{\prime }}⟩={\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{v{w}^{\prime }}{\gamma }_{l_1}-⟨C_l^{v{f}^{\prime }}⟩={\displaystyle \sum _{l_1}}{M}_{l{l}_1}^{v{w}^{\prime }}{\gamma }_{l_1},\end{array}$$(83)

where the second equality assumes that the expectation of f′ vanishes.

Whether to correlate n or δ is a well-known question for real-space estimators of galaxy–galaxy lensing (Joachimi et al. 2021). Formally, we can construct a direct estimator of the signal from either the number density n using expectation (81),

$$\begin{array}{c}\hfill {\widehat{\gamma }}_l^n={\displaystyle \sum _{l_1}}{(M^{\overline{n}{w}^{\prime }})}_{l{l}_1}^{-1}{C}_l^{n{f}^{\prime }},\end{array}$$(84)

or from the density contrast δ using expectation (83) and the definition (42) of the mean density,

$$\begin{array}{cc}\hfill {\widehat{\gamma }}_l^{\delta }& ={\displaystyle \sum _{l_1}}{(M^{v{w}^{\prime }})}_{l{l}_1}^{-1}\frac{C_l^{n{f}^{\prime }}-C_l^{\overline{n}{f}^{\prime }}}{{\overline{n}}_0}\hfill \\ \hfill & ={\displaystyle \sum _{l_1}}{(M^{\overline{n}{w}^{\prime }})}_{l{l}_1}^{-1}{C}_l^{n{f}^{\prime }}-{\displaystyle \sum _{l_1}}{(M^{\overline{n}{w}^{\prime }})}_{l{l}_1}^{-1}{C}_l^{\overline{n}{f}^{\prime }}.\hfill \end{array}$$(85)

In real space, the mean number density $\overline{n}$ corresponds to an equivalent distribution of uniform random points (‘randoms’); the inverse mixing matrix in both ${\widehat{\gamma }}_l^n$ and ${\widehat{\gamma }}_l^{\delta }$ thus corresponds to a normalisation by weighted pairs of randoms and the observed positions of the field. The difference between ${\widehat{\gamma }}_l^n$ and ${\widehat{\gamma }}_l^{\delta }$ is the second term in the estimator (85), which corresponds to correlations between randoms and field values. It was shown by Singh et al. (2017) that the estimator (85) has lower variance relative to the estimator (84), particularly on large scales, because it suppresses covariance terms that couple with the survey mask. This suppression also increases the accuracy of covariance estimation via resampling techniques, as these modify the effective survey window in the subsampling compared to the original survey. Moreover, the subtraction of correlations around random points can subtract residual additive systematics in the signal. As in the case of angular clustering, we therefore generally prefer the density contrast δ instead of the number density n to measure cross-correlations.

5. Finite resolution maps

The preceding sections demonstrate how we can obtain angular power spectra from discrete observations, and how we can relate their expectations to the intrinsic two-point statistics of the observed point processes or random fields. We now turn to the practical task of computing the angular power spectra. Retracing our steps, we find that this can be done in one of two ways:

• (i)

  Compute C_l^nn′, C_l^ff′, and C_l^nf′ directly using their respective expressions (30), (33), and (34).

• (ii)

  Compute n_lm and f_lm from their analytical expansions (29) and (32), then compute the angular power spectra (9) from the spherical harmonic coefficients.

If N is the number of observations, and l_max is the highest angular mode number of interest, then the former method has a runtime complexity of O(N²l_max), i.e. quadratic in the number of observations, which is the same as for real-space estimators. The complexity of the latter method, however, is O(Nl_max²), and it is hence favourable when l_max ≪ N. We therefore generally want to obtain angular power spectra C_l from their constituent coefficients a_lm.

The a_lm computed from the sums (29) and (32) potentially still contain more information than we need: if the observed points are sufficiently dense, they probe scales beyond our desired scale of l_max. We can then reduce the computational complexity further by introducing a spatial binning of the points – or, in other words, by making a map.

Map-making consists of two separate but related parts. The first is sampling, so that spherical functions are represented by their values in a finite set of basis points on the sphere. The number and locations of the basis points are chosen such that it is possible to accurately recover angular modes up to some given l_max from the spherical harmonic expansion (2). Several sampling schemes for that purpose have been proposed; commonly used in astronomy are, e.g. the scheme of Driscoll & Healy (1994), schemes with exact spherical harmonic transforms for band-limited functions (Huffenberger & Wandelt 2010; McEwen & Wiaux 2011), and HEALPix (Górski et al. 2005).

Sampling a random set of point masses, such as the spherical functions (28) and (31) we construct from our discrete observations, with a fixed set of basis points will result in a map that is almost surely zero everywhere. The second part of map-making is hence the collection of function values (i.e. observed points) over a finite region around each sampling point. This is achieved using spherical convolution, which we define below. The area over which observations are collected is, at least in principle, entirely independent of the sampling scheme. Naturally, we want every observation to be counted, in which case this area must be large enough to cover the spaces between sampling points. On the other hand, the area should also not be much larger than necessary, or we needlessly degrade the angular power spectra that we wish to measure. In practice, there is hence always a close match between the convolution and the sampling scheme.

5.1. Spherical convolution

Convolution is a mathematical operation that produces a new function F from a given function f and convolution kernel K. The value of the convolution in a point is obtained by making this point the origin of a local copy of K and computing the integral of f weighted by that kernel. Convolution is therefore not a local operation. And since spin-weighted spherical functions are always defined relative to a local coordinate frame (Boyle 2016), it follows that convolution on the sphere has to explicitly take this non-local nature into account.

As an illustrative example, consider the following situation, where the dot marks the centre of a small, essentially flat patch of the sphere, and the arrows indicate the complex argument of a local spin-weighted function of constant magnitude:

$$\begin{array}{c}\hfill 2pt\begin{array}{ccc}& \downarrow & \\ \to & \mathbf{·}& \leftarrow \\ & \uparrow \end{array}\end{array}$$

Intuitively, the sum of the function values should be zero. Not accounting for coordinate frame effects, this is indeed the case if the dot is near the equator. But if the dot marks the north pole, all arrows point north, and naive summation produces an incorrect result.

For a spherical convolution that treats non-zero spin weights in the correct manner, we define the convolution F of a spherical function f and a symmetric convolution kernel K as

$$\begin{array}{c}\hfill F(\widehat{u})={\displaystyle \int f}(\widehat{u}\prime ){\mathrm{e}}^{\mathrm{i}S\alpha }K(\theta ){\mathrm{e}}^{-\mathrm{i}s{\alpha }^{\prime }}\mathrm{d}\widehat{u}\prime ,\end{array}$$(86)

where the angles θ, α, α′ are the separation and relative orientation of $\widehat{u}$ and $\widehat{u}\prime$ as in the spherical harmonic addition theorem (8). Here, s is the spin weight of the convolved function f, and S is the spin weight of the convolution F. We can choose S freely: under a rotation of γ in $\widehat{u}$, the angle α in the convolution (86) transforms as α ↦ α − γ, so that F indeed picks up the phase factor e^−iS​γ of a function with spin weight S. Overall, the convolution (86) is equivalent to the directional spin-weighted spherical convolution of McEwen et al. (2015) with a symmetric kernel, and reduces to the usual spherical convolution of scalar functions when the spin weight is zero (e.g. Wandelt & Górski 2001).

Most importantly, the definition (86) of spherical convolution yields a useful convolution theorem for spherical harmonic expansions. The convolution kernel K is a function of separation, similar to an angular correlation function, so that we can apply the expansion (5) into Wigner d functions,

$$\begin{array}{c}\hfill K(\theta )={\displaystyle \sum _l}\frac{2l+1}{4\pi }{K}_l{d}_{Ss}^l(\theta ),\end{array}$$(87)

where the coefficients K_l of the expansion are given by (6),

$$\begin{array}{c}\hfill K_l=2\pi {\displaystyle {\int }_0^{\pi }}K(\theta )d_{Ss}^l(\theta )sin(\theta )\mathrm{d}\theta .\end{array}$$(88)

Inserting the expansion (87) into definition (86) and substituting the addition theorem (8) yields the integral (3) for the coefficients f_lm in the spherical harmonic expansion (2) of f. We thus obtain the desired spherical harmonic convolution theorem,

$$\begin{array}{c}\hfill F_{\mathit{lm}}=K_l{f}_{\mathit{lm}},\end{array}$$(89)

where the coefficients F_lm of the convolution F are the product of the coefficients K_l of the convolution kernel K and the coefficients f_lm of the convolved function f.

Carrying out the convolution (86) requires computing the phase factors e^iS​α and e^−isα′ at each point, which can often be done efficiently (see Appendix C). However, if the support of the convolution kernel K is sufficiently small, the local geometry is close to flat, and α ≈ α′. In that case, the phase factors reduce to unity if we choose a convolution with S = s that does not change the spin weight of the function.

It remains to find a tractable convolution kernel K. For a function f with s = 0, the natural choice is a spherical disc of some chosen angular size β > 0. However, for s ≠ 0, the same kernel does not produce analytically tractable coefficients (88) for the convolution theorem. We hence propose a modified convolution kernel that works for any spin weight S = s ≥ 0, and reduces to a spherical disc if s = 0,

$$\begin{array}{c}\hfill K(\theta )=\{\begin{array}{cc}[cos\frac{\theta }{2}{]}^{2s}\hfill & \mathrm{if}\theta \le \beta ,\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\end{array}$$(90)

where β is the angular size (i.e. radius) of the kernel. The effective area of the convolution kernel is⁸

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_K=2\pi {\displaystyle {\int }_0^{\pi }}K(\theta )sin(\theta )\mathrm{d}\theta =4\pi \frac{1-[cos\frac{\beta }{2}{]}^{2s+2}}{s+1}·\end{array}$$(91)

The choice of kernel (90) is firstly motivated by the fact that there is a known expression for its coefficients (88)⁹,

$$\begin{array}{c}\hfill K_l=\{\begin{array}{cc}\frac{4\pi [sin\frac{\beta }{2}][cos\frac{\beta }{2}{]}^{2s+1}{d}_{ss+1}^l(\beta )}{\sqrt{l(l+1)-s(s+1)}}\hfill & \mathrm{if}l>s,\hfill \\ [3pt]4\pi \frac{1-[cos\frac{\beta }{2}{]}^{4s+2}}{2s+1}\hfill & \mathrm{if}l=s.\hfill \end{array}\end{array}$$(92)

Secondly, for small kernel sizes β ≲ 1°, the convolution kernel (90) is essentially a flat spherical disc even when the spin weight is s = 2 (Fig. 2), in which case the coefficients (92) for s = 0 and s = 2 become essentially the same (Fig. 3). This makes the specific kernel (90) a good practical choice for maps when the resolution is below the degree-scale.

In summary, the convolution (86) means that we can create finite resolution maps of the point-mass like number density (28) or field (31) by picking a sampling scheme and for each grid point summing each observed point with the weight K(θ) given by the convolution kernel K, omitting the phase factors in the convolution (86) if the resolution allows it. For a convolution kernel such as (90) with small angular size β ≪ π, an alternative method is to reverse the order of operations, and find all grid points closer than β for each observed point. This can result in vast performance improvements, particularly if the grid points can be queried efficiently, e.g. when using Cartesian or HEALPix grids.

Once maps are created, it suffices to compute their spherical harmonic coefficients F_lm, and reconstruct the coefficients f_lm of the spherical function from the convolution theorem (89),

$$\begin{array}{c}\hfill f_{\mathit{lm}}=\frac{1}{K_l}{F}_{\mathit{lm}}.\end{array}$$(93)

Naturally, this is only possible when K_l ≠ 0, which limits the angular mode numbers l that can be recovered for a given convolution kernel. However, if the deconvolution (93) is possible for all numbers l ≤ l_max, we are readily able to compute the angular power spectrum (9) of f and f′ from the finite-resolution maps F and F′. In cases where the deconvolution is impossible or undesirable, we can instead use the convolution theorem (89) to model the effect of the convolution on the expected angular power spectra. In practice, this can be done at no computational cost, by absorbing the convolution kernel into the mixing matrix (26).

5.2. HEALPix pseudo-convolution

Even for small kernels, the radius search required by the spherical convolution comes at a non-negligible computational cost. Given the number of galaxies observed by Euclid, this cost quickly becomes prohibitive, unless faster, specialised algorithms can be found. However, when using the HEALPix grid for sampling, we can alternatively follow the standard procedure of simply summing the points in each HEALPix pixel. But even though all HEALPix pixels have the same area, this operation is not a true spherical convolution, due to the slightly varying pixel shapes (Górski et al. 2005). Nevertheless, summation over HEALPix pixels does obey the convolution theorem (89) approximately, and the normalised coefficients Ω_K⁻¹K_l are known as the HEALPix pixel window function (Fig. 4). As it turns out, this pseudo-convolution can be adequate for Euclid analysis, which we demonstrate in Sect. 6.

There is, however, one fundamental difference between a true spherical convolution and HEALPix pseudo-convolution. According to the convolution theorem (89), the convolution kernel K_l is imprinted on all spherical harmonic coefficients of a map, and consequently on the angular power spectrum (9),

$$\begin{array}{c}\hfill C_l^{F{F}^{\prime }}=K_l^2{C}_l^{f{f}^{\prime }}.\end{array}$$(94)

In particular, it follows that the convolution kernel should also affect additive bias terms such as A^nn′ in the spectrum (37) of point processes, or A^ff′ in the expected spectrum (59) of random fields. These bias terms will therefore no longer be constant after a true convolution. This is not the case for the HEALPix pseudo-convolution (Fig. 5), since the non-overlapping HEALPix pixels cannot imprint structure, such as the convolution kernel, below the pixel scale. In fact, a similar effect occurs when sampling a true convolution so sparsely that the convolution kernel areas no longer overlap.

There is a practical consequence of this difference between true convolution and HEALPix pseudo-convolution. For a map created by true convolution, the deconvolution (93) turns an additive bias term in the angular power spectrum back into a constant. For a HEALPix map, the same deconvolution of the pixel window function turns an additive bias A into an l-dependent bias A/K_l². When comparing HEALPix spectra and their expectations, the additive bias must therefore either be subtracted from the measured spectra before deconvolution of the pixel window function, or the correct l-dependent bias must be used, e.g. in expectations (49), (52), and (59).

5.3. Maps

We can now define the maps we make. For simplicity, we always call one value of a finite map a pixel with area Ω_K, with the understanding that this may refer either to an actual HEALPix pixel or to the kernel of a true convolution.

The convolution of the number density n with definition (28) is the map N of number counts in each pixel. The convolution of $\overline{n}$ is the map $\overline{N}$ of mean number counts; we write it as $\overline{N}={\overline{N}}_{0}V$ using the mean number of points per pixel ${\overline{N}}_0={\mathrm{\Omega }}_K{\overline{n}}_0$ and a map V that is the convolution of v divided by Ω_K. We call V the ‘visibility map’, since it is the pixel-averaged equivalent of the visibility v, with pixel values between 0 and 1. Since the convolution (86) is a linear operation, the expectation (44) translates from the number densities to the number count maps, $⟨N⟩=\overline{N}$.

To isolate the clustering signal in the number count map N, we define a map Δ for the density contrast (50),

$$\begin{array}{c}\hfill \mathrm{\Delta }=\frac{N-\overline{N}}{{\overline{N}}_0}·\end{array}$$(95)

It is clear that Δ is the convolution of δ, but divided by the pixel area, so that the numerical values of Δ have the correct, intuitive scale where −1 means “empty space”. Deconvolution of Δ must therefore be carried out with the normalised convolution kernel Ω_K⁻¹ K_l.

For a field f (e.g. cosmic shear), we compute the map F by summing the weighted field values w_k f_k in each pixel, and dividing the result by a constant mean pixel weight ${\overline{W}}_0$. Specifically, we choose ${\overline{W}}_0$ to be the mean weight divided by the mean visibility, as computed from the maps, which makes ${\overline{W}}_0$ relatively insensitive to the survey footprint and systematics. The map F is hence the convolution of the function f with definition (31), divided by ${\overline{W}}_0$ to remove explicit dependencies on pixel area and overall weight factors. Similarly, we compute the weight map W as the total weight in each pixel, i.e. the convolution of w with definition (54), divided by ${\overline{W}}_0$.

Since ${\overline{W}}_0$ contains a factor of the pixel area, deconvolution of F and W is carried out with the normalised convolution kernel Ω_K⁻¹ K_l. Furthermore, since both F and W are scaled identically, the resulting mixing matrix automatically applies the correct factors of ${\overline{W}}_0$ to the expected angular power spectra. However, we do need to account for the scaling by ${\overline{W}}_0$ when computing any additive bias terms.

In particular, we do not average the field values in each pixel by dividing the map F by the map W. For Euclid, the resolution of our maps is such that about half of all observed pixels contain fewer than two observed values. Using a weighted average would simply divide out the given weights in these pixels, resulting in an unweighted cosmic shear map (see, e.g. Hikage et al. 2011, 2019; Nicola et al. 2021). This is clear when looking at the spherical functions f and w with definitions (31) and (54), respectively: dividing f by w where both are non-zero is equivalent to using unit weights in f.

6. Validation

In the preceding sections, we derive the overall theory of angular power spectra from discrete sets of observations, their expectations, and ways to efficiently compute spectra from maps. We now turn to the validation of our findings. One part of this are the explicit assumptions that we have made throughout; these represent specific scientific questions that are partially the subject of active research in their own right, and we do not investigate their validity here.

In what follows, we validate our specific implementation of the methodology described above. This is a publicly available code called Heracles¹, developed within the Euclid Science Ground Segment. It contains routines for catalogue reading, map-making, spherical harmonic transforms, angular power spectra, and mixing matrices. The code can be used as a Python library, e.g. for data exploration in a notebook interface, or via a standalone command-line interface, e.g. for batch data processing. In particular, the code also contains an implementation of the discrete angular power spectrum methodology, which is based on fast, non-uniform computation of spherical harmonics (Reinecke et al. 2023) as implemented in the ducc package¹⁰.

To validate the performance of Heracles, we carry out the following series of tests:

• (i)

  We estimate the mean density of galaxies from the visible sky fraction in the various Euclid data releases,

• (ii)

  we test if the phase factors in the spherical convolution can be neglected for map-based spectra,

• (iii)

  we assess the overall accuracy of our measurements, and

• (iv)

  we apply the methodology in a data processing setting that mimics the first Euclid data release.

Where simulations are created, we generally employ the same flat ΛCDM cosmology as Euclid’s Flagship simulation (Euclid Collaboration: Castander et al. 2025), with parameter values Ω_m = 0.319, Ω_b = 0.049, A_s = 2.1 × 10⁻⁹, n_s = 0.96, and h = 0.67.

6.1. Mean density estimation

Constructing the density contrast (50) requires knowledge of the mean density of galaxies over the entire sky, which we must estimate from the visible sky fraction. If our estimate is inaccurate, we bias the angular power spectrum C_l^δδ′ in a non-trivial manner with respect to the expectation (52). The problem in estimating the mean density accurately is that the visible part of the sky might be particularly over- or underdense compared to the true mean. The probability of this depends on the area of the observed sky, as well as the typical size of large-scale density fluctuations, and hence the clustering of points. This problem is closely related to the integral constraint for real-space estimators.

To test the impact on Euclid observations, we generate 10 000 lognormal realisations of a galaxy distribution with a linear galaxy bias (Tessore et al. 2023). To account for the redshift evolution of galaxy clustering, we test redshifts z = 0.5, 1.0, 2.0, using a redshift-dependent bias that was fitted to the Euclid Flagship simulation (Euclid Collaboration: Lepori et al. 2022). Since the error in the mean density is a function of visible sky fraction, we further use a representative footprint (Euclid Collaboration: Scaramella et al. 2022) for each of the three Euclid data releases DR1 (1 year, 6% sky coverage), DR2 (3 years, 18% sky coverage), and DR3 (6 years, 36% sky coverage). The results show that the expected relative error in the mean density is at the per mille level for all data releases and redshifts, with a scatter that stays below the per cent level for DR2 and beyond (Fig. 6).

6.2. Phase factors

To test whether or not we can neglect the phase factors in the convolution (86), we need to quantify their impact on maps at the required resolution for Euclid. Consider a fixed pixel located at a position ${\widehat{u}}_0$. By neglecting the phase factors, the convolution (86) is approximated as

$$\begin{array}{c}\hfill F({\widehat{u}}_0)\approx {\displaystyle \int f}(\widehat{u}\prime )K(\theta )\mathrm{d}\widehat{u}\prime .\end{array}$$(96)

To cancel the phase factors, the approximation effectively applies a position-dependent multiplicative bias e^{−i (S​α − sα′)} to the function f over the pixel area, where the angles α and α′ are taken with respect to ${\widehat{u}}_0$. Under a rotation of γ′ in $\widehat{u}\prime$, the angle α′ transforms as α′↦α′−γ′; as a function of $\widehat{u}\prime$, the bias therefore has a spin weight of s.

For cosmic shear, we can make maps of this bias, using the expressions from Appendix C and setting S = s = 2. The result is shown in Fig. 7 for HEALPix maps with resolution parameters NSIDE = 1, 4, 16, 64. Parallel transport along a meridian has phase factors of unity, so that the phase factor bias is effectively a function of azimuthal distance from the pixel centre. The pixel resolution starts to fall below the degree-scale at NSIDE = 64, and the phase factor bias becomes small, due to the essentially flat geometry of the pixels.

Overall, we expect no impact from neglected phase factors for HEALPix maps with resolution parameter NSIDE ≳ 1024. If necessary, the phase factor bias could be mitigated even further by choosing the coordinate system such that the poles fall into a masked region, e.g. the galactic plane.

6.3. Accuracy

We now validate the results we derive in Sect. 4 with simulations. To characterise the accuracy of our measurement, these require many realisations that, for Euclid, cover a significant fraction of the sky. Usually, lognormal simulations would be the method of choice here; however, the transformations involved in sampling lognormal fields are not exact (Tessore et al. 2023), which adds an element of uncertainty to the validation. Gaussian random fields can be simulated with exactly prescribed two-point statistics; however, for realistic angular power spectra and values of linear galaxy bias, the realised density fields almost surely contain regions where the number density becomes negative. Instead, we use squared Gaussian random fields, which is a toy model we develop in Appendix E that can be simulated exactly and remains physically valid everywhere.

To accurately quantify the effects of map-based measurements, we require simulations that are not themselves affected by pixel effects. We therefore simulate the fields not in real space, but via their spherical harmonic expansion (2). We then sample points using a rejection sampling scheme that accepts or rejects points with a probability that is proportional to the simulated density field, evaluated in each sampled point from the spherical harmonic expansion. Instead of a survey footprint map, we draw points from a spherical cap of 2500 deg², matching the anticipated area of Euclid’s DR1, and located at the centre of the Euclid Flagship simulation (Euclid Collaboration: Castander et al. 2025). The weak lensing fields are subsequently evaluated at the sampled positions from their spherical harmonic expansions, without any intermediary interpolation.

We then generate 10 000 realisations of these simulations. To show results for both auto- and cross-spectra, we simulate two Gaussian tomographic redshift bins centred on representative redshifts of z = 0.5, 1.0 with a width of σ_z = 0.125. To simulate galaxy clustering, we use realistic redshift-dependent galaxy bias (Euclid Collaboration: Lepori et al. 2022), with a Euclid-like galaxy density of 2 galaxies per arcmin² in each tomographic bin. For each simulated galaxy, we store its position, as well as its observed ellipticity from weak lensing with a random intrinsic ellipticity drawn from a hyperbolic normal distribution (Tessore et al. 2023), using a per-component standard deviation σ_ϵ = 0.26. Instead of using the correct weak lensing action (56), we simply sum the intrinsic ellipticity and the shear γ from weak lensing, since we might otherwise pick up biases due to the reduced shear approximation (Deshpande et al. 2020). Finally, we give every galaxy a random shear weight from a log-uniform distribution between 10⁻² and 10², to simulate a dynamic range that should exceed any real shear measurement method.

For each simulation, we measure the angular power spectra of the reconstructed density field δ, the cosmic shear E- and B-mode (Appendix D), and their cross-correlations in the two tomographic bins for angular modes up to l = 3000. In addition, we measure the angular power spectra of the visibility V and the shear weight W in each bin for modes up to l = 9000, from which we compute the mixing matrices (25) truncated at l = 3000 and l₁ = 6000. By construction, our simulations are band-limited at l = 6000, so that these truncated mixing matrices should contain all non-zero entries, and yield exact expectations for our measurements.

We thus obtain measured and expected angular power spectra for all combinations of probes across the two tomographic bins: angular clustering, cosmic shear, and galaxy–galaxy lensing. To reduce noise, we average the spectra over 32 angular bins with logarithmic spacing between l = 10 and l = 3000. We then compute the mean of the bias ΔC_l = C_l − ⟨C_l⟩ between measurements and expectations, which we scale relative to the standard deviation σ_l of each measurement over the set of realisations.

We carry out the measurements for each simulation using (i) the exact angular power spectra computed from discrete sets of points, and (ii) map-based angular power spectrum from HEALPix maps with resolution parameter NSIDE = 4096. The results are shown in Fig. 8. For the exact, discrete spectra, we find agreement at the 1% level relative to the standard deviation, except for effects that can be ascribed to source–lens clustering (i.e. Assumption 6). For the HEALPix-based spectra, the results show an additional per-cent-level bias in angular clustering at small scales. Further testing reveals that the small-scale HEALPix bias has a dependency on the location of the survey footprint, and we can hence ascribe it to the pseudo-convolution with varying pixel shapes. The dependence of the recovered angular power spectra on pixel shapes is more thoroughly explored elsewhere (Hall & Tessore, in prep.).

To demonstrate that the bias in Fig. 8 is in fact source–lens clustering, we run a second set of simulations where the positions of shears are distributed according to an independent (but identically clustered) realisation of large-scale structure. When source–lens clustering is thus taken into account, the discrete spectra show a relative bias that is consistent with zero at the 1% level for all probes.

Overall, we find that the HEALPix-based measurements are only marginally biased with respect to the exact discrete angular power spectra. In light of the lower computational cost, we therefore adopt this method for Euclid’s DR1 analysis, which will enable faster turnaround in the data processing. However, since accuracy is a function of survey area and galaxy density, this may no longer be the case for subsequent data releases.

6.4. Applicability to Euclid DR1

To demonstrate that we have a viable pipeline for Euclid’s first data release (DR1), we process a realistic DR1-like data volume. To do so, we select galaxies contained in the provisional northern DR1 footprint (Euclid Collaboration: Scaramella et al. 2022) from the Euclid Flagship simulation (Euclid Collaboration: Castander et al. 2025), obtained from CosmoHub (Tallada et al. 2020; Carretero et al. 2017). The Euclid data processing pipeline aims to support up to 13 tomographic redshift bins (Euclid Collaboration: Mellier et al. 2025), and since the number of spectra and mixing matrices increases quadratically with the number of tomographic bins, we want to ensure compliance with such a setting. Using the simulated photometric redshifts, we therefore bin galaxies into 13 equi-populated tomographic redshift bins. We then measure all 780 possible auto- and cross-correlations between galaxy positions and cosmic shear E- and B-modes. To compare the measurements with expectations, we further compute mixing matrices for all spectra from the simulated visibility and shear weight maps. For this test, we apply the HEALPix-based methodology, with resolution parameter NSIDE = 4096 and maximum angular mode number l = 5000 for all probes, which exceeds the “optimistic” forecast of scale cuts (Euclid Collaboration: Blanchard et al. 2020).

The results are shown in Figs. A.1 and A.2. To compute the expected spectra, we obtain theoretical full-sky predictions with the Cosmology Likelihood for Observables in Euclid code, CLOE (Euclid Collaboration: Joudaki et al., in prep.), using the implemented prescription for photometric harmonic-space observables (Euclid Collaboration: Cardone et al., in prep.). We use HMCode2020 (Mead et al. 2021) to model the non-linear matter power spectrum, as provided in the public code CAMB (Challinor & Lewis 2011). For the background cosmology, we use the same parameter values as in the Euclid Flagship simulation (Euclid Collaboration: Castander et al. 2025). To compute the observables, we use the simulated redshift distributions n(z), and a linear galaxy bias measured from the Euclid Flagship simulation (Euclid Collaboration: Lepori et al. 2022). Since we are only testing the feasibility of the data processing here, we do not perform any fine-tuning of the non-linear modelling, or take into account systematic effects such as magnification bias or intrinsic alignments. This is visible, e.g. in the angular clustering at small scales l ≳ 1000. Nevertheless, we obtain a level of agreement between measurements and expectations that is in line with previous results (Euclid Collaboration: Castander et al. 2025).

Processing the data from DR1-like catalogues to all 780 angular power spectra is very fast: obtaining the spherical harmonic expansions of the fields and weights in one tomographic bin takes around three minutes of wall-clock time, and all tomographic bins can be processed in parallel. The subsequent computation of angular power spectra from all combinations of spherical harmonic coefficients of the fields has negligible cost. Computation of the mixing matrices from the spherical harmonic coefficients of the weights is a more resource-intensive operation, taking around 75 CPU core-hours in total; however, all mixing matrices can be computed in parallel as necessary. Overall, we therefore expect no significant impact on Euclid’s DR1 processing from the measurement of angular power spectra.

7. Discussion and conclusion

We derived a complete framework to obtain the exact measurements and expectations for the angular power spectra from discrete sets of data. Starting from an exact map-free formalism, we find new results such as exact non-stochastic expressions for the additive (noise) biases (36) and (60) for angular clustering and random fields, respectively. Furthermore, by explicitly tracking what assumptions enter our measurements, we are able to separate the methodological accuracy of our results (i.e. when all assumptions are true) from the true accuracy of the results, for example in cases of intrinsic alignments and source–lens clustering, which need to be treated at the level of theoretical predictions. When validating our results on simulations, we find that discrete angular power spectra can achieve biases of less than 1% with respect to their standard deviation in a Euclid DR1-like setting. Overall, we are therefore confident in our ability to measure angular power spectra to the very high level of accuracy required to achieve Euclid’s ambitious science goals.

Using the theory of spherical convolution (86), we can connect our exact theory with the standard practice of measuring angular power spectra from HEALPix maps, for example. Conceptually, this is a step away from the picture in which the observed maps are ‘noisy tracers’ of an underlying continuous field, such as the true galaxy density or the true cosmic shear. In practice, the main difference between the exact map-based formalism and standard practice is that observations are summed over one map pixel (i.e. convolution kernel), but are not averaged. This side-steps common issues with map-based spectra (e.g. empty pixels) or the fact that pixels containing just one single observation point are effectively unweighted after averaging.

When analysing Euclid DR1-like simulations, we find that HEALPix-based spectra can perform at a level of accuracy that is comparable to the discrete spectra. An exception is the case of angular galaxy clustering, where the HEALPix pseudo-convolution, due to varying pixel shapes, introduces a position-dependent bias at the per cent level for angular mode numbers l at or above the resolution parameter NSIDE. However, since this bias is both small and mitigable by relevant analysis choices (i.e. map resolution and scale cuts), we plan to employ the map-based methodology for Euclid’s first data release.

We make our implementation of the methodology presented here available in the form of a code called Heracles¹. This code, originally created for 3 × 2 pt data processing in the Euclid Science Ground Segment, was designed from the ground up to be user-friendly and widely applicable to any given probe and survey, and will be maintained for public use.

---

Acknowledgments

We thank Kevin Wolz and David Alonso for their collaboration during the completion of this manuscript. We also thank Ben Wandelt for a number of excellent comments, and the anonymous referee for improving our manuscript. NT thanks Martin Reinecke for pointing out the possibility of implementing the discrete angular power spectrum computation using ducc. This work was supported by UK Space Agency grants ST/W002574/1 and ST/X00208X/1. AH acknowledges support through a Royal Society University Research Fellowship. GCH acknowledges support through the ESA research fellowship programme. This work has made use of CosmoHub, developed by PIC (maintained by IFAE and CIEMAT) in collaboration with ICE-CSIC. CosmoHub received funding from the Spanish government (MCIN/AEI/10.13039/501100011033), the EU NextGeneration/PRTR (PRTR-C17.I1), and the Generalitat de Catalunya. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft-und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid website (https://www.euclid-ec.org/).

References

Appendix A: Additional figures

Appendix B: Alternative estimators for angular clustering

In Sect. 4.1, we measure the two-point statistics for point processes using the particular choice of density contrast (50). Here, we consider a number of alternative choices.

Firstly, we can trivially replace the mean number density $\overline{n}$ by a catalogue of random points (‘randoms’). Given the definition (42) of the mean number density, these randoms must be distributed according to the visibility v (Baleato Lizancos & White 2024).

Secondly, in full correspondence to real-space methods, we can directly construct estimators of the angular power spectrum 𝔴_l from the expectation (49). For example, using a formal inverse ${(M^{\overline{n}\overline{n}\prime })}^{-1}$ of the mixing matrix, we can construct the estimator

$$\begin{array}{c}\hfill {\widehat{\mathfrak{w}}}_l^{\mathrm{N}}={\displaystyle \sum _{l^{\prime }}}{(M^{\overline{n}{\overline{n}}^{\prime }})}_{l{l}^{\prime }}^{-1}[C_{l^{\prime }}^{n{n}^{\prime }}-C_{l^{\prime }}^{\overline{n}{\overline{n}}^{\prime }}-A^{n{n}^{\prime }}].\end{array}$$(B.1)

In the taxonomy of Kerscher et al. (2000), this corresponds to the “natural” real-space estimator (DD − RR)/RR. Furthermore, by expectation (44), we have $⟨C_l^{n\overline{n}\prime }⟩=C_l^{\overline{n}\overline{n}\prime }$, and we can hence construct a more advanced estimator

$$\begin{array}{c}\hfill {\widehat{\mathfrak{w}}}_l^{\mathrm{LS}}={\displaystyle \sum _{l^{\prime }}}{(M^{\overline{n}{\overline{n}}^{\prime }})}_{l{l}^{\prime }}^{-1}[C_{l^{\prime }}^{n{n}^{\prime }}-C_{l^{\prime }}^{n{\overline{n}}^{\prime }}-C_{l^{\prime }}^{\overline{n}{n}^{\prime }}+C_{l^{\prime }}^{\overline{n}{\overline{n}}^{\prime }}-A^{n{n}^{\prime }}].\end{array}$$(B.2)

This is the harmonic-space equivalent of the Landy & Szalay (1993) estimator (DD − DR − RD + RR)/RR.

The estimators (B.1) and (B.2) both rely on inversion of the mixing matrix. We can similarly construct a partial-sky variant of the natural estimator in harmonic space,

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathfrak{w}}}_l^{\mathrm{N}}=C_l^{n{n}^{\prime }}-C_l^{\overline{n}{\overline{n}}^{\prime }}-A^{n{n}^{\prime }},\end{array}$$(B.3)

as well as a partial-sky variant of the harmonic-space Landy-Szalay estimator,

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathfrak{w}}}_l^{\mathrm{LS}}=C_l^{n{n}^{\prime }}-C_l^{n{\overline{n}}^{\prime }}-C_l^{\overline{n}{n}^{\prime }}+C_l^{\overline{n}{\overline{n}}^{\prime }}-A^{n{n}^{\prime }}.\end{array}$$(B.4)

The respective expectation of both partial-sky estimators is the product of mixing matrix and full-sky expectation. In particular, the partial-sky Landy–Szalay estimator (B.4) is essentially the same as the measured angular power spectrum (51) of the density contrast (50).

Lastly, to see why we normalise the density contrast (50) by a constant ${\overline{n}}_0$, consider an alternative definition of the density contrast with an arbitrary normalisation function q,

$$\begin{array}{c}\hfill {\delta }_q(\widehat{u})=\frac{n(\widehat{u})-\overline{n}(\widehat{u})}{{\overline{n}}_{0}q(\widehat{u})}.\end{array}$$(B.5)

It follows from the definitions of the number density (28) and mean number density (42) that δ_q is equivalent to δ under a change of weights $w_k\mapsto w_k/q({\widehat{u}}_k)$ and, consequently, a change of visibility $v(\widehat{u})\mapsto v(\widehat{u})/q(\widehat{u})$. Defining the density contrast (50) with a different normalisation therefore effectively replaces the given set of weights with a different set of weights.

Appendix C: Relative orientation on the sphere

To obtain the angles α and α′ that describe the relative orientation of points $\widehat{u}$ and $\widehat{u}\prime$ on the sphere, it suffices to solve the spherical triangle shown in Fig. C.1 (for more information, see Hall & Tessore in prep.),

$$\begin{array}{cc}\hfill \alpha & =arctan\frac{sin\vartheta \prime sin(\phi -\phi \prime )}{sin\vartheta cos\vartheta \prime -cos\vartheta sin\vartheta \prime cos(\phi -\phi \prime )},\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill \alpha \prime & =-arctan\frac{sin\vartheta sin(\phi -\phi \prime )}{sin\vartheta \prime cos\vartheta -cos\vartheta \prime sin\vartheta cos(\phi -\phi \prime )}.\hfill \end{array}$$(C.2)

The respective numerators and denominators are written here such that their signs yield the correct quadrant for the inverse tangent.

The same angles can be expressed in terms of the components of the unit vectors $\widehat{u}=\{x,y,z\}$ and $\widehat{u}\prime =\{x\prime ,y\prime ,z\prime \}$ as

$$\begin{array}{cc}\hfill \alpha & =arctan\frac{y{x}^{\prime }-x{y}^{\prime }}{z^{\prime }-zcos\theta },\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill {\alpha }^{\prime }& =-arctan\frac{y{x}^{\prime }-x{y}^{\prime }}{z-z^{\prime }cos\theta },\hfill \end{array}$$(C.4)

with $cos\theta =\widehat{u}·\widehat{u}\prime$. This form is often useful in applications where points are available as vectors, since it requires no additional trigonometric operations. In fact, for α = arctan(q/p), we find a familiar expression for the spin-2 phase factors that appear, e.g. in the spherical harmonic addition theorem (8),

$$\begin{array}{c}\hfill {\mathrm{e}}^{2\mathrm{i}\alpha }=\frac{p^2-q^2+2\mathrm{i}pq}{p^2+q^2}.\end{array}$$(C.5)

The phase factors can hence be computed entirely in terms of the vector components of $\widehat{u}$ and $\widehat{u}\prime$.

Appendix D: Decomposition into E- and B-modes

If the spherical function f is a complex-valued random field, the two-point statistics of f and a second, not necessarily distinct, field f′ are not fully characterised by the expected angular correlations (16) alone. Like for any complex random variable, we also require the associated pseudo-correlation ⟨C^f*​f′(θ)⟩, i.e. the correlation of the complex conjugated random field f^* and f′,

$$\begin{array}{c}\hfill ⟨f(\widehat{u})f^{\prime }({\widehat{u}}^{\prime })⟩={\mathrm{e}}^{\mathrm{i}s\alpha }⟨C^{f^{\ast }{f}^{\prime }}(\theta )⟩{\mathrm{e}}^{\mathrm{i}{s}^{\prime }{\alpha }^{\prime }},\end{array}$$(D.1)

where we have used the fact that the spin weight of f^* is −s. The same information is contained in the pseudo-spectrum¹¹

$$\begin{array}{c}\hfill ⟨{(f^{\ast })}_{\mathit{lm}}^{\ast }{f}_{l^{\prime }{m}^{\prime }}^{\prime }⟩={\delta }_{l{l}^{\prime }}^{\mathrm{K}}{\delta }_{m{m}^{\prime }}^{\mathrm{K}}⟨C_l^{f^{\ast }{f}^{\prime }}⟩,\end{array}$$(D.2)

which is merely expectation (13) applied to f^* and f′.

Instead of using spectra and pseudo-spectra, it is often more convenient to work with a different decomposition of the harmonic-space two-point statistics, namely that into E- and B-modes (Zaldarriaga & Seljak 1997). For a spherical function f with spin weight s, the respective E- and B-modes are defined as linear combinations of the spherical harmonic coefficients of f and f^*,

$$\begin{array}{cc}\hfill E_{\mathit{lm}}& =-\frac{f_{\mathit{lm}}+{(-1)}^s{(f^{\ast })}_{\mathit{lm}}}{2},\hfill \end{array}$$(D.3)

$$\begin{array}{cc}\hfill B_{\mathit{lm}}& =-\frac{f_{\mathit{lm}}-{(-1)}^s{(f^{\ast })}_{\mathit{lm}}}{2\mathrm{i}},\hfill \end{array}$$(D.4)

where the overall negative sign is the convention adopted by HEALPix. The E- and B-mode spectra are then obtained by using E_lm and B_lm in the angular power spectrum (9). Since the coefficients (D.3) and (D.4) are linear combinations of f_lm and (f^*)_lm, it is clear that the resulting E- and B-mode spectra are linear combinations of the spectra C_l^ff′, C_l^f*​f′, etc. of the constituent fields and their complex conjugates.

To compute, e.g. the expectation for partial-sky E- and B-mode spectra, it therefore suffices to apply the mixing matrix formalism described above to the individual spectra of the fields, and express the result in terms of the full-sky E- and B-mode spectra. In doing so, one finds that the mixing matrices also introduce mixing between E- and B-modes (Brown et al. 2005).

Appendix E: Squared normal fields

In this section we consider a Gaussian random field X on the sphere that is transformed into a random field Y by an arbitrary function t,

$$\begin{array}{c}\hfill Y(\widehat{u})=t(X(\widehat{u})).\end{array}$$(E.1)

As shown by Tessore et al. (2023), a band-limited angular power spectrum G_l for X does not generally result in a band-limited angular power spectrum C_l for Y. In practice, we are not generally able to construct a spectrum G_l that, after transformation, reproduces a desired spectrum C_l exactly. For validation, we now try and identify a special case where that is possible. In particular, we look for a transformation with the following two characteristics:

• i)

  The transformed field Y is bounded from below, so that we can simulate density contrasts δ for angular clustering that respect the physical bound δ ≥ −1.

• ii)

  The transformed field Y has a band-limited spectrum.

Lognormal fields satisfy the first criterium, but lack a strictly band-limited spectrum. However, it turns out that both criteria are fulfilled by squaring a Gaussian random field.

Let X be a normal random variable with zero mean and variance σ² ≤ 1, and let $a=\sqrt{1-{\sigma }^2}$. Define the random variable Y as a quadratic transformation of X, resulting in a scaled and shifted non-central χ² random variable with 1 degree of freedom,

$$\begin{array}{c}\hfill Y=\lambda [{(X-a)}^2-1],\end{array}$$(E.2)

where λ > 0 is the scale parameter of the distribution, which also fixes the minimum value of Y.¹² A straightforward calculation shows that Y has zero mean and variance

$$\begin{array}{c}\hfill ⟨Y^2⟩=2{\lambda }^2{\sigma }^2(2-{\sigma }^2).\end{array}$$(E.3)

Inserting σ² = 1 − a², the variance can equivalently be expressed in terms of a,

$$\begin{array}{c}\hfill ⟨Y^2⟩=2{\lambda }^2(1-a^4),\end{array}$$(E.4)

and the value of a can hence be obtained from the variance of the transformed random variable,

$$\begin{array}{c}\hfill a=(1-\frac{⟨Y^2⟩}{2{\lambda }^2}{)}^{\frac{1}{4}}.\end{array}$$(E.5)

The transformation (E.2) is therefore readily obtained in either direction. Setting λ = 1, the standardised probability distribution function of Y is

$$\begin{array}{c}\hfill f_Y(y;a)=\frac{exp(-\frac{a^2+y+1}{2(1-a^2)})cosh(\frac{a\sqrt{y+1}}{1-a^2})}{\sqrt{2\pi (1-a^2)(y+1)}}.\end{array}$$(E.6)

The distribution approaches normality for small variances, and becomes more skewed as the variance increases, similar to the lognormal distribution (Fig. E.1).

We then apply the transformation (E.2) pointwise to a pair X and X′ of jointly homogeneous Gaussian random fields on the sphere. By expectation (16), there is a correlation function G such that $⟨X(\widehat{u})X\prime (\widehat{u}\prime )⟩=G(\theta )$. It can be shown that the transformed fields Y and Y′ are also jointly homogeneous (Tessore et al. 2023), and there is hence an angular correlation function C such that $⟨Y(\widehat{u})Y\prime (\widehat{u}\prime )⟩=C(\theta )$. Using the transformation (E.2), we can compute C(θ) in terms of G(θ),

$$\begin{array}{c}\hfill C(\theta )=2\lambda \lambda \prime G(\theta )(G(\theta )+2aa\prime ).\end{array}$$(E.7)

By completing the square, we also obtain the inverse relation,

$$\begin{array}{c}\hfill G(\theta )=\sqrt{\frac{C(\theta )}{2\lambda \lambda \prime }+{(aa\prime )}^2}-aa\prime .\end{array}$$(E.8)

Furthermore, relation (E.7) is readily transformed to harmonic space using expectation (24),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill C_l& =2\lambda {\lambda }^{\prime }{\displaystyle \sum _{l_1{l}_2}}\frac{(2l_1+1)(2l_2+1)}{4\pi }{\left(\begin{array}{ccc}{l}_1& l_2& l\\ 0& 0& 0\end{array}\right)}^2{G}_{l_1}{G}_{l_2}\hfill \\ \hfill & +4\lambda {\lambda }^{\prime }a{a}^{\prime }{G}_l.\hfill \end{array}\end{array}$$(E.9)

Since the triangle condition l ≤ l₁ + l₂ is imposed on the sum by the Wigner 3j symbols, it follows that C_l is indeed band-limited if G_l is band-limited, at twice the angular mode number. In this case, the non-linear solver for G_l proposed by Tessore et al. (2023) can produce an essentially exact transformation.

For our simulations, the variance of each random field, and hence a, is determined by its spectrum. It remains to fix the value of λ. For angular clustering, we simulate the density contrast δ, for which we set λ = 1. For cosmic shear, we simulate the convergence field κ, for which we set λ using the fitting formula of Hilbert et al. (2011).

All Figures