CURLING - I. The Influence of Point-like Image Approximation on the Outcomes of Cluster Strong Lens Modeling

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 modeling 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 $\Omega_{\rm 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 modeling. This approach reduces the bias in magnification from 46.2 per cent to 0.09 per cent for $\mu \sim 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.


INTRODUCTION
Galaxy clusters, with total masses of 10 14 ∼ 10 15 M ⊙ , 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 (Liang et al. 2014;Yuan et al. 2012; Mahler ★ E-mail: hyshan@shao.ac.cn † E-mail: nan.li@nao.cas.cnet al. 2018;Caminha et al. 2017).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 (Welch et al. 2022;Mahler et al. 2019;Treu et al. 2015;Mercurio et al. 2021), the characteristics of dark matter (Massey et al. 2015;Harvey et al. 2015;Caminha et al. 2019;Jauzac et al. 2019;Bergamini et al. 2023;He et al. 2020), and enabling precise constraints of cosmological parameters (Gilmore & Natarajan 2009;Jullo et al. 2010;Acebron et al. 2017;Caminha et al. 2021;Grillo et al. 2018Grillo et al. , 2020)).
Accurate lens modeling 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 sub-stantial 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 modeling.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 modeling (D'Aloisio & Natarajan 2011a; 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 (Tully & Fisher 1977;Faber & Jackson 1976), the scatter of these galaxy scaling relations can result in biased modeling outputs (Bergamini et al. 2021;Granata et al. 2022).The accuracy of lens modeling is also sensitive to the measurements of multiple image redshifts, as demonstrated by Acebron et al. (2017) who analyzed 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 modeling 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;Nightingale et al. 2018;Birrer & Amara 2018;He et al. 2023).However, due to the computational limitations, current algorithms for modeling 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 modeling 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 modeling.For instance, by incorporating the data from Multi-Unit Spectroscopic Explorer (MUSE) on the Very Large Telescope (VLT) (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 enables the spectroscopic confirmation of multiple images, thereby mitigating the risks of misidentifying image families and enhancing the robustness of modeling.Consequently, the errors arising from the modeling 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 modeling 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 CDM cosmology (Ω Λ,0 = 0.7, Ω m = 0.3, ℎ = 0.7) with a constant dark energy equation of state parameter  =  0 = −1.

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.

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-fit model parameter sets and the galaxy catalogs1 of MACSJ0416.1-2403(Bergamini et al. 2021), Abell 2744 (Bergamini et al. 2023), MACSJ1206.2-0847,and Abell S1063 (Bergamini et al. 2019), to describe the mass distribution of the four mock clusters.These referenced models, built with highquality and extensive spectro-photometric 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: a) dark matter haloes associated with the large-scale, smooth clumps, b) external clumps representing surrounding substructures within the cluster's strong lensing field, c) the brightest cluster galaxies (BCGs), and d) 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 (Kassiola & Kovner 1993;Elíasdóttir et al. 2007, PIEMD) 1. Model parameters for the four cluster lens systems in the simulation. represents the cluster redshift, N s is the number of lensed sources put in the lens system, and  s,toy denotes the source redshift distribution used in the analyses.Note that in each system, only three of the lensed sources are used for lens modeling due to computational constraints.Within a lensed system, each line provides the parameters for a specific component, where Δ and Δ  are the respective relative position to the cluster centre; e is the ellipticity of the clump calculated as e = 1 -/,  and  are the semi-major and semi-minor axis;  is the position angle; r core , r cut , and  are the core radius, cut radius and velocity dispersion; m ref ,  ref core , r ref cut ,  ref are parameters for computing the scaling relations of member galaxies in Eq.2.
density is expressed as: where  0 represents the central density of the profile,  core and  cut are the core radius and the truncation radius, respectively.Within the regions where  core <  <  cut , the profile behaves as  ∼  −2 , while it transitions to  ∼  −4 in the outer regions where  >  cut .
To define the potential of a PIEMD profile in Lenstool (Jullo et al. 2007), we utilize a parameter set consisting of { , , , , r core , r cut ,  }, where ,  describe the centre of the potential,  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  ≃ 0.2 − 0.5 (Robertson et al. 2020), and the strongest cluster lenses peak at  ≃ 0.4 (Fox et al. 2022), we set the redshifts of the four mock clusters to  = [0.35,0.4, 0.45, 0.5].
We summarize the properties of the mock clusters in Table 1.

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 (Bergamini et al. 2023;Mahler et al. 2023;Zitrin et al. 2013;Caminha et al. 2022), 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 2 arcsec 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 modeling (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 Figure 1.Note that the redshift distribution of HUDF galaxies starts from  = 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  = 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 ′′ within a field of 180 ′′ × 180 ′′ .We note that some of the images are missed due to grid resolution in Lenstool while performing lens modeling.Accordingly, we exclude the sources that are found by Lenstool with less than one image from the catalog.Therefore, the number of lensed sources in each lensing system is slightly different, we summarize the number as N s in Table 1.
In addition to determining the positions of the sources, for the novel modeling 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érsic profile, with a random ellipticity  and position angle .The size and Sésic index of the source are given based on their evolution with redshift (Mowla et al. 2019).Therefore, the final range of semi-major axis is from 0.3 to 0.8 arcsec, and the range of Sésic index varies from 1.0 for sources at higher redshifts to 4.5 for lower redshift sources.The superimposed surface brightness distribution  obs within the field-of-view, summed from each source, then becomes the observable in our pixelated modeling 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.Nonetheless, in the context of comparing traditional lens modeling and pixelated modeling, 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.

METHODOLOGY
In this section, we describe the methods of analysing the simulations generated in Section 2. Firstly, 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 modeling methods, encompassing both the traditional lens modeling and a novel pixelated algorithm proposed in this work.

Lensing theory
The lens equation relates the position of source  at redshift  s to the position of image  at redshift   : where ì  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  = 1 2 ΔΨ = 1 2 (Ψ 11 + Ψ 22 ) as a function of the effective lensing potential Ψ, to delineate the magnification of images, the shear vectors  1 = 1 2 (Ψ 11 − Ψ 22 ) and  2 = Ψ 12 = Ψ 21 , where the complex shear  =  1 +  2 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  = Σ Σ crit , where is the critical surface mass density of the Universe, one can also calculate the cluster mass by summing the convergence in an aperture.

Strong lensing cosmological constraints
We further write down the reduced deflection angle in the lens equation Eq. 3 as mann-Robertson-Walker cosmology is given by: where H 0 = 100 ℎ km s −1 is the Hubble constant at present day, Ω m and Ω X are the corresponding densities of total matter and dark energy, and  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 ,  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   ,   1 , and   2 are the lens redshift and two distinct source redshifts,  represents the corresponding angular diameter distances.
Alternatively, for the strong lensing system with only one family of multiple images, as proposed in Golse et al. (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.

Lens modeling
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 modeling process.To perform lens modeling, 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 modeling 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 favor 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  2  for a single multiple image family  on the image plane is defined as follows: where the position of the observed image  ì   obs is a combination of the position of the simulated image with a random positional error of 0.1 ′′ , ì   (p) is the position of image  predicted by the model with parameter set {p},    denotes the positional error of this image, and the overall  2 is obtained by summing over all image systems.We assume a smaller global Gaussian error of    = 0.1 ′′ for the position of each multiple image in space-based measurements (D'Aloisio & Natarajan 2011a) .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.

The pixelated modeling algorithm
Given the significant computational complexity involved in modeling the observed multiple images pixel by pixel, current approaches often simplify the analysis by treating the multiple images as point sources.As depicted in Figure 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 nonuniform 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.Nonetheless, as we look forward to the surveys from the next generation telescopes, such as the Stage-IV surveys, strong lens modeling 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 modeling.
The steps in our method closely follow those of Lenstool, but with a key difference in the calculation of  2 .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   obs is the observed surface brightness on pixel  of the total image,   (p) is the model-predicted surface brightness on pixel  with model parameters p, and   is the error on pixel .Similar to the implementation in Lenstool, the observed surface brightness is modeled as the sum of the ideal brightness distribution and a noise field.The noise map is tailored for a JWST setup, incorporating a Total area as a function of the difference between the recovered and fiducial magnification maps.Left: comparison between Lenstool (blue) and pixelated modeling (pink) with the best-fit result defined as model with minimum  2 ; right: same with left but with the best-fit result defined as medians of the parameters.
readout noise of 15.77  − , a zero-point magnitude of 28.0, and a total exposure time of 7, 537 seconds 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, e.g., charge transfer efficiency, hot pixels, residuals in flat fields, are mitigated before the modeling process.The calculation of  2 is performed by integrating over all the pixels in the entire field of 72 ′′ × 72 ′′ , with a pixel size of 0.18 ′′ .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érsic profile, and we keep these parameters fixed while modeling.The Probability Distribution Function (PDF) of the lens parameters is optimized using the Monte-Carlo-Markov-Chain (MCMC) algorithm implemented in the emcee (Foreman-Mackey et al. 2013) package.

RESULTS
Before analyzing 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 (Abell 2744-like) from Lenstool (blue) and pixelated modeling (pink) in Figure 3.To ensure a comprehensive comparison, we also include the modeling 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 10000 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 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 modeling 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 modeling, although we notice a slight deviation from a Gaussian distribution for the constraint on r 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 modeling to traditional modeling.This is attributed to the fact that pixelated modeling 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 modeling 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.

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 (Okabe et al. 2020;Newman et al. 2013).As a way to compare the results from traditional and pixelated modeling methods, we present the recovered mass profiles in Figure 4.In this figure, we plot the enclosed mass obtained using the two modeling methods (blue shaded region: Lenstool, red: pixelated modeling) 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 modeling are more precise compared to traditional modeling.This is expected because pixelated modeling incorporates more information by using the intensity of each pixel as a constraint, allowing for a larger degree of freedom.

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. 2023;Bradley et al. 2023), HFF (Lotz et al. 2017;Bouwens et al. 2022;Yang et al. 2022;Fox et al. 2022), and RELICS (Salmon et al. 2020;Coe et al. 2019;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, i.e., Lenstool, and the pixelated method, we need to establish a definition for the best-fit 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 probability distribution function (PDF), often referred to as the Maximum A Posteriori (MAP) estimate.In Lenstool, this is represented as the parameter set with minimum  2 .However, the MAP estimate may not always be the best summary of the posterior estimation.In some cases, the best-fit 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-fit model designated as the one with the minimum  2 , and as the medians of the samples.The magnification maps illustrated in this subsection are obtained using the best-fit models based on Eq. 6, where we assume a source redshift of  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.
Figure 5 provides an illustration of the difference between the best-fit (medians) and the fiducial magnification maps obtained with Lenstool (left) and pixelated modeling (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 −2 .However, it is worth noting that in regions close to critical curves in the left panel, the difference can exceed 10 4 .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 (Meena et al. 2023;Kelly et al. 2022), or faint galaxies from the , where the corresponding values of magnification on pixels with  fid are calculated as  best in the recovered magnification map.The panels and color coding are the same as in Figure 6.epoch of reionization (EoR) (Yan et al. 2023;Yue et al. 2014), 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 2 .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 modeling 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 modeling using Lenstool.
We calculate the total area representing different magnification differences, denoted as  = | best −  fid |, within the four cases: using Lenstool or pixelated modeling as the fitting algorithm, and preferring minimum  2 or medians as the best fit.These results are shown in Figure 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 modeling.
We proceed with a detailed pixel-by-pixel comparison of the magnification maps.To generate the conditional distribution P( X | ref ) shown in Figure 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  fid > 1000.Consistent with the observations from Figure 5 and 6, we find that pixels with higher magnification, particularly those near critical curves, exhibit inaccuracies regardless of whether the minimum  2 or median values are used as the best-fit 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

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 2011b;Acebron et al. 2017;Magaña et al. 2018).In this subsection, we analyze 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 modeling algorithms, in Figure 8.To derive the PDFs P  (cosmo) for the cosmological parameters from each cluster , we marginalize over the mass parameters within the corresponding MCMC chains.We observe that despite the absence of systematic errors in the modeling 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 Ω m,fid = 0.3, while the median values of  lie in the region of lower values compared to the fiducial  fid = −1.0.
Given the complexity involved in the process of cluster strong lens modeling, 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, P total (cosmo) = (ΔΩ m = 0.04, Δ = 0.16 in the ideal case) demonstrates its competitiveness with other available cosmological probes (CMB, ΔΩ m = 0.04, Δ = 0.27, Planck Collaboration et al. ( 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  = −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 modeling 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 modeling.This correspondence can be attributed to the comprehensive consideration of the intensity of every pixel in the pixelated modeling 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  = −1.0 that are expected to be unbiased.
Based on the above analyses, we emphasize the importance of developing lens modeling 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.

DISCUSSION AND CONCLUSION
With the designed large field of view and deep imaging capabilities, stage-IV surveys, such as Euclid, CSST, 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 modeling, 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 Abell 2744 at redshift  = 0.4 to investigate the mass map obtained from lens modeling, as illustrated in Figure 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 Figure 5 to Figure 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 modeling and demonstrate the improvement achieved with the new approach.Next, we focus on the cosmological constraints obtained from cluster strong lenses.As shown in Figure 8, we find that strong lensing clusters can provide powerful constraints compared to other tools like Cosmic Microwave Background (CMB Planck Collaboration et al. 2020; Hinshaw et al. 2013;Komatsu et al. 2011), Baryonic Acoustic Oscillation (BAO Alam et al. 2017;Abbott et al. 2022), weak gravitational lensing (Troxel et al. 2018;Joudaki et al. 2018;Bocquet et al. 2019), and others.Nonetheless, 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 modeling requires careful treatments.Our results have demonstrated that the Gaussian-approximated likelihood employed in traditional modeling methods introduces bias to the modeling results, particularly in the recovered magnification maps and cosmological constraints.Furthermore, additional caution should be exercised during the modeling process.For instance, the proper way for model optimization is completed on the image plane, which compares each observed multiple image to its corresponding modelpredicted one iteratively and is quite time-consuming.In order to improve efficiency, some studies (Jullo et al. 2010;Gilmore & Natarajan 2009;Acebron et al. 2017) have opted to perform optimization on the source plane by minimizing the discrepancy among modelpredicted source positions without the need for an additional step of re-lensing.However, this alternative approach can lead to a biased optimization, favoring 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 modeling 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 modeling approach to solve the bias induced in traditional modeling.The approach is coded on the basis of Lenstool, but instead of computing the  2 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, modeling the simulated A2744-like cluster with Lenstool requires 9 hours and 48 minutes, whereas the same system with the proposed pixelated method takes 25 hours and 34 minutes.Furthermore, considering all image families would result in a modeling 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 (Pasha & Miller 2023;Campagne et al. 2023;Phan et al. 2019).
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 modeling 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 (Li et al. 2019;D'Aloisio & Natarajan 2011a;Bayliss et al. 2014;Kuchner et al. 2017), re-quiring proper modeling.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.

Figure 1 .
Figure1.HUDF source redshift distribution.The redshifts of the lensed sources in our sample are assigned according to the HUDF photometric redshift catalog.We exclude the redshifts below  = 1 due to their low lensing efficiency.

Figure 2 .
Figure 2.An illustration of the different appearances between the arc-like image (rainbow), generated by the great potential of the intervening lens, and the point-like image (black circle) under the assumption of the observed multiple images being compact point sources.

Figure 3 .
Figure 3. Posteriors obtained with Lenstool (blue), pixelated modeling (pink) and lenstool-emcee (green) for the parameters of cluster main halo and cosmology.The plotted contours are 1  and 2  levels.Truths are indicated as the grey dashed lines.

Figure 4 .
Figure 4. Circularly averaged mass profiles of the cluster recovered with blue: Lenstool and red: pixelated modeling.Shaded regions are obtained from the randomly selected 2000 realizations in the MCMC outputs.The input mass profile is indicated as the black line.

Figure 5 .
Figure 5. Log-scale difference between the recovered and fiducial magnification maps at  s = 10 with Lenstool (left panel) and pixelated modeling (right panel), with critical curves overlaid as black dots.Note: Regions with little difference (< 10 −7 ) are plotted as white areas since they become Not-a-Number (NaN) in log scale.

Figure 7 .
Figure 7. Conditional distribution P ( best | fid ) (top) with the error of best-fit recovered magnification map (bottom), where the corresponding values of magnification on pixels with  fid are calculated as  best in the recovered magnification map.The panels and color coding are the same as in Figure 6.

Figure 8 .
Figure8.Cosmological constraints derived from the four simulated clusters (Left: obtained with Lenstool, right: obtained with the pixelated modeling).Each colored point denotes the difference between the median values with 1- errors of the cosmological parameters Ω m and  from an individual cluster system and the fiducial values, the input cosmology is marked with the black dashed lines.The grey star represents the result from the stack of the four clusters.
profile.The PIEMD