© The Authors 2025

1. Introduction

Strong gravitational lensing by galaxy clusters can be used as a powerful and versatile tool in astrophysics and cosmology. For instance, the sensitivity of gravitational lensing to the total mass distribution of an object allows us to probe its dark matter (DM) content (e.g. Grillo et al. 2015; Schuldt et al. 2019; Limousin et al. 2022; Wang et al. 2022; Cha & Jee 2023). Properties of high redshift galaxies that would otherwise be undetected are revealed thanks to the lensing magnification effect that turns gravitational lenses into ‘cosmic telescopes’ (e.g. Coe et al. 2012; Alavi et al. 2016; Acebron et al. 2018; Bouwens et al. 2022; Castellano et al. 2023).

In recent years, supernovae (SNe) strongly lensed by galaxy clusters have gained importance as a probe of cosmological parameters such as the Hubble constant, H₀, which is the current expansion rate of the Universe. Its value has been measured from the cosmic microwave background with the Planck satellite (Planck Collaboration VI 2020, H₀ = 67.4 ± 0.5 km s⁻¹ Mpc⁻¹), and through the type Ia SN distance ladder calibrated with Cepheid (Riess et al. 2022, H₀ = 73.0 ± 1.0 km s⁻¹ Mpc⁻¹) or in combination with stars from the tip of the red giant branch and J-region asymptotic giant branch stars (Freedman et al. 2025, H₀ = 69.96 ± 1.05 (stat) ± 1.12 (sys) km s⁻¹ Mpc⁻¹). The increased precision of these measurements from Planck Collaboration VI (2020) and Riess et al. (2022) has led to a ∼5σ tension. This tension could reveal unknown systematic effects in the measurements (e.g. Efstathiou 2020; Yeung & Chu 2022; Freedman & Madore 2023; Riess et al. 2024) or could be evidence for new physics beyond the standard flat ΛCDM (Valentino et al. 2021) model, which is a spatially flat universe consisting of dark energy described by the cosmological constant Λ and cold dark matter (CDM).

An independent way of measuring H₀ is through time-delay cosmography. This method, first proposed by Refsdal (1964), makes use of the variability of objects such as quasars or SNe, which are observed multiple times through the strong lensing effect. Because of the different path lengths the light takes to get to the multiple image positions and the differences in the gravitational potential, the variability will be observed with a time delay between the images. The time delays are measured from light curves obtained from monitoring by the COSmological MOnitoring of GRAvItational Lenses (COSMOGRAIL; Courbin et al. 2004; Vuissoz et al. 2008; Courbin et al. 2011; Tewes et al. 2013; Millon et al. 2020a,b) collaboration for lensed quasars, for example. The time delays are related to the lensing potential, ψ, of the lens, and the so-called time-delay distance, D_Δt, which is proportional to 1/H₀ (Suyu et al. 2010).

