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

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). 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. 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. 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. 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.

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.