# Bayesian mass mapping with weak lensing data using karmma – validation with
simulations and application to Dark Energy Survey Year 3 data

Supranta S. Boruah<EMAIL_ADDRESS>Department of Astronomy and Steward
Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85719, USA
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia,
PA 19104, USA  Pier Fiedorowicz  Department of Physics, University of Arizona,
1118 E. Fourth Street, Tucson, AZ, 85721, USA  Lawrence Livermore National
Laboratory, Livermore, CA 94550, USA  Eduardo Rozo  Department of Physics,
University of Arizona, 1118 E. Fourth Street, Tucson, AZ, 85721, USA

###### Abstract

We update the field-level inference code karmma to enable tomographic forward-
modelling of shear maps. Our code assumes a lognormal prior on the convergence
field, and properly accounts for the cross-covariance in the lensing signal
across tomographic source bins. We use mock weak lensing data from $N$-body
simulations to validate our mass-mapping forward model by comparing our
posterior maps to the input convergence fields. We find that karmma produces
more accurate reconstructions than traditional mass-mapping algorithms.
Moreover, the karmma posteriors reproduce all statistical properties of the
input density field we tested — one- and two-point functions, and the peak and
void number counts — with $\approx 10\%$ accuracy. Our posteriors exhibit a
small bias that increases with decreasing source redshift, but these biases
are small compared to the statistical uncertainties of current (DES) cosmic
shear surveys. Finally, we apply karmma to Dark Energy Survey Year 3 (DES-Y3)
weak lensing data, and verify that the two point shear correlation function
$\xi_{+}$ is well fit by the correlation function of the reconstructed
convergence field. This is a non-trivial test that traditional mass mapping
algorithms fail. The code is publicly available at
https://github.com/Supranta/KaRMMa.git. karmma DES-Y3 mass maps are publicly
available at https://zenodo.org/records/10672062.

##  I Introduction

Most standard analyses of cosmological data use 2-point summary statistics
such as correlation functions and power spectra for extracting information
from these data sets [e.g, 1, 2]. However, we lose valuable non-Gaussian
information while compressing the full data set to 2-point summary statistics.
Field-based inference, where one extracts information from the full field, is
emerging as a powerful alternative for optimally extracting information from
cosmological data sets across a broad variety of probes, including CMB lensing
[3, 4, 5], integrated Sachs-Wolfe (ISW) effect [6], galaxy clustering [7, 8],
peculiar velocity [9, 10, 11] and weak lensing [12, 13, 14].

In this paper, we focus on field-based inference of weak lensing data. Such an
analysis requires an accurate forward model of the observed shear field. Once
we have an accurate forward model for the weak lensing data, we can use it to
reconstruct the unknown dark matter distribution of the Universe. Indeed,
reconstructed mass maps from weak lensing surveys are crucial for extracting
non-Gaussian information [e.g, 15], for studying voids [16], and for cross-
correlation studies [17, 18]. As such, they are an essential data product for
Stage-IV weak lensing surveys. However, the reconstruction of the unseen dark
matter distribution from weak lensing shear data constitutes a non-trivial
inverse problem. While standard inversion methods such as Kaiser-Squires
inversion [19] are widely used, they suffer a variety of problems. For
example, these methods lead to biased reconstruction in the presence of a
survey mask. Furthermore, in the presence of shape noise the reconstructed
mass maps do not have the expected statistical properties. For instance, the
two-point functions of the recovered maps are inconsistent with direct
measurements of the two point shear correlation function (see Figure 10).

Forward modelled mass-map reconstruction naturally address these difficulties.
For example, by forward-modelling the convergence field in the masked region,
we get mass maps free from masking effects [20, 21]. Likewise, depending on
the accuracy of the forward model, we are able to reconstruct mass maps with
the correct 2-point and one-point function.

A number of different weak lensing field-based analysis methods have been
proposed in the literature, though none so far have been applied to data.
Running approximate numerical structure formation simulations is likely to be
the most accurate forward model. This is the approach followed by [22, 13],
who used Lagrangian Perturbation Theory (LPT) to forward model the cosmic
shear data within the borg Bayesian forward modeling framework [23, 7].
However, this approach is computationally expensive and analyses are often
limited to small areas with a limited resolution.

Because of these difficulties, a number of alternative field-based inference
methods have tried to directly model the observed 2-dimensional projected
fields, thereby bypassing the need to generate a three-dimensional matter
density field. For example, [24, 25] developed a Bayesian hierarchical model
that simultaneously inferred the mass distribution and power spectrum of the
weak lensing maps in tomographic bins. However, their method assumed
Gaussianity of the field and only worked on flat sky. [26, 27] proposed a
similar algorithm, almanac, which performs wide-field mass mapping on a
sphere, but still assumes Gaussianity. Similarly, [16] sampled mass maps with
DES-Y3 data using the constrained realization approach [28, 29, 30] assuming
Gaussianity of the convergence field.

Since the late-time density field is highly non-Gaussian, it is important to
include the non-Gaussianity of the density/convergence field in the forward
model. To that end, [21, 31] introduced the karmma algorithm. In karmma , the
convergence field is modeled as a lognormal random field, enabling us to by-
pass the need to run expensive numerical simulations to reconstruct the matter
density field, while simultaneously improving upon the Gaussian assumption
made in other works (e.g. almanac). More recently [14] introduced miko, a
flat-sky field-level inference code that also utilizes the same lognormal
prior for the convergence field. The lognormal model has been shown to be a
good approximation of the late-time non-Gaussian convergence field [32, 33].
Indeed, in our earlier works we demonstrated that the karmma posteriors from
mock observations of simulated shear maps accurately reproduced a broad range
of statistical properties of the corresponding input convergence fields.
However, the algorithm had a severe limitation: it was only able to forward
model a single tomographic source bin at a time. Here, we extend the karmma
algorithm using the methodology of [12] to enable tomographic mass mapping on
a sphere while modeling the non-Gaussianity of the convergence field. We
validate the updated code by running it on mock weak lensing data based on
$N$-body simulations. We find that the mass maps reconstructed with karmma
have smaller errors than traditional mass mapping approaches. Further, the
statistical properties of the posterior maps are consistent with those from
numerical simulations. After validating our method with mock simulations, we
apply karmma to the Dark Energy Survey (DES) Year 3 (Y3) weak lensing data. We
make our maps publicly available at https://zenodo.org/records/10672062. To
the best of our knowledge, these are the first publicly-available Bayesian
mass maps reconstructed with weak lensing data on the curved sky.