The method of time-delay cosmography has been applied successfully to lensed quasars (e.g. Suyu et al. 2010, 2014; Wong et al. 2017, 2024; Jee et al. 2019; Birrer et al. 2019; Chen et al. 2019; Rusu et al. 2020; Shajib et al. 2020) in the H₀ lenses in COSMOGRAIL’s Wellspring (H0LiCOW; Suyu et al. 2017), COSMOGRAIL, STRong lensing Insights into the Dark Energy Survey (STRIDES; Treu et al. 2018), and the Strong lensing at High Angular Resolution Program (SHARP; Lagattuta et al. 2012; Chen et al. 2022) collaborations. The H₀ inferences of six individual lenses were combined by the H0LiCOW collaboration for a joint constraint on H₀, yielding $H_0=73.3_{-1.8}^{+1.7}{\mathrm{km}\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ (Wong et al. 2020). The Time-Delay COSMOgraphy (TDCOSMO) collaboration showed that relaxing the assumptions on the mass model and using only stellar velocity dispersion measurements of the lens galaxies to constrain the shape of the mass profile results in $H_0=74.5_{-6.1}^{+5.6}{\mathrm{km}\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, i.e. the uncertainty on H₀ increases from 2% to 8% but the median value does not change significantly (Birrer et al. 2020).

Lensed SNe can be used in a similar way to determine the Hubble constant. The advantages of SNe over quasars for H₀ inference are that their light curves change drastically on timescales of weeks, thereby providing precise time delays with less monitoring (Pierel & Rodney 2019; Pierel et al. 2021; Huber et al. 2022; Huber & Suyu 2024). Furthermore, SNe, in contrast to quasars, will fade away, making the modelling of the light distribution of the lens galaxy or the arc more precise (Ding et al. 2021). In the case of (at least) SNe Ia, chromatic microlensing effects (due to stars and compact objects in the lensing galaxy) on the light curves are mitigated in the early parts of colour curves when microlensing is achromatic due to similar spatial intensity distributions across wavelengths in the early phases of SNe (e.g. Goldstein et al. 2018; Huber et al. 2021). For sources lensed by galaxy clusters, the time delays are of the order of months to years and can be measured with precisions of 1–3% (Fohlmeister et al. 2013; Dahle et al. 2015; Kelly et al. 2016a, 2023a). From the eight lensed SN systems discovered so far (six cluster-scale lenses and two galaxy-scale lenses), two had sufficiently long time-delay measurements for time-delay cosmography, both lensed by galaxy clusters.

The first H₀ measurement from a lensed SN came from images of SN Refsdal (Kelly et al. 2015, 2016b) at z = 1.49 in the Hubble Frontier Field cluster MACS J1149.5+2223. It was classified as a type II SN based on the shape of the light curve and from spectroscopy (Kelly et al. 2016b). With the measured time delays and combination of eight lens models, Kelly et al. (2023b) inferred $H_0=64.8_{-4.3}^{+4.4}{\mathrm{km}\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, and $H_0=66.6_{-3.3}^{+4.1}{\mathrm{km}\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ using the two models most consistent with the observations, which is roughly a 6% precision on H₀. By strictly following their blind modelling methodology presented and frozen in Grillo et al. (2016), Grillo et al. (2024) measured $H_0=65.1_{-3.4}^{+3.5}{\mathrm{km}\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ with SN Refsdal using more relaxed assumptions on the background cosmology. The H₀ measurement from Grillo et al. (2024) is robust with respect to cosmological model assumptions as a result of the presence of multiple lensed background sources at different redshifts being lensed by the cluster, in addition to SN Refsdal that provides time-delay constraints.

After SN Refsdal, a SN called ‘SN H0pe’ was detected in the JWST imaging of the PLCK G165.7+67.0 (G165) galaxy cluster field (z = 0.35), which was triply imaged as a result of strong lensing by G165 (PID 1176; Windhorst et al. 2023; Frye et al. 2023, 2024). Follow-up JWST NIRCam and NIRSpec observations (PID 4446) confirmed the SN to be a type Ia at z = 1.783 (Frye et al. 2024; Chen et al. 2024; Polletta et al. 2023). Using the three images of SN H0pe, the two relative time delays and three absolute magnifications were measured by Pierel et al. (2024a) using photometry and Chen et al. (2024) using spectroscopy. The values for these five observables were also predicted by seven different lens modelling approaches all incorporating identical lensing evidence and blinded from each other and from the time-delay measurement team. A joint fitting of the five predicted observables to the measured values yielded $H_0=75.4_{-5.5}^{+8.1}$ km s⁻¹ Mpc⁻¹ (Pascale et al. 2025). This is only the second measurement of H₀ from a lensed SN after SN Refsdal, and the first using a standardizable candle (a type Ia SN). We refer to the above references, especially Pascale et al. (2025), for more details on the SN H0pe analysis. The measurements from SN Refsdal and SN H0pe are consistent with each other at the ∼1.5 − 2σ level.

With the recent discovery of the lensed SN Encore (Pierel et al. 2024b), a third system suitable for time-delay cosmography is added to the sample. SN Encore is a type Ia SN at redshift z = 1.95 that is lensed into multiple images by the galaxy cluster MACS J0138−2155 (Newman et al. 2018a,b). The same source galaxy was previously host to another likely type Ia SN, called SN Requiem, which was discovered in archival Hubble Space Telescope (HST) data taken in 2016 (Rodney et al. 2021). Since SN Requiem was discovered after the SN had faded, it could not be used to obtain a precise H₀ measurement, but Rodney et al. (2021) predicted the appearance of a future image in the year ∼2037, the detection of which will yield a time delay with percent-level uncertainty. In contrast, SN Encore was discovered while the SN was still visible, and Pierel et al. (2024b) and Granata et al. (2025) were able to acquire additional observations of SN Encore and the galaxy cluster. Ongoing time-delay measurements of SN Encore (Pierel et al., in prep.) and the giant lensing arcs of the SN host galaxy provide an excellent opportunity to use this cluster-scale lens as a probe for H₀. In order to infer the Hubble constant, an accurate model of the lensing galaxy cluster MACS J0138−2155 is needed in addition to the time delays. Furthermore, the mass model could be used to predict the properties of future images of SN Encore and SN Requiem, which would facilitate their detections.

In this paper, we present the data processing approach and have collected all the necessary ingredients for building a cluster mass model. In particular, we measured the photometry of galaxies in the field of the cluster and combined them with newly acquired spectroscopic observations presented in Granata et al. (2025) in order to identify galaxies that are members of the galaxy cluster. This dataset, together with the spectroscopically identified lensed background sources, was shared with seven independent modelling teams using a variety of modelling software. The cluster mass modelling by each team was performed independently in a blind analysis where the findings of each team were not communicated until all teams finalized their mass models. The comparison of the mass modelling results from the blind analysis is presented in an upcoming paper (Suyu et al., in prep.). In this paper, we present the modelling details of the GLEE¹ team, which is composed of the first six authors. The GLEE model encapsulates the dark matter and galaxy mass components parametrically, and employs multi-lens-plane ray tracing. In order to maintain a blind analysis throughout our work, the first part of this paper draft (Sects. 1–4), which contains the input data, was shared throughout this work with all modelling teams as the teams checked and all agreed on the same input data for building their mass models. The rest of the paper draft (Sects. 5–7) that specifically describes the GLEE mass model was only accessible to the GLEE team before unblinding, and was subsequently shared with other teams after unblinding. The model results from our blind analysis are presented here without further modification. The H₀ constraint from our mass modelling will be presented in Pierel et al. (in prep.).

The outline of the paper is as follows. We present the observational data used in our analysis in Sect. 2 and the photometric measurements of galaxies in Sect. 3. Through the spectroscopic data presented in Granata et al. (2025) and the photometric measurements, we identified galaxies that are members of the galaxy cluster, additional galaxies along the line-of-sight that are potentially significant for the cluster mass model, and background source galaxies that are strongly lensed by the cluster in Sect. 4. The methods and details of our mass modelling are presented in Sect. 5 and the results, including predictions of SNe Encore and Requiem’s positions and magnifications, are in Sect. 6. We summarize in Sect. 7.

Throughout the paper, magnitudes are reported in the AB magnitude system and parameter estimates are given by the median of its one-dimensional marginalized posterior probability density function. The quoted uncertainties show the 16th and 84th percentiles (corresponding to a 68% credible interval). Since the time-delay measurements are not yet available, our mass models in this work adopt a flat ΛCDM cosmology with $H_0=70\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ and Ω_M = 1 − Ω_Λ = 0.3, especially to predict the properties of SN Encore and SN Requiem.

2. Observations

MACS J0138−2155 (shown in Fig. 1) is a galaxy cluster at redshift z_d = 0.336, discovered within the MAssive Cluster Survey (MACS; Ebeling et al. 2001; Repp & Ebeling 2018). This cluster acts as a strong gravitational lens and magnifies the light of a near-infrared (NIR) source galaxy MRG−M0138 at a redshift z_s = 1.949 (Newman et al. 2018b) into two tangential arcs and one radial arc, which makes this system one of the brightest extragalactic objects in the NIR (Newman et al. 2018a,b). The bright central region of the source galaxy is lensed into at least five multiple images located in the giant tangential arc (∼20″ wide) to the south of the cluster centre, the tangential arc (∼15″ wide) west of the brightest cluster galaxy (BCG), and the radial arc north-east of the BCG. The source galaxy is the host of SN Encore (Pierel et al. 2024b), which is itself lensed into multiple images. Three images 1a, 1b, and 1c of SN Encore have been observed in the two tangential arcs (see white circles in Fig. 1, although image 1c has low signal-to-noise ratio and is not clearly visible²). A fourth image 1d is expected in the future in the radial arc, and a possible central image 1e (fifth image, also in the future with the longest time delay) would be in the central region of the bright BCG (dashed white circles in Fig. 1). To build the total cluster mass model we use imaging data from the Hubble Space Telescope (HST) and the James Webb Space Telescope (JWST), as well as spectroscopic observations from Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010), which are described in this section.

Fig. 1. Colour image of MACS J0138−2155 from combining JWST/NIRCam and HST/WFC3 data. The RGB image is constructed from the filters F105W+F115W+F125W (blue), F150W+F160W+F200W (green), and F277W+F356W+F444W (red). The observed images of SN Encore are shown in solid white circles, the expected positions of images D and possibly E are in dashed white circles. The observed and predicted images for SN Requiem are shown in blue following the same pattern. The images were drizzled to a pixel scale of 0.04″/pix. (Image Credit: STScI, A. Koekemoer, T. Li).

2.1. HST and JWST imaging

The first visit of HST was in 2016 (Proposal ID 14496, PI: Newman) when ACS and WFC3 images were taken. In the scope of the REQUIEM galaxy survey (Akhshik et al. 2020, 2023), 6 orbits of WFC3 observations were performed in 2019 (Proposal ID 15663, PI: Akhshik). SN Encore was discovered in JWST data taken in November 2023 (Proposal ID 2345, PI: Newman), which led to follow-up observations in December 2023 and January 2024 with HST (Proposal ID 16264, PI: Pierel) and JWST (Proposal ID 6549, PI: Pierel). We summarize the HST and JWST observations in Table 1 and Table 2, respectively. The newly obtained HST and JWST imaging data and reduction are described in detail in Pierel et al. (2024b).

All available imaging data were combined into mosaics with an updated version of the HST ‘mosaicdrizzle’ pipeline (Koekemoer et al. 2011). This produced separate epoch mosaics, as well as full depth mosaics where all epochs per filter were combined. Each mosaic product contains multiple extensions, including the science image and the error image (providing the 1σ uncertainties of the science image pixels). The absolute astrometry of all mosaics is directly aligned to Gaia-DR3 (Prusti et al. 2016; Vallenari et al. 2023). In Table 3, we summarize all imaging data products that are used in our analysis. The filters were chosen such that a large wavelength range is covered with preference for filters with a large total exposure time. The final data products are drizzled to a pixel size of 0.04″/pixel for all the filters in order to balance between the different native pixel sizes of the different detectors and to sample sufficiently the point-spread functions (PSFs).

2.2. MUSE spectroscopy and redshift measurements

MACS J0138−2155 was the target of two programs carried out by the integral-field spectrograph MUSE (Bacon et al. 2010), at the Very Large Telescope of the European Southern Observatory. The first ∼49 minutes were obtained within the program ID 0103.A-0777 (September 2019, PI: Edge), later complemented by additional 2.9 hours of data collected after the detection of SN Encore in the Target of Opportunity program ID 110.23PS (December 2023, PI: Suyu). The average seeing across the observations was approximately 0.8″. The data reduction pipeline adopted for the creation of the final data-cube is described in detail in Granata et al. (2025).

The MUSE spectroscopy was used in Granata et al. (2025) as the basis for the construction of a redshift catalogue containing 107 reliable redshift measurements of extragalactic objects. They identified 50 cluster members in the redshift range z = 0.324 − 0.349, that yielded a cluster redshift of z_d = 0.336. Finally, they performed a blind search of spectral features appearing simultaneously at different locations across the data cube, which allowed for the identification of two additional background sources systems 5 and 6 with secure redshifts that are strongly lensed by the galaxy cluster, and an additional source with a tentative redshift measurement that is strongly lensed by a cluster member (this source is not included in our final model). In total, there are four independent background source galaxies with secure spectroscopic redshifts from z = 0.767 to 3.420 that are strongly lensed into multiple image systems.

3. Photometry and Sérsic profiles for objects in MACS J0138−2155

Accurate photometry of galaxies in the cluster field of view is important to photometrically select cluster members for the mass model, as well as to scale their mass parameters with a scaling relation for the modelling. Therefore, we fitted the surface brightness of several objects in the field using Sérsic profiles (Sérsic 1963) described in Sect. 3.1. The Sérsic profile fitting gives accurate photometry, especially in crowded fields of galaxy clusters. We subtracted the extended light of the BCG to avoid biases in the photometry of nearby galaxies. This procedure is described in Sect. 3.2, and the photometry of other objects in the field in Sect. 3.3.

3.1. Sérsic profiles

To measure the magnitudes of the galaxies, we fitted their two-dimensional surface brightness (SB) with Sérsic profiles (Sérsic 1963) in all filters listed in Table 3. The Sérsic profile is parametrized as

$$\begin{array}{c}\hfill I_{\mathrm{S}}(x,y)=A_{\mathrm{S}}exp[-k\{{\left(\frac{\sqrt{{(x-x_{\mathrm{S}})}^2+{\left(\frac{y-y_{\mathrm{S}}}{q_{\text{S}}}\right)}^2}}{r_{\mathrm{eff}}}\right)}^{\frac{1}{n_{\mathrm{s}}}}-1\}],\end{array}$$(1)

where A_S is the amplitude, x_S and y_S are the centroid coordinates, q_S is the axis ratio (0 < q_S ≤ 1), and n_S the Sérsic index. The value for the normalization constant k is set such that r_eff is the half-light radius in the direction of the semi-major axis (k is thus not a freely variable parameter). The light distribution is rotated by a position angle ϕ_S, which is measured east of north (where ϕ_S = 0° corresponds to the major axis being aligned with the northern direction).

We constructed the point-spread function (PSF) from multiple stars in the field with the software STARRED (Michalewicz et al. 2023; Millon et al. 2024). We manually chose bright, unsaturated stars in the field. The number of stars used depends on the field of view (FOV), but we used at least three stars for each filter.

The galaxy SB can be modelled as a combination of multiple Sérsic components to obtain SB distribution that fit better to the observations. In that case, the total SB is

$$\begin{array}{c}\hfill I_{\mathrm{S},\mathrm{tot}}(x,y)={\displaystyle \sum _{i=1}^{N_{\mathrm{sersic}}}}{I}_{\mathrm{S},i}(x,y),\end{array}$$(2)

where N_sersic is the number of Sérsic components.

The Sérsic parameters are optimized to fit to the observed intensity value I_i^obs of pixel i by minimizing the χ² function of the SB

$$\begin{array}{c}\hfill {\chi }_{\mathrm{SB}}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{p}}}}\frac{|I_i^{\mathrm{obs}}-{(\mathrm{PSF}\otimes {\mathrm{I}}_{\mathrm{S},\mathrm{tot}})}_i{|}^2}{{\sigma }_{\mathrm{tot},i}^2},\end{array}$$(3)

where N_p is the number of data pixels in the cutout, ${\sigma }_{\mathrm{tot},\mathit{i}}$ is the uncertainty of pixel i. The symbol ⊗ represents the convolution of the PSF with the predicted Sérsic intensity I_S, tot.

3.2. BCG subtraction

Brightest cluster galaxies are among the most massive and bright galaxies in the Universe. It is thus important to subtract off its light to avoid bias in the photometry of nearby cluster members and contamination of the giant arcs. In addition, clusters often exhibit intracluster light (ICL; Contini 2021; Montes 2022), an extended and diffuse light envelope caused by unbound stars in the gravitational well of the DM halo. While BCGs can often be well fit with two Sérsic profiles to account for the inner and outer parts of the galaxy separately (Gonzalez et al. 2005; Donzelli et al. 2011; Chu et al. 2021), the disentanglement of BCG and ICL light is not straightforward (Kluge et al. 2020, 2021). The ICL often has to be modelled as a separate component (e.g. Joo & Jee 2023). We fitted the BCG light with three Sérsic profiles to also account for possible ICL contribution using the modelling software GLEE (Suyu & Halkola 2010; Suyu et al. 2012). The light centroids are linked among the multiple Sérsic profiles. We found that we fitted the extended light distribution at the centre of the cluster slightly better (in terms of χ_SB² and by a visual check of the modelling residuals) with three Sérsic profiles instead of two. In Fig. 2, we show, as an example, the observed image, the BCG model, and the BCG subtracted image in the HST F160W filter, which is the filter we later used to obtain the cluster member scaling relations given the well-understood photometric calibrations of HST (Sect. 5.1). At this wavelength, the BCG has a magnitude of m_BCG, F160W = 15.30 with a best-fit axis ratio of q_BCG, F160W = 0.58 and position angle ϕ_BCG, F160W = 53.2° (east of north), measured from its second brightness moments. The central region of the BCG exhibits prominent dust features, which leads to strong residuals in the central ∼1″, but overall the extended light is subtracted well. The subtraction also reveals more details of the radial arc and even a central image.

Fig. 2. Results from modelling the BCG light with Sérsic profiles in the HST F160W filter. The left panel shows the observed image, the middle panel shows our BCG light model, and the right panel is the BCG subtracted image. The images cover the central 50″ × 50″ around the BCG.

3.3. Sérsic fitting of galaxies

We modelled the surface brightness of objects in the HST/JWST fields of the cluster with the software Morphofit (Tortorelli et al. 2018, 2023; Tortorelli & Mercurio 2023), which fits to galaxies in an automated and parallelized way. In Morphofit, the initial object detection and estimation of Sérsic parameter values is done with SExtractor (Bertin et al. 1996) and the final Sérsic fit is conducted with the modelling software Galfit (Peng et al. 2011). The Sérsic parametrization is the same as in Eq. (1). Morphofit will create cutouts of all objects and model them in parallel on multiple cores. We selected objects from the SExtractor catalogue with F160W magnitudes MAG_AUTO < 23.5, i.e., all objects with an automatically detected F160W magnitude fainter than MAG_AUTO = 23.5 are rejected to limit the number of objects in the mass model by including only the most massive ones. In a lens cluster at a similar redshift, Bergamini et al. (2023) identified a threshold magnitude of m_F160W ≤ 21; the inclusion of a fainter object (fainter than this threshold) in the mass model generally has no significant impact. Nevertheless we chose a uniform, inclusive approach and included all cluster members within the FOV that are above our magnitude threshold.

Initially, we modelled each galaxy with two Sérsic profiles with forced photometry (i.e., linking the centroid coordinates) among the multiple filters. Models with a single Sérsic component showed insufficient light subtraction, especially in the brighter central parts. Adding a second Sérsic component accounts for the bulge and disc decomposition that we observe in many of the galaxies in the field. For simplicity we decided to uniformly model all galaxies with two Sérsic profiles. After a visual quality check, we added a third Sérsic profile to the models that still show an excess of light in the residuals that could potentially impact the photometry. Other bright objects, which are sufficiently close, are fit simultaneously to avoid bias in the photometry. We fitted a total of 482 galaxies (in addition to the BCG) across all bands and estimated their total integrated AB magnitudes. All fitted objects are labelled uniquely with a photometric ID. The magnitudes for the cluster members selected in Sect. 4.1 and the background galaxy (see Sect. 4.2) are summarized in Table A.1, which can be downloaded online³.

In Fig. 3, we show as examples the modelled surface brightness and residuals for several objects close to the BCG. The three blended objects 92, 103, and 106 in the first row of the plot highlight the necessity of fitting Sérsic profiles to the galaxies instead of using automatic aperture photometry from SExtractor alone. We could deblend these objects and got reliable magnitudes for the individual galaxies, although there are some extended residual features, possibly due to tidal interactions of the galaxies. The second row shows the massive background galaxy at a redshift of z_bg = 0.371, and the third row a bright cluster member just south of the BCG.

Fig. 3. Examples for Sérsic fits to objects in the field in the JWST F200W filter. From left to right: data image, model image, and residual image. The photometric IDs of the objects are shown in blue on the data image. Object 82 is a bright background galaxy at a redshift zbg = 0.371, and the other galaxies are cluster members close to the BCG.

4. Deflector galaxies and multiply lensed sources

In our lens mass model of MACS J0138−2155, we included as deflectors spectroscopically confirmed and photometrically selected cluster members. In total, we found 84 objects at the cluster redshift of z_d = 0.336, including the BCG. We describe our cluster member selection in Sect. 4.1. In addition to the mass at the cluster redshift, we included two line-of-sight (LOS) perturber galaxies, one in the foreground (z_fg = 0.309) and one in the background (z_bg = 0.371). We describe the LOS objects in Sect. 4.2. Our image-position constraints, summarized in Sect. 4.3, consist of 23 multiple images from eight image systems at four different source redshifts, selected from JWST imaging and MUSE spectroscopic data.

4.1. Cluster members

Since MUSE covers only the central 1 × 1 arcmin² region of MACS J0138−2155 and some cluster members are located outside this area, we used the photometric information of the galaxies described in Sect. 3 to identify additional members. In particular, we used our spectroscopic catalogue to identify spectroscopic cluster members, and used both colour-magnitude and colour-colour diagrams to further select cluster members photometrically to complement the spectroscopic members.

We used the sample of 50 spectroscopically confirmed cluster galaxies in the redshift range z = 0.324 − 0.349 identified by Granata et al. (2025) as the basis for our catalogue of members. These galaxies are marked by cyan squares on the JWST image in Fig. 4. After excluding the three jellyfish galaxies⁴, we used the remaining spectroscopic objects to define the red sequence in the colour-magnitude diagram shown in Fig. 5. We used the F555W−F814W colour from HST, given that these two filters roughly bracket the rest-frame 4000 $\mathrm{\AA }$ break of the cluster members, providing a tighter red sequence that is helpful for selecting additional members photometrically. The spectroscopically confirmed cluster members are shown as blue points, whereas spectroscopically confirmed non-members are red. Galaxies with photometric measurements from Sect. 3 and without spectroscopic redshifts are in grey.

Fig. 4. Colour image of MACS J0138−2155 built from JWST NIRCam filters (red: F277W, F356W, F444W; green: F150W, F200W; blue: F115W). Spectroscopically confirmed objects are marked with squares, coloured in cyan for cluster members, and blue and yellow for the foreground and background objects, respectively. Objects selected from the colour-magnitude fit are marked as red circles, and green triangles are the colour-colour selections.

We fitted a linear relation to the spectroscopic cluster members (blue points) through a least-squares fit and obtained the 1σ uncertainty on this linear relation (the blue line and the shaded blue region). We performed a σ-clipping (at 2σ) to exclude outliers in spectroscopic members, which are marked by red crosses⁵. The BCG, which is more than three magnitudes brighter than the other cluster members in this cluster, was clipped in this process. We then repeated the fitting and σ-clipping until no additional spectroscopic cluster members are detected/excluded as outliers, i.e., all remaining spectroscopic cluster members lie within 2σ of the linear fit (within twice the thickness of the light-shaded blue 1σ region).

All the galaxies without spectroscopy, but with colours and magnitudes lying within 1σ of the red sequence (in shaded blue), were selected as photometric cluster members and are marked by yellow plus signs in Fig. 5. In total, 28 galaxies were selected as cluster members with the colour-magnitude method. We however excluded four galaxies that are marked by red minus signs in the red sequence since we cannot readily use their photometry in the scaling relations of elliptical galaxies – one of them is a spiral galaxy, one is outside the F160W FOV without F160W photometry, and two of them have unreliable photometry from the elliptical Sérsic fit (one is located at the edge of the WFC3 chip and one has non-elliptical morphology). The resulting 24 colour-magnitude selected cluster members are marked by red circles on the JWST image shown in Fig. 4.

Fig. 5. Colour-magnitude diagram in the field of MACS J0138−2155. Spectroscopically confirmed cluster members and non-members are indicated as blue and red points, respectively. Sources with no spectroscopic confirmation are shown in grey. The blue line and region indicate the fitted red sequence and the 1σ scatter. Red crosses (X) mark cluster members that do not follow the red sequence and yellow plus signs indicate the photometrically selected members that are included in our lens model. Objects inside the red sequence with a red minus sign are not included as photometric cluster members due to the lack of, or unreliable, F160W photometry.

We supplemented the colour-magnitude selection of cluster members by using additional colour-colour selection. Following Frye et al. (2024), we used the F277W−F444W and F814W−F150W colours, where we substituted the JWST F090W filter in Frye et al. (2024) with the HST F814W filter due to their similar wavelengths and the lack of JWST F090W observations for SN Encore. Figure 6 shows the colour-colour diagram on two different scales, where the right-hand panel is a zoom-in of the left-hand panel. The symbols are the same as in Fig. 5, where red points are spectroscopic non-members, blue points are spectroscopic members, grey circles are objects without spectroscopic redshifts, red crosses are spectroscopic members that are not within the red-sequence region, yellow plus signs are photometric members (grey) that lie within the red sequence and were selected as cluster member, and red minus signs are objects not included as cluster members due to unreliable/missing F160W photometry. The sizes/transparency of the filled grey circles are related to the brightness of the galaxies, where brighter galaxies in F160W are represented by bigger and darker dots (as shown by the colour bar). Open grey circles are galaxies without F160W photometry. In the right-hand panel, we define a region with −0.80 < F277W − F444W < −0.68 and 0.75 < F814W − F150W < 1.3 to include photometric galaxies (filled grey circles) as cluster members, some of which are on the brighter side.

Fig. 6. Colour-colour diagram in the field of MACS J0138−2155, where the right panel is a zoom-in of the left panel. Spectroscopically confirmed cluster members and non-members are indicated in blue and red, respectively. Sources with no spectroscopic confirmation are shown in grey, with size and shade related to the galaxy brightness as shown in the colour bar. Red crosses (X) mark spectroscopic cluster members that do not follow the red sequence and yellow plus signs indicate the colour-magnitude selected members (from Fig. 5) that are included in our lens model. Photometric galaxies inside the blue shaded region in the right-hand panel are selected as cluster members via our colour-colour selection.

We tuned this colour-colour selection region (light blue shaded region) to include the brighter photometric members while excluding the spectroscopic non-members. In total, we selected 32 photometric galaxies via the colour-colour diagram, of which 14 are additional cluster members that were not previously selected through the colour-magnitude diagram. However, two of the 14 galaxies do not have HST F160W photometry as they are outside the F160W field of view, and two other galaxies are located very close to the edge of the FOV. We excluded these four galaxies (marked by red minus signs in Fig. 6) from the model since we cannot assign mass to them via scaling relations without the F160W photometry. Their effect on the mass model is likely to be minimal given their faintness compared to other cluster members and their large distance to the cluster core. The 28 colour-colour selected galaxies are marked by triangles in Fig. 4.

As visible in Fig. 4, our spectroscopic, colour-magnitude, and colour-colour selection of cluster members include significantly massive early-type galaxies that are crucial for incorporating into the lens mass model, especially those galaxies near the cluster centre where we have lensed multiple image systems. The union of the colour-magnitude and colour-colour selections ensures a higher completeness in the cluster members, especially for including a few of the massive cluster member that are missed in one of the selections. The objects that are close to but outside of the borders of the selection regions (either in colour-magnitude in Fig. 5 or colour-colour in Fig. 6) are either faint or substantially distant from the cluster centre, and their lensing effects are expected to be negligible. In total, we included 50 spectroscopically selected, and 34 photometrically selected cluster members in our model.

4.2. LOS objects

We checked for potentially significant line-of-sight objects for cluster mass modelling, i.e. galaxies near arcs or that appear to be massive. We spectroscopically confirmed a small foreground galaxy close to the giant arc at a redshift of z_fg = 0.309 (RA, Dec = 24.5159632°, −21.9300802°), and a massive background galaxy at z_bg = 0.371 (RA, Dec = 24.5136457°, −21.9244117°), about 8″ northwest of the BCG. We included both into our mass model given their proximity to the giant arcs or the high mass of the background galaxy (based on its stellar velocity dispersion measurement). The final selection of galaxies to be included in the model is shown in Fig. 4. This includes cluster members, as well as LOS objects.

4.3. Strongly lensed background sources as constraints

We securely identified eight multiple image systems from four independent background galaxies, with a total of 23 images in the redshift range of 0.767 ≤ z ≤ 3.420, which are summarized in Table 4. These multiple image systems were voted on and agreed upon by eight groups of researchers⁶, as described in Suyu et al. (in prep.). Each multiple image system was labelled numerically, and the individual lensed images of each system were labelled alphabetically, in decreasing order of brightness. The multiple images of SN Encore (JWST) and SN Requiem (HST) are Systems 1 and 2, respectively. System 3 is associated with the corresponding emissions from the centre of the SN host galaxy at z = 1.949. The positions were identified from the JWST images. Systems 4.1, 4.2, and 4.3 correspond to [O II] emissions at z_s, 4 = 0.767 identified in the MUSE data, where different clumps of emission are lensed into multiple image systems; since these emissions likely come from the same galaxy, we labelled them with the same starting number 4. The image positions of Systems 4.1 and 4.3 are based on MUSE data, whereas System 4.2 is based on JWST images. Additionally, two Ly-α emission regions were found in MUSE at a redshift of z_s, 5 = 3.152 and z_s, 6 = 3.420, which we labelled as Systems 5 and 6, respectively; these two Ly-α sources do not have clear detections in the JWST imaging. We assigned elliptical errors on their positions based on the spatial extent of the images in the data (in either JWST or MUSE data, depending on the data in which the image positions weredetected/measured).

In addition to these eight multiple image systems that were designated as “gold” systems, there is one more lensed image system that was designated as “silver” (where we refer to Suyu et al. (in prep.) for the criteria of gold and silver systems). The silver system corresponds to the spectroscopic ID 1755 (Granata et al. 2025), with a “likely” instead of a “secure” redshift measurement coming from two blended background sources. Furthermore, the source is primarily lensed by a single cluster member rather than the full cluster, so the lensing configuration is mainly informative on the total mass distribution of the specific cluster member rather than the global cluster mass distribution. This specific cluster member is ≳17″ away from all the SN images and therefore has negligible impact on the SN image properties (of both SN Encore and SN Requiem). Therefore, we did not include this silver system into our mass modelling. Figure 7 shows an overview of the eight multiple image systems and their redshifts.

Fig. 7. Positions of the eight multiple image systems (the ‘gold’ systems) used to constrain our mass model. We show the positions in arcseconds relative to the BCG at (RA, Dec) = (24.5157084°, −21.9254717°).

5. Lens mass models

Having selected the cluster members, identified important LOS objects, and determined the multiple image systems, we constructed our cluster lens mass model. We describe the parametrization of the different mass components in Sect. 5.1, and present seven different GLEE mass models in Sect. 5.2 in order to quantify systematic uncertainties due to mass modelling assumptions. We then describe the methods of modelling MACS J0138−2155 with the multiple image systems as constraints in Sect. 5.3. The general modelling methodology is very similar to that developed by Grillo et al. (2015) and Grillo et al. (2016) for two lens galaxy clusters of the Hubble Frontier Fields (Lotz et al. 2017).

5.1. Mass parametrizations

5.1.1. Galaxy mass distributions

We modelled the cluster members and LOS galaxies with a truncated dual pseudo-isothermal elliptical mass distribution (dPIE; Elíasdóttir et al. 2007; Suyu & Halkola 2010) with vanishing core radius. Their dimensionless projected surface mass density, i.e., convergence, is

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{dPIE}}(x,y){|}_{{\mathrm{z}}_{\mathrm{s}}=\infty }=\frac{{\theta }_{\mathrm{E},\infty }}{2}(\frac{1}{R_{\mathrm{em}}}-\frac{1}{\sqrt{R_{\mathrm{em}}^2+r_{\mathrm{t}}^2}}),\end{array}$$(4)

