CURLING – I. The influence of point-like image approximation on the outcomes of cluster strong lens modelling

ABSTRACT

Cluster-scale strong lensing is a powerful tool for exploring the properties of dark matter and constraining cosmological models. However, due to the complex parameter space, pixelized strong lens modelling in galaxy clusters is computationally expensive, leading to the point-source approximation of strongly lensed extended images, potentially introducing systematic biases. Herein, as the first paper of the ClUsteR strong Lens modelIng for the Next-Generation observations (CURLING) program, we use lensing ray-tracing simulations to quantify the biases and uncertainties arising from the point-like image approximation for JWST-like observations. Our results indicate that the approximation works well for reconstructing the total cluster mass distribution, but can bias the magnification measurements near critical curves and the constraints on the cosmological parameters, the total matter density of the universe Ω_m, and dark energy equation of state parameter w. To mitigate the biases, we propose incorporating the extended surface brightness distribution of lensed sources into the modelling. This approach reduces the bias in magnification from 46.2 per cent to 0.09 per cent for μ ∼ 1000. Furthermore, the median values of cosmological parameters align more closely with the fiducial model. In addition to the improved accuracy, we also demonstrate that the constraining power can be substantially enhanced. In conclusion, it is necessary to model cluster-scale strong lenses with pixelized multiple images, especially for estimating the intrinsic luminosity of highly magnified sources and accurate cosmography in the era of high-precision observations.

1 INTRODUCTION

Galaxy clusters, with total masses of |$\rm 10^{14}\sim 10^{15} \, {\rm M}_{\odot }$|⁠, represent the most massive self-gravitationally bound structures in the Universe. As predicted by Einstein’s general relativity (GR), these massive clusters exhibit gravitational lensing effects due to their dense mass distribution. This phenomenon leads to the deflection and strong lensing of light rays from background sources, such as distant galaxies, quasars, and even star candidates, resulting in the formation of multiple images or giant arcs (Yuan et al. 2012; Liang et al. 2014; Caminha et al. 2017; Mahler et al. 2018). Gravitational lensing by galaxy clusters has evolved into a powerful tool for studying various aspects of astrophysics and cosmology, for instance, providing insights into the origin and evolution of galaxies (Treu et al. 2015; Mahler et al. 2019; Mercurio et al. 2021; Welch et al. 2022), the characteristics of dark matter (Harvey et al. 2015; Massey et al. 2015; Caminha et al. 2019; Jauzac et al. 2019; He et al. 2020; Bergamini et al. 2023, Lin et al. 2023), and enabling precise constraints of cosmological parameters (Gilmore & Natarajan 2009; Jullo et al. 2010; Acebron et al. 2017; Grillo et al. 2018, 2020; Caminha et al. 2022a, Wagner et al. 2019).

Accurate lens modelling is crucial for achieving reliable scientific results. However, this task is challenging due to the intricate nature of strong lensing clusters. These clusters, characterized by the substantial lensing cross-section, consist of multiple clumpy dark matter haloes and numerous substructures in a dynamically perturbed state. Over the past few decades, efforts have been made to address various sources of systematic errors associated with lens modelling. For example, errors can arise in parametric methods, where the assumption of analytical models will inevitably deviate from the true mass distribution (Acebron et al. 2017; Meneghetti et al. 2017; Raney et al. 2020). The presence of large-scale structures along the line of sight of the lens system can impact lens modelling (D’Aloisio & Natarajan 2011; Chirivì et al. 2018), but is often overlooked in most of the lensing analyses. Given the limitations on accurately measuring the magnitude and velocity dispersion of each cluster galaxy, it is not possible to depict each mass component in a cluster lens, instead galaxy scaling relations are adopted on the basis of mass-to-light relations (Faber & Jackson 1976; Tully & Fisher 1977), the scatter of these galaxy scaling relations can result in biased modelling outputs (Bergamini et al. 2021; Granata et al. 2022). The accuracy of lens modelling is also sensitive to the measurements of multiple image redshifts, as demonstrated by Acebron et al. (2017) who analysed the difference on cosmological constraints between using spectroscopic and photometric data.

In addition to the systematic errors as mentioned earlier, another source of error in strong lens modelling arises from the positional error associated with multiple images. During the optimization of lens parameters, accurate likelihood computation involves utilizing a pixelated approach to describe the image light profile and taking into account the signal-to-noise ratio of each pixel (Dye & Warren 2005; Suyu et al. 2006; Kneib & Natarajan 2011; Birrer & Amara 2018; Nightingale, Dye & Massey 2018; He et al. 2023, Wagner et al. 2017, Lin et al. 2022). However, due to the computational limitations, current algorithms for modelling strong lenses at cluster scale simplify the observed multiple images as point sources, assuming a Gaussian distribution for the associated errors. This assumption relies on certain conditions, such as the compactness of the source galaxy, and the smoothness of the surface brightness profile. Given that the morphology of a lensed object is typically more complex than a simple Gaussian shape, relying on this simplification in lens modelling can introduce notable biases. These biases can subsequently impact the accuracy of the derived lensing properties.

