Weak-lensing Mass Reconstruction of Galaxy Clusters with a Convolutional Neural Network - II: Application to Next-Generation Wide-Field Surveys

In this section, we provide a brief overview of the WL theory. For more details on gravitational lensing, readers can refer to review papers [8, 25, e.g.,].

The distortion of source galaxies in the WL regime can be expressed as $\mathbf{x^{\prime}}=\mathbf{Ax}$, where $\mathbf{x}$ ($\mathbf{x^{\prime}}$) is the position in the source (image) plane. The matrix $\mathbf{A}$ is expressed as:

$$\mathbf{A}=(1-\kappa)\begin{pmatrix}1-g_{1}&-g_{2}\\ -g_{2}&1+g_{1}\end{pmatrix}.$$

(1)

In Equation 1, $\kappa$ denotes the convergence, and $g_{1(2)}$ is the first (second) component of the reduced shear $g=({g_{1}}^{2}+{g_{2}}^{2})^{1/2}$. The reduced shear $g$ is computed from the relation $g=\gamma/(1-\kappa)$, where $\gamma$ is the shear. The convergence $\kappa$ is defined as:

$$\kappa=\frac{\Sigma}{\Sigma_{c}},$$

(2)

where $\Sigma$ is the surface mass density, and the critical surface mass density $\Sigma_{c}$ is given by:

$$\Sigma_{c}=\frac{c^{2}D_{s}}{4\pi GD_{d}D_{ds}}.$$

(3)

Here, $c$ is the speed of light, $D_{d(s)}$ is the angular diameter distance to the lens (source), and $D_{ds}$ is the angular diameter distance between the lens and the source. $\gamma$ and $\kappa$ are related by the following convolution:

$$\bm{\gamma}\mathbf{(x)}=\frac{1}{\pi}\int\mathbf{D}(\mathbf{x}-\mathbf{x^{\prime}})\kappa(\mathbf{x^{\prime}})d\mathbf{x^{\prime}},$$

(4)

where the convolution kernel $\mathbf{D}$ is:

$$\mathbf{D}\mathbf{(x)}=-\frac{1}{(x_{1}-\mathbf{i}x_{2})^{2}}.$$

(5)

We follow the procedure outlined in H21 for generating the training data. We use the publicly available convergence map from the MassiveNuS cosmological simulation [44]¹¹ 1 http://columbialensing.org. From the convergence map datasets at five different source redshifts ($z=0.5,1.0,1.5,2.0,~\rm and~2.5$), we select the dataset created for the source plane at $z=2.5$. Note that this $z=2.5$ source redshift of the convergence map is different from the source redshift that we assume for creating the WL galaxy shape catalog described below. In H21, we selected the convergence map for the source plane at z=1.5; however, we find that using the dataset for the source plane at $z=2.5$ improves the reconstruction quality. Due to the higher lensing efficiency, the convergence map at $z=2.5$ contains richer structures. One potential concern is the risk of bias, as the source plane at $z=2.5$ is significantly higher than the typical value ($z\sim 1$) achieved with ground-based data. However, our CNN model is designed to learn the relationship between the shear catalog and the corresponding underlying mass distribution, rather than the characteristic convergence pattern at a specific source redshift. Consequently, we find that this particular choice of source redshift for our training dataset does not introduce a significant bias in the mass reconstruction. The dataset consists of 10,000 convergence maps covering a field size of $3\fdg 5\times 3\fdg 5$ with a resolution of $512\times 512$ ($\sim 0\farcm 410$ per grid cell).

When the nominal 10-year mission is completed, the LSST is expected to cover the entire southern sky with a depth of 27 mag. From the Subaru Hyper-Suprime Cam (HSC) archival images with the expected median seeing ($\raise 0.81805pt\hbox{$\scriptstyle\sim$}0.7\arcsec$) of the LSST, we find that the average source number density is $\scriptstyle\sim$33 ${\rm arcmin}^{-2}$. Since the Subaru/HSC imaging data share many properties with the Rubin Observatory, we adopt this number density as the average source density for our mock data [30]. The noise in the shape of each source galaxy is comprised of the shape dispersion and measurement components. As in H21, we adopt $\sigma_{e}=0.24$ as the intrinsic shape dispersion. For the measurement error, we utilize the relation between galaxy magnitude and measurement error determined from an analysis of Subaru/HSC data in HyeongHan et al. [30].

We divide the 10000 convergence maps into 7000 for training, 2800 for validation, and 200 for testing. Data augmentation is employed to increase the number of training data by applying four rotations and two-axis flips, but we do not use random crops as done in H21. The resulting number of convergence maps for training is $7000\times(4+2)=42000$. Within the 3$\fdg$5 x 3$\fdg$5 field, there are no preferred locations for galaxy clusters unlike in H21.

