MARS: A New Maximum-entropy-regularized Strong Lensing Mass Reconstruction Method

MARS is based on a free-form source plane minimization with maximum-entropy regularization. Here, we briefly summarize the algorithm unique to our method. For general details of the SL mass reconstruction, readers are referred to excellent review papers (e.g., Kochanek 2006; Kneib & Natarajan 2011).

Gravitational lensing is a nonlinear mapping that transforms the source position β to the image position θ via the following lens equation:

$$\begin{eqnarray}&&{\boldsymbol{\beta }}={\boldsymbol{\theta }}-{\boldsymbol{\alpha }}({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 1 }$$

where α is the deflection angle. One can compute the deflection angle α via either the derivative of the deflection potential Ψ:

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}={\rm{\nabla }}{\rm{\Psi }}\end{eqnarray} \tag{ 2 }$$

or the convolution of the convergence κ:

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}({\boldsymbol{\theta }})=\displaystyle \frac{1}{\pi }\int \kappa ({{\boldsymbol{\theta }}}^{{\prime} })\displaystyle \frac{{\boldsymbol{\theta }}-{{\boldsymbol{\theta }}}^{{\prime} }}{| {\boldsymbol{\theta }}-{{\boldsymbol{\theta }}}^{{\prime} }{| }^{2}}{d}^{2}{\theta }^{{\prime} }.\end{eqnarray} \tag{ 3 }$$

The convergence κ is the dimensionless surface mass density defined as

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{{\rm{\Sigma }}}{{{\rm{\Sigma }}}_{c}},\end{eqnarray} \tag{ 4 }$$

where Σ (Σ_c ) is the (critical) surface mass density. The Σ_c is defined as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{c}=\displaystyle \frac{{c}^{2}{D}_{s}}{4\pi {{GD}}_{{ds}}{D}_{d}},\end{eqnarray} \tag{ 5 }$$

where D_ds is the angular diameter distance between the source and lens, D_s is to the source, and D_d is to the lens.

The first method (obtaining the deflection angle α by taking the derivative of the lensing potential Ψ; Equation (2)) has been advocated over the second one (Equation (3)) by some authors (e.g., Seitz et al. 1998; Bradač et al. 2005; Jee et al. 2007), who claim that the second method leads to biased results because the mass distribution outside the reconstruction field can affect the evaluation of the deflection angle (Equation (3)). However, this weakness of the second method can be overcome by simply and sufficiently extending the reconstruction field (as we demonstrate in Section 4.1). In this paper, we choose the second method because our experiment shows that the function minimization is much more stable than when we use the first method, although the field extension increases the number of free parameters by a significant factor.

Both κ and α are functions of the source redshift. Because our sources are at different redshifts, we first evaluate them at a reference redshift z_f and then scale the deflection angle for each source according to its redshift as follows:

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}({z}_{s})={\boldsymbol{\alpha }}({z}_{f})W({z}_{f},{z}_{s}),\end{eqnarray} \tag{ 6 }$$

where W(z_f , z_s ) is the cosmological weight:

$$\begin{eqnarray}&&W({z}_{f},{z}_{s})=\displaystyle \frac{D({z}_{f})D({z}_{l},{z}_{s})}{D({z}_{s})D({z}_{f},{z}_{l})}.\end{eqnarray} \tag{ 7 }$$

D(z) is the angular diameter distance to the redshift z while D(z_A , z_B ) is between the redshift z_A and z_B . In this paper, we set z_f = 9.

All multiple images from each source should converge to the same positions in the source plane when "delensed," and we minimize their scatters:

$$\begin{eqnarray}&&{\chi }^{2}=\sum _{i=1}^{I}\sum _{j=1}^{J}\sum _{k=1}^{K}\displaystyle \frac{{({{\boldsymbol{\theta }}}_{i,j,k}-{{\boldsymbol{\alpha }}}_{i,j,k}(z)-{{\boldsymbol{\beta }}}_{k})}^{2}}{{{\sigma }_{i}}^{2}},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{\beta }}}_{k}=\displaystyle \frac{1}{J}\sum _{j=1}^{J}\left({{\boldsymbol{\theta }}}_{i,j,k}-{{\boldsymbol{\alpha }}}_{i,j,k}(z)\right),\end{eqnarray} \tag{ 9 }$$

