Comparing Mass Mapping Reconstruction Methods with Minkowski Functionals

Abstract

Using higher-order statistics to capture cosmological information from weak lensing surveys often requires a transformation of observed shear to a measurement of the convergence signal. This inverse problem is complicated by noise and boundary effects, and various reconstruction methods have been developed to implement the process. Here we evaluate the retention of signal information of four such methods: Kaiser-Squires, Wiener filter, DarkMappy, and DeepMass. We use the higher order statistics Minkowski functionals to determine which method best reconstructs the original convergence with efficiency and precision. We find DeepMass produces the tightest constraints on cosmological parameters, while Kaiser-Squires, Wiener filter, and DarkMappy are similar at a smoothing scale of 3.5 arcmin. We also study the MF inaccuracy caused by inappropriate training sets in the DeepMass method and find it to be large compared to the errors, underlining the importance of selecting appropriate training cosmologies.

1 Introduction

Weak gravitational lensing uses the distortion of the shapes and/or size of galaxies as a probe of the matter distribution along the line of sight. Maps quantifying lensing can be generated by measuring the observed shapes of galaxies, utilizing this shear effect from gravitational lensing (Hikage et al., 2019; Hamana et al., 2020; Asgari et al., 2021; Amon et al., 2022; Secco et al., 2022). Statistics measured from weak gravitational lensing are effective probes of the evolution of large-scale structure. Galaxy surveys, such as the Dark Energy Survey (DES; Flaugher, 2005), Hyper Supreme-Cam (HSC; Aihara et al., 2017) survey, and Kilo-Degree Survey (KiDS; Kuijken et al., 2019), constructed detailed galaxy catalogues of tens to hundreds of millions of observed galaxies (Gatti et al., 2021a; Giblin et al., 2021; Aihara et al., 2022). The Legacy Survey Space and Time at the Rubin Observatory (LSST; LSST Science Collaboration et al., 2009; LSST Dark Energy Science Collaboration, 2012) and the space telescopes Euclid (Scaramella et al., 2021) and Nancy Grace Roman (Spergel et al., 2015) will go further in depth and area, and so improve our understanding of dark matter and dark energy.

Weak lensing causes source galaxies shapes to be magnified by the convergence field and distorted by the shear field. While only the shear signal can be directly inferred from galaxy shape observations, convergence can conveniently be represented as an integral of the projected matter density along the line of sight and can be derived from shear. Without noise or boundary effects, the observed shear field can be used to measure the convergence field trivially, but this is not the case for real data. We need to solve an inverse problem, working backwards from observational shear measurements to infer the underlying convergence in the presence of these effects in real data, which becomes complex.

Various non-parametric mass mapping reconstruction methods have been developed to address this problem and take into account noise and masks in the data to solve the inverse problem. Linear approaches like Kaiser Squires (KS) (Kaiser & Squires, 1993), Wiener filter (WF) (Lahav et al., 1994), Simon et al. (2009), and VanderPlas et al. (2011) and non-linear approaches like GLIMPSE (Leonard et al., 2014), Glimpse2D (Lanusse et al., 2016), KS+ (Pires et al., 2020), DeepMass (Jeffrey et al., 2020), MCAlens (Starck et al., 2021), DarkMappy (Price et al., 2021), KaRMMa (Fiedorowicz et al., 2022), DeepPosterior (Remy et al., 2022), SKS+ (Kansal, 2023), and others (Marshall et al., 2002; Simon et al., 2011b; Leonard et al., 2012; Alsing et al., 2015; Schneider et al., 2017; Horowitz et al., 2019; Porqueres et al., 2021) have all been shown to be effective. The original method, Kaiser-Squires (KS) is a linear inversion of the shear map into a convergence map, but it cannot take noise or boundary effects into account (Kaiser & Squires, 1993). Wiener filtering is also a linear filter, but unlike KS, it can transform noisy data well (Lahav et al., 1994). Sparsity methods like GLIPMSE and DarkMappy are non-linear and able to reconstruct information in the non-Gaussian regime (Leonard et al., 2014; Price et al., 2021). Some methods like DarkMappy and DeepPosterior are able to approximate the uncertainty that comes from reconstruction. Other methods including DeepPosterior have employed machine learning as a useful tool in denoising convergence signal (Jeffrey et al., 2020; Shirasaki et al., 2021; Remy et al., 2022). In this paper we explore four reconstruction methods: Kaiser Squires, Wiener filter, DarkMappy, and DeepMass.

Measuring statistics from the best reconstructed convergence maps is expected to lead to tighter constraints on cosmological parameters. Two-point statistics, which describe the correlation between pairs of data points, can constrain these parameters with high precision (Heymans et al., 2021; Miyatake et al., 2021; Abbott et al., 2022). However, statistics beyond two-point are required to fully probe non-Gaussian structure on small scales. Higher order statistics including N-point such as the one-point probability distribution function (Liu & Madhavacheril, 2019; Barthelemy et al., 2020; Thiele et al., 2020; Boyle et al., 2021) and three-point correlation functions (Schneider & Lombardi, 2003; Takada & Jain, 2004; Gong et al., 2023), moments (Van Waerbeke et al., 2013; Petri et al., 2015; Vicinanza et al., 2016; Chang et al., 2018; Peel et al., 2018; Vicinanza et al., 2018; Gatti et al., 2020, 2021b; Porth & Smith, 2021), topological descriptors such as Betti numbers (Feldbrugge et al., 2019; Heydenreich et al., 2021; Parroni et al., 2021) and peak counts (Marian et al., 2009; Kratochvil et al., 2010; Liu et al., 2015; Kacprzak et al., 2016; Shan et al., 2017; Martinet et al., 2018; Peel et al., 2018; Ajani et al., 2020; Harnois-Déraps et al., 2021; Zürcher et al., 2022), scattering transform coefficients (Cheng et al., 2020; Cheng & Ménard, 2021; Gatti et al., 2023), and field level inference (Porqueres et al., 2021; Boruah et al., 2022) can measure information in the non-Gaussian regime in weak lensing maps.