This paper is structured as follows: in section II, we present the simulations
and the weak lensing data used in this work. In section III, we present the
main features of karmma . Then in section IV, we validate karmma with mock
weak lensing simulations before presenting the Bayesian mass maps produced
from DES-Y3 data in section V. Finally, we conclude in section VI. In Appendix
A, we show the results of our tests with lognormal mocks.

##  II Data

###  II.1 Mock maps from simulations

In this section, we describe the mock weak lensing catalogues used in section
IV to test the performance of our code. We use the publicly available
simulations of [34, hereafter T17]. T17 produced a suite of $108$ mock weak
lensing catalogues by ray-tracing $N$-body simulations run with cosmological
parameters $\Omega_{m}=0.279$, $\Omega_{b}=0.046$, $h=0.7$, $\sigma_{8}=0.82$
and $n_{s}=0.97$. Each mock catalogue consists of full-sky shear and
convergence maps on redshift shells separated by $150~{}h^{-1}$ Mpc. We use
healpix maps at a resolution of $N_{\text{side}}=4096$ for this work.

To generate simulated shear maps, we adopt the survey properties from the
DES-Y3 analysis [35, 36]. Specifically, we use four tomographic redshift bins.
The source redshift distribution of each bin is taken from [37], shown here in
Figure 1. The effective number density in each of the tomographic bins is
$[1.476,1.479,1.484,1.461]$ arcmin-2. Finally, we adopt the shape noise
estimate from the DES-Y3 shear catalog, $\sigma_{\epsilon}=0.261$ [38].

We produce tomographic maps of the quantity, $\mathcal{A}$, (here
$\mathcal{A}$ can be either the convergence or shear map) by taking a weighted
average of the maps at various redshift shells with a given redshift
distribution, $n(z)$, via

$\mathcal{A}^{i}=\frac{\sum_{j=1}^{N_{\text{shells}}}n^{i}(z_{j})\mathcal{A}_{\text{shell}}^{j}}{\sum_{j=1}^{N_{\text{shells}}}n^{i}(z_{j})},$
(1)

where $n^{i}$ is the redshift distribution of the $i$-th tomographic bin and
$\mathcal{A}^{j}_{\text{shell}}$ is the map of the quantity $\mathcal{A}$ on
the $j$-th redshift shell. Note that the finite redshift resolution leads to a
systematic error of $\mathcal{O}(\lesssim 5\%)$ on the power spectrum of the
convergence field. Consequently, for the purposes of testing our code, we use
the power spectrum measured from the T17 simulations, rather than a theory
power spectrum computed using Boltzmann codes such as camb [39] or class [40].
We change this in section V while applying to DES-Y3 data.

Our ‘observed’ shear map is a noisy realization of the true shear.
Specifically, we produce our noisy shear maps at a higher resolution of
$N_{\rm side}=1024$ and then downgrade to a resolution of $N_{\rm side}=256$
to mimic the impact of pixelization. The number of sources in each sky pixel
is a Poisson random draw from the mean source density of the appropriate
tomographic bin in the higher resolution maps. Note that this ignores source
clustering which is an important systematic effect while considering non-
Gaussian information [41]. We add to each shear component a Gaussian shape
noise with zero mean and standard deviation, $\sigma_{\epsilon}/\sqrt{N_{p}}$,
where, $N_{p}$ is the number of source galaxies sampled in that pixel. A
resolution of $N_{\text{side}}=256$ corresponds to a pixel size with angular
resolution of $\sim 13$ arcmin.

###  II.2 DES-Y3 weak lensing data

Figure 1: The source redshift distribution, $n(z)$, for the DES-Y3 used in
this work.

In section V, we run karmma on the Dark Energy Survey Year 3 (DES-Y3) weak
lensing data to produce Bayesian mass maps from weak lensing data. We use the
DES-Y3 metacalibration [42, 43] shape catalogue [38]. These galaxies were
divided into the $4$ tomographic bins according to their estimated photometric
redshifts [37]. Then we construct a healpix map of estimated shear for each of
the tomographic bin by taking the weighted average over the galaxy
ellipticities, $\epsilon$, in each pixel as,

$\gamma^{\alpha}_{p}=\frac{\sum_{i\in p}w_{i}\epsilon_{i,\alpha}}{\sum_{i\in
p}w_{i}},\alpha=1,2,$ (2)

where $w_{i}$ are the per-galaxy inverse variance weights, $p$ is the pixel
label, and $\alpha$ denotes the two polarizations. The estimated variance of
each shear component in the pixel $p$ is given as [44, 45],

$\sigma^{2}_{\epsilon,p}=\frac{\sum_{i\in
p}w^{2}_{i}(e^{2}_{i,1}+e^{2}_{i,2})/2}{(\sum_{i\in p}w_{i})^{2}}.$ (3)

Figure 2: Illustration of the various steps involved in the karmma forward
model. Here, SHT refers to the Spherical harmonic transform, KS refers to
Kaiser-Squires. Note that each of the steps illustrated above are
differentiable, thus allowing us to compute the gradient with respect to the
input variables – a necessary feature for using HMC. The maximum $\ell$ value
of the spherical harmonic modes and the healpix map resolution of the maps is
shown at each step. See section III for details.  Figure 3: An illustration of
the properties of the karmma posteriors for the first tomographic redshift bin
in our simulations. The input convergence map used to generate our mock
observations is shown in the top left panel. The centre and right panels on
the top row shows two randomly chosen karmma map samples. The mean map and the
corresponding pixel-by-pixel RMS uncertainty are shown in the bottom left and
bottom centre panels, respectively. The bottom right panel shows the signal to
noise ratio computed from the sample mean and standard deviation.

##  III Mass map inference with KaRMMa

karmma models the convergence field as realizations of a lognormal random
field. The observed shears are then forward-modeled from the resulting
convergence fields. We determine the posterior distribution of convergence
maps by sampling the posterior using Hamiltonian Monte Carlo (HMC) methods. In
this way, rather than providing a single “best-fit” convergence, we recover
the full posterior distribution of convergence maps consistent with the input
shear data. These posterior quantify the uncertainty in the reconstructed mass
maps. We describe this algorithm in further detail below.

