Three-Dimensional Reconstruction of Weak-Lensing Mass Maps
with a Sparsity Prior. II. Weighing Triaxial Cluster Halos

Following the definition of Jing & Suto (2002), we adopt a halo density profile of

where

Here, $\delta_{\mathrm{ce}}$ is the concentration parameter of the halo. $\rho_{\mathrm{crit}}(z)$ is the critical density of the universe at redshift $z$. We denote the scale radius of the triaxial halo as $R_{0}$, and $a,b$ and $c$ are scaling factors that describes the shape of the halo. We set $b=c=1$ in the rest of the analysis in this paper, for simplicity and normalization.

From equation (8), we see that the case of an isotropic halo model with $\alpha=1$ reproduces the NFW halo profile. Various literature uphold different values of $\alpha$ ranging between $\alpha=1$ and $\alpha=1.5$ (see e.g., Navarro et al. 1997, Moore et al. 1999, Oguri et al. 2003). Therefore, we perform our subsequent analyses both for halos with $\alpha=1$ and $\alpha=1.5$.

Jing & Suto (2002) also define a length scale $R_{e}$ such that $R_{e}/r_{\text{vir}}=0.45$  ($r_{\text{vir}}$ is the virial radius) and $\frac{R_{e}}{R_{0}}=c_{e}$ is the concentration parameter. The average density within an ellipsoid of $R_{e}$ is

where, from Oguri et al. (2001a), we have

Assuming that $\Omega_{m}+\Omega_{\Lambda}=1$ , $\omega_{\text{vir}}=1/\Omega_{\text{vir}}-1$ , and the density parameter $\Omega_{\text{vir}}$ at the redshift of virialization $z_{\text{vir}}$ is (Oguri et al., 2001b)

In this paper, we take the empirical relation between the concentration parameter ($c\equiv\frac{r_{\mathrm{vir}}}{R_{0}}$), the viral mass ($M_{\mathrm{vir}}$) and redshift ($z$) of the halo:

where $A=6.02$, $B=-0.12$, $C=0.16$, and $M_{\odot}$ is the solar mass (Child et al., 2018).

To describe rotations of a traxial halo, we introduce two coordinate systems: The $(x,y,z)$ system is the dark matter halo’s system, with its origin at the halo’s center and $z$ axis lies along the major principle axis of the halo. The $(x^{\prime},y^{\prime},z^{\prime})$ coordinate represent the observer’s coordinate system, with its origin also set at the halo’s center. We define $(\theta,\phi)$ to be the polar coordinates of the line-of-sight direction of the observer with the halo’s long axis as the $z$-axis. Just like Oguri et al. (2003), we note that the $x^{\prime}$-axis lies in the $x$-$y$ plane. We also define $\zeta$ such that

where the primed coordinate $(x^{\prime},y^{\prime})$ represent normalized observer’s coordinate. Fig. 3 is an illustrative plot that shows the relationship between the coordinates. Then we get an expression for the lensing convergence $\kappa$ as:

where

with $f=\sin^{2}\theta(\frac{c^{2}}{a^{2}}\cos^{2}\phi+\frac{c^{2}}{b^{2}}\sin^{2}\phi)+\cos^{2}\theta$, and

The subscripts, “TNFW”, corresponds to “Triaxial NFW”. $\Sigma_{\text{crit}}$ is the lensing critical surface mass density defined as

where $D_{\mathrm{OL}},D_{\mathrm{OS}}$, and $D_{\mathrm{LS}}$ are the angular diameter distances from the observer to the lens plane, from the observer to the source plane, and from the lens plane to the source plane, respectively.

Once we have an expression for $\kappa$, we may follow Keeton (2001) to calculate the shear field. Note that although in $\alpha=1$ case an analytical solution can be yield, an analytical solution does not exist for a density profile with $\alpha=1.5$ . Differing from equation (4) in Li et al. (2021), we did not adopt truncation at the viral radius to facility numerical computation of shear fields.

We use halos with different triaxial shapes to produce shear fields for simulation. This is represented with a wide range of ellipticity, defined as $\frac{a}{c}$ (see equation 9). However, to reduce dimensionality and computational time during the mass map reconstruction, the dictionaries are prepared with isotropic halos with $a\kern-2.79999pt=\kern-2.79999ptb\kern-2.79999pt=\kern-2.79999ptc\kern-2.79999pt=\kern-2.79999pt1$ described in the previous section. The shear field are measured at redshift of $z\kern-2.79999pt=\kern-2.79999pt0.05,0.36,0.47,0.56,0.69,0.80,0.91,1.03,1.22,1.50,$ and $2.50$.

