Three-Dimensional Reconstruction of Weak Lensing Mass Maps with a Sparsity Prior. I. Cluster Detection

We propose a novel method to reconstruct high-resolution three-dimensional mass maps using data from photometric weak-lensing surveys. We apply an adaptive LASSO algorithm to perform a sparsity-based reconstruction on the assumption that the underlying cosmic density field is represented by a sum of Navarro-Frenk-White halos. We generate realistic mock galaxy shape catalogues by considering the shear distortions from isolated halos for the configurations matched to Subaru Hyper Suprime-Cam Survey with its photometric redshift estimates. We show that the adaptive method significantly reduces line-of-sight smearing that is caused by the correlation between the lensing kernels at different redshifts. Lensing clusters with lower mass limits of $10^{14.0} h^{-1}M_{\odot}$, $10^{14.7} h^{-1}M_{\odot}$, $10^{15.0} h^{-1}M_{\odot}$ can be detected with 1.5-$\sigma$ confidence at the low ($z<0.3$), median ($0.3\leq z<0.6$) and high ($0.6\leq z<0.85$) redshifts, respectively, with an average false detection rate of 0.022 deg$^{-2}$. The estimated redshifts of the detected clusters are systematically lower than the true values by $\Delta z \sim 0.03$ for halos at $z\leq 0.4$, but the relative redshift bias is below $0.5\%$ for clusters at $0.4<z\leq 0.85$. The standard deviation of the redshift estimation is $0.092$. Our method enables direct three-dimensional cluster detection with accurate redshift estimates.


INTRODUCTION
Weak gravitational lensing causes small, coherent distortion of the shapes of distant galaxies. Information on the foreground mass distribution is imprinted in the distorted galaxy images, and thus weak-lensing offers a direct physical probe into the mass distribution in our universe, including both visible matter and invisible dark matter (see Kilbinger 2015;Mandelbaum 2018, for recent reviews). Ongoing and future observational programs such as the Subaru Hyper Suprime-Cam survey (HSC; Aihara et al. 2018), the Kilo-Degree Survey  (Laureijs et al. 2011), and the Nancy Grace Roman Space Telescope (Spergel et al. 2015) are aimed at studying largescale mass distribution with high precision.
Statistics of density peaks in two-dimensional (2D) and three-dimensional (3D) mass maps can be used as a powerful cosmological probe (Jain & Van Waerbeke 2000;Fan et al. 2010;Lin et al. 2016). One can detect massive clusters of galaxies by identifying high signal-tonoise (S/N) peaks in a mass map without any reference to the mass-to-light ratio (Schneider 1996;Hamana et al. 2004).
2D density reconstruction techniques recover integration of the projected mass along the line of sight, and have been extensively studied so far (Kaiser & Squires 1993;Lanusse et al. 2016;Price et al. 2020) and applied to large-scale surveys (Oguri et al. 2018;Chang et al. 2018;Jeffrey et al. 2018). Cluster detection and identification from 2D mass maps has been applied to wide-field weak-lensing surveys (Shan et al. 2012;Miyazaki et al. 2018a;Hamana et al. 2020). Cross-matching with another cluster catalogs (e.g., an optically selected one) is usually needed to infer physical quantities such as mass and to estimate redshift for individual clusters located in 2D mass maps.
In principle, one can directly reconstruct 3D mass distributions by using photometric redshift (photo-z) information of the source galaxies (Hu & Keeton 2002;Bacon & Taylor 2003;Massey et al. 2007;Simon et al. 2009;VanderPlas et al. 2011). Unfortunately, these methods either do not have enough spatial resolution to identify individual clusters, or suffer from smearing along the line of sight. These are critical obstacles that need to be overcome for practical searches of clusters in 3D mass maps. Alternatively, Hennawi & Spergel (2005) propose to perform a maximum-likelihood detection of clusters, by convolving tomographic shear measurements with 3D filters that match the tangential shears induced by multiscale Navarro-Frenk-White (NFW) halos. Their method can be used effectively to detect clusters, but does not fully reconstruct wide-field mass distributions.
In the present paper, we develop a novel method for high-resolution 3D reconstruction. We model a given 3D density field as a sum of the NFW (Navarro et al. 1997) basis 'atoms' that are setup in comoving coordinates. A basis atom is defined by a 2D NFW surface density profile on the transverse plane and one-dimensional Dirac delta function in the line of sight direction. We apply the adaptive LASSO algorithm (Zou 2006) to find a sparse solution for a pixelized map. We examine the performance of cluster detection using the reconstructed mass maps. To this end, we apply shear distortions generated by isolated halos using realistic HSC-like galaxy shapes with photo-z estimates.
The rest of the paper is organized as follows. In Section 2, we propose the new method for 3-D density map reconstruction. In Section 3, we study the cluster detection from the reconstructed mass map using isolated halo simulations with the HSC observational condition. In Section 4, we summarize and discuss the future development of the method.