and I is the total number of sources, J is the total number of multiple images for each source, and K is the total number of "knots" for each multiple image. Here, "knots" refer to the positions of distinctive morphological features within a source if identifiable, which also must converge if the mass reconstruction is flawless. In this paper, we do not utilize the information (1) because our synthetic SL data do not provide them and (2) because A1689 has only a few sources with clear knots. Nevertheless, as explained below, we set up virtual knots as a convenient way to evaluate local magnifications. The notation σ_i in Equation (8) is the maximum scatter of the multiple-image positions from the ith source in the source plane and is given by

$$\begin{eqnarray}&&{\sigma }_{i}=\sqrt{{\sigma }_{x,i}^{2}+{\sigma }_{y,i}^{2}},\end{eqnarray} \tag{ 10 }$$

where σ_x,i and σ_y,i are the largest differences in β along the x- and y-axes, respectively. If we had used a fixed value for σ_i (for example, σ_i = 0001), the mass reconstruction would have been biased toward high-magnification results because a higher-magnification model can make the χ² term smaller. Using the current scheme prevents the converging model from favoring a high-magnification solution. Although one can consider a full lensing Jacobian matrix and compute the local magnification with it, in this paper, we implement this by simply creating two virtual knots 02 from the source center in such a way that the two vectors connecting the knots and the center are perpendicular to each other. Figure 1 illustrates the concept of our virtual knots in the image plane and the way that we evaluate σ_x,i and σ_y,i in the source plane.

Zoom In Zoom Out Reset image size

The following cross-entropy ⁵ is used to regularize our mass reconstruction to achieve smooth results and prevent overfitting:

$$\begin{eqnarray}&&R=\sum _{m,n}\left(p(m,n)-\kappa (m,n)+\kappa (m,n)\mathrm{ln}\displaystyle \frac{\kappa (m,n)}{p(m,n)}\right),\end{eqnarray} \tag{ 11 }$$

where p(m, n) and κ(m, n) are the prior and convergence at the position (m, n), respectively. This regularization is different from the one in Jee et al. (2007), where the the first two terms in Equation (11) are dropped and p and κ are normalized in such a way that the summation of each becomes unity. We find that this modification reduces the artificial flattening of the mass peak.

Our final target function to minimize is then

$$\begin{eqnarray}&&f={\chi }^{2}+{rR},\end{eqnarray} \tag{ 12 }$$

where r is a regularization control parameter to adjust the balance between χ² and R. Choosing too large values leads to oversmoothed results whereas too small values cause overfitting. In the latter case, the resulting mass map is not uniquely determined and shows a number of (unphysical) small-scale features. Because Equation (12), composed of the log-likelihood and prior terms, is based on the Bayesian principle, one may suggest a value, which makes the χ² value per component approach unity (e.g., Seitz et al. 1998). However, in our case, because σ_i (Equation (10)) in the χ² term (Equation (8)) is not a ∼68% error estimate but a maximum dispersion, this scheme is not directly applicable. Instead, we examine the relation between r and the resulting mean scatter in the source plane and choose the value in such a way that the mean scatter becomes ∼0001. We discuss the mass map dependence on the r parameter in Section 5.3.

Our mass reconstruction is performed in two steps. In the first one, we embed a 50 × 50 target reconstruction field in the center of a 70 × 70 grid, where the 10 pixel margin surrounding the target field is introduced to include the impact of the mass outside the 50 × 50 target reconstruction field. After the minimization of the function f in the first step, we double the resolution (140 × 140) and refine the result from the initial 70 × 70 solution. As explained below, our mass reconstruction is computationally intensive. This two-step implementation is found to speed up the convergence substantially compared to the case where we start with the 140 × 140 grid from the beginning. We considered different sizes for the margin and found that larger margins do not change the final result significantly. Because the observed multiple images do not lie exactly on the lattice point of the grid, we use bicubic interpolation to compute α at the image positions.

We start the function minimization from a flat prior. If the minimization is finished for the given prior, we restart the minimization after updating the prior with the new convergence map obtained in the previous minimization. This prior needs to be slightly smoothed to prevent overfitting. We applied a Gaussian convolution kernel whose smoothing scale is a σ = 0.6 (1.2) grid cell for the 70 × 70 (140 × 140) grid. This minimization cycle is repeated until there is no improvement.