Another higher-order statistic, Minkowski functionals (MFs) are morphological descriptors that describe the topology of a continuous field (Minkowski, 1903). MFs are unbiased, have low variance, and can have additional resilience to some systematic uncertainties, making them effective probes of the underlying dark matter distribution (Zürcher et al., 2021). MFs have been applied to Cosmic Microwave Background data (Schmalzing & Gorski, 1998; Eriksen et al., 2004; Hikage et al., 2006; Buchert et al., 2017; Collaboration et al., 2020; Hamann & Kang, 2023; Chingangbam & Rahman, 2023), 2D convergence maps (Mecke et al., 1993; Kratochvil et al., 2012; Petri et al., 2013; Vicinanza et al., 2019; Parroni et al., 2020; Grewal et al., 2022; Euclid Collaboration et al., 2023), 2D density fields (Grewal et al., 2022), and 3D density fields from spectroscopic data (Hikage et al., 2003; Wiegand & Eisenstein, 2017; Sullivan et al., 2019; Appleby et al., 2022; Liu et al., 2022). Like other higher-order statistics, MFs produce tighter constraints on cosmological parameters because they are sensitive to small-scale structure.

In this paper we use MFs as a test of different convergence reconstruction methods, using them to quantify the information (particularly on small scales) that the methods can reconstruct in maps. We compare the four different mass mapping methods KS, WF, DarkMappy, and DeepMass, applying them to simulated images and using MFs to constrain $w$CDM parameters and thus computing figures of merit to describe method power.

In Section 2 we review weak lensing formalism and the mass mapping reconstruction methods, in Section 3 we describe the simulations, observables, and evaluation metrics, in Section 4 we compare the cosmological constraints from the four reconstruction methods, and we conclude in Section 5.

2 Formalism

2.1 Convergence and Shear

The magnification and distortion of a lensed galaxy image is represented to first order by the Jacobian matrix

$$\mathcal{A}(\boldsymbol{\theta})=\left(\delta_{ij}-\frac{\partial^{2}\psi(\boldsymbol{\theta})}{\partial\theta_{i}\partial\theta_{j}}\right)=\begin{pmatrix}1-\kappa-\gamma_{1}&-\gamma_{2}\\ -\gamma_{2}&1-\kappa+\gamma_{1}\end{pmatrix},$$ (1)

where $\theta$ is the observed angular position, $\psi$ is the lensing potential, the shear $\gamma\equiv\gamma_{1}+i\gamma_{2}$, and convergence is $\kappa$ (Bartelmann & Schneider, 2001). Convergence is a measure of the matter distribution along the line of sight. It is given by a projection of the overdensity $\delta$:

$$\kappa(\mathbf{\hat{\theta}})=\frac{3\mathrm{H}_{0}^{2}\Omega_{\mathrm{m}}}{2c^{2}}\int_{0}^{\chi_{\mathrm{lim}}}d\chi\frac{q(\chi)}{a(\chi)}f_{K}(\chi)\delta(f_{K}(\chi)\mathbf{\hat{\theta}},\chi),$$ (2)

where $\mathrm{H}_{0}$ is the Hubble constant, $\Omega_{\mathrm{m}}$ is the matter density, $c$ is the speed of light, $\chi$ is comoving distance, $\chi_{\mathrm{lim}}$ is the horizon distance, $a$ is the scale factor, and the lens efficiency $q$ is defined as

$$q(\chi)=\int_{0}^{\chi_{\mathrm{lim}}}n(\chi^{\prime})d\chi^{\prime}\frac{f_{K}(\chi^{\prime}-\chi)}{f_{K}(\chi^{\prime})}.$$ (3)

$n(\chi)\mathrm{d}\chi$ is the source galaxy distribution and $f_{K}$ is the comoving angular distance (Kilbinger, 2015).

While convergence is harder to observe directly, shear can be measured from galaxy shape observations. The convergence $\kappa$ is related to the shear $\gamma$ via the lensing potential $\psi$. Rearranging the following equations for $\kappa$ and $\gamma$ and solving for $\psi$,

$$\kappa=\frac{1}{2}\Delta\psi,$$ (4)

$$\gamma_{1}=\frac{1}{2}\left(\partial_{1}^{2}\psi-\partial_{2}^{2}\psi\right);\gamma_{2}=\partial_{1}\partial_{2}\psi,$$ (5)

and then performing a Fourier transform in the flat sky regime, the convergence can be expressed as

$$\tilde{\kappa}=\frac{k_{1}^{2}-k_{2}^{2}}{k^{2}}\tilde{\gamma_{1}}+\frac{2k_{1}k_{2}}{k^{2}}\tilde{\gamma_{2}}.$$ (6)

$k_{1}$ and $k_{2}$ are spatial frequencies in Fourier space ($k^{2}=k_{1}^{2}+k_{2}^{2}$), and $\tilde{\gamma_{1}}$ and $\tilde{\gamma_{2}}$ are the two orthogonal components of shear distortion (Remy et al., 2022). As the expectation value of matter overdensity $\delta$ is zero, the first moment of shear and convergence will also be zero. Known as mass sheet degeneracy, this feature means that at $k$ = 0, the relationship between convergence and shear is undefined. It is important to note that this transformation is only optimal when performed on pure signal in the absence of noise and missing data.

In the presence of noise and gaps from masking, we must go beyond this simple treatment. In this more realistic regime, various algorithms attempt the inversion from shear to convergence in different ways. We review several of these methods here.

2.2 Kaiser-Squires

Kaiser-Squires uses a direct linear inversion of the shear in Fourier space to obtain the reconstructed convergence field (Kaiser & Squires, 1993). Scalar perturbations in matter density are captured by the E mode gradient-like patterns, and the tensor perturbation counterparts are captured by the curl-like B mode. The shear-convergence relation can be expanded to give both the E and B modes:

$$\tilde{\kappa}=\tilde{\kappa_{E}}+i\tilde{\kappa_{B}}=\left(\frac{k_{1}^{2}-k_{2}^{2}}{k^{2}}+i\frac{2k_{1}k_{2}}{k^{2}}\right)(\tilde{\gamma_{1}}+i\tilde{\gamma_{2}}).$$ (7)

We can define the forward model in Fourier space $\mathbf{P}$ as

$$\textbf{{P}}=\left(\frac{k_{1}^{2}-k_{2}^{2}}{k^{2}}+i\frac{2k_{1}k_{2}}{k^{2}}\right)$$ (8)

(Kaiser & Squires, 1993). We note that we do not expect weak lensing to generate B modes, which come from systematics and higher-order lensing effects.

KS disregards noise and boundary effects, which means it requires continuous fields to reconstruct maps with masks (Van Waerbeke et al., 2013). In common with other approaches, we use a Gaussian kernel to fill in the gaps and smooth the field, though this leads to loss of structure on small scales and suppresses peaks in the convergence (Jeffrey et al., 2018). This also leads to leakage of the E and B modes of the reconstructed convergence field (Remy et al., 2022).

2.3 Wiener filter

A Wiener filter (WF) is a linear inversion of signal and noise together, aiming to suppress the latter. The WF method uses the expected power spectrum of the convergence field at a fiducial cosmology and of the noise (Simon et al., 2011a). In Fourier space:

$$\tilde{\kappa}=\textbf{SP}^{\dagger}\left[\textbf{{PSP}}^{\dagger}+\textbf{{N}}\right]^{-1}\tilde{\gamma},$$ (9)

where S is the covariance matrix of power spectrum signal and N is the noise covariance matrix (Lahav et al., 1994; Zaroubi et al., 1995). However, in filtering out the noise, WF smooths out small-scales features in the map more generally. While WF can reconstruct the field with precision on large scales, it does not typically retain small-scale information in convergence maps as well (Simon et al., 2011a; Horowitz et al., 2019).

2.4 DarkMappy

DarkMappy was primarily designed for cluster analysis; here we apply the method to a wide-field convergence. DarkMappy models the convergence field using wavelets (Price et al., 2021). It uses a sparsity approach to minimise the number of wavelets used and varies their coefficients in a hierarchical Bayesian approach, which enables an approximation of the signal uncertainty that comes from reconstruction. As it is a non-linear approach, we expect it to measure small-scale effects smoothed out by KS and WF.

A wavelet dictionary is a collection of wavelets, which are oscillatory functions with zero mean that are localised in Fourier and real space (Mallat, 2012). A wavelet dictionary $\mathbf{\Psi}$ can sparsely represent the signal $\kappa$:

$$\kappa=\mathbf{\Psi}\alpha,$$ (10)

where only a subset of non-zero coefficients $\alpha$ is needed (Leonard et al., 2014).

Unlike the previous methods, DarkMappy fits the data in real space. The forward model operation in Eq. (8) is transformed into a measurement operator $\mathbf{\Phi}$ using a Fourier transform $\mathbf{F}$:

$$\mathbf{\Phi}=\mathbf{F}^{-1}\mathbf{PF}.$$ (11)

$\mathbf{\Phi}$ maps convergence onto shear in the measurement equation, where the total shear is then given by:

$$\gamma=\mathbf{\Phi}\kappa+n,$$ (12)

where $n$ is the noise contaminating the signal. Using a wavelet-based, sparsity-promoting, $l_{1}$-norm prior, the DarkMappy model for the reconstructed map is then given by the maximum a posteriori solution:

$$\kappa^{\mathrm{map}}=\underset{\kappa}{\mathrm{argmin}}\biggl{\{}\mu\parallel\mathbf{\Psi}^{\dagger}\kappa\parallel_{1}+\frac{1}{2}\parallel\Sigma^{-\frac{1}{2}}(\mathbf{\Phi}\kappa-\gamma)\parallel_{2}^{2}\biggr{\}}.$$ (13)

The first term is the prior, which imposes sparsity, and the second term is the likelihood of the data. $\mu$ regularises the weighting between the two terms and $\Sigma$ is the covariance of the measured shear (Price et al., 2021).

The search process for the MAP uses knock-out hypothesis testing of the posterior to differentiate between original signal and reconstruction effects (Cai et al., 2018). This algorithm tests whether features (e.g. peaks) are significant enough to be included in the MAP model using a statistical test that compares $\kappa$ models that include them to those using $\kappa$ smoothed by surrounding pixels.

2.5 DeepMass

DeepMass is a convolutional neural network (CNN) that aims to denoise noisy maps. The network learns parameters $\Theta$ in convolution layers to reconstruct convergence maps from noisy shear maps:

$$\kappa=\mathcal{F}_{\Theta}(\gamma),$$ (14)

where $\mathcal{F}_{\Theta}(\gamma)$ is the posterior estimate of convergence (Jeffrey et al., 2020). Convergence and shear training data is used to minimise a mean-square-error (MSE) cost function $J$ as the network finds $\Theta$:

$$J(\Theta)=\parallel\mathcal{F}_{\Theta}(\gamma)-\kappa\parallel_{2}^{2}.$$ (15)

Conceptually, the ‘true’ output noise-free convergence maps $\boldsymbol{\kappa}$ are drawn from a prior distribution $\mathrm{P}(\kappa)$, and the corresponding input noisy shear maps are drawn from the likelihood $\mathrm{P}(\gamma|\kappa)$ (Jeffrey et al., 2020). The reconstructed convergence $\kappa$ can then be approximated as

$$\kappa=\mathcal{F}_{\Theta}(\boldsymbol{\gamma})=\int\kappa\mathrm{P}(\kappa|\gamma)\textrm{d}\kappa.$$ (16)