Table 1 presents the architecture of our new CNN model. The input layer is the $3\times 512\times 512$ array consisting of the three ($\epsilon_{1},\epsilon_{2}$, and $\Delta\epsilon$) channels. The first two channels are the average ellipticity components within the grid cell, while the last channel is the average ellipticity measurement error. Figure 1 presents an example of the three-channel input layers. Unlike in H21, we did not smooth the input layers with a Gaussian kernel²² 2 H21 needed to smooth the input layers because the source density per pixel was much lower than in the current setup. . To mitigate overfitting associated with empty pixels, we applied smoothing to the input layers. However, in this study, we did not smooth the input layers because most pixels contain sufficient ellipticity information.

Our CNN model is composed of 8 two-dimensional (2D) convolutional layers, whereas the previous model used 11 2D convolutional layers. Instead of increasing the depth of our CNN model, we improved it by employing different filter sizes in the initial layers and concatenating them. This approach allows our CNN model to capture both small and large-scale structures with shallower layers than H21, making it more compact and efficient. Additionally, unlike in H21, we maintain the input layer size ($512\times 512$) throughout the CNN operations to minimize resolution loss. The number of free parameters in our CNN model is $5.6\times 10^{6}$.

Similarly to the approach in H21, a skip connection is added after each combination of the convolution and transposed convolution layers. This measure is motivated by the residual neural network [24, ResNet;], which is shown to improve the overall performance of the CNN model by reducing the problem of vanishing gradients. The activation function for the convolutional layer is the Leaky Rectified Linear Activation (LReLU) function, which has a slope of 0.2 when $x<0$. To prevent overfitting, we applied Dropout layers [56] after the Concatenate and Multiply layers, randomly excluding neurons during training. In this study, we excluded 25 percent of neurons per convolution layer. For training our CNN model, the Adam optimizer [40] was used with a learning rate of $10^{-5}$. We ran 1000 epochs and adopted gradient clipping for the optimizer [62] to prevent exploding gradients.

In H21, we utilized a weighted mean-squared-error loss function inspired by the focal loss [43]. This approach emphasized the importance of higher $\kappa$ values. The loss function in H21 is given by:

$$L=\sum_{\bm{x}}\omega_{f}(\bm{x})[\kappa_{pred}(\bm{x})-\kappa_{truth}(\bm{x})]^{2},$$

(6)

where $\omega_{f}(\bm{x})$ is a weight function defined as:

$$\omega_{f}(\bm{x})=1+\frac{|\kappa_{truth}|}{{\rm max}(\kappa_{truth})}.$$

(7)

In this study, we aim to reconstruct both high-contrast individual clusters and more diffuse large-scale structures. The ratio of high ($\gtrsim$0.2) to low ($\lesssim$0.1) $\kappa$ pixels in our WF dataset is lower than that in H21, rendering the original weighting scheme less effective for our purposes. Thus, we have modified the loss function to better suit the broader range of structures as follows:

$$L=\sum_{\bm{x}}{\omega_{f}(\bm{x})}^{6}[\kappa_{pred}(\bm{x})-\kappa_{truth}(\bm{x})]^{2},$$

(8)

where $\omega_{f}(\bm{x})$ is the weight defined in Equation 7.

This revision amplifies the influence of high-$\kappa$ pixels more than the original weighting scheme. The choice of the exponent 6 is empirically determined to provide optimal balance and training stability, enhancing the model’s sensitivity to high-$\kappa$ pixels while preventing overfitting. By increasing the weight of high-$\kappa$ pixels, the model is better tailored to accurately predict these pixels while still maintaining its ability to reconstruct low-$\kappa$ pixels without overfitting.

Figure 2 compares the mass reconstruction results from the current CNN, H21, and KS93 methods for five randomly selected fields from the test dataset. The $\kappa$ values from our CNN model range from $\raise 0.81805pt\hbox{$\scriptstyle\sim$}-0.05$ to $\raise 0.81805pt\hbox{$\scriptstyle\sim$}0.3$, consistent with the true convergence maps. Notably, the reconstructed maps from our CNN model exhibit good agreement in features with the true $\kappa$ values for both low ($\lesssim$0.1) and high ($\gtrsim$0.2) $\kappa$ values, marking an improvement over the results from the H21 model. In contrast, the KS93 reconstructions tend to significantly overestimate the $\kappa$ values mainly due to the mass-sheet degeneracy. We also found that the residuals between the truth and reconstructed mass maps from our CNN model do not present any significant systematic pattern.

The model from H21 tends to smooth out small-scale structures, a trend discernible across all samples in Figure 2. We attribute this to the reduction of the input channels in the first convolutional layer. While this approach was introduced to enhance training efficiency, the reduction makes it difficult to capture small-scale structures. The KS93 results show correlations with the large-scale structures of the true mass distribution but significantly smooth out finer details. Of course, one can consider decreasing the kernel size for KS93. However, when we decreased the kernel size to restore the compactness of the highest peaks, the result showed numerous spurious fluctuations across the field. This is not surprising because the reduced kernel size is inadequate in low density regions.