METHOD
The lensing shear field γ observed from background galaxy images is related to the foreground density contrast field δ = ρ/ρ − 1 through a linear transform: where is the shear measurement error due to the random orientation of galaxy shapes (shape noise) and the sky variance (photon noise). The matrix operator T = P·Q includes the physical lensing effect denoted by a matrix operator Q and the systematic effects in observations represented by a matrix operator P. The latter includes smoothing of the shear field in the transverse plane and photometric redshift uncertainties.
In this section, we first introduce the NFW dictionary to model the density contrast field (Section 2.1). Then we explain the weak-lensing operator Q in Section 2.2, whereas the systematics operator P is introduced in Section 2.3. Reconstruction of the density contrast field is performed by solving a linear problem of Equation (1) or of the equivalent Equation (10) in its practical form introduced later in this section. To this end, we devise an adaptive LASSO algorithm that achieves high-resolution reconstruction in Section 2.4.

Model dictionary
In order to reconstruct high-resolution, high-S/N mass maps, we incorporate prior information on the density contrast field into the reconstruction by modeling the density field as a sum of basis atoms in a "dictionary": where Φ is the matrix operator transforming from the projection coefficient vector x to the density contrast field δ. The column vectors of Φ are the basis "atoms" of the model dictionary and denoted as φ s . There have been a few studies that adopt different dictionaries for 3D weak-lensing map reconstructions: (i) Simon et al. (2009) perform a reconstruction in Fourier space, which is equivalent to representing the density contrast field with sinusoidal functions. However, sinusoidal functions are not localized in configuration space, and the density contrast field is not sparse in Fourier space. Therefore, sparsity priors cannot be directly applied to this model dictionary for high-resolution mass map reconstructions. (ii) Leonard et al. (2014) model the density contrast field with starlets (Starck et al. 2015). However, the starlet dictionary does not account for the angular scale difference at different lens redshifts, and is not specifically designed to model the clumpy mass distribution in the universe. It is worth exploring other dictionaries for our purpose of weak-lensing mass reconstruction.
The normalized 2D profiles of the smoothed basis "atoms". The pixel size is 1 . The leftmost column is the point-mass atom, and the other columns show the NFW atoms with different scale radii (in pixels) as indicated. The upper (lower) panels show the basis atoms in Fourier (configuration) space. We smooth the 2D profiles using a Gaussian kernel with a 1.5 pixel width.

1D normalized profile
Point mass NFW r s = 1 NFW r s = 2 NFW r s = 4 NFW r s = 8 Figure 2. The normalized one-dimensional (1D) profiles of smoothed basis atoms, which are slices of the corresponding 2D profiles (shown in Figure 1) at y = 0. The scale radii are in pixels.
In the standard cosmological model, dark matter is concentrated in roughly spherical "halos", which have the NFW density profile (Navarro et al. 1997). Motivated by this fact, we generate a model dictionary using NFW atoms with N typical scale radii r s (s = 1, ..., N ) in comoving coordinates. An atom has the NFW surface density profile on the transverse plane (Takada & Jain 2003) and the Dirac δ function in the line of sight direction. Following Leonard et al. (2014), we neglect the size of halos along the line of sight since the resolution scale of the reconstruction is much larger than the halos.
The multiscale NFW atom is defined as where r θ is the projection of the comoving position on the transverse plane, (4) f = 1/[ln(1 + c) − c/(1 + c)], and c is the halo concentration. For simplicity, we fix c = 4 for the NFW atoms in our dictionary. In the present paper, we adopt a hard truncation on the NFW profile at a radius equals cr s . Studying the influence of different truncation forms (Oguri & Hamana 2011) on the mass map reconstruction is left to our future work. Compared to other model dictionaries, our NFW dictionary is motivated by a physical consideration on the clumpy mass distribution in the universe. Furthermore, the multiscale NFW atoms are set up in comoving coordinates, with an account of the scale difference in the angular coordinates for halos at different lens redshifts. The corresponding NFW atom in the angular separation coordinates is where χ(z) is the comoving distance to redshift z.
For the NFW dictionary, the transform from the projection coefficient vector to the density contrast field of Equation (2) is written as To simplify the notation, we compress the projection coefficient vectors on the basis atoms with multiple scales into a single column vector: and write the dictionary transform operator into a row vector: For an additional test and comparison, we also construct a dictionary with the point-mass atoms which are represented by the 3D Dirac function as The 2D profiles of the NFW atoms and the point-mass atom on the transverse plane are shown in Figure 1. The corresponding 1D profiles are plotted in Figure 2. The profiles are smoothed with a Gaussian kernel and pixelized on linearly spaced grids. Details on the smoothing and pixelizing operations are described in Section 2.3.2. We define the forward transform operator A = P·Q·Φ where P represents systematic effects and Q represents the physical lensing effect. With Equations (1) and (2), we write the transform from the coefficient vector x to the observed lensing shear field as  Top panel: Normalized lensing kernels as a function of source redshift with lens redshifts fixed. The solid lines are the kernels for the source galaxies with precise spectroscopic redshifts, whereas the dashed lines are for source redshifts with HSC-like photometric redshift errors. Bottom panel: Correlation matrix between the lensing kernels of different lens redshifts. We normalize the diagonal terms to 1 . The color map is the correlation matrix for spectroscopic redshift. We also compare the result for the lensing kernel of spectroscopic (photometric) redshift by solid (dashed) contours at levels 0.7, 0.85, and 0.98 .