where the 2D elliptical mass radius with respect to their centre is

$$\begin{array}{c}\hfill R_{\mathrm{em}}(x,y)=\sqrt{\frac{x^2}{{(1+e)}^2}+\frac{y^2}{{(1-e)}^2}},\end{array}$$(5)

and the ellipticity, e, is related to the axis ratio, q, by

$$\begin{array}{c}\hfill e=\frac{1-q}{1+q}·\end{array}$$(6)

The parameter θ_E, ∞ is the Einstein radius for a source at redshift z_s = ∞, and r_t is the truncation radius. The mass distribution is then rotated by its position angle ϕ. We fixed the coordinates of each object on the image plane to the observed positions from the HST F160W image (from the Sérsic fit described in Sect. 3). For simplicity, we assumed all the cluster members and LOS perturbers to be spherical, with the exception of the BCG for which we considered different scenarios, a spherical BCG and an elliptical BCG, given the substantial mass of the BCG.

To reduce the number of free parameters in our model (such that the number of free parameters does not exceed the number of observables) and increase the computational efficiency of our modelling, we adopted scaling relations for the Einstein and truncation radii of the cluster members⁷, i.e., we scaled them with respect to those of a reference galaxy. The choice of the reference galaxy does not impact the result of the modelling. In our case, we chose to scale the values of θ_E, ∞ and r_t with power laws with respect to the luminosity of the member using object ID_phot = 115 as a reference (m_{F160W,  ID115} = 17.99), so that θ_E, ∞ and r_t for the i-th member are