In practice, the network is trained on pairs of truth convergence maps and noisy WF convergence maps, which are reconstructed from shear maps (Jeffrey et al., 2020). The model is two-dimensional, so that multiple models need to be built for data with more than one redshift bin. DeepMass utilises a U-Net architecture with a contracting path (encoder), reducing the dimensions of the input data while extracting high-level features, and an expanding path (decoder), which complements the contracting path by restoring the full spatial dimension and generating a high-resolution output (Ronneberger et al., 2015). It also uses average pooling, where each region of the input data is replaced with its mean value. As the resolution decreases, the model captures more comprehensive physical features in the convergence map by considering a wider area and reducing the dimensionality of the data (Géron, 2019). While DeepMass can measure small-scale structure in convergence maps, the model does not retrieve uncertainties in the reconstruction (Jeffrey et al., 2020; Remy et al., 2022). As with all CNNs, DeepMass risks overfitting to map structure and/or creating fake artefacts in the reconstructed signal.

2.6 Minkowski functionals

Minkowski functionals are field integrals that characterise the topological properties of continuous fields (Minkowski, 1903; Zürcher et al., 2021). In this work we measure these mathematical descriptors from a convergence field. MFs measure the properties of excursion sets of this field, which are given by the region of the field above a given threshold (Parroni et al., 2020). The first three functionals quantify area, perimeter, and mean curvature of an excursion set and are defined as

$$V_{0}(t)=\frac{1}{A}\int\limits_{A}\Xi(\kappa(\mathbf{x})-t)\mathrm{d}^{2}\mathbf{x},$$ (17)

$$V_{1}(t)=\frac{1}{4A}\int\limits_{A}\delta(\kappa(\mathbf{x})-t)\sqrt{\kappa_{x}^{2}+\kappa_{y}^{2}}\mathrm{d}^{2}\mathbf{x},$$ (18)

$$V_{2}(t)=\frac{1}{2\pi A}\int\limits_{A}\delta(\kappa(\mathbf{x})-t)\\ \left(\frac{2\kappa_{x}\kappa_{y}\kappa_{xy}-\kappa_{x}^{2}\kappa_{yy}-\kappa_{y}^{2}\kappa_{xx}}{\kappa_{x}^{2}+\kappa_{y}^{2}}\right)\mathrm{d}^{2}\mathbf{x},$$ (19)

where $\mathbf{x}=(x,y)$ is the location in the field, $t$ is a chosen threshold, $A$ is the total area of the map, $x$ and $y$ are flat sky coordinates, $\kappa(\textbf{x})$ is the field value in two dimensions, and $\kappa_{x}$, $\kappa_{y}$, $\kappa_{xy}$, $\kappa_{yy}$, and $\kappa_{xx}$ are derivatives of the field (Minkowski, 1903; Petri et al., 2013). The derivatives are also evaluated at position $\mathbf{x}$; this is dropped in Eq. (18) and (19) for brevity.

The Heaviside function $\Xi$ in Eq. (17) identifies the area of the field region above the threshold; similarly, the Dirac delta functions in Eq. (18) and (19) select the regions where the height matches the threshold to measure the perimeter and connectivity of the field as the curvature of the boundary, respectively (Minkowski, 1903).

3 Methodology

3.1 Simulations

To generate our convergence maps we use CosmoGridV1, which comprises a suite of full-sky lightcone simulations with nside = 2048 and shells that allow for probing of multiple redshifts, up to z $<$ 3.5 (Kacprzak et al., 2022). The simulations were run with PKDGRAV3, a high performance self-gravitating astrophysical N-body treecode (Potter et al., 2016). CosmoGridV1 varies the cosmological parameters $\Omega_{\mathrm{m}}$, $\sigma_{8}$, $w_{0}$, $n_{\mathrm{s}}$, $\Omega_{\mathrm{b}}$, and $\mathrm{H}_{0}$ and baryonic parameters $M^{0}_{\mathrm{c}}$ and $\nu$, but uses fixed neutrino masses, each with $m_{v}$ = 0.02eV. To make realistic lensing maps, baryonic effects are included using a shell-based baryon correction model at the map level. Particle count maps are adjusted by a 2D displacement function measured from halo catalogues (Kacprzak et al., 2022). Following Schneider et al. (2019), the mass dependence of the gas profile $M_{\mathrm{c}}$ is broken down into the amplitude $M_{\mathrm{c}}^{0}$ and redshift dependence parameter $\nu$:

$$M_{\mathrm{c}}=M_{\mathrm{c}}^{0}(1+z)^{\nu}.$$ (20)

CosmoGridV1 provides 200 simulations at a fiducial cosmology with different seeds for the initial conditions, accompanied by maps with paired matching seeds with two steps, one above and one below for each of the cosmological parameters, shown in Table 1. For each cosmology, we cut out 50,000 15 degree $\times$ 15 degree patches with random positions and rotations from the 200 simulations, removing values at the poles to avoid distortion. The resulting convergence maps have resolution 256 $\times$ 256 pixels with pixel size 3.5 arcmin.

The CosmoGridV1 convergence maps have only positive values; here we subtract the mean after reconstruction before measuring observables. The mass-sheet degeneracy means the zero-point is arbitrary when going from shear back to convergence.

	Fiducial	$\Delta$ Fiducial
$\Omega_{\mathrm{m}}$	0.26	$\pm$ 0.01
$\sigma_{8}$	0.84	$\pm$ 0.015
$w_{0}$	-1	$\pm$ 0.05
$n_{\mathrm{s}}$	0.9649	$\pm$ 0.02
$\Omega_{\mathrm{b}}$	0.0493	$\pm$ 0.001
$\mathrm{H}_{0}$	67.3	$\pm$ 2.0
$M_{\mathrm{c}}^{0}$	13.82	$\pm$ 0.1
$\nu$	0.0	$\pm$ 0.1

Table 1: CosmoGridV1 Parameter Shift Values. $\mathrm{H}_{0}$ has units $\mathrm{kms^{-1}Mpc^{-1}}$.

In this work, we do our analysis with 5,000 map patches, a balance between efficiency and precision found by requiring convergence of the KS contours described below, where there is less than a 5% shift in contours. For DeepMass training and testing, we use 50,000 random patches from 200 maps to avoid overfitting at the given cosmology. However, all observables are only measured from 5,000 maps.

3.1.1 Redshift

