The Open Journal of Astrophysics

# Almanac: Weak Lensing power spectra and map inference on the masked sphere

Arthur Loureiro1,2,3 Lorne Whiteway2 Elena Sellentin4,5 Javier Silva
Lafaurie4,5 Andrew H. Jaffe3 Alan F. Heavens3 1 Institute for Astronomy,
University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ,
UK  2 Department of Physics and Astronomy, University College London, Gower
Street, London WC1E 6BT, UK  3 Astrophysics Group and Imperial Centre for
Inference and Cosmology (ICIC), Blackett Laboratory, Imperial College London,
Prince Consort Road, London SW7 2AZ, UK  4 Mathematical Institute, Leiden
University, Snellius Gebouw, Niels Bohrweg 1, NL-2333 CA Leiden, The
Netherlands  5 Leiden Observatory, Leiden University, Oort Gebouw, Niels
Bohrweg 2, NL-2333 CA Leiden, The Netherlands<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a field-based signal extraction of weak lensing from noisy
observations on the curved and masked sky. We test the analysis on a simulated
Euclid-like survey, using a Euclid-like mask and noise level. To make optimal
use of the information available in such a galaxy survey, we present a
Bayesian method for inferring the angular power spectra of the weak lensing
fields, together with an inference of the noise-cleaned tomographic weak
lensing shear and convergence (projected mass) maps. The latter can be used
for field-level inference with the aim of extracting cosmological parameter
information including non-gaussianity of cosmic fields. We jointly infer all-
sky $E$-mode and $B$-mode tomographic auto- and cross-power spectra from the
masked sky, and potentially parity-violating $EB$-mode power spectra, up to a
maximum multipole of $\ell_{\rm max}=2048$. We use Hamiltonian Monte Carlo
sampling, inferring simultaneously the power spectra and denoised maps with a
total of $\sim 16.8$ million free parameters. The main output and natural
outcome is the set of samples of the posterior, which does not suffer from
leakage of power from $E$ to $B$ unless reduced to point estimates. However,
such point estimates of the power spectra, the mean and most likely maps, and
their variances and covariances, can be computed if desired.

##  1 Introduction

Cosmic shear analysis has developed over the past few decades into a powerful
cosmological probe for studying the large-scale structure of the Universe
(Bartelmann & Schneider, 2001; Munshi et al., 2008; Kilbinger, 2015). Photons
from distant galaxies travel towards us along geodesic paths that are
distorted by the gravitational potential of dark and baryonic matter along the
way. The result is a coherent distortion in the observed shapes of background
galaxies: a weak gravitational lensing effect that shears and magnifies
objects (Bartelmann & Schneider, 2001). The shape distortions are sensitive to
the total matter content of the Universe and its clustering (Hall, 2021), and
so this effect makes cosmic shear an effective probe for studying the growth
of structure and the geometric evolution of the Universe (Heavens, 2009).
Measuring how the gravitational potential along the line-of-sight affects the
ellipticities of galaxies in photometric galaxy surveys also allows us to map
the invisible distribution of dark matter around galaxies (Kaiser & Squires,
1993; Marshall et al., 2002; Lanusse et al., 2016; Jeffrey et al., 2021;
Wallis et al., 2022) and to understand the structure of the cosmic web (Rossi,
2012; Joachimi et al., 2015; Codis et al., 2018).

Recent years have seen an impressive increase in both the quality and quantity
of weak lensing data and analysis techniques. Current photometric surveys such
as the Kilo Degree Survey (KiDS; Heymans et al. 2021; Asgari et al. 2021;
Loureiro et al. 2022), the Dark Energy Survey (DES; DES Collaboration et al.
2022; Amon et al. 2022; Secco et al. 2022), and the Subaru Hyper Supreme-Cam
survey (HSC; Mandelbaum et al. 2018; Hikage et al. 2019) have, together,
revolutionised cosmology by using cosmic shear to challenge its standard
model. Contemporary weak lensing experiments infer a growth of structure
parameter ($S_{8}$) that is in tension with the value inferred from the cosmic
microwave background (CMB) (Heymans et al., 2013; MacCrann et al., 2015; Lemos
et al., 2021; Amon & Efstathiou, 2022).

New data from forthcoming cosmic shear surveys such as Euclid (Laureijs et
al., 2011) and LSST (LSST Science Collaboration et al., 2009), combined with
improved data analysis techniques, could shed light on this tension and help
us better understand the geometry and growth of structure in the Universe.
These future weak lensing experiments will measure the shapes of billions of
galaxies across large portions of the sky with unprecedented accuracy and
precision. The question ‘how to extract optimally the cosmic shear information
from these data?’ is hotly debated in the field (Kilbinger et al., 2014;
Jarvis, 2015; Alonso et al., 2019; Nicola et al., 2021; Asgari et al., 2021;
Porqueres et al., 2021, 2022).

The vastly increased precision expected from Euclid and LSST raises the
question of whether the traditional approach of summary statistics such as
correlation functions and pseudo-power spectrum estimates are accurate enough
for robust inference. There is some evidence that shear correlation functions
fail (Hartlap et al., 2009; Sellentin & Heavens, 2018; Sellentin et al.,
2018), while pseudo-Cℓ estimates are close to gaussian (Upham et al., 2022).
Avoiding having to assume the gaussianity of such summaries is a major benefit
of field-level inference approaches, whose goal is to sample correctly from
the likelihood of the (close to) raw field-level data

There are two general approaches to field-level inference of cosmic shear: one
firmly based in a cosmological model with a physical gravity model, the other
(the approach taken in this paper) more agnostic, inferring the underlying
shear maps and their two-point statistics in a cosmology-independent way with
minimal assumptions.

In an example of the former, model-based approach, Porqueres et al. (2022)
showed that field-level inference can strongly lift the standard ‘banana’
cosmic shear parameter degeneracy (the degeneracy between the amplitude of
matter density fluctuations $\sigma_{8}$ and the total matter density
$\Omega_{\rm m}$). By extending the Borg framework (Jasche & Wandelt, 2013;
Jasche & Lavaux, 2019), which includes a full gravity model, sampling of the
initial gaussian density field is performed and forward modelled to the
tomographic shear fields, where the likelihood is applied. This is a very
powerful framework, but it ties the inference to a particular cosmological
model.

The approach of this work — as embodied in our algorithm, called Almanac — is
complementary, inferring the shear maps and the power spectra of the
(denoised, unmasked) fields in a cosmology-independent way. The samples of
maps and power spectra produced by Almanac are not tied to a particular model
(beyond a few basic symmetry assumptions), and can subsequently be used for
parameter inference (although the goal of this paper is simply to draw map and
power spectra samples). The main challenge is efficient scaling with the
highest multipole probed, so we concentrate on a relatively simple 2-bin
tomographic example. More details of Almanac are presented in a companion
paper (Sellentin, Loureiro et al., in prep).

The paper is structured as follows: weak lensing theory in section 2 is
followed by the hierarchical data model in section 3, and the Cholesky
reparameterisation in section 4. Simulations are described in section 5, and
results and convergence tests presented in section 6. We give our conclusions
in section 7.

##  2 Weak lensing theory

Galaxies are biased tracers of the dark matter density field. By contrast,
weak distortions in the observed shapes of background galaxies (caused by the
deflection of light by foreground structures) are an excellent tracer of the
projected matter field. This section outlines the formalism behind using this
weak lensing as a cosmological probe. Further details can be found in reviews
by Bartelmann & Schneider (2001) and Kilbinger (2015).

The usual approximation of the lens equation maps the unperturbed two-
dimensional source angular position $\mathbf{\theta}_{s}$ to the observed
position $\mathbf{\theta}_{\text{obs}}$ via a mapping ${\sf{A}}_{ij}$:
$\theta_{s,i}\approx{\sf{A}}_{ij}\theta_{\text{obs},j}$. This mapping can be
expressed in terms of the effective lensing potential ${\tilde{\Psi}}{}$ as

${\sf{A}}_{ij}=\delta_{ij}-\partial_{i}\partial_{j}{\tilde{\Psi}}{},$ (1)

where ${\tilde{\Psi}}{}$ is a weighted projection along the line-of-sight of
the three-dimensional Newtonian potential ${\Psi}{}$ (Kaiser, 1992):

${\tilde{\Psi}}{}(\chi_{\rm
s},{\hat{n}}{})=2\int_{0}^{\chi_{s}}d\chi\frac{f_{K}(\chi_{\rm
s}-\chi)}{f_{K}(\chi)f_{K}(\chi_{s})}{\Psi}{}(\chi,{\hat{n}}{})\,.$ (2)

Here $\chi_{\rm s}$ is the comoving distance to the lensing source, and
$f_{K}$ is related to the curvature $K$ of the universe (for a flat universe
$f_{K=0}(\chi)=\chi$, where $\chi$ is a comoving distance variable). We also
define the tomographic-bin averaged quantity

${\tilde{\Psi}}{}({\hat{n}})=\int d\chi\
\bar{n}(\chi){\tilde{\Psi}}{}(\chi,{\hat{n}})$ (3)

where $\bar{n}(\chi)$ is the normalised density of sources (averaged over the
sky); ${\tilde{\Psi}}{}$ is a scalar (i.e., spin-0) field on the sky.

The mapping in Eq. (1) can also be written in terms of the spin-0 lensing
convergence field $\kappa$ and the complex spin-2 shear field
$\bm{\gamma}=\gamma_{1}+\textrm{i}\mkern 1.0mu{}\gamma_{2}$:

$\displaystyle\mathsf{A}=\left(\begin{array}[]{cc}1-\kappa-\gamma_{1}&-\gamma_{2}\\\
-\gamma_{2}&1-\kappa+\gamma_{1}\end{array}\right).$ (6)

The potentially observable quantities $\kappa$ and $\gamma$ are related to the
lensing potential ${\tilde{\Psi}}{}$ by

$\kappa=\frac{1}{4}(\eth\bar{\eth}+\bar{\eth}\eth){\tilde{\Psi}}{}$ (7)

and

$\bm{\gamma}=\gamma_{1}+\textrm{i}\mkern
1.0mu{}\gamma_{2}=\frac{1}{2}\eth\eth{\tilde{\Psi}}{},$ (8)