$$\begin{array}{cc}\hfill {\theta }_{\mathrm{E},\infty ,i}& ={\theta }_{\mathrm{E},\infty }^{\mathrm{ref}}{\left(\frac{L_i}{L_{\mathrm{ID}115}}\right)}^{2\alpha },\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill r_{\mathrm{t},i}& =r_{\mathrm{t}}^{\mathrm{ref}}{\left(\frac{L_i}{L_{\mathrm{ID}115}}\right)}^{\beta },\hfill \end{array}$$(8)

where L_i and L_ID115 are measured in the HST band F160W, and α and β are power-law indices that we determine through properties of cluster members, particularly velocity-dispersion measurements.

Cluster-scale strong lensing is sensitive to the total mass of the member galaxies, but unless galaxy-scale lensing systems are observed (Granata et al. 2023; Galan et al. 2024), the reconstruction of their mass structure is affected by a degeneracy between the values of the Einstein and truncation radii of the members, which affects our predictions on their compactness. Several recent works (Bergamini et al. 2019, 2021; Granata et al. 2022) have shown that this degeneracy can be broken with the introduction of independent kinematic priors on Eq. (7).

The value of θ_E, ∞ is linked to the velocity dispersion parameter of the dPIE profile (Elíasdóttir et al. 2007) by

$$\begin{array}{c}\hfill {\theta }_{\mathrm{E},\infty }=4\pi {\left(\frac{{\sigma }_{\mathrm{dPIE}}}{c}\right)}^2.\end{array}$$(9)

For a vanishing core radius, σ_dPIE is well approximated by the value of the line-of-sight stellar velocity dispersion measured within a small aperture (smaller than the truncation radius, Bergamini et al. 2019). We know that a power-law relation (also known as the Faber-Jackson law, Faber & Jackson 1976) holds between the measured stellar velocity dispersion of elliptical galaxies and their total luminosity. Granata et al. (2025) provided us with a sample of 13 passive cluster galaxies in MACS J0138−2155 with reliable velocity-dispersion measurements within their effective radii. We can thus calibrate the Faber-Jackson law for the cluster members

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{v},i}={\sigma }_{\mathrm{v}}^{\mathrm{ref}}{\left(\frac{L_i}{L_{\mathrm{ID}115}}\right)}^{\alpha },\end{array}$$(10)

where α was defined in Eq. (7). From Eq. (9), we derive ${\theta }_{\mathrm{E},\infty }^{\mathrm{ref}}=4\pi {\left(\frac{{\sigma }_v^{\mathrm{ref}}}{c}\right)}^2$. Granata et al. (2025) fit Eq. (10) adopting a Bayesian technique analogous to Bergamini et al. (2019). They find $\alpha =0.{25}_{-0.05}^{+0.05}$ and ${\sigma }_{\mathrm{v}}^{\mathrm{ref}}={206}_{-13}^{+14}\mathrm{km}{\mathrm{s}}^{-1}$, which is compatible with the stellar velocity dispersion value for the reference member ID 115 (211.1 ± 3.1 km s⁻¹); the scatter is $\mathrm{\Delta }{\sigma }_v={25}_{-4}^{+6}\mathrm{km}{\mathrm{s}}^{-1}$ for the Faber-Jackson relation, which we used to estimate the mean and standard deviation of the Gaussian prior on θ_E, ∞^ref. For the second scaling relation in Eq. (8), we chose β = 0.7, so that $M_{\mathrm{tot},i}/L_i\propto L_i^{0.2}$, consistent with the observed Fundamental Plane relation.

By varying only two parameters, i.e. the Einstein and truncation radii of our reference galaxy, we determined the total mass distribution of all cluster members that follow the scaling relations. In practice, we excluded the BCG and object ID_phot = 116 from the scaling relations and used Gaussian priors on their Einstein radii from their measured velocity dispersions. The BCG is too bright to follow the scaling relation; the galaxy with ID_phot= 116 is close to the critical curve of image systems 4.1, 4.2 and 4.3 and the predicted image multiplicity and positions are very sensitive to the mass distribution of this cluster member. While we assumed that the BCG is spherical in most of our mass models, we also considered an elliptical BCG in one of our mass models. Furthermore, the jellyfish galaxies are free to vary outside the scaling relation as well, with their truncation radii fixed to 5″ and a uniform prior on their Einstein radii (see Tables B.1 and C.1 for details on the assigned priors on the mass components).

For the two LOS galaxies (see Sect. 4.2), we fixed their centroids to the observed light distributions. The foreground galaxy does not have reliable photometry, and given its proximity to the southern giant arc, we allowed both its Einstein and truncation radii to vary. The background galaxy has a measured velocity dispersion, which we used to impose a prior on its Einstein radius.

5.1.2. Cluster DM halo mass distributions

We included two cluster DM (dark matter) haloes in our model, since we found that a single cluster DM halo cannot reproduce the positions of the eight multiple image systems well. We found that the primary DM halo is more centrally concentrated and massive, whereas the secondary DM halo is extended and cored. The multiple image systems require that the centroid of the secondary DM halo to be offset from that of the primary DM halo. As will be seen in Sect. 6, the secondary halo is a perturbation of the primary halo to allow for non-elliptical halo mass distribution. Hereafter, we refer to the secondary halo as a perturbative halo that alters the primary halo structure, as opposed to a separate and physically present DM halo. For the primary DM halo, we considered three possible profiles: (1) pseudoisothermal elliptical mass distribution (PIEMD; Kassiola & Kovner 1993), (2) softened power-law elliptical mass distribution (SPEMD; Barkana 1999), and (3) elliptical Navarro Frenk & White profile (NFW; Navarro et al. 1997; Oguri 2021). For the perturbative DM halo, we modeled it with a PIEMD.