Because the number of free parameters is extremely large (4900 and 19,600 for the 70 × 70 and 140 × 140 grids, respectively) and the function is highly nonlinear, caution is needed when we select a minimization algorithm. We choose the L-BFGS-B algorithm (Byrd et al. 1995; Zhu et al. 1997), which is based on the quasi-Newton method. The algorithm does not require a full Hessian matrix to determine descent directions, and this is a significant merit in our study because the full Hessian matrix evaluation is prohibitively expensive.

It is important to understand the limitation of the SL mass reconstruction in order to interpret the result quantitatively and compare it with other studies. In the LTM approach, it is straightforward to compute the parameter uncertainties because the number of free parameters is small and the system is in general overdetermined. However, because the quoted uncertainties are valid only within the LTM assumptions, one must include systematic errors due to the chosen LTM models when assessing the validity of the result.

In the non-LTM or free-form method, computation of the parameter uncertainties is nontrivial. In particular, if the minimization function does not contain any regularization term, the system is highly underdetermined and thus the solution is not unique. To bypass the difficulty, some authors perform many minimization runs and quote the standard deviation of these solutions. However, because the mean is not one of the valid solutions, it is questionable whether the standard deviation can be used as a reliable uncertainty proxy.

Our MARS free-form method provides a unique solution because the maximum-entropy term makes the system overdetermined. Obviously, this is a significant advantage over other free-from methods and allows us to estimate the parameter (i.e., mass cell) uncertainties from the Hessian matrix. Although the matrix is extremely large (70² × 70² and 140² × 140² for the 70 × 70 and 140 × 140 grids, respectively), the computation is needed only once after the convergence is reached. Also, because evaluation of each matrix element is independent, computation of the entire matrix elements can be easily parallelized. This Hessian-based uncertainty estimation is also used in Jee et al. (2007). However, their parameters are potentials and thus complex error propagation using all covariances was needed to obtain mass uncertainties. Our Hessian matrix is a square matrix whose elements are the second derivatives ⁶ of the function f (Equation (12)). Under the assumption that the posterior of each κ value follows the Gaussian distribution, we can use the inverse of the Hessian matrix to obtain the covariance of κ.

The covariance matrix needs to be rescaled if the χ² value (Equation (8)) per multiple image departs from unity when the convergence is reached. We adopted σ = 002 to be our fiducial mean source plane scatter and multiplied the per-image χ² value to our covariance matrix. In Section 5.4, we discuss how the uncertainty depends on our choice for the mean source plane scatter.

We note that the uncertainty map computed in this way has some limitations. First, it is derived under the assumption that the error distribution is Gaussian. Also, the values are valid only for the given prior. Finally, our uncertainty map does not include systematic errors due to the regularization and other method-specific artifacts.

To test our SL reconstruction method, we use the two mock clusters from M17, whose SL data are publicly available. ⁷ During the review process of the current article, we were informed that the "true" deflection field data released by M17 are not self-consistent with their multiple-image catalogs. Our independent investigation reveals that indeed the scatters in the source plane measured with their "true" deflection field data are as large as ∼01 for both clusters. Note that the results in this paper are obtained with the same multiple-image catalogs available to the M17 participants. Here we provide only a brief summary on how these two mock clusters are generated. For details, we refer readers to M17.

The mock galaxy cluster Ares at z = 0.5 is created with the semianalytic code MOKA (Giocoli et al. 2012). Ares consists of three components: (1) the main dark matter halo, (2) cluster members, and (3) the brightest cluster galaxy (BCG), which all follow analytic profiles. The SL data are produced under the flat ΛCDM model with Ω_m = 0.272 and H₀ = 70.4 km s⁻¹ Mpc⁻¹. The projected mass distribution shows many sharp "spikes" at the location of galaxies due to the cuspy centers of the analytic profile halos.