Weak gravitational lensing
The lensing transform operator can be expressed as K(z l , z s ) is the lensing kernel between lens redshift z l and source redshift z s (Bartelmann & Schneider 2001), which is given by where E(z) is the Hubble parameter as a function of redshift, in units of H 0 .
is the Kaiser-Squares kernel (Kaiser & Squires 1993), which decays proportional to the inverse-square of the distance. Here θ 1,2 are the two components of the angular position vector θ.
The top panel of Figure 3 shows the lensing kernels as a function of source redshift for lenses at different lens redshifts. The bottom panel of Figure 3 shows the correlation between the lensing kernels. Each lensing kernel is highly nonlocal, and the lensing kernels of different lens redshifts are strongly correlated as can be seen in the bottom panel. These inherent properties of the lensing kernels result in strong correlation between the column vectors that constitute the forward transform matrix A. As will be discussed in Section 2.4, the strong correlation makes it challenging to reconstruct mass maps with high-resolution in the line of sight direction.

Systematics
Shear measurement deviates from the true, physical shear owing to a variety of systematic effects in real observations. In this section, we discuss the influence of several major systematics on the lensing shear measurement, and describe the corresponding transform operator by decomposing into three parts of R (photometric redshift uncertainties), W (smoothing), and M (masking).

Photometric redshift Uncertainty
Photometric redshifts (photo-z) are estimated with a limited number of broad optical and near-infrared bands in the current generation weak-lensing surveys. For example, 9 bands are used for KIDS+VIKING survey (Hildebrandt et al. 2020), 5 bands for DES survey and HSC survey. Unlike the high-precision spectroscopic redshift (spect-z) estimation, the photo-z estimation suffers from considerable statistical uncertainties.  are divided into ten bins according to the photo-z bestfit estimates. Figure 5 shows the average probability density function (PDF) for galaxies in each redshift bin.
In the presence of photo-z uncertainties, a source galaxy with a best-fit photo-z estimate z s has a posterior probability P (z|z s ) of being actually located at redshift z. This means that there is a possibility P (z|z s ) that the galaxy image is actually distorted by the shear at redshift z. Therefore, the photo-z uncertainty smears the lensing kernel statistically, which we model by a smearing operator R = dzP (z|z s ) .
(14) Figure 3 shows the lensing kernels and their correlations for source redshifts with photo-z uncertainties demonstrated in Figure 5. Compared to the lensing kernels of spect-z that converge to zero at source redshifts lower than the lens redshift, the lensing kernels of photo-z do not converge to zero at the low source redshifts. This is simply because the galaxies with photo-z estimated lower than the lens redshifts may actually be located at higher redshifts. In addition, the photo-z uncertainty only slightly increases the correlations between lensing kernels at different lens plane.