Observational efforts hold the potential to significantly get rid of the systematic errors and improve the accuracy of current lens modelling. For instance, by incorporating the data from Multi-Unit Spectroscopic Explorer on the Very Large Telescope (Bacon et al. 2012) of each member galaxy, it becomes possible to estimate galaxy masses with greater accuracy than by solely relying on scaling-relation models. Such data also enable the spectroscopic confirmation of multiple images, thereby mitigating the risks of misidentifying image families and enhancing the robustness of modelling. Consequently, the errors arising from the modelling process are likely to become more pronounced.

To investigate the impact of assuming multiple images as point sources on cluster strong lensing analyses, we employ simulated cluster strong lensing systems to study their model-derived mass profiles, magnification maps, and cosmological parameters. The paper is organized as follows. In Section 2, we provide an overview of the simulated cluster lenses. We briefly introduce the theory of strong lensing analyses, including the generation of mass map and magnification map, the constraint on cosmology, and the modelling procedure in Section 3. Our main results are presented in Section 4. We conclude with a discussion and summary of our work in Section 5. Throughout the paper, we consider a flat w-cold dark matter cosmology (Ω_{Λ, 0} = 0.7, Ω_m = 0.3, and h  = 0.7) with a constant dark energy equation of state parameter w = w₀ = −1.

2 LENSING RAY-TRACING SIMULATIONS

In this section, we outline the process of setting up the simulation of cluster lens systems. In the sense that stacking the constraints from individual clusters can be an effective approach to yield competitive estimates of the values of the cosmological parameters, we include four lens systems in the simulation, each consisting of two parts: the galaxy cluster acting as the strong lens, and the background sources within the cluster field that are lensed into multiple images.

2.1 Cluster lenses

To fully exploit the potential of cluster strong gravitational lensing effect, it is crucial to select massive clusters with higher lensing efficiency as the targets. In light of this, we choose to build the clusters based on the strong lensing clusters that have been extensively studied as part of the Hubble Frontier Field (HFF) program (Lotz et al. 2017). We replicate the publicly available best-fitting model parameter sets and the galaxy catalogues¹ of MACSJ0416.1−2403 (Bergamini et al. 2021), A2744 (Bergamini et al. 2023), MACSJ1206.2-0847, and AS1063 (Bergamini et al. 2019), to describe the mass distribution of the four mock clusters. These referenced models, built with high-quality and extensive spectrophotometric data, along with the largest number of secure multiple images identified so far, provide robust mass distributions for depicting clusters with great strong lensing strengths.

The mass of a mock cluster is thus composed of several components: (i) dark matter haloes associated with the large-scale, smooth clumps, (ii) external clumps representing surrounding substructures within the cluster’s strong lensing field, (iii) the brightest cluster galaxies (BCGs), and (iv) galaxy-scale subhaloes corresponding to the member galaxies within the cluster.

Each mass component in the cluster is constructed using an analytic pseudo-isothermal elliptical mass distribution (PIEMD, Kassiola & Kovner 1993; Elíasdóttir et al. 2007) profile. The PIEMD density is expressed as:

where ρ₀ represents the central density of the profile, r_core and r_cut are the core radius and the truncation radius, respectively. Within the regions where r_core < r < r_cut, the profile behaves as ρ ∼ r⁻², while it transitions to ρ ∼ r⁻⁴ in the outer regions where r > r_cut. To define the potential of a PIEMD profile in lenstool (Jullo et al. 2007), we utilize a parameter set consisting of { x, y, e, θ, |$r_{\mathrm{ core}}$|⁠, |$r_{\mathrm{ cut}}$|⁠, σ }, where x, y describe the centre of the potential, e and θ define its ellipticity and position angle, σ represents the velocity dispersion. The masses of the galaxy-scale subhaloes are determined based on the scaling relations (Granata et al. 2022), that are related to their respective galaxy magnitudes, following the equations:

we fix the parameters of slopes α = 4.0 and β = 4.0.

Given that the typical redshift of a strong lensing cluster is z ≃ 0.2–0.5 (Robertson et al. 2020), and the strongest cluster lenses peak at z ≃ 0.4 (Fox et al. 2022), we set the redshifts of the four mock clusters to z  = [0.35, 0.4, 0.45, 0.5].

We summarize the properties of the mock clusters in Table 1.

2.2 Background sources

In order to simulate the multiple images lensed by the cluster in each system, it is important to not only define the properties of the lens, but also characterize the population of background sources, including their number counts, redshifts, and spatial distribution.