where $\eth$ and $\bar{\eth}$ are the spin-raising and spin-lowering
differential operators (see, e.g., Castro et al. (2005) for details).

An arbitrary spin-$s$ function $f$ on the sphere may be represented in the
basis of spin-$s$ spherical harmonic functions ${}_{s}Y_{\ell m}$:

$f({\hat{n}})=\sum_{\ell m}f_{\ell m}\ {}_{s}Y_{\ell m}({\hat{n}})$ (9)

where

$f_{\ell m}=\int d\Omega\ f({\hat{n}})\ {}_{s}Y_{\ell m}^{*}({\hat{n}}).$ (10)

In this basis Eqs. (7) and (8) become

$\displaystyle\kappa_{\ell m}$
$\displaystyle=-\frac{1}{2}\ell(\ell+1)\,{\tilde{\Psi}}{}_{\ell
m}\quad\textrm{and}$ (11) $\displaystyle\gamma_{\ell m}$
$\displaystyle=\frac{1}{2}\sqrt{(\ell-1)\ell(\ell+1)(\ell+2)}\,{\tilde{\Psi}}{}_{\ell
m}\,,$ (12)

respectively, and hence

$\gamma_{\ell m}=-\sqrt{\frac{(\ell-1)(\ell+2)}{\ell(\ell+1)}}\kappa_{\ell
m}.$ (13)

The shear field may also be decomposed into $E$\- and $B$-modes, with
coefficients

$E_{\ell m}=-\frac{1}{2}\int{\rm
d}\Omega\left[\gamma({\hat{n}}{})\,_{+2}Y^{*}_{\ell
m}({\hat{n}}{})+\gamma^{*}({\hat{n}}{})\,_{-2}Y^{*}_{\ell
m}({\hat{n}}{})\right]$ (14)

and

$B_{\ell m}=\frac{\textrm{i}\mkern 1.0mu{}}{2}\int{\rm
d}\Omega\left[\gamma({\hat{n}}{})\,_{+2}Y^{*}_{\ell
m}({\hat{n}}{})-\gamma^{*}({\hat{n}}{})\,_{-2}Y^{*}_{\ell
m}({\hat{n}}{})\right]\,.$ (15)

The integrals are taken over the full sphere. The advantage of this
decomposition is that to a good approximation, cosmic shear produces only
$E$-modes (e.g. Schneider, 2005; Bartelmann & Maturi, 2017), and the presence
of $B$-modes can indicate systematic errors in the data (e.g. Heymans et al.,
2006; Asgari et al., 2019). These definitions use the $H=1$ sign convention
(in the notation of Lewis (2005)), as is appropriate for HEALPix; one
consequence is that $\gamma=-(E+\textrm{i}\mkern 1.0mu{}B)$.

Galaxy surveys directly observe neither shear nor convergence, observing
instead galaxy shape ellipticities $\bm{\epsilon}$ (and optionally sizes
and/or fluxes). These shapes are noisy estimates of the shear field in the
weak lensing regime ($\kappa\ll 1$). In this regime one has, for a large
number of galaxies,
$\langle\bm{\epsilon}\rangle\approx\bm{\gamma}/(1-\kappa)\approx\bm{\gamma}$.
Both shape noise (due to the random intrinsic shapes of galaxies) and
measurement noise (a smaller effect than shape noise) are taken into account
in the standard deviation $\sigma_{\epsilon}$ of the distribution of the
ellipticities; this distribution is usually assumed to be approximately
gaussian for a sufficiently large area density of sources (Joachimi & Bridle,
2010; Van Waerbeke et al., 2013; Alsing et al., 2016; Chang et al., 2018).
Intrinsic alignments can also introduce correlations of ellipticities (e.g.,
Heavens et al., 2000; Catelan et al., 2001; Kirk et al., 2015; Kiessling et
al., 2015; Joachimi et al., 2015; Blazek et al., 2015), and this effect will
be included in the inferred power spectra.

$\pi({\sf{C}})$${\sf{C}}_{\ell}$$\mathcal{G}(\bm{a}|{\sf{C}}_{\ell})$$\ell$$\bm{a}$$\mathcal{L}(\bm{d}|\bm{a},{\sf{N}})$${\sf{N}}$$\bm{d}$
Figure 1:  Directed acyclic graph representing the Bayesian hierarchical model
in Eq. (16) for maps and angular power spectra inference. Begin by drawing the
angular power spectra ${\sf{C}}_{\ell}$ from a prior distribution
$\Pi({\sf{C}}_{\ell})$. Then draw a set of spherical harmonics coefficients
for the shear fields $\bm{a}$, conditional upon the angular power spectra.
Finally draw a set of noisy tomographic cosmic shear fields given the signal,
$\bm{a}$, and the noise covariance, ${\sf{N}}$, based on the realisations from
the previous steps.

##  3 The posterior and HMC sampling

###  3.1 Joint posterior of fields and their covariance

The posterior from which we draw samples consists of the full sky maps in
harmonic space $\bm{a}$ and their power spectra ${\sf{C}}$, given the data
$\bm{d}$ and the noise covariance ${\sf{N}}$. Fig. 1 shows a directed acyclic
graph for the model presented in this section. Here $\bm{a}$ and $\bm{d}$ are
concatenations of the spherical harmonic coefficients and the field values
from different tomographic bins into a single map parameter vector and a
single data vector. The power spectrum ${\sf{C}}$ contains auto- and cross-
power spectra across tomographic bins, with $E$\- and $B$-mode correlations,
as well as the parity-violating $EB$-mode cross-correlations. For more
details, see the companion paper (Sellentin, Loureiro et al., in prep). The
posterior is

$\mathcal{P}({\sf{C}},\bm{a}|\bm{d},{\sf{N}})\propto\mathcal{L}(\bm{d}|\bm{a},{\sf{N}})\mathcal{G}(\bm{a}|{\sf{C}})\pi({\sf{C}})\,.$
(16)

The first factor is the probability at field level:

$\mathcal{L}(\bm{d}|\bm{a},{\sf{N}})\propto\exp\left(-\frac{1}{2}(\bm{d}-{\sf{Y}}\bm{a})^{\rm
T}{\sf{N}}^{-1}(\bm{d}-{\sf{Y}}\bm{a})\right)\,.$ (17)

The second factor assumes a gaussian prior for the spherical harmonic
coefficients:

$\mathcal{G}(\bm{a}|{\sf{C}})=\frac{1}{\sqrt{|2\pi{\sf{C}}|}}\exp\left(-\frac{1}{2}\bm{a}^{{\textrm{T}}{}}{\sf{C}}^{-1}\bm{a}\right).$
(18)

This factor merits some discussion. The aim of the method is to infer the
power spectra and the maps in a cosmology-independent way, making no
assumptions about the gravity model. Thus, in contrast to BORG-WL (Porqueres
et al., 2021, 2022), we do not sample the initial gaussian field and then
forward model. Instead, we want the least informative prior for the
coefficients, for a given power spectrum, with no reliance on a model-
dependent knowledge (in principle) of the statistics of the shear fields. This
least informative, maximum-entropy prior is a gaussian (Jaynes, 2003), which
ensures that we are not using any model-dependent information. This does not
imply that the fields themselves are realisations of a zero-mean gaussian with
isotropic covariance, nor that a map created from the spherical harmonic
coefficients of a given posterior sample is a realisation of such a process,
since these are conditioned on the data. This is similar to the known
behaviour of Wiener filters: where the data are highly constraining, the
posterior will have small variance around the data, but where they are noisy,
the variance will be dominated by the prior information about the signal.

The last factor in the posterior is a prior on the angular power spectra. In
all cases we set the prior to zero outside the set of positive definite
matrices, enforced by the coordinate choices described in Section 4. We then
assign a constant prior, $\pi({\sf{C}})=\text{const}$, but note that this
choice is essentially cosmetic. Any use of these results for further analysis
will require a prior on cosmological parameters $\theta$, corresponding to a
delta-function distribution
$\pi(\theta,{\sf{C}})=\pi(\theta)\delta({\sf{C}}-{\sf{C}}[\theta])$. In this
case the constant prior is equivalent to sampling from the likelihood function
itself, which can then be modulated by such a cosmological prior.

More generally, we could choose the prior to be a power of the determinant,
$\pi({\sf{C}})=|{\sf{C}}|^{q}$, enforcing different information about the
scaling of the spectra. A Jeffreys prior corresponds to $q=-(N_{p}+1)/2$
(where ${\sf{C}}_{\ell}$ has $N_{p}\times N_{p}$ entries), as discussed in
Sellentin & Heavens (2016); a frequency-matching choice of $q$ can also be
made (Percival et al., 2022) if desired. All of these would lead to different
detailed inferences, such as regarding the significance of non-zero power of a
particular ${\sf{C}}_{\ell}$.

###  3.2 Hamiltonian Monte Carlo

Our aim is to infer the denoised maps and power spectra from simulated input
data, using the Bayesian hierarchical model depicted in Fig. 1. Only a few
samplers can handle such a high-dimensional parameter space. Gibbs sampling
(Geman & Geman, 1984; Casella & George, 1992) as used by Wandelt et al.
(2004); Eriksen et al. (2004); Larson et al. (2007); Alsing et al. (2016);
Alsing et al. (2017); Colombo et al. (2022) can generate samples using the
conditional distributions $\mathcal{P}(\bm{a}|{\sf{C}},\bm{d})$ and
$\mathcal{P}({\sf{C}}|\bm{a},\bm{d})$; however, this method is probably too
inefficient to sample the full posterior distribution in Eq. (16) for upcoming
weak lensing analyses.

Instead, we take an approach similar to that used by most field-level
inference works (Jasche & Kitaura, 2010; Jasche & Wandelt, 2013; Leclercq et
al., 2015; Jasche & Lavaux, 2019; Porqueres et al., 2021, 2022): we implement
a tuned version of the HMC algorithm. This section outlines the key details;
see the Almanac core paper (Sellentin, Loureiro et al., in prep) for more
information.

In HMC (Hanson, 2001; Hajian, 2007; Neal, 2011) we view the negative logarithm
of the posterior density as a potential energy