Smoothing
The source galaxies are not uniformly distributed in the sky, with substantial fluctuations in the number density. We first smooth the lensing shear measurements, and pixelize the smoothed shear field onto a regular grid. After these procedures, the fast Fourier transform (FFT) can be directly conducted on the transverse plane in each redshift bin to compute the convolution operation in Equation (1). Another benefit of smoothing is that it reduces bias arising from the aliasing effect in the pixelization process since the smoothing kernel reduces the amplitude of the shear signal at high frequency.
We convolve the lensing shear measured from galaxy images with a 3D smoothing kernel W ( θ, z). The expectation of the lensing shear field after smoothing is where z i and θ i refer to the photometric redshift, and transverse position of the ith galaxy in the galaxy sample, respectively. γ L refers to the shear field before the smoothing.
The smoothing kernel W ( θ, z) can be decomposed into a transverse component W t ( θ) and a line of sight component W l (z) as We use an isotropic 2D Gaussian kernel and a 1D tophat kernel to smooth the shear field in the transverse plane and in the line of sight direction, respectively. These components are given by where we set β = 1.5 arcmin in this paper. By definition, the smoothing kernel is normalized as Assuming that the galaxy number distribution varies slowly at the smoothing scale, the smoothed shear can be written into where W is the smoothing operator defined as We note that another widely used scheme is to average the shear measurements in each pixel. Such a scenario is equivalent to resampling the shear field smoothed with a 3D top-hat kernel with the same scale as the pixels.
After smoothing, we pixelize the shear field on N x × N y × N z grids, where N x and N y are the number of pixels of the two orthogonal axes of the transverse plane, and N z is the number of pixels in the line of sight direction. We denote by γ α the smoothed shear measurements recorded on the pixel with index α, where 1 ≤ α ≤ N x × N y × N z . The grids on the transverse planes are equally spaced with a pixel size of 1 . In the line of sight direction, we set binning with equal galaxy number as shown in Figure 4.
Similarly, we pixelize the projection coefficient vector x onto an N x × N y × N l grid. In this paper, the projection coefficient vector is pixelized in equal spacing ranging from redshift 0.01 to redshift 0.85. Here, we use N l to denote the number of the lens planes and x β to denote the projection field element with index where the last N is the number of NFW dictionary frames representing different physical scale radii. The corresponding pixelized elements of the forward transform matrix A is denoted as A αβ .
In order to see the effect of smoothing clearly, we plot the histograms of the statistical shear measurement errors due to shape noise and sky variance for the galaxies in tract 9347 of the first-year HSC shear catalog ). In Figure 6, we compare the error on an individual galaxy basis and that of individual pixels after the smoothing and pixelization procedures described in this section. The standard deviation of the individual pixel errors is much smaller than for the galaxies owing to the smoothing. We also note that, even though the shear measurement error for the galaxies does not follow a Gaussian distribution, the error after smoothing and pixelization is well described by a Gaussian distribution, which is simply explained by the central limit theorem.

Masking
In real observations, shear measurements can be performed in a finite region of the sky, and the boundaries often have complicated geometries. Moreover, many isolated sub-regions near the bright stars are masked out.
We define the masking window function as where n sm is the smoothed galaxy number density. We define the masking operator as where δ D ( r) is 3D Dirac delta function.
Taking into account the systematics discussed in previous Sections 2.3.1-2.3.3, the systematic operator is

Density map reconstruction
The projection coefficients can be estimated by optimizing a penalized loss function. An estimator is generally defined aŝ where Σ (γ − Ax) 2 2 is the l 2 norm of residuals weighted by the inverse of the covariance matrix Σ of the shear measurement error, which measures the difference between the prediction and the data. C(x) is the penalty term measuring the deviation of the coefficient estimate x from the prior assumption. The estimation with the "penalty" term prefers parameters that are able to describe the observation with a specified prior information. The penalty parameter λ adjusts the relative weight between the data and the prior assumption in the optimization process.
There have been a few studies that adopt different penalties for 3D weak-lensing map reconstructions: (i) Simon et al. (2009) propose to use the Wiener filter, which is also known as l 2 ridge penalty (C = x 2 2 ), to find a penalized solution in Fourier space. Oguri et al. (2018) apply the method of Simon et al. (2009) to the first-year data of the Hyper Suprime-Cam Survey (Aihara et al. 2018). It is found that the density maps reconstructed by the method suffer from significant line of sight smearing with standard deviation of σ z = 0.2 − 0.3. (ii) Leonard et al. (2014) use GLIMPSE algorithm, which adopts a derivative version of the l 1 LASSO penalty (C = x 1 1 ) to find a sparse solution in the starlet dictionary (Starck et al. 2015). GLIMPSE reduces the smearing by adopting a "greedy" coordinate descent algorithm that forces the structure to grow only on the most related lens redshift plane.
In this paper, we use another derivative version of the LASSO penalty -the adaptive LASSO penalty. We first normalize the column vectors of the forward transform matrix in Section 2.5. We then introduce the loss function with the adaptive LASSO penalty in Section 2.5.1. Finally, we find the minimum of the loss function in Section 2.5.2 with the FISTA algorithm (Beck & Teboulle 2009).