First, to address the number of multiply lensed sources in a cluster field, we rely on the deeply resolved cluster images obtained from HFF and JWST (Zitrin et al. 2013; Caminha et al. 2022b; Bergamini et al. 2023; Mahler et al. 2023), which have depths of up to 30 AB magnitude, similar to the expected depths of next-generation surveys for studying strong lensing clusters. These surveys have been demonstrated to have the ability of observing up to hundreds of multiple image families. However, in addition to the identification of multiple images under the consideration of observational depths, we also emphasize the importance of extensive spectroscopic follow-up observations at the cores of galaxy clusters. This is crucial for ensuring robust measurement of the source redshifts, which is essential for accurate lens modelling (Grillo et al. 2015; Bergamini et al. 2021). Based on these criteria, we assume that a cluster contains more than 20 image families that can be identified and spectroscopically verified.

The source redshift distribution is generated based on the photometric redshift catalog of the Hubble Ultra-Deep Field (HUDF, Beckwith et al. 2006), as shown in Fig. 1. Note that the redshift distribution of HUDF galaxies starts from z = 0.025, however, in our simulation, we exclude the redshifts below a certain threshold due to their low lensing efficiency. As a result, the sources in our sample fall within the redshift range of z = 1.0–8.4.

In each lens system, we use the same redshift distribution for the sources, while their positions vary due to the diverse mass distributions of the cluster lenses. The source positions are determined using the ray-tracing algorithm outlined in Li et al. (2016), which involves numerically solving the lens equation. With the mass distribution of the lens as introduced in Section 2.1 ready, next is to calculate the deflection angle map and build the lens equation. The final step is to map the locations of pixels on the image plane to the source plane, by placing the source on the source plane, light rays of the source are mapped from the source plane onto the image plane again to generate the lensed image. A source position is determined provided that more than one lensed image is generated. The ray-tracing and source-locating are performed with a pixel size of 0.05 arcsec within a field of 180 arcsec × 180 arcsec. We note that some of the images are missed due to grid resolution in lenstool while performing lens modelling. Accordingly, we exclude the sources that are found by lenstool with less than one image from the catalogue. Therefore, the number of lensed sources in each lensing system is slightly different, we summarize the number as |$N_{\mathrm{ s}}$| in Table 1.

In addition to determining the positions of the sources, for the novel modelling method proposed in this work, we also compute the extended source surface brightness distribution with the above ray-tracing algorithm. For simplicity, we assume that the light of each source follows a S|$\rm \acute{e}$|rsic profile, with a random ellipticity q and position angle θ. The size and S|$\rm \acute{e}$|sic index of the source are given based on their evolution with redshift (Mowla et al. 2019). Therefore, the final range of semimajor axis is from 0.3 to 0.8 arcsec, and the range of S|$\rm \acute{e}$|sic index varies from 1.0 for sources at higher redshifts to 4.5 for lower redshift sources. The superimposed surface brightness distribution S_obs within the field of view, summed from each source, then becomes the observable in our pixelated modelling method.

We have constructed the cluster strong lens system in a comprehensive manner, taking into account all the image families expected to be observed in next-generation surveys. None the less, in the context of comparing traditional lens and pixelated modelling, which can be computationally expensive, we randomly select only three of these image families as constraints for the following analyses. The redshifts of these selected sources are marked in Table 1.

3 METHODOLOGY

In this section, we describe the methods of analysing the simulations generated in Section 2. First, we introduce the concept of measuring the cluster mass profile and magnification map based on the lensing formalism. We also discuss the capability of cluster strong lenses as a probe for constraining cosmological parameters. Subsequently, we provide a comprehensive explanation of our modelling methods, encompassing both the traditional lens modelling and a novel pixelated algorithm proposed in this work.

3.1 Lensing theory

The lens equation relates the position of source $\beta$ at redshift z_s to the position of image θ at redshift z_l:

where |$\boldsymbol{\alpha }$| is the reduced deflection angle as a result of lensing potential. The mapping from unlensed coordinates on the source plane to that on the image plane can be described by the Jacobian matrix,

defining the convergence |$\kappa = \frac{1}{2} \Delta \Psi = \frac{1}{2}(\Psi _{11} + \Psi _{22})$| as a function of the effective lensing potential Ψ, to delineate the magnification of images, the shear vectors |$\gamma _1 = \frac{1}{2}(\Psi _{11}-\Psi _{22})$| and γ₂ = Ψ₁₂ = Ψ₂₁, where the complex shear γ = γ₁ + iγ₂ describes the stretching of images, we can write the Jacobian matrix as

The magnification of the lens is then given by the inverse of the determinant of the Jacobian matrix,

tangential and radial magnification factors are defined as