Our lognormal assumption states that the convergence field $\kappa$ is a non-
linear transform of a Gaussian random field $y$ where

$\kappa=e^{y}-\lambda.$ (4)

The parameter $\lambda$ is called the shift parameter. At the resolutions of
this study ($N_{\rm side}=256$, or $\approx 13\ {\rm arcmin}$), a multivariate
lognormal model provides a good description of the convergence field [32, 33].

Figure 4: Same as Figure 3, but for the fourth tomographic bin. Compared to
Figure 3, the structures are produced with a higher SNR in the higher redshift
bins.

Here, we use the lognormal formalism developed in [21, 12]. Specifically, we
sample the spherical harmonics (or Fourier wave modes) of the $y$ maps.
Furthermore, following [12] we introduce new reparameterized variables,
${\bm{x}}$, defined via

$y_{lm}={\mathbf{L}}(l)x_{lm},$ (5)

where, $y_{lm}$ is the spherical harmonics of the $y$ variable defined in
equation (4) and ${\mathbf{L}}(l)$ is the Cholesky decomposition of a
$n_{\text{bin}}\times n_{\text{bin}}$ matrix constructed from the power
spectrum $C(l)$ at a wave mode $l$ in various tomographic bins
($n_{\text{bin}}$ being the number of tomographic bins).111Note that the
definition of the unit variance variables differ from the definition in [12]
where the $x_{lm}$ variables were related to the $y_{lm}$ through the
eigenvalue decomposition of $C(l)$.  This reparameterization de-correlates the
variables in different redshift bins such that all the elements of ${\bm{x}}$
have unit variance and no cross-correlation under the prior.222Note that
posterior samples of $x$ can be correlated due to the likelihood. As shown in
[12], this redefinition makes the sampling of mass maps more efficient. The
power spectrum and the shift parameters are computed at the T17 cosmological
parameters using their simulations.

Following [31], we filter our maps at $\ell_{\text{max}}=2N_{\text{side}}$ to
avoid aliasing. Note that the low-pass filter is applied on the $\kappa$ field
and not on the $y$ (or $x$) fields. Also note that because of the filtering,
the resulting $\kappa$ field is no longer a lognormal field.

Once we have the convergence field, the shear field is constructed using the
Kaiser-Squires (KS) relation. We perform our inference on the sphere, thus
requiring the curved sky version of the KS relation; the spherical harmonic
coefficients of the shear field are related to the spherical harmonic
coefficients of the $\kappa$ field via [46]

$\tilde{\gamma}_{lm}=\sqrt{\frac{(l-1)(l+2)}{l(l+1)}}\kappa_{lm}.$ (6)

The shear field is then obtained as a sum over the spin-2 spherical harmonics,

$\gamma({\bm{\theta}})=\sum_{lm}\tilde{\gamma}_{lm}~{}{}_{2}Y_{lm}({\bm{\theta}}).$
(7)

Since the convergence field is low-pass filtered, the shear field is also low-
pass filtered. By contrast, the observed shear field is not filtered, which
leads to some amount of model mis-specification. However, we have found that
given the resolution of our maps, this slight model mis-specification is
easily preferable to the aliasing that occurs when we do not filter the field
[31].

Finally, we use pixelized shear maps to perform the inference, so we must
account for the pixel window function in our forward model. The pixelization
operation is implemented in the harmonic space by multiplying each spherical
harmonic mode by the spherical transform of the pixel window function,

$\tilde{\gamma}_{lm}\rightarrow
W^{N_{\text{side}}}_{p}(l)\tilde{\gamma}_{lm},$ (8)

where, $W^{N_{\text{side}}}_{p}$ is the healpix pixel window function at
resolution of $N_{\text{side}}$. While in principle the pixel window function
for each healpix pixel is different, an average pixel window function can be
computed using the healpy function healpy.sphtfunc.pixwin. The need for
including the pixel window function in the forward model has previously been
noted in [21, 14].

Given a $\kappa$ map, the observed shear map is modeled as a noisy realization
of the true shear map, where the noise is set by the shape noise of the
survey. In particular, the noise is uncorrelated from pixel to pixel, so that

$\log
P(\gamma_{\text{obs}}|{\bm{x}})=\frac{\chi^{2}}{2}=\sum_{i=1}^{N_{\text{pix}}}\frac{[\gamma^{i}({\bm{x}})-\gamma^{i}_{\text{obs}}]^{2}}{2\sigma^{2}_{\epsilon,i}}.$
(9)

The summation is over all the pixels within the survey mask.

Given an ‘observed’ shear map, the posterior of ${\bm{x}}$ is simply

$P({\bm{x}}|\gamma_{\text{obs}})\propto P(\gamma_{\text{obs
}}|{\bm{x}})P({\bm{x}}).$ (10)

Because ${\bm{x}}$ are unit variance variables, the prior, $P({\bm{x}})$ is a
unit variance normal distribution.

The various steps involved in our forward model are illustrated in Figure 2.
Note that all the intermediate steps are differentiable allowing us to take
the gradient of the posterior with respect to the variables ${\bm{x}}$, as is
required to sample using HMC. We implement our forward model in the python
package pytorch [47] and use the probabilistic programming language pyro [48]
to sample the mass-maps. We use the No U-Turn Sampler [NUTS, 49] to sample
${\bm{x}}$ (or equivalently ${\bm{\kappa}}$).

We emphasize that in the current implementation, the transformation from the
$x$ variables to convergence are conditioned on the power spectrum and shift
parameter characterizing the lognormal prior. Both of these are cosmology
dependent, and therefore one must specify a cosmological model to sample from
these maps. One can, in principle, simultaneously sample the cosmological
parameters characterizing the prior as part of the inference process as
described in [12]. In this way, karmma can produce joint cosmology and mass
map posteriors. However, we have not yet characterized the shift parameter
$\lambda$ as a function of cosmology, nor have we validated the resulting
cosmological inference pipeline. For this reason, we postpone this
simultaneous inference process to future work.