The convergence for a PIEMD is

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{PIEMD}}(x,y){|}_{{\mathrm{z}}_{\mathrm{s}}=\infty }=\frac{{\theta }_{\mathrm{E},\infty }}{2\sqrt{R_{\mathrm{em}}^2+r_{\mathrm{c}}^2}},\end{array}$$(11)

where R_em is the elliptical mass radius (Eq. (5)), e is the ellipticity (Eq. (6)), θ_E, ∞ is the lens Einstein radius for a source at redshift infinity and r_c is the core radius. The mass distribution is then rotated by its position angle ϕ, and centred at (x_DM, y_DM). All six parameters are free to vary within flat priors.

The convergence of a SPEMD is

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{SPEMD}}(x,y)=E{(u^2+s^2)}^{-\gamma },\end{array}$$(12)

where u² = x² + y²/q², q is the axis ratio, s = 2r_c/(1 + q), r_c is the core radius, γ is the radial power-law slope (γ = 0.5 corresponds to an isothermal profile), and E is the strength/amplitude of the lens. The mass distribution is centred at (x_DM, y_DM) and has a position angle ϕ. The SPEMD profile therefore has seven free parameters, and we adopt flat priors for these parameters.

The three-dimensional mass density distribution of the NFW profile is

$$\begin{array}{c}\hfill {\rho }_{\text{NFW}}(r)=\frac{{\rho }_{\mathrm{s}}}{(r/r_{\mathrm{s}}){(1+r/r_{\mathrm{s}})}^2},\end{array}$$(13)

where r is the radial coordinate, ρ_s is the characteristic density, and r_s is the scale radius. This can be projected to obtain the convergence (e.g. Golse & Kneib 2002). The elliptical NFW convergence has six free parameters with flat priors: centroid (x_DM, y_DM), axis ratio (q), position angle (ϕ), normalization (κ₀), and a scale radius (r_s). We followed Oguri (2021) for calculating the lensing deflection angles and potentials of elliptical NFW convergence, and used the implementation from Wang et al. (in prep.).

5.1.3. Mass sheet

Since the mass-sheet degeneracy (MST; Falco et al. 1985; Gorenstein et al. 1988) is one of the dominant sources of uncertainty in time-delay cosmography, we also considered adding a constant sheet of mass parameterized by κ₀ in our mass model. The value of κ₀ is not constrained by a single background source, since any value of κ₀ can produce a new mass distribution that fits the image positions equally well. This new mass distribution consists of adding κ₀ at the cluster redshift while reducing the mass of the other components such as dark matter haloes and galaxies, and produces the same predictions of lensed image positions and relative magnifications of a single source (e.g. Schneider 2019). In our case with four sources at different redshifts, κ₀ can no longer take an arbitrary value and still fit to the observed image positions. The model with the mass sheet therefore allows us to quantify model uncertainties due to the mass-sheet degeneracy (as first shown in Grillo et al. 2020).

5.2. Mass models

We considered seven different GLEE mass models of the galaxy cluster in order to quantify uncertainties due to mass modelling assumptions. We summarize the seven models in Table B.1, including the number of variable mass model parameters and also indicating the priors that are imposed on some of the parameters based on the kinematic measurements of Granata et al. (2025).

The first four models explore the variety of plausible cluster DM halo mass distributions given its dominant role in reproducing the observed image positions. In particular, the primary DM halo is described by a cored-isothermal profile (i.e., PIEMD) in the iso_halo model, a power-law profile (i.e., SPEMD) in the PL_halo model, and a NFW profile in the NFW_halo model. The iso_halo+sheet model is parameterized in the same way as the iso_halo model with the addition of a mass sheet with the amplitude κ₀ as an extra free parameter. The next model, iso_halo+ell_BCG, is similar to the iso_halo model except that the BCG is described by an elliptical isothermal profile instead of a circular one. This assesses the impact of assuming the BCG to be circular. All of these five models have deflectors at different redshifts, given the presence of the foreground and background LOS galaxies, and employ multi-plane lens modelling (e.g. Wong et al. 2017; Chirivì et al. 2018; Wang et al. 2022; Acebron et al. 2024; Schuldt et al. 2024). These models also use the elliptical positional uncertainties listed in Table 4.

The next single_plane model is similar to the iso_halo model but has the foreground galaxy at z_fg = 0.309 (near the southern giant arc) and the bright background galaxy at z_bg = 0.371 placed at the cluster redshift of z_d = 0.336, which reduces the model to a single deflector plane. This helps to assess whether it is important to use multi-plane lensing for this cluster.

The last ‘baseline’ model is considered in order to facilitate direct comparison of modelling results from the various independent modelling teams, which are presented in Suyu et al. (in prep.). Since the teams use different modelling software with varying capabilities, not all software have multi-plane lensing or the usage of elliptical positional uncertainties (instead of circular uncertainties). Therefore, the teams agreed to build the so-called baseline model with all the deflectors assumed to be at the cluster redshift and all the positional uncertainties to be circular, which is a setup that all the teams can perform. For each multiple image position, the size of the circular uncertainty (σ^circ) is the geometric mean of the semi-major (σ^major) and semi-minor (σ^minor) axes of the elliptical uncertainty tabulated in Table 4, i.e., ${\sigma }^{\mathrm{circ}}=\sqrt{{\sigma }^{\mathrm{major}}{\sigma }^{\mathrm{minor}}}$. Such a baseline model will enable a more direct comparison of the modelling software and methodology. Each team can then build their ‘ultimate’ models, relaxing assumptions from the baseline model that they would use for cosmographic analysis.

Having presented the various mass model parameterizations that we have considered, we describe next our procedure to constrain the model parameter values based on the image positions of the multiple image systems.

5.3. Mass modelling with lensed image positions

Our lens mass modelling of MACS J0138−2155 was performed with GLEE (Suyu & Halkola 2010; Suyu et al. 2012). We first optimized the model with the multiple image positions described in Sect. 4.3 by minimizing the angular separation between the observed, θ^obs, and model-predicted, θ^pred, image positions on the image plane. The model parameters η (vector of length N_par, the number of lens parameters) are varied to minimize the χ² function

$$\begin{array}{c}\hfill {\chi }_{\mathrm{im}}^2={\displaystyle \sum _{j=1}^{N_{\mathrm{sys}}}\sum _{i=1}^{N_{\mathrm{im}}^j}}{\left(\mathrm{\Delta }{\mathit{\theta }}_{\mathit{ij}}\right)}^{\mathrm{T}}\mathsf{R}{{\mathsf{C}}_{\mathsf{im}}}^{-1}{\mathsf{R}}^{\mathrm{T}}\left(\mathrm{\Delta }{\mathit{\theta }}_{\mathit{ij}}\right),\end{array}$$(14)

where N_sys is the number of multiple-image systems, N_im^j is the number of multiple images in image system j, and the image positional offset is

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{\theta }}_{\mathit{ij}}={\mathit{\theta }}_{\mathit{ij}}^{\mathrm{obs}}-{\mathit{\theta }}_{\mathit{ij}}^{\mathrm{pred}}(\mathit{\eta },{\mathit{\beta }}_j^{\mathrm{mod}}).\end{array}$$(15)

The two-dimensional vectors ${\mathit{\theta }}_{\mathit{ij}}^{\mathrm{obs}}$ and ${\mathit{\theta }}_{\mathit{ij}}^{\mathrm{pred}}$ are the observed and the predicted positions for image i in image system j, respectively. The predicted image positions ${\mathit{\theta }}_{\mathit{ij}}^{\mathrm{pred}}$ are computed given the lens mass parameters η and the modelled source position ${\mathit{\beta }}_j^{\mathrm{mod}}$ of system j. We obtain ${\mathit{\beta }}_j^{\mathrm{mod}}$ by ray tracing the multiple lensed image positions in system j back to the source plane through the lens equation, and taking the magnification-weighted average of these corresponding source positions. The rotation matrix R, its transpose R^T and covariance matrix C_im of the image positions depict elliptical positional uncertainties. Specifically, for an error ellipse of an image position i of the j^th multiple image system with semi-major and semi-minor axes of ${\sigma }_{\mathit{ij}}^{\mathrm{major}}$ and ${\sigma }_{\mathit{ij}}^{\mathrm{minor}}$, where the axis ${\sigma }_{\mathit{ij}}^{\mathrm{major}}$ is rotated counter-clockwise from the positive y-axis by the position angle ${\varphi }_{\mathit{ij}}^{\mathrm{im}}$, the matrices are expressed by

$$\begin{array}{c}\hfill \mathsf{R}=\left(\begin{array}{cc}cos{\varphi }_{\mathit{ij}}^{\mathrm{im}}& -sin{\varphi }_{\mathit{ij}}^{\mathrm{im}}\\ sin{\varphi }_{\mathit{ij}}^{\mathrm{im}}& cos{\varphi }_{\mathit{ij}}^{\mathrm{im}}\end{array}\right),\end{array}$$(16)

and

$$\begin{array}{c}\hfill {\mathsf{C}}_{\mathrm{im}}=\left(\begin{array}{cc}{\left({\sigma }_{\mathit{ij}}^{\mathrm{minor}}\right)}^2& 0\\ 0& {\left({\sigma }_{\mathit{ij}}^{\mathrm{major}}\right)}^2\end{array}\right).\end{array}$$(17)

For the special case where we have circular positional uncertainties with ${\sigma }_{\mathit{ij}}^{\mathrm{major}}={\sigma }_{\mathit{ij}}^{\mathrm{minor}}={\sigma }_{\mathit{ij}}^{\text{circ}}$, the matrix product $\mathsf{R}{{\mathsf{C}}_{\mathsf{im}}}^{-1}{\mathsf{R}}^{\mathrm{T}}$ is the identity matrix multiplied by ${\sigma }_{\mathit{ij}}^{-2}$, and Eq. (14) reduces to the familiar

$$\begin{array}{c}\hfill {\chi }_{\mathrm{im},\mathrm{circ}}^2={\displaystyle \sum _{j=1}^{N_{\mathrm{sys}}}\sum _{i=1}^{N_{\mathrm{im}}^j}}\frac{|\mathrm{\Delta }{\mathit{\theta }}_{\mathit{ij}}{|}^2}{{\sigma }_{\mathit{ij}}^2}·\end{array}$$(18)

To sample the parameter space, we initially used EMCEE, an ensemble sampler using parallel runs of Markov chain Monte Carlo (MCMC) chains (Foreman-Mackey et al. 2013) to obtain a sampling covariance matrix for more efficient sampling and optimization with simulated annealing (Kirkpatrick et al. 1983) and Metropolis-Hastings MCMC chains.