Normalization
The l 2 norm of the ith column vector of A weighted by the inverse of the noise covariance matrix is defined as The column vectors have different weighted l 2 norms. Since the gradient descent algorithm, which will be used to solve Equation (24) in Section 2.5.2, takes each column vector equally, we normalize the column vectors before performing the density map reconstruction to boost the convergence speed of the gradient descent iteration. The normalized forward transform matrix and projection parameters are given by

Adaptive LASSO
The LASSO algorithm uses l 1 norm of the projection coefficient vector to regularize the modeling. The estimator is defined aŝ where • 1 1 and • 2 2 refer to the l 1 norm and l 2 norm, respectively, and λ is the penalty parameter for the LASSO estimation.
The LASSO algorithm searches and selects the parameters that are relevant to the measurements, and simultaneously estimates the values of the selected parameters. It has been shown by Zou (2006) that when the column vectors of the forward transform matrix A are highly correlated, the algorithm cannot select the relevant atoms from the dictionary consistently. In addition, the estimated parameters are often biased owing to the shrinkage in the LASSO regression. We note that, for the density map reconstruction problem here, the column vectors are highly correlated even in the absence of photo-z uncertainties since the lensing kernels for lenses at different redshifts overlap significantly, i.e. highly correlated as shown in Figure 3. Therefore, the LASSO algorithm cannot precisely determine the consistent mass distribution in redshift, and the reconstructed map suffers from smearing in the line of sight direction even in the absence of noises. Figure 7 shows an example reconstruction results for a single halo with mass M 200 = 10 15 h −1 M at redshift 0.35 3 . The shear measurement and photo-z uncertainties are not included in this simulation. We find significant smearing of the mass distribution with the LASSO algorithm (left panel of Figure 7).
To overcome the problem, Zou (2006) proposes the adaptive LASSO algorithm, which uses adaptive weights to penalize different projection coefficients in the l 1 penalty. The adaptive LASSO algorithm performs a two-step process. In the first step, the standard (nonadaptive) LASSO is used to estimate the parameters. Let us denote the preliminary estimation asx LASSO . In the second step, the preliminary estimate is used to calculate the nonnegative weight vector for penalization aŝ where we set the hyperparameter τ to 2 . The adaptive LASSO estimator is then given bŷ where "•" refers to the element-wise product. λ ada is the penalty parameter for the adaptive LASSO, which does not need to be the same as the penalty parameter for the preliminary LASSO estimation λ. The adaptive weights enhance the shrinkage in the soft thresholding for the coefficients with smaller amplitudes, whereas the weights suppress the shrinkage for the coefficients with larger amplitudes.
To simplify the equations in the following, we rewrite the loss function with the Einstein notation as and the quadruple term in the loss function as 2.5.2. FISTA Beck & Teboulle (2009) propose the Fast Iterative Soft Thresholding Algorithm (FISTA) to solve the LASSO problem. Since the loss functions of the LASSO and the adaptive LASSO differ only in their penalty terms, FISTA is also applicable to solve the adaptive LASSO problem in a straightforward manner. In this paper, we apply FISTA to solve both the preliminary LASSO estimation and the adaptive LASSO estimation.
We first explain the LASSO preliminary estimation. The coefficients are initialized as x (1) i = 0. According to FISTA, we iteratively update the projection coefficient vector x. Taking the nth iteration as an example, a temporary update is first calculated as where ST is the soft thresholding function defined as The soft thresholding is a part of the LASSO algorithm, which selects the modes with amplitude greater than λ, and shrinks the selected estimation by λ. The coefficient µ is the step size of the gradient descent iteration, and ∂ i G(x (n) ) is the ith element of the gradient vector of G at point x (n) : where Re {•} returns the real part of the input vector. FISTA requires an additional update using a weighted average between x (n+1) and x (n) as t (n+1) = 1 + 1 + 4(t (n) ) 2 2 , where the relative weight is initialized as t (1) = 1.
FISTA converges as long as the gradient descent step size µ satisfies where the operation • returns the spectrum norm of the input matrix. The spectral norm is estimated by simulating a large number of random vectors with l 2 norms equal one with different realizations. The matrix operator A † · Σ −1 · A is subsequently applied to each random vector to yield a corresponding transformed vector. The spectral norm is determined as the maximum l 2 norm of the transformed vectors. As summarized in the below, we first initialize the projection coefficients as zero, and use FISTA to find the global minimum of the LASSO loss function. Thusobtained global minimum is the preliminary one. We then use the preliminary LASSO estimation to weight the coefficients and calculate the adaptive LASSO loss function. Finally, we set the preliminary LASSO estimation as the "warm" start of the adaptive LASSO estimation, and FISTA again to find the global minimum of the adaptive LASSO loss function, which is our final solution.