In our analysis, we use approximately DES Y3-like tomographic redshifts from the CosmoGridV1 shells, modelled as a Smail-type distribution:

$$n(z)=z^{\alpha}\mathrm{exp}\left[-\left(\frac{z}{z_{0}}\right)^{\beta}\right],$$ (21)

(Smail et al., 1995; Kacprzak et al., 2022). The survey parameter values can be seen in Table 2 and the final distributions can be seen in Figure 1.

Bin	$\alpha$	$\beta$	$z_{0}$
1	1.99	1.44	0.20
2	3.46	2.34	0.39
3	6.03	3.60	0.66
4	3.53	4.49	1.03

Table 2: Smail-type distribution parameters for DES Y3-like redshift bins (Fischbacher et al., 2022)

Figure 1: The DES Y3-like $n_{z}$ redshift distributions from CosmoGridV1 used in this work (Kacprzak et al., 2022).

3.1.2 Noise and Mask

At our chosen resolution our maps are in the low noise regime. We only take shape noise into account, disregarding any other observational effects. To obtain the noise standard deviation per pixel, we use

$$\sigma=\frac{\sigma_{e}}{\sqrt{\frac{n_{\mathrm{gal}}}{n_{\mathrm{bin}}}p^{2}}},$$ (22)

where the shape dispersion $\sigma_{e}$ is 0.26 per galaxy, $n_{\mathrm{gal}}$ is a DES Y3-like number density (10 galaxies/arcmin${}^{2}$), $n_{\mathrm{bin}}$ is the four DES Y3-like redshift bins, and $p$ is the pixel size in arcmin. With 256${}^{2}$ total pixels and 15 $\times$ 15 square degree maps, the pixel size $p$ is 3.5 arcmin and the noise level $\sigma$ is 0.047 per pixel. The map noise is computed as an uncorrelated random Gaussian field with mean 0 and $\sigma$ from Eq. (22).

Figure 2: Mask applied to a 15 $\times$ 15 degree flat sky convergence map of redshift bin 0 at fiducial cosmology. The colourbar shows convergence of each pixel. We note that in practice the map is applied to shear map, only visualised here on convergence data.

We apply a mask to our maps by hand to approximate the structure of real masks, covering approximately 2.8% of the map, as seen in Figure 2.

The noise and mask are added at the shear step. A Gaussian kernel with a standard deviation of one pixel has to be applied to the masked shear map before reconstruction. This smoothing fills in gaps and decreases the resolution of the maps. This is a crucial step for Kaiser Squires, which only works for continuous fields.

3.1.3 Applications of Reconstruction Methods

For the Wiener filter reconstruction calculation of the signal matrix $\mathbf{S}$ in Eq. (9), we compute the 2D flat-sky power spectrum $P(k)$ with CCL (Chisari et al., 2019) using the cosmology of the fiducial convergence map and incorporating baryonic effects using HMCode (Mead et al., 2015; McCarthy et al., 2022).

Splitting $k$ into its components $k_{1}$ and $k_{2}$ in Eq. (6), $P(k)$ can be found in 2D space. We relate the flat-sky 2D wavenumber $k$ to the angular multipole $\ell$ as

$$\ell=\frac{2\pi}{bs}k,$$ (23)

where $b$ is the box width in pixels and $s$ is the pixel size in radians. Using CCL, we then compute the angular power spectrum $C_{\ell}$ at the corresponding cosmology and perform a scalar transformation using

$$P(k)=\frac{b}{s}C_{\ell}.$$ (24)

The noise matrix $\mathbf{N}$ in Eq. (9) is a uniform matrix with value $\sigma$ from Eq. (22). It describes the statistics of the noise in the convergence map.

We train and test the DeepMass model on realisations from a single (fiducial) cosmology and apply it to the nearby non-fiducial cosmologies. To test potential systematic uncertainties from this approach, we will also compare it to models trained directly on the non-fiducial data.

DeepMass does not include a padding functionality in the model, so the output map does not have the same dimensions as the input map. As such, we pad input data to reduce edge effects on the signal. To account for any correlation in signal between different redshifts, we have extended DeepMass to a 3D model across all four redshift bins. We note sporadic issues where DeepMass yields NaN values. We then rerun the model at the same training cosmology.

3.2 Observables

The functional integrals from Eq. (17), (18), and (19) are converted to sums over the number of pixels $N$, so we can calculate MFs from the 2D reconstructed CosmoGridV1 convergence fields:

$$V_{0}(t_{j})=\frac{1}{N}\sum_{i}\Theta(\kappa(\textbf{x}_{i})-t),$$ (25)

$$V_{1}(t_{j})=\frac{1}{4N}\sum_{i}\Delta(\kappa(\textbf{x}_{i})-t)\sqrt{\kappa_{x}^{2}+\kappa_{y}^{2}},$$ (26)

$$V_{2}(t_{j})=\frac{1}{2\pi N}\sum_{i}\Delta(\kappa(\textbf{x}_{i})-t)\\ \left(\frac{2\kappa_{x}\kappa_{y}\kappa_{xy}-\kappa_{x}^{2}\kappa_{yy}-\kappa_{y}^{2}\kappa_{xx}}{\kappa_{x}^{2}+\kappa_{y}^{2}}\right),$$ (27)

where $\Delta$ is 1 when $\kappa(\textbf{x}_{i})$ is between the thresholds $t_{j}$ and $t_{j+1}$ and 0 outside the range. Following Grewal et al. (2022), we evenly space the thresholds from $\mu-3\sigma$ to $\mu+3\sigma$, where $\mu$ is 0 and $\sigma$ is the field value standard deviation for each redshift in each cosmology. We differ slightly here, as in that work we chose thresholds for each individual map realisation. We note the thresholds are dynamic across reconstruction methods; we expect this to slightly improve our posterior contours compared to other approaches.

3.3 Evaluating the Posterior

We measure MFs from the pixels in the reconstructed posterior distributions of each method, then build a Fisher matrix $\mathrm{F}$ to evaluate the constraining power on cosmological parameters:

$$\mathrm{F}_{ij}=\sum_{mn}\frac{\partial X_{m}}{\partial\theta_{i}}\mathrm{C}^{-1}_{mn}\frac{\partial X_{n}}{\partial\theta_{j}},$$ (28)

where $X_{m}$ is the measured MF, $\theta_{i}$ is the cosmological parameter, and the covariance matrix $C_{mn}$ is $\langle(X_{m}-\overline{X}_{m})(X_{n}-\overline{X}_{n})\rangle$. The covariance is obtained from simulations at the fiducial cosmology. The Fisher derivatives are taken from the slope of the line of best fit through the three cosmologies in Table 1: minus parameter step size, fiducial, and plus parameter step size. The paired initial condition seeds decrease the noise in this calculation. To match DES Y3 constraining power, we divide the covariance matrix by the square root of the ratio of DES sky area (5000 sq deg) to our map size (15$\times$15 sq deg); this scale factor is 4.714 (Gatti et al., 2021b).

The correlation matrix derived from the covariance matrix is shown in Figure 3. While $V_{0}$ has some correlation with itself, it is anti-correlated with $V_{1}$ and $V_{2}$, which are highly correlated with each other. There is also strong correlation across redshift bins. We note the difference between this matrix and our previous MF correlation matrix in Grewal et al. (2022) can be attributed to the thresholding methodology.

Figure 3: Correlation coefficient matrix of all MFs, which have been concatenated from all four redshift bins. There is strong correlation between the MFs across redshift. The boxes represent one functional, where there are three per bin (denoted by a thicker black line). The MFs have been measured from a true convergence map at fiducial cosmology.

With the center set at the fiducial cosmology, we follow the process in Coe (2009) to generate the ellipse parameters and build the contours in Section 4. We marginalise over the other parameters by inverting the Fisher matrix before removing them.

3.4 Figures of Merit

We use two figures of merit (FOM) to measure the constraining power of different reconstruction methods: the lensing amplitude parameter $S_{8}$ and the inverse of the area of the $\Omega_{\mathrm{m}}$ - $\sigma_{8}$ contour ellipse. $S_{8}\equiv\sigma_{8}\left(\frac{\Omega_{\mathrm{m}}}{0.3}\right)^{\alpha}$ (where $\alpha$=0.5) is perpendicular to the lensing degeneracy direction for $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ for typical survey redshifts. It is thus the minor axis of the $\Omega_{\mathrm{m}}$ - $\sigma_{8}$ contours, and is usually the best constrained parameter combination.

Following Euclid Collaboration et al. (2023), we can transform the Fisher matrix in $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ into the matrix for $\Omega_{\mathrm{m}}$ - $S_{8}$ with

$$\mathrm{F}(\Omega_{\mathrm{m}},\mathrm{S}_{8})=\mathrm{M^{T}}\mathrm{F}(\Omega_{\mathrm{m}},\sigma_{8})\mathrm{M},$$ (29)

where

$$\mathrm{M}=\begin{bmatrix}\dfrac{\partial\Omega_{\mathrm{m}}}{\partial\Omega_{\mathrm{m}}}&\dfrac{\partial\Omega_{\mathrm{m}}}{\partial\mathrm{S}_{8}}\\[10.0pt] \dfrac{\partial\sigma_{8}}{\partial\Omega_{\mathrm{m}}}&\dfrac{\partial\sigma_{8}}{\partial\mathrm{S}_{8}}\end{bmatrix}=\begin{bmatrix}1&0\\ -0.5\dfrac{\sigma_{8}}{\Omega_{\mathrm{m}}}&\left(\dfrac{0.3}{\Omega_{\mathrm{m}}}\right)^{0.5}\end{bmatrix}.$$ (30)

In this work, we also aim to study the uncertainty in the DeepMass reconstruction. We quantify it with $\chi^{2}$ per pixel:

$$\chi^{2}=\frac{1}{N_{\mathrm{pixels}}}\sum\left(\frac{X_{\mathrm{recon}}-X_{\mathrm{true}}}{\sigma}\right)^{2},$$ (31)

using the $\sigma$ defined in Eq. (22). $X_{\mathrm{true}}$ is the true MF observable and $X_{\mathrm{recon}}$ is the reconstructed MF observable, fiducial or non-fiducial.

4 Results

Figure 4: Visualising the output for true maps, true smoothed maps, and the four different reconstruction methods. 15 $\times$ 15 degree flat sky convergence maps of redshift bin 3 at fiducial cosmology. The shared colourbar shows convergence of each pixel. A Gaussian kernel with a standard deviation of one arcmin has been applied to the true map to achieve the smoothed map.

Figure 4 shows the true, true smoothed, and reconstructed convergence maps for the four methods: Kaiser Squires, Wiener filter, DarkMappy, and DeepMass. As expected, peaks in the KS and WF maps have been smoothed out because small scales are noise-dominated; however, DarkMappy looks similarly noise-dominated despite it being a non-linear method. DeepMass appears most like the truth convergence map, which is evidence that it has denoised the convergence signal; this is consistent with predictions. While DeepMass takes more time to train the model than the other methods do to run, once the model is generated, DeepMass reconstructs convergence maps very quickly.

Since MFs are not linear, their total value does not come from summing the MFs from noise and signal. Here we have a separate analysis with signal and noise to study the effects of the latter. The Minkowski functional outputs for true, true-plus-noise, KS, WF, DarkMappy, and DeepMass are compared in Figure 5. The MF curves have the same shape, which is the typical shape for a convergence field. The distributions are not expected to be exactly the same, and here we see differing amplitudes for different methods corresponding to the noise level in the field. We do not expect KS to reduce the noise in any way, but some of the noise in KS has been smoothed out. As such, KS is not as noisy as true-plus-noise. Figure 5 shows KS, WF, and DarkMappy are comparable, while DeepMass has the lowest variation in amplitude due to noise suppression after smoothing.