For each mass model in Table B.1, we ran multiple iterations of Metropolis-Hastings MCMC sampling. At each iteration, we updated the sampling covariance matrix based on the posterior distribution of the previous iteration. We adjusted the step size of our MCMC sampling to obtain an acceptance rate of ∼0.25, following Dunkley et al. (2005) for efficient sampling. In our final iteration, we ran two MCMC chains of 2 × 10⁶ steps each, where the two chains had different random number seeds. The two separate chains enabled us to test the statistical accuracy of our results.

In total, we have 23 multiple image positions from eight systems in Table 4, giving us 46 constraints from both the x and y directions. Furthermore, we have three constraints on the Einstein radii of three galaxies (BCG, cluster member with ID_phot of 116, and background LOS galaxy) in our mass model based on their measured velocity dispersions; the Einstein radii of these three galaxies cannot be reconstructed well from the multiple image positions and are mostly constrained by their measured velocity dispersions. We have a prior on the reference velocity dispersion in the Faber-Jackson relation, but we do not count this as an additional constraint since the prior range is large and this reference velocity dispersion is determined mostly by the multiple image positions. The number of mass parameters, N_par, for each of our mass models ranges from 25 to 27, as tabulated in Table B.1. The eight multiple image systems additionally require eight source positions, thus 16 x and y values. Therefore, the number of degrees of freedom (D.O.F) is D.O.F = 46 + 3 − N_par − 16 = 33 − N_par, which ranges from 6 to 8, depending on the mass model.

The SN host galaxy is lensed into giant arcs covering thousands of image intensity pixels, which provide even more constraints on the cluster mass model. However, modelling such giant arcs is computationally challenging (e.g. Acebron et al. 2024), and we defer this to future work.

6. Modelling results and predictions for fixed background cosmology

In this section we present the mass modelling results for the flat ΛCDM cosmological model with $H_0=70\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ and Ω_m = 1 − Ω_Λ = 0.3. By fixing the background cosmological model, we can quantify uncertainties associated with the mass modelling before using our mass model for cosmography.

We first present the goodness of fit of the seven GLEE mass models and their generic properties in terms of the cluster mass distribution in Sect. 6.1. We then present our ultimate model in Sect. 6.2 as a combination of multiple mass model parameterizations that fit the data comparably well. In Sect. 6.3, we present our baseline model for direct model comparison with other modelling teams. In Sect. 6.4, we predict the image positions and magnifications for both SN Encore and SN Requiem from our ultimate model. We refrain from presenting the predicted time delays of SN Encore from our model in this paper in order to keep our cosmological analysis blinded. By keeping our time-delay predictions hidden at this stage from the collaborators (Pierel et al., in prep.) who are currently measuring the time delays, their time-delay measurements will not be influenced by our model predictions. The time-delay predictions will be presented in Suyu et al. (in prep.), which will remain hidden from the time-delay analysis group until the time-delay measurements are finalized. In Sect. 6.5, we present the convergence and shear values at the SN Encore image positions for future microlensing studies.

6.1. Goodness of fit of the mass models

Fig. 8. Critical curves of the multiple models for a source at redshift z = 1.949, which is the redshift of SN Encore (and SN Requiem). The panels show the iso_halo model in blue (all panels), PL_halo in red (left panel), iso_halo+sheet in green (left panel), iso_halo+ell_BCG in brown (middle pannel), single_plane in magenta (right panel), and baseline in orange (right panel). The best-fit model was used for each curve. The critical curves of the different models are overall similar, especially near the southern giant arc and the radial arc (of the lensed host galaxy of SN Encore), and hence the overlapping critical lines.

In Table B.1, we list the χ² value of the most probable (highest posterior) mass distribution of each model in fitting to the multiple image positions. All models apart from the NFW_halo have the most-probable χ_im, MP² ≲ 7. Given that the corresponding D.O.F are between 6 and 8, we achieved a reduced χ² ∼ 1 for these models, showing a good fit to the observables. The NFW_halo model has its most-probable χ_im, MP² of 23 that is substantially higher than all other models, demonstrating that the cluster main dark matter halo is better represented by a cored isothermal profile. Given this result, we did not further consider the NFW_halo model in this work.

When comparing the mass model parameter constraints and especially the predicted SN image positions, magnifications and time delays between the two MCMC chains, we find good agreement, with the model predicted values typically agreeing within 0.5σ. Given the large number of mass model parameters, including a few parameters such as the truncation radii that are poorly constrained, the sampling of the full posterior distribution needs a long MCMC chain. This is the reason why the two chains of 2 × 10⁶ in length do not provide perfectly matching parameter constraints.

Focusing specifically on the predicted time delay between images SN1a and SN1b of SN Encore (see Fig. 7) for cosmographic analysis based on this time-delay pair (with anticipated time-delay measurement in the near future), the offset of the median time delay values between the two chains is within ∼1% for half of the six models (not considering the NFW_halo model), although the offset could reach up to ∼4% in the cases of the PL_halo and baseline models. Since we anticipate that the upcoming time-delay measurement uncertainty will be ≳5% based on existing data, our parameter uncertainties of ≲4% is currently tolerable. In contrast to the shortest time delay between SN1a and SN1b, the longer time delays of the other lensed images have more accurate predictions; all models have < 1% difference between the two chains for the delay pair of SN1a and (future) SN1d. Similarly, the model predictions are accurate to within 1% for the corresponding pair (2d–2a) of SN Requiem images. Therefore, with a future time-delay measurement of the reappearance of SN Requiem or SN Encore, the uncertainty due to statistical fluctuations in the sampling of our model is subpercent.

Given the overall consistency in the two chains for each mass model, we picked one of the two chains to present our results. All six models (i.e., excluding the NFW_halo model that does not fit well to the observed data) show similar global properties for the galaxy cluster, with a concentrated primary DM halo centred within 3″ of the BCG and an extended perturbative DM halo that is located ∼30″ southwest of the BCG in a region with cluster members (Fig. 4). The perturbative DM halo has a large core radius of ∼30″ for the four multi-lens-plane models, whereas it has a core radius of ∼20″ for the single_plane model and ∼50″ for the baseline model. The cluster halo is dominated by the primary halo component, and the perturbative halo adds more mass near the region of the giant tangential arcs, making the overall cluster halo non-elliptical. This finding is in agreement with the independent mass model of Acebron et al. (2025), who showed that the perturbation likely originates from the mass concentration located ∼140″ southeast of the BCG. The need for such additional perturbative halo components to go beyond elliptical haloes (which may be overly simplistic) and fit strong-lensing constraints has also been found in studies of other clusters (e.g. Limousin et al. 2025).

The modelled source positions of image systems 4.1 and 4.3 (that are identified from the MUSE data, see Table 4 and Fig. 7) of the most probable models are closely located on the source plane of the O II emitter. For the iso_halo+sheet model, the source positions of systems 4.1 and 4.3 are located within 0.11″; for the other five models (iso_halo, PL_halo, iso_halo+ell_BCG, single_plane, baseline), the source positions are within 0.05″. Therefore, even though systems 4.1 and 4.3 were treated as two separate image systems of four and two multiple images, respectively, they are in fact likely a single source system with six images. The treatment of a 6-image system as two separate 4-image and 2-image systems would provide two fewer constraints on the mass model, but would not lead to biases on the mass model parameter values and H₀. A more detailed study of the systems 4.1 and 4.3 is deferred to future work.

6.2. Ultimate model

The four multi-plane lens models, iso_halo, PL_halo, iso_halo+sheet, and iso_halo+ell_BCG, all fit well the observed image positions with an image χ_im, MP² ≤ 6.3, and a reduced χ_im, MP² ∼ 1. The different parameterizations enable us to incorporate systematic uncertainties from model assumptions. Despite the different mass parameterizations and even considering a model that incorporates explicitly a mass sheet to account for the mass-sheet degeneracy (which is expected to be the dominant mass-modelling uncertainty), we obtained globally similar cluster mass properties. The critical curves for the different models shown in Fig. 8 (left and middle panels) are nearly identical for these four mass models. The predicted magnifications and time delays of both SNe Requiem and Encore from the four different models mostly agree within their estimated 1σ uncertainties, although for a few cases the magnification and delay predictions from the different models (obtained from samples of the final MCMC chain, which we describe in more detail in Sect. 6.4) are more discrepant but still within 2σ.

We list the parameters of each model and their inferred values in Table C.1. The parameter values for the galaxies (cluster members, foreground and background galaxies) agree within their 1σ uncertainties across the four models (iso_halo, PL_halo, iso_halo+sheet and iso_halo+ell_BCG in columns 4–7). The reference galaxy has a lower velocity dispersion (∼170 − 180 km s⁻¹) compared to the prior (206 ± 25 km s⁻¹), indicating that the cluster mass models require on average slightly lower velocity dispersions for the cluster members compared to the Faber-Jackson fit, although the preferred values from the lens mass model are consistent within ∼1σ of the Faber-Jackson relation. In the case of the iso_halo+ell_BCG where the BCG is elliptical (instead of circular as in other models), the axis ratio of the BCG is ∼0.9 which is nearly circular, and thus its position angle has broad uncertainties. This shows that a circular BCG mass distribution is a good approximation for fitting to the current data of multiple image positions. The primary DM halo parameters also agree well across the four different models, and its position angle is consistent with the observed position angle of the BCG light profile. The PL_halo model with a radially variable DM density profile has a resulting radial profile slope of γ = 0.51 ± 0.03, which agrees well with the other isothermal halo models (with γ = 0.5) within the uncertainties. The centroid of the primary DM halo is also closely aligned with that of the BCG (within ∼1.5″). In contrast, the perturbative DM halo centroid is ∼30″ southwest of the BCG, and has a large core radius. Furthermore, there is degeneracy in the perturbative halo component and the mass sheet in the iso_halo+sheet model, where the slightly positive convergence sheet of ∼0.09 contributes to the extended DM halo, resulting in a smaller Einstein radius of the perturbative DM halo compared to the other three mass models. In summary, Table C.1 shows a high level of consistency across the different models, despite the different mass-model parameterizations.

As an illustration, we show a breakdown of the contributions to χ_im, MP² from each multiple image system for the iso_halo model. For the SN Encore (system 1), SN Requiem (system 2), and their host galaxy centroid (system 3), we are able to reproduce their positions within a root-mean-square (rms) scatter of ≲0.03″ for each system, which is comparable to the positional uncertainties. We are also able to reproduce the image positions of system 4.2 ([O II] arc with positions from the JWST image) with a rms of 0.02″. The remaining image systems have larger rms offsets between the predicted and observed image positions, ranging from 0.21″ to 0.44″, because of the larger observed positional uncertainties from the MUSE data (as listed in Table 4). Overall, we fit all 23 identified multiple image positions within their observational uncertainties and with an rms scatter of 0.24″ for the iso_halo model. The results for the PL_halo, iso_halo+sheet and iso_halo+ell_BCG models are similar, given their similar values of χ_im, MP².

