# Euclid preparation. Angular power spectra from discrete observations

N Tessore, B Joachimi, A Loureiro, A Hall, G Cañas-Herrera, I Tutusaus, N Jeffrey, K Naidoo, J.D Mcewen, A Amara, et al.

► To cite this version:

N Tessore, B Joachimi, A Loureiro, A Hall, G Cañas-Herrera, et al. Euclid preparation. Angular power spectra from discrete observations. Astron. Astrophys., 2025, 694, pp.A141. 10.1051/0004-6361/202452018. [arXiv:2004.06361]

HAL Id: hal-04685451

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

Submitted on 9 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Euclid preparation

## LIX. Angular power spectra from discrete observations

Euclid Collaboration: N. Tessore¹,*, B. Joachimi¹, A. Loureiro²,³, A. Hall⁴, G. Cañas-Herrera⁵,⁶, I. Tutusaus⁷, N. Jeffrey¹, K. Naidoo¹, J. D. McEwen⁸, A. Amara⁹, S. Andreon¹⁰, N. Auricchio¹¹, C. Baccigalupi¹²,¹³,¹⁴,¹⁵, M. Baldi¹⁶,¹¹,¹⁷, S. Bardelli¹¹, F. Bernardeau¹⁸,¹⁹, D. Bonino²⁰, E. Branchini²¹,²²,¹⁰, M. Brescia²³,²⁴,²⁵, J. Brinchmann²⁶,²⁷, A. Caillat²⁸, S. Camera²⁹,³⁰,²⁰, V. Capobianco²⁰, C. Carbone³¹, V. F. Cardone³²,³³, J. Carretero³⁴,³⁵, S. Casas³⁶,³⁷, M. Castellano³², G. Castignani¹¹, S. Cavuoti²⁴,²⁵, A. Cimatti³⁸, C. Colodro-Conde³⁹, G. Congedo⁴, C. J. Conselice⁴⁰, L. Conversi⁴¹,⁴², Y. Copin⁴³, F. Courbin⁴⁴,⁴⁵,⁴⁶, H. M. Courtois⁴⁷, M. Cropper⁸, A. Da Silva⁴⁸,⁴⁹, H. Degaudenzi⁵⁰, G. De Lucia¹³, J. Dinis⁴⁸,⁴⁹, F. Dubath⁵⁰, C. A. J. Duncan⁴⁰, X. Dupac⁴², S. Dusini⁵¹, M. Farina⁵², S. Farrens⁵³, F. Faustini⁵⁴,³², S. Ferriol⁴³, M. Frailis¹³, E. Franceschi¹¹, M. Fumana³¹, S. Galeotta¹³, W. Gillard⁵⁵, B. Gillis⁴, C. Giocoli¹¹,⁵⁶, P. Gómez-Alvarez⁵⁷,⁴², A. Grazian⁵⁸, F. Grupp⁵⁹,⁶⁰, L. Guzzo⁶¹,¹⁰, S. V. H. Haugan⁶², H. Hoekstra⁶³, W. Holmes⁶⁴, F. Hormuth⁶⁵, A. Hornstrup⁶⁶,⁶⁷, P. Hudelot¹⁹, K. Jahnke⁶⁸, M. Jhabvala⁶⁹, E. Keihänen⁷⁰, S. Kermiche⁵⁵, A. Kiessling⁶⁴, B. Kubik⁴³, M. Kümmel⁶⁰, M. Kunz⁷¹, H. Kurki-Suonio⁷²,⁷³, S. Ligori²⁰, P. B. Lilje⁶², V. Lindholm⁷²,⁷³, I. Lloro⁷⁴, G. Mainetti⁷⁵, E. Maiorano¹¹, O. Mansutti¹³, O. Marggraf⁷⁶, M. Martinelli³²,³³, N. Martinet²⁸, F. Marulli⁷⁷,¹¹,¹⁷, R. Massey⁷⁸, E. Medinaceli¹¹, S. Mei⁷⁹, M. Melchior⁸⁰, Y. Mellier⁸¹,¹⁹, M. Meneghetti¹¹,¹⁷, E. Merlin³², G. Meylan⁴⁴, J. J. Mohr⁶⁰,⁵⁹, M. Moresco⁷⁷,¹¹, B. Morin⁵³, L. Moscardini⁷⁷,¹¹,¹⁷, E. Munari¹³,¹², R. Nakajima⁷⁶, S.-M. Niemi⁵, C. Padilla⁸², S. Paltani⁵⁰, F. Pasian¹³, K. Pedersen⁸³, W. J. Percival⁸⁴,⁸⁵,⁸⁶, V. Pettorino⁵, S. Pires⁵³, G. Polenta⁵⁴, M. Poncet⁸⁷, L. A. Popa⁸⁸, F. Raison⁵⁹, A. Renzi⁸⁹,⁵¹, J. Rhodes⁶⁴, G. Riccio²⁴, E. Romelli¹³, M. Roncarelli¹¹, E. Rossetti¹⁶, R. Saglia⁶⁰,⁵⁹, Z. Sakr⁹⁰,⁷,⁹¹, A. G. Sánchez⁵⁹, D. Sapone⁹², B. Sartoris⁶⁰,¹³, M. Schirmer⁶⁸, P. Schneider⁷⁶, T. Schrabback⁹³, A. Secroun⁵⁵, G. Seidel⁶⁸, M. Seiffert⁶⁴, S. Serrano⁹⁴,⁹⁵,⁹⁶, C. Sirignano⁸⁹,⁵¹, G. Sirri¹⁷, L. Stanco⁵¹, J. Steinwagner⁵⁹, P. Tallada-Crespi³⁴,³⁵, A. N. Taylor⁴, I. Tereno⁴⁸,⁹⁷, R. Toledo-Moreo⁹⁸, F. Torradeflot³⁵,³⁴, L. Valenziano¹¹,⁹⁹, T. Vassallo⁶⁰,¹³, Y. Wang¹⁰⁰, J. Weller⁶⁰,⁵⁹, G. Zamorani¹¹, E. Zucca¹¹, A. Biviano¹³,¹², M. Bolzonella¹¹, A. Boucaud⁷⁹, E. Bozzo⁵⁰, C. Burigana¹⁰¹,⁹⁹, M. Calabrese¹⁰²,³¹, D. Di Ferdinando¹⁷, J. A. Escartin Vigo⁵⁹, F. Finelli¹¹,⁹⁹, J. Gracia-Carpio⁵⁹, S. Matthew⁴, N. Mauri³⁸,¹⁷, A. Pezzotta⁵⁹, M. Pöntinen⁷², V. Scottez⁸¹,¹⁰³, A. Spurio Mancini¹⁰⁴,⁸, M. Tenti¹⁷, M. Viel¹²,¹³,¹⁵,¹⁴,¹⁵⁰, M. Wiesmann⁶², Y. Akrami¹⁰⁶,¹⁰⁷, S. Anselmi⁵¹,⁸⁹,¹⁰⁸, M. Archidiacono⁶¹,¹⁰⁹, F. Atrio-Barandela¹¹⁰, A. Balaguera-Antolinez³⁹,¹¹¹, M. Ballardini¹¹²,¹¹,¹¹³, D. Benielli⁵⁵, A. Blanchard⁷, L. Blot¹¹⁴,¹⁰⁸, H. Böhringer⁵⁹,¹¹⁵,¹¹⁶, S. Borgani¹¹⁷,¹²,¹³,¹⁴, S. Bruton¹¹⁸, R. Cabanac⁷, A. Calabro³², B. Camacho Quevedo⁹⁴,⁹⁶, A. Cappi¹¹,¹¹⁹, F. Caro³², C. S. Carvalho⁹⁷, T. Castro¹³,¹⁴,¹²,¹⁰⁵, K. C. Chambers¹²⁰, A. R. Cooray¹²¹, S. de la Torre²⁸, G. Desprez¹²², A. Díaz-Sánchez¹²³, S. Di Domizio²¹,²², H. Dole¹²⁴, S. Escoffier⁵⁵, A. G. Ferrari¹⁷, P. G. Ferreira¹²⁵, I. Ferrero⁶², A. Finoguenov⁷², A. Fontana³², F. Fornari⁹⁹, L. Gabarra¹²⁵, K. Ganga⁷⁹, J. García-Bellido¹⁰⁶, T. Gasparetto¹³, E. Gaztanaga⁹⁶,⁹⁴,³⁷, F. Giacomini¹⁷, F. Gianotti¹¹, G. Gozalias¹²⁶, C. M. Gutierrez¹²⁷,¹²⁸, W. G. Hartley⁵⁰, H. Hildebrandt¹²⁹, J. Hjorth⁸³, A. Jimenez Muñoz¹³⁰, S. Joudaki³⁷, J. J. E. Kajava¹³¹,¹³², V. Kansal¹³³,¹³⁴, D. Karagiannis¹³⁵,¹³⁶, C. C. Kirkpatrick⁷⁰, S. Kruk⁴², F. Lacasa¹³⁷,¹²⁴, M. Lattanzi¹¹³, A. M. C. Le Brun¹⁰⁸, J. Le Graet⁵⁵, L. Legrand¹³⁸, J. Lesgourgues³⁶, T. I. Liaudat¹³⁹, J. Macias-Perez¹³⁰, M. Magliocchetti⁵², F. Mannucci¹⁴⁰, R. Maoli¹⁴¹,³², J. Martín-Fleitas¹⁴², C. J. A. P. Martins¹⁴³,²⁶, L. Maurin¹²⁴, R. B. Metcalf⁷⁷,¹¹, M. Miluzio⁴²,¹⁴⁴, P. Monaco¹¹⁷,¹³,¹⁴,¹², A. Montoro⁹⁶,⁹⁴, C. Moretti¹⁵,¹⁰⁵,¹³,¹²,¹⁴, G. Morgante¹¹, C. Murray⁷⁹, S. Nadathur³⁷, N. A. Walton¹⁴⁵, L. Patrizii¹⁷, V. Popa⁸⁸, D. Potter¹⁴⁶, P. Reimberg⁸¹, I. Risso¹⁴⁷,¹⁴⁸, P.-F. Rocci¹²⁴, R. P. Rollins⁴, M. Sahlén¹⁴⁹, E. Sarpa¹⁵,¹⁰⁵,¹⁴, A. Schneider¹⁴⁶, M. Sereno¹¹,¹⁷, P. Simon⁷⁶, K. Tanidis¹²⁵, C. Tao⁵⁵, G. Testera²², R. Teyssier¹⁵⁰, S. Toft⁶⁷,¹⁵¹, S. Tosi²¹,²², A. Troja⁸⁹,⁵¹, M. Tucci⁵⁰, C. Valieri¹⁷, J. Valiviita⁷²,⁷³, D. Vergani¹¹, G. Verza¹⁵²,¹⁵³, P. Vielzeuf⁵⁵, M. L. Brown⁴⁰, and E. Sellentin¹⁵⁴,⁶³

