© The Authors 2026

1. Introduction

Massive objects, such as galaxies and galaxy clusters, can cause the strong gravitational lensing effect; namely, background sources are observed multiple times. Since this effect depends on the total mass of the lens, it is a powerful tool for dark matter (DM; e.g., Grillo et al. 2015; Schuldt et al. 2019; Meneghetti et al. 2020; Wang et al. 2022) and galaxy evolution (e.g., Annunziatella et al. 2017; Mercurio et al. 2021) studies. Furthermore, thanks to the lensing magnification, it enables the detection and investigation of otherwise too-faint, high-redshift sources (e.g., Vanzella et al. 2021, 2022; Meštrić et al. 2022, 2023; Stiavelli et al. 2023; Morishita et al. 2024).

As proposed by Refsdal (1964), the values of the time delays between the multiple images of a strongly lensed supernova (SN) can be used to measure the value of the Hubble constant, H₀, which describes the current expansion rate of the Universe. This technique recently gained significant attention, given the current discrepancy between the measurements from the cosmic microwave background from Planck, resulting in H₀ = (67.4±0.5) km s⁻¹ Mpc⁻¹ (Planck Collaboration VI 2020), and the Cepheid distance ladder approach used by the SH0ES team, providing H₀ = (73.0±1.0) km s⁻¹ Mpc⁻¹ (Riess et al. 2022). The discrepancy raises the question of unknown systematics (e.g., Freedman & Madore 2024; Riess et al. 2024) or new physics (e.g., Di Valentino et al. 2021), and time-delay cosmography (TDC) could play a major role in its clarification.

The first-discovered spatially resolved gravitationally lensed SN (glSN; Kelly et al. 2015, 2016) was imaged multiple times by the massive galaxy cluster MACS J1149.5+2223 (hereafter MACS1149; see Fig. 1 for the cluster core), which lies at a redshift of z_c = 0.5422 (Grillo et al. 2016). The glSN, called SN Refsdal and located at z_SN = 1.489, was lensed into six multiple images, while one image (called SY) appeared before the discovery, so only four time delays could be measured (Rodney et al. 2021; Kelly et al. 2023a). These time delays, together with the observed positions of several multiple images of different background sources, were used by Kelly et al. (2023b) to measure the value of the Hubble constant, resulting in $H_0=(64.8_{-4.3}^{+4.4})\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$, from the combination of eight different mass models, and $H_0=(66.6_{-3.3}^{+4.1})\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$, by combining their two best-fitting models (Oguri-a^⋆ and Grillo-g from Treu et al. 2016). By directly including the measured time delays in the frozen Grillo-g mass model (Grillo et al. 2016, hereafter G16), and by assuming a more general background cosmological model, Grillo et al. (2024) blindly obtained $H_0=(65.1_{-3.4}^{+3.5})\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$, accounting for both the statistical and systematic uncertainties.

Fig. 1. Color-composite image (F435W for blue, F606W + F814W for green, and F105W + F125W + F140W + F160W for red) of the cluster core with the SN Refsdal and its host. Shown are the cluster members (red triangles for photometrically selected members and red circles for spectroscopically confirmed members) and the multiple image systems (blue squares for SN Refsdal and the host and blue pentagons otherwise). We indicate the regions (in orange) used for the SN host reconstruction and extended image model.

More recently, two other glSNe, SN H0pe, which is strongly lensed by the galaxy cluster PLCK G165.7+67.0 (Frye et al. 2023), and SN Encore (Pierel et al. 2024), located in the same host galaxy as the SN Requiem (Rodney et al. 2021) behind the cluster MACS J0138.0−2155, were discovered. A first estimate of the Hubble constant using the SN H0pe time delays was presented by Pascale et al. (2025), resulting in a (post-blind) value of $H_0=(71.8_{-8.1}^{+9.2})\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$, while the analysis with seven different mass models for the SN Encore system resulted in $H_0=(66.9_{-8.1}^{+11.2})\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$ (Pierel et al. 2026; Suyu et al. 2025).

Although significant effort has been made to find more glSNe and improve the H₀ measurements, the relative uncertainty on the value of H₀ is still relatively high (∼5–6% for SN Refsdal and 12% for SN H0pe) when compared with that on the time delays. The uncertainty on the reconstructed total mass distribution of the lens used for TDC contributes, together with that on the measured time delays, to the total error budget of the H₀ value. In particular, for SN Refsdal, the cluster total mass model uncertainty (∼5%) strongly dominates over the time delay uncertainty (∼1%). This reflects the necessity of improving the cluster’s total mass reconstruction and of better understanding its statistical and systematic uncertainties (e.g., Grillo et al. 2020). So far, all the considered cluster total mass models for glSNe rely on the identified multiple image positions of different background sources (e.g., Treu et al. 2016; Kelly et al. 2023b; Acebron et al. 2025; Ertl et al. 2025; Pascale et al. 2025), which represents the state-of-the-art “point-like” approach on cluster scale. However, on galaxy scale, where so far only multiply-imaged quasars have been exploited for TDC (e.g., Suyu et al. 2010; Wong et al. 2017; Birrer et al. 2019; Chen et al. 2019; Shajib et al. 2020; Wong et al. 2024; Birrer 2025), modeling the extended surface-brightness distribution of the quasar’s host galaxy is common, and this has been shown to provide significantly more information than when just using the quasar’s point-like images. Lately, some remarkable effort has been made in the direction of the extended source reconstructions on group-scale systems (Ding et al. 2025; Bolamperti et al. 2024; Wang et al. 2022), while so far, solely Acebron et al. (2024) has successfully pushed forward this extremely time-consuming approach on a cluster-scale system, specifically on the galaxy cluster SDSS J1029+2623. This cluster contains a strongly lensed quasar with a measured time delay and is therefore well suited for TDC. In this work, we followed a similar approach to Acebron et al. (2024) and relied on the same modeling software, GLEE (Suyu & Halkola 2010; Suyu et al. 2012), while focusing on weighting schemes and the impact on the statistical uncertainty rather than the extracted source properties. In contrast, the approach and the aim of Şengül et al. (2023) are notably different, as the authors focused on distortions of giant arcs to detect smaller perturbers and measure their masses (such as DM components or supermassive black holes), did not constrain the total mass of the cluster using a direct reconstruction of the arc as we did here.