In Fig. 9, we show the effective convergence (dimensionless surface mass density) for the SN redshift (z = 1.949) from our multi-lens-plane iso_halo model that accounts for the three lens redshift planes (z_fg = 0.309 of the foreground galaxy, z_d = 0.336 of the galaxy cluster, and z_bg = 0.371 of the background galaxy). The effective convergence, κ_eff, is obtained by taking $\frac{1}{2}\mathbf{\nabla }·{\mathit{\alpha }}_{\mathrm{tot}}$, where α_tot is the total (scaled) deflection angle from summing up the deflection angles $\widehat{\mathit{\alpha }}$ in all lens planes to the z = 1.949 source plane: ${\mathit{\alpha }}_{\mathrm{tot}}={\displaystyle \sum _{i=1}^{N-1}}\frac{D_{\mathit{iN}}}{D_N}{\widehat{\mathit{\alpha }}\mathbf{(}\mathit{\theta }}_i)$ with the source plane as the N^th plane (e.g. Gavazzi et al. 2008; Chirivì et al. 2018). For our specific case with three lens redshift planes, the source plane is thus at N = 4. The contours of the effective convergence show that most of the mass is in the primary dark matter halo given the nearly concentric elliptical contours. The contours broaden slightly beyond the tangential arcs in the southwest direction from the BCG, due to the presence of the perturbative dark matter halo. Despite the seemingly small amount of convergence associated with the perturbative halo, its presence is necessary to fit to the observed image positions with a reduced χ_im² ∼ 1. Without this perturbative halo, the χ_im² was an order of magnitude higher.

Fig. 9. Effective convergence, κeff, of the iso_halo model for the SN Encore redshift (z = 1.949). The total κeff distribution (in grey) near the cluster core (inner ∼10″ region) is dominated by the primary DM halo and the BCG, and κeff follows roughly the observed BCG and galaxy intensity distribution. The extension of the contour level with κeff = 0.6 towards the southwest part of the BCG is due to the presence of the perturbative DM halo located around (xDM2, yDM2)∼(20″, −20″), for which we show its position in blue with its 1σ uncertainty. The green contours show the κeff of only the main DM halo, the red contours show κeff for the main and perturbative halo combined, which highlights the deformation of the overall κeff due to the perturbative halo.

We obtained our ultimate model based on the lensed image positions by combining these four multi-lens-plane models (iso_halo, PL_halo, iso_halo+sheet, and iso_halo+ell_BCG) with equal weights. Our ultimate model therefore incorporates systematic uncertainties due to mass model parameterizations in predicting the properties (e.g. positions, magnifications and time delays) of SN Encore and SN Requiem, and also for measuring cosmological parameters in future studies.

6.3. single_plane and baseline models

We consider a scenario where we place the LOS galaxies at the same redshift as the galaxy cluster. This will enable us to assess the impact of assuming a single lens plane in the mass modelling of the MACS J0138−2155 cluster, especially since most of the cluster-lens models to date adopt a single lens plane. Furthermore, the positional uncertainties of multiple lensed images are often considered to be circular rather than elliptical. Circular uncertainties are sufficient approximations especially when the rms offset between the predicted and observed image positions is substantially larger than the nominal positional uncertainties, signaling that additional origins of uncertainties associated with mass model assumptions dominate over the positional uncertainty. Since several of our mass models in Table B.1 fit well to the observed image positions within the elliptical positional uncertainties, such that the usages of elliptical versus circular uncertainties could lead to differences, we assess here the difference in modelling results between elliptical and circular positional uncertainties.

We thus consider two specific models: (1) the single_plane model, which takes the iso_halo model and assumes the foreground galaxy at z_fg = 0.309 and background galaxy at z_bg = 0.371 are both at the cluster redshift z_d = 0.336 instead; and (2) the baseline model, which is the same as the single_plane model but with circular positional uncertainties σ^circ that is obtained as the geometric mean of the semi-major and semi-minor axes of the elliptical positional uncertainty, i.e., ${\sigma }^{\mathrm{circ}}=\sqrt{{\sigma }^{\mathrm{major}}{\sigma }^{\mathrm{minor}}}$. As noted in Sect. 5.2, the baseline model with a single lens plane and circular positional uncertainties is the setup that all the seven independent modelling teams in Suyu et al. (in prep.) have agreed to construct in order to facilitate a direct model comparison.

We show the critical curves of the single_plane (magenta) and baseline (orange) models in comparison to that of the iso_halo (blue) model in the right-hand panel of Fig. 8. The critical curves again overlap well, even at the locations near the foreground galaxy (at (x_fg, y_fg)∼(−1″, −17″)) and the background galaxy (at (x_bkg, y_bkg)∼(7″, 4″)) that are approximated at the cluster redshift in the single_plane and baseline models. There is a noticeable difference in the critical curves near the jellyfish galaxy at (x_JF1, y_JF1)∼(19″, −13″) where the single_plane and baseline models have a lower resulting Einstein radius for this jellyfish galaxy. Nonetheless, since the multiple lensed image positions and the giant arcs are not near this critical curve, the difference in the critical curve is expected to have minimal impact on the model properties of SN Encore and SN Requiem.

We tabulate the model parameter values of the baseline model in Table C.1 in the last column. In comparison to the values of iso_halo, we see that almost all the parameter values agree within their 1σ uncertainties, except for the Einstein radius of the first jellyfish galaxy that is almost ∼2σ lower in the baseline model. We find that adopting single-plane lensing with elliptical positional uncertainties, i.e. the single_plane model, shifts the predicted magnification and time delays by ∼1σ for most predictions, and up to ∼2σ for a few predictions, relative to the iso_halo model. Curiously, the change of elliptical to circular positional uncertainties in the baseline reduces the shift relative to the iso_halo to within ∼1σ. While it is merely a coincidence that the biases induced by the two assumptions of single-plane lens and circular-positional uncertainties cancel each other to some extent, this finding suggests that the baseline model provides an approximation of the galaxy cluster MACS J0138−2155 mass distribution that is accurate to ∼1σ of the mass modelling uncertainties in comparison to our ultimate model in Sect. 6.2. In Table 5, we show in columns 4 and 5 the breakdown of χ_im, MP² and rms scatter for the most-probable baseline model. The values are similar to those of the iso_halo model, with the same resulting rms scatter of 0.24″ between all the observed and predicted image positions.

6.4. Predicted positions and magnifications of SN Encore and SN Requiem

To predict the multiple image positions of a given model, we first map the individual observed image positions back to the source plane through the multi-plane lens equation. We then obtain the model source position as the magnification-weighted average of the mapped source positions. Given the model source position, we then solve the non-linear lens equation to obtain the model-predicted image positions. We also compute the lensing magnifications at the model-predicted image positions.

For our ultimate model, which is the combination of the four multi-plane mass models, we combined the respective chains, each of length 2 × 10⁶, to obtain a final chain length of 8 × 10⁶. We then thinned this chain by a factor of 1000, reaching a final thinned chain of 8000 samples. For each sample in this chain, we then solved the lens equation to obtain the predicted image positions and magnifications. For the baseline model, we similarly thinned it by a factor of 1000 and use the thinned chain to predict the image positions/magnifications.

We find that all samples in the chains predict at least four images for either SN Encore and SN Requiem. For SN Encore, the four predicted images are near the 1a, 1b, 1c and 1d locations marked on Fig. 1. For SN Requiem, the four images that are persistently predicted by the models are similar in configuration to that of SN Encore; we denote them as images 2a, 2b, 2c (matching their corresponding observed images) and 2d. Images 1d and 2d are predicted near the radial arc, and will appear in the future given their model time delays of multiple years relative to the first images 1a and 2a, respectively. While some mass models in the MCMC chain predict a fifth central image of SN Encore or SN Requiem near the core of the BCG, this image is not always predicted. We denote these images as 1e for SN Encore and 2e for SN Requiem, and compute the percentage of models in the chains that predict these fifth images. Furthermore, some samples in the chain produce additional (future) images of SNe Encore and Requiem near the jellyfish galaxies along the radial arc. Similarly, some samples in the chain produce additional images of SN Encore formed by the foreground galaxy (at (x_fg, y_fg) listed in Table C.1) near image 1b. Since these additional images are not persistently predicted across all samples in the chain and are sensitive to the mass distributions of the jellyfish galaxies or the foreground galaxy, which are simplistically approximated as spherical in our model, we defer to future work for investigating in detail these non-persistent images. Future detections of such additional images from the jellyfish galaxies, if any, will help place additional constraints on the mass distribution of inner cluster core and the jellyfish galaxies.

In Table 6 we list the predicted image positions and magnifications of both SN Encore and SN Requiem based on our ultimate model. The predicted positions of 1a, 1b, 2a, 2b, and 2c agree well with the observed positions in Table 4 within both the observational uncertaintes in Table 4 and the model uncertainties in Table 6, which are in line with the image position rms and χ_im² computed in Table 5. The model positional uncertainties on these image positions are on the order of tens of milliarcseconds, as shown in Table 6. The model positional uncertainties of undetected/future images (1c, 1d, 1e, 2d, 2e) are larger, up to ∼0.7″.

We list the predictions from the baseline model in Table D.1. The predictions from the ultimate model and the baseline models agree well overall within ∼1σ of the model uncertainties. This shows that the baseline model is a sufficient approximation for modelling SN Encore and SN Requiem, based on the current image positions and kinematic constraints. Nonetheless, we caution that there is a larger shift in the model predictions when comparing the single_plane model (with elliptical positional uncertainties) to the ultimate model. Therefore, to err on the cautious side, we advocate for using the ultimate model from multi-plane lensing for future cosmographic analysis, particularly in measuring the value of H₀.

The magnification values of the multiple images predicted from both our ultimate and baseline models are substantially higher than the values predicted by Newman et al. (2018a) and Rodney et al. (2021) for the SN host galaxy and SN Requiem, by a factor up to ∼6. In contrast, our predicted magnifications of SN Requiem agree with those of Acebron et al. (2025) within the estimated uncertainties, and we refer to Acebron et al. (2025) for a more detailed discussion.

6.5. Convergence and shear of SN Encore

To enable future microlensing studies of SN Encore, we list in Table 7 the effective convergence and shear values at the positions of SN Encore multiple images from the four multi-plane models that constituted the ultimate model. For the images 1a and 1b, we evaluated the convergence/shear at the observed image positions (in Table 4); for image 1c, that is not clearly visible, and the future image 1d, we used the median of the predicted image positions of the ultimate model (in Table 6). Since images 1e is not always predicted to be present in our mass models, we do not list its convergence and shear here.