To consider variation of halo profiles as found in recent high-resolution $N$-body simulations (Navarro et al., 2010), we include the ability to simulate shear field produced by triaxial halos (Jing & Suto, 2002) with both (i) the NFW radial profile with $\alpha=1$ in equation (8) (Navarro et al., 1997) and (ii) the cuspy NFW radial profile with $\alpha=1.5$ (Jing & Suto, 2000).

We account for statistical uncertainties in shape estimation from galaxy intrinsic shape noise and measurement error due to image noise, calculated using the first-year shear catalog of the HSC first-year data (Mandelbaum et al., 2017). Specifically, we utilized the formulation of Shirasaki et al. (2019), where we have

In the above expression $\epsilon^{\mathrm{int}}$ represents the per-component intrinsic shape error, and $\epsilon^{\mathrm{mea}}$ represents the per-component shape measurement error. Also, $e^{\mathrm{ran}}=e^{\mathrm{obs}}e^{i\phi}$, where $\epsilon^{\mathrm{obs}}$ is the distortion of some individual galaxy and $e^{i\phi}$ serves to rotate the observed shape by some random angle $\phi$. $N_{1}$ and $N_{2}$ are random numbers drawn from a Gaussian centered at $0$ with a standard deviation of $\sigma_{\mathrm{e}}$. $\epsilon_{\mathrm{rms}}$ is the root-mean-square of the intrinsic galaxy shape for each shape component. $\sigma_{e}$ is the standard deviation of the shape measurement error due to image noise for each shape component. Note, $\epsilon_{\mathrm{rms}}$ and $\sigma_{\mathrm{e}}$ are estimated from image simulations at single galaxy level using realistic galaxy image simulations (Mandelbaum et al., 2018a).

We assume that the multiplicative and additive biases in the shear catalog is fully corrected in this paper; therefore, we have an expression for the observed shear:

where $\mathcal{R}=1-\left\langle\epsilon_{\mathrm{rms}}^{2}\right\rangle$. We can then substitute

to get a mock shear field.

For realistic noisy tests, we adopt realistic HSC-like galaxy number density ($\sim$20 arcmin^-2) (Mandelbaum et al., 2018b; Li et al., 2022) when producing shear fields in order to test the performance of our algorithm with noisy setup. However, for the noiseless tests in section 3.3, we adopt an extremely high galaxy number density (2000 arcmin^-2) to suppress the random noise in the sub-pixel galaxy distribution.

The performance of the our mass mapping algorithm may depend on the regularization parameter $\lambda$ for the noisy case. We present the how SPLINV behaves differently with different $\lambda$ values in this section.

To determine the sparsity parameter $\lambda$ in equation (5) that optimizes reconstruction results, we perform various single halo reconstruction of shear field produced by an isotropic halo of mass $10^{14.6},10^{14.8}$, and $10^{15.0}$ respectively, at an intermediate redshift ($z=0.2425$) in the center of in a $48\times 48$ pixelized grid covering $98\,\text{arcmin}\times 98\,\text{arcmin}$ of sky area. Examining the density plots (Figs. 4—6) for different values of $\lambda$, we find that while “the best” $\lambda$ parameter potentially exists for each case, a smaller $\lambda$ value tend to make SPLINV to provide a smaller mass estimation than those provided by SPLINV with a larger $\lambda$, most likely due to a smaller $\lambda$ relaxes the sparsity condition as can be seen in equation  (5).

We conclude that a relatively higher $\lambda$ should more strongly enforce the sparsity condition, while effective making a cutoff for false detections with small masses. For this same reason, larger $\lambda$ reduce the probability of detecting halos with small masses.

We find that the optimal value of $\lambda$ depends on both the mass and the redshift of the halo, as reconstruction of a halo with larger SNR (higher mass and lower redshift) prefers a larger $\lambda$. From this, we conclude that we should find an optimized $\lambda$ for interval targeted detection mass/redshift, and then recursively apply SPLINV to detect galaxy halo in each mass/redshift interval. For detecting halos with relatively smaller masses ($~{}10^{14.6}~{}\textup{M}_{\odot}$), authors recommend using $\lambda=2$ and for halos with large masses $~{}10^{15.2}\textup{M}_{\odot}$, we recommend using $\lambda=4$. We leave the study on the optimal setup of $\lambda$ using realistic ray-tracing simulations (Takahashi et al., 2017) to future works. More specifically, what our findings can be concluded as:

1. (i)

   For simulations with small halo masses ($10^{14.6}M_{\odot}\leq M\leq 10^{14.8}\,M_{\odot}$), we find a $\gtrapprox 18\%$ positive mass estimation bias exists for reconstruction with $\lambda>=2.5$. This is possibly due to a form of Eddington Bias (e.g., see Kelly (2007) and Eddington (1913)), where only halo’s shear signal boosted by noises gets detected which are then confused with halo with a large mass.

2. (ii)

   For halos with larger mass ($M\geq 10^{15.0}M_{\odot}$) we do not find significant mass bias with $\lambda\geq 3$ due to the high SNR of these halos.

3. (iii)

   For halos with larger mass, SPLINV slightly underestimate halo masses with small $\lambda$ ($\lambda\leq 2.5$). This is possibly because the strong signal of higher mass halos may be cause the sparsity condition of SPLINV to fail and be construed by our algorithm as caused by multiple halos.

Because we found SPLINV with $\lambda=2$ performs well with smaller mass halos (which are more abundant in the universe) and only suffers overestimation slightly, we set $\lambda=2$ as fiducial setup and put results with $\lambda=4$ into Appendix (A).

In this section, we present reconstruction results for triaxial NFW ($\alpha=1$) and cuspy NFW ($\alpha=1.5$) halos.

We simulate $100$ halos with $10$ different shapes (from $\frac{a}{c}=1$ to $\frac{a}{c}=0.5$) and $10$ different redshifts (from $z=0.0625$ to $z=0.4675$), and for each halo, we generate $500$ realizations of observational noise as described in Section 3.2.2.

In Figs. 7 and 8, we show the estimated relative mass biases for the two types of halo with ellipticity ($\frac{a}{c}$) ranging from $1$ to $0.5$ reconstructed with dictionary generated numerically with the same mass and concentration parameter, but with $\frac{a}{c}=1$ (isotropic). At lower redshifts, corresponding to stronger lensing signal, the choice of $\lambda=2$ gives halo mass estimations with error less than $10\%$ or even better. Observing the first panel in Figs. 7 and 8, we again see the effect that, because the shear produced by the underlying halo was too small, even with $\lambda=2$, only the halos whose shear was boosted by noise gets picked up by our algorithm, resulting in an overestimation. However, the non-monotonous pattern in the second panel of Figs. 7 and 8 indicates that one could potentially optimize the value of $\lambda$ for each halo at each redshift. We also see that for reconstruction of more massive halos results in an underestimation of masses. This could be caused by the fact that a smaller $\lambda$ enforces a weaker sparsity condition and SPLINV may confuse the large signal due to the massive halo as signal generated by two separate halos. Another trend we find is that, the smaller cuspy NFW halos are generally more significantly affected by noise and hence will tend to have larger mass estimation bias.

Additionally, we find that there is a small $\frac{a}{c}$ dependence on the estimated mass bias. However, this dependence is not nearly as strong as that in the noiseless case, which tells us that we should focus on optimizing the value of $\lambda$ or other detection strategies before we try to include other parameters (like the triaxiality of halo models or its rotation) that complicate our dictionary space.

Here we present the results of redshift estimations from noisy reconstructions in Figs. 9 and 10. We note that the slight overestimation of redshift for halos with low redshifts in the figures is due to the discrete nature and the lower boundary of the redshift bins: there cannot be an underestimation for redshifts for these halos. Other than this, we observe that the amplitude of the relative redshift estimation bias is consistently below $5\%$, with no significant dependence on the shape ($\frac{a}{c}$ value) of the halo.

In the previous sections, we focus on the cases where halo model used for reconstruction is the same as those used to create the shear field. In this section, we study the potential model bias due to the systematic difference between halo models in the universe and those used in our dictionary. Following the previous sections, we are using isotropic models in our model dictionary. In this section, we study the mass and redshift estimation under the condition that the dictionary used for construction does not match the underlying halo in the simulation that produces the shear field.

In Fig. 11, we show the effect of systematic error due to the models used for mass map reconstruction being different. More specifically, we show the result of estimating mass of a cuspy halo with the assumption that the underlying mass field is composed of NFW halos. Comparing Fig. 11 with Fig. 5, we find that although we used the “wrong” dictionary in Fig. 11, the reconstructed result still resembles that in Fig. 5.