Finally, the superiority of our new CNN model can also be evaluated by the loss $L$ (see the value in the lower right corner of each panel in Figure 2) evaluated with Equation 8. The $L$ value from our new CNN results spans the range $200-210$, which is $10-20\%$ lower than those from the H21 results. The KS93 result shows significantly higher $L$ values.

Following the qualitative visual analysis in §III.1, we present a quantitative evaluation by directly comparing convergence values. Figure 3 displays pixel-by-pixel correlations between the truth and reconstructed $\kappa$ maps. We find that the current CNN results exhibit a stronger correlation, with slopes of $\raise 0.81805pt\hbox{$\scriptstyle\sim$}0.47$ and $\raise 0.81805pt\hbox{$\scriptstyle\sim$}0.70$ for the H21 and current models, respectively. In contrast, the slope deviates substantially from unity for the KS93 mass reconstruction.

Figure 4 compares the probability distribution function (PDF) of $\kappa$ values. Once again, it is evident that the distribution derived from the current CNN model matches the truth more closely than the H21 model. In particular, the agreement is significantly improved in the recovery of high $\kappa$ values ($\kappa>0.1$). The distribution is considerably wider in the KS93 reconstruction.

In §III.1 and §III.2, we evaluated the fidelity of our CNN mass reconstruction by comparing 2D mass maps, pixel-by-pixel values, and $\kappa$ distributions. Here, we assess the similarity in structures between the reconstructed map and truth by comparing their $\kappa-\kappa$ auto-correlation functions.

To calculate the auto-correlation function, we used TreeCorr ³³ 3 https://github.com/rmjarvis/TreeCorr [34, 33] on 200 mass maps from the truth and the reconstruction. Figure 5 displays the results. We find that the correlation is well recovered at both small and large scales. The discrepancy at $\lesssim 2$ pixels ($\lesssim 0.8$ arcmin) is due to the inevitable resolution loss in the reconstruction caused by the finite source density.

One of the key merits of future WL observations is their ability to detect galaxy clusters solely based on projected mass maps. This method offers significant advantages over traditional approaches such as galaxy overdensity mapping, Sunyaev-Zeldovich effect observations, and X-ray surveys because the signals directly correlate with their masses.

We assess the cluster detection capability of our CNN model by the following. First, we identify cluster candidates as local peaks on both the truth and reconstructed maps. We then compute projected masses within a 1 Mpc radius for each identified cluster and evaluate the detection performance by comparing the projected mass and completeness. Clusters from the truth and CNN results were matched based on a distance criterion, which we set to 6 pixels ($\sim 2\farcm 4$).

Figure 6 compares the truth and CNN reconstructed projected cluster masses. The slope from the linear regression (dashed magenta line) is very close to the one-to-one correlation (solid black line), indicating that the convergence values obtained directly from the CNN maps serve as accurate mass proxies.

In Figure 7, we present the completeness of cluster detection, defined as the ratio of the number of detections in the reconstructed mass maps to those in the true convergence maps, as a function of the cluster mass. All three methods achieve $100\%$ completeness for clusters with masses larger than $\raise 0.81805pt\hbox{$\scriptstyle\sim$}3.5\times 10^{14}M_{\odot}$. Our CNN model maintains a completeness of $\raise 0.81805pt\hbox{$\scriptstyle\sim$}75\%$ down to $\raise 0.81805pt\hbox{$\scriptstyle\sim$}10^{14}M_{\odot}$. The KS93 completeness is much lower, yielding $\raise 0.81805pt\hbox{$\scriptstyle\sim$}35\%$ at $\raise 0.81805pt\hbox{$\scriptstyle\sim$}10^{14}M_{\odot}$. The completeness measured from the H21 result is $\raise 0.81805pt\hbox{$\scriptstyle\sim$}10$% lower than that of the current version in the $\lesssim 2.5\times 10^{14}M_{\odot}$ regime.

In the upcoming era of the Vera C. Rubin Observatory, the Roman Space Telescope, and the Euclid mission, WF WL mass reconstructions will be performed routinely. By leveraging the unique property of WL, which eliminates the need for dynamical assumptions (i.e., it is only influenced by the gravitational potential), our CNN model would be an effective and independent tool for detecting galaxy cluster candidates. This capability enables a range of studies, such as the cluster mass function reconstruction [47, 18, 11] and cosmology [61, 59, 28]. In addition, the independently obtained galaxy cluster candidate catalogs will synergize with other cluster surveys across various wavelengths [10, 41, 15].