${\psi}{}(\bm{a},{\sf{C}}|\bm{d})=-\ln\mathcal{P}(\bm{a},{\sf{C}}|\bm{d}).$
(19)

We then augment the parameter space with momentum variables drawn from a
gaussian distribution, whose negative logarithm can be viewed as kinetic
energy. A Hamiltonian trajectory of the resulting dynamical system is then
computed using a leapfrog integrator. At the end of a trajectory, a
Metropolis-Hastings step (Metropolis et al., 1953; Hastings, 1970) decides if
the trajectory is accepted (this corrects for any non-conservation of energy
by the leapfrog integrator). If it is accepted, the new sample is the
trajectory’s end point; if it is rejected, the new sample is the trajectory’s
start point. In either case, the momenta are resampled and another trajectory
is computed (starting at the new sample). By iteration, the sampler builds up
a Monte Carlo Markov Chain.

Step sizes and other parameters of the leapfrog integrator must be carefully
tuned; see the companion paper (Sellentin, Loureiro et al., in prep) for
details of this tuning. The sampling is challenging, since the shape of the
posterior in the $\bm{a}\times{\sf{C}}$ plane exhibits a ‘stingray’ geometry
(Neal, 2000; Betancourt & Girolami, 2013) that results in large correlation
lengths and slow convergence; a reparametrisation of the parameter space, as
discussed in the following section, can improve convergence.

##  4 Cholesky decomposition coordinates

(a)

(b)

Figure 2: A toy example illustrating the geometry of the potential density in
(a) the matrix-log parameterisation $\\{\bm{a},\sf{G}\\}$ and in (b) the
Cholesky decomposition parameterisation $\\{\bm{x},{\sf{K}}\\}$. This toy
example shows a slice in the map $\times$ angular power spectra dimensions for
the negative log-posterior (potential) illustrated by the DAG in Fig. 1 using
different coordinate systems. The potential density varies from low (pink) to
high (blue) and values have been omitted since the normalisation is arbitrary.
The matrix-log parameterisation has a stingray shape, also known as Neal’s
Funnel (Neal, 2000); this posterior shape is known to be difficult to sample
from.

The particular challenge of tomographic weak lensing is that the covariance
matrix, although block-diagonal in $\ell$ and $m$, is dense for the $E$ and
$B$ correlations across $N_{\rm b}$ tomographic bins. Thus there are dense
sub-block matrices of size $2N_{\rm b}\times 2N_{\rm b}$, with strong
correlations between neighbouring bins. We cannot straightforwardly sample the
elements of ${\sf{C}}$ since the HMC random walk can easily wander out of the
subset of positive-definite matrices. To handle this we instead explored
sampling either a) the matrix logarithm ${\sf{G}}$ (where
${\sf{C}}=\exp({\sf{G}})$) or b) the ‘diagonal-log’ ${\sf{K}}$ (see Eq. (22))
of the Cholesky factor ${\sf{L}}$ (where ${\sf{C}}={\sf{L}}{\sf{L}}^{\rm T}$
with ${\sf{L}}$ lower-triangular) of the covariance matrix. Figure 2 shows,
for a toy example, the geometry of the potential (the negative of the
logarithm of the posterior density) using both parameterisations. For the
strong correlations inherent in cosmic shear, the Cholesky decomposition was
found to lead to chains with shorter correlation lengths (indeed, the matrix
log formalism failed badly for this application). See below for further
details of the Cholesky reparametrisation here.

The exact shape of the posterior density of course depends on the coordinates
used. When parameterising using natural coordinates $\\{\bm{a},{\sf{C}}\\}$
the posterior density has an irksome shape in the vicinity of
$\bm{a}=0,{\sf{C}}=0$. At such points there is poor agreement between $\bm{a}$
and the data (this reduces the posterior density) but excellent agreement
between $\bm{a}$ and ${\sf{C}}$ (this increases the posterior density). The
former factor is bounded, while the latter is unbounded as we approach
${\sf{C}}=0$. Thus there is an asymptote: the posterior density becomes
infinite at $\bm{a}=0,{\sf{C}}=0$ (note though that this behaviour affects
only a small volume in parameter space; the integrated probability at the
asymptote is tiny). The asymptotic behaviour contrasts, e.g., with a gaussian
distribution, where the density reaches a finite maximum (a ‘peak’) within the
interior of the parameter space. If we reparameterise to coordinates
$\\{\bm{a},{\sf{G}}=\ln({\sf{C}})\\}$ then the same asymptotic behaviour
occurs, only now at ${\sf{G}}$ near minus infinity. However when ${\sf{G}}$ is
large negative the posterior density falls off very quickly when the $\bm{a}$
move away from zero. Thus in these coordinates the posterior density shape
will resemble a funnel, or a stingray (Fig. 2(a)). Such an extreme shape is
difficult to sample from using an HMC sampler. A remedy can be to find a
different parameterisation, in order to avoid the occurrence of the stingray
(or, as we do, to move it to a less bothersome location in parameter space).
For alternative approaches, see Jewell et al. (2009), Racine et al. (2016),
Millea et al. (2021).

In an analysis of a single $a_{\ell m}$ and a single $C_{\ell}$, the ratio
$x=a_{\ell m}/\sqrt{C_{\ell}}$ will always be of unit variance. Therefore,
although it may not solve all problems, switching to $x$ can avoid the
occurrence of a stingray for small $C_{\ell}$. As we analyze many multipoles
and power spectra jointly, the multivariate analogue is to flatten all
$\bm{a}_{\ell m}$ modes of the same $\ell$ jointly. This can be achieved by
Cholesky decomposing the covariance matrices ${\sf{C}}_{\ell}$.

Each matrix ${\sf{C}}_{\ell}$
($\ell_{\textrm{min}}\leq\ell\leq\ell_{\textrm{max}}$) is relatively small
(for example, dimension $3$ for a CMB $TEB$-analysis, and dimension $2N_{\rm
b}$ for a cosmic shear analysis with $N_{\rm b}$ redshift bins and both $E$\-
and $B$-modes) and hence can be decomposed at reasonable numerical cost.

In the rest of this section we fix $\ell$ (between $\ell_{\rm min}$ and
$\ell_{\rm max}$) and for notational simplicity we drop the dependence on
$\ell$ (so for example $\bm{a}$ now refers to $a_{\ell m}$ for varying $m$ and
this fixed $\ell$, while ${\sf{C}}$ refers to what was previously called
${\sf{C}}_{\ell}$).

Each matrix ${\sf{C}}$ is relatively small (for example, dimension $3$ for a
CMB $TEB$-analysis, and dimension $2N_{\rm b}$ for a cosmic shear analysis
with $N_{\rm b}$ redshift bins and both $E$\- and $B$-modes) and hence can be
decomposed at reasonable numerical cost.

For positive definite ${\sf{C}}$ the Cholesky factor ${\sf{L}}$ is the unique
lower-triangular matrix with positive elements on the diagonal satisfying
${\sf{C}}={\sf{L}}{\sf{L}}^{{\textrm{T}}{}}$. Given $\bm{a}$ define $\bm{x}$
via

$\bm{x}={\sf{L}}^{-1}\bm{a};$ (20)

then if $\bm{a}$ has mean zero and covariance ${\sf{C}}$ then $\bm{x}$ will
have mean zero and unit covariance

$\langle\bm{x}\bm{x}^{{\textrm{T}}{}}\rangle=\mathbb{I}.$ (21)

We define the diagonal-log ${\sf{K}}$ of ${\sf{L}}$ to be