CLUSTER DETECTION
We simulate weak-lensing shear induced by NFW halos with various masses and redshifts. The generated shear fields are used to distort the HSC mock galaxy shapes with different realizations of shear measurement and photo-z uncertainties (Section 3.1).
We then test our algorithm using the mock catalogs with varying the penalty parameter in our model (Section 3.2). We also compare our method that uses the NFW dictionary with one using the point mass dictionary (Section 3.3).

Simulations
We use halos with a variety of masses and redshifts in the range 10 14 h −1 M < M < 10 15 h −1 M , and 0.05 < z < 0.85, respectively. We divide the parameter space into eight redshift bins and eight mass bins with equal separation. We randomly shift the input halo redshift and halo mass from the bin center by a small amount in order to avoid repeatedly sampling at the exact same halo mass and redshift.
The concentration parameter c h of a NFW halo is determined as a function of the halo mass (M 200 ) and redshift (z h ) according to Ragagnin et al. (2019): The weak-lensing shear fields of the NFW halos are calculated according to Takada & Jain (2003). The shear distortions are applied to one hundred realizations of galaxy catalogs with the HSC-like shear measurement error and photo-z uncertainty.
The mock galaxy catalogs are generated using the HSC S16A shear catalog . We use the galaxies in a 1 square degree field at the center of tract 9347 (Aihara et al. 2018). The average galaxy number density in this region is 22.94 arcmin −2 . The positions of galaxies are randomized to distribute homogeneously in the one-square degree stamp. We randomly assign redshift for each galaxy following the MLZ photo-z probability distribution function (Tanaka et al. 2018).
We simulate the HSC-like shear measurement error due to shape noise and sky variance with different realizations by randomly rotating the galaxies in the shear catalog. The shear measurement error on individual galaxy level and individual pixel level are demonstrated in Figure 6. The standard deviation map of the noise is demonstrated in Figure 8.