respectively. Regions with λ_t = 0 and λ_r = 0 are defined as critical curves in the lens plane, and caustics in the corresponding source plane, in which the magnifications are ideally infinite. These areas are particularly advantageous for observing distant and faint galaxies.

As the definition of convergence |$\kappa = \frac{\Sigma }{\Sigma _{\rm crit}}$|⁠, where |$\Sigma _{\rm crit} = \frac{c^2}{4\pi G} \frac{D_{OS}}{D_{LS}D_{OL}}$| is the critical surface mass density of the Universe, one can also calculate the cluster mass by summing the convergence in an aperture.

3.2 Strong lensing cosmological constraints

We further write down the reduced deflection angle in the lens equation equation (3) as

wherein c is the speed of light, and ϕ denotes the projected Newtonian lens potential. The angular positions of multiple images on the source and lens plane are interconnected through the lens potential and cosmic distances D_OL, D_LS, and D_OS. The angular diameter distance between redshift z₁ and z₂ in a flat Friedmann–Robertson–Walker cosmology is given by:

where |${\rm H}_{0} = 100\ h\ \rm km\ s^{-1}$| is the Hubble constant at present day, Ω_m and Ω_X are the corresponding densities of total matter and dark energy, and w_X is the parameter for the equation of state of dark energy.

The above equations demonstrate the capacity of cluster strong lenses to constrain the cosmological parameters {Ω_m, w_X}, although the dependence of cosmology is entangled with the cluster mass. In order to disentangle the degeneracy between cluster mass distribution and cosmology, more than two families of multiple images have to be provided, via the family ratio Ξ:

where z_L, |$z_{S_{1}}$|⁠, and |$z_{S_{2}}$| are the lens redshift and two distinct source redshifts, D represents the corresponding angular diameter distances. Alternatively, for the strong lensing system with only one family of multiple images, as proposed in Golse, Kneib & Soucail (2002), an additional prior on the mass distribution obtained from observations independent of strong lensing should be introduced to break the degeneracy between cluster mass and cosmology.

3.3 Lens modelling

In the context of providing predictions for the impact of point-like multiple images on future surveys, we exclude other systematic errors, and therefore adopt the same form of mass distribution as that used in our simulated clusters during the modelling process. To perform lens modelling, we utilize the Bayesian Markov Chain Monte Carlo (MCMC) sampler, which is implemented in the publicly available lenstool software (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009). This modelling package allows us to sample the parameters related to both the cluster mass distribution and cosmology. As proposed in Kneib & Natarajan (2011), optimizing the model on the source plane can introduce biases due to the magnification factor, which tends to favour flat density profiles and large ellipticities. Hence, to improve parameter accuracy, we choose to perform the optimization on the image plane throughout the work, even though it requires more computational time since it involves an extra step of relensing, and models that generate different image configurations will be rejected (Kneib & Natarajan 2011). The calculation of |$\chi _{i}^{2}$| for a single multiple image family i on the image plane is defined as follows:

where the position of the observed image j|$\boldsymbol{\theta } _{\rm obs}^{j}$| is a combination of the position of the simulated image with a random positional error of |$\rm 0.1^{\prime \prime }$|⁠, |$\boldsymbol{\theta }^{j}({\rm p})$| is the position of image j predicted by the model with parameter set {p}, σ_ij denotes the positional error of this image, and the overall χ² is obtained by summing over all image systems. We assume a smaller global Gaussian error of σ_ij = 0.1 arcsec for the position of each multiple image in space-based measurements (D’Aloisio & Natarajan 2011). Furthermore, we make the assumption that all available multiple images can be measured with spectroscopic redshifts, thus we do not treat the redshifts of multiple images as free parameters throughout the entire process.

3.4 The pixelated modelling algorithm

Given the significant computational complexity involved in modelling the observed multiple images pixel by pixel, current approaches often simplify the analysis by treating the multiple images as point sources. As depicted in Fig. 2, the likelihood is uniform in each direction under the assumption that the multiple image is treated as a point source. However, this assumption is only approximate due to the arc-like shapes produced by the lensing potential, resulting in non-uniform likelihood in each direction. As a result, defining the arc-like shape error solely based on a Gaussian function can lead to biased constraints.

The bias introduced by the point-like multiple image assumption can be buried under various unresolved systematic errors arising from the limitations of observations. None the less, as we look forward to the surveys from the next generation telescopes, such as the Stage-IV surveys, strong lens modelling of galaxy clusters becomes increasingly critical for a deeper understanding of the Universe. In this context, it is imperative that we thoroughly comprehend and address all types of systematic errors. To mitigate such bias, this work introduces an alternative approach: pixelated lens modelling.