Schuldt et al. (2024, hereafter S24) recently published an enhanced strong lens mass model of MACS1149, incorporating new VLT/MUSE observations and following the state-of-the-art point-like strong lens modeling approach. Exploiting the new spectroscopic data, and thanks to a reanalysis of previous observations, this model includes 17 previously unknown multiple images from six different background sources, extending the covered source redshift range from z = 3.703 (system 14 in G16) to z = 5.983. The availability of many multiple images at various redshifts gives the significant advantage of reducing the mass-sheet degeneracy, which is one of the main sources of systematic uncertainties in gravitational lensing studies. This newly presented point-like model resulted in a root-mean-square (RMS) value of the distances between the observed and model-predicted positions of the multiple images of 0$\stackrel{″}{.}$39, with only 34 free parameters.

Building upon this point-like mass model, we present the first extended model of the surface brightness of the host galaxy of SN Refsdal, which is so far the best studied glSNe. We incorporated 77 000 Hubble Space Telescope (HST) pixels of the SN host galaxy (i.e., the pixels within the orange contours in Fig. 1), which translates into an increase of the number of observables by two orders of magnitude, while the underlying parametrization of the cluster total mass remains the same as in S24. This approach also allowed us to reconstruct the delensed surface-brightness distribution of the SN host, which is a spiral galaxy at z_SN = 1.49. With fixed cosmology, we carried out a detailed analysis of the effect of this model extension on the statistical errors on the cluster’s total mass and on the model-predicted time delays. Moreover, we tested different weighting schemes between the point-like and extended image contributions to the value of the total χ². These insights directly expanded the work presented by Acebron et al. (2024) adopting a single weighting scheme and will pave the way toward a new generation of cluster mass models that will enable more precise measurements of the value of the Hubble constant and of other quantities such as the magnification factors of background sources. In particular, with the ongoing Euclid (Laureijs et al. 2011) and Legacy Survey of Space and Time at the Rubin Observatory (Ivezic et al. 2008) observational campaigns, the number of detected glSNe will significantly increase (e.g., Arendse et al. 2024), and complete lensing analyses for TDC will become even more crucial.

The paper is organized as follows. We introduce the exploited observables and the parametrization of the strong lens mass model in Sect. 2, along with a brief discussion of the optimization procedure for an extended source. In Sect. 3 we show our novel results and compare them with those from our previous point-like model from S24. Section 4 summarizes the impact of the new modeling on the derived quantities, such as the total mass of the cluster and the predicted time delays. We conclude our work in Sect. 5.

Throughout this paper, we assume a flat Λ cold dark matter (ΛCDM) cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_M = 1 − Ω_DE = 0.3. In this cosmology, 1′′ corresponds to 6.36 kpc at the cluster redshift, z_c = 0.5422.

2. Lens mass model

The presented mass model was developed with the GLEE (Suyu & Halkola 2010; Suyu et al. 2012) modeling software, while we also made use of glee_tools.py (Schuldt et al. 2023). GLEE allowed us to obtain the total mass distribution of the cluster with a parametric approach, which we describe further in Sect. 2.1. The values and errors of the model parameters are typically constrained by the identified multiple image positions of the strongly lensed sources (see Sect. 2.2), but GLEE also allows one to exploit the extended images of lensed sources in the mass model as constraints (e.g., Suyu & Halkola 2010; Acebron et al. 2024). This additionally enables the reconstruction of the (unlensed) surface-brightness distribution of those sources by minimizing the difference between the model-predicted and observed flux values on the regions covered by the multiple images. This procedure is detailed in Sect. 2.3.

2.1. Mass parametrization

The presented extended-image mass model is based on the point-like mass model of S24 and assumes the exact same cluster total mass parametrization to allow a fair comparison to the model from S24, considered as reference mass model. In short, the overall cluster DM distribution is described by three pseudo-isothermal elliptical mass distribution (Kassiola & Kovner 1993) profiles, two of which are centered in the cluster core. Each of them is characterized by six parameters: the coordinates of the center, x_H and y_H; the axis ratio, q_H; its orientation, ϕ_H; the halo strength (or Einstein radius), θ_H; and the core radius, r_c,H. Moreover, we include, as in S24, a fourth circular (q_H ≡ 1) halo associated with the galaxy-group halo located in the north of the cluster. These cluster halos describe mostly the DM component of the galaxy cluster, which hereafter we refer to as the DM halo.

Furthermore, we included the 308 securely identified cluster members from S24, of which 195 (63%) are spectroscopically confirmed. Following G16, the mass distribution of two cluster members located close to the SN Refsdal host are described by a dual pseudo-isothermal elliptical (dPIE) mass distribution (Elíasdóttir et al. 2007; Suyu & Halkola 2010) with vanishing core radius. All of the remaining 306 cluster members are represented through a circular (q ≡ 1) dPIE profile, whose truncation and Einstein radii are scaled according to their magnitudes in the HST F160W band. While these assumptions result in an approximation of the cluster’s total mass distribution, these mass profiles and scaling relations are commonly adopted for cluster mass models. The cluster member distribution in the cluster core is shown in Fig. 1 with red triangles if photometrically identified, and red circles if spectroscopically confirmed. The distribution of all cluster members is displayed in Fig. 1 of S24.