Figure 5: Comparison of the residual errors in a randomly chose karmma mass
map sample (blue) and Kaiser-Squires mass map (red) created with the same
data. The different panels show the results in different tomographic bins. As
we can see from the figure, the residuals of the karmma mass maps are much
lower than that of the Kaiser-Squires mass maps, showing the better quality of
the karmma mass maps. The improvement is especially noteworthy for the lowest
redshift bin, where the signal-to-noise is the lowest.  Figure 6: Fractional
error between the two-point functions of the input simulation maps and those
estimated from the karmma posteriors. The bottom-left triangles shows the
correlation functions between the different tomographic bins, while the upper-
right triangles shows the corresponding power spectra computed using
pseudo-$C(l)$s. Each panel corresponds to a different cross-bin combination,
as labeled. The blue lines show the difference in the map statistics averaged
across the full 108 simulated maps. The grey bands shows the mean error in one
map as estimated using the karmma posteriors, while the blue band shows the
error on the mean estimated using all 108 simulations. The difference between
the simulated maps and our posteriors increases with decreasing redshift and
smaller scales, and is typically $\lesssim 5\%$. Figure 7: Comparison of non-
Gaussian statistics of the karmma maps (blue) and of the mock simulations (red
dashed). The different rows corresponds to the three non-Gaussian statistics
computed here – the PDF of $\kappa$ values (top panels), the peak counts
(middle panels) and the void counts (bottom panels). The different columns
show these statistics in the $4$ tomographic bins. As we can see, all these
non-Gaussian summary statistics are well-approximated by karmma at high
redshift, whereas at low-redshift, these summary statistics are significantly
biased due to the breakdown of the lognormal model.

##  IV Validating KaRMMa with mock simulations

We test the performance of karmma using the simulations presented in section
II.1. The multivariate lognormal prior requires the power spectrum and the
shift parameters as its input. We measure these quantities in the simulations,
averaging over all 108 simulations, and use these as the inputs to the
lognormal prior.333We remind the reader that the use of the simulated power
spectrum is necessary due to the small differences between the power spectrum
in the simulations and theory power spectra. When applied to data, we replace
the power spectrum by the appropriate theory power spectrum for our fiducial
cosmology. We then run karmma on all $108$ mock data sets created from the T17
data. This allows us to check how well the lognormal forward model performs on
mock weak lensing data from $N$-body simulations. Each karmma run requires
$\sim\mathcal{O}(50-100)$ CPU core hours to produce $100$ independent samples
on AMD Zen2 processors.

Figures 3 and 4 illustrate the properties of the karmma posteriors when run on
our simulated data set. The figures correspond to the mass maps from
tomographic redshift bins 1 and 4 respectively. The input mass map used to
generate the simulation is shown on the top left panel. The centre and right
panels on the top row show two randomly chosen samples for our posterior
distribution of maps. The mean and standard deviation of the posteriors are
shown in the bottom left and bottom centre panels respectively. The bottom
right panels shows the signal-to-noise ratio defined as,

$\text{SNR}=\frac{\text{Mean}(\kappa)}{\sigma[\kappa]}.$ (11)

By comparing the input convergence map to the karmma posterior and mean map,
we can see that karmma recovers the dark matter distribution of the input map.
We also see that the higher redshift bins have a higher signal-to-noise ratio,
as expected.

Next we assess the quality of the reconstructed mass maps by computing the
pixel-by-pixel errors in the output mass maps. The error between the
convergence value in the $i$-th pixel of a reconstructed karmma sample
(labelled $s$) and the simulation is

$\Delta\kappa_{i}=\kappa^{i}_{\text{KARMMA},s}-\kappa^{i}_{\text{sim}}.$ (12)

In Figure 5, we compare the histogram of the convergence residuals for two
different maps: 1) a randomly chosen map from our posterior sample; and 2) a
traditional Kaiser–Squires mass map. Both maps are generated using the same
input shear map. We find the errors in the karmma mass maps are always much
lower than those in Kaiser–Squires maps. The standard deviation of the karmma
residuals are $85\%$, $77\%$, $66\%$ and $60\%$ lower than the standard
deviations of the KS residuals in the $4$ tomographic bins respectively. The
improvement is especially noticeable for the lowest redshift bin, where the
signal-to-noise is the lowest. As shown in [12], adequately accounting for the
covariance between tomographic bins dramatically improves the signal-to-noise
of the $\kappa$ reconstruction at low redshift by allowing us to extract
information about low redshift structures contained in the high redshift
tomographic bins.

As we now demonstrate, one of the principal advantages of karmma is that the
posterior maps faithfully reproduce the statistical properties of the input
convergence field. Figure 6 shows the difference between the mean two-point
function in our posterior maps and that of the input convergence map in both
real and harmonic space. For the latter, we compute pseudo-$C(l)$s using the
publicly available code namaster [50]. The blue line shows the mean
difference, averaged across all 108 simulations. The grey bands show the mean
error in the posterior for one simulation, while the blue bands show the error
on the mean. The difference in the two point functions between the input
simulated maps and our reconstructed maps is $\lesssim 5\%$, with the
difference increasing with decreasing redshift. These biases arise because of
the failure of the lognormal model to accurately model non-linear structure
formation. In particular, Appendix A demonstrates that we do not observe any
such biases when we test karmma using log-normal convergence maps as input.

Apart from the 2-point function, the lognormal model also captures many non-
Gaussian features of the convergence field. Figure 7 compares three such
statistics, namely the PDF of $\kappa$ values in each pixel (i.e. the one-
point function), and the peak and void counts of the $\kappa$ maps. The peaks
(voids) of the map are defined as the number of pixels which have the highest
(lowest) $\kappa$ value among all its neighbors. These statistics are well-
known to be sensitive to the non-Gaussian information of the convergence field
[51, 52]. We see that the non-Gaussian statistics of the karmma maps are close
to those in the simulations, particularly for the high redshift bins. However,
the differences between the simulated maps and our posteriors increases with
decreasing redshift, being particularly large for the lowest redshift bin.
This is not surprising, as it is at this redshift that non-linear structure
formation is most important. We note that because the signal-to-noise of the
convergence in an individual pixel is low ($S/N\lesssim 1$), the one point of
the posteriors is nearly identical to that of the prior.

##  V Bayesian mass maps with DES-Y3 data

Figure 8: The mean (top) and standard deviation (bottom) of the karmma
posterior maps for DES-Y3 data in each of the 4 tomographic bins, as labeled.

We run karmma on pixelized shear maps produced from the DES-Y3 shape catalogs
as detailed in section II.2. The input $C(l)$ for this run is computed for the
T17 cosmological parameters using pyccl [53]. We use the shift parameters
computed from the T17 simulations for these karmma runs. We produce a chain
with $500$ independent samples. The maps showing the mean and standard
deviation of the posteriors in each of the four tomographic bins is shown in
Figure 8.The figure highlights the fact that karmma produces full posterior
distributions of the mass maps, which is particularly important if these maps
are to be used for cross-correlation studies [e.g. 54].