The steps in our method closely follow those of lenstool, but with a key difference in the calculation of χ². Instead of comparing the observed and model-predicted positions of multiple images, we compute the corresponding extended surface brightness of every image and sum them to create a total image. The comparison between models and observations is then made on a pixel-by-pixel basis,

where |$S_{\rm obs}^{i}$| is the observed surface brightness on pixel i of the total image, Sⁱ(p) is the model-predicted surface brightness on pixel i with model parameters p, and σ_i is the error on pixel i. Similar to the implementation in lenstool, the observed surface brightness is modelled as the sum of the ideal brightness distribution and a noise field. The noise map is tailored for a JWST setup, incorporating a readout noise of 15.77 e⁻, a zero-point magnitude of 28.0, and a total exposure time of 7537 s across nine exposures. This setup is designed to simulate the noise contributions from both the sky and read noise within the aperture, while we assume that the non-linear effects, for example, charge transfer efficiency, hot pixels, residuals in flat fields, are mitigated before the modelling process. The calculation of χ² is performed by integrating over all the pixels in the entire field of 72 arcsec × 72 arcsec, with a pixel size of 0.18 arcsec. We note that this resolution is currently not at the JWST level, which is considered due to computational costs. As introduced in Section 2.2, to reduce complexity, we assume that the light distribution of each source follows a same S|$\rm \acute{e}$|rsic profile, and we keep these parameters fixed while modelling. The probability distribution function (PDF) of the lens parameters is optimized using the MCMC algorithm implemented in the emcee (Foreman-Mackey et al. 2013) package.

4 RESULTS

Before analysing the impact of assuming multiple images as point sources on the recovery of cluster mass profiles, magnification maps, and cosmography, we first visualize the parameter PDFs of the simulated cluster (A2744-like) from lenstool (blue) and pixelated modelling (pink) in Fig. 3. To ensure a comprehensive comparison, we also include the modelling results from lenstool with emcee as the sampler, referred to as ‘lenstool-emcee’ (green). This additional analysis aims to address the possibility that the investigated bias is induced by the choice of sampling algorithm. Note that the cluster model is simplified throughout the entire work. Specifically, we vary the parameters of the first main halo and cosmology while keeping all other parameters fixed. Additionally, the number of image families used as constraints is reduced to three. These reductions in the number of constraints and free parameters during optimization are intended to ensure computational efficiency. In our analysis, we generate the posteriors by running 100 walkers for 10 000 steps, while discarding the first 3000 steps as part of the ‘burn-in’ phase.

Our analysis reveals that parameters obtained with lenstool and lenstool-emcee exhibit similar degeneracies, the constraint on the main halo’s truncation radius |$r_{\mathrm{ cut}}$| with lenstool-emcee is even less accurate compared to that with lenstool. Therefore, it is reasonable to exclude the possibility that the observed bias is solely due to the sampling algorithm. On the other hand, we find a significant improvement in the posteriors obtained through our pixelated modelling approach. The distribution of each parameter from lenstool is deviated from a Gaussian distribution centred on the input value, as indicated by the grey dashed line. In contrast, such a deviation is not apparent in the parameters derived from pixelated modelling, although we notice a slight deviation from a Gaussian distribution for the constraint on |$r_{\mathrm{ cut}}$| at large values. This can be understood as constraints from strong lensing primarily lying in the central region of a cluster field. In addition, the parameter distributions are narrower when comparing pixelated modelling to traditional modelling. This is attributed to the fact that pixelated modelling employs information from every pixel in the image, as opposed to relying on a handful of multiple image positions. However, cautions should be taken that since we are in the ideal case, discussing the accuracy of lens modelling under the assumption of point-like multiple images, it is worth noting that this analysis does not take into account any sources of systematic errors from the observational side. Therefore, the exact level of precision can only be determined by conducting more realistic simulations.

4.1 Mass profiles

Investigating the mass distribution of strong lensing clusters is fundamental for understanding the nature of dark matter and the formation and evolution of galaxy clusters (Newman et al. 2013; Okabe et al. 2020). As a way to compare the results from traditional and pixelated modelling methods, we present the recovered mass profiles in Fig. 4. In this figure, we plot the enclosed mass obtained using the two modelling methods (blue-shaded region: lenstool, red: pixelated modelling) as a function of its radius R, where the input mass distribution is compared with the masses recovered with 2000 randomly selected realizations.

We find that both methods can provide accurate constraints on the cluster mass profile, given that the same form of cluster mass model with the simulation is used for fitting. This suggests that the investigated bias has little impact on the recovery of a strong lensing cluster mass profile. We also observe that the results from pixelated modelling are more precise compared to traditional modelling. This is expected because pixelated modelling incorporates more information by using the intensity of each pixel as a constraint, allowing for a larger degree of freedom.

4.2 Magnification maps