NFW atoms
In this section, we test the performance of our algorithm by adopting models where the matter density field is represented by multiscale NFW atoms. The dictionary is constructed with three frames of different NFW scale radii in comoving coordinate: 0.12 h −1 Mpc, 0.24 h −1 Mpc, and 0.36 h −1 Mpc. The truncation radii are set to four times the scale radii for the atoms in the dictionary, i.e., we assume c h = 4. Note that each frame of our dictionary fixes the scale radius in comoving coordinates and thus the NFW atoms have different angular sizes when placed at different redshift.
We test the algorithm with varying the penalty parameter for the LASSO estimation with λ = 3.5, 4.0, and 5.0. The corresponding penalty parameters for the final adaptive LASSO estimations are set to λ ad = λ τ +1 . Here, both the LASSO estimation and our adaptive LASSO estimation select the pixels with signal-to-noise ratios greater than λ in each gradient descent iteration, and the local density is estimated for the selected pixels with a shrinkage of the estimation amplitude. The LASSO estimation shrinks the density amplitudes by λ for every selected pixels, whereas the adaptive LASSO estimation suppresses the shrinkage if the preliminary estimation for the pixel is greater than λ, by downweighting their penalties, and otherwise it enhances the shrinkage.
We'd like to reconstruct with the resolution limit set by the Gaussian smoothing kernel with a standard deviation of 1. 5 and the 1 pixel scale as described in Sections 2.3.2. In addition, we smooth the reconstructed density field with the same Gaussian kernel in each lens redshift plane. Figure 9 shows the 3D density maps reconstructed with different penalty parameters for a halo with M 200 = 10 15.02 h −1 M at redshift 0.164. The corresponding pixel value histograms are shown in Figure 9. Our adaptive LASSO algorithm assigns zero to a fraction of the reconstructed pixels that are mostly noises, while retaining strong signals. The right panels demonstrate a case where almost all of the pixels with pure noise are set to zero when a large penalty parameteris adopted. It is important to note that the reconstructed density maps are not compromised by line of sight smearing, i.e. the reconstructed lens is localized in redshift. Following Leonard et al. (2014), we normalize the detected peaks in the lth (l = 1...20) lens redshift plane to account for the peak amplitude difference arising from the difference in the norm of the lensing kernels in different redshift bins: where the normalization matrix is defined as where K(z l , z s ) is the lensing kernel. In Figure 10, we show the histograms of the normalized peaks with different penalty parameters. There, we stack the histograms from 100 realizations of all the halos sampled in the redshift-mass plane. In addition, we generate 1000 realizations of pure noise catalogs and perform the reconstruction using the noise catalogs in order to examine the noise properties. The histograms of the normalized peaks detected from the pure noise catalogs are shown in Figure 10 along with the best-fit Gaussian functions of the noise peak histograms. We find that the number counts including both true and false peaks decrease as the penalty parameter λ increases. In addition, the standard deviation of noise peaks decreases as λ increases. As a result, with λ = 5.0, we find a clearer excess in the positive peak counts compared with the noise peak histograms, especially at the high density contrast. This is expected because a higher penalty parameter prefers a sparser solution; more peaks originating from the noise are removed than those from real clusters at the high density contrast.
We demonstrate the position offsets of the detected halos compared to the input positions in Figure 11. The left panel shows the 2D number histograms stacked over all the simulations as a function of the angular distance and the redshift offsets. We see clear clustering of the peaks close to the position of the input halo on the stacked number histogram. For each simulation, we identify positive peaks closest to the input position (in the pixel unit). If a closest peak is located inside the region denoted with the dashed box in the left panel of Figure 11, we regard it as a "true" peak detection.
Other identified peaks, which include both positive and negative peaks, are judged to be false detection.
The right panel of Figure 11 shows the estimated redshift of the "true" detections averaged over the simulations (with different noise realizations) for each halo as a function of the input redshift. As shown, the estimated redshifts are slightly lower than the ground truth by ∆z ∼ 0.03 for halos in the low-redshift range (z ≤ 0.4). For halos at 0.4 < z ≤ 0.85, the relative redshift bias is below 0.5%. The standard deviation of the redshift estimation is 0.092.
In order to reduce the number of false detections, we select peaks as galaxy clusters if the peak values are greater than an empirically determined threshold. The threshold is set in units of the standard deviation of the noise peaks. We use different detection thresholds (1.5σ and 3.0σ) to detect galaxy clusters from the mass maps reconstructed with different penalty parameters. Figures 12 and 13 show the detection rates for halos in the redshift-mass plane for different penalty parameters (λ = 3.5 and λ = 5.0). The corresponding number of false detections per square degree as a function of detection threshold are also demonstrated in the figures.
As shown, the false peak density is successfully reduced for relatively large detection threshold, but the detection rate of halo also decreases. After a few experiments, we have decided to set the detection threshold to 1.5σ and set the penalty parameter λ to 5.0 since the combination suppresses the false detection to 0.022 while keeping a high halo detection rate. In summary, our method is able to detect halos with minimal mass of 10 14.0 h −1 M , 10 14.7 h −1 M , 10 15.0 h −1 M for the low (z < 0.3), median (0.3 ≤ z < 0.6) and high (0.6 ≤ z < 0.85) redshift ranges, respectively.
Using the detection rate measured from our simulations, we are able to predict the number density of detected clusters by assuming the halo mass function of Tinker et al. (2008). We use HMF (Murray et al. 2013), an open-source package 4 , to calculate the halo mass function. The predicted halo detection number density for the setup λ = 5 and 1.5σ detection threshold is shown in Figure 14. The resulting cluster number density is 0.49 deg −2 , which is much higher than the false detection rate of 0.022 deg −2 . This cluster number density corresponds to 78.4 clusters for the first-year HSC shear catalog ) with a survey area of ∼ 160 deg 2 . The expected number of detection is slightly higher than the number of 2D cluster detections (63 detected clusters) for the first year HSC shear catalog (Miyazaki et al. 2018b). Furthermore, our 3D detection method provides an accurate redshift estimation for individual clusters. In contrast, the redshift information is not provided from the 2D mass map reconstruction.