Figure 9 compares the mean of our posterior maps for the lowest tomographic
redshift bin to the publicly available mass maps released by the DES
collaboration [55]. Specifically, we compare our mean mass map to their Null
B-mode prior Kaiser-Squires method, Wiener filter method, and the wavelet
transform based glimpse method [56]. We find that karmma can resolve smaller
structures than any of the other methods. This improvement is driven by the
fact that karmma self-consistently accounts for the covariance across all
tomographic redshift bins, which, as noted earlier, allows us to extract
information on the low-redshift density field from the high-redshift shear
data [12].

In Figure 10, we compare the correlation function and the pseudo $C(l)$ of the
DES-Y3 mass maps reconstructed using karmma and the same publicly available
mass maps used in Figure 9. In that figure, we also show the correlation
functions and the pseudo-$C(l)$s measured by the DES collaboration with black
dots.444DES-Y3 data vectors are available at
https://des.ncsa.illinois.edu/releases/y3a2/Y3key-products and the
pseudo-$C(l)$ measurements are taken from [45]. As we can see from the figure,
barring karmma , none of the other methods recover the expected 2-point
functions in the mass maps. This is unsurprising. Indeed, the fact that best
point-estimates of the convergence map are expected to be biased further
highlights the importance of being able to properly sample the posterior
distribution of the maps.

We run 2 additional karmma analysis with different priors to test the
sensitivity of the recovered mass maps to the input cosmological parameters.
We use the best fit cosmological parameters from Planck Collaboration
cosmological analysis [57] and DES-Y3 3$\times$2 pt analysis [1]. For each of
these cosmologies, we recompute the power spectrum and shift parameters. The
shift parameter at different cosmologies is computed by rescaling the shift
parameter in the T17 simulations with the ratio of the shift parameter
predicted using cosmomentum [58]. That is, we only rely on analytic
predictions for the shift parameter for the purposes of rescaling the shift
parameter from our fiducial cosmology to a new cosmology. The results of these
runs are shown in Figure 11, where we plot the fractional difference of the
recovered correlation function and pseudo-$C(\ell)$’s for the karmma runs with
different cosmological parameters with respect to the T17 2-point functions.
As we can see from the figure, the power spectrum on large scales is largely
insensitive to the input prior. On large scales, the signal-to noise is high
and therefore the large-scale modes are determined by the likelihood. By
contrast, the signal-to-noise ratio on small scales is low, which renders the
posterior at small scales sensitive to the adopted prior.

In Figure 12 we compare the distribution of $\chi^{2}$ values computed using
equation (9) for each map in the karmma posterior. Interestingly, we see that
the karmma maps produced using a Planck cosmology provide a better fit for the
DES-Y3 data than the DES-Y3 or T17 cosmological parameters. However, we
caution against interpreting this result as evidence for a Planck cosmology
since we did not account for observational systematics in our analysis.
Rather, we want to highlight the difference in the likelihood between the
Planck and DES-Y3 posteriors suggests that a field-level analysis of the
DES-Y3 data may have the potential to distinguish between these different
cosmologies.

Figure 9: Comparison of the mass map for the lowest tomographic redshift bin
in the DES-Y3 data. Each panel corresponds to a different algorithm: karmma
mean map (top left), Wiener filtering (top right), Kaiser-Squires with null
B-mode prior (bottom left), and glimpse (bottom right). The fact that karmma
adequately models the covariance between tomographic redshift bins enable us
to robustly recover structure on scales significantly smaller than those
resolved in other algorithms. Figure 10: Comparison of the correlation
function (bottom left) and the pseudo-$C(\ell)$ (top right) of mass maps
created from DES-Y3 data with the same data vectors computed directly from the
DES-Y3 shape catalogue (black). The blue shaded region shows the $95\%$
confidence interval of the summary statistics of karmma posteriors. The
correlation function and pseudo-$C(\ell)$’s of DES-Y3 mass maps created with
different methods, namely Null-B prior KS (red dash dotted), Wiener filter
(green dashed), glimpse (cyan dotted), are also shown in the figure. As we can
see from the figure, none of the mass mapping methods except karmma reproduce
the correlation functions of the DES-Y3 weak lensing data.  Figure 11: The
fractional difference in the recovered correlation function and pseudo
$C(\ell)$’s for karmma runs with different cosmological parameters with
respect to the theory predictions from T17 cosmology. The grey shaded region
shows the $95\%$ confidence interval for the karmma run with the fiducial
cosmological parameters of T17. The blue (red) solid line shows the mean of
the karmma posterior with Planck (DES-Y3) cosmological parameters. The blue
(red) dashed lines in the $C(\ell)$ plot shows the input $C(\ell)$ prior for
the corresponding cosmological parameters. As we can see from the figure, the
$C(\ell)$’s on large scales are not substantially impacted by the input
cosmological parameters. However, the small-scale power depends on the input
cosmological parameters. Figure 12: $\chi^{2}$ distribution from the karmma
posteriors with different cosmological parameters. The posteriors with DES-Y3
cosmological parameters (red vertical lines) have the highest $\chi^{2}$,
followed by the T17 (gray histograms) and P19 (blue diagonal lines). The fact
that these distributions are clearly differentiated suggests a field-level
cosmological analysis of the DES shear maps may be able to distinguish between
the Planck and DES-Y3 cosmologies.

##  VI Conclusion

The karmma algorithm [21, 31] models the convergence field as a lognormal
random field. By doing so, karmma can perform a forward-modelled Bayesian
reconstruction of the density field of the Universe as constrained by cosmic
shear data. In this paper, we have extended the karmma framework so that it
can simultaneously forward model multiple tomographic redshift bins while
fully accounting for their expected covariance. This enhancement significantly
improves the signal-to-noise of the reconstructed maps, particuarly in the low
signal-to-noise regime (See section IV and [12]). We validated our method on
simulated cosmic shear data generated using $N$-body simulations, and
demonstrated that the karmma posteriors accurately reproduce a variety of
Gaussian and non-Gaussian statistics of the input simulated maps (see Figures
6 and 7). However, we do find evidence of small $\lesssim 5\%-10\%$ biases
which increase with decreasing redshifts. These biases arise due to the
failure of the lognormal model to correctly capture non-linear structure
formation. These results suggest that further improvements require improving
the field-level prior assumed for the convergence field. Some options for such
an improvement include analytic extensions such as the double log model of [8]
or using generative artificial intelligence models [59, 60, 61, 62].