The magnification of strong lensing clusters has become a primary scientific objective in projects such as JWST (Atek et al. 2023; Bradač et al. 2024; Bradley et al. 2023), HFF (Lotz et al. 2017; Bouwens et al. 2022; Fox et al. 2022; Yang et al. 2022), and Reionization Lens- ing Cluster Survey (RELICS, Coe et al. 2019; Salmon et al. 2020; Neufeld et al. 2022), providing valuable opportunities for observing and studying high-redshift objects. Therefore, it is essential to investigate the robustness of the recovered magnification maps.

Before comparing the results between the traditional method, that is, lenstool, and the pixelated method, we need to establish a definition for the best-fitting lens model based on the output from MCMC sampling. In some literature (Acebron et al. 2019; Remolina González et al. 2021), it is defined as the model parameters that maximize the PDF, often referred to as the maximum a posteriori (MAP) estimate. In lenstool, this is represented as the parameter set with minimum χ². However, the MAP estimate may not always be the best summary of the posterior estimation. In some cases, the best-fitting model is defined using the mean, median, or mode of the samples (Caminha et al. 2016; Raney et al. 2020). In this work, we present the results with both the best-fitting model designated as the one with the minimum χ², and as the medians of the samples. The magnification maps illustrated in this subsection are obtained using the best-fitting models based on equation (6), where we assume a source redshift of z_s = 10. This redshift is chosen as it represents a typical value for discovering and studying high-redshift lensed objects. We utilize the absolute value of magnification throughout the analysis.

Fig. 5 provides an illustration of the difference between the best-fitting (medians) and the fiducial magnification maps obtained with lenstool (left) and pixelated modelling (right) at the core of the cluster field, where critical curves are also superimposed as black dots for reference. To enhance clarity, the difference is normalized by the fiducial magnification map. In most regions of the two panels, the magnification difference is below 10⁻². However, it is worth noting that in regions close to critical curves in the left panel, the difference can exceed 10⁴. These regions near critical curves are of particular interest as they exhibit extreme magnification (up to thousands), that are expected to host individual high-redshift stars (Kelly et al. 2022; Meena et al. 2023), or faint galaxies from the epoch of reionization (Yue et al. 2014; Yan et al. 2023), offering valuable insights into various outstanding questions such as the galaxy formation. Even in the cluster core regions that are not close to the critical curves, there are also large areas with a magnification difference |μ_best − μ_fid|/μ_fid > 10². Therefore, the substantial disparity observed between the recovered and fiducial magnification maps raises concerns about the potential misinterpretation of the properties of lensed objects in these critical regions. On the other hand, we notice that pixelated modelling offers improved accuracy in predicting magnification maps, as shown in the right panel, especially in regions near critical curves. This helps reduce the likelihood of extreme errors that are more prevalent in traditional modelling using lenstool.

We calculate the total area representing different magnification differences, denoted as δμ = |μ_best − μ_fid|, within the four cases: using lenstool or pixelated modelling as the fitting algorithm, and preferring minimum χ² or medians as the best fit. These results are shown in Fig. 6. For both criteria of best fit, we find that lenstool fails to produce magnification maps that are as accurate as those obtained by pixelated modelling.

We proceed with a detailed pixel-by-pixel comparison of the magnification maps. To generate the conditional distribution P(⁠|$\rm \mu _{X}|\mu _{ref}$|⁠) shown in Fig. 7, we select all the pixels with a specific magnification value μ_fid in the fiducial map (e.g. μ_fid = 1) and calculate the median value of these corresponding pixels μ_best in the recovered map. This process is repeated for magnification values spanning the range from μ_fid = 0 to > 1000. Consistent with the observations from Figs 5 and 6, we find that pixels with higher magnification, particularly those near critical curves, exhibit inaccuracies regardless of whether the minimum χ² or median values are used as the best-fitting model. These discrepancies suggest that the properties of lensed sources identified in these regions should be regarded with caution, as they may not be reliable. The biases in the lenstool recovered maps become more pronounced as the magnification μ increases. At μ = 10, the biases are |$\rm |\mu _{best}-\mu _{fid}|/\mu _{fid} = 0.005$| (lenstool) and |$\rm |\mu _{best}-\mu _{fid}|/\mu _{fid} = 0.0007$| (pixelated modelling), while the biases reach |$\rm |\mu _{best}-\mu _{fid}|/\mu _{fid} = 0.4624$| (lenstool) and |$\rm |\mu _{best}-\mu _{fid}|/\mu _{fid} = 0.0009$| (pixelated modelling) at μ = 1000, indicating that only with pixelated modelling can we obtain a magnification with accuracy.

4.3 Cosmological constraints