The most probable convergence and shear values from the four models mostly agree within 0.02 and all within 0.04, showing broad consistency in the results from the different mass model parameterizations. For microlensing studies, the differences in convergence and shear across the models are small in comparison to the less well-constrained smooth matter mass fraction (which is the convergence due to the smooth DM distribution in fraction of the total convergence composed of DM and baryons) and also the stochasticity of stellar distributions. The use of any of the convergence and shear values listed in Table 7 for microlensing studies should therefore suffice. Microlensing can be significant especially for image 1b, that has the foreground galaxy located within 2″, and image 1d, that has a jellyfish galaxy stripped across the neighbouring area. Future microlensing studies of this system should account for the effects of these nearby galaxies.

7. Summary

We constructed new lensing models of the galaxy cluster MACS J0138−2155 based on new JWST imaging and MUSE spectroscopy. The analysis of the spectroscopic dataset is performed by Granata et al. (2025), which we combined with the photometric analysis here to obtain the key ingredients for building our cluster total mass models. We summarize the main results as follows.

• We fitted nearly 500 galaxies in the FOV of the HST and JWST imaging data using Sérsic profiles and we obtained their photometries across six JWST and five HST filters.

• We used these galaxy photometries to select 34 additional cluster members based on colour-magnitude and colour-colour criteria, which complemented the spectroscopic sample of 50 cluster members and resulted in a total of 84 cluster members. The photometric catalogue of the cluster members is publicly released along with this paper.

• In addition to the cluster, there is one foreground galaxy (near the southern giant arc) and one massive background galaxy, which we included explicitly in our multi-lens-plane mass model.

• We identified the image positions of eight multiple image systems with a total of 23 images from four distinct spectroscopic redshifts in a range of 0.767 ≤ z ≤ 3.420, based on the combination of JWST, HST, and MUSE data.

• We considered seven different lens mass models with GLEE, exploring a range of DM halo profiles, the inclusion of a mass sheet, the impact of multi-plane versus single-plane modelling, and the effect of having elliptical versus circular BCG. With two DM halo components (one main DM halo and a second, perturbative component) and most of the cluster galaxies that follow the Faber-Jackson relation, six of the seven models fit to the observed image positions with a reduced χ_im² ∼ 1. Cored-isothermal DM haloes could fit well to the observables, whereas the NFW DM halo results in a reduced χ_im² that is approximately four times higher.

• We combined the four multi-lens-plane mass models (iso_halo, PL_halo, iso_halo+sheet and iso_halo+ell_BCG, which have remarkably consistent results) to form our ultimate model for predicting observables of SN Encore and SN Requiem, and also for future cosmographic analyses.

• We considered a baseline model with the two LOS galaxies approximated at the cluster redshift and with circularized positional uncertainties of the image positions. The baseline model is designed to facilitate the direct comparison of independent mass models from different teams (Suyu et al., in prep.). We found the baseline model provides a good approximation to our ultimate model.

• We predicted the positions and magnifications of SN Encore and SN Requiem, including ones that will only appear in the future. Observations of these predicted images in the future will provide a unique opportunity to test our models.

• We provide the convergence and shear values at the positions of the SN Encore multiple images to enable future microlensing studies.

The exciting discovery of SN Encore in the galaxy that previously hosted SN Requiem gives us an excellent opportunity to study cosmology, particularly to constrain the Hubble constant, and the SN host galaxy. Our work in building an accurate total mass model of the cluster is crucial as it lays the foundation for enabling state-of-the-art cosmological and astrophysical studies of this cluster. The dataset presented in this paper also forms the basis for independent cluster-mass modelling of this cluster. The comparison of seven independent mass models based on this unique dataset will be presented in a forthcoming publication.

Data availability

Table A.1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/702/A157

Acknowledgments

We thank Jake Summers and Jordan D’Silva for their help in reducing the JWST images, and the anonymous referee for the helpful feedback. This paper is based in part on observations with the NASA/ESA Hubble Space Telescope and James Webb Space Telescope obtained from the Mikulski Archive for Space Telescopes at STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555 for HST, and under NASA contract NAS 5-03127 for JWST. The specific observations used in this work can be accessed via DOI:https://doi.org/10.17909/snj9-an10. SE and SHS thank the Max Planck Society for support through the Max Planck Fellowship for SHS. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (LENSNOVA: grant agreement No 771776). This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. SS has received funding from the European Union’s Horizon 2022 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 101105167 – FASTIDIoUS. We acknowledge financial support through grants PRIN-MIUR 2017WSCC32 and 2020SKSTHZ. AA acknowledges financial support through the Beatriz Galindo programme. AA and JMD acknowledge support from the project PID2022-138896NB-C51 (MCIU/AEI/MINECO/FEDER, UE), Ministerio de Ciencia, Investigación y Universidades. YF is supported by JSPS KAKENHI Grant Numbers JP22K21349 and JP23K13149. MJJ and SC acknowledge support for the current research from the National Research Foundation (NRF) of Korea under the programs 2022R1A2C1003130, RS-2023-00219959, and RS-2024-00413036. Support was provided to JDRP through program HST-GO-16264. JDRP is supported by NASA through an Einstein Fellowship grant No. HF2-51541.001, awarded by STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. AZ acknowledges support by Grant No. 2020750 from the United States-Israel Binational Science Foundation (BSF) and Grant No. 2109066 from the United States National Science Foundation (NSF), and by the Israel Science Foundation Grant No. 864/23.

References

Appendix A: Measured photometry of cluster members in MACS J0138−2155

In Table A.1, we list the photometric ID, the position in RA and dec, and the magnitudes in all filters listed in Table of the objects included in the model. The full table is available at the CDS. The BCG is labeled as index 0. If an object is not modeled in a specific band (e.g. if it is not in the field of view) its magnitude is marked as −. The reported magnitudes are AB magnitudes.

Appendix B: Lens models

In Table B.1, we provide an overview of the lens mass models. We describe the different models and list their number of mass parameters (N_par), as well as their χ_im, MP² values.

Appendix C: Model parameters

In Table C.1, we list the final model parameter values for the four multi-plane lens models (iso_halo, PL_halo, iso_halo+sheet and iso_halo+ell_BCG) and the single-plane baseline model.

Appendix D: Predicted image positions and magnifications of supernovae from the baseline model

In Table D.1 we list the predicted image positions and magnifications of SN Encore and SN Requiem from the baseline model where the LOS galaxies are approximated to be at the cluster redshift plane, and positional uncertainties of lensed images are circularized (see Sect. 6.3 for more details). Comparing to the ultimate model predictions in Table 6 from multi-plane lens modelling, the predictions from the baseline model are mostly consistent within the 1σ modelling uncertainties. This shows that the baseline model provides a good approximation to the final best model for the current sets of image positions and kinematic data.

All Tables

All Figures

Fig. 1. Colour image of MACS J0138−2155 from combining JWST/NIRCam and HST/WFC3 data. The RGB image is constructed from the filters F105W+F115W+F125W (blue), F150W+F160W+F200W (green), and F277W+F356W+F444W (red). The observed images of SN Encore are shown in solid white circles, the expected positions of images D and possibly E are in dashed white circles. The observed and predicted images for SN Requiem are shown in blue following the same pattern. The images were drizzled to a pixel scale of 0.04″/pix. (Image Credit: STScI, A. Koekemoer, T. Li).

In the text

	Fig. 2. Results from modelling the BCG light with Sérsic profiles in the HST F160W filter. The left panel shows the observed image, the middle panel shows our BCG light model, and the right panel is the BCG subtracted image. The images cover the central 50″ × 50″ around the BCG.
In the text

	Fig. 3. Examples for Sérsic fits to objects in the field in the JWST F200W filter. From left to right: data image, model image, and residual image. The photometric IDs of the objects are shown in blue on the data image. Object 82 is a bright background galaxy at a redshift zbg = 0.371, and the other galaxies are cluster members close to the BCG.
In the text

Fig. 4. Colour image of MACS J0138−2155 built from JWST NIRCam filters (red: F277W, F356W, F444W; green: F150W, F200W; blue: F115W). Spectroscopically confirmed objects are marked with squares, coloured in cyan for cluster members, and blue and yellow for the foreground and background objects, respectively. Objects selected from the colour-magnitude fit are marked as red circles, and green triangles are the colour-colour selections.

In the text

Fig. 5. Colour-magnitude diagram in the field of MACS J0138−2155. Spectroscopically confirmed cluster members and non-members are indicated as blue and red points, respectively. Sources with no spectroscopic confirmation are shown in grey. The blue line and region indicate the fitted red sequence and the 1σ scatter. Red crosses (X) mark cluster members that do not follow the red sequence and yellow plus signs indicate the photometrically selected members that are included in our lens model. Objects inside the red sequence with a red minus sign are not included as photometric cluster members due to the lack of, or unreliable, F160W photometry.

In the text

Fig. 6. Colour-colour diagram in the field of MACS J0138−2155, where the right panel is a zoom-in of the left panel. Spectroscopically confirmed cluster members and non-members are indicated in blue and red, respectively. Sources with no spectroscopic confirmation are shown in grey, with size and shade related to the galaxy brightness as shown in the colour bar. Red crosses (X) mark spectroscopic cluster members that do not follow the red sequence and yellow plus signs indicate the colour-magnitude selected members (from Fig. 5) that are included in our lens model. Photometric galaxies inside the blue shaded region in the right-hand panel are selected as cluster members via our colour-colour selection.

In the text

	Fig. 7. Positions of the eight multiple image systems (the ‘gold’ systems) used to constrain our mass model. We show the positions in arcseconds relative to the BCG at (RA, Dec) = (24.5157084°, −21.9254717°).
In the text

Fig. 8. Critical curves of the multiple models for a source at redshift z = 1.949, which is the redshift of SN Encore (and SN Requiem). The panels show the iso_halo model in blue (all panels), PL_halo in red (left panel), iso_halo+sheet in green (left panel), iso_halo+ell_BCG in brown (middle pannel), single_plane in magenta (right panel), and baseline in orange (right panel). The best-fit model was used for each curve. The critical curves of the different models are overall similar, especially near the southern giant arc and the radial arc (of the lensed host galaxy of SN Encore), and hence the overlapping critical lines.

In the text

Fig. 9. Effective convergence, κeff, of the iso_halo model for the SN Encore redshift (z = 1.949). The total κeff distribution (in grey) near the cluster core (inner ∼10″ region) is dominated by the primary DM halo and the BCG, and κeff follows roughly the observed BCG and galaxy intensity distribution. The extension of the contour level with κeff = 0.6 towards the southwest part of the BCG is due to the presence of the perturbative DM halo located around (xDM2, yDM2)∼(20″, −20″), for which we show its position in blue with its 1σ uncertainty. The green contours show the κeff of only the main DM halo, the red contours show κeff for the main and perturbative halo combined, which highlights the deformation of the overall κeff due to the perturbative halo.

In the text