${\sf{K}}_{\alpha\beta}=\left\\{\begin{array}[]{ll}\ln({\sf{L}}_{\alpha\beta})&\mbox{if
}\alpha=\beta\,,\\\
{\sf{L}}_{\alpha\beta}&\mbox{otherwise.}\end{array}\right.$ (22)

Then $\\{\bm{x},{\sf{K}}\\}$ provides an alternative parameterisation of
$\\{\bm{a},{\sf{C}}\\}$; note that an arbitrary ${\sf{K}}$ will generate a
positive-definite ${\sf{C}}$.

###  4.1 Jacobian for Cholesky coordinates

Reparameterising from $\\{\bm{a},{\sf{C}}\\}$ to $\\{\bm{x},{\sf{K}}\\}$
introduces a Jacobian factor in the posterior density. Now

$p(\bm{x},{\sf{K}})\ \mathrm{d}\bm{x}\ \mathrm{d}{\sf{K}}=p(\bm{a},{\sf{C}})\
\mathrm{d}\bm{a}\ \mathrm{d}{\sf{C}}$ (23)

and so the Jacobian factor is

$|{\sf{J}}|=\begin{vmatrix}\partial\bm{a}/\partial\bm{x}&\partial\bm{a}/\partial{\sf{K}}\\\
\partial{\sf{C}}/\partial\bm{x}&\partial{\sf{C}}/\partial{\sf{K}}\\\
\end{vmatrix}=|\partial\bm{a}/\partial\bm{x}|\
|\partial{\sf{C}}/\partial{\sf{K}}|;$ (24)

the second equality follows from ${\partial{\sf{C}}}/{\partial\bm{x}}=0$
(since it is evaluated at constant ${\sf{K}}$). We have
$\bm{a}={\sf{L}}\bm{x}$ and so

$|\partial\bm{a}/\partial\bm{x}|=|{\sf{L}}|=\prod_{\alpha=1}^{n}{\sf{L}}_{\alpha\alpha}.$
(25)

The determinant of the transformation from ${\sf{C}}$ to ${\sf{L}}$ is (Gupta,
1999):

$|\partial{\sf{C}}/\partial{\sf{L}}|=2^{n}\prod_{\alpha=1}^{n}{\sf{L}}_{\alpha\alpha}^{n+1-\alpha},$
(26)

while

$|\partial{\sf{L}}/\partial{\sf{K}}|=\prod_{\alpha=1}^{n}\frac{\partial{\sf{L}}_{\alpha\alpha}}{\partial{\sf{K}}_{\alpha\alpha}}=\prod_{\alpha=1}^{n}{\sf{L}}_{\alpha\alpha}.$
(27)

Combining these results yields

$|{\sf{J}}|=2^{n}\prod_{\alpha=1}^{n}{\sf{L}}_{\alpha\alpha}^{n+3-\alpha}=2^{n}\exp\left[\sum_{\alpha=1}^{n}(n+3-\alpha){\sf{K}}_{\alpha\alpha}\right]\
.$ (28)

###  4.2 Negative Log Posterior in Cholesky coordinates

In Cholesky coordinates $\\{\bm{x},{\sf{K}}\\}$ the negative logarithm
${\psi}{}$ of the posterior (Eq. (16)) becomes:

$\begin{split}{\psi}{}(\bm{x},{\sf{K}})&=\frac{1}{2}(\bm{d}-{\sf{Y}}{\sf{L}}\bm{x})^{\rm
T}{\sf{N}}^{-1}(\bm{d}-{\sf{Y}}{\sf{L}}\bm{x})+\frac{\bm{x}^{\rm
T}\bm{x}}{2}\\\ &+\sum_{\alpha=1}^{n}(\alpha-2-n-2q){\sf{K}}_{\alpha\alpha}\
.\end{split}$ (29)

Here we have included the Jacobian term (Eq. (28)), we have used for example
$\bm{a}={\sf{L}}\bm{x}$ and $|{\sf{C}}|=|{\sf{L}}|^{2}$, and we have silently
dropped irrelevant additive constants. The final summand includes an implied
sum over $\ell$.

This new coordinate system creates a posterior shape that can be explored more
efficiently by a HMC sampler than the shape of the posterior in the
$\\{\bm{a},\ln({\sf{C}})\\}$ coordinates. The data products outputted by the
sampler are transformed back to the usual $\\{\bm{a},{\sf{C}}\\}$-parameters
that a cosmologist would expect.

###  4.3 Gradients in Cholesky coordinates

The HMC sampler requires for its leapfrog routine the derivatives of the
potential (i.e., the negative logarithm of the posterior) ${\psi}{}$ with
respect to both $\bm{x}$ and ${\sf{K}}$ parameters. Straightforward
calculations show that the derivative with respect to $\bm{x}$ is

$\frac{\partial{\psi}{}}{\partial\bm{x}_{\alpha}}=-[{\sf{L}}^{{\textrm{T}}{}}{\sf{Y}}^{{\textrm{T}}{}}{\sf{N}}^{-1}(\bm{d}-{\sf{Y}}{\sf{L}}\bm{x})]_{\alpha}+\bm{x}_{\alpha}\,$
(30)

and that the derivative with respect to ${\sf{K}}$ is

$\begin{split}\frac{\partial{\psi}{}}{\partial{\sf{K}}_{\alpha\beta}}=&\Big{(}-(\delta_{\alpha\beta}^{K}({\sf{L}}_{\alpha\beta}-1)+1)\\\
&\,\times\bm{x}_{\beta}[{\sf{Y}}^{{\textrm{T}}{}}{\sf{N}}^{-1}(\bm{d}-{\sf{Y}}{\sf{L}}\bm{x})]_{\alpha}\Big{)}\\\
&\,+\delta_{\alpha\beta}^{K}(\alpha-2-n-2q)\ ;\end{split}$ (31)

here $\delta^{K}$ is the Kronecker delta.

###  4.4 Hessian in Cholesky coordinates

The Hessian of the posterior can be a good approximation for the mass matrix
of an HMC sampler. If one does not want to (or cannot) use the full Hessian,
one could still set the leapfrog step size to be the inverse square root of
the diagonal of the Hessian; Almanac sets the initial step sizes in this way
(although the step sizes are later tuned to better values). The diagonal
elements of the Hessian are

$\frac{\partial^{2}{\psi}{}}{\partial\bm{x}_{\alpha}\partial\bm{x}_{\alpha}}=[{\sf{L}}^{{\textrm{T}}{}}{\sf{Y}}^{{\textrm{T}}{}}{\sf{N}}^{-1}{\sf{Y}}{\sf{L}}]_{\alpha\alpha}+1$
(32)

and

$\begin{split}\frac{\partial^{2}{\psi}{}}{\partial{\sf{K}}_{\alpha\beta}\partial{\sf{K}}_{\alpha\beta}}=&(\delta_{\alpha\beta}^{K}({\sf{L}}_{\alpha\beta}^{2}-1)+1)\bm{x}_{\beta}^{2}[{\sf{Y}}^{{\textrm{T}}{}}{\sf{N}}^{-1}{\sf{Y}}]_{\alpha\alpha}\\\
&\,-\delta_{\alpha\beta}^{K}{\sf{L}}_{\alpha\beta}\bm{x}_{\beta}[{\sf{Y}}^{{\textrm{T}}{}}{\sf{N}}^{-1}(\bm{d}-{\sf{Y}}{\sf{L}}\bm{x})]_{\alpha}\
.\end{split}$ (33)

We may estimate the diagonal terms of
$[{\sf{Y}}^{{\textrm{T}}{}}{\sf{N}}^{-1}{\sf{Y}}]_{\alpha\alpha}$ using Monte
Carlo integration, as described in the Appendix of Sellentin, Loureiro et al.,
in prep. This technique allows us to calculate Eq. 32 and the first term in
Eq. 33; however, the second term in Eq. 33 is challenging to calculate
efficiently and we therefore estimate it by numerical differentiation.

We have suppressed in this outline some details that are required for
practical implementation. The values in ${\sf{L}}$, when applied to $\bm{x}$,
must be used repeatedly ($2\ell+1$ times); this introduces $2\ell+1$ factors
into the log determinant of the Jacobian (which then cancels out of the
negative log posterior). It also introduces a sum over $2\ell+1$ terms in the
${\sf{K}}$ gradient, and a double sum in the ${\sf{K}}$ hessian. We must also
account for the variance being split between the real and imaginary parts of
$a_{\ell m}$ except when $m=0$; this introduces factors of $\sqrt{2}/2$ when
$m\neq 0$.

###  4.5 Reduced shear

The quantity that controls the shape distortion is not strictly the shear, but
rather the reduced shear (e.g. Bartelmann & Maturi, 2017),
$\gamma/(1-\kappa)$. For the tests in this paper we have not included reduced
shear. In the future we intend to include reduced shear in the posterior
calculation. This will complicate the calculation of the derivatives and we
propose not including reduced shear in the derivatives calculation. The small
(percent-level) errors in the resulting trajectories will then be dealt with
completely by the Metropolis-Hastings acceptance/rejection step, which also
deals with inaccuracies in the leapfrog integrator.

##  5 Simulation

Figure 3: Setup for the Euclid-like cosmic shear survey simulation, showing
the input angular power spectra to the GRF simulation for the $E$-modes (green
line) and $B$-modes (orange line), the mask angular power spectra (dashed
black line), the noise power spectra (yellow dot-dashed line), and the data
realisation angular power spectra for $E$\- and $B$-modes (pink and blue dots,
respectively). The bottom left figure shows the angular mask used in the
simulation.

In order to demonstrate Almanac’s application to the next generation of cosmic
shear experiments, we used gaussian random fields (GRF, Bardeen et al. 1986;
Hoffman & Ribak 1992) to simulate a Euclid-like cosmic shear survey. Euclid is
an excellent showcase for Almanac’s performance and application due to sky
coverage, complex mask geometry, and expected data quality. Since Almanac
infers all scales available from a given multipole range to properly forward-
model the underlying field, galaxy surveys with large sky coverage such as
Euclid and LSST are ideal data sets for the methodology described in this
work.

The Euclid photometric sample will likely be divided into ten redshift
tomographic bins with an equal effective galaxy number density of $\bar{n}=3$
galaxies/arcmin2 for each bin (Laureijs et al., 2011; Euclid Collaboration et
al., 2020; Tutusaus et al., 2020; Taylor et al., 2021; Euclid Collaboration et
al., 2022). We arbitrarily select two of these ten tomographic bins ($0.65\leq
z<0.79$ and $1.54\leq z<1.83$) to define the weak lensing tracers for the
underlying theoretical angular power spectra. The fiducial ${\sf{C}}_{\ell}$
for the $E$-modes used in the simulation were calculated using CCL (Chisari et
al., 2019) with a Planck 2018 fiducial cosmology (Planck Collaboration et al.,
2020)111$\Omega_{\rm cdm}=0.27$; $\Omega_{\rm b}=0.045$; $\Omega_{\rm k}=0$;
$n_{\rm s}=0.96$; $\sigma_{8}=0.83$; $h=0.67$., in a range between $\ell_{\rm
min}=4$ and $\ell_{\max}=2048$. Meanwhile, the $B$-mode power spectra were set
to $10^{-15}$ and the ${\sf{C}}_{\ell}^{E_{i}B_{j}}$ were set to $10^{-20}$,
both constant across the multipole range used for the $E$-mode power spectra.

The survey angular geometry mask was created by masking out the galactic plane
and the ecliptic based on the description in Euclid Collaboration et al.
(2022). We then randomly masked pixels using a Poisson distribution to mimic
the effects of stars in the field. The resulting mask is shown in the bottom
left of Fig. 3.

The last element in the simulation is shape noise, taken to be
$\sigma_{\epsilon}=0.27$. Note (see Fig. 3) that the $B$-modes are almost six
orders of magnitude below the noise power spectra.

The noise covariance in real (pixel) space ${\sf{N}}$ is taken to be a
diagonal matrix. It has $2N_{b}N_{\rm pix}$ diagonal elements (the factor of
$2$ corresponding to the two components of shear); for a given shear
component, tomographic bin $i$, and pixel $p$ the corresponding entry in
${\sf{N}}$ is ${\sigma_{\epsilon}^{2}}/{(2\bar{n}_{i})}$ if $p$ is unmasked
and infinity otherwise. Here $\bar{n}_{i}$ is the lensing source number per
pixel in bin $i$.

The simulations are performed using a HEALPix (Gorski et al., 1999) resolution
of $N_{\rm side}=1024$ ($\Delta\Omega_{\rm pix}=11.8$ arcmin2) by sampling
correlated $\bm{a}$ from a multivariate gaussian distribution,
$\bm{a}\sim\mathcal{N}(0;{\sf{C}}_{\ell})$, and then transforming the signal
maps into real space. Next, using the aforementioned galaxy density and shape
noise, we add gaussian noise to the signal maps with a variance of
$\sigma^{2}=\sigma_{\epsilon}^{2}/(2\bar{n}_{i})$ per component. In a final
step, we apply the mask shown in Fig. 3 to all the galaxy ellipticity data
maps.

Although a more realistic N-body simulation with ray-tracing would include
higher-order effects into the tests, a GRF simulation allows for more control
over the quantities we are probing and testing. It lets us quickly simulate
correlated tomographic cosmic shear data, allowing full control and
understanding of the effect of signal, noise, and mask; a GRF is the most
reliable way to have full control of the simulation’s underlying angular power
spectrum and noise. This approach also allows us to keep fixed the
simulation’s signal realisation while turning off shape noise and mask effects
to produce the ground truth. More complex simulations should not pose any
difficulties for the Almanac analysis presented in Sec. 6 as the method is
sensitive to the statistical properties of the input data and makes no
assumptions on the underlying physics of the analysed data maps beyond the
spin, noise, and statistical isotropy.

##  6 Results

We apply the Almanac BHM described in Sec. 3 to the Euclid-like simulation
described in Sec. 5 probing a multipole range $4\leq\ell\leq 2048$ with
$N_{\rm side}=1024$. Here we describe the main outputs from this analysis: the
marginalised angular power spectrum and the inferred shear and convergence
full sky maps. We further demonstrate other benefits of having the full joint
posteriors of maps and angular power spectra, and we discuss implications for
future cosmic shear surveys (such as our approach to $EB$ leakage, and the
recovery of large scale modes).

We keep the $B$-mode auto-power spectra in the modelling, with a small signal
amplitude in the simulated data; although in the standard model of cosmology
these modes are expected to be negligible for cosmic shear, we want the
ability to detect systematic contamination that could appear through the
$B$-modes (Heymans et al., 2006; Asgari et al., 2019).

The full posterior contains all the spin-2 spherical harmonics for both
redshift tomographic bins and their angular power spectra; this results in a
total of 16,813,990 dimensions (of which only 20,450 are for the angular power
spectra, ${\sf{C}}_{\ell}$).

We ran two simultaneous chains using different starting points and different
random seeds. Each chain contains $N_{\rm samples}\approx 2.25\times 10^{5}$
samples (after burn-in and tuning) and each had an acceptance ratio of $\sim
86\%$. For both chains we used $N_{\rm tune}=2\times 10^{4}$ and $N_{\rm lf-
tune}=5\times 10^{3}$ samples to perform the step size tuning described in
Sellentin, Loureiro et al., in prep. Unless stated otherwise, all the results
presented in this section come from combining the chains (once the chains were
individually converged — see Sec. 6.3 for a discussion of the convergence
diagnostics for chains used in this analysis).

###  6.1 Inferred Angular Power Spectra

Figure 4: Inferred marginalised angular power spectra for the $E$\- and
$B$-modes. The blue (resp. light blue) hashed area is the 68.3% (resp. 95%)
C.I. The starting point for one of the chains is shown in light grey. The red
line is the ground truth given to the input simulation shown in Fig. 3. The
credible intervals for the $B$-mode auto-power spectra extend close to the
lower bound of zero. Figure 5: Similar to Fig. 4 but showing the parity
violating modes instead. Figure 6: Marginalised one-dimensional posterior
distributions for $E$-modes of different multipoles
($\ell=8,\,10,\,15,\,150,\,1600$). As expected, marginalised posteriors for
large scale modes are significantly non-gaussian while those for small scale
modes are much closer to gaussian (due to the higher number of spherical
harmonic coefficients available).

For each element of the power spectra we calculate the one-dimensional
marginalised 68.3% and 95% credible intervals (C.I.) from the sampled
posterior. The results are shown in Fig. 4 for the $E$\- and $B$-modes in both
tomographic bins and in Fig. 5 for the parity violating modes. The natural
outcome of the analysis is the joint posterior of the power spectra and the
maps. This high-dimensional distribution is not simple in shape, so we do not
show point estimates (e.g., mean, mode), as they do not necessarily act as
useful summaries of the posterior. The apparent noisy spikiness in the
marginalised angular power spectra shown in Figs. 4 and 5 is expected since
the marginalised multipoles are independent and since we are not binning our
results in bandpowers.

Fig. 4 demonstrates Almanac’s ability to infer the underlying angular power
spectrum for cosmic shear data on the sphere up to scales relevant for the
current and next generation of surveys. Although not a requirement for a
Bayesian analysis, the marginalised $E$-mode ${\sf{C}}_{\ell}$s agree well
with the input power spectra given to the Euclid-like simulations within the
68.3% C.I. for all multipoles in the analysed range, including the largest
scales available in the simulation. Our prior for ${\sf{C}}$ excludes non-
positive-definite matrices, and hence all diagonal elements of ${\sf{C}}$ must
be positive; in particular, the $B$-mode auto-power spectra cannot be zero or
negative. The low signal-to-noise ratio of the B-mode observations implies
that the marginalised $C_{\ell}^{BB}$ posterior will often peak at positive
values even when no such signal is present. This is evident from the form of
the map-marginalised posterior (Bond et al., 1998). Moreover, sampling near
the $C_{\ell}^{BB}=0$ boundary is challenging due to the funnel shape of the
posterior in the full $(a_{\ell m},C_{\ell})$ parameter space. This region
comprises a very small fraction of the full parameter volume. In most cases,
the credible intervals extend to or very near the lower bound. Both the
$B$-mode upper bound and the parity-violating modes are efficient probes for
systematic contamination that can be calculated in a fully Bayesian way with
the measurements produced by Almanac.

We show several examples of marginalised one-dimensional posterior
distributions of $E$-mode multipoles in Fig. 6. This plot illustrates that the
method can properly infer scales that would usually be considered too large
(given the fraction of the sky made available by the survey’s mask), since the
data in the unmasked regions, in conjunction with the prior, constrain modes
on all scales.

A particular issue when using point estimators of two-point statistics for
weak lensing data is the occurrence of ambiguous modes due to the survey’s
mask and anisotropic noise. Commonly referred to as ‘$EB$ leakage’ (Bunn et
al., 2003; Bunn, 2008; Leistedt et al., 2017; Liu et al., 2019; Nicola et al.,
2021), this effect causes high (low) $B$-modes to be confused with low (high)
$E$-modes, yielding a negative correlation between these modes on the
marginalised posterior distribution for a given multipole. Since we are not
working with estimators of the $E$\- and $B$-mode power spectra, instead
inferring the full posterior of the underlying whole sky fields and their
angular power spectra, such leakage of power from $E$ to $B$ does not arise.
What may generically happen is that the posterior samples of $E$\- and
$B$-modes are correlated. To check for this effect, we calculate the
$r$-correlation factor for each pair of auto power spectra samples
$\\{C^{E_{i},E_{i}}_{\ell=L},C^{B_{i},B_{i}}_{\ell=L}\\}$ for a given
multipole $L$. The $r$-correlation factors for both tomographic bins are
plotted in Fig. 7, which shows a rather small level of correlation. Although
we find no significant correlations between $E$\- and $B$-modes in our
analysis, we outline here another advantage of a Bayesian approach. Since we
have the full posterior of angular power spectra, modes that demonstrate
significant correlation can be removed by marginalising. We highlight,
however, that since the full posterior distribution is known, this
marginalisation would only be necessary if one were to take point estimates of
the marginalised individual modes of the $E$-mode angular power spectra.

Figure 7: Correlations between the samples of same multipoles between auto-
power spectra of $E$\- and $B$-modes. For each pair of $C_{\ell}^{EE}$ and
$C_{\ell}^{BB}$ highly anti-correlated modes indicate that lower (higher)
$B$-modes are degenerate with higher (lower) $E$-modes. We find no evidence
for this effect in our analysis.

###  6.2 Inferred Shear and Convergence Maps

Figure 8: Maps for the first redshift tomographic bin ($0.65\leq z<0.79$). Row
1: simulated galaxy ellipticity data; row 2: its underlying cosmic shear
signal (fiducial); row 3: recovered mean shear components; row 4: errors for
the components of the recovered shear; row 5: typical convergence (left) and
its signed S/N ratio (right). The fiducial maps contain the same signal
realisation as the simulated maps but without the shape noise and without the
mask. A ‘typical’ map is a randomly selected (post burn-in) sample. The mean
maps recovered most of the underlying structure shown in the fiducial maps in
areas where data is present and are able to recover some of the large-scale
information under the mask. The uncertainty on the recovered shear maps
reflects the data missing due to the survey geometry. Figure 9: Same as Fig. 8
but for the second tomographic bin ($1.54\leq z<1.83$). Figure 10:
Reconstructed projected gravitational potential (or lensing potential),
${\tilde{\Psi}}{}(\hat{n})$, from typical Almanac maps for both tomographic
bins. Top: the ground-truth lensing potential obtained from the fiducial
maskless noiseless simulations (second row on Fig. 8 and 9). Bottom: projected
lensing potential reconstructed from the simulated data (first row on Fig. 8
and 9) using a typical Almanac map. Here, we overlay the mask used in the
simulations for visualisation purposes.  Figure 11: Consistency check between
the inferred shear maps and the simulated signal map showing the $E$-mode
angular power spectra for a typical map. In red we show the angular power
spectra for the signal realisation in the simulations, while in grey we show
the pseudo-${\sf{C}}_{\ell}$ for the typical map. Figure 12: The residuals
from the samples of the shear component maps obtained by Almanac compared to
the underlying noiseless maskless signal map. The residuals are expected to be
distributed as a standard normal (and we therefore show such a standard normal
for comparison).

We next consider the spherical harmonic coefficients; these make up the larger
subset of the joint set of parameters. In this section we analyse the
posterior samples when marginalised over the angular power spectra. For each
$\bm{a}$ sample in the chains, we make a spherical harmonic transformation,
following the formalism stated in Sec. 2, to obtain real-space maps for shear
and convergence. The $\bm{a}$ samples are quite large, so to reduce file size
(and hence increase speed) only one $\bm{a}$ sample in 800 is saved; this
thinning yields 577 samples across both chains.

The main results for the inferred shear and convergence maps and their
uncertainties are shown in Fig. 8 (first bin, $0.65\leq z<0.79$) and Fig. 9
(second bin, $1.54\leq z<1.83$). In these figures, the data maps (first row)
show the simulated galaxy ellipticities (here data is missing due to the
survey geometry mask, and the maps contain shape noise according to the
description in Sec. 5). The fiducial maps (second row) show the same simulated
maps but without the shape noise and survey mask (i.e., they contain only the
underlying cosmic shear signal). The mean maps (third row) and standard
deviation maps (fourth row) are calculated from the statistics of the $\bm{a}$
samples transformed into real space. Finally, the ‘typical’ convergence map
(fifth row, left) is transformed from the $E^{\rm typical}_{\ell m}$ of a
randomly selected (post burn-in) sample (compare Eq. (13)):

$\kappa^{\rm typical}_{\ell m}=-\sqrt{\frac{\ell(\ell+1)}{(\ell-1)(\ell+2)}}\
E^{\rm typical}_{\ell m}$ (34)

while the signed S/N ratio map (fifth row, right) is defined as

$(S/N)_{\kappa}\equiv\frac{\kappa^{\rm
typical}(\theta,\phi)}{\sigma_{\kappa}(\theta,\phi)}$ (35)

where $\sigma_{\kappa}$ is the standard deviation of the convergence maps at
the relevant pixel.

Comparing the mean shear maps to the fiducial map, we see that Almanac samples
can properly recover the underlying structure of shears where the input galaxy
ellipticity maps have observations. For the masked areas, where observations
are not present, the mean maps display suppression of structure for small
scales. This is expected from the averaging of the samples. Each sample is
characterized by a power spectrum that contains small-scale modes, so the
sampled maps have fine structures even in the masked regions. However, in the
masked regions, the small-scale structure is not well constrained, so
different samples have different structures there, and these tend to average
out when the mean map is computed. The masked regions also have higher
posterior uncertainty, and this is reflected in the uncertainty maps shown in
the fourth row of Figs. 8 and 9.

Finally, we can use Eqs. (11) and (12) to obtain typical maps for the
tomographic projected lensing potential (or the projected Newtonian
gravitational potential, ${\tilde{\Psi}}{}$, as defined by Eq. (2)) for the
full sky. These maps are shown in Fig. 10, where they are compared to the
ground-truth obtained from the fiducial simulations (same underlying signal as
the data but full-sky and noiseless). We overlay the mask (white semi-
transparent region in the maps); this allows us to observe that Almanac’s
inferred gravitational field is in excellent agreement with fiducial field in
unmasked regions, while, as expected, the inferred map displays an overall
similar structure under the mask but with different locations of voids and
walls than that found in the fiducial map.

We perform two main consistency checks in the recovered maps, one for the
typical maps and another for the mean shear maps. The first test (see Fig. 11)
compares the underlying signal realisation from the fiducial maskless and
noiseless simulation with the pseudo-${\sf{C}}_{\ell}$ estimate of the typical
shear maps. We see that it is a typical map that best reflects the underlying
two point statistics. As expected, the agreement between the two is better for
the largest scales, while the typical map scatters around the fiducial value
for the smaller scales (due to the missing information resulting from the
mask). By contrast, the mean map will be over-smoothed in regions of low
signal-to-noise ratio. One could also consider the statistics of the maximum a
posteriori map, but the funnel shape of the posterior forces this to $a_{\ell
m}=0$, and in fact the small parameter volume in the neighbourhood of the
maximum means that even the maximum a posteriori sample is far from this
point. The second test (see Fig. 12) checks the distribution of residuals
between the shear maps samples and the fiducial shear map

$\mathcal{R}(\gamma_{i})\equiv\frac{1}{\sigma_{\gamma_{i}}}\left(\bar{\gamma}_{i}-\gamma^{\rm
fiducial}_{i}\right)$ (36)

where $\bar{\gamma}_{i}$ is the mean map for the samples of shear component i;
these are expected to follow a standard normal distribution. The figure shows
the distributions of the residuals for all recovered shear component maps, and
a standard normal for comparison. The small difference seen for the lower
redshift tomographic bins can be explained by their lower signal-to-noise
affecting the ability of the method to accurately estimate these maps.
However, the difference is below 1-$\sigma$ and does not raise concerns about
the maps shown in Figs. 8 and 9.

###  6.3 Convergence diagnostics

In this section we describe the convergence diagnostics used to assess the
validity of the chains (and hence the validity of the results shown in the
previous sections). Assessing convergence of a high-dimensional Monte Carlo
chain is a notoriously difficult problem. We examined several convergence
diagnostics on the individual chains, probing different aspects of the
convergence: a trace of the log posterior, the Fraction of Missing
Information, the correlation length of the chains, Hanson’s statistic, and the
Gelman-Rubin test. Further details are in Sellentin, Loureiro et al., in prep.

Fig. 13 shows the negative log posterior trace and its histogram for the
combined chains. The trace (upper plot) shows no obvious correlations, and
both chains show visually similar behaviour. This first test demonstrates that
both chains are consistently probing the same underlying distribution. The
histogram of posterior values (bottom plot) also shows a well-behaved
distribution, suggesting that burn-in was properly removed for both chains and
that both chains burned-in to the same region of the underlying posterior.

(a)

(b)

Figure 13: Diagnostics for the log posterior for the two chains (both
individually and combined). (a) The log posterior trace for the combined
(i.e., concatenated) chain (blue region is chain 1, orange region is chain 2).
(b) Combined (pink) and individual chain (blue and orange) histograms of the
log posterior values.

As a second test, we calculate the Fraction of Missing Information (FMI;
Betancourt 2016) 222This statistic has also been called the Bayesian Fraction
of Missing Information (BFMI); however, the statistic is not specifically
Bayesian and hence the shorter name is preferred.. This diagnostic assesses if
the sampler is correctly exploring the full range of total energy submanifolds
within phase space. It is defined as (Betancourt, 2016, 2017)

$\mathrm{FMI}=\frac{\sum_{n=1}^{N}(E_{n}-E_{n-1})^{2}}{\sum_{n=0}^{N}(E_{n}-\bar{E})^{2}},$
(37)

where $E_{n}$ is the total energy of the Hamiltonian system at the nth sample
in the chain and $\bar{E}$ is the mean energy of the Hamiltonian system for
the whole chain. In our experience with verifying against other convergence
diagnostics (correlation lengths, acceptance ratio, and the Hanson statistic),
an FMI $>0.7$, coupled with similar histograms in Fig. 14, indicates that the
sampler’s hyperparameters and the coordinate system used are such that the
sampler is properly exploring the posterior distribution. At the same time, an
FMI $<0.7$ is associated with poor convergence. Helpfully, the FMI can
diagnose bad behaviour with fewer samples than other statistics. For the
chains presented in Sec. 6, we achieve a FMI of 0.87 for chain 1 and 0.84 for
chain 2.

The requirement that FMI should be close to unity can be visualised by looking
at the histograms of $E=E_{n}-\bar{E}$ and $\Delta E=E_{n}-E_{n-1}$ shown in
Fig. 14 for both chains individually. Ideally the histograms of $P(\Delta E)$
and $P(E)$ would be identical, as this would demonstrate that the sampler has
correctly explored the full range of available energies within phase space.
Even with a very low S/N for the $B$-modes, Fig. 14 shows that the coordinate
transformation discussed in Sec. 4 allows the HMC sampler to properly explore
phase space (Betancourt, 2016, 2017).

As a third test we compute the correlation length of individual chains; this
is defined, for each parameter, to be the length at which the auto-correlation
of the parameter values in the chain drops below a threshold of $0.1$ for the
first time. We find the chains have a median correlation length of $3042$ and
$2585$, respectively, which corresponds to having a total of $165$ fully
independent samples in the chains. Both these tests demonstrate significant
convergence for the individual chains, meaning that the two chains can be
combined for more robust constraints.

As a fourth test we examine Hanson’s convergence statistic (Hanson, 2001;
Taylor et al., 2008),

${\cal
H}_{i}=\frac{\frac{1}{3N}\sum_{n}(x_{n,i}-\bar{x}_{i})^{3}\nabla_{x_{i}}{\psi}{}}{\frac{1}{N}\sum_{n=1}^{N}(x_{n,i}-\bar{x}_{i})^{2}}\,,$
(38)

where $x_{n,i}$ is the $n$th sample of the $i$th parameter. Our $\cal H$
equals that in Hanson (2001) and is the reciprocal of that in Taylor et al.
(2008). For a distribution with exponentially decaying tails, ${\cal H}_{i}$
is expected to tend to unity for all dimensions $i$, if convergence is
reached.

Taylor et al. (2008) states that in practice, ${\cal H}_{i}$ between 0.8 and
1.2 suggests good convergence in the chains. Fig. 15 shows this diagnostic
applied to the spherical harmonic dimensions of the posterior; this
demonstrates good convergence.

(a)

(b)

Figure 14: Fraction of Missing Information histograms for the two chains. The
two histograms display a similar distribution for both chains.

(a)

(b)

Figure 15: Hanson diagnostics for the spherical harmonic coefficients. The
quoted tolerance limits suggest good convergence (Taylor et al., 2008).

Finally, since we have two independent chains with different starting points,
we can perform a Gelman-Rubin test to further assess convergence. The results
of this test are shown in Fig. 16 for each of the angular power spectrum
dimensions. The median is very close to unity (1.0034) with a small tail to
higher values. Although excellent convergence is not achieved for all
dimensions, the overall distribution of marginalised angular power spectra
demonstrates a satisfactory convergence with $R<1.1$ for the vast majority of
parameters and $R<1.4$ when considering all dimensions. Hence, we are
confident that the chains have converged to the correct posterior
distribution.

Figure 16: Distribution of Gelman-Rubin $R$ statistics, for the power spectra,
for $\ell_{\rm max}=2048$. A very small number of $R$ values are larger than
the range shown.

##  7 Discussion

This work presents a Bayesian hierarchical model for joint inference of
angular power spectra and maps applied to simulated Stage IV tomographic
cosmic shear data distributed on the sphere. The method hierarchically models
the fields and their covariances, with a harmonic space parameterisation and a
Cholesky decomposition of the power spectra.

This Bayesian approach has the advantage of producing samples of the posterior
distribution of the power spectra and map coefficients, thereby avoiding the
simplifying assumption of gaussianity that is a feature of summary statistics
(and that may or may not be accurate). The posterior distribution has more
than 16.8 million parameters, and is thus a challenge for a standard Bayesian
approach to joint inference. Such a high-dimensional posterior requires an
efficient sampler; we have developed a Hamiltonian Monte Carlo sampler with
tuning stages to optimise the step size and the acceptance ratio (Sellentin,
Loureiro et al., in prep). To deal with the low S/N intrinsically found in the
$B$-modes (expected to be zero for the standard cosmological model and usual
extensions), we performed a coordinate transformation to both the spherical
harmonic coefficients, by re-scaling them to their variance, and to the
field’s covariance in harmonic space — described in Eqs. (20) and (22).

Although past works have applied similar methods to cosmic shear data (Alsing
et al., 2016; Alsing et al., 2017) and CMB lensing (Millea et al., 2021), this
is the first time this methodology has been developed for a complete lensing
sky, without the flat-sky approximation. Moving away from the flat-sky
approximation is crucial for analysing the next generation of cosmic shear
surveys due to the significant increase in sky coverage compared to the
previous generations (Lemos et al., 2017; Tutusaus et al., 2020; Matthewson &
Durrer, 2021). For the method presented here, we can use the increase in sky
coverage to our advantage: as more modes are available in the observed data,
Almanac can robustly infer the missing data under the mask and its angular
power spectra.

To test the methodology, we simulated a Euclid-like survey using gaussian
simulations to control the simulation’s signal, noise, and underlying angular
power spectra. The simulations contained two tomographic bins with the
expected Euclid mask, galaxy density, and shape noise. We then analysed the
cosmic shear simulated data with two chains using different starting points
over the multipole range $4\leq\ell\leq 2048$. From this analysis, we inferred
the angular power spectra of cosmic shear $E$\- and $B$-modes from the
simulated data. Here, we highlight the highly non-gaussian posterior
distributions shown in Fig. 6 for the largest scales and how, even for scales
larger than the ones available in data, we can still accurately infer the
angular power spectra of the $E$-modes. Although the largest scales in the
angular power spectra of cosmic shear carry little cosmological information,
these scales allow us to reconstruct the inferred denoised shear maps, as
shown in Sect. 6.2. We note, however, that such scales are extremely valuable
for galaxy clustering studies with primordial non-gaussianities (Andrews et
al., 2022; Baumann & Green, 2022). We will study the effects of the non-
gaussianities on the angular power spectra in a future work.

Since the method presented here also infers the spherical harmonic
coefficients for the underlying signal, we can convert these samples into
recovered signal maps for cosmic shear, convergence, and lensing potential.
The maps we recover in each sample contain information inferred under areas
where the telescope has not collected data. When naively analysing typical
maps in the chains, we can recover a significant amount of information, but
not sufficient to precisely recover the full signal for small scales (missing
due to the mask). Nevertheless, in future works we will explore how the mass
maps we recovered, under a gaussian signal assumption, compare to other
mainstream mass mapping methods in the literature (Kaiser & Squires, 1993;
Marshall et al., 2002; Lanusse et al., 2016; Chang et al., 2018; Price et al.,
2021; Jeffrey et al., 2021).

Differently from (pseudo) angular power spectrum estimators (Alonso et al.,
2019; Nicola et al., 2021; Chon et al., 2004) or correlation function
estimators (Jarvis et al., 2004; Kilbinger et al., 2014), here we are probing
the full posterior distribution of full-sky $E$\- and $B$-modes
${\sf{C}}_{\ell}$. This allows us to inspect for systematic contamination in
parity-violating and $B$-modes. For point estimates, the mask can make it
difficult to disambiguate between $E$-mode power and $B$-mode power. In some
analyses this is described as ‘$EB$ leakage’ for which various solutions are
proposed (Schneider et al., 2010; Leistedt et al., 2017; Liu et al., 2019;
Ghosh et al., 2021). With Almanac, on the other hand, we have the full
posterior and such leakage does not occur. The ambiguity may simply show up as
correlation between the $E$ and $B$ mode power spectra (and wider error bars
on each power spectrum element individually), and we typically find the
marginal distributions of $B$-mode power spectra to have a mode at zero, as we
expect. Only if these posteriors are converted to point estimates would one
expect such $EB$ leakage to occur.

In comparison with other BHMs for field-level inference, such as Borg (Jasche
& Wandelt, 2013; Jasche & Lavaux, 2019; Porqueres et al., 2021, 2022), Almanac
has the advantage of being agnostic to the underlying gravitational model, so
the map and power spectra samples may be used to test many cosmological
models. However, without a specific gravity model included, Almanac does not
sample the initial conditions of the Universe. Other approaches use a log-
normal prior, but do not yet analyse full-sky tomography (Fiedorowicz et al.,
2022; Boruah et al., 2022).

In terms of computational resources for future Euclid and LSST applications,
the method is limited by the spherical harmonic transformations performed
using libsharp2 (Reinecke & Seljebotn, 2013; Ishioka, 2018). For the set-up we
presented in this work, we are able to produce 100 samples from the posterior
distribution in less than 40 minutes on an Apple M1 Pro laptop333Note,
however, that the analysis was not performed on a laptop, but on local cluster
computers.. Future investigations of parallelising the posterior calculations
over redshift tomographic bins could potentially provide a significant speed-
up in the calculations.

Future works will investigate extensions of Almanac to non-gaussian
likelihoods and to the extraction of cosmological parameters via different
approaches, including naively using the marginalised samples and their
covariances for inference. We will also explore the use of the Almanac to
galaxy clustering and to joint analysis of cosmic shear, galaxy clustering,
and CMB.

## Acknowledgments

We thank the support staff of Leiden University’s ALICE High Performance
Computing infrastructure and UCL’s Hypatia cluster (especially Edd Edmondson),
Martin Reinecke for assistance with libsharp (Reinecke & Seljebotn, 2013), and
Malak Olamaie, Sreekumar Balan, Mike Hobson, Joris Bierkens and Justin Alsing
for useful discussions. Almanac was created by making significant
modifications to a software package originally written by Sreekumar Balan and
we record our gratitude to him for making this code available to us. The work
we present here would not have been possible without the following packages
and software: NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020),
Matplotlib (Hunter, 2007), HEALPix (Gorski et al., 1999), libsharp 2 (Ishioka,
2018), GetDist (Lewis, 2019), Pandas (Wes McKinney, 2010), AstroPy (Astropy
Collaboration et al., 2018), and Jupyter Lab (Kluyver et al., 2016). This work
used computing equipment funded by the Research Capital Investment Fund (RCIF)
provided by UKRI, and was partially funded by the UCL Cosmoparticle
Initiative.