Since it was initially proposed by Link & Pierce (1998), cluster strong lensing has been extensively studied to demonstrate its feasibility on cosmography, as well as to explore related statistical and systematic errors (Gilmore & Natarajan 2009; D’Aloisio & Natarajan 2011; Acebron et al. 2017; Magaña et al. 2018). In this subsection, we analyse the cosmological constraints from our simulated lens systems, and investigate the impact of the assumption of point-like multiple images on this issue.

We present the cosmological constraints of the four clusters in the simulation, obtained from both modelling algorithms, in Fig. 8. To derive the PDFs P_i(cosmo) for the cosmological parameters from each cluster i, we marginalize over the mass parameters within the corresponding MCMC chains. We observe that despite the absence of systematic errors in the modelling process, the posterior from each individual system is not perfectly aligned with the fiducial cosmology (marked as the black dashed lines). The median values of Ω_m are more likely to be larger than the fiducial |$\rm \Omega _{m, fid} = 0.3$|⁠, while the median values of w lie in the region of lower values compared to the fiducial w_fid = −1.0.

Given the complexity involved in the process of cluster strong lens modelling, a stacked cosmological analysis has been proposed to enhance the figure of merit as well as mitigate potential systematic errors among different clusters. Therefore, we also incorporate the combining constraint from the four clusters. The combined PDF is derived by taking the product of PDFs from all the four systems,

In the left panel obtained with lenstool, the cosmological constraint arising from a stack of four HFF-like cluster lens systems (ΔΩ_m = 0.04 and Δw = 0.16 in the ideal case) demonstrates its competitiveness with other available cosmological probes (cosmic microwave background, CMB, ΔΩ_m = 0.04, Δw = 0.27; Planck Collaboration VI 2020). Combining the systems seems to help eliminate the deviation from the fiducial values within errors, we obtain the cosmological constraints Ω_m = 0.28 ± 0.04 and w = −0.97 ± 0.16, exhibiting slight deviations but are consistent with the fiducial cosmology within the 1σ confidence level. However, we argue that this biased cosmological constraint is an inevitable issue when the Gaussian light profile assumption is used in modelling methods for multiple images. It is important to recognize that the direction of deviation from the truth in each lens system is not linear. Therefore, the cosmological constraints by combining several systems can be even more biased according to the selection bias of strong lensing clusters.

In turn, the parameter PDFs of each lens system follow a more typical Gaussian distribution, with medians that closely align with the input values when using pixelated modelling. This correspondence can be attributed to the comprehensive consideration of the intensity of every pixel in the pixelated modelling approach, eliminating the need for the Gaussian light profile approximation for multiple images. By removing this source of systematic error, we can obtain cosmological constraints with median values of Ω_m = 0.30 and w = −1.0 that are expected to be unbiased.

Based on the above analyses, we emphasize the importance of developing lens modelling algorithms that operate on pixel scale. Such algorithms have the potential to alleviate the deviation from the truth, and provide more accurate and reliable cosmological constraints. By incorporating more information, specifically treating each pixel as a constraint, we can obtain more stringent cosmological contours compared to traditional algorithms that rely solely on the positions of images. Consequently, even with a limited number of observed clusters, we can provide competitive constraints through cluster strong lensing.

5 DISCUSSION AND CONCLUSION

With the designed large field of view and deep imaging capabilities, stage-IV surveys, such as Euclid, China Survey Space Telescope (CSST, 2011), and JWST, will have the ability to observe numerous multiple images and giant arcs in cluster fields. This presents a valuable opportunity to gain deeper insights into systematic errors from the observational side and holds great promise in the cluster strong lensing field. In this context, the assumption of multiple images to be point sources and approximating their light profiles as a two-dimensional Gaussian distribution in traditional modelling, buried in other systematic errors, will emerge as a non-negligible issue. In this work, we employ a series of simulated cluster strong lenses to assess the impact of the above bias and propose a possible solution to address it.

In the first part, we take advantage of the cluster lens system similar to A2744 at redshift z = 0.4 to investigate the mass map obtained from lens modelling, as illustrated in Fig. 4. We find that the cluster mass profile can be accurately recovered using lenstool, implying that the bias in discussion has little influence. Regarding the recovery of the magnification map, as illustrated in Figs 5–7, we observe that the magnification map of the targeted cluster can be measured accurately in most regions of the strong lensing field, if valuable constraints from the data are provided. However, in regions close to critical curves, we find that the magnification difference is still too large to robustly study the properties of high-redshift objects. To address this issue, we exploit a pixelated approach to model the lens system, which takes into account the light distribution on each pixel of the image. We compare the performance of the pixelated approach to traditional lens modelling and demonstrate the improvement achieved with the new approach. Next, we focus on the cosmological constraints obtained from cluster strong lenses. As shown in Fig. 8, we find that strong lensing clusters can provide powerful constraints compared to other tools like CMB (Komatsu et al. 2011; Hinshaw et al. 2013; Planck Collaboration VI 2020), baryonic acoustic oscillation (Alam et al. 2017; Abbott et al. 2022), weak gravitational lensing (Joudaki et al. 2018; Troxel et al. 2018; Bocquet et al. 2019), and others. None the less, the bias is also introduced into the cosmological constraint thanks to the Gaussian-like approximation of multiple images. The pixelated approach also demonstrates enhancement in the same figure, where the peak of parameter PDF align with the true value. Moreover, since the pixelated approach utilizes more information, we can expect competing cosmological constraints even with a single cluster lens system.