Figure 5: The first three MFs, $V_{0}$, $V_{1}$, and $V_{2}$, calculated from fiducial convergence maps. The curves shown are for true signal and true-plus-noise, along with the four reconstructed signals visualised in Figure 4: Kaiser Squires, Wiener filter, DarkMappy, and DeepMass. The $t$ values on the curve correspond to an excursion set, and the three MFs on the $y$-axis describe the perimeter, area, and curvature of the excursion set of the convergence.

We marginalise over all cosmological and baryonic parameters in Table 1 in our analysis except for the two used in our Fisher matrix. Figure 6 shows the contour comparison in the reconstruction methods for $\Omega_{\mathrm{m}}$ - $\sigma_{8}$, $\Omega_{\mathrm{m}}$ - $S_{8}$, and $\Omega_{\mathrm{m}}$ - $w_{0}$. The contour directions are typical for convergence constraints. Figures of merit are shown in Table 3 with the standard deviation of $S_{8}$ in the first column and the inverse of the $\Omega_{\mathrm{m}}$ - $\sigma_{8}$ contour area in the second column.

Method	$\sigma_{S_{8}}$	1/$E_{A}$
True	0.00392	28328.11
True + Mask	0.00387	19445.44
True + Noise	0.00819	3378.91
Kaiser Squires	0.00848	3484.66
Wiener Filter	0.00700	3994.15
DarkMappy	0.00940	2967.18
Deepmass	0.00409	6606.95

Table 3: Figures of merit for different reconstruction methods. $E_{A}$ is the ellipse area of the $\Omega_{\mathrm{m}}$ - $\sigma_{8}$ contour.

Figure 6: $w$CDM Fisher constraints for true signal, true-plus-mask, true-plus-noise, KS reconstruction, WF reconstructions, DarkMappy reconstruction, and DeepMass reconstruction. The top left plot shows constraints on $\Omega_{\mathrm{m}}$ vs $w_{0}$, the bottom left plot shows constraints on $\Omega_{\mathrm{m}}$ vs $\sigma_{8}$, and the bottom right plot shows constraints on $\Omega_{\mathrm{m}}$ vs $S_{8}$. KS, WF, and DarkMappy are comparable at these scales, while DeepMass shows significantly tighter constraints, implying the model has effectively denoised the signal.

DeepMass outperforms the other reconstruction methods by roughly 50$\%$ and has constraints closest to the true contour. KS, WF, and DarkMappy have similar constraints within a 30% interval of each other, and are comparable to true-plus-noise, which shows adding noise doubles the size of the contour. We note the choice of threshold range can affect the size of the contours. As mentioned previously, KS is not expected to perform as well as true-plus-noise, but some of the noise has been smoothed out. We can see in Figure 6 that the step of masking the data has a much smaller impact than the step of adding noise to the fields, the latter doubling the size of the parameter constraints. We use standard value 0.5 for $\alpha$ in the $\sigma_{S_{8}}$ calculation, which corresponds to the width perpendicular to the direction of the true-plus-mask contour. We also include the second column for inverse ellipse area to more completely quantify constraining power of cosmological parameters, and we see that the true map has the largest value, followed by true-plus-mask. On the whole, $\Omega_{\mathrm{m}}$ is better constrained than other cosmological parameters. We can see the same trends in the $\Omega_{\mathrm{m}}$ - $S_{8}$ and $\Omega_{\mathrm{m}}$ - $w_{0}$ planes. We investigate baryonic parameters in Appendix A.

Using Eq. (31), we compare $\chi^{2}$ values per pixel to measure the effect of training DeepMass at a (slightly) wrong cosmology. Figure 7 shows a comparison of DeepMass outputs for the fiducial and non-fiducial models (-$\Delta\Omega_{\mathrm{m}}$) applied to non-fiducial data (-$\Delta\Omega_{\mathrm{m}}$). Each histogram displays the $\chi^{2}$ values across 5000 output maps. The incorrect model histogram illustrates the distinction between the fiducial model and true data with -$\Delta\Omega_{\mathrm{m}}$, while the correct model histogram depicts differences for the non-fiducial model on the same input data. As expected, the histogram for the correct model has a lower average $\chi^{2}$ value. The quantification of uncertainty from the reconstruction is determined as the difference between the means of the two histograms: 3.8%. While this seems relatively small, it leads to changes in the MFs comparable to cosmologically-induced changes, as seen in Figure 8. As this can only worsen the constraining power of DeepMass, which is already close to the ideal constraints from truth data, it cannot have had a significant effect on our results here, as we focus on the reconstruction precision (not accuracy). However, this demonstrates that the selection of training model cosmologies is critical for future analysis accuracy. Tests in Jeffrey et al. (2020) used the truth cosmologies for optimal training, but analyses of real data will need thoughtful and perhaps iteratively-chosen training samples and careful emulators to interpolate between them.

Figure 7: Normalised $\chi^{2}$/pixel comparison of DeepMass reconstruction using the correct model versus the incorrect model. The correct model (purple) gives a lower $\chi^{2}$ value than the incorrect model (orange), indicating that choice of model is important. We measure a 3.8% uncertainty in the cosmology choice in DeepMass training, which has been calculated as the percent change in the difference in the means of the two $\chi^{2}$ histograms.

Figure 8: MFs for the two DeepMass outputs in Figure 7. Errorbars are shown on the output from the correct model.

5 Conclusion

We want to use the most precise reconstructed convergence maps in order to best capture the original signal. In this way, measuring statistics from such detailed reconstructed maps enables improved measurement of cosmological parameters. In this paper we investigate four non-parametric mass mapping methods developed to reconstruct convergence from shear: Kaiser Squires, Wiener filter, DarkMappy, and DeepMass. We smooth masked and noisy shear maps and reconstruct convergence maps using the different methods. KS requires smoothing because it can only take a continuous field as input. We then measure Minkowski functionals from the reconstructed convergence maps and perform a Fisher analysis to make contours. We compare constraints on $\Omega_{\mathrm{m}}$, $\sigma_{8}$, $S_{8}$, and $w_{0}$ (see Appendix A for baryonic parameter analysis).