2.2. Multiple image positions

While the mass parametrization (see Sect. 2.1) was kept the same, we considered four different combinations and weightings of the observables, as summarized in Table 1. In the first two models, hereafter referred to as Sw and Ew, we included all multiple images from S24, namely 63 multiple images belonging to 18 knots of the SN Refsdal host (z_SN = 1.489) and 43 multiple images from 16 other strongly lensed sources, covering the redshift range between 1.24 and 5.98. In the other two models, hereafter referred to as Ewo and Swo, we excluded the identified multiple-image positions of the SN Refsdal and of its host¹. The observed positions of the SN Refsdal and the identified multiple images of the SN host are shown in Fig. 1 (blue squares). We refer to G16 and S24 for their exact coordinates and redshifts.

Since the cluster’s total mass distribution is relatively complex and the image multiplicity high, we adopted the expression of S24 for the image-position χ² term, namely

$$\begin{array}{c}\hfill {\chi }_{}^2\text{img}={\displaystyle \sum _{j=1}^{N^{\text{fam}}}}\{\begin{array}{cc}{\displaystyle {\sum }_{k=1}^{N_j^{\text{obs}}}}\frac{{(x_{j,k}^{\text{obs}}-x_{j,k}^{\text{pred}})}^2}{\mathrm{\Delta }{x}_{j,k}^2}+\frac{{(y_{j,k}^{\text{obs}}-y_{j,k}^{\text{pred}})}^2}{\mathrm{\Delta }{y}_{j,k}^2}& \text{if}{N}_j^{\text{obs}}\le N_j^{\text{pred}}\\ \infty & \text{otherwise}\end{array},\end{array}$$(1)

which ensures the correct image multiplicity for all sources. The mass model relying only on the point-like image positions was obtained by minimizing this term, and keeping it here allowed us to draw a direct comparison to the mass model presented by S24. In this equation, N^fam describes the number of families, N_j^obs and N_j^pred the number of observed and predicted image positions of family j, respectively; x^obs and y^obs the observed x and y coordinates, respectively; and x^pred and y^pred are the predicted x and y coordinates, while Δ describes the image position uncertainty of the corresponding image position. In the next section, we provide further details on the parameter optimization.

2.3. Extended image model and source surface-brightness reconstruction

In addition to the image positions described in Sect. 2.2, we included the flux values of the pixels of the SN host in our mass model as observables. Since we only wanted to model and reconstruct the surface-brightness distribution of the SN host, we followed the approach of G16 and created a difference image between the HST filters F435W and F606W. In detail, we used a representative subset of spectroscopically confirmed cluster members to determine a scaling factor for the F606W band to best subtract the light of these cluster members, resulting in an image with only the (lensed) background sources. The best aperture size to estimate the scaling factor was determined in such a way that the remaining fluxes within 1″ around the selected cluster members were minimal.

To create error maps, which were not directly available for the HST images, we approximated them for both filters using glee_tools.py (Schuldt et al. 2023). In this approximation, the error for pixel j of filter k is defined as

$$\begin{array}{c}\hfill {\sigma }_{\text{tot,}k,j}^2=\{\begin{array}{cc}{\sigma }_{\text{poisson,}j}^2+{\sigma }_{\text{bkgr,}j}^2& \text{if}{\sigma }_{\text{bkgr,}j}\le {\sigma }_{\text{poisson,}j}\hfill \\ {\sigma }_{\text{bkgr,}j}^2& \text{otherwise}\hfill \end{array},\end{array}$$(2)

with the Poisson noise (Hasinoff 2012)

$$\begin{array}{c}\hfill {\sigma }_{\text{poisson,}j}=\sqrt{\frac{|d_j|}{t_j}},\end{array}$$(3)

which we obtained from the observed pixel intensity, d_j (in e⁻-counts per second) and the exposure time, t_j, of pixel j, and σ_bkgr,j approximated as a constant value equal to the standard deviation from an empty region. After generating both error maps separately, we combined them via

$$\begin{array}{c}\hfill {\sigma }_{\text{tot,}j}=\sqrt{{\sigma }_{\text{tot, F435W,}j}^2+{\sigma }_{\text{tot, F606W,}j}^2},\end{array}$$(4)

without different weightings of the individual error maps, as we applied a weighting to the final error map (see Table 1).

In addition, we needed to define which pixels of the difference image belong to the SN host and should be taken into account during the modeling. For this, we created an arcmask (see, e.g., Ertl et al. 2023, for further details). The considered area is shown with orange contours in Fig. 1, and contains N_pix = 77 000 pixels with a pixel size of 0$\stackrel{″}{.}$03.

In each iteration of the parameter optimization, all pixels within the arcmask are mapped from the observed image plane onto the source plane, using the given mass model through the ray-tracing equation. In this way, the source was reconstructed on a grid of pixels, in our case with 100 × 100 pixels. Since we usually reconstruct astrophysical sources, which typically show quite smooth surface-brightness distributions, a regularization term was applied to penalize very irregular surface-brightness reconstructions. The regularization strength was adapted for each source, which also allowed the reconstruction of relatively clumpy spiral galaxies as the SN host. Specifically, we adopted the curvature form of the regularization term implemented in GLEE, and refer the reader to Suyu et al. (2006) for further details on the definition and procedure.

Furthermore, in each optimization step, the pixel grid of the reconstructed source was mapped back into the image plane, predicting the lensed SN host’s intensity distribution. By comparing, for each pixel, the predicted flux values, d^pred, a vector of length N_pix, with the observed ones, d^obs, a vector of the same length, we were able to define a chi-square term as

$$\begin{array}{c}\hfill {\chi }_{}^2\text{ext}={({\mathit{d}}^{\text{obs}}-{\mathit{d}}^{\text{pred}})}^T{C}_{\text{D}}^{-1}({\mathit{d}}^{\text{obs}}-{\mathit{d}}^{\text{pred}}),\end{array}$$(5)

where C_D is the image covariance matrix described in Suyu et al. (2006). This led to a very significant increase of the runtime and of the required random access memory for the matrix multiplication.

Assigning smaller uncertainties to the knots of the SN host instead of a full extended-image mass model as described here is possible and leads to a better reconstruction of the host. However, the extended image model has several further advantages. First, lowering the uncertainties on the multiple images works only for extended sources with significant luminosity structures, such as the SN host galaxy we have here, but not for those with relatively smooth surface-brightness distributions, such as the quasar host reconstructed by Acebron et al. (2024), as one cannot identify many knots there. Second, all these knots – 63 multiple images from 18 different knots in MACS1149 – need to be properly matched and located, while the pixels used for an extended reconstruction can easily be defined. Third, the lowered uncertainties would only take the identified positions into account, while with our reconstruction we simultaneously made use of the pixel intensities, which have to match each other on a large number of pixels, affecting the deflection angle and magnification values in that region. Fourth, it enables studies of the background source and allows the measurement of extracted quantities such as the half-light radius of the reconstructed source (Acebron et al. 2024).

Furthermore, this approach can easily be extended to simultaneously model multiband observations of the same source, or to sources at different redshifts. For instance, Schuldt et al. (2019) presented mass models of the galaxy-scale system called the Cosmic Horseshoe (Belokurov et al. 2007) by modeling the main arc in multiple filters and the radial arc simultaneously.

We then optimized the mass parameter values using simulated annealing and Markov chain Monte Carlo techniques (Dunkley et al. 2005; Foreman-Mackey et al. 2013), as well as emcee (Foreman-Mackey et al. 2013), by minimizing the total chi-square value, χ² = χ_img² + χ_ext². However, since the extended image provided almost three orders of magnitude more observables at a specific redshift in our case (the SN host redshift z_SN = 1.489) than the point-like images, we tested different weightings of the two chi-square terms. This was achieved by a constant scaling of the combined error map affecting χ_ext² and a constant scaling of the image position uncertainties influencing the χ_img² term. We remark that the precise weighting should depend on the specific scientific case to address. Raising the contribution of the multiple images from sources at different redshifts helps to break the mass-sheet degeneracy (Falco et al. 1985; Schneider & Sluse 2013), while more weight on the extended images leads to a better reconstruction of the SN host. Since it is not obvious which weighting is preferred for TDC and its impact on the reconstructed mass model parameters and the quantities derived from them, as for the time delays, we provide an extensive comparison of the different weightings in Sect. 3.

In detail, as already noted in Sect. 2.2, we considered models with (Sw and Ew ) and without (Swo and Ewo) the multiple images of the SN Refsdal and its host, as these image positions were measured from the HST images, similarly to the flux values used in the extended-image model. We also explicitly note that the SN image positions were excluded in the Swo and Ewo models. Furthermore, since we adjusted all the uncertainties to obtain a reduced χ² ∼ 1 for the correct estimate of uncertainties of the model parameters, we considered, in both cases, two different weightings of the χ² terms. For the models Ew and Ewo, we fixed the image-position uncertainties to those from S24, which means that the χ_img² contribution is negligible as long as the image multiplicity is correct. In case of a wrong image multiplicity, χ_img² and consequently χ² are both infinity, and the model is thus rejected, regardless of the extended-image reconstruction (i.e., irrespective of the χ_ext² value). In other words, in these two models we gave full weight to the extended image, which is the preferred choice when one is particularly interested in studying the source properties. Additionally, we considered the models Sw and Swo, where we gave a similar weight (∼50%) to the two χ² terms by adjusting the error map and the image position uncertainties. In these cases, we gave more weight to the sources at different redshifts than in the previous two models, which was crucial for breaking the mass-sheet degeneracy. After obtaining a good fit for each of the four models considered, we ran final emcee chains for 96 days each with ten cores, leading to final chain lengths of more than 600 000 iterations per model. Due to the long runtime, we decided to use an approximate point spread function, that, according to our tests with shorter chains and a point spread function properly constructed from stars in the field, led to the same conclusions overall. The absolute inferred mass parameter values and the details of the source reconstruction (irrelevant in this work) differed only slightly.

3. Modeling results and comparison to the point-like model

As described above, we built four different strong lensing models of the MACS1149 cluster, taking into account the extended surface-brightness distribution of the SN host, and compared them to the state-of-the-art point-like model presented by S24. An overview of all models is presented in Table 1, where we give the corresponding χ² values of the optimized mass models. Since the difference in the models is the weighting of the different χ² terms, we report the values for the image positions (χ_img²) and extended image (χ_ext²) separately. Furthermore, because the associated uncertainties differ for the different models (see Sect. 2), we also report chi-square values calculated adopting the same uncertainties, namely the uncertainties from S24 for the multiple images, resulting in χ_img,S24², and the original error map for the extended image obtained from Eq. (4), resulting in χ_ext,orig², as these values can be properly compared across the five models.

By looking at the χ_ext,orig² values, we find an expected significant improvement in the reconstructed image, when compared to the point-like model, and we can see that the models Ew and Ewo models outperform the Sw and Swo models. Moreover, it is also expected that the χ_img,S24² value increases as we give less weight to this term. However, considering that the Ew and Ewo models have full weight on the extended image and that the multiple image positions of the SN host and other sources at different redshifts contribute a negligible amount, it is noteworthy that these image positions are nonetheless very well reconstructed. This is shown in the RMS values that increased by only 0$\stackrel{″}{.}$05 for the Sw and Swo models and by a factor of approximately two for Ew and Ewo, compared to the point-like model, which has its full weight on the multiple image positions. We note that the reported RMS values were always calculated including the image positions of the SN host (labelled RMS(w) in Table 1) to allow a direct and fair comparison between the different models. We remind the reader again that we kept the mass model parametrization exactly as in S24 to allow a fair comparison to that model, despite the high number of newly introduced observables, so a slightly higher RMS is expected. However, in a future work more flexible mass models could be explored to better reproduce the now incorporated 77,000 HST pixel values. Moreover, we note that an RMS value of 0$\stackrel{″}{.}$73, which is the highest value we obtained with the model flexibility of S24, is comparable to state-of-the-art point-like models from other clusters (e.g., Bergamini et al. 2019; Richard et al. 2021; Furtak et al. 2023), where the full weight is on the multiple image positions.

The improvement in the image reconstruction is visualized in Fig. 2, where we compare Ew, resulting in the lowest χ_ext,orig² value, with the point-like model S24. Although the Ew model still shows some residuals, the reconstructed source has significantly more structure than that from the point-like model. We highlight the refinement going from S24 to Ew in Fig. 3, where we show the corresponding normalized residuals as absolute ratio (left) and difference (right). While the majority of the 77 000 pixels show similar normalized residuals (green), we see a clear trend toward lower values (blue) for the Ew model, corresponding to a better image reconstruction. It is remarkable that the normalized residuals of more than 200 pixels improved by more than 1.5σ compared to the reconstruction with the point-like model, while the normalized residuals for only four pixels increased by 1.5σ. This reflects the power of the extended image model and translates into a refined cluster-mass reconstruction around the SN host, which is crucial for TDC and studies of the physical properties of the source galaxy.

Fig. 2. Surface-brightness reconstruction from model Ew (top row), showing the best reconstruction, and from the point-like model S24 (bottom row), which was not optimized on the extended image. We show, from left to right, the observed image, the model-predicted image, the normalized residuals, and the reconstructed source (the SN host galaxy). We illustrate the sizes of the images and of the reconstructed sources with green bars of 5″and 1″, respectively. We again note that the Ew model is optimized with an approximate point spread function (see Sect. 2.3 for details), which slightly affects the final appearances of the reconstructed image and source.

Fig. 3. Comparison of normalized residuals (NRes) of the model Ew and the point-like model by S24, as an absolute ratio (left) and absolute difference (right). The bottom panels show the 1D histograms, and the top panels show the corresponding 2D image color-coded in the same way as the histograms.

By adding the surface-brightness distribution of the SN host galaxy to the other point-like multiple images, we obtain both a significant improvement in the extended image reconstruction and unprecedented precision in the inferred values of the individual mass model parameters. This is illustrated in Fig. 4, where we show the relative errors of the 34 free parameters for the five models, Ew, Ewo, Sw, Swo, and the point-like model. The statistical uncertainties of all model parameters, which were estimated by extracting 100 random iterations from the final sampling chains, are reduced by between one and two orders of magnitude. Solely the uncertainties on the three varying redshifts of sources 6, 7, and 8 (see G16 and S24) are comparable among models S24, Ew, and Ewo. This can be explained by the fact that the adopted image-position uncertainty is the same for these three models, while for models Sw and Swo it is much lower (see Table 1). As a consequence, a small change of the redshifts in the Sw and Swo models results in a high χ_img² value, such that these models were penalized. On the other hand, the redshift values of models Ew and Ewo can vary more broadly, similarly to those in S24. We note that the relative error is, for some parameters, slightly misleading. For instance, the center of the first DM halo is very close to the origin of the selected coordinate system (see G16 and S24), so its relative error always appears relatively high. However, since this holds for all five models, the comparison between the models, which is the aim of the plot, is fair and clearly demonstrates that we obtain a significant improvement in precision, regardless of the chosen weighting and inclusion or exclusion of the multiple images of the SN host.

Fig. 4. Relative errors of the 34 free parameters (left) and of the magnification factors at the model-predicted positions of the point-like multiple images (right) from the five different models (indicated by different colors and symbols). For the definition of the individual model parameters, see Sects. 2.1 and S24. While the 106 multiple images are, for simplicity, not labeled, we highlight those belonging to the SN and its host (and thus excluded in models Ewo and Swo; white background), those with variable redshift (since no spectroscopic redshift is known; light gray), those belonging to the northern galaxy group (see Fig. 2 of S24; medium gray), and those from additional sources in the cluster core (dark gray). We observe a significant reduction in the statistical uncertainties of the model parameters (left) and magnifications (right) when we include the 77 000 HST pixels as direct constraints in the strong lensing model.

4. Predicted cumulative total mass profile, time delays, magnifications, and multiband source reconstruction

After presenting the improvement on the mass parameter values, we considered the impact on some derived quantities. First, we computed the cumulative total mass profile of the cluster MACS1149, as shown in Fig. 5. We compared our four new extended-image and previous point-like models and found good agreement among them. We obtained a projected total mass integrated within 500 kpc of (6.31±0.07) × 10¹⁴M_⊙, with the 1σ errors approximately one order of magnitude smaller when the extended surface brightness is included in the models. We note that the quoted uncertainties are only the statistical uncertainties, while the systematic uncertainties should be assessed with different cluster’s total mass parametrizations.

Fig. 5. Cumulative projected total mass profile M(< R) of MACS1149 as function of the radius, R, measured from the brightest cluster galaxy (see G16. The plotted bands correspond to the 1σ intervals obtained from 100 random iterations of the final sampling chains. We compared our four new extended image models with the point-like model from S24 and found significantly smaller statistical uncertainties for the extended image models. The positions of the multiple images are marked with small dashes on the lower x-axis and color-coded according to their redshifts. The bar lengths reflect the positional uncertainties (a small dash corresponds to an HST identification, while a long dash corresponds to a MUSE-only identification; see S24 for details). We further show the positions of the 77 000 HST image pixels used in the extended image models with a gray histogram, where the height (in log-scale) corresponds to the number of pixels per radial bin.

Second, we compared the model-predicted values of the time delays between multiple images of the SN observed in MACS1149. Because all our models were obtained with a fixed cosmology (i.e., flat ΛCDM; see Sect. 1), we only considered the four presented models (Ew, Ewo, Sw, and Swo) and the point-like models from S24 and G16, which were built under the same cosmological assumptions. Figure 6 shows the median values and 1σ errors of the four predicted time delays for all six models, obtained from 100 different iterations of the final sampling chains. A comparison with the measured time delays would not be fair, given the fixed values of the cosmological parameters in our lensing models. The updated measurement of the values of the cosmological parameters, including the Hubble constant, by exploiting the extended image models exceeds the scope of this study and is thus left for a future work.

Fig. 6. Median values and 1σ uncertainties of time delays obtained at the model-predicted positions of the multiple images of the SN with fixed (flat ΛCDM) cosmology. The four new extended image models and the two previous point-like models by S24 and G16 are shown with different colors. We notice the significantly lower statistical uncertainties for the extended image models compared to the point-like models. We note again that the extended image models are constructed with an approximate point spread function (see Sect. 2.3 for details), which slightly affects the exact values of the predicted time delays.

Although the predicted values of the point-like models S24 and G16 differ slightly, they are consistent, given the 1σ uncertainties. The uncertainty on Δt_SX:S1, which is the longest time delay and thus the most relevant one for cosmological applications, obtained with the mass model from S24 than with that from G16. This can probably be ascribed to the inclusion of the SX observed position in S24, which is different from G16. The latter mass model was developed before the detection of the SN at the SX position (i.e., this position was a blind prediction of that model). Considering this, the predicted value and error of the Δt_SX:S1 time delay with model G16 represent a remarkable result.

The model-predicted time-delay median values of the extended image models are all well within the 1σ intervals of S24 and G16, and, regardless of the weighting, the obtained statistical uncertainties are more than an order of magnitude smaller. Although a lower uncertainty is expected, given the improved reconstruction of the cluster’s total mass (see Fig. 4), the observed effect is remarkable. The remarkably small error bars are responsible for the fact that only models Ew and Ewo, and Sw and Swo, respectively, are consistent within 1σ. Consequently, the model-predicted time delays mildly depend on the weighting, i.e., whether we reconstruct the SN host’s surface brightness, and thus the local cluster’s total mass distribution, in the best way (models Ew and Ewo), or if we give 50% of the weight to the other multiple images, which are at different redshifts, to reduce the mass-sheet degeneracy. We remind the reader that the quoted uncertainties only refer to the statistical errors, and the small differences among the four models can be used to quantify a first source of systematic uncertainties. Until now, with point-like models, the systematic uncertainties were found to be significantly smaller than the statistical ones (Grillo et al. 2020). Given the increased precision with the source surface-brightness reconstruction, systematic effects will become relevant.

It is noteworthy that the statistical uncertainty on the time-delay value Δt_SX:S1 predicted with point-like models is ∼5%, exceeding the error on the measured time delay, which is on the order of 1.5% (Kelly et al. 2023a). Both these terms contribute with similar weights to the resulting precision on the inferred value of H₀. Thanks to the extended image models, we achieved a sub-percent statistical uncertainty on the model-predicted time delays. This translates into significantly more precise measurements of the values of the cosmological parameters.

Another strong application of cluster lensing, including MACS1149 itself (e.g., Stiavelli et al. 2023; Morishita et al. 2024), is the study of high-redshift sources that are amplified thanks to the lensing magnification effect. Our extended image models are expected to impact this research field in two ways. As already shown in Fig. 2, we obtained a delensed version of the SN host galaxy. However, while we used a difference image from two HST bands to minimize the light contribution from all cluster members for the modeling, we can also use a given extended image model to delens the images in a single band. As a demonstration, we used the final Ew model, achieving the lowest χ_ext,orig², to reconstruct the delensed SN host galaxy from the individual HST images. In detail, we obtained the reconstruction in three HST bands, F814W, F606W, and F435W, and in Fig. 7 we show them as a color image. We modified the mapped pixels slightly, mainly the regions of the four SN images, to exclude the light contribution from nearby cluster members and the SN images (slightly visible in F435W). While the extended image modeling is computationally very expensive, the source reconstruction is not (this is roughly equivalent to a single model iteration), and thus it can be easily done in multiple bands and with images from other telescopes.

Fig. 7. Color composite image of the SN host reconstruction from the Hubble Frontier Field HST bands F814W (red), F606W (green), and F435W (blue), using the Ew model. The mapped region was slightly modified compared to that used in the model optimization (shown in Fig. 2) to minimize light contamination from some cluster members and the SN images. This figure underscores our ability to directly model the SN host in multiple bands. The image has an extension of 2$\stackrel{″}{.}$05 and 2$\stackrel{″}{.}$27 along the x and y axes, which correspond to 17.84 kpc and 19.75 kpc, respectively, at a redshift of zSN = 1.49.

The extended image modeling and the delensing of background sources is only possible for extended ones, such as the SN host with 77 000 HST pixels. This technique is not directly applicable for compact objects, extending only over a small amount of pixels. Therefore, resulting in an additional application scenario, we show the predicted magnifications of the 106 multiple images, which is a crucial measurement for studying background sources, in Fig. 4 (see Grillo et al. 2016, and S24 for the exact coordinates). As we see from Fig. 4 (right panel), where we show the relative uncertainties on the magnification, we found an improvement by ∼1 order of magnitude between the classical point-like model and the four extended image models. As expected, this holds specifically for the multiple images from the SN host (white background in Fig. 4). Therefore, we expect this model approach to impact the high-redshift studies in general. It is understandable that models Ew and Ewo preform better than Sw and Swo on the multiple image positions of the SN host, while the performance on the other images is overall worse. This holds in particular for the images with photometric redshifts (light gray band in Fig. 4), as these images add very little information to Ew and Ewo, while, because of the reduced image position uncertainty, they contribute significantly to models Sw and Swo.

It is also interesting to see the behavior of the northern galaxy group, where we lack a significant improvement for most of the seven multiple images (medium gray background in Fig. 4). This is likely due to the angular distance to the SN host galaxy, where the additional observables are located. In particular, models Ew and Ewo lack any significant constraints in this region, which underscores that the choice of weighting depends on the project’s scope and can be used to mitigate this effect. However, we did not find a strong dependency on the weighting for the overall predicted quantities (see Figs. 5 and 6). Furthermore, we explicitly note that all extracted quantities (total cluster mass, time delays, etc.) fall well within the previous statistical uncertainties from our point-like mass model S24, while the differences between our new mass models reflect systematic uncertainties that became relevant. These systematic uncertainties are well below the statistical uncertainties of the mass model from S24 (see Fig. 6), which is in agreement with previous findings by Treu et al. (2016) and Grillo et al. (2020), for example.

5. Summary and conclusion

We obtained the first total mass model of the strong lensing cluster MACS1149 by combining 106 identified multiple images and the extended surface-brightness distribution of the SN host galaxy. The cluster mass was modeled with parameterized profiles, following the point-like model from S24 to enable a direct comparison, and composed of three elliptical DM halos, a spherical DM halo associated with the galaxy group in the north (see Fig. 2 in S24), two elliptical profiles for two cluster members close to the SN host, and a spherical profile for each of the 306 additional cluster members, where their Einstein and truncation radii were scaled according to given scaling relations. Additionally, three background sources without spectroscopic redshift were included with a variable redshift value. In total, we optimized 34 free parameters (see Fig. 4). While this seems a small number of parameters, including a few hundred mass components and several dozen background sources already leads to a relatively long runtime, when optimizing a model on the image plane. We faced significant computational challenges by incorporating 77 000 HST pixels in the modeling. Furthermore, the random access memory (RAM) increased remarkably, from ∼1GB/core to around 180GB/core, making a parallelized optimization extremely difficult.

Since we considered the flux values of 77 000 HST pixels in addition to the 212 multiple-image positions (x and y coordinates for 106 images) as observables, we examined, for the first time, different combinations of the two groups of observables. Namely, we tested the effect of the weighting with four different models, fixing the model parametrization as in S24. In detail, we analyzed a first version effectively with full weight on the extended image of the SN host, and a second version where we rescaled the image position uncertainties and the surface-brightness error map so that both terms contributed similarly to the total chi-square value. Furthermore, for the two previous cases, we built a version with all identified image positions, and one excluding those of the SN host, as these were considered as already included in the HST pixels. Regardless of the weighting and of the included or excluded SN host image positions, the values of the model parameters were measured with precision between one and two orders of magnitude higher than in the point-like model (see Fig. 4). We also found slightly lower statistical uncertainties for models with equal weights (i.e., models Sw and Swo). As expected, we obtained a better reconstruction of the SN host with the models adopting full weight on the extended images (see Fig. 3 and Table 1).

Going beyond the simple comparison of the parameter values, we quantified the impact on some derived quantities of the models. First, given the powerful use of gravitational lensing to measure the total mass of a lens, we obtained the projected cumulative total mass profile of the lensing cluster MACS1149 (see Fig. 5) with remarkably small errors. Then, since the cluster produced multiple images of SN Refsdal, we predicted the values of their time delays and compared them with those from Grillo et al. (2016) and S24 within a fixed cosmological model. In particular, we did not find any bias on the predicted time delays using the different weighting and showed that the statistical uncertainty, associated with the remaining uncertainty in the cluster’s total mass reconstruction, decreased from ∼5%, for the point-like model, to less than 1% for the extended image models. This will translate into significantly more precise estimates of the values of the cosmological parameters. We thus reached a regime in which the systematic errors need to be quantified, as we can see from the minor differences between the four different mass model predictions. Moreover, we measured the magnification factors for all 106 identified multiple images with the four new mass models and the point-like model S24 and compared their relative errors (see Fig. 4). In all cases, we found a significant improvement in precision when taking the surface-brightness distribution of the SN host into account.

Furthermore, we reconstructed the SN host galaxy in multiple HST bands and show, for illustration purposes, a color-composite image in Fig. 7. Thanks to extended image models, such as those developed in this work, precise source reconstructions will be possible in the future, using the HST, James Webb Space Telescope (JWST), or other observations. This will likely significantly impact high-redshift source studies. While it is clear that for this specific science case the best source reconstruction is preferred (i.e., model Ew or Ewo), the other weighting (i.e., model Sw and Swo) is likely more effective in reducing the mass-sheet degeneracy thanks to the multiple images at different redshifts, and it is thus more suitable for cosmological applications. This shows that the preferred weighting might depend on the specific scientific goal.

In conclusion, we present the potentialities of extended strong lensing models in the Hubble Frontier Field galaxy cluster MACS1149, showing the first multiple-imaged and spatially resolved supernova. We demonstrate the impact on several scientific applications, from TDC to high-redshift studies. With an improvement of more than one order of magnitude in the statistical uncertainties of the model parameters, this work paves the way for a new generation of mass models of cluster lenses exploiting high-resolution imaging data from the HST, JWST, and Euclid telescope. Accelerated lensing codes such as Herculens (Galan et al. 2022), Gravity.jl (Lombardi 2024), or GPU-based GLEE and GLAD (Wang et al. 2025) will play a fundamental role in reducing the computational time and thus making analyses of this kind feasible.

Acknowledgments

We thank the referee for constructive comments that helped to improve the clarity of the presented work. SS has received funding from the European Union’s Horizon 2022 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101105167 – FASTIDIoUS. We acknowledge financial support through grants PRIN-MIUR 2017WSCC32 and 2020SKSTHZ. PB acknowledges financial support through grant PRIN-MIUR 2020SKSTHZ and support from the Italian Space Agency (ASI) through contract “Euclid - Phase E”, INAF Grants “The Big-Data era of cluster lensing” and “Probing Dark Matter and Galaxy Formation in Galaxy Clusters through Strong Gravitational Lensing”. This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. This work uses the following software packages: Astropy (Astropy Collaboration 2013, 2018), Emcee (Foreman-Mackey et al. 2013), GLEE (Suyu & Halkola 2010; Suyu et al. 2012), glee_tools.py (Schuldt et al. 2023), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011; Harris et al. 2020), Python (Van Rossum & Drake 2009), Scipy (Virtanen et al. 2020).

References

All Tables

All Figures

Fig. 1. Color-composite image (F435W for blue, F606W + F814W for green, and F105W + F125W + F140W + F160W for red) of the cluster core with the SN Refsdal and its host. Shown are the cluster members (red triangles for photometrically selected members and red circles for spectroscopically confirmed members) and the multiple image systems (blue squares for SN Refsdal and the host and blue pentagons otherwise). We indicate the regions (in orange) used for the SN host reconstruction and extended image model.

In the text

Fig. 2. Surface-brightness reconstruction from model Ew (top row), showing the best reconstruction, and from the point-like model S24 (bottom row), which was not optimized on the extended image. We show, from left to right, the observed image, the model-predicted image, the normalized residuals, and the reconstructed source (the SN host galaxy). We illustrate the sizes of the images and of the reconstructed sources with green bars of 5″and 1″, respectively. We again note that the Ew model is optimized with an approximate point spread function (see Sect. 2.3 for details), which slightly affects the final appearances of the reconstructed image and source.

In the text

	Fig. 3. Comparison of normalized residuals (NRes) of the model Ew and the point-like model by S24, as an absolute ratio (left) and absolute difference (right). The bottom panels show the 1D histograms, and the top panels show the corresponding 2D image color-coded in the same way as the histograms.
In the text

Fig. 4. Relative errors of the 34 free parameters (left) and of the magnification factors at the model-predicted positions of the point-like multiple images (right) from the five different models (indicated by different colors and symbols). For the definition of the individual model parameters, see Sects. 2.1 and S24. While the 106 multiple images are, for simplicity, not labeled, we highlight those belonging to the SN and its host (and thus excluded in models Ewo and Swo; white background), those with variable redshift (since no spectroscopic redshift is known; light gray), those belonging to the northern galaxy group (see Fig. 2 of S24; medium gray), and those from additional sources in the cluster core (dark gray). We observe a significant reduction in the statistical uncertainties of the model parameters (left) and magnifications (right) when we include the 77 000 HST pixels as direct constraints in the strong lensing model.

In the text

Fig. 5. Cumulative projected total mass profile M(< R) of MACS1149 as function of the radius, R, measured from the brightest cluster galaxy (see G16. The plotted bands correspond to the 1σ intervals obtained from 100 random iterations of the final sampling chains. We compared our four new extended image models with the point-like model from S24 and found significantly smaller statistical uncertainties for the extended image models. The positions of the multiple images are marked with small dashes on the lower x-axis and color-coded according to their redshifts. The bar lengths reflect the positional uncertainties (a small dash corresponds to an HST identification, while a long dash corresponds to a MUSE-only identification; see S24 for details). We further show the positions of the 77 000 HST image pixels used in the extended image models with a gray histogram, where the height (in log-scale) corresponds to the number of pixels per radial bin.

In the text

Fig. 6. Median values and 1σ uncertainties of time delays obtained at the model-predicted positions of the multiple images of the SN with fixed (flat ΛCDM) cosmology. The four new extended image models and the two previous point-like models by S24 and G16 are shown with different colors. We notice the significantly lower statistical uncertainties for the extended image models compared to the point-like models. We note again that the extended image models are constructed with an approximate point spread function (see Sect. 2.3 for details), which slightly affects the exact values of the predicted time delays.

In the text

Fig. 7. Color composite image of the SN host reconstruction from the Hubble Frontier Field HST bands F814W (red), F606W (green), and F435W (blue), using the Ew model. The mapped region was slightly modified compared to that used in the model optimization (shown in Fig. 2) to minimize light contamination from some cluster members and the SN images. This figure underscores our ability to directly model the SN host in multiple bands. The image has an extension of 2$\stackrel{″}{.}$05 and 2$\stackrel{″}{.}$27 along the x and y axes, which correspond to 17.84 kpc and 19.75 kpc, respectively, at a redshift of zSN = 1.49.

In the text