*(Affiliations can be found after the references)*

Received 27 August 2024 / Accepted 8 November 2024

<footnote>* Corresponding author; n.tessore@ucl.ac.uk</footnote>
ABSTRACT

In this paper we present the framework for measuring angular power spectra in the *Euclid* mission. The observables in galaxy surveys, such as galaxy clustering and cosmic shear, are not continuous fields, but discrete sets of data, obtained only at the positions of galaxies. We show how to compute the angular power spectra of such discrete data sets, without treating observations as maps of an underlying continuous field that is overlaid with a noise component. This formalism allows us to compute the exact theoretical expectations for our measured spectra, under a number of assumptions that we track explicitly. In particular, we obtain exact expressions for the additive biases ('shot noise') in angular galaxy clustering and cosmic shear. For efficient practical computations, we introduce a spin-weighted spherical convolution with a well-defined convolution theorem, which allows us to apply exact theoretical predictions to finite-resolution maps, including *HEALPix*. When validating our methodology, we find that our measurements are biased by less than 1% of their statistical uncertainty in simulations of *Euclid*'s first data release.

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

1. Introduction

The photometric survey of the *Euclid* mission (Laureijs et al.
2011; Euclid Collaboration: Mellier et al. 2025) will infer cos-
mology using correlations between the observed angular posi-
tions of galaxies (angular galaxy clustering), their observed
shapes (cosmic shear), and the cross-correlation between posi-
tions 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 sur-
veys 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 mea-
sured in terms of real angular separation on the sky. Con-
versely, in harmonic-space methods, observations first undergo
a spherical harmonic transform before two-point statistics are
extracted. Examples of real-space methods include angular cor-
relation functions (Peebles 1973; Schneider et al. 2002), COSE-
BIs (Schneider et al. 2010), and band powers (Schneider et al.
2002), while examples of harmonic-space methods include var-
ious 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 effec-
tively slightly different probes of the same underlying informa-
tion. 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 angu-
lar 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 observ-
ables of the CMB are continuous temperature and polarisation
fields, for which maps are created by carefully planned obser-
vations that are slightly oversampled with respect to the instru-
ment’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’ con-
tribution 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 ambi-
tious science goals (Euclid Collaboration: Mellier et al. 2025),
we expect about half of the observed pixels in the result-
ing maps to be empty. Our approach for *Euclid* is there-
fore 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 with-
out 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 imple-
mentation 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 subse-
quently by Wolz et al. (2024).

In addition to practical computation, the discrete angular
power spectra offer an additional advantage on the theoreti-
cal side: not having to assume the existence of intermediary
maps with resolution-dependent noise greatly simplifies theoret-
ical predictions for the measured spectra. Apart from a number
of scientific assumptions, which we track and call out explic-
itly, this approach allows us to obtain an exact theory for the
expectations of our measurements. For angular galaxy cluster-
ing, we find that the additive shot noise bias is not random, but
is a known number, and we obtain an expected galaxy cluster-
ing signal that depends directly on the angular correlation func-
tion w(θ), as originally defined by Peebles (1973), instead of the
two-point statistics of an ancillary galaxy density field. For cos-
mic 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 cos-
mic shear field. In light of the stringent requirements on admis-
sible 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<fnref n="1" />. 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 $\hat{u}$ and $\hat{u}'$. A spherical function $f$ has spin weight $s$ if the function value $f(\hat{u})$ transforms under a rotation $\gamma$ of the coordinate frame in $\hat{u}$ as