Accurate lens modelling requires careful treatments. Our results have demonstrated that the Gaussian-approximated likelihood employed in traditional modelling methods introduces bias to the modelling results, particularly in the recovered magnification maps and cosmological constraints. Furthermore, additional caution should be exercised during the modelling process. For instance, the proper way for model optimization is completed on the image plane, which compares each observed multiple image to its corresponding model-predicted one iteratively and is quite time-consuming. In order to improve efficiency, some studies (Gilmore & Natarajan 2009; Jullo et al. 2010; Acebron et al. 2017) have opted to perform optimization on the source plane by minimizing the discrepancy among model-predicted source positions without the need for an additional step of re-lensing. However, this alternative approach can lead to a biased optimization, favouring mass models with flat density profiles and large ellipticities (Kneib & Natarajan 2011). Based on these considerations, we maintain our approach of fitting the clusters on the image plane throughout the paper, to ensure the attainment of the most accurate results.

The accurate modelling of strong lensing clusters has always posed a challenging task due to the complex structures of the clusters and the large number of strongly lensed multiple images that serve as constraints. In this paper, we propose the idea of creating a new modelling approach to solve the bias induced in traditional modelling. The approach is coded on the basis of lenstool, but instead of computing the χ² and corresponding likelihood function using the positions of multiple images, we prefer to fit the image plane pixel by pixel. Even with the currently available algorithms, the process of fully exploring the parameter space and obtaining results from MCMC or other sampling methods can be time-consuming, often taking days to weeks to fit on the image plane. In our work, modelling the simulated A2744-like cluster with lenstool requires 9 h and 48 min, whereas the same system with the proposed pixelated method takes 25 h and 34 min. Furthermore, considering all image families would result in a modelling time of months for a cluster lens. Hence, any opportunity to accelerate this process is necessary. To address this challenge, one potential solution is to upgrade the algorithm using GPU-based acceleration. Promising tools for this upgrade include jax (Bradbury et al. 2018) and numpyro (Bingham et al. 2018), which can expedite both the processes of computing the deflection field and the sampling. These novel tools have the potential to significantly enhance speed and have already been applied to multiple purposes in astrophysics (Phan, Pradhan & Jankowiak 2019; Campagne et al. 2023; Pasha & Miller 2023).

Our current work primarily relies on analytical simulations to comprehend the bias due to point-like image approximation and test the feasibility of our proposed modelling approach. However, it must be noted that these simulations may not fully capture the intricate details of real clusters in the Universe. Factors such as structures along the line of sight and the surrounding environment can influence the characteristics of strong lensing signals (D’Aloisio & Natarajan 2011; Bayliss et al. 2014; Kuchner et al. 2017; Li et al. 2019), requiring proper modelling. Additionally, adopting parametric models only might underfit a cluster-scale strong lensing system with the structures mentioned above, so it is imperative to consider a combination of parametric and free-form approaches (Beauchesne et al. 2021), or incorporate local potential correction methods (Vegetti & Koopmans 2009; Galan et al. 2022). Our future research endeavors will be dedicated to addressing the aforestated issues.

ACKNOWLEDGEMENTS

The authors thank the anonymous referee who provided useful suggestions that improved this manuscript. We acknowledge the support by National Key R&D Program of China no. 2022YFF0503403 and the Ministry of Science and Technology of China (grant nos 2020SKA0110100). This work is supported by China Manned Space Project (nos CMS-CSST-2021-A01, CMS-CSST-2021-A04, CMS-CSST-2021-A07, and CMS-CSST-2023-A03). HYS acknowledges the support from NSFC of China under grant 11973070, Key Research Program of Frontier Sciences, CAS, grant no. ZDBS-LY-7013, and Program of Shanghai Academic/Technology Research Leader. RL acknowledges the support of National Nature Science Foundation of China (nos 11988101 and 12022306), CAS Project for Young Scientists in Basic Research (no. YSBR-062). AA has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 101024195 ROSEAU. EJ acknowledges the support of the Centre National d’Etude Spatiale (CNES).

DATA AVAILABILITY

The data underlying this article will be shared on a reasonable request to the authors.

Footnotes

References