The other mock galaxy cluster Hera at z = 0.507 is generated from a high-resolution (m_DM = 10⁸ h⁻¹ M_⊙) zoom-in resimulation of the initial low-resolution cosmological simulation of Planelles et al. (2014), which assumes a flat ΛCDM with Ω_m = 0.24, Ω_b = 0.04, σ₈ = 0.8, and H₀ = 72 km s⁻¹ Mpc⁻¹. Compared to Ares, the galaxy-scale structures are somewhat smoothed because of the finite resolution of the simulation.

Ares and Hera differ greatly in mass. Ares comprises two halos, whose virial masses are 1.32 × 10¹⁵ h⁻¹ M_⊙ and 8.8 × 10¹⁴ h⁻¹ M_⊙ for the southeastern and northwestern halos, respectively. The projected distance between the two halos is ∼400 h⁻¹ kpc. On the other hand, the total mass of Hera is 9.4 × 10¹⁴ h⁻¹ M_⊙. This cluster is also bimodal as in the case of Ares, comprising the eastern and western halos. The projected separation is much smaller (∼130 h⁻¹ kpc) than that of Ares.

Because Ares is more massive than Hera, the angular extent of the multiple-image distribution is also larger. For Ares, we choose the central 200'' × 200'' region as our reconstruction field, which contains 242 multiple images from 85 sources. For Hera, we reconstruct the central 100'' × 100'' region and there are 65 multiple images from 19 source galaxies. We display the color-composite images of Ares (Hera) and its multiple-image distribution in Figure 2 (3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The galaxy cluster A1689 at z = 0.18 has been the subject of a number of SL studies (e.g., Broadhurst et al. 2005; Zekser et al. 2006; Limousin et al. 2007; Coe et al. 2010; Diego et al. 2015). We compiled multiple images from the following three studies: Broadhurst et al. (2005, hereafter B05), Limousin et al. (2007, hereafter L07), and Coe et al. (2010, hereafter C10). The multiple images of A1689 used in the current study are listed in Table 1, which contains 109 multiple images from 34 source galaxies. We classify them into the "gold," "silver," and "gold candidate" samples. The gold images have spectroscopic redshifts. The silver images have only photometric redshifts, which all three papers agree to be bona fide multiple images. The gold candidates are the multiple images whose spectroscopic redshifts were unknown to the authors in the three studies, but later measured by Bina et al. (2016, hereafter B16). There are 73 multiple images from 22 source galaxies in the gold sample, 31 multiple images from 10 source galaxies in the silver sample, and 5 multiple images from 2 source galaxies in the gold candidate sample. We display the multiple-image distribution used in this study in Figure 4, which shows the central 200'' × 200'' region, where we perform the SL mass reconstruction. We follow the numbering scheme of B05 while adopting the coordinates of C10.

Zoom In Zoom Out Reset image size

4.1.1. Convergence Map Comparison

Figure 5 shows the truth, reconstructed convergence, and uncertainty maps of Ares and Hera. Our reconstruction closely traces the main clumps of the truth. In the case of Ares, the reconstructed convergence map reveals the two main clumps located in the NW and SE regions. The small galaxy-scale spikes are missing in our mass reconstruction because their features are smaller than our gird cell. Instead, our convergence shows small "wrinkles," which we believe are artifacts due to the resolution limit (for example, it is difficult to make two multiple images within a grid cell converge in the source plane).

Zoom In Zoom Out Reset image size

The reconstructed mass map of Hera also closely follows the truth for the main clumps. We find that the artifacts mentioned for the Ares reconstruction are substantially reduced here because the true mass distribution, on which the SL data are based, is much smoother. In Hera, one noticeable discrepancy between our reconstruction and the truth is the location of the eastern peak, which is shifted south by ∼10 kpc with respect to the truth. We attribute the offset to the combined effect of the relatively insufficient number of multiple images around this eastern peak and the regularization.

Our uncertainty map shows that the uncertainty varies from σ_κ ∼ 0.1 (0.035) to ∼0.45 (0.11) for Ares (Hera). The uncertainty distribution is consistent with our expectation that in general the error is large where the convergence is high. Also, the error is somewhat reduced at the location of the multiple images.

M17 present residual [(κ − κ_true)/κ_true] convergence maps for the nine different SL mass reconstruction algorithms, which were used in their SL modeling comparison project. M17 show that in general the LTM models perform better for both Ares and Hera. This is somewhat expected (1) because the two mock clusters were created under the LTM assumption and (2) because analytic profiles are reasonably good priors in the outer regions where there are no multiple images. In M17, the performance difference between the LTM and free-form methods is greater for Ares, which comprises cuspy parametric halos. For comparison with the results in M17, we also display our residual convergence maps with the same color table and range in Figure 6. The residual map for Ares shows that our mass reconstruction quality is similar to those of the LTM models in the region where the SL data are present. Compared with GRALE, which leads the free-form method in performance in M17, our mass reconstruction gives smaller residual convergence values over a larger area with less large-scale features. In the case of Hera (right panel of Figure 6), our reconstruction is superior to any of the M17 results in terms of the amplitude and uniformity of the residuals (see Figure 10 in M17). All LTM residuals possess galaxy-scale spikes, which we attribute to mismatches between the assumed and true halo shapes. All free-form methods suffer from large-scale residuals.

Zoom In Zoom Out Reset image size

As for the mass estimation of galaxy-scale substructures, MARS cannot constrain them because the current grid-based representation cannot distinguish the projected mass belonging to the cluster halo from those bound to the substructures. In addition, these galaxy-scale substructures are smaller than our grid cell. M17 shows that some LTM-based methods outperform free-from algorithms for the substructure mass recovery, which is not surprising when we consider the fact that the LTM methods explicitly include the galaxy-scale halos into their models.

4.1.2. Radial Profile Comparison

We compare the radial convergence and cumulative mass profiles in Figures 7–10. M17 state that the SL data are expected to provide tight constraints at 20'' < r < 60'' (10'' < r < 30'') for Ares (Hera). Thus, we choose the radial ranges for comparison accordingly. Because both Ares and Hera comprise two main halos, we measure the radial profiles twice for each cluster using the two centers of the main clumps. Note that M17 use only the stronger peak (the southeastern peak for Ares and the western peak for Hera) for each cluster to measure the profiles. Overall, both radial convergence and cumulative mass profiles are in excellent agreement with the truth. In the study of M17, the three best-performing LTM methods produce the mass profiles for Ares, which differ from the truth by <2%, while most free-form methods depart from the truth by 5%–15%. Our result (see the right panel of Figure 8) shows that the agreement of the cumulative mass profile in the 20'' < r < 60'' radial range is at the ≲1% level. When we choose the northwestern peak as the reference (note that this measurement is not presented in M17), the difference starts to increase at r > 30 rapidly and reaches ∼10% at r = 60'' (see the left panel of Figure 8). Because the northeastern peak has a smaller (≲30'') Einstein radius, this behavior is not surprising. In the case of Hera, the mass profiles agree with the truth within ≲1% for the radial range 10'' < r < 30'' regardless of the center choice. The best LTM methods in M17 show agreements at the ∼3% level. The free-form results in M17 differ from the truth by 5%–10%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1.3. Magnification Comparison

In Figure 11, we display the magnification maps for Ares and Hera. Also shown in Figure 11 are the true magnification maps that we computed from the publicly available convergence and shear data. The overall shape and size of the reconstructed critical curves are well recovered, although the largest discrepancies in magnification are found near the critical curves as shown in the residual maps (Figure 12). As noted by M17, the critical curves define the regions where the magnification diverges and thus small differences in the mass reconstruction can cause large magnification errors near them. M17 report that some free-form approach such as the Bradac-Hoag, Coe, and Diego-multire methods produce highly fluctuating irregular critical lines, which are not observed in the truth. As shown, our entropy-based regularization enables us to obtain relatively smooth critical curves. Following M17, we also compute the correlation between our and the truth magnification maps and display the result in Figure 13, which indicates that the accuracy (i.e., departure of the median from the perfect correlation) is within ∼10% in the range 1 < μ < 30 whereas the precision (i.e., distance between the 25th and 75th percentiles) is ∼20% at μ = 10. Note that some details in our magnification comparison methods may be different from M17. For example, we do not use exactly the same region for comparison. With this caveat, we report that our prediction on the median magnification is accurate within 10% for both Ares and Hera. In the case of Ares, we report that our magnification reconstruction results are significantly better than some (Bradac-Hoag, Diego-multires, Lam, and Coe) of the non-LTM methods while the Diego-reggrid and GRALE methods show a similar performance to ours. Clearly, most LTM methods outperform MARS for the Ares magnification reconstruction. In the case of Hera, the performance of MARS is similar to those of GLAFIC and Diego-overfit, which are the best-performing methods in the LTM and free-form algorithms, respectively. Because the two mock clusters are generated under the LTM assumption, there exists a strong correlation between dark matter and luminous galaxies in the mock data. As noted by M17, it is possible that LTM approaches can benefit from the use of this extra information. M17 show that most LTM methods perform better in the case of Ares than Hera, which may hint at this possibility because Ares has tighter dark-matter–galaxy correlations. In addition, the critical curves from the LTM methods are much smoother by construction, which is a significant advantage because even small irregular features around the critical curves in non-LTM methods critically contribute to the degradation of the magnification correlations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we focus mainly on the convergence, uncertainty, magnification, and radial profile. In Appendix C and D, we present our reconstruction of the multiple images in the source plane and the giant arcs in the image plane.

4.2.1. Convergence and Uncertainty Maps

We display the reconstructed κ distribution and uncertainty map of A1689 in Figure 14. Also, Figure 15 shows the κ contour overlaid on the color-composite image of the cluster. Readers are reminded that here we only present the result obtained with the gold sample. We discuss how the result changes when we include the gold candidate or silver samples in Appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mass map shows that the projected mass distribution of A1689 is characterized by the dominant mass peak centered on the BCG at (α, δ) ∼ (13^h11^m29^s, −1°20'30'') and the secondary, but much weaker, mass peak centered on the bright galaxy at (α, δ) ∼ (13^h11^m32^s, −1°20'00''). This feature has also been observed in some previous LTM-based studies (e.g., Broadhurst et al. 2005; Halkola et al. 2006; Zekser et al. 2006; Limousin et al. 2007). Some free-form methods do not constrain the secondary peak mentioned above as the characteristic feature of A1689. For example, Coe et al. (2010) obtained a mass map with several substructures whose projected densities are comparable and do not identify our secondary clump as a distinguishing feature of A1689. Diego et al. (2005) presented a result from averaging over 1000 solutions. Although their result shows a weak secondary substructure near our mass peak, it is offset by ∼20'', which is surprising if true within the LTM paradigm. We note that in their revised analysis (Diego et al. 2015), this secondary peak coincides with the galaxy.

Given that our mass map is the solution that successfully delenses all multiple images to converge in the source plane with a negligible scatter of ∼0001, its smoothness is remarkable for a free-from result. As already mentioned, the mass reconstructions from Coe et al. (2010) and Diego et al. (2005) suffer from high-amplitude small-scale fluctuations. Mohammed et al. (2014) produce a smooth mass reconstruction of A1689 using the GRALE method, which does not assume LTM. However, because the authors still assume that the cluster consists of Plummer spheres, this smoothness is achieved by design.

We also present the multiple-image scatters in the image plane in Figure 16. Except for the two largest scatters at ∼0025 and ∼0073, we find that most multiple images in the A1689 field have a small (<001) scatter. We note that these scatters are evaluated without including false multiple images, which would produce catastrophically large scatters. Because MARS is designed to minimize scatters only in the source plane, the model cannot prevent false multiple images. This weakness is shared by the SL reconstruction algorithms that employ only source plane minimization.

Zoom In Zoom Out Reset image size

Our scrutiny of the mass map indicates that the mass map contains several "pinches" or "quadrupoles" characterized by alternating over- and underdense regions. We believe that they are artifacts caused by the resolution limitation and happen where compact halos generate multiple images densely distributed in small areas. As discussed in Section 4.1, the mass reconstruction for Ares produces many such features while no such artifact is seen in the result for Hera. The uncertainty map (middle panel of Figure 14) somewhat resembles the mass reconstruction (top panel of Figure 14), which is similar to our mock-cluster cases (Section 4.1). Near the central mass peak, this convergence uncertainty (δκ) is high (δκ ∼ 0.8) whereas near the field boundary it is low (0.3 ≲ δ κ ≲ 0.4). When rescaled by the local convergence, the relative (δ κ/κ) error map (bottom panel) shows the opposite trend, with the estimate being high (δ κ/κ ∼ 1) near the field edges, which is consistent with the relative difference map shown for the mock cluster (Figure 6). Note that the uncertainties are relatively low at the location of the multiple images.

4.2.2. A1689 Magnification

Figure 17 presents our critical curves of A1689. The global morphology comprising the inner (radial) r ∼ 15'' and outer (tangential) r ∼ 40'' critical curves is similar to the findings in previous studies. However, details differ greatly among the studies, including ours. Therefore, despite the fact that A1689 is one of the most studied SL clusters, the community still needs to converge on the magnification details and thus the interpretation of high-z galaxies found behind the cluster. We repeat that because the relation between mass and magnification distributions are highly nonlinear, even small differences in the convergence lead to large changes in the magnification.

Zoom In Zoom Out Reset image size

As seen in the mock-cluster results, the smoothness of our critical curve is noteworthy when compared with other free-form results. For example, because no regularization is imposed in Coe et al. (2010), their critical curves show highly oscillatory features and vary widely among different solutions, which, however, are all capable of converging all used multiple-image positions in the source plane. In Diego et al. (2005), only the critical curve estimated from the mean (out of 1000) convergence map is shown.

4.2.3. A1689 Radial Profile

The experiments with the mock-cluster SL data (Section 3.1) show that our mass reconstruction recovers the cluster radial profiles with superb (≲1%) accuracy in the regions where multiple images are present. Thus, if all multiple images are real in the gold sample, we expect that in the current study a similarly high-fidelity mass profile can be obtained for A1689 in the SL regime.

Our test with the mock-cluster data shows that our mass reconstruction recovers well the central profile down to r ∼ 4'' (∼25 kpc). At smaller radii, the density becomes nonnegligibly underestimated because of the limitations in resolution and regularization. Considering the plate scale and resolution of our mass reconstruction for A1689, we believe that the smallest radius at which we can restore the density without severe bias is ∼12 kpc. Given the extent of the multiple images, the mass density is expected to be constrained by the data out to ∼240 kpc.

We display the projected radial density profile of A1689 in the top panel of Figure 18. Also plotted is the best-fit (c₂₀₀ = 5.53 ± 0.77, and ${r}_{s}={538}_{-100}^{+90}$ kpc) NFW profile. This NFW profile is a reasonable approximation of the overall shape of our radial profile in the SL regime, although we observe some indications of the densities predicted by the best-fit NFW model being 1σ higher at r ≲ 15 kpc and lower at r ≳ 200 kpc. Because there exist several multiple images at r ≲ 15 kpc, it is unlikely that the shallower (than the best-fit NFW prediction) profile is entirely due to the aforementioned artifact, but is somewhat constrained by the data. The "shoulder" (see the arrow in the top panel) at r ∼ 160 kpc coincides with the location of the secondary mass peak at $(\alpha ,\delta )=({13}^{{\rm{h}}}{11}^{{\rm{m}}}{32}^{{\rm{s}}},-1^\circ 20^{\prime} 00^{\prime\prime} )$, which indicates that our best-fit NFW model might be nonnegligibly affected by the presence of the substructure. The cumulative projected mass profile (the bottom panel of Figure 18) shows that the aperture mass of A1689 is ∼4 × 10¹⁴ M_⊙ at r = 237 kpc, which is in good agreement with the previous studies (e.g., Broadhurst et al. 2005; Halkola et al. 2006; Limousin et al. 2007). At face value, the best-fit concentration value c₂₀₀ = 5.53 ± 0.77 is slightly lower than the results in some previous studies. However, when we limit our comparison to those who used only the SL data (e.g., Broadhurst et al. 2005; Limousin et al. 2007; Coe et al. 2010), our measurement is consistent at the ∼1σ level.

Zoom In Zoom Out Reset image size

As mentioned in Section 2.1, our final reconstruction is achieved in two steps (i.e., initial reconstruction in the 50 × 50 grid followed by a second reconstruction in the 100 × 100 grid ⁸ ). If significant systematic differences are found between the two results, this might indicate that the resolution in the second step is insufficient and the result should be further improved with even a higher-resolution grid. Here we demonstrate that although the larger grid provides a higher resolution in some cases, no significant systematic difference is found.

Figure 19 compares the 50 × 50 and 100 × 100 versions side by side for the three clusters studies in this paper. It is clear that the larger grid improves the resolution of the resulting mass map. In addition, the 100 × 100 versions reduce quadrupole-like artifacts mentioned in Section 4.2.1. In particular, the artifact reduction is most significant in Ares, whose truth map possesses many sharp spikes, which cannot be sufficiently represented even with the 100 × 100 resolution.

Zoom In Zoom Out Reset image size

We also notice that the convergence values at the mass peaks are slightly reduced in the higher-resolution version. We speculate that this might arise from our two-step mass reconstruction method, which obtains the lower-resolution version first and uses the smoothed version as a prior to the second reconstruction. However, the differences make only a negligible impact on the radial mass profiles.

Therefore, we conclude that although the higher-resolution grid improves the resolution and reduces some artifacts, the change in grid resolution does not lead to any significant systematic differences that lead to changes in scientific interpretation.

In order to show that our algorithm provides quasi-unique solutions, we demonstrate that the mass reconstruction does not depend on our initial guesses. In Figure 20, we present the A1689 mass reconstruction results when we start the function optimization from four different initial conditions: two flat mass sheets and two NFW profiles.

Zoom In Zoom Out Reset image size

The 2D mass map comparison shows that the agreement among the four results is excellent in the region constrained by multiple images (to guide eyes, we plot the locations of the multiple images with black dots). The mass map in the outskirt where no SL data exist somewhat depends on the initial conditions, which however is not surprising because no constraints are provided there.

For a more quantitative analysis, we compare the radial profiles from these four mass maps in the bottom panel of Figure 20. As indicated by the 2D comparison, no significant difference is found out to r ∼ 200 kpc. Beyond this radius, some small systematic differences are noticeable. In Figure 21, we display the residual mass maps.

Zoom In Zoom Out Reset image size

In summary, we find that MARS achieves a quasi-unique solution for the region constrained by the SL data. We note that the optimization requires a much longer time if the initial setups are far from the final result (e.g., flat initial conditions).

Equation (12) consists of the χ² and regularization terms. The regularization term is to prevent overfitting and stabilize the convergence to a quasi-unique solution. In order to illustrate the dependence of the mass map on the regularization control parameter r, we display the three (r = 0, 0.1, 1, and 10) cases in Figure 22. The r = 1 case corresponds to our fiducial result, which leads to a mean source plane scatter of ∼10⁻³''. The r = 0.1 (10) results yield a mean source plane scatter of ∼10⁻⁴'' (∼10⁻²''). The differences among the r = 0.1, 1, and 10 cases are small, but noticeable. As expected, the smaller (larger) r value produces a less (more) smoothed mass map.

Zoom In Zoom Out Reset image size

The r = 0 mass map is obtained without any entropy prior and it is not surprising that the resulting mass map shows a number of unphysical, small-scale fluctuations. Because the system is underdetermined, the solution is not unique and the result displayed in Figure 22 is just one of many possible solutions.

As explained in Section 2.2, we estimate the mass reconstruction errors with the Hessian matrix by assuming that the posteriors follow Gaussian distributions. Another important assumption is the size of the multiple-image position scatters in the source plane. For our fiducial result, we assume σ = 002 after considering the finite resolution of our mass grid. Because it is difficult to define/estimate the exact value for this scatter from the first principle, here we examine the dependence of the mass reconstruction error estimation on the size of the scatter.

In Figure 23, we show the uncertainty maps of Hera and Ares for the three additional cases: σ = 0005, 005, and 05. We find that the relative scales across the field are not affected by the assumed scatter values, which is somewhat expected because changing the scatter is equivalent to scaling the likelihood with a simple multiplicative factor. In terms of the absolute values, we notice significant changes. For Ares, the decrease in error is nearly an order of magnitude when we change σ = 0005 to σ = 005. For the case of σ = 05, the decrease is about 50% relative to the σ = 005 case. The sensitivity is reduced for Hera, and we attribute this to the difference in the number of multiple images.

Zoom In Zoom Out Reset image size

In Figure 24, we show the relative uncertainty (δ κ/κ) maps of Hera and Ares for the three cases. Similarly to Figure 23, the uncertainty values scale with the size of the source plane scatter.

Zoom In Zoom Out Reset image size