$$f(\hat{u}) \mapsto e^{-i\gamma s} f(\hat{u}). \qquad (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 $s Y_{lm}$:

$$f(\hat{u}) = \sum_{lm} f_{lm} s Y_{lm}(\hat{u}). \qquad (2)$$

Here and in the following, sums always extend over all admissible values $l \ge |s|$ and $-l \le m \le l$. The coefficients $f_{lm}$ of the expansion are obtained by integration against the spherical harmonics,

$$f_{lm} = \int f(\hat{u}) [s Y_{lm}(\hat{u})]^* d\hat{u}, \qquad (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} \equiv {}_0Y_{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'}(\theta)$ as the correlation of $f(\hat{u})$ and $f'(\hat{u}')$ over all points $\hat{u}, \hat{u}'$ on the sphere separated by the angle $\theta$,

$$C^{ff'}(\theta) = \frac{1}{8\pi^2} \iint_{\hat{u}, \hat{u}' = \cos\theta} e^{i\alpha\alpha'} f^*(\hat{u}) f'(\hat{u}') e^{-i\alpha'\alpha} d\hat{u} d\hat{u}'. \qquad (4)$$

Here, if the spin weights $s$ and $s'$ are non-zero, the angles $\alpha$ and $\alpha'$, defined in Appendix C, describe a rotation of the respective coordinate frames in $\hat{u}$ and $\hat{u}'$ 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),

$$C^{ff'}(\theta) = \sum_l \frac{2l + 1}{4\pi} C_l^{ff'} d_{ss'}^l(\theta). \qquad (5)$$

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

$$C_l^{ff'} = 2\pi \int_0^\pi C^{ff'}(\theta) d_{ss'}^l(\theta) \sin(\theta) d\theta. \qquad (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'}(\theta)$ in the angular power spectrum (6) by its definition (4),

$$C_l^{ff'} = \frac{1}{4\pi} \iint f^*(\hat{u}) f'(\hat{u}') e^{i\alpha\alpha'} d_{ss'}^l(\theta) e^{-i\alpha'\alpha} d\hat{u} d\hat{u}', \qquad (7)$$

<footnote n="1">¹ https://github.com/heracles-ec/heracles</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Overview of the methodology presented below. We apply the formalism for discrete angular power spectra in three distinct ways: (i) Exact spherical harmonic coefficients can be computed from the discrete data, without the use of maps. In turn, angular power spectra can be computed from combinations of spherical harmonic coefficients. (ii) The angular power spectra themselves can be computed from the discrete data. This is inefficient for practical computation, but makes it possible to obtain exact expressions for the expected spectra. (iii) The discrete data can be turned into maps, and subsequently into spherical harmonic coefficients by means of a spherical harmonic transform. This can yield the same results as the discrete transformation.</figcaption></figure>

where the angles $\alpha, \theta, \alpha'$ still depend on $\hat{u}$ and $\hat{u}'$, but we now have two unrestricted integrals over the sphere. Inserting the spherical harmonic addition theorem

$$ \sum_m s Y_{lm}(\hat{u}) s' Y_{lm}^*(\hat{u}') = \frac{2l+1}{4\pi} e^{is\alpha} d_{ss'}^l(\theta) e^{-is'\alpha'} \quad (8) $$

into the integrand in definition (7), the two integrals decouple, and reduce to the spherical harmonic coefficients $f_{lm}$ and $f'_{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,

$$ C_l^{ff'} = \frac{1}{2l+1} \sum_m f_{lm}^* f_{lm'} \quad (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 $\theta$. 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<fnref n="2" />.

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

$$ f_{lm} \mapsto \sum_{\mu} f_{l\mu} D_{\mu m}^{l}(R), \quad (10) $$

where $D_{\mu m}^l$ is the Wigner $D$ function. The importance of homogeneous random fields on the sphere is closely related to this

<footnote n="2">² Generally speaking, homogeneity is invariance under translations, whereas invariance under rotations is isotropy. The translations of the sphere are rotations, and any rotation of the sphere can be achieved by three translations, so that homogeneity and isotropy are equivalent.</footnote>
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_{lm}^* f_{l'm'}^*$ of modes from the spherical harmonic expansion (2) of a pair of functions $f$ and $f'$. Under the rotation (10), the product transforms as

$$f_{lm}^* f_{l'm'}^* \mapsto \sum_{\mu\mu'} f_{l\mu}^* f_{l'\mu'}^* D_{\mu m}^{l\mu}(R) D_{\mu'm'}^{l'm'}(R). \quad (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 $\langle \cdot \rangle$, we find

$$\langle f_{lm}^* f_{l'm'}^* \rangle = \sum_{\mu\mu'} \langle f_{l\mu}^* f_{l'\mu'}^* \rangle D_{\mu m}^{l\mu}(R) D_{\mu'm'}^{l'm'}(R). \quad (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

$$\langle f_{lm}^* f_{l'm'}^* \rangle = \delta_{ll'}^K \delta_{mm'}^K \langle C_l^{ff'} \rangle, \quad (13)$$

where $\delta_{ll'}^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 $\langle C_l \rangle$ of random fields, in which case the expression (9) is used as an estimator for $\langle C_l \rangle$. However, we prefer to think of the sum (9) as the actual, realised, observable angular power spectrum, and $\langle C_l \rangle$ as its expectation over realisations.

Having obtained the two-point expectation (13) in harmonic space, its equivalent $\langle f^*(\hat{u}) f'(\hat{u}') \rangle$ 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):

$$\langle f^*(\hat{u}) f'(\hat{u}') \rangle = \frac{2l+1}{4\pi} \sum_l \langle C_l^{ff'} \rangle e^{-is\alpha} d_{ss'}^l(\theta) e^{is'\alpha'}. \quad (14)$$

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

$$\langle C_{l}^{ff'}(\theta) \rangle = \frac{2l+1}{4\pi} \sum_{l'} \langle C_{l}^{ff'} \rangle d_{ss'}^{l}(\theta), \quad (15)$$

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

$$\langle f^*(\hat{u}) f'(\hat{u}') \rangle = e^{-is\alpha} \langle C^{ff'}(\theta) \rangle e^{is'\alpha'}. \quad (16)$$

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

$$\langle C_l^{ff'} \rangle = 2\pi \int_0^\pi \langle C^{ff'}(\theta) \rangle d_{ss'}^l(\theta) \sin(\theta) d\theta, \quad (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$,

$$f(\hat{u}) = g(\hat{u}) w(\hat{u}). \quad (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_1, s_2$ are the respective spin weights of $f, g, w$, it follows that $s = s_1 + s_2$ 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'$,

$$C_{l'm'}^{ff'}(\theta) = \frac{1}{8\pi^2} \iint_{\hat{u}, \hat{u}' = \cos\theta} \left[ e^{is_1\alpha} g^*(\hat{u}) g'(\hat{u}') e^{-is_1\alpha'} \right] \times \left[ e^{is_2\alpha} w^*(\hat{u}) w'(\hat{u}') e^{-is_2\alpha'} \right] d\hat{u} d\hat{u}'. \quad (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'$,

$$e^{is_1\alpha} \langle g^*(\hat{u}) g'(\hat{u}') \rangle e^{-is_1\alpha'} = \langle C^{gg'}(\theta) \rangle. \quad (20)$$

We can factor $\langle C^{gg'}(\theta) \rangle$ out of the integral, which reduces to the angular correlation function (4) of $w$ and $w'$,

$$\langle C_{l'm'}^{ff'}(\theta) \rangle = \langle C^{gg'}(\theta) \rangle C^{ww'}(\theta). \quad (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),

$$\langle C_l^{ff'} \rangle = 2\pi \int_0^\pi \langle C^{gg'}(\theta) \rangle C^{ww'}(\theta) d_{ss'}^l(\theta) \sin(\theta) d\theta. \quad (22)$$

We then expand the angular correlations $\langle C^{gg'}(\theta) \rangle$ and $\langle C^{ww'}(\theta) \rangle$ in the integral using relation (5). Since $s = s_1 + s_2$ and $s' = s'_1 + s'_2$, we obtain Gaunt's integral for the $d$ functions (Edmonds 1960),

$$\frac{1}{2} \int_0^\pi d_{s_1 s'_1}^{l_1}(\theta) d_{s_2 s'_2}^{l_2}(\theta) d_{s s'}^{l}(\theta) \sin(\theta) d\theta = (-1)^{s-s'} \begin{pmatrix} l_1 & l_2 & l \\ s_1 & s_2 & -s \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l \\ s'_1 & s'_2 & -s' \end{pmatrix}, \quad (23)$$

where the right-hand side contains the Wigner $3j$ symbols. The result expresses the expected angular power spectrum $\langle C_l^{ff'} \rangle$ in terms of the angular power spectra $\langle C_l^{gg'} \rangle$ and $\langle C_l^{ww'} \rangle$,

$$\langle C_l^{ff'} \rangle = (-1)^{s-s'} \sum_{l_1 l_2} \frac{(2l_1 + 1)(2l_2 + 1)}{4\pi} \langle C_{l_1}^{gg'} \rangle C_{l_2}^{ww'} \\ \times \begin{pmatrix} l_1 & l_2 & l \\ s_1 & s_2 & -s \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l \\ s'_1 & s'_2 & -s' \end{pmatrix}. \quad (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 $\langle C_{l}^{ff'} \rangle$ as a function of the expected angular power spectrum $\langle C_{l}^{gg'} \rangle$ of the underlying random fields. In that case, the convolution (24) can be separated into a linear operator containing the sum over $l_2$,

$$M_{ll_1}^{ww'} = (-1)^{s-s'} \sum_{l_2} \frac{(2l_1 + 1)(2l_2 + 1)}{4\pi} C_{l_2}^{ww'} \times \begin{pmatrix} l_1 & l_2 & l \\ s_1 & s_2 & -s \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l \\ s'_1 & s'_2 & -s' \end{pmatrix}, \quad (25)$$

which can subsequently be applied to any given $\langle C_{l_1}^{gg'} \rangle$,

$$\langle C_{l}^{ff'} \rangle = \sum_{l_1} M_{ll_1}^{ww'} \langle C_{l_1}^{gg'} \rangle. \quad (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_1 w_2$ is the product of spherical functions $w_1$ and $w_2$, then

$$M^{ww'} \neq M^{w_1 w'_1} M^{w_2 w'_2}, \quad (27)$$

except for special cases. The reason is a lack of homogeneity: the random field $w_2 g$ that yields the mixing matrix $\bar{M}^{w_2 w'_2}$ is no longer homogeneous, and the product $w_1 (w_2 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 $\hat{u}_k$, $k = 1, 2, \dots$, on the sphere, as well as a set of weights $w_k^3$. 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 $\delta^D$,

$$n(\hat{u}) = \sum_k w_k \delta^D(\hat{u} - \hat{u}_k), \quad (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,

$$n_{lm} = \sum_k w_k Y_{lm}^*(\hat{u}_k). \quad (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 $\hat{u}'_k$ 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

$$C_l^{nn'} = \frac{1}{4\pi} \sum_{kk'} w_k w_{k'} P_l(\cos \theta_{kk'}), \quad (30)$$

where $P_l = d_{00}^l$ is the Legendre polynomial, and $\theta_{kk'}$ is the angle between $\hat{u}_k$ and $\hat{u}_{k'}$. 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, \dots$, which are observed at points $\hat{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 $\delta^D$,

$$f(\hat{u}) = \sum_k g_k w_k \delta^D(\hat{u} - \hat{u}_k), \quad (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_1$ and $s_2$ of $g$ and $w$: if a rotation $\gamma$ of the sphere in $\hat{u}_k$ transforms $g_k$ into $e^{-i s_1 \gamma} g_k$ and $w_k$ into $e^{-i s_2 \gamma} w_k$, the function value $f(\hat{u}_k)$ transforms into $e^{-i s \gamma} f(\hat{u}_k)$.

<footnote n="3">³ In what follows, we always implicitly assume that points have spin-0 weights, since that is the only practically relevant case. However, our results generalise straightforwardly to the spin-weighted case. For unweighted observations, the weights are set to unity.</footnote>
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,

$$f_{lm} = \sum_k g_k w_k s Y_{lm}^*(\hat{u}_k). \qquad (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

$$C_l^{ff'} = \frac{1}{4\pi} \sum_{kk'} g_k^* g_{k'} w_k w_{k'} e^{i s \alpha_{kk'}} d_{ss'}^l (\theta_{kk'}) e^{-i s' \alpha_{kk'}}, \qquad (33)$$

where the angles $\alpha_{kk'}$, $\theta_{kk'}$, $\alpha'_{kk'}$ are defined for each pair of points $\hat{u}_k, \hat{u}_{k'}$ 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,

$$C_l^{nf'} = \frac{1}{4\pi} \sum_{kk'} g_{k'}^* w_k w_{k'} d_{0s'}^l (\theta_{kk'}) e^{-i s' \alpha_{kk'}}, \qquad (34)$$

where $s'$ is the spin weight of $f'$, and the angles $\theta_{kk'}$, $\alpha'_{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 \equiv 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 \neq k'$, meaning $\hat{u}_k$ and $\hat{u}_{k'}$ are not the same observed point) and over the set of degenerate pairs of identical points (denoted by $k = k'$, meaning $\hat{u}_k$ and $\hat{u}_{k'}$ are the same observed point),

$$C_l^{nn'} = \frac{1}{4\pi} \left\{ \sum_{k \neq k'} + \sum_{k = k'} \right\} w_k w_{k'} P_l(\cos \theta_{kk'}). \qquad (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 $\theta_{kk'} = 0$ for $k \neq 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

$$A^{nn'} = \frac{1}{4\pi} \sum_{k \neq k'} w_k w_{k'}. \qquad (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

$$C_l^{nn'} = \frac{1}{4\pi} \sum_{k \neq k'} w_k w_{k'} P_l(\cos \theta_{kk'}) + A^{nn'}, \qquad (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,

$$\langle C_l^{nn'} - A^{nn'} \rangle = \frac{1}{4\pi} \left\{ \sum_{k \neq k'} w_k w_{k'} P_l(\cos \theta_{kk'}) \right\}. \qquad (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 $\hat{u}_k$ and $\hat{u}_{k'}$ 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,

$$ \langle C_l^{nn'} - A^{nn'} \rangle = \frac{1}{4\pi} \left\langle \sum_{\text{all } k \neq k'} w_k w_{k'} P_l(\cos \theta_{kk'}) \right\rangle. \quad (39) $$

Naturally, we have no knowledge about the unobserved points $\hat{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 $\hat{u}_k$, by making

**Assumption 1.** There exist functions $v$ and $v'$ that describe the expected weight conditional on the observed position,

$$ \langle w_k | \hat{u}_k \rangle = v(\hat{u}_k), \quad (40) $$

and similarly $\langle w'_{k'} | \hat{u}_{k'} \rangle = v'(\hat{u}_{k'})$.

We call $v$ and $v'$ the (weighted) 'visibility' of the respective observation; for unit weights, the value $v(\hat{u})$ is a number between 0 and 1 that describes the a priori probability that a point sampled in a given position $\hat{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,

$$ \langle C_l^{nn'} - A^{nn'} \rangle = \frac{1}{4\pi} \left\langle \sum_{\text{all } k \neq k'} v(\hat{u}_k) v'(\hat{u}_{k'}) P_l(\cos \theta_{kk'}) \right\rangle. \quad (41) $$

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

**Assumption 2.** 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$^4$, and hence terms in the sum. Introducing functions $\bar{n}$ and $\bar{n}'$ with

$$ \bar{n}(\hat{u}) = \frac{N}{4\pi} v(\hat{u}), \quad (42) $$

and similarly $\bar{n}'(\hat{u}') = (4\pi)^{-1} N' v'(\hat{u}')$, the expectation (41) is

$$ \langle C_l^{nn'} - A^{nn'} \rangle = 4\pi \langle \bar{n}(\hat{u}) \bar{n}'(\hat{u}') P_l(\cos \theta) \rangle, \quad (43) $$

where $\hat{u}, \hat{u}'$ is a pair of random points, and $\theta$ is the angular separation between them. The functions $\bar{n}$ and $\bar{n}'$ 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(\hat{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.** 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),

$$ \langle n(\hat{u}) \rangle = \bar{n}(\hat{u}). \quad (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 $\bar{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,

$$ \langle \bar{n}(\hat{u}) \bar{n}'(\hat{u}') P_l(\cos \theta) \rangle = \iint \bar{n}(\hat{u}) \bar{n}'(\hat{u}') P_l(\cos \theta) p(\hat{u}, \hat{u}') d\hat{u} d\hat{u}', \quad (45) $$

with $p(\hat{u}, \hat{u}')$ $d\hat{u} d\hat{u}'$ the a priori probability of the point process to generate a pair of points in $d\hat{u} d\hat{u}'$. 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.** 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(\hat{u}, \hat{u}')$ in the integral (45) depends only on the angular distance $\theta$ between the pair of points $\hat{u}$ and $\hat{u}'$. It can be written as (Landy & Szalay 1993; Peebles 1973)

$$ p(\hat{u}, \hat{u}') = \frac{1 + w(\theta)}{(4\pi)^2}, \quad (46) $$

<footnote>$^4$ For simplicity, we use $NN'$ for the number of pairs here, while the true number of pairs might be slightly different, e.g. $N(N-1)$ for an auto-correlation. If the difference is significant, one can introduce a pair count correction factor.</footnote>
where $w$ is the expected angular correlation function of density fluctuations in the observed point processes, which describes the clustering of points<fnref n="5" />.

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

$$ \langle C_{l}^{nn'} - A^{nn'} \rangle = \frac{1}{4\pi} \iint \bar{n}(\hat{u}) \bar{n}'(\hat{u}') [1 + w(\theta)] P_{l}(\cos \theta) d\hat{u} d\hat{u}'. \quad (47) $$

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

$$ \langle C_{l}^{nn'} - A^{nn'} \rangle = 2\pi \int_{0}^{\pi} C^{\bar{n}\bar{n}'}(\theta) [1 + w(\theta)] P_{l}(\cos \theta) \sin(\theta) d\theta. \quad (48) $$

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

$$ \langle C_{l}^{nn'} - A^{nn'} \rangle = C_{l}^{\bar{n}\bar{n}'} + \sum_{l_1} M_{ll_1}^{\bar{n}\bar{n}'} w_{l_1}, \quad (49) $$

where $M_{ll_1}^{\bar{n}\bar{n}'}$ is the mixing matrix (25) due to the mean density functions $\bar{n}$ and $\bar{n}'$, and $w_l$ is the angular power spectrum of the point processes, obtained from the intrinsic angular correlation function $w$ 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 $w_l$. However, the signal is doubly contaminated when the mean densities $\bar{n}$ and $\bar{n}'$ contain systematic variations, by both the angular power spectrum $C_{l}^{\bar{n}\bar{n}'}$ and by the associated mixing matrix $M_{ll_1}^{\bar{n}\bar{n}'}$. To remove these contaminations, we can directly manipulate expression (49) until it yields an estimator for $w_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 $w_l$ in the expectation (49). We directly construct spherical functions $\delta$ and $\delta'$ for the number density contrast of the observed points,

$$ \delta(\hat{u}) = \frac{n(\hat{u}) - \bar{n}(\hat{u})}{\bar{n}_0}, \quad (50) $$

and equivalently for $\delta'(\hat{u}')$, where $\bar{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 $\bar{n}$; for the alternative case (see Appendix B). The density contrast (50) is hence a linear combination of the spherical functions $n$ and $\bar{n}$, and it follows that the angular power spectrum of $\delta$ and $\delta'$ is

$$ C_{l}^{\delta\delta'} = \frac{C_{l}^{nn'} - C_{l}^{\bar{n}\bar{n}'} - C_{l}^{\bar{n}\bar{n}'} + C_{l}^{\bar{n}\bar{n}'} \bar{n}_{0}\bar{n}_{0}'}{\bar{n}_{0}\bar{n}_{0}'} \quad (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 $\langle C_{l}^{\delta\delta'} \rangle$ of the angular power spectrum (51) is readily computed using expressions (42), (44), and (49),

$$ \langle C_{l}^{\delta\delta'} - A^{\delta\delta'} \rangle = \sum_{l_1} M_{ll_1}^{\nu\nu'} w_{l_1}, \quad (52) $$

where $M^{\nu\nu'}$ is the mixing matrix for the visibilities $v$ and $v'$, and $A^{\delta\delta'} = (\bar{n}_0\bar{n}_{0'})^{-1}A^{nn'}$ 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(\hat{u}_k)$ of $g$ in the observed points, we can use the defining property of the delta function to factor $g(\hat{u})$ out of the sum in definition (31),

$$ f(\hat{u}) = g(\hat{u}) \sum_{k} w_k \delta^{\mathrm{D}}(\hat{u} - \hat{u}_k). \quad (53) $$

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

$$ w(\hat{u}) = \sum_{k} w_k \delta^{\mathrm{D}}(\hat{u} - \hat{u}_k). \quad (54) $$

We can therefore write $f(\hat{u}) = g(\hat{u})w(\hat{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 $\hat{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'(\hat{u}') = g'(\hat{u}')w'(\hat{u}')$. To this end, we firstly require

**Assumption 5.** 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.** The distribution of observed values $g_k$ is conditional on the observed positions $\hat{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 = gw$ and $f' = g'w'$, while $w$ and $w'$ are considered fixed

<footnote n="5">⁵ The angular correlation function $w$ is not to be confused with the weight $w$.</footnote>
functions. We can hence use the mixing matrix formalism (26) to obtain the expected angular power spectrum of $f$ and $f'$,

$$ \langle C_{l}^{ff'} \rangle = \sum_{l_1} M_{ll_1}^{ww'} \langle C_{l_1}^{gg'} \rangle, \qquad (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 $\epsilon_k$, which probes the cosmic shear field through the effect of weak gravitational lensing on the intrinsic galaxy shapes (e.g. Bartelmann & Schneider 2001),

$$ \epsilon_k = \frac{\epsilon_k^i + g_k}{1 + g_k^* \epsilon_k^i}, \qquad (56) $$

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

$$ \langle \epsilon_k | g_k \rangle = g_k. \qquad (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,

$$ \langle |\epsilon_k|^2 \rangle = \langle |g_k|^2 \rangle + \left\langle |\epsilon_k^i|^2 \frac{(1 - |g_k|^2)^2}{1 - |g_k|^2 |\epsilon_k^i|^2} \right\rangle. \qquad (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.** 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,

$$ \langle C_{l}^{ff'} \rangle = \sum_{l_1} M_{ll_1}^{ww'} \langle C_{l_1}^{gg'} \rangle + A^{ff'}, \qquad (59) $$

where $A^{ff'}$ is the additive bias due to the noise variance $\sigma_{kk'}^2$ from degenerate pairs of identical objects (denoted as before by

$k \equiv k')^6,$

$$ A^{ff'} = \delta_{ss'}^{K} \frac{1}{4\pi} \sum_{k \equiv k'} w_k w_{k'} \sigma_{kk'}^2, \qquad (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 $\sigma_{kk'}^2$ due to shape noise for each observed value $g_k$ or $g_{k'}$. In that situation, we can construct an estimate $\mathcal{A}^{ff'}$ of the additive bias from the variance of the noisy observations (Nicola et al. 2021),

$$ \mathcal{A}^{ff'} = \delta_{ss'}^{K} \frac{1}{4\pi} \sum_{k \equiv k'} w_k w_{k'} \epsilon_k^* \epsilon_{k'}^*. \qquad (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,

$$ \langle \mathcal{A}^{ff'} \rangle = A^{ff'} + \delta_{ss'}^{K} \frac{1}{4\pi} \sum_{k \equiv k'} w_k w_{k'} \langle C^{gg'}(0) \rangle, \qquad (62) $$

where the expected zero-lag angular correlation $\langle C^{gg'}(0) \rangle$ is the intrinsic variance of the random fields $g$ and $g'$. Subtracting $\mathcal{A}^{ff'}$ from the measured angular power spectrum $C_{l}^{ff'}$ and taking the expectation using expressions (59) and (62), we obtain

$$ \langle C_{l}^{ff'} - \mathcal{A}^{ff'} \rangle = \sum_{l_1} M_{ll_1}^{ww'} \langle C_{l_1}^{gg'} \rangle - \delta_{ss'}^{K} \frac{1}{4\pi} \sum_{k \equiv k'} w_k w_{k'} \langle C^{gg'}(0) \rangle. \qquad (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 $\langle C^{gg'}(0) \rangle$ by a sum over the expected angular power spectrum,

$$ \langle C^{gg'}(0) \rangle = \delta_{s_1 s'_1}^{K} \sum_{l} \frac{2l + 1}{4\pi} \langle C_{l}^{gg'} \rangle, \qquad (64) $$

with $s_1$ and $s'_1$ the respective spin weights of $g$ and $g'$, as above. The expectation (63) is therefore equivalent to

$$ \langle C_{l}^{ff'} - \mathcal{A}^{ff'} \rangle = \sum_{l_1} M_{ll_1}^{ww'} \langle C_{l_1}^{gg'} \rangle, \qquad (65) $$

where we have introduced the reduced mixing matrix

$$ M_{ll_1}^{ww'} = M_{ll_1}^{ww'} - \delta_{ss'}^{K} \delta_{s_1 s'_1}^{K} \frac{2l_1 + 1}{4\pi} \frac{1}{4\pi} \sum_{k \equiv k'} w_k w_{k'}. \qquad (66) $$

In the expectation (65), the bias introduced by $\mathcal{A}^{ff'}$ is thus completely absorbed into $M_{ll_1}^{ww'}$.

As it turns out, the reduced mixing matrix has a much simpler interpretation than the definition (66) suggests. Consider the

<footnote n="6">6 Here, degenerate pairs refer to observations of the same random field value $g_k \equiv g_{k'}$. Apart from auto-correlations, such pairs also arise, e.g. for cosmic shear when one set of galaxies is observed with two different shape measurement methods, where it may be the case that $g_k \equiv g_{k'}$ but $\epsilon_k \neq \epsilon_{k'}$ and $w_k \neq w_{k'}$.</footnote>
angular power spectrum (33) for the pair of weight functions $w$ and $w'$ with respective spin weights $s_2$ and $s'_2$. Following expression (35), we split $C_l^{ww'}$ into contributions from true pairs of distinct points ($k \neq 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'$,

$$A^{ww'} = \delta_{s_2 s'_2}^{K} \frac{1}{4\pi} \sum_{k \equiv k'} w_k w'_{k'} \quad (67)$$

Since $s = s_1 + s_2$ and $s' = s'_1 + s_2$, we can substitute $A^{ww'}$ for the sum in expression (66),

$$M_{l_1}^{ww'} = M_{l_1}^{ww'} - \delta_{ss'}^{K} \delta_{s_1 s'_1}^{K} \frac{2l_1 + 1}{4\pi} A^{ww'} \quad (68)$$

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

$$\delta_{ss'}^{K} \delta_{s_1 s'_1}^{K} = \sum_{l_2 \sigma} (2l_2 + 1) \begin{pmatrix} l & l_1 & l_2 \\ -s & s_1 & \sigma \end{pmatrix} \begin{pmatrix} l & l_1 & l_2 \\ -s' & s'_1 & \sigma \end{pmatrix}. \quad (69)$$

Using the fact that $A^{ww'}$ vanishes unless $s_2 = s'_2$, an equivalent way to write expression (68) is therefore

$$M_{l_1}^{ww'} = M_{l_1}^{ww'} - \delta_{ss'}^{K} \sum_{l_2} \frac{(2l_1 + 1)(2l_2 + 1)}{4\pi} A^{ww'} \times \begin{pmatrix} l_1 & l_2 & l \\ s_1 & s_2 & -s \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l \\ s'_1 & s'_2 & -s' \end{pmatrix}. \quad (70)$$

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

$$M_{l_1}^{ww'} = (-1)^{s-s'} \sum_{l_2} \frac{(2l_1 + 1)(2l_2 + 1)}{4\pi} (C_{l_2}^{ww'} - \delta_{ss'}^{K} A^{ww'}) \times \begin{pmatrix} l_1 & l_2 & l \\ s_1 & s_2 & -s \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l \\ s'_1 & s'_2 & -s' \end{pmatrix}. \quad (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 $\sigma^2$ is a delta-like contribution to the expected angular correlation function of the random fields,

$$\langle C^{gg'}(\theta) \rangle \mapsto \langle C^{gg'}(\theta) \rangle + \sigma^2 \delta^D(\cos \theta - \cos 0). \quad (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

$$C_l^{ff'} - \mathcal{A}^{ff'} \iff C^{ff'}(\theta) - C^{ff'}(0) \delta^D(\cos \theta - \cos 0) \quad (73)$$

for the random fields, and

$$C_l^{ww'} - A^{ww'} \iff C^{ww'}(\theta) - C^{ww'}(0) \delta^D(\cos \theta - \cos 0) \quad (74)$$

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

$$\langle C^{ff'}(\theta) - C^{ff'}(0) \delta^D(\cos \theta - \cos 0) \rangle = [\langle C^{gg'}(\theta) \rangle + \sigma^2 \delta^D(\cos \theta - \cos 0)] \times [C^{ww'}(\theta) - C^{ww'}(0) \delta^D(\cos \theta - \cos 0)], \quad (75)$$

where we can evaluate the right-hand side for all $\theta \ge 0$ without knowing the value of $\sigma^2$.

### 4.3. Cross-correlations, galaxy–galaxy lensing

The final case of interest is the cross-correlation of points $\hat{u}_k$ generated by a point process, and observed values $g'_{k'} = g(\hat{u}_{k'})$ 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 $\hat{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(\hat{u}_k)$ under Assumption 1,

$$\langle C_l^{nf'} \rangle = \frac{1}{4\pi} \left\langle \sum_{\text{all } k} \sum_{k'} g'_{k'} v(\hat{u}_k) w'_{k'} d_{0s'}^l(\theta_{kk'}) e^{-is'\alpha'} \right\rangle. \quad (76)$$

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

**Assumption 8.** 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 $\bar{n}$. Using the definition (54) of the weight function $w'$, we may also replace the remaining sum over $k'$ by an integral over $\hat{u}'$,

$$\langle C_l^{nf'} \rangle = \int \langle g'(\hat{u}') \bar{n}(\hat{u}) w'(\hat{u}') d_{0s'}^l(\theta) e^{-is'\alpha'} \rangle d\hat{u}', \quad (77)$$
where the angles $\theta$ and $\alpha'$ now describe the relative orientation between the random point $\hat{u}$ and $\hat{u}'$.

Using Assumption 6, we can factor the weight $w'(\hat{u}')$ out of the integral in expectation (77). Furthermore, using Assumption 3, the expectation is conditional on the mean number density $\bar{n}$, and only the position $\hat{u}$ in $\bar{n}(\hat{u})$ is random. The remaining expectation in Eq. (77) therefore reduces to the random point $\hat{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 $\hat{u}$ is

$$ \langle g'(\hat{u}') e^{-is'\alpha'} | \hat{u} \rangle = \gamma(\theta), \quad (78) $$

where $\gamma$ is the expected angular cross-correlation function$^7$. For galaxy–galaxy lensing, the expected correlation is more commonly written in terms of a tangential component $\gamma_t$ and cross-component $\gamma_\times$ as $\gamma(\theta) = \gamma_t(\theta) + i \gamma_\times(\theta)$.

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

$$ \langle C_l^{nf'} \rangle = \frac{1}{4\pi} \iint \bar{n}(\hat{u}) w'(\hat{u}') \gamma(\theta) d_{0s'}^l(\theta) d\hat{u} d\hat{u}'. \quad (79) $$

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

$$ \langle C_l^{nf'} \rangle = 2\pi \int_0^{\pi} C^{\bar{n}w'}(\theta) \gamma(\theta) d_{0s'}^l(\theta) \sin(\theta) d\theta, \quad (80) $$

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

$$ \langle C_l^{nf'} \rangle = \sum_{l_1} M_{ll_1}^{\bar{n}w'} \gamma_{l_1}, \quad (81) $$

where $\gamma_l$ is the angular power spectrum associated with the angular cross-correlation function $\gamma$.

Overall, we therefore obtain the intuitively clear result that the expected angular power spectrum is given by the intrinsic spectrum $\gamma_l$ for point process and random field, modulated by a mixing matrix coming from the mean number density $\bar{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 $\delta$ instead of the number density $n$. Since the density contrast (50) is linear in $n$, the angular power spectrum of $\delta$ and $f'$ is

$$ C_l^{\delta f'} = \frac{C_l^{nf'} - C_l^{\bar{n}f'}}{\bar{n}_0}. \quad (82) $$

The expectation $\langle C_l^{\delta f'} \rangle$ then follows immediately from the definition (42) of the mean density and expectation (81),

$$ \langle C_l^{\delta f'} \rangle = \sum_{l_1} M_{ll_1}^{vv'} \gamma_{l_1} - \langle C_l^{vf'} \rangle = \sum_{l_1} M_{ll_1}^{vv'} \gamma_{l_1}, \quad (83) $$

where the second equality assumes that the expectation of $f'$ vanishes.

Whether to correlate $n$ or $\delta$ 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),

$$ \hat{\gamma}_l^n = \sum_{l_1} (M^{\bar{n}w'})^{-1}_{ll_1} C_l^{nf'}, \quad (84) $$

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

$$ \begin{aligned} \hat{\gamma}_l^\delta &= \sum_{l_1} (M^{vv'})^{-1}_{ll_1} \frac{C_l^{nf'} - C_l^{\bar{n}f'}}{\bar{n}_0} \\ &= \sum_{l_1} (M^{\bar{n}w'})^{-1}_{ll_1} C_l^{nf'} - \sum_{l_1} (M^{\bar{n}w'})^{-1}_{ll_1} C_l^{\bar{n}f'}. \end{aligned} \quad (85) $$

In real space, the mean number density $\bar{n}$ corresponds to an equivalent distribution of uniform random points ('randoms'); the inverse mixing matrix in both $\hat{\gamma}_l^n$ and $\hat{\gamma}_l^\delta$ thus corresponds to a normalisation by weighted pairs of randoms and the observed positions of the field. The difference between $\hat{\gamma}_l^n$ and $\hat{\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 $\delta$ 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^2l_{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}^2)$, and it is hence favourable when $l_{max} \ll 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

<footnote n="7">⁷ Not to be confused with the shear in gravitational lensing.</footnote>
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:

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

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

$$F(\hat{u}) = \int f(\hat{u}') e^{iS\alpha} K(\theta) e^{-iS\alpha'} d\hat{u}', \quad (86)$$

where the angles $\theta, \alpha, \alpha'$ are the separation and relative orientation of $\hat{u}$ and $\hat{u}'$ 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 $\gamma$ in $\hat{u}$, the angle $\alpha$ in the convolution (86) transforms as $\alpha \mapsto \alpha - \gamma$, so that $F$ indeed picks up the phase factor $e^{-iS\gamma}$ 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,

$$K(\theta) = \sum_{l} \frac{2l + 1}{4\pi} K_l d_{Ss}^{l}(\theta), \quad (87)$$

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

$$K_l = 2\pi \int_0^{\pi} K(\theta) d_{Ss}^l(\theta) \sin(\theta) d\theta. \quad (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,

$$F_{lm} = K_l f_{lm}, \quad (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\alpha}$ and $e^{-iS\alpha'}$ 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 $\alpha \approx \alpha'$. 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 $\beta > 0$. However, for $s \neq 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 \ge 0$, and reduces to a spherical disc if $s = 0$,

$$K(\theta) = \begin{cases} [\cos \frac{\theta}{2}]^{2s} & \text{if } \theta \le \beta, \\ 0 & \text{otherwise,} \end{cases} \quad (90)$$

where $\beta$ is the angular size (i.e. radius) of the kernel. The effective area of the convolution kernel is<fnref n="8" />

$$\Omega_K = 2\pi \int_0^{\pi} K(\theta) \sin(\theta) d\theta = 4\pi \frac{1 - [\cos \frac{\beta}{2}]^{2s+2}}{s+1}. \quad (91)$$

The choice of kernel (90) is firstly motivated by the fact that there is a known expression for its coefficients (88)<fnref n="9" />,

$$K_l = \begin{cases} \frac{4\pi \left[ \sin \frac{\beta}{2} \right] \left[ \cos \frac{\beta}{2} \right]^{2s+1} d_{ss+1}^{l}(\beta)}{\sqrt{l(l+1)-s(s+1)}} & \text{if } l > s, \\ \frac{4\pi \left[ 1 - \left[ \cos \frac{\beta}{2} \right]^{4s+2} \right]}{2s+1} & \text{if } l = s. \end{cases} \quad (92)$$

Secondly, for small kernel sizes $\beta \le 1^\circ$, 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

<footnote n="8">⁸ Not to be confused with the cosmological curvature matter density.</footnote>

<footnote n="9">⁹ The convolution kernel (90) and coefficients (92) follow from the integral (4.11.9) of Varshalovich et al. (1988), which in fact yields a more general, spin-changing convolution kernel with $S \neq s$.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2. Convolution kernel (90) with no cut-off for spin weight $s = 0$ (black) and spin weight $s = 2$ (red). For sufficiently small angles the convolution kernel becomes indistinguishable from a flat spherical disc even in the spin-weighted case.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Normalised coefficients of the convolution kernel (90) for spin weights $s = 0$ (black), $s = 2$ (red), and kernel sizes $\beta = 10^{\circ}$ (top), $\beta = 1^{\circ}$ (middle), $\beta = 1$ arcmin (bottom). There is excellent agreement between the coefficients for $s = 2$ and $s = 0$, except in the case of a large kernel ($\beta = 10^{\circ}$) at large angular scales ($l < 10$).</figcaption></figure>

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(\theta)$ 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 $\beta \ll \pi$, an alternative method is to reverse the order of operations, and find all grid points closer than $\beta$ 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),

$$f_{lm} = \frac{1}{K_l} F_{lm}. \qquad (93)$$

Naturally, this is only possible when $K_l \neq 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 \le 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 $\Omega_K^{-1} 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),

$$C_l^{ff'} = K_l^2 C_l^{ff'}. \qquad (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^2$. 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).
<figure><img src="image_6.png" /><figcaption>Fig. 4. HEALPix pixel window function (red) for resolution parameters NSIDE = 256 (top), NSIDE = 1024 (middle), NSIDE = 4096 (bottom). Also shown is the convolution kernel of a spherical disc with the same pixel area (black). The HEALPix pixel window function is only provided up to $l = 4$ NSIDE, where it starts to fall below the kernel of the spherical disc.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Angular power spectrum of 1000 000 uniform random points, computed using the HEALPix pseudo-convolution (red) and a true convolution with spherical discs of the same area (black). Both results agree with their respective expectation (dashed). The convolution kernel is only imprinted on the true convolution.</figcaption></figure>

### 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 $\Omega_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 $\tilde{n}$ is the map $\tilde{N}$ of mean number counts; we write it as $\tilde{N} = \tilde{N}_0 V$ using the mean number of points per pixel $\tilde{N}_0 = \Omega_K \tilde{n}_0$ and a map $V$ that is the convolution of $v$ divided by $\Omega_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, $\langle N \rangle = \tilde{N}$.

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

$$ \Delta = \frac{N - \tilde{N}}{\tilde{N}_0}. \qquad (95) $$

It is clear that $\Delta$ is the convolution of $\delta$, but divided by the pixel area, so that the numerical values of $\Delta$ have the correct, intuitive scale where $-1$ means “empty space”. Deconvolution of $\Delta$ must therefore be carried out with the normalised convolution kernel $\Omega_K^{-1} 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 $\tilde{W}_0$. Specifically, we choose $\tilde{W}_0$ to be the mean weight divided by the mean visibility, as computed from the maps, which makes $\tilde{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 $\tilde{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 $\tilde{W}_0$.

Since $\tilde{W}_0$ contains a factor of the pixel area, deconvolution of $F$ and $W$ is carried out with the normalised convolution kernel

$\Omega_K^{-1} K_l$. Furthermore, since both $F$ and $W$ are scaled identically, the resulting mixing matrix automatically applies the correct factors of $\tilde{W}_0$ to the expected angular power spectra. However, we do need to account for the scaling by $\tilde{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<fnref n="10" />.

<footnote n="10">¹⁰ https://gitlab.mpcdf.mpg.de/mtr/ducc</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Simulated distribution of the relative error when estimating the mean density $\bar{n}$ from part of the sky, using the footprint of Euclid DR1 (black), DR2 (red), and DR3 (blue) for redshifts $z = 0.5$ (top), $z = 1.0$ (middle), and $z = 2.0$ (bottom).</figcaption></figure>

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 $\Lambda$CDM cosmology as Euclid's Flagship simulation (Euclid Collaboration: Castander et al. 2025), with parameter values $\Omega_m = 0.319$, $\Omega_b = 0.049$, $A_s = 2.1 \times 10^{-9}$, $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_{\hat{l}}^{\hat{\delta}\hat{\delta}}$ 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

<figure><img src="image_9.png" /><figcaption>Fig. 7. Argument of the phase factor bias $e^{-2i(\alpha-\alpha')}$ for HEALPix maps with NSIDE = 1, 4, 16, 64, shown in orthographic projection with the north pole at the centre and the equator at the border. The angles $\alpha$ and $\alpha'$ are defined in Appendix C. For NSIDE = 64, the resolution is at the degree-scale, and the phase factors are close to unity.</figcaption></figure>

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 $\hat{u}_0$. By neglecting the phase factors, the convolution (86) is approximated as

$$F(\hat{u}_0) \approx \int f(\hat{u}') K(\theta) d\hat{u}'. \qquad (96)$$

To cancel the phase factors, the approximation effectively applies a position-dependent multiplicative bias $e^{-i(\alpha-\alpha')}$ to the function $f$ over the pixel area, where the angles $\alpha$ and $\alpha'$ are taken with respect to $\hat{u}_0$. Under a rotation of $\gamma'$ in $\hat{u}'$, the angle $\alpha'$ transforms as $\alpha' \mapsto \alpha' - \gamma'$; as a function of $\hat{u}'$, 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 $\ge$ 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
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Bias relative to statistical uncertainty between measured and expected angular power spectra from 10 000 simulations with a Euclid DR1-like set-up. Shown are spectra for combinations of galaxy density $\delta$ and cosmic shear E- and B-mode in two tomographic bins. The discrete angular power spectra show very good agreement between measurement and expectation (black), except for effects that can be ascribed to source-lens clustering. The same effects are visible for map-based spectra from HEALPix; in addition, these also show a small residual in angular clustering due to the pseudo-convolution of HEALPix pixels (red). When source-lens clustering is taken into account (blue), the relative bias of the discrete spectra is consistent with zero at the 1% level (error bars). Points of the three data sets are slightly offset for better visibility.</figcaption></figure>

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$^2$, 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 interme-
diary 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 $\sigma_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$^2$ 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 $\sigma_e = 0.26$. Instead of using the correct weak lensing action (56), we simply sum the intrinsic ellipticity and the shear $\gamma$ 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^{-2}$ and $10^2$, 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 $\delta$, 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_1 = 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 $\Delta C_l = C_l - \langle C_l \rangle$ between measurements and expectations, which we scale relative to the standard deviation $\sigma_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 \gtrsim 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 con-
nect our exact theory with the standard practice of measuring
angular power spectra from *HEALPix* maps, for example. Con-
ceptually, 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 for-
malism and standard practice is that observations are summed
over one map pixel (i.e. convolution kernel), but are not aver-
aged. This side-steps common issues with map-based spectra
(e.g. empty pixels) or the fact that pixels containing just one sin-
gle 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<fnref n="1" />. 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.

Acknowledgements. 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.11), 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 Nederlandse Onderzoeksg school 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 (www.euclid-ec.org).

References

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

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

Alsing, J., Heavens, A., Jaffe, A. H., et al. 2016, MNRAS, 455, 4452

Baleato Lizancos, A., & White, M. 2024, JCAP, 05, 010

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

Boyle, M. 2016, J. Math. Phys., 57, 092504

Brown, M. L., Castro, P. G., & Taylor, A. N. 2005, MNRAS, 360, 1262

Carretero, J., Tallada, P., Casals, J., et al. 2017, Proceedings of the European Physical Society Conference on High Energy Physics, 488

Challinor, A., & Lewis, A. 2011, Phys. Rev. D, 84, 043516

Deshpande, A. C., Kitching, T. D., Cardone, V. F., et al. 2020, A&A, 636, A95

Driscoll, J. R., & Healy, D. M. 1994, Adv. Appl. Math., 15, 202

Dupac, X., & Tauber, J. 2005, A&A, 430, 363

Edmonds, A. R. 1960, Angular Momentum in Quantum Mechanics (Princeton: Princeton University Press)

Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191

Euclid Collaboration (Lepori, F., et al.) 2022, A&A, 662, A93

Euclid Collaboration (Scaramella, R., et al.) 2022, A&A, 662, A112

Euclid Collaboration (Mellier, Y., et al.) 2025, A&A, in press, https://doi.org/10.1051/0004-6361/202450810

Euclid Collaboration (Castander, F. J., et al.) 2025, A&A, in press, https://doi.org/10.1051/0004-6361/202450853

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

Heavens, A. F., & Taylor, A. N. 1995, MNRAS, 275, 483

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

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

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

Hilbert, S., Hartlap, J., & Schneider, P. 2011, A&A, 536, A85

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

Huffenberger, K. M., & Wandelt, B. D. 2010, ApJS, 189, 255

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

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

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

Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. D, 55, 7368

Kerscher, M., Szapudi, I., & Szalay, A. S. 2000, ApJ, 535, L13

Kiessling, A., Cacciato, M., Joachimi, B., et al. 2015, Space Sci. Rev., 193, 67

Kirk, D., Brown, M. L., Hoekstra, H., et al. 2015, Space Sci. Rev., 193, 139

Kitching, T. D., & Deshpande, A. C. 2022, Open J. Astrophys., 5, 6

Kitching, T., Deshpande, A., & Taylor, P. 2021, Open J. Astrophys., 4, 17

Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193]

Linke, L., Unruh, S., Wittje, A., et al. 2024, ArXiv e-prints [arXiv:2407.09810]

Loureiro, A., Whiteway, L., Sellentin, E., et al. 2023, Open J. Astrophys., 6, 6

Maraio, A., Hall, A., & Taylor, A. 2023, MNRAS, 520, 4836

McEwen, J. D., & Wiaux, Y. 2011, IEEE Trans. Signal Process., 59, 5876

McEwen, J. D., Leistedt, B., Büttner, M., Peiris, H. V., & Wiaux, Y. 2015, ArXiv e-prints [arXiv:1509.06749]

Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401

More, S., Sugiyama, S., Miyatake, H., et al. 2023, Phys. Rev. D, 108, 123520

Nicola, A., García-García, C., Alonso, D., et al. 2021, JCAP, 03, 067

Peebles, P. J. E. 1973, ApJ, 185, 413

Percival, W. J., Burkey, D., Heavens, A., et al. 2004, MNRAS, 353, 1201

Reinecke, M., Belkner, S., & Carron, J. 2023, A&A, 678, A165

Rodríguez-Monroy, M., Weaverdyck, N., Elvin-Poole, J., et al. 2022, MNRAS, 511, 2665

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

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

Seitz, C., & Schneider, P. 1997, A&A, 318, 687

Sellentin, E., Loureiro, A., Whiteway, L., et al. 2023, Open J. Astrophys., 6, 31

Singh, S., Mandelbaum, R., Seljak, U., Slosar, A., & Vazquez Gonzalez, J. 2017, MNRAS, 471, 3827

Tadros, H., Ballinger, W. E., Taylor, A. N., et al. 1999, MNRAS, 305, 527

Tallada, P., Carretero, J., Casals, J., et al. 2020, Astron. Comput., 32, 100391

Tegmark, M. 1997, Phys. Rev. D, 55, 5895

Tegmark, M., & de Oliveira-Costa, A. 2001, Phys. Rev. D, 64, 063001

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

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

Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1988, Quantum Theory of Angular Momentum (World Scientific)

Wandelt, B. D., & Górski, K. M. 2001, Phys. Rev. D, 63, 123002

Wandelt, B. D., Hivon, E., & Górski, K. M. 2001, Phys. Rev. D, 64, 083003

Wolz, K., Alonso, D., & Nicola, A. 2025, J. Cosmol. Astropart. Phys., 2025, 028

Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830

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

<footnote n="2">² Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm SE-106 91, Sweden</footnote>
³ Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK

⁴ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

⁵ European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

⁶ Institute Lorentz, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

⁷ Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France

⁸ Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

⁹ School of Mathematics and Physics, University of Surrey, Guildford, Surrey GU2 7XH, UK

¹⁰ INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy

¹¹ INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy

¹² IFPU, Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151 Trieste, Italy

¹³ INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy

¹⁴ INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste, TS, Italy

¹⁵ SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, TS, Italy

¹⁶ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy

¹⁷ INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

¹⁸ Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay, 91191 Gif-sur-Yvette Cedex, France

¹⁹ Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis Boulevard Arago, 75014 Paris, France

²⁰ INAF-Osservatorio Astrofísico di Torino, Via Osservatorio 20, 10025 Pino Torinese, TO, Italy

²¹ Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146 Genova, Italy

²² INFN-Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy

²³ Department of Physics “E. Pancini”, University Federico II, Via Cinthia 6, 80126 Napoli, Italy

²⁴ INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

²⁵ INFN Section of Naples, Via Cinthia 6, 80126 Napoli, Italy

²⁶ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal

²⁷ Faculdade de Ciências da Universidade do Porto, Rua do Campo de Alegre, 4150-007 Porto, Portugal

²⁸ Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

²⁹ Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy

³⁰ INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

³¹ INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy

³² INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzo Catone, Italy

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

³⁴ Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain

³⁵ Port d'Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra, (Barcelona), Spain

³⁶ Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany

³⁷ Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK

³⁸ Dipartimento di Fisica e Astronomia “Augusto Righi” - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

³⁹ Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, 38204 San Cristóbal de La Laguna, Tenerife, Spain

⁴⁰ Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK

⁴¹ European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy

⁴² ESAC/ESA, Camino Bajo del Castillo, s/n, Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain

⁴³ Université Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, Villeurbanne F-69100, France

⁴⁴ Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

⁴⁵ Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, 08028 Barcelona, Spain

⁴⁶ Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig de Lluís Companys 23, 08010 Barcelona, Spain

⁴⁷ UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 Rue Enrico Fermi, 69622 Villeurbanne, France

⁴⁸ Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande PT1749-016, Lisboa, Portugal

⁴⁹ Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal

⁵⁰ Department of Astronomy, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland

⁵¹ INFN-Padova, Via Marzolo 8, 35131 Padova, Italy

⁵² INAF-Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere, 100, 00100 Roma, Italy

⁵³ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

⁵⁴ Space Science Data Center, Italian Space Agency, Via del Politecnico snc, 00133 Roma, Italy

⁵⁵ Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

⁵⁶ Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irneio 46, 40126 Bologna, Italy

⁵⁷ FRACTAL S.L.N.E., Calle Tulipán 2, Portal 13 1A, 28231 Las Rozas de Madrid, Spain

⁵⁸ INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy

⁵⁹ Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany

⁶⁰ Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany

⁶¹ Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy

⁶² Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern 0315, Oslo, Norway

⁶³ Leiden Observatory, Leiden University, Einsteinweg 55, 2333 CC Leiden, The Netherlands

⁶⁴ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

⁶⁵ Felix Hormuth Engineering, Goethestr. 17, 69181 Leimen, Germany

⁶⁶ Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark

⁶⁷ Cosmic Dawn Center (DAWN), Denmark

⁶⁸ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

⁶⁹ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

⁷⁰ Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland

⁷¹ Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 Quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland

⁷² Department of Physics, PO Box 64, 00014 University of Helsinki, Finland
73 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland

74 NOVA Optical Infrared Instrumentation Group at ASTRON, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands

75 Centre de Calcul de l'IN2P3/CNRS, 21 Avenue Pierre de Coubertin, 69627 Villeurbanne Cedex, France

76 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany

77 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Via Piero Gobetti 93/2, 40129 Bologna, Italy

78 Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK

79 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

80 University of Applied Sciences and Arts of Northwestern Switzerland, School of Engineering, 5210 Windisch, Switzerland

81 Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France

82 Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra, (Barcelona), Spain

83 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 155, 2200 Copenhagen, Denmark

84 Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

85 Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

86 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada

87 Centre National d'Etudes Spatiales – Centre Spatial de Toulouse, 18 Avenue Edouard Belin, 31401 Toulouse Cedex 9, France

88 Institute of Space Science, Str. Atomistilor, Nr. 409 Măgurele, Ilfov 077125, Romania

89 Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, Via Marzolo 8, 35131 Padova, Italy

90 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany

91 Université St Joseph; Faculty of Sciences, Beirut, Lebanon

92 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile

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

94 Institut d'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain

95 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain

96 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain

97 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal

98 Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain

99 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy

100 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

101 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy

102 Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA), Loc. Lignan 39, I-11020 Nus, (Aosta Valley), Italy

103 Junia, EPA Department, 41 Bd Vauban, 59800 Lille, France

104 Department of Physics, Royal Holloway, University of London, London TW20 0EX, UK

105 ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy

106 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain

107 CERCA/ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA

108 Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France

109 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy

110 Departamento de Física Fundamental, Universidad de Salamanca, Plaza de la Merced s/n, 37008 Salamanca, Spain

111 Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain

112 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

113 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

114 Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

115 Ludwig-Maximilians-University, Schellingstrasse 4, 80799 Munich, Germany

116 Max-Planck-Institut für Physik, Boltzmannstr. 8, 85748 Garching, Germany

117 Dipartimento di Fisica – Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy

118 Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA

119 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France

120 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

121 Department of Physics & Astronomy, University of California Irvine, Irvine, CA 92697, USA

122 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia B3H 3C3, Canada

123 Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain

124 Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405, Orsay, France

125 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK

126 Department of Computer Science, Aalto University, PO Box 15400, Espoo FI-00 076, Finland

127 Instituto de Astrofísica de Canarias, c/ Via Lactea s/n, La Laguna E-38200, Spain

128 Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco Sanchez, La Laguna E-38200, Spain

129 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany

130 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000 Grenoble, France

131 Department of Physics and Astronomy, Vesilinnantie 5, 20014 University of Turku, Finland

132 Serco for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692, Madrid, Spain

133 ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia

134 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Hawthorn, Victoria 3122, Australia

135 School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK

136 Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa

137 Université Libre de Bruxelles (ULB), Service de Physique Théorique CP225, Boulevard du Triophe, 1050 Bruxelles, Belgium
138 ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil

139 IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette Cedex, France

140 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

141 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy

142 Aurora Technology for European Space Agency (ESA), Camino Bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692, Madrid, Spain

143 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal

144 HE Space for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692, Madrid, Spain

145 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

146 Department of Astrophysics, University of Zurich, Winterthur-erstrasse 190, 8057 Zurich, Switzerland

147 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy

148 INFN-Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy

149 Theoretical Astrophysics, Department of Physics and Astronomy, Uppsala University, Box 515, 751 20 Uppsala, Sweden

150 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA

151 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark

152 Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA

153 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010 New York, NY, USA

154 Mathematical Institute, University of Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands
Appendix A: Additional figures

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Angular power spectra (red) for angular clustering (upper triangle) and cosmic shear (lower triangle) in the Euclid Flagship simulation with a DR1-like footprint. For cosmic shear, the B-mode spectrum due to mode mixing is shown in blue. Also shown is the expectation for each spectrum (dashed), as computed from the respective mixing matrices. All spectra are binned into 32 angular bins with logarithmic spacing between $l = 10$ and $l = 5000$. The y-axis changes to linear scaling when passing through the origin.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. A.2. Same as Fig. A.1, but for galaxy–galaxy lensing.</figcaption></figure>
## 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 $\bar{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 $w_l$ from the expectation (49). For example, using a formal inverse $(M^{\bar{n}\bar{n}})^{-1}$ of the mixing matrix, we can construct the estimator

$$ \hat{w}_l^N = \sum_{l'} (M^{\bar{n}\bar{n}'})^{-1}_{ll'} [C_{l'}^{nn'} - C_{l'}^{\bar{n}\bar{n}'} - A^{nn'}]. \quad (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 $\langle C_{l}^{nn'} \rangle = C_{l}^{\bar{n}\bar{n}'}$, and we can hence construct a more advanced estimator

$$ \hat{w}_l^{LS} = \sum_{l'} (M^{\bar{n}\bar{n}'})^{-1}_{ll'} [C_{l'}^{nn'} - C_{l'}^{\bar{n}\bar{n}'} - C_{l'}^{\bar{n}\bar{n}'} + C_{l'}^{\bar{n}\bar{n}'} - A^{nn'}]. \quad (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,

$$ \tilde{w}_l^N = C_{l}^{nn'} - C_{l}^{\bar{n}\bar{n}'} - A^{nn'}, \quad (B.3) $$

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

$$ \tilde{w}_l^{LS} = C_{l}^{nn'} - C_{l}^{\bar{n}\bar{n}'} - C_{l}^{\bar{n}\bar{n}'} + C_{l}^{\bar{n}\bar{n}'} - A^{nn'}. \quad (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 $\bar{n}_0$, consider an alternative definition of the density contrast with an arbitrary normalisation function $q$,

$$ \delta_q(\hat{u}) = \frac{n(\hat{u}) - \bar{n}(\hat{u})}{\bar{n}_0 q(\hat{u})}. \quad (B.5) $$

It follows from the definitions of the number density (28) and mean number density (42) that $\delta_q$ is equivalent to $\delta$ under a change of weights $w_k \mapsto w_k/q(\hat{u}_k)$ and, consequently, a change of visibility $v(\hat{u}) \mapsto v(\hat{u})/q(\hat{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 $\alpha$ and $\alpha'$ that describe the relative orientation of points $\hat{u}$ and $\hat{u}'$ on the sphere, it suffices to solve the

<figure><img src="image_13.png" /><figcaption>Fig. C.1. Angles $\alpha, \theta, \alpha'$ that describe the relative orientation between two points $\hat{u}$ and $\hat{u}'$ on the sphere can be obtained from the spherical triangle between the north pole, $\hat{u}$, and $\hat{u}'$, with $\Delta\varphi = \varphi - \varphi'$.</figcaption></figure>

spherical triangle shown in Fig. C.1 (for more information, see Hall & Tessore in prep.),

$$ \alpha = \arctan \frac{\sin \vartheta' \sin(\varphi - \varphi')}{\sin \vartheta \cos \vartheta' - \cos \vartheta \sin \vartheta' \cos(\varphi - \varphi')}, \quad (C.1) $$

$$ \alpha' = -\arctan \frac{\sin \vartheta \sin(\varphi - \varphi')}{\sin \vartheta' \cos \vartheta - \cos \vartheta' \sin \vartheta \cos(\varphi - \varphi')}, \quad (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 $\hat{u} = \{x, y, z\}$ and $\hat{u}' = \{x', y', z'\}$ as

$$ \alpha = \arctan \frac{yx' - xy'}{z' - z \cos \theta'}, \quad (C.3) $$

$$ \alpha' = -\arctan \frac{yx' - xy'}{z - z' \cos \theta'}, \quad (C.4) $$

with $\cos \theta = \hat{u} \cdot \hat{u}'$. This form is often useful in applications where points are available as vectors, since it requires no additional trigonometric operations. In fact, for $\alpha = \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),

$$ e^{2i\alpha} = \frac{p^2 - q^2 + 2i pq}{p^2 + q^2}. \quad (C.5) $$

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

## 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 $\langle C^{ff'}(\theta) \rangle$, i.e. the correlation of the complex conjugated random field $f^*$ and $f'$,

$$ \langle f(\hat{u}) f'(\hat{u}') \rangle = e^{i s \alpha} \langle C^{f'*f'}(\theta) \rangle e^{i s' \alpha'}, \qquad (\text{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<fnref n="11" />

$$ \langle (f^*)_{lm}^* f'_{lm'} \rangle = \delta_{ll'}^K \delta_{mm'}^K \langle C_l^{f'*f'} \rangle, \qquad (\text{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^*$,

$$ E_{lm} = - \frac{f_{lm} + (-1)^s (f^*)_{lm}}{2}, \qquad (\text{D.3}) $$

$$ B_{lm} = - \frac{f_{lm} - (-1)^s (f^*)_{lm}}{2i}, \qquad (\text{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$,

$$ Y(\hat{u}) = t(X(\hat{u})). \qquad (\text{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 $\delta$ for angular clustering that respect the physical bound $\delta \ge -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.

<figure><img src="image_14.png" /><figcaption>Fig. E.1. Shape of squared normal distributions (solid) and lognormal distributions with the same variance (dashed). For small variances, the distribution is approximately normal (black). For large variances, the distribution saturates at its lower bound (blue). In the intermediate regime, the distribution has a roughly lognormal shape (red).</figcaption></figure>

Let $X$ be a normal random variable with zero mean and variance $\sigma^2 \le 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 $\chi^2$ random variable with 1 degree of freedom,

$$ Y = \lambda [(X - a)^2 - 1], \qquad (\text{E.2}) $$

where $\lambda > 0$ is the scale parameter of the distribution, which also fixes the minimum value of $Y$.<fnref n="12" /> A straightforward calculation shows that $Y$ has zero mean and variance

$$ \langle Y^2 \rangle = 2\lambda^2 \sigma^2 (2 - \sigma^2). \qquad (\text{E.3}) $$

Inserting $\sigma^2 = 1 - a^2$, the variance can equivalently be expressed in terms of $a$,

$$ \langle Y^2 \rangle = 2\lambda^2 (1 - a^4), \qquad (\text{E.4}) $$

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

$$ a = \left(1 - \frac{\langle Y^2 \rangle}{2\lambda^2}\right)^{\frac{1}{4}}. \qquad (\text{E.5}) $$

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

$$ f_Y(y; a) = \frac{\exp\left(-\frac{a^2+y+1}{2(1-a^2)}\right) \cosh\left(\frac{a\sqrt{y+1}}{1-a^2}\right)}{\sqrt{2\pi}(1-a^2)(y+1)}. \qquad (\text{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

<footnote n="11">¹¹ The prefix ‘pseudo-’ is used here in the statistical sense, and not to be confused with meaning ‘partial sky’, for which it is unfortunately sometimes also used.</footnote>

<footnote n="12">¹² For that reason, the scale parameter $\lambda$ of a lognormal random variable is commonly called the “shift” parameter (Tessore et al. 2023).</footnote>
the sphere. By expectation (16), there is a correlation function $G$ such that $\langle X(\hat{u}) X'(\hat{u}') \rangle = 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 $\langle Y(\hat{u}) Y'(\hat{u}') \rangle = C(\theta)$. Using the transformation (E.2), we can compute $C(\theta)$ in terms of $G(\theta)$,

$$C(\theta) = 2\lambda\lambda' G(\theta) (G(\theta) + 2aa'). \qquad (E.7)$$

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

$$G(\theta) = \sqrt{\frac{C(\theta)}{2\lambda\lambda'} + (aa')^2 - aa'}. \qquad (E.8)$$

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

$$C_l = 2\lambda\lambda' \sum_{l_1 l_2} \frac{(2l_1 + 1)(2l_2 + 1)}{4\pi} \begin{pmatrix} l_1 & l_2 & l \\ 0 & 0 & 0 \end{pmatrix}^2 G_{l_1} G_{l_2} \\ + 4\lambda\lambda' aa' G_l. \qquad (E.9)$$

Since the triangle condition $l \le l_1 + l_2$ 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 $\lambda$. For angular clustering, we simulate the density contrast $\delta$, for which we set $\lambda = 1$. For cosmic shear, we simulate the convergence field $\kappa$, for which we set $\lambda$ using the fitting formula of Hilbert et al. (2011).