Following our simulation tests, we applied karmma to DES-Y3 weak lensing data
to produce the first Bayesian forward-modelled tomographic mass maps from
Stage-III weak lensing data on a sphere. These maps are also the first to
self-consistently account for the covariance between tomographic redshift bins
as part of the mass map reconstruction algorithm, which in turn improved the
signal-to-noise of the resulting maps, particularly at low redshifts. We show
that these maps have the correct theoretically expected statistical properties
such as the $2$-point function. We make these maps publicly available at
https://zenodo.org/records/10672062.

An important limitation in our work is the relatively coarse angular
resolution of our lensing maps ($\approx 13\ {\rm arcmin}$). Since the
information on small-scales is highly non-Gaussian, one might expect that
field-level inference will become even more useful at higher resolution.
However, this leads to increased computational demands due to: i) spherical
harmonic transforms at a higher resolution; and ii) a larger parameter space.
Significantly improving resolution over that achieved in this work will
require improved computational spherical harmonics code and faster sampling
methods. Improved spherical harmonic transform codes have been produced using
GPU acceleration [e.g, 63] or spherical Fourier neural operators [64].
Improvement in sampling can be achieved using improved sampling methods such
as microcanonical Hamiltonian Monte-Carlo methods [65, 66].

Despite this limitation, the inclusion of tomographic mass mapping on a sphere
represents an important step towards enabling a full field-level cosmological
analysis of cosmic shear data. In addition to varying cosmology, our future
work will focus on incorporating weak lensing systematics as part of the
inference process. The field-level approach is particularly well-suited for
incorporating some of these systematics. For example, it has been shown that
the non-Gaussian clustering of galaxies can improve the inference of
photometric redshifts [67, 68]. Likewise, Bayesian photo-$z$ inference
algorithms [69, 70, 71, 72] can naturally be extended to incorporate density
field inference to further improve photometric redshift estimation. Similarly,
we will need to model intrinsic alignments at the field level in order to
properly extract information using our forward modeling framework [73].
Fortunately, we see no reason why any of these challenges should be
prohibitive, suggesting that field-based inference of cosmological parameters
from cosmic shear data could soon be realized.

###### Acknowledgements.

We thank Marco Gatti for providing the DES-Y3 shape catalogue, and for useful
comment on the manuscript that helped improved the clarity of our
presentation. The computation presented here was performed on the High
Performance Computing (HPC) resources supported by the University of Arizona
TRIF, UITS, and Research, Innovation, and Impact (RII) and maintained by the
UArizona Research Technologies department. SSB was supported by the Department
of Energy Cosmic Frontier program, grant DE-SC0020215, and NSF grant 2009401.
ER’s work is supported is supported by NSF grant 2009401. ER also receives
funding from DOE grant DE-SC0009913 and NSF grant 2206688.

##  Appendix A Tests with lognormal mocks

Figure 13: Fractional error in the correlation function (left),
pseudo-$C(\ell)$ (centre) and the 1-pt PDF (right) of the first redshift bin
of the karmma posteriors for runs on lognormal mocks. As we can see, we do not
see the large biases in these summary statistics that we saw for the karmma
runs on $N$-body simulations in section IV. Thus, it shows that the biases
seen in section IV and Figures 6, 7 arises due to the model mis-specification
due to the lognormal prior.

In Figure 13, we show the fractional error in the correlation function,
pseudo-$C(\ell)$’s and the 1-pt PDF in the karmma posterior maps on runs with
lognormal mocks. For the ease of visualization, we only show the results for
the first redshift bins, where we saw in section IV that the uncertainty was
the largest. As we can see from the figure, we do not see the
$\sim\mathcal{O}(5\%)$ bias in the 2-pt functions or the even larger bias in
the 1-pt PDF. This demonstrates that the biases we see in our runs with
$N$-body simulations arise because of model mis-specification due to the
assumed lognormal prior.