In the evaluation of machine learning (ML) models, a null test serves as a critical validation step to ensure that the model does not produce false positive signals from random noise fields. Since H21 placed clusters always near the center of the field in the training stage, a moderate level of bias was observed for null input data.

In order to investigate the same issue for the current model, we generated 1000 mock WL observations from null fields (i.e., $\kappa=0$) with the same observational properties (i.e., including measurement and shape noise) as the training data and reconstructed the corresponding mass maps. Figure 8 shows the $\kappa$ distributions from these mass maps. The $\kappa$ values range between -0.04 and 0.06, and the distribution is positively skewed, with a mode around $-0.01$. However, compared to H21, the distribution is relatively symmetric, and the bias level is very low. This improvement is somewhat expected, as the current training dataset is produced from random cosmological ray tracing fields. Our null test provides confidence in the robustness of the model and its ability to produce reliable mass reconstructions under realistic observational conditions.

Traditional WL mass reconstructions are affected by boundary and finite-field effects due to the lack of information beyond the mass reconstruction field. In §III.1, our CNN results do not exhibit these artifacts for square boundaries.

Real observations often have irregular observational footprints, which can introduce biases in mass reconstruction. To investigate these effects under more realistic conditions, we reconstruct $\kappa$ maps with masked input data. We adopt the observational footprint of the Subaru/HSC survey of the Coma cluster presented in HyeongHan et al. [30] as our mask. Figure 9 shows the resulting mass reconstruction.

While all three methods can restore the overall mass distributions, we find that the current CNN results are the most satisfactory. Again, the KS93 maps show highly inaccurate convergence values. The H21 result is significantly better that the KS93 result; However, its convergence values are systematically lower than the true values. This experiment demonstrates that the current CNN model can be applied to real observations, where complex boundary shapes are present.

While ML has shown remarkable promise in various fields, it is not without its pitfalls. A significant challenge arises when models that perform exceptionally well on training data fail to generalize to real-world applications. ML models are often sensitive to the quality and representativeness of the training data. If the training set does not adequately reflect the diversity and complexity of real-world scenarios, the model may become biased or overly simplistic. This lack of robustness can result in poor performance in practical applications, where conditions may vary significantly from those encountered during training. Additionally, changes in the underlying data distribution, known as dataset shift, can further exacerbate these issues, rendering a once-effective model ineffective.

To ensure that the current CNN model excels not only on training datasets but also on real-world data, we experimented with actual WF WL observations. We present the mass reconstruction of the Coma cluster in the left panel of Figure 10. We used the WL shear catalog from [30], which measured WL signals from Subaru/Hyper-Suprime Cam imaging data covering approximately $3\fdg 5\times 3\fdg 5$ centered on its two BCGs. The Coma cluster, being a very massive cluster, may serve as an excellent real-world test case, as its characteristics are expected to differ significantly from those of the average cluster in our training dataset.

As the true mass distribution in the Coma field is unknown, we cannot assess its fidelity using the same metrics presented in §III. Instead, we provide comparisons with previous results as follows. First, we find that the mass estimate derived from our convergence field is consistent with previous studies. Fitting an NFW profile provides $M_{200}=7.4_{-3.1}^{2.4}\times 10^{14}M_{\sun}$, which aligns with literature values [21, 48, e.g.,]. Second, the reconstruction reveals a number of significant substructures, some of which were also reported by previous low-resolution WL studies. Third, the mass distribution is well traced by those of the Coma cluster galaxies and intracluster light [48, 36, e.g.,].

Most notably, the CNN mass map reveals some of its intracluster filaments reported by [30]. Intracluster filaments are the terminal ends of intercluster filaments that stretch over tens of megaparsecs. Due to their low contrast, their detection has been considered elusive. [30] identified the Coma cluster’s intracluster filaments not through mass reconstruction, but using two novel methods: the matched-filter technique and the shear peak statistic [45]. Here, we independently detect the intracluster filament of the Coma based on the current CNN mass map as follows. First, we ran DisPerSE⁴⁴ 4 https://www2.iap.fr/users/sousbie/web/html/indexd41d.html [53, 54] on the CNN mass map. To set the threshold cut for filament detection, we adopt a persistence above the 3-sigma level, where persistence indicates the absolute difference between the saddle and maximum critical points. The solid pink lines in the right panel of Figure 10 depict the resulting filaments. The three main branches originating from the cluster center is clearly detected and aligns well with the result from [30]. Second, we measured the azimuthal excess of the convergence with respect to the background value. Figure 11 shows that this method too reveals the three main filaments identified by DisPerSE, which is also consistent with the findings of [30] (compare with Figure 2 of their paper). The consistency in intracluster filament detection showcases the utility of CNN-based WL mass reconstruction for identifying large-scale structures in future WL surveys, including the detection of cosmic webs. We refer readers to HyeongHan et al. [30] for further scientific details on the Coma cluster.