## References

  * Alonso et al. (2019) Alonso D., Sanchez J., Slosar A., LSST Dark Energy Science Collaboration 2019, MNRAS, 484, 4127
  * Alsing et al. (2016) Alsing J., Heavens A., Jaffe A. H., Kiessling A., Wandelt B., Hoffmann T., 2016, MNRAS, 455, 4452
  * Alsing et al. (2017) Alsing J., Heavens A., Jaffe A. H., 2017, MNRAS, 466, 3272
  * Amon & Efstathiou (2022) Amon A., Efstathiou G., 2022, MNRAS, 516, 5355
  * Amon et al. (2022) Amon A., et al., 2022, Phys. Rev. D, 105, 023514
  * Andrews et al. (2022) Andrews A., Jasche J., Lavaux G., Schmidt F., 2022, arXiv e-prints, p. arXiv:2203.08838
  * Asgari et al. (2019) Asgari M., et al., 2019, A&A, 624, A134
  * Asgari et al. (2021) Asgari M., et al., 2021, A&A, 645, A104
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Bardeen et al. (1986) Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15
  * Bartelmann & Maturi (2017) Bartelmann M., Maturi M., 2017, Scholarpedia, 12, 32440
  * Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
  * Baumann & Green (2022) Baumann D., Green D., 2022, J. Cosmology Astropart. Phys, 2022, 061
  * Betancourt (2016) Betancourt M., 2016, arXiv e-prints, p. arXiv:1604.00695
  * Betancourt (2017) Betancourt M., 2017, arXiv e-prints, p. arXiv:1701.02434
  * Betancourt & Girolami (2013) Betancourt M. J., Girolami M., 2013, arXiv e-prints, p. arXiv:1312.0906
  * Blazek et al. (2015) Blazek J., Vlah Z., Seljak U., 2015, J. Cosmology Astropart. Phys, 2015, 015
  * Bond et al. (1998) Bond J. R., Jaffe A. H., Knox L., 1998, Phys. Rev. D, 57, 2117
  * Boruah et al. (2022) Boruah S. S., Rozo E., Fiedorowicz P., 2022, MNRAS, 516, 4111
  * Bunn (2008) Bunn E. F., 2008, arXiv e-prints, p. arXiv:0811.0111
  * Bunn et al. (2003) Bunn E. F., Zaldarriaga M., Tegmark M., de Oliveira-Costa A., 2003, Phys. Rev. D, 67, 023501
  * Casella & George (1992) Casella G., George E. I., 1992, The American Statistician, 46, 167
  * Castro et al. (2005) Castro P. G., Heavens A. F., Kitching T. D., 2005, Phys. Rev. D, 72, 023516
  * Catelan et al. (2001) Catelan P., Kamionkowski M., Blandford R. D., 2001, MNRAS, 320, L7
  * Chang et al. (2018) Chang C., et al., 2018, MNRAS, 475, 3165
  * Chisari et al. (2019) Chisari N. E., et al., 2019, ApJS, 242, 2
  * Chon et al. (2004) Chon G., Challinor A., Prunet S., Hivon E., Szapudi I., 2004, MNRAS, 350, 914
  * Codis et al. (2018) Codis S., Jindal A., Chisari N. E., Vibert D., Dubois Y., Pichon C., Devriendt J., 2018, MNRAS, 481, 4753
  * Colombo et al. (2022) Colombo L. P. L., et al., 2022, arXiv, p. 2208.14276
  * DES Collaboration et al. (2022) DES Collaboration et al., 2022, Phys. Rev. D, 105, 023520
  * Eriksen et al. (2004) Eriksen H. K., et al., 2004, ApJS, 155, 227
  * Euclid Collaboration et al. (2020) Euclid Collaboration et al., 2020, A&A, 642, A191
  * Euclid Collaboration et al. (2022) Euclid Collaboration et al., 2022, A&A, 662, A112
  * Fiedorowicz et al. (2022) Fiedorowicz P., Rozo E., Boruah S. S., 2022, arXiv e-prints, p. arXiv:2210.12280
  * Geman & Geman (1984) Geman S., Geman D., 1984, IEEE - Transactions on Pattern Analysis and Machine Intelligence, pp 721–741
  * Ghosh et al. (2021) Ghosh S., Delabrouille J., Zhao W., Santos L., 2021, J. Cosmology Astropart. Phys, 2021, 036
  * Gorski et al. (1999) Gorski K. M., Wandelt B. D., Hansen F. K., Hivon E., Banday A. J., 1999, arXiv e-prints, pp astro–ph/9905275
  * Gupta (1999) Gupta A.K.and Nagar D., 1999, Matrix Variate Distributions. Chapman & Hall/CRC, London
  * Hajian (2007) Hajian A., 2007, Phys. Rev. D, 75, 083525
  * Hall (2021) Hall A., 2021, MNRAS, 505, 4935
  * Hanson (2001) Hanson K. M., 2001, in Sonka M., Hanson K. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4322, Medical Imaging 2001: Image Processing. pp 456–467, doi:10.1117/12.431119
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hartlap et al. (2009) Hartlap J., Schrabback T., Simon P., Schneider P., 2009, A&A, 504, 689
  * Hastings (1970) Hastings W. K., 1970, Biometrika, 57, 97
  * Heavens (2009) Heavens A., 2009, Nuclear Physics B Proceedings Supplements, 194, 76
  * Heavens et al. (2000) Heavens A., Refregier A., Heymans C., 2000, MNRAS, 319, 649
  * Heymans et al. (2006) Heymans C., White M., Heavens A., Vale C., van Waerbeke L., 2006, MNRAS, 371, 750
  * Heymans et al. (2013) Heymans C., et al., 2013, MNRAS, 432, 2433
  * Heymans et al. (2021) Heymans C., et al., 2021, A&A, 646, A140
  * Hikage et al. (2019) Hikage C., et al., 2019, PASJ, 71, 43
  * Hoffman & Ribak (1992) Hoffman Y., Ribak E., 1992, ApJ, 384, 448
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ishioka (2018) Ishioka K., 2018, Journal of the Meteorological Society of Japan. Ser. II, 96, 241
  * Jarvis (2015) Jarvis M., 2015, TreeCorr: Two-point correlation functions, Astrophysics Source Code Library, record ascl:1508.007 (ascl:1508.007)
  * Jarvis et al. (2004) Jarvis M., Bernstein G., Jain B., 2004, MNRAS, 352, 338
  * Jasche & Kitaura (2010) Jasche J., Kitaura F. S., 2010, MNRAS, 407, 29
  * Jasche & Lavaux (2019) Jasche J., Lavaux G., 2019, A&A, 625, A64
  * Jasche & Wandelt (2013) Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894
  * Jaynes (2003) Jaynes E. T., 2003, Probability theory: The logic of science. Cambridge University Press, Cambridge
  * Jeffrey et al. (2021) Jeffrey N., et al., 2021, MNRAS, 505, 4626
  * Jewell et al. (2009) Jewell J. B., Eriksen H. K., Wandelt B. D., O’Dwyer I. J., Huey G., Górski K. M., 2009, ApJ, 697, 258
  * Joachimi & Bridle (2010) Joachimi B., Bridle S. L., 2010, A&A, 523, A1
  * Joachimi et al. (2015) Joachimi B., et al., 2015, Space Sci. Rev., 193, 1
  * Kaiser (1992) Kaiser N., 1992, ApJ, 388, 272
  * Kaiser & Squires (1993) Kaiser N., Squires G., 1993, ApJ, 404, 441
  * Kiessling et al. (2015) Kiessling A., et al., 2015, Space Sci. Rev., 193, 67
  * Kilbinger (2015) Kilbinger M., 2015, Reports on Progress in Physics, 78, 086901
  * Kilbinger et al. (2014) Kilbinger M., Bonnett C., Coupon J., 2014, athena: Tree code for second-order correlation functions, Astrophysics Source Code Library, record ascl:1402.026 (ascl:1402.026)
  * Kirk et al. (2015) Kirk D., et al., 2015, Space Sci. Rev., 193, 139
  * Kluyver et al. (2016) Kluyver T., et al., 2016, in Loizides F., Scmidt B., eds, Positioning and Power in Academic Publishing: Players, Agents and Agendas. IOS Press, Netherlands, pp 87–90, https://eprints.soton.ac.uk/403913/
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lanusse et al. (2016) Lanusse F., Starck J. L., Leonard A., Pires S., 2016, A&A, 591, A2
  * Larson et al. (2007) Larson D. L., Eriksen H. K., Wandelt B. D., Górski K. M., Huey G., Jewell J. B., O’Dwyer I. J., 2007, ApJ, 656, 653
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Leclercq et al. (2015) Leclercq F., Jasche J., Wandelt B., 2015, J. Cosmology Astropart. Phys, 2015, 015
  * Leistedt et al. (2017) Leistedt B., McEwen J. D., Büttner M., Peiris H. V., 2017, MNRAS, 466, 3728
  * Lemos et al. (2017) Lemos P., Challinor A., Efstathiou G., 2017, J. Cosmology Astropart. Phys, 2017, 014
  * Lemos et al. (2021) Lemos P., et al., 2021, MNRAS, 505, 6179
  * Lewis (2005) Lewis A., 2005, Phys. Rev. D, 71, 083008
  * Lewis (2019) Lewis A., 2019, arXiv e-prints, p. arXiv:1910.13970
  * Liu et al. (2019) Liu H., Creswell J., von Hausegger S., Naselsky P., 2019, Phys. Rev. D, 100, 023538
  * Loureiro et al. (2022) Loureiro A., et al., 2022, A&A, 665, A56
  * MacCrann et al. (2015) MacCrann N., Zuntz J., Bridle S., Jain B., Becker M. R., 2015, MNRAS, 451, 2877
  * Mandelbaum et al. (2018) Mandelbaum R., et al., 2018, PASJ, 70, S25
  * Marshall et al. (2002) Marshall P. J., Hobson M. P., Gull S. F., Bridle S. L., 2002, MNRAS, 335, 1037
  * Matthewson & Durrer (2021) Matthewson W. L., Durrer R., 2021, J. Cosmology Astropart. Phys, 2021, 027
  * Metropolis et al. (1953) Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087
  * Millea et al. (2021) Millea M., et al., 2021, ApJ, 922, 259
  * Munshi et al. (2008) Munshi D., Valageas P., van Waerbeke L., Heavens A., 2008, Phys. Rep., 462, 67
  * Neal (2000) Neal R. M., 2000, arXiv e-prints, p. physics/0009028
  * Neal (2011) Neal R., 2011, in Brooks S., Gelman A., Jones G. L., Meng X.-L., eds, , Handbook of Markov Chain Monte Carlo. Chapman & Hall/CRC, pp 113–162, doi:10.1201/b10905
  * Nicola et al. (2021) Nicola A., García-García C., Alonso D., Dunkley J., Ferreira P. G., Slosar A., Spergel D. N., 2021, J. Cosmology Astropart. Phys, 2021, 067
  * Percival et al. (2022) Percival W. J., Friedrich O., Sellentin E., Heavens A., 2022, MNRAS, 510, 3207
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Porqueres et al. (2021) Porqueres N., Heavens A., Mortlock D., Lavaux G., 2021, MNRAS, 502, 3035
  * Porqueres et al. (2022) Porqueres N., Heavens A., Mortlock D., Lavaux G., 2022, MNRAS, 509, 3194
  * Price et al. (2021) Price M. A., McEwen J. D., Pratley L., Kitching T. D., 2021, MNRAS, 500, 5436
  * Racine et al. (2016) Racine B., Jewell J. B., Eriksen H. K., Wehus I. K., 2016, ApJ, 820, 31
  * Reinecke & Seljebotn (2013) Reinecke M., Seljebotn D. S., 2013, A&A, 554, A112
  * Rossi (2012) Rossi G., 2012, MNRAS, 421, 296
  * Schneider (2005) Schneider P., 2005, arXiv e-prints, pp astro–ph/0509252
  * Schneider et al. (2010) Schneider P., Eifler T., Krause E., 2010, A&A, 520, A116
  * Secco et al. (2022) Secco L. F., et al., 2022, Phys. Rev. D, 105, 023515
  * Sellentin & Heavens (2016) Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132
  * Sellentin & Heavens (2018) Sellentin E., Heavens A. F., 2018, MNRAS, 473, 2355
  * Sellentin et al. (2018) Sellentin E., Heymans C., Harnois-Déraps J., 2018, MNRAS, 477, 4879
  * Taylor et al. (2008) Taylor J. F., Ashdown M. A. J., Hobson M. P., 2008, MNRAS, 389, 1284
  * Taylor et al. (2021) Taylor P. L., et al., 2021, The Open Journal of Astrophysics, 4, 6
  * Tutusaus et al. (2020) Tutusaus I., et al., 2020, A&A, 643, A70
  * Upham et al. (2022) Upham R. E., et al., 2022, A&A, 660, A114
  * Van Waerbeke et al. (2013) Van Waerbeke L., et al., 2013, MNRAS, 433, 3373
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wallis et al. (2022) Wallis C. G. R., Price M. A., McEwen J. D., Kitching T. D., Leistedt B., Plouviez A., 2022, MNRAS, 509, 4480
  * Wandelt et al. (2004) Wandelt B. D., Larson D. L., Lakshminarayanan A., 2004, Phys. Rev. D, 70, 083511
  * Wes McKinney (2010) Wes McKinney 2010, in Stéfan van der Walt Jarrod Millman eds, Proceedings of the 9th Python in Science Conference. pp 56 – 61, doi:10.25080/Majora-92bf1922-00a