## References

  * DES Collaboration [2022] DES Collaboration, Phys. Rev. D 105, 023520 (2022), arXiv:2105.13549 [astro-ph.CO] .
  * eBOSS Collaboration [2021] eBOSS Collaboration, Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 [astro-ph.CO] .
  * Millea  _et al._ [2019] M. Millea, E. Anderes, and B. D. Wandelt, Phys. Rev. D 100, 023509 (2019), arXiv:1708.06753 [astro-ph.CO] .
  * Millea  _et al._ [2020] M. Millea, E. Anderes, and B. D. Wandelt, Phys. Rev. D 102, 123542 (2020), arXiv:2002.00965 [astro-ph.CO] .
  * Millea  _et al._ [2021] M. Millea _et al._ , Astrophys. J. 922, 259 (2021), arXiv:2012.01709 [astro-ph.CO] .
  * Zhou and Dodelson [2023] A. J. Zhou and S. Dodelson, Phys. Rev. D 108, 083506 (2023), arXiv:2304.01387 [astro-ph.CO] .
  * Jasche and Lavaux [2019] J. Jasche and G. Lavaux, A&A 625, A64 (2019), arXiv:1806.11117 [astro-ph.CO] .
  * Porth  _et al._ [2023] L. Porth, G. M. Bernstein, R. E. Smith, and A. J. Lee, MNRAS 518, 3344 (2023), arXiv:2111.13702 [astro-ph.CO] .
  * Boruah  _et al._ [2022a] S. S. Boruah, G. Lavaux, and M. J. Hudson, MNRAS 517, 4529 (2022a), arXiv:2111.15535 [astro-ph.CO] .
  * Prideaux-Ghee  _et al._ [2023] J. Prideaux-Ghee, F. Leclercq, G. Lavaux, A. Heavens, and J. Jasche, MNRAS 518, 4191 (2023), arXiv:2204.00023 [astro-ph.CO] .
  * Bayer  _et al._ [2023a] A. E. Bayer, C. Modi, and S. Ferraro, JCAP 2023, 046 (2023a), arXiv:2210.15649 [astro-ph.CO] .
  * Boruah  _et al._ [2022b] S. S. Boruah, E. Rozo, and P. Fiedorowicz, MNRAS 516, 4111 (2022b).
  * Porqueres  _et al._ [2022] N. Porqueres, A. Heavens, D. Mortlock, and G. Lavaux, MNRAS 509, 3194 (2022), arXiv:2108.04825 [astro-ph.CO] .
  * Junzhe Zhou  _et al._ [2023] A. Junzhe Zhou, X. Li, S. Dodelson, and R. Mandelbaum, arXiv e-prints , arXiv:2312.08934 (2023), arXiv:2312.08934 [astro-ph.CO] .
  * Euclid Collaboration [2023] Euclid Collaboration, arXiv e-prints , arXiv:2301.12890 (2023), arXiv:2301.12890 [astro-ph.CO] .
  * Kovács  _et al._ [2022] A. Kovács _et al._ , MNRAS 510, 216 (2022), arXiv:2112.07699 [astro-ph.CO] .
  * Liu and Hill [2015] J. Liu and J. C. Hill, Phys. Rev. D 92, 063517 (2015), arXiv:1504.05598 [astro-ph.CO] .
  * Hojjati  _et al._ [2017] A. Hojjati _et al._ , MNRAS 471, 1565 (2017), arXiv:1608.07581 [astro-ph.CO] .
  * Kaiser and Squires [1993] N. Kaiser and G. Squires, Astrophys. J. 404, 441 (1993).
  * Mawdsley and DES Collaboration [2020] B. Mawdsley and DES Collaboration, MNRAS 493, 5662 (2020), arXiv:1905.12682 [astro-ph.CO] .
  * Fiedorowicz  _et al._ [2022a] P. Fiedorowicz, E. Rozo, S. S. Boruah, C. Chang, and M. Gatti, MNRAS 512, 73 (2022a), arXiv:2105.14699 [astro-ph.CO] .
  * Porqueres  _et al._ [2021] N. Porqueres, A. Heavens, D. Mortlock, and G. Lavaux, MNRAS 502, 3035 (2021), arXiv:2011.07722 [astro-ph.CO] .
  * Jasche and Wandelt [2013] J. Jasche and B. D. Wandelt, MNRAS 432, 894 (2013), arXiv:1203.3639 [astro-ph.CO] .
  * Alsing  _et al._ [2016] J. Alsing, A. Heavens, A. H. Jaffe, A. Kiessling, B. Wandelt, and T. Hoffmann, MNRAS 455, 4452 (2016), arXiv:1505.07840 [astro-ph.CO] .
  * Alsing  _et al._ [2017] J. Alsing, A. Heavens, and A. H. Jaffe, MNRAS 466, 3272 (2017), arXiv:1607.00008 [astro-ph.CO] .
  * Loureiro  _et al._ [2023] A. Loureiro, L. Whiteaway, E. Sellentin, J. S. Lafaurie, A. H. Jaffe, and A. F. Heavens, The Open Journal of Astrophysics 6, 6 (2023), arXiv:2210.13260 [astro-ph.CO] .
  * Sellentin  _et al._ [2023] E. Sellentin, A. Loureiro, L. Whiteway, J. S. Lafaurie, S. T. Balan, M. Olamaie, A. H. Jaffe, and A. F. Heavens, arXiv e-prints , arXiv:2305.16134 (2023), arXiv:2305.16134 [astro-ph.CO] .
  * Hoffman and Ribak [1991] Y. Hoffman and E. Ribak, ApJL 380, L5 (1991).
  * Zaroubi  _et al._ [1995] S. Zaroubi, Y. Hoffman, K. B. Fisher, and O. Lahav, Astrophys. J. 449, 446 (1995), arXiv:astro-ph/9410080 [astro-ph] .
  * Kodi Ramanah  _et al._ [2019] D. Kodi Ramanah, G. Lavaux, and B. D. Wandelt, MNRAS 490, 947 (2019), arXiv:1906.10704 [astro-ph.CO] .
  * Fiedorowicz  _et al._ [2022b] P. Fiedorowicz, E. Rozo, and S. S. Boruah, arXiv e-prints , arXiv:2210.12280 (2022b), arXiv:2210.12280 [astro-ph.CO] .
  * Xavier  _et al._ [2016] H. S. Xavier, F. B. Abdalla, and B. Joachimi, MNRAS 459, 3693 (2016), arXiv:1602.08503 [astro-ph.CO] .
  * Clerkin  _et al._ [2017] L. Clerkin _et al._ , MNRAS 466, 1444 (2017), arXiv:1605.02036 [astro-ph.CO] .
  * Takahashi  _et al._ [2017] R. Takahashi, T. Hamana, M. Shirasaki, T. Namikawa, T. Nishimichi, K. Osato, and K. Shiroyama, Astrophys. J. 850, 24 (2017), arXiv:1706.01472 [astro-ph.CO] .
  * Amon and DES Collaboration [2022] A. Amon and DES Collaboration, Phys. Rev. D 105, 023514 (2022), arXiv:2105.13543 [astro-ph.CO] .
  * Secco and DES Collaboration [2022] L. F. Secco and DES Collaboration, Phys. Rev. D 105, 023515 (2022), arXiv:2105.13544 [astro-ph.CO] .
  * Myles  _et al._ [2021] J. Myles _et al._ , MNRAS 505, 4249 (2021), arXiv:2012.08566 [astro-ph.CO] .
  * Gatti  _et al._ [2021] M. Gatti _et al._ , MNRAS 504, 4312 (2021), arXiv:2011.03408 [astro-ph.CO] .
  * Lewis  _et al._ [2000] A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), arXiv:astro-ph/9911177 [astro-ph] .
  * Lesgourgues [2011] J. Lesgourgues, arXiv e-prints , arXiv:1104.2934 (2011), arXiv:1104.2934 [astro-ph.CO] .
  * Gatti  _et al._ [2024] M. Gatti, , and DES Collaboration, MNRAS 527, L115 (2024), arXiv:2307.13860 [astro-ph.CO] .
  * Huff and Mandelbaum [2017] E. Huff and R. Mandelbaum, arXiv e-prints , arXiv:1702.02600 (2017), arXiv:1702.02600 [astro-ph.CO] .
  * Sheldon and Huff [2017] E. S. Sheldon and E. M. Huff, Astrophys. J. 841, 24 (2017), arXiv:1702.02601 [astro-ph.CO] .
  * Nicola  _et al._ [2021] A. Nicola, C. García-García, D. Alonso, J. Dunkley, P. G. Ferreira, A. Slosar, and D. N. Spergel, JCAP 2021, 067 (2021), arXiv:2010.09717 [astro-ph.CO] .
  * Doux and DES Collaboration [2022] C. Doux and DES Collaboration, MNRAS 515, 1942 (2022), arXiv:2203.07128 [astro-ph.CO] .
  * Wallis  _et al._ [2022] C. G. R. Wallis, M. A. Price, J. D. McEwen, T. D. Kitching, B. Leistedt, and A. Plouviez, MNRAS 509, 4480 (2022), arXiv:1703.09233 [astro-ph.CO] .
  * Paszke  _et al._ [2019] A. Paszke _et al._ , in  _Advances in Neural Information Processing Systems 32_ (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Bingham  _et al._ [2019] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. A. Szerlip, P. Horsfall, and N. D. Goodman, J. Mach. Learn. Res. 20, 28:1 (2019).
  * Hoffman and Gelman [2011] M. D. Hoffman and A. Gelman, arXiv e-prints , arXiv:1111.4246 (2011), arXiv:1111.4246 [stat.CO] .
  * Alonso  _et al._ [2019] D. Alonso, J. Sanchez, A. Slosar, and LSST Dark Energy Science Collaboration, MNRAS 484, 4127 (2019), arXiv:1809.09603 [astro-ph.CO] .
  * Liu  _et al._ [2015] J. Liu, A. Petri, Z. Haiman, L. Hui, J. M. Kratochvil, and M. May, Phys. Rev. D 91, 063507 (2015), arXiv:1412.0757 [astro-ph.CO] .
  * Davies  _et al._ [2021] C. T. Davies, M. Cautun, B. Giblin, B. Li, J. Harnois-Déraps, and Y.-C. Cai, MNRAS 507, 2267 (2021), arXiv:2010.11954 [astro-ph.CO] .
  * Chisari and LSST Dark Energy Science Collaboration [2019] N. E. Chisari and LSST Dark Energy Science Collaboration, ApJS 242, 2 (2019), arXiv:1812.05995 [astro-ph.CO] .
  * Nguyen  _et al._ [2020] N.-M. Nguyen, J. Jasche, G. Lavaux, and F. Schmidt, JCAP 2020, 011 (2020), arXiv:2007.13721 [astro-ph.CO] .
  * Jeffrey  _et al._ [2021] N. Jeffrey, M. Gatti, and DES Collaboration, MNRAS 505, 4626 (2021), arXiv:2105.13539 [astro-ph.CO] .
  * Lanusse  _et al._ [2016] F. Lanusse, J. L. Starck, A. Leonard, and S. Pires, A&A 591, A2 (2016), arXiv:1603.01599 [astro-ph.CO] .
  * Planck Collaboration [2020] Planck Collaboration, A&A 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Friedrich  _et al._ [2020] O. Friedrich, C. Uhlemann, F. Villaescusa-Navarro, T. Baldauf, M. Manera, and T. Nishimichi, MNRAS 498, 464 (2020), arXiv:1912.06621 [astro-ph.CO] .
  * Fiedorowicz  _et al._ [2022c] P. Fiedorowicz, E. Rozo, S. Boruah, W. Coulton, S. Ho, and G. Fabbian, in  _Machine Learning for Astrophysics_ (2022) p. 8.
  * Yiu  _et al._ [2022] T. W. H. Yiu, J. Fluri, and T. Kacprzak, JCAP 2022, 013 (2022), arXiv:2112.12741 [astro-ph.CO] .
  * Dai and Seljak [2023] B. Dai and U. Seljak, arXiv e-prints , arXiv:2306.04689 (2023), arXiv:2306.04689 [astro-ph.CO] .
  * Shirasaki and Ikeda [2023] M. Shirasaki and S. Ikeda, arXiv e-prints , arXiv:2310.17141 (2023), arXiv:2310.17141 [astro-ph.CO] .
  * Tian  _et al._ [2022] C. Tian, S. Li, and H. Liu, ApJS 263, 21 (2022), arXiv:2208.10154 [astro-ph.IM] .
  * Bonev  _et al._ [2023] B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar, arXiv e-prints , arXiv:2306.03838 (2023), arXiv:2306.03838 [cs.LG] .
  * Robnik  _et al._ [2022] J. Robnik, G. B. De Luca, E. Silverstein, and U. Seljak, arXiv e-prints , arXiv:2212.08549 (2022), arXiv:2212.08549 [stat.CO] .
  * Bayer  _et al._ [2023b] A. E. Bayer, U. Seljak, and C. Modi, arXiv e-prints , arXiv:2307.09504 (2023b), arXiv:2307.09504 [astro-ph.CO] .
  * Jasche and Wandelt [2012] J. Jasche and B. D. Wandelt, MNRAS 425, 1042 (2012), arXiv:1106.2757 [astro-ph.CO] .
  * Tsaprazi  _et al._ [2023] E. Tsaprazi, J. Jasche, G. Lavaux, and F. Leclercq, arXiv e-prints , arXiv:2301.03581 (2023), arXiv:2301.03581 [astro-ph.CO] .
  * Leistedt  _et al._ [2016] B. Leistedt, D. J. Mortlock, and H. V. Peiris, MNRAS 460, 4258 (2016), arXiv:1602.05960 [astro-ph.CO] .
  * Leistedt  _et al._ [2023] B. Leistedt, J. Alsing, H. Peiris, D. Mortlock, and J. Leja, ApJS 264, 23 (2023), arXiv:2207.07673 [astro-ph.IM] .
  * Sánchez and Bernstein [2019] C. Sánchez and G. M. Bernstein, MNRAS 483, 2801 (2019), arXiv:1807.11873 [astro-ph.CO] .
  * Alarcon  _et al._ [2020] A. Alarcon, C. Sánchez, G. M. Bernstein, and E. Gaztañaga, MNRAS 498, 2614 (2020), arXiv:1910.07127 [astro-ph.CO] .
  * Tsaprazi  _et al._ [2022] E. Tsaprazi, N.-M. Nguyen, J. Jasche, F. Schmidt, and G. Lavaux, JCAP 2022, 003 (2022), arXiv:2112.04484 [astro-ph.CO] .