Next, we study whether the just using isotropic halo models in our dictionary affect our ability to reconstruct highly anisotropic halos. Comparing Fig. 12 with the left panel of Fig. 5, we see good agreement with reconstruction using isotropic halo model and we rotate a highly anisotropic halo with $\frac{a}{c}=0.5$ in the polar direction for $\theta=30^{\circ},60^{\circ},90^{\circ}$. A set of illustrative plots is shown in Fig. 13. The results from this section and the previous one indicate that, even when the true halo that constitutes the $\kappa$ map of the universe may be anisotropic, one may still recover the underlying mass map using isotropic models.

To test the performance of our algorithm in realistic multi-halo cases, we consider the following simulation set-up. We use the same parameters adopted in the previous sections but with a sky covering $256$ arcmin $\times 256$ arcmin area, corresponding to 128 pixels $\times$ 128 pixels. The center of the stamp is set to $(\text{ra},\text{dec})=(0^{\circ},0^{\circ})$ . The three cases are:

1. (i)

   The first halo has $M=10^{14.8}~{}M_{\odot}$ at $z=0.1975$ with $\text{ra}=4440^{\prime\prime}$ and $\text{dec}=5520^{\prime\prime}\,$, and the second halo has $M=10^{14.6}~{}{M}_{\odot}$ at $z=0.2875$ with $\text{ra}=-4440^{\prime\prime}$ and $\text{dec}=-4440^{\prime\prime}$ . The distance are chosen to be far enough so that, even in the noisy simulations, there is little chance that two of the halos are falsely detected as one.

2. (ii)

   The first halo has $M=10^{14.8}~{}\textup{M}_{\odot}$ at $z=0.2875$ with $\text{ra}=4440^{\prime\prime}$ and $\text{dec}=5520^{\prime\prime}$. The second halo has $M=10^{15.2}\textup{M}_{\odot}$ at $z=0.3775$ with $\text{ra}=-4440^{\prime\prime}$ and $\text{dec}=-4440^{\prime\prime}$. The third halo has $M=10^{14.6}~{}M_{\odot}$ at $z=0.1975$ with $\text{ra}=0^{\prime\prime}$ and $\text{dec}=0^{\prime\prime}$ . The distance are chosen to be far enough so that, even in the noisy simulations, there is little chance that two of the halos are falsely detected as one.

3. (iii)

   Same with (ii) but we first perform a reconstruction with dictionary containing one halo with $M=10^{15.2}~{}M_{\odot}$ (with a higher lambda) first, and then subtract the shear field produced by a realistic reconstruction result, containing information on estimated mass and redshift.

The halo models used for the simulations and the reconstructions are both isotropic in this section.

In Fig. 14, we present result of case (i) with dictionary composed of halos with $M=10^{14.6}\textup{M}_{\odot}$ and with $M=10^{14.8}$ with $\lambda=2$. While the estimation for the halo with $M=10^{14.6}\textup{M}_{\odot}$ is accurate, we see an overestimation of halo mass for the halo with $M=10^{14.8}\textup{M}_{\odot}$. This is probably due to the fact that we chose a generally small $\lambda$ for the halo with this mass and a SPLINV confuses the shear field produced by the $10^{14.8}\textup{M}_{\odot}$ halo with that produced by the $10^{14.6}\textup{M}_{\odot}$, but with a much larger mass to match the strength of shear field. The slight decrease in detection probability is probability due to increase in parameter space caused by one additional available choice of atom which causes the gradient descent algorithm harder to converge.

In Fig. 15 and Fig. 16, we show reconstruction results of the two smaller mass halo in case (ii) (a good detection on the more massive halo can always be achieved with a high value of $\lambda$). This is done in 2 ways: in the first method, we used $\lambda=2$ and with dictionary composed of halos with $M=10^{14.6}\textup{M}_{\odot}$ and with $M=10^{14.8}~{}M_{\odot}$ without modifying the shear field. In the second method we first perform some detections of the larger mass halo. Randomly select a set of mass and redshift estimation (we used $M\approx 10^{15.11}M_{\odot}$ and $z=0.2875$ to produce the result), and then proceed with reconstruction of the remaining two halos with $\lambda=2$ and with dictionary composed of halos with $M=10^{14.6}\textup{M}_{\odot}$ and with $M=10^{14.8}~{}M_{\odot}$.

We find that there is no significant difference in performance whether we subtract the shear field produced by the large halo or not. However, the underestimation in halo with $M=10^{14.8}~{}M_{\odot}$ and overestimation in halo with $M=10^{14.6}~{}M_{\odot}$ is still present. This result shows that, if we focus on detecting smaller mass halos, we may safely use dictionaries of those smaller mass halos without worrying about the shear field produced by large mass halos to interfere with our detection, keeping in mind that any anomalous large mass estimation may be caused by some large halo.