Point mass atoms
We perform an additional test by substituting the default NFW dictionary with point-mass. This test may indicate some certain limitation of our method when applied to a case with very compact (point) objects, which however is unlikely to happen in actual lensing observations.
We set the penalty parameter for the preliminary LASSO to λ = 3.5 and 5.0. Pramanik & Zhang (2020) Figure 11. The left panel shows the stacked 2D distribution of the deviations of detected peak positions from the centers of the corresponding input halos. The x-axis is for the deviated distance in the transverse plane, and the y-axis is for the deviation in redshift. In each simulation, the positive peak inside the dashed black box with the minimal offset (in the pixel unit) from the position of the input halo is taken as "true" detection. The right panel shows the redshift deviation of detected peaks. The x-axis is the input halo redshifts, and the y-axis is the redshift of the detected peak. The cross-points denote the average of the detected peaks for each halo over different noise realizations, and the error-bars indicate the uncertainties of the averages. The deep gray area indicates relative redshift bias less than 0.05, and the light gray area for relative redshift bias less than 0.5. These results are based on our reconstruction with the NFW dictionary with λ = 3.5 . propose to incorporate group information into different adaptive LASSO penalty weights by setting the weights for projection coefficients in a group to the average of the adaptive weights in the same group. We assume that the neighboring pixels in the same redshift plane belong to the same structural group (e.g., galaxy cluster and void), and smooth the amplitude of the preliminary LASSO estimation in each lens redshift plane with a top-hat filter of comoving diameter r c = 0.25 h −1 Mpc. Let us denote the amplitude of the preliminary LASSO estimation with the smoothing as x LASSO sm , and adopt the penalty weights given byŵ = 1/ x LASSO τ sm . Figure 15 shows the reconstruction result with the point-mass dictionary. Interestingly, several "discrete" masses are assigned to at different redshift bins in the neighboring region of the true halo center. In contrast, The expected number density of detected clusters per square degree as a function of halo redshift (xaxis) and halo mass (y-axis). The number density in total is 0.49 deg −2 .
as we have seen in Figure 9, the NFW dictionary manages to recover a consistent mass distribution. The problem of the point-mass dictionary originates from the fact that the profile of the point-mass atom in the transverse plane is much more compact than the profile of the input halo, especially when placed at low redshifts.

SUMMARY
We have developed a novel method to generate highresolution 3D density maps from weak-lensing shear measurement with photometric redshift information. A key improvement over previous similar methods is that we represent a 3D density field by a collection of NFW atoms with different physical sizes. With a prior assumption that the clumpy mass distribution is sparse in 3D, we reconstruct the density map using the adaptive LASSO algorithm (Zou 2006). We show that adopting the standard LASSO algorithm results in significant smearing of structure in the line of sight direction even in the absence of galaxy shape noise and photometric redshift uncertainties. Our adaptive LASSO algorithm efficiently reduces the smearing of structure.
We have examined the performance of cluster detection with the reconstructed 3D mass maps using mock catalogs that apply shear distortions from isolated halos to galaxies with HSC-like shapes and photo-z uncertainties. Under the realistic conditions, our method is able to detect halo with minimal mass limits of 10 14.0 h −1 M , 10 14.7 h −1 M , 10 15.0 h −1 M at low (z < 0.3), median (0.3 ≤ z < 0.6) and high (0.6 ≤ z < 0.85) redshifts, respectively, with an average false detection of 0.022 deg −2 . The estimated redshifts of the clusters detected in the reconstructed mass maps are slightly lower than the true redshift by ∆z ∼ 0.03 for halos at low redshifts (z ≤ 0.4). The relative redshift bias is below 0.5% for halos at 0.4 < z ≤ 0.85, and the standard deviation of the redshift estimation is 0.092.
We will apply the newly developed 3D mass map reconstruction technique to the wide-field data from HSC survey (e.g., Mandelbaum et al. 2018;Li et al. 2020) and perform galaxy cluster detection in our future work.  Figure 15. The density maps reconstructed from the mock galaxy shape catalog with the point-mass dictionary. The penalty parameter is λ = 3.5 (left) and λ = 5.0 (right). The axes and color maps have the same meanings as Figure 7. The input halo mass is M200 = 10 15.02 h −1 M , and its redshift is z = 0.164.