The contour plots for different combinations of cosmological parameters all show the same results for different reconstruction methods. DeepMass outperforms the rest of the methods by a significant degree, demonstrating the model is a successful denoiser and an effective tool for reconstructing convergence from shear. Beyond a 50$\%$ improvement in constraining power, the model is also computationally efficient.

We also measure the uncertainty from using different cosmologies in the DeepMass model training. By comparing $\chi^{2}$ values of true data MFs and MFs from data trained on the correct model or the same data trained on an incorrect model, we can quantify the impact of cosmological choices on model performance with the difference in average $\chi^{2}$ values. This is about 3.8%, indicating the choice of training model cosmology has an impact on reconstruction precision.

This difference aligns closely with observed variations in average MF DeepMass outputs from data trained on the correct or incorrect model. Both effects are comparable with the change in cosmology for $\Omega_{m}$. These findings demonstrate a strong dependence of DeepMass on the cosmology of the training dataset, emphasising the need for the consideration of input cosmology in model training to ensure accurate and reliable results.

The remaining methods, KS, WF, and DarkMappy are comparable at the scales we use here. DarkMappy was expected to outperform KS and WF. One reason it did not could be that it was designed for cluster analysis and may be tuned for such scenarios rather than our wide-field tests; alternative choices of its parameters and wavelet dictionary could improve its performance. Another could be the resolution of the field. In this paper we have used CosmoGridV1 simulations with DES Y3-like redshift bins and noise, which have led to a more limited resolution. We expect an improvement in resolution with less noisy survey conditions like for LSST or in future LSST-like forecasts, as well as with newer simulations.

Additionally, the choice of threshold levels can influence contour sizes. Here we use dynamic thresholds, where the range is set by the mean and standard deviation of the field for each cosmology, redshift, and reconstruction method. Using fixed thresholds more carefully tuned to each case has less variability. This can lead to smaller contours, but the thresholds are less consistently useful for different inputs like redshift bins, reconstruction methods, and simulation conditions like noise level.

We use Minkowski functionals to evaluate the different reconstruction methods, but there are many other statistics that could be used. Other higher order statistics may yield different results and show more distinction between the contours. Likewise, we have chosen four reconstruction methods here, representative of the different approaches to mass mapping. Other methods may be as or more precise.

In future analysis, using different (or multiple statistics) could lead to better information about how the different methods perform. Additionally, using lower noise conditions like LSST and other reconstruction methods could lead to stronger differentiation among methods.

We finally note that mass mapping for higher-order statistics is an ideal problem for a community challenge like the ones presented in Zuntz et al. (2021) and Leonard et al. (2023). This would both focus authors of reconstruction methods on useful metrics for upcoming surveys and incentivise them to optimise parameters and choices for a specific scenario.

6 Acknowledgements

We thank Tomasz Kacprzak for help with the CosmoGridV1 simulations. NG thanks Harry Rendell-Bhatti for the useful discussions. TT acknowledges funding from the Swiss National Science Foundation under the Ambizione project PZ00P2_193352. Results in this paper made use of many software packages, including Numpy, Scipy, and CCL (Harris et al., 2020; Virtanen et al., 2020; Chisari et al., 2019).

Appendix A Effects of Baryons

A.1 Comparing Fisher calculations for Baryonic Parameter Constraints

Here we explore the effects of baryons in our simulations and analysis. Figure 9 shows contour plots for baryonic parameters. The left plot shows the $M^{0}_{\mathrm{c}}-\nu$ plane, and the right plot shows the $\sigma_{8}-M^{0}_{\mathrm{c}}$ plane. The results are unexpected, as the true maps should have the tightest constraints. We would expect smoothing and noise in the reconstruction process to lead to a loss of information, especially for small-scales where we expect to see baryonic effects. These contours are not reliable. To investigate this discrepancy, we study the baryonic effects at each step of the analysis.

Figure 9: Contours for different baryonic parameter combinations. The same four reconstruction methods are compared, which produce similar constraints, with DeepMass performing the best.

A.2 Baryonic non-fiducial Convergence Maps

We first compare the true baryonic and fiducial convergence maps to visualise the effect of a change in the baryonic parameters $M_{\mathrm{c}}^{0}$ and $\nu$. Figure 10 shows the difference between the non-fiducial $\pm\Delta M_{\mathrm{c}}^{0}$ and $\pm\Delta\nu$ baryonic convergence maps and the fiducial convergence maps. We can see that the change in signal is very small compared to the original convergence value. This is roughly two orders of magnitude smaller than the change for nonbaryonic parameters.

Figure 10: The difference in the non-fiducial and fiducial convergence maps for baryon parameters $M_{\mathrm{c}}^{0}$ and $\nu$.

A.3 Minkowski functionals for Baryonic Parameters

We next compare the Minkowski functional observables measured from the baryonic and fiducial convergence maps. In Figures 11 and 12, the functionals at different baryon cosmologies are shown to be nearly identical. This is not the case for the cosmological parameters we constrain in this paper; an example of variation in MF observables for $\Omega_{\mathrm{m}}$ is shown in Figure 13. Here, any difference in the MFs falls within the noise regime. This may be because the $\Delta$ parameter step size is too small or the pixel size is too large to visualise or measure baryonic effects. These changes in step size are small enough that numerical noise dominates true Fisher matrix contours.

Figure 11: MFs at different values of $\nu$ (top row). Difference in MFs of $\Delta\nu$ and fiducial maps (bottom row).

Figure 12: MFs at different values of $M_{\mathrm{c}}^{0}$ (top row). Difference in MFs of $\Delta M_{\mathrm{c}}^{0}$ and fiducial maps (bottom row).

Appendix B Minkowski functionals for Non-Baryonic Parameters

For non-baryonic parameters, we see the expected change in the Minkowski functionals that corresponds to a change in cosmology. Figure 13 shows MFs for $\pm\Delta\Omega_{\mathrm{m}}$ as well as the fiducial value for $\Omega_{\mathrm{m}}$.

Figure 13: MFs at different values of $\Omega_{\mathrm{m}}$; the change is comparable to the variation from training cosmology in Figure 8.