© The Authors 2025

1. Introduction

Strong gravitational lensing of cluster-sized lenses has become an important tool in astrophysics. It provides large magnifications that allow one to study faint galaxies and sub-kiloparsec regions at high redshifts that would otherwise be inaccessible with current telescopes. It thus offers an insight into the formation of first galaxies and the sources of reionisation (e.g. Hashimoto et al. 2018; Vanzella et al. 2023; Welch et al. 2023; Asada et al. 2023; Strait et al. 2023; Mowla et al. 2024; Adamo et al. 2024; Bradač et al. 2024; Fujimoto et al. 2024). In extreme cases, it even enables studies of individual stars over cosmological distances (e.g. Welch et al. 2022; Meena et al. 2023; Fudamoto et al. 2025; Furtak et al. 2024). Since gravitational lensing allows one to estimate mass distribution independently of the hydrostatic equilibrium and dynamical equilibrium of cluster galaxies, it is also an essential tool for investigating the interplay between baryons and dark matter on different scales and characterising properties of merging and out-of-equilibrium clusters (see Natarajan et al. 2024 for a review). Finally, since lensing depends on the geometry of the Universe, it can also be utilised to constrain values of cosmological parameters (e.g. Refsdal 1964; Jullo et al. 2010; Magaña et al. 2018; Caminha et al. 2022a).

The various applications of strong lensing rely on an accurate reconstruction of the mass distribution, which in turn relies on the quality and abundance of data, namely, positions and accurate redshifts of multiple images of strongly lensed background sources. This first requires high-resolution imaging in several photometric bands, which enable precise position measurements and multiple-image system identification based on colours and morphologies. It also requires spectroscopic redshift measurements, which are essential for confirming multiple images and reducing systematic errors (e.g. Johnson & Sharon 2016). Until recently, the progress in the field of strong lensing in galaxy clusters was spearheaded by the Hubble Space Telescope (HST) through its large field of view and high-resolution imaging in the visible wavelengths and complemented by ground-based spectroscopic observations, for example, with Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2012) at the Very Large Telescope (VLT). The James Webb Space Telescope (JWST), launched in 2021, expanded those capabilities to near-infrared wavelengths, paving the way for a new generation of lens models (e.g. Caminha et al. 2022b; Mahler et al. 2023; Furtak et al. 2023; Diego et al. 2024a; Gledhill et al. 2024).

The CAnadian NIRISS Unbiased Cluster Survey (CANUCS) is a ∼200 h survey that aims to use the unprecedented capabilities of JWST in combination with strong gravitational lensing in galaxy clusters to study high redshift galaxies. With three of the four instruments on board JWST, CANUCS observed five massive clusters known to be effective gravitational lenses, including MACS J0416.1-2403 (hereafter MACS0416).

MACS0416 is a massive galaxy cluster at z = 0.397 discovered by the MAssive Cluster Survey (MACS; Ebeling et al. 2010). It is a cluster merger featuring a two-peaked X-ray surface brightness distribution (Mann & Ebeling 2012) coinciding with a bimodal mass distribution and is likely in a pre-collisional phase (Balestra et al. 2016). Due to its elongated mass distribution and large Einstein radius, MACS0416 is a very effective gravitational lens with an extraordinary number of multiple images (e.g. Zitrin et al. 2013; Balestra et al. 2016; Hoag et al. 2016; Caminha et al. 2017; Diego et al. 2024a). It has also sparked interest due to its abundance of high magnification transient events, enabling the studies of stars at z ∼ 1 (Rodney et al. 2018; Chen et al. 2019; Kaurov et al. 2019; Yan et al. 2023; Diego et al. 2024a). Since its discovery, the cluster has been targeted by several cluster surveys. It was selected as a part of the Cluster Lensing and Supernova survey with Hubble (CLASH; Postman et al. 2012), which provided constraints for the first detailed mass model by Zitrin et al. (2013). Due to its essential lensing properties, MACS0416 was among the six clusters observed as a part of the Hubble Frontier Fields (HFF) program (Lotz et al. 2017, 2014), which provided deep (≲28.5 mag) imaging data in seven HST ACS/WFC3 bands of the central cluster region. The HFF imaging data has since been used in numerous lensing analyses (e.g. Jauzac et al. 2014, 2015; Johnson et al. 2014; Diego et al. 2015; Kawamata et al. 2016; Sebesta et al. 2016). The extended cluster region was also observed by The Beyond Ultra-deep Frontier Fields and Legacy Observations (BUFFALO) program (Steinhardt et al. 2020), and it was included in the Flashlights program (Kelly et al. 2019).

MACS0416 has also been targeted by several spectroscopic follow-up observations. The cluster was observed with VLT’s VIsible MultiObject Spectrograph (VIMOS) as a part of the CLASH-VLT spectroscopic campaign (presented in Balestra et al. 2016), resulting in a lens model constrained with 30 spectroscopically confirmed multiple images belonging to ten systems (Grillo et al. 2015; see our Table 1). This spectroscopic dataset contributed to the development of Richard et al. (2014), Johnson et al. (2014) and Jauzac et al. (2014) lens models. The HST/WFC3 infrared grism spectroscopy obtained by the Grism Lens-Amplified Survey from Space (GLASS; Treu et al. 2015; Schmidt et al. 2014) led to a non-parametric lens model and a dataset of 30 spectroscopic multiple images of 15 background sources (Hoag et al. 2016). A substantial improvement of the multiple-image catalogue was facilitated by several VLT/MUSE observations. The southwest region of the cluster was observed with MUSE as a part of the programme 094.A-0525 (PI: Bauer) with 11 h of exposure time. The northeast region has been targeted by the programs GTO 094.A-0115B (PI: Richard) and 0100.A-0764 (PI: Vanzella, Vanzella et al. 2021), reaching a total integration time of 17.1 h. This makes MACS0416 a cluster with some of the deepest MUSE observations to date.

Leveraging the first MUSE observations, Caminha et al. (2017) presented a catalogue of 102 spectroscopically confirmed multiple images from 37 sources. Their model was improved by including a hot X-ray emitting gas component (Bonamigo et al. 2017, 2018) and information on galaxy kinematics (Bergamini et al. 2019). The dataset of multiple images was further expanded by Bergamini et al. (2021) by leveraging the data presented in Vanzella et al. (2021). Richard et al. (2021) published a MUSE catalogue of 198 multiple images from 71 systems. This number was further increased to 237 spectroscopically confirmed multiple images by Bergamini et al. (2023, hereafter Be23). Another catalogue of 214 multiple images was compiled by Diego et al. (2024b) by re-evaluating systems from existing literature. Before the onset of JWST, these catalogues presented the most extensive dataset of spectroscopically confirmed multiple images for this cluster, and they served as a starting point for constructing the initial catalogue used in our search for additional multiple-image systems. Apart from the aforementioned works, MUSE datasets have been used for several other lens reconstructions (Vega-Ferrero et al. 2019; Kawamata et al. 2018; Limousin et al. 2022; Cha & Jee 2023; Perera et al. 2025).

MACS0416 was recently observed with the JWST Near Infrared Camera (NIRCam) by the CANUCS survey and the Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS) project (Windhorst et al. 2023), covering the wavelength range from 0.8 to 5 μm. Three epochs of observations of the PEARLS program with an additional CANUCS epoch have enabled studies of transient events in highly magnified regions of lensed galaxies (e.g. Yan et al. 2023; Diego et al. 2023). During the preparation of this paper, the PEARLS collaboration published a new lens model constrained by leveraging the new JWST imaging data (Diego et al. 2024a, henceforth Di24). With 343 multiple-image candidates, belonging to 119 multiple-image systems, the Di24 catalogue represents the largest dataset of its kind to date. It contains all previously known systems with MUSE spectroscopic redshift as well as 41 additional candidates without spectroscopic confirmation. We used Di24 catalogue to supplement our dataset, which we further expanded with our imaging and spectroscopic data. For easier comparison of our work with theirs, we cross-matched the catalogues and adopted the names in Di24 of the multiple-image systems.

In addition to NIRCam imaging, the CANUCS program also includes Near Infrared Imager and Slitless Spectrograph (NIRISS; Doyon et al. 2012) wide-field slitless spectroscopy in F115W, F150W, and F200W bands as well as Near Infrared Spectrograph (NIRSpec) multi-object prism spectroscopy. The new spectroscopic data obtained with JWST are the basis of this work.

In Table 1, we list catalogues from several strong lensing studies with the number of multiple images and multiple-image systems. The combination of having a number of strong lensing studies and the largest catalogue of reliable spectroscopic multiple-image systems available (303 in this work) also make MACS0416 interesting for investigating the evolution of the lens models due to an increasing abundance of data. A recent study by Perera et al. (2024) has provided a comparison of the mass reconstructions using different datasets (including the one presented in this work).

Throughout this work we assume a flat ΛCMD cosmology with Ω_Λ = 0.7, Ω_m = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹. At the cluster redshift z = 0.396, a projected distance of 1″ corresponds to a physical scale of 5.340 kpc. Magnitudes are given in the AB system (Oke & Gunn 1983).

2. Data

2.1. Imaging and photometry

In this work we use CANUCS NIRCam observations in filters F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W with 6.4 ks exposure time per filter. In addition, we use the archival HST/ACS imaging data in F435W, F606W, and F814W filters and HST/WFC3 data in F105W, F110W, F125W, F140W, F160W filters from the HFF (Lotz et al. 2017, 2014) and CLASH (Postman et al. 2012) programs. The data reduction and production of the photometric catalogues follows the procedure outlined in Noirot et al. (2023) and Asada et al. (2024). The NIRCam and WFC3 images were first processed with the official STScI JWST pipeline and Grizli (Brammer 2023) and drizzled on the same pixel scale (40 mas/pixel) and aligned with Gaia DR3 astrometry (Gaia Collaboration 2016, 2023). The diffuse cluster light and the bright cluster galaxies were modelled and removed following the procedure described in Martis et al. (2024). The bright cluster galaxy subtracted images were then PSF-homogenised to match the PSF in the F444W filter. Source detection was performed on the χ_mean detection image (Drlica-Wagner et al. 2018). Photometric catalogues were produced with the Photutils package (Bradley et al. 2023) using Kron (Kron 1980) and fixed apertures of different sizes. Photometric redshifts were computed with EAZY-py (Brammer et al. 2008) using the latest standard templates (tweak_fsps_QSF_12_v3) and templates from Larson et al. (2023). EAZY-py also returns the “risk” value z_phot_risk, as defined in Tanaka et al. (2018), which we use as an indicator of the photometric redshift reliability. In this work, we only show photometric redshifts with z_phot_risk< 0.2 (Sect. 3.1).

2.2. NIRISS grism spectroscopy

The cluster was observed with NIRISS using wide-field slitless spectroscopy mode (Willott et al. 2022) in F115W, F150W and F200W filters and two orthogonal grisms (GR150C and GR150R, R ∼ 150) with 9.7 ks exposure time per configuration. Data reduction follows the procedure described in Noirot et al. (2023) and includes the initial processing and source modelling with Grizli. The contamination of the cluster galaxies is modelled using the isophotal models of the bright cluster galaxies (Sect. 2.1) and subtracted from the spectra (the procedure is described in Estrada-Carpenter et al. 2024). We use our Photutils catalogues derived from NIRCam images to obtain the positions of sources and their spectra. The continuum of each source is then modelled using an iterative polynomial fit. The contamination from other sources is subtracted from the spectrum of each source in the catalogue. The redshift is fitted using the procedure described in Noirot et al. (2023). Each source is fitted with three different methods. The first method fits the spectra using a narrow spectral range (1.03–1.26 μm in F115W, 1.35–1.65 μm in F150W, 1.79–2.20 μm in F200W) to avoid the less sensitive regions at the edges of each filter response curve. The second method uses the full spectral range to reliably fit sources with emission lines with wavelengths close to the edges of the filter response curves. The third method is a joint fit of the spectro-photometric data with Grizli using all available JWST NIRISS, NIRCam and HST filters. The obtained grism redshifts z_niriss are prioritised in that order and are selected upon visual inspection of each spectral fit and assessment of its quality. For some sources, none of the three fits returned a reliable redshift. Those sources were refitted with custom modification of the data range (e.g. by removing configurations with visible unsubtracted contamination) or by applying narrower redshift priors and then visually reinspected. The grism spectra used in this work are provided in Appendix B together with details of the individual fits.

While Grizli provides a redshift probability density function, its width likely underestimates the true redshift uncertainty (see Noirot et al. 2023). To account for the possible systematic errors we compare grism redshifts of known multiple images with their MUSE redshifts z_muse from Richard et al. (2021) catalogue. We find that all grism redshifts apart from the two outliers K8 and K5 (see Sect. 3.2) match MUSE redshift within δz_niriss = ±0.005(1 + z_spec) which we use as an error estimate. This error encompasses a small systematic offset – we found that our NIRISS redshifts z_niriss may be, on average, underestimated by ∼0.002(1 + z_spec) when compared with z_muse. While the issue will be further investigated before the release of the CANUCS NIRISS redshift catalogues, we note that this difference is of negligible consequence for the strong lensing analysis.

In Sect. 3.2 (see Fig. 2), we also utilise spatially resolved emission line maps to match sub-components of multiple images. Maps are obtained as described in Estrada-Carpenter et al. (2024). To this end, the extended source was first segmented based on the pixel signal-to-noise ratio. All regions were then fitted simultaneously with priors obtained with broadband photometry and by using a linear combination of spectral templates. The emission line maps were obtained by forward modelling the galaxy with the best-fit model, excluding the line of interest, subtracting it, and dithering the residuals.

2.3. NIRSpec prism spectroscopy

Following the NIRCam and NIRISS observations, their reduction and subsequent target selection, MACS0416 was observed with NIRSpec using the PRISM/CLEAR disperser with nominal resolving power ∼100 and the approximate wavelength coverage between 0.6 and 5.3 μm. Targets were selected for a variety of science goals, including multiple-image system candidates for this project (see Sect. 3.1), and were observed through the Micro-Shutter Assembly (MSA) with 2.9 ks exposures for each MSA configuration. The cluster was observed with three MSA configurations. The data is then processed using the STScI JWST pipeline and the msaexp package (Brammer 2022). This includes mitigating the snowball residuals, masking out pixels approaching saturation, standard wavelength calibration, flat-field and path-loss corrections, photometric calibration and the extraction of 1-D spectra (e.g. see Withers et al. 2023, and Desprez et al. 2024 for a more detailed outline of the procedure). The redshifts z_nirspec were obtained using msaexp. In Appendix C we show the NIRSpec spectra used in this work and describe cases that required distinct treatment.

The NIRSpec redshift uncertainty was estimated by comparing NIRSpec and MUSE redshift from Richard et al. (2021), similarly to Sect. 2.2. We found that the difference between MUSE and NIRSpec redshift is below 0.002(1 + z_spec) for more than 68% of multiple images. We use this as an error estimate, which takes into account potential small systematic deviations (e.g. Bunker et al. 2024). As a test, we also estimated the uncertainty by fitting a Gaussian profile to the strongest spectral line in spectra of several objects and obtained the line position uncertainty. We find it to be consistent with our error estimate.

3. Catalogue of multiple images

In this section we describe our catalogue of strong lensing constraints and the methods used to obtain them. Based on reliability, we divide our catalogue in several categories. The gold category, which is used for lens model optimisation, is the most secure and encompasses reliable images from existing catalogues with MUSE redshift z_muse (see Sect. 3.2) and new spectroscopic systems, which we trust due to multiple spectroscopic redshift measurements with NIRISS or NIRSpec, morphology, photometric redshift, or/and lensing configuration. It contains 303 images from 111 multiple-image systems. We also provide a catalogue of less reliable sources which may be of interest for future studies. These sources are divided into silver or bronze categories with intermediate and low degrees of confidence, respectively. We also add a fourth category of quartz images. Those images appear in previous works with measured MUSE redshifts but are excluded from our gold category because we cannot reliably identify or confirm them with our NIRCam imaging (see Sect. 3.2). The total number of multiple images in all catalogues is 415, belonging to 150 systems (124 with spectroscopic redshifts). This makes our catalogue the largest catalogue of multiple images in a galaxy cluster field to date. The number of secure spectroscopic multiple images in our gold category transcends the number obtained with MUSE data (Be23) by 28%. The catalogues of spectroscopic lensing constraints from previous studies are summarised in Table 1. In Fig. 1, we show the HST & JWST image of the central region of MACS0416, overlaid with our lensing catalogue.

Fig. 1. RGB image of the central region of MACS0416 with multiple images. For the RGB image, we used several JWST/NIRCam and HST/WCS filters. F277W, F356W, F410M, and F444W are shown in red; F115W, F150W, and F200W are in green; and F814W, F606W, F435W, and F090W are in blue. Newly identified systems or multiply imaged clumps belonging to the gold, silver, or bronze category are shown in green, cyan, or magenta, respectively. Known systems and clumps from Di24 and Be23 are depicted in gold, silver, or bronze, with a colour corresponding to their category in this work. The bronze category also includes quartz images. Squares and diamonds represent images with new spectroscopic redshift, obtained with NIRISS and NIRSpec, respectively. Empty circles represent images for which we could not measure spectroscopic redshift directly but have new spectroscopic redshifts inferred from other images of the system. Stars represent newly identified candidate systems for which we could not obtain spectroscopic redshift. The red solid line represents the tangential (outermost) critical curve from our lensing model for redshift z = 9.

Multiple images in our catalogue are labelled as follows. New CANUCS systems are denoted with the prefix C. Previously known systems are denoted with the prefix K and the identifier from Di24, allowing for a straightforward cross-matching of the catalogues. The lowercase letter following the main ID is the label of the multiply imaged clump. We note that each clump is considered a separate system when reporting the total system count. The decimal digit represents the image ID. For instance, image K14b.3 is the third image of the second (b) clump of a galaxy with ID 14 in Di24 catalogue.

We assign systems with at least one spectroscopic redshift measurement, a system redshift z_sys, which is used for lens modelling. If MUSE redshift is available, we set z_sys = z_muse due to its higher spectral resolution compared to NIRSpec and NIRISS (except in systems K5 and K8 where we update the spectroscopic redshift; see Sect. 3.2). If z_muse is unavailable, we use z_niriss, z_nirspec or the mean of the two if both are available, due to their comparable spectral resolutions. We note that the redshift differences between different measurements are small and have a negligible impact on the lensing analysis.

Gold catalogue in Lenstool format, as well as the full catalogue with indicated gold, silver and bronze categories are available on CANUCS website¹ and will be available on Mikulski Archive for Space Telescopes (MAST)² with DOI 10.17909/ph4n-6n76 together with other CANUCS data products and lens models. Any relevant changes to the previous lensing catalogues are listed and described in the following sub-sections.

3.1. New multiple images

The search for new multiple images and multiple-image systems was performed by leveraging the strong lensing model from Be23, which was constrained with over a hundred spectroscopically confirmed multiple images, uniformly distributed across the cluster field and covering a wide redshift range 1 ≲ z ≲ 6. The model reproduced positions of their multiple images with root-mean-square distance ${\mathrm{\Delta }}_{\mathrm{r}\mathrm{m}\mathrm{s}}=0\stackrel{″}{.}43$. We used Be23 best model parameters to predict the positions of potential counter images of all sources in our photometric catalogues brighter than the 28.5 magnitude in F200W filter, using Lenstool and Eazy redshifts. The magnitude limit was applied to reduce computational and manual inspection time. It is motivated by the fact that identifying multiple-image candidates based on colour and morphology similarity becomes increasingly difficult at fainter magnitudes. Each source with predicted multiple images was then visually inspected. We examined a region around the predicted counter images in several RGB combinations of the HST and JWST filters, trying to identify sources with similar colours, morphology or photometric redshift. The visually inspected region around each source was broad, spanning over ∼15″ – much broader than Δ_rms of the model. This strategy was adopted to prevent any potential lensing model inaccuracies from biasing our dataset, as well as to account for errors in photometric redshifts. The regions around the source and image predictions were also checked for any neighbouring sources with similar configurations of multiple images, not predicted by the lens model either due to catastrophic redshift errors (e.g. with redshift below cluster redshift) or due to them not being included in our photometric catalogues. We note that a vast majority of multiple-image systems in this cluster follow a similar configuration, forming a sequence of three multiple images oriented perpendicularly to the elongated cluster morphology.

After compiling a list of candidates for multiple-image systems, we inspected their grism spectra and their fitted redshifts (see Sect. 2.2). The strategy for obtaining grism redshifts differed from the strategy used to obtain CANUCS grism redshift catalogues (Noirot et al., in prep.), where each spectrum is evaluated independently. We jointly evaluated the spectra of each multiple-image system and used all of them to assess the reliability of individual redshifts or spectral features. For some images, we could identify a potential emission line in a single orientation only (for instance, if the other orientation was too contaminated). If the redshift obtained from such emission lines matched the spectroscopic redshift of other images within redshift uncertainty, we considered the candidate spectroscopically confirmed. This approach allowed us to extract useful information from several spectra that would, on their own, prove inconclusive (see Appendix B). The newly found CANUCS systems with NIRISS spectroscopic redshifts are listed in Table 2.

After the initial search for multiple images we compiled the list of targets for spectroscopic follow-up observations with NIRSpec. The candidates were assigned priority values based on their reliability and unavailable grism redshifts. We note that our degree of confidence in some systems was high enough (e.g. due to colour and morphology similarity) that they did not necessarily require spectroscopic confirmation of all images. In such cases, a single spectroscopic redshift was enough to include the system in the gold category of multiple images with z_spec. Due to the limited possible configurations of the three MSA masks and the large number of follow-up candidates, not all of our lensing candidates were observed. On the other hand, we obtained NIRSpec spectra of some images that were initially not identified as multiple-image candidates (e.g. K100.1 at z_nirspec = 7.24 ± 0.02).

After its release, we used the multiple-image catalogue from the PEARLS collaboration (Di24), which leveraged JWST imaging data, to assess the effectiveness of our multiple-image identification method. As mentioned above, we cross-matched the catalogues and renamed our candidates to match the Di24 catalogues for a straightforward comparison. We independently found 20 multiple-image system candidates from Di24. For 11 of them, we also obtained spectroscopic redshifts, which are listed in Table A.1 with their class. Then, we visually inspected the remaining 22 system candidates from Di24, which we previously missed, and evaluated them, following the same procedure as before. We included systems K100, K99 and image K115a.3 in the gold catalogue. For systems K100 and K99 we obtained spectroscopic redshifts (see Table A.1). We also added 1 new silver system and 4 new bronze systems. For consistency, we did not include other system candidates from Di24 in our catalogue as they were too faint or otherwise could not be confirmed with the same degree of confidence as used for our candidates. Regardless, those candidates may still be interesting for future studies, so we direct the reader to Di24 for a full list. We also note that our spectroscopic redshifts of the Di24 systems are in good agreement with Di24 geometric redshift, predicted by their lens model – with relative differences below a few percent for most systems (see Fig. B.1 in Di24).

In addition, we also found 15 systems not reported in Di24. For systems C1, C2, C3 and C4, we obtained grism redshifts and included C1, C2 and C3 in the gold category and C4 in the silver category. We provide a list of new multiple-image systems with spectroscopic redshifts in Table 2. For candidates from C5 to C15, we did not obtain spectroscopic redshifts. In Table A.2, we provide their photometric redshifts and class. The class of the new multiple images was assigned upon qualitative assessment of the similarity of their colours in several HST and JWST RGB combinations (in Appendix B, we show cutouts of the new images), morphological features, their lensing configuration and photometric and spectroscopic redshifts.

In total, we identified eight gold, nine silver, and 21 bronze new multiple images. We also applied our method of finding multiple images to the other cluster of the CANUCS survey (e.g. Gledhill et al. 2024, Desprez et al. in prep.).

3.2. Known systems with MUSE redshift

Our initial catalogue of secure multiple images was compiled from Be23, Diego et al. (2024b), and Richard et al. (2021) catalogues, which were all based on HST imaging and VLT/MUSE observations, described in Sect. 1. A small correction ($\sim 0\stackrel{″}{.}2$) was applied to positions from the previous catalogues to match the positions in our JWST images. We adopted the system names from Di24. We also include additional multiply lensed clumps from Be23, which provide additional information on image deformation and are particularly useful for constraining critical curves (e.g. Bergamini et al. 2021; Grillo et al. 2016). We keep the same order of clumps as in Be23 while indexing them with letters instead of numerals. For example, system 12.6 from Be23 corresponds to the clump K1f where the name K1 indicates the galaxy with ID 1 in Di24. We visually inspected JWST data of all systems and assigned them to the gold or quartz category. The quartz systems, which were considered reliable in Di24, but are excluded from the gold category in this work, include systems K27, K48, K51, K56, K58, K59, K60, K61, K75 and multiple images K1a.3, K17.3, K25.2, K45.3, K57.3. and K77.4. Some of those sources were excluded from the gold category because they are not detected or are extremely faint (e.g. K51, K61, K75) or obscured by a foreground galaxy (K27.1). We also downgraded clump c of system K62 (corresponding to images 21.3b and 21.3c in Be23) to the quartz category due to faintness. We note that the candidates in the quartz system were not rejected by our study; we excluded them from the gold category since our data alone did not provide a sufficient degree of confidence. In some cases Be23 and Di24 select different candidates for the same image or report significantly different positions ($\gtrsim 0\stackrel{″}{.}3$). In such cases, we picked a more reliable candidate based on similar morphology or clear detection and added it to the quartz catalogue (K45.3, K58.1, K60.1, K61.3). We excluded the alternative candidates from all our catalogues (see Be23 and Di24 for a full list of candidates). The two images of system K59 at redshift z = 4.530 have different positions in the Be23 and Di24 catalogues, with K59.2 candidates separated by more than 3″. Despite clear detection of the Di24 candidate and absence of the Be23 candidate in the NIRCam data, we opted to include the latter in the quartz catalogue, as it more closely matches the MUSE detection (see Richard et al. 2021). We also excluded the K59.3 candidate from all catalogues, as it was detected as a possible counter image of Di24 candidates.

In a few spectroscopically confirmed extended systems we were not able to unambiguously identify multiply imaged features with our imaging data. An example is system K1 at redshift z = 0.94 (Hoag et al. 2016; Caminha et al. 2017), which contains a large extended arc. The arc is close to the northern brightest cluster galaxy (BCG), which increases the probability of microlensing events. Hence, it has been studied extensively due to its abundance of caustic-crossing transients (Chen et al. 2019; Kaurov et al. 2019; Yan et al. 2023). The counterimage of the arc (K1a.3) is spectroscopically confirmed (Caminha et al. 2017) and is lensed by the cluster member with ID 8971 (ID from Be23; see Sect. 4). However, as it is difficult to securely identify exact counterparts of individual clumps in the arc with HST and JWST images, we decided to include the counterimage in the bronze category using the position from Di24 (the image was excluded on the same grounds from Bergamini et al. 2021 and Be23 lens models).

We find a similar situation in system K8, which comprises a very red arc and a spiral galaxy at redshift z = 2.04. The arc has also been investigated by Di24 as it hosted two extremely magnified star candidates. The system has a strong H_α emission line in the F200W NIRISS grism channel, enabling a reliable spectroscopic redshift measurement with NIRISS in both the arc and counterimage. A part of the arc has also been observed with NIRSpec, confirming the spectroscopic redshift of z_nirspec = 2.046 ± 0.001 (see Table A.3 and Figs. B.3 and C.2). The updated redshift differs from the previously measured MUSE redshift (z_muse = 1.953; see Table B.4 in Richard et al. 2021) by 0.09, which is above our redshift uncertainty. Furthermore, we identified the counterpart of the brightest clump of the arc in the counterimage. To this end, we examined several RGB combinations of JWST and HST filters, two of which are shown in Fig. 2. Additionally, we produced the map of H_α emission line flux from F200W grism spectrum (spectrum is shown in Fig. B.3 and the emission line map in the right panel of Fig. 2). The brightest clumps a.2 and a.3 in the arc show clear H_α emission. The emission is expected to be higher than its counterpart a.1 in the spiral image due to higher magnifications in the arc. This rules out two of the brightest H_α emitting clumps above and below the galactic centre. Based on colour similarity and faint H_α emission, we identify the clump on the bottom left relative to the galactic centre to be the most likely candidate for clump a.1. We verified that given the shape of the caustic and the positions of bright clumps in the source plane, the selected a.1 clump is the farthest below the caustic line; in other words, it is one of the most likely to be multiply imaged. Unlike K8a.1, we could not securely identify fainter counter images K8b.1 and K8c.1 to clumps K8b and K8c in the arc.

Fig. 2. Multiply imaged clumps a, b, and c of system K8. The lower panels show the extended arc, with two counter images of the spiral galaxy shown in the upper panels. The two panels on the left show two RGB combinations of NIRCam images. In the right panel, we show Hα emission line flux obtained from the F200W grism spectra (see Appendix B), which was used to identify image K8a.1.

We did a minor correction to the spectroscopic redshift of system K5. Each of the three images of system K5 shows clear [O III] λλ4959, 5007 (blended at the grism resolution), [O II] λ3727 and H_α emission lines in the grism spectrum (see Fig. B.3), yielding spectroscopic redshift of z_niriss = 1.84 ± 0.01, which differs from z_muse = 1.8178 (Richard et al. 2021) by more than our redshift uncertainty estimate.

We also spectroscopically confirmed image K62b.3, which was reported in Di24 (62.c) as a candidate. We measured z_nirspec = 5.14 ± 0.01 and included it in the gold category with system redshift set to z_muse = 5.106. In system K26 we identified an alternative candidate to K26.3 from Di24. The grism spectrum (see Fig. B.3) of our candidate shows the [O III] emission, yielding z_niriss = 2.92 ± 0.02. This agrees with z_muse = 2.9259; hence, we included our candidate in the gold catalogue. All relevant updates to known systems with MUSE redshifts are summarised in Table A.3. We also modified some of the known systems without spectroscopic redshift from Di24. We provide our alternative candidate K81.3 to image 81.c from Di24 as well as a new counter image K96.3 of system K96. Both image candidates were found using the method described in Sect. 3.1 and are listed in Table A.3.

In some lensed sources with resolved structure in our imaging data, we included different features as separate multiple-image systems. This includes several features in highly magnified arcs, included in Be23 as well as additional new families K8b, K8c, K22b, K49b, K55b, K64b, K85b and K115b. Out of those, K8b, K8c and K115b are situated in strongly lensed arcs and are useful for locally constraining the critical curves. Others belong to well-separated images and provide additional information on their orientation and local deformations.

K33 and K77 (14.1 and 14.2 in Be23) are a unique lensing system consisting of two compact low luminosity star-forming stellar complexes at redshift 3.221 (Vanzella et al. 2017). The sources are first multiply lensed by the global cluster profile and then further split into several images by two cluster members, resulting in a system of seven multiple-image candidates (shown in Fig. 3). The systems are included as in Di24 with a few modifications. We excluded K77.4 (as given in Di24) since it is indistinguishable from K33.4. We added our K77.4 candidate, which belongs to the same multiple image as K33.5. However, due to its uncertain position, we include it in the quartz category. We included system candidate 110 from Di24 to system K33 as images K33.6 and K33.7. Of the two, we include K33.6 to our gold catalogue, following Be23. For image K33.7, which has not been included as a reliable candidate in previous studies, we provide a tentative spectroscopic confirmation based on the [O III] emission line at around 2.1 μm in the grism spectra, shown in Fig. 3 (see also the NIRSpec spectrum of K33.1 in Fig. C.2). The emission can also be seen in other images of the system. This makes system K33 a very likely system with 7 multiple images. In this work, we add K33.7 to the silver category and do not use it to constrain the model.

Fig. 3. Upper panel: RGB cutouts (with the same filter combinations as used in Fig. 1) of cluster members 8971 and 5068 with surrounding multiple-image systems K33, K77, and K111. The inset cutout shows the counterimage K33.4. Gold, silver or bronze colours indicate the category (reliability) of multiple images. Dotted and dashed lines with arrows represent the orientation of 150C and 150R grism spectrum, respectively. They are also indicated on the bottom two panels. The dotted line indicating the 150C spectrum orientation extends to galaxy 8971, positioned 6″ to the northeast of galaxy 5068. Second panel: RGB cutouts with removed cluster galaxies, enabling a better colour comparison. Third panel: F200W 150C grism spectrum of galaxy 5068 showing potential [O III] emission lines observed at 2.1 μm from K33.6 and K33.7, and emission from K33.2, K77.1, and K33.5 on the left. The wavelength calibration matches the spectrum of galaxy 5068. Fourth panel: F200W 150R grism spectrum of galaxy 5068, indicating the [O III] emission of K33.6 and K33.7.

In Table 3 we summarise the differences between our and Di24 catalogues. We indicate sources that are missing in either catalogue and sources that were upgraded to or downgraded from the most reliable categories.

4. Lens modelling

The strong lensing model was optimised using the parametric strong lens modelling software Lenstool³ (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009), which uses Markov chain Monte Carlo sampler BayeSys (Skilling 2004) to sample the posterior probability of the model parameters ξ. The posterior probability is computed as

$$\begin{array}{c}\hfill P(\xi |r^o)=P(r^o|\xi )P(\xi ),\end{array}$$(1)

where the prior probability P(ξ) is the product of either uniform or Gaussian probabilities for each parameter. The likelihood P(r^o|ξ) given the positions of N multiple images r^o, each with uncertainty σ_p is defined as

$$\begin{array}{c}\hfill P(r^o|\xi )={\left(2\pi {\sigma }_p^2\right)}^{-\frac{N}{2}}exp\left(-\frac{{\chi }^2(r^o,\xi )}{2}\right).\end{array}$$(2)

The χ² is calculated from the distance between the observed positions r^o and model predicted positions r^m:

$$\begin{array}{c}\hfill {\chi }^2(r^o,\xi )={\displaystyle \sum _{j=1}^{n_{\mathrm{sys}}}}{\displaystyle \sum _{i=1}^{n_{\mathrm{im}}^j}}\frac{||r_{i,j}^o-r_{i,j}^m(\xi ){||}^2}{{\sigma }_p^2},\end{array}$$(3)

where the sum goes over n_sys multiple-image systems, each containing n_im^j multiple images. For optimisation, Lenstool uses 10 interlinked chains to avoid local minima and varies the likelihood influence during the burn-in phase to control the convergence speed (see Jullo et al. 2007).

To evaluate the accuracy of the lens model, we used the root-mean-square distance from the observed positions of multiple images and positions predicted by the best model with parameters ξ_best:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{rms}}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N}||r_i^o-r_i^m({\xi }_{\mathrm{best}}){||}^2}.\end{array}$$(4)

In this work, we use the mass distribution parameterisation from Be23, which was introduced as the best-performing lens model (LM-4HALOS) in Bergamini et al. (2021). It describes the total (baryonic and dark) matter distribution as a superposition of several elliptical halos with dual pseudo isothermal elliptical (dPIE) density profile (Limousin et al. 2005; Elíasdóttir et al. 2007). Each dPIE halo is parametrised with 8 parameters: coordinates x and y of the centre, redshift z (set to 0.397 for all halos), ellipticity e, orientation θ, core radius r_core, cut radius r_cut and the Lenstool fiducial velocity dispersion σ_lt, which is related to central velocity dispersion σ₀ of the dPIE profile as ${\sigma }_{\mathrm{lt}}=\sqrt{2/3}{\sigma }_0$. Positions x and y in this work are given in arcsec relative to the reference position α = 64.038142°, δ = −24.067472°, which coincides with northern BCG. The total mass of each dPIE halo is computed as

$$\begin{array}{c}\hfill M_{\text{tot}}=\frac{\pi {\sigma }_0^2{r}_{\text{cut}}}{G}·\end{array}$$(5)

The global mass component in our lensing model features 4 large-scale cluster halos. Two are centred on the northeastern and southwestern BCGs with 30″ and 15″ wide uniform position priors, respectively (H1 and H2 in Fig. 4). A second halo is added to the southern region with a wider position prior as a second-order correction to the southern mass component (H4). Following previous studies (Caminha et al. 2017; Bergamini et al. 2019, 2021; Be23), another halo (H3) was added to the northeast of the northern BCG with fewer free parameters (the halo is kept circular without free ellipticity and orientation), which serves as a correction to the northeastern cluster component and accounts for a small overdensity of cluster galaxies. All four cluster components have free σ_lt and r_core, while we fix their r_cut to a large value (2000″), as it cannot be constrained by the multiple images in the inner cluster regions.

Fig. 4. Relative difference between the convergence (κ) from the updated strong lensing model with JWST data and the lens model from Be23, which uses the same parameterisation. The positions of four main cluster halos (H1–H4) and the two galaxies not included in the scaling relations (galaxy 8971 and the foreground galaxy) are marked with red pluses. Positions of the northern and southern BCGs are indicated with blue circles. Green crosses represent new multiple images and old multiple images with new spectroscopic redshifts, not included in the Be23 nor the Di24 catalogues. Black solid and dashed lines represent the tangential (outermost) critical curve for sources at redshift 9, derived from the new and Be23 lensing models, respectively.

All prior ranges of the free parameters are listed in Table 4 and match the values from Be23. Only one modification was required in the position prior of the H3 cluster halo, which had to be narrowed down to −39″ < x < −25″ and 2″ < y < 16″. This modification proved necessary as the model with extended prior range (−55″ < x < −25″, 0″ < y < 30″; see Be23) had problems with convergence after the addition of new spectroscopic systems – the model placed the 3^rd cluster halo in a local minimum of χ² outside of the region covered with multiple images. After narrowing down the prior range, we found that the new solution, consistent with the best-fit position from Be23, yielded significantly better χ² and Δ_rms. We also note that the chosen prior range is nonetheless much larger than the 68% posterior confidence interval (see Table 5). This also further demonstrated the need for the additional mass component H3 in the overdensity of cluster galaxies to the northeast of the main cluster component (H1).

The model also includes four dPIE clumps describing the cluster gas distribution, derived from X-ray Chandra observations (Bonamigo et al. 2018, 2017). They are included as fixed clumps without free parameters, which is justified by a smaller set of assumptions when deriving the gas component mass distribution and smaller statistical errors compared to other components (Bonamigo et al. 2017). The parameters of the four clumps describing the cluster gas can be found in Table 1 of Be23.

Smaller sub-halos, which are associated with cluster members, were modelled as circular dPIE profiles with a vanishing core radius r_core. The cut radius r_cut,i and velocity dispersion σ_lt,i of each cluster member i scale with luminosity L_i of each galaxy were calculated according to the following relations (Jorgensen et al. 1996; Natarajan & Kneib 1997):

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{lt},i}& ={\sigma }_{\mathrm{lt}}^{\mathrm{ref}}{\left(\frac{L_i}{L_{\mathrm{ref}}}\right)}^{\alpha },\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill r_{\mathrm{cut},i}& =r_{\mathrm{cut}}^{\mathrm{ref}}{\left(\frac{L_i}{L_{\mathrm{ref}}}\right)}^{{\beta }_{\mathrm{cut}}}.\hfill \end{array}$$(7)

In this work, we use the cluster member catalogue with magnitudes from Bergamini et al. (2021), comprising 213 galaxies, 171 of which have spectroscopic confirmation. The positions of galaxies were first shifted to match the astrometry of CANUCS images (the correction was small, $\sim 0\stackrel{″}{.}2$, and was also applied to the other mass components). The slopes of the scaling relations were derived in Bergamini et al. (2021), following the procedure described in Bergamini et al. (2019). They leverage measured velocity dispersions of cluster members with MUSE spectroscopy to obtain α = 0.3, and derive β_cut = 0.6 by assuming that the sub-halo mass M_tot,i scales with luminosity as $M_{\mathrm{tot},\mathit{i}}\propto L_i^{0.2}$. As in Be23 and Bergamini et al. (2021), we assume a uniform prior between 1″ and 50″ for $r_{\mathrm{cut}}^{\mathrm{ref}}$, and a Gaussian prior for ${\sigma }_{\mathrm{lt}}^{\mathrm{ref}}$ centred at 248 km s⁻¹ with standard deviation of 28 km s⁻¹. The Gaussian prior was determined with MUSE spectroscopy as described in Bergamini et al. (2019) and helps break the degeneracy between ${\sigma }_{\mathrm{lt}}^{\mathrm{ref}}$ and $r_{\mathrm{cut}}^{\mathrm{ref}}$, allowing for more accurate characterisation of the sizes of the cluster members. Luminosity L_i of the cluster members was estimated with HST F160W magnitudes from Bergamini et al. (2021), and for L_ref, the reference magnitude of 17.02 was chosen. As a test, we optimised a lens model of the upper half of the cluster (to reduce the dataset size and computational time). We used measured ellipticities instead of circular sub-halos, as well as using NIRCam F090W magnitudes in place of F160W magnitudes. Neither case resulted in a χ² improvement compared to an identical run with the old cluster member catalogue; hence, we did not implement these changes in the final model. We leave the refinement of the cluster member catalogue with new JWST data for a future study, for instance, with new NIRISS spectroscopic catalogues of the full cluster field (Noirot et al., in prep.).

Two galaxy-scale halos were not modelled with scaling relations (Eqs. (6) and (7)). A foreground galaxy (z = 0.112), which perturbs the lensing potential south of the H4 cluster halo, was modelled as a circular dPIE halo at the cluster redshift with free σ_lt and r_cut and vanishing r_core. Cluster member 8971, which produces the galaxy-galaxy lensing of systems K77 and K33 (left panel of Fig. 3, discussed in Sect. 3.2), was included as an elliptical dPIE halo with free ellipticity e, orientation θ, σ_lt and r_cut with vanishing r_core. See Table 4 for the prior values.

The parameter optimisation was done in two steps. For the first run of the model, we used the position uncertainty ${\sigma }_p=0\stackrel{″}{.}3$. We note that this is larger than the actual position uncertainty of the multiple images from the HST and JWST data. However, the actual accuracy (∼40 mas) is beyond the current lens modelling capabilities, which do not consider contributions such as from the line-of-sight mass distribution or the scatter around the mass-to-light relation of the cluster galaxies (e.g. Jullo et al. 2010). The limitations of the lens model were further reflected in high χ² values after the first run. After obtaining the minimum of the χ², we rescaled the position uncertainties σ_p by a constant factor so that the minimal χ² was approximately equal to the number of degrees of freedom (D.o.F):

$$\begin{array}{c}\hfill \mathrm{D}.\mathrm{o}.\mathrm{F}=2N-2n_{\mathrm{sys}}-N_{\mathrm{par}}.\end{array}$$(8)

The term N_par in Eq. (8) represents the number of free parameters of the mass distribution (30 in this work), and 2n_sys is equal to the number of source positions of multiple-image systems, which represent an additional 222 free parameters. With 2N = 606 measured coordinate positions of multiple images in our gold catalogue, we obtained D.o.F = 354. After rescaling σ_p so that χ²/D.o.F ≈ 1, we reran the model. The position uncertainty used in the final optimisation was 0$\stackrel{″}{.}$49. This rescaling is standard practice in strong lensing studies (e.g. Be23, Grillo et al. 2016), and it serves as an attempt to account for systematic inaccuracies of the lens model. We emphasise that rescaling the uncertainty does not influence the best-fit model; rather, it provides better-approximated uncertainties of the model parameters.

After the 83 000 burn-in steps, during which the χ² value reached a stable value showing no significant average variation over large number of steps, we sampled the posterior with 40 000 samples. The optimisation required over two weeks of computational time on 130 CPU cores. The best-fit model resulted in ${\mathrm{\Delta }}_{\mathrm{r}\mathrm{m}\mathrm{s}}=0\stackrel{″}{.}52$. We do not use other fit quality criteria such as the BIC or AIC (e.g. see Bergamini et al. 2021) as the comparison of different model parameterisations is beyond the scope of this work. We also note the dependence of those criteria on the adopted σ_p, which does not reflect the uncertainty of the measurements in our case. In Table 5, we provide the median values and the 68% confidence intervals of the lensing model parameters. We note that the median parameter values may differ from the best model values, which we provide separately as a Lenstool parameter file on MAST. However, it enables a more straightforward comparison with the analogous table in Be23.

We also note that similar to the Bergamini et al. (2021) model, our model also predicts the multiple-image K33.7, for which we provided a tentative spectroscopic confirmation in Sect. 3.2 and included it in the silver category. K33.7 is lensed by a cluster member with ID 5068, which is modelled with scaling relations (Eqs. (6) and (7)).

5. Discussion

5.1. Comparisson with Bergamini et al. (2023)

In Fig. 4 we show the positions of the main mass components as well the difference of convergence (κ) maps and the critical curves at redshift z = 9 between our and Be23 best-fit models. Comparing our mass distribution with Be23, we can see that the relative κ difference between the models is less than a few percent over most of the inner cluster regions. The median difference |κ_new − κ_old|/κ_new in the regions displayed in Fig. 4 is only 1.6%. We note that the close agreement between our density maps is, to a large degree, the consequence of using the same lens modelling method and parameterisation as Be23. A recent comparison of different MACS0416 lens models (including our, Be23 and Di24 models) by Perera et al. (2024) shows that the median relative density differences can reach ∼20% if different modelling methods are used. The overall agreement between our and Be23 mass distributions is also reflected in the positions of the critical curves. As the critical curves trace the regions of extreme magnification, their accurate description is vital for studying high redshift objects. In Fig. 4, we show the tangential critical curve for sources at z = 9 to illustrate that the high redshift critical curves derived from our model align well with those from Be23.

The differences between our and Be23 density maps and critical curves are more drastic in the southern part of the cluster. The foreground galaxy in our best-fit model is significantly more massive than in Be23 (Fig. 4), which follows from our higher best-fit velocity dispersion σ_lt = 136 km s⁻¹ (Table 5) compared to only 77 km s⁻¹ in Be23 best-fit model (see the σ_lt dependence of the dPIE mass in Eq. (5)). However, we note that our median σ_lt of the foreground galaxy is significantly lower with large relative uncertainties (${\sigma }_{\mathrm{lt}}={99}_{-32}^{+27}{\text{km}}^{}$). The poorly constrained mass distribution around the foreground galaxy is unsurprising as there are very few constraints in the southern-most cluster regions (see Fig. 1).

Fig. 5. Top: Projected mass Mr in the circular aperture with radius r centred at the northern (red) or the southern BCG (yellow) computed from the best-fit model (solid line) with indicated 68% confidence intervals. The profiles from the best-fit Be23 models are shown in green and blue. The dashed lines represent the cumulative mass Mr contained in the cluster member halos (modelled with the scaling relations). Small vertical bars represent the distances of multiple images from the nothern or southern BCG. Middle: Difference between the total cumulative mass Mr from our model and the Be23 best-fit model relative to our Mr. The confidence intervals are computed from our Mr uncertainty. Bottom: Relative difference between the cumulative mass contained in the cluster member halos from our model and that of Be23.

In Fig. 5 we show the projected total mass M_r enclosed in the circular aperture with radius r (in kpc) from the northern or southern BCG as well as the mass contained in the cluster member halos. The total mass profiles from our model show an excellent agreement with the Be23 best-fit model, with a relative difference of less than a few percent in the radial range containing multiple images. This is consistent with the literature showing that the circularly averaged total mass profiles are well constrained by lens models – even if different modelling techniques are employed (Raney et al. 2020; Meneghetti et al. 2017). However, we find a larger difference when comparing the mass contained in the cluster member halos. Our best-fit scaling relation normalisation values ${\sigma }_{\mathrm{lt}}^{\mathrm{r}\mathrm{e}\mathrm{f}}=206{\text{km}}^{}$ and $r_{\mathrm{cut}}^{\mathrm{r}\mathrm{e}\mathrm{f}}=16.0″$ (Table 5) result in cluster members with 36% higher masses but more extended profiles with lower central densities than in Be23 (${\sigma }_{\mathrm{lt}}^{\mathrm{r}\mathrm{e}\mathrm{f}}=217{\text{km}}^{}$, $r_{\mathrm{cut}}^{\mathrm{r}\mathrm{e}\mathrm{f}}=10.5″$ ). This can be seen in Fig. 4 and in the last panel of Fig. 5 where the higher total masses of our galaxies and lower mass enclosed in the inner 30 kpc around the BCGs result in a rise of the cumulative mass difference with r. Significant cluster member mass differences despite the matching total mass profiles indicate a known degeneracy between the cluster galaxies and the smooth mass component in parametric lens modelling (e.g. Limousin et al. 2016). Whereas the total mass profiles already seem to be well constrained with 88 multiple-image systems in Be23, the cluster member component shows larger variations upon the addition of new systems. We find this despite the Gaussian prior on ${\sigma }_{\mathrm{lt}}^{\mathrm{r}\mathrm{e}\mathrm{f}}$, which is added to reduce the sub-halo mass uncertainty (Bergamini et al. 2019).

We also compared the magnifications of multiple images computed with Di24, Be23 and our lens model. For a meaningful comparison of the lens models, we used the positions and redshifts of spectroscopic multiple images given in Table B from Di24 (while accounting for a slight shift in astrometry between different works). We found that the rms relative magnification difference 2(μ₁ − μ₂)/(μ₁ + μ₂) between our and Be23 model is 10% for images with μ < 10 (in our model) and 35% for μ > 10. The larger relative differences at higher μ are a consequence of rapid magnification variations near the critical curves, which makes those images susceptible to small differences between lens models. Despite different dataset sizes, the magnification differences between our two lens models (both constrained with Lenstool) are relatively small. When comparing our magnifications with those provided by Di24, which used a different modelling code (WSLAP+; Diego et al. 2005), the relative differences are larger, reaching 41% for images with μ < 10 and 85% for μ > 10. For 80% of the multiple images, Di24 magnifications are higher than our best model μ. Large μ differences when using different modelling methods (which transcend the statistical uncertainties given by individual models) are a known issue that has been investigated with the Hubble Frontier Lens models (Priewe et al. 2017; Raney et al. 2020). As our models differ from Di24 not only in the methods but also in the set of constraints, we cannot meaningfully investigate the impact of the modelling techniques on magnifications in this work.

When comparing Δ_rms we find that our ${\mathrm{\Delta }}_{\mathrm{r}\mathrm{m}\mathrm{s}}=0\stackrel{″}{.}52$ is larger than ${\mathrm{\Delta }}_{\mathrm{r}\mathrm{m}\mathrm{s}}=0\stackrel{″}{.}43$ in Be23 by ∼20%. The increase is not surprising, considering that the newly derived mass distribution accommodating the extended dataset (with > 20% more multiple-image systems) differs from the one that minimised the χ² (and Δ_rms) of the old dataset using the same model parametrisation.

5.2. The effectiveness of NIRISS for strong lensing studies

In Fig. 6, we show the redshift and F200W magnitude distribution of multiple images with spectroscopic system redshifts. Different markers and colours indicate whether system redshifts of multiple images could be measured with JWST (NIRISS or NIRSpec), MUSE or both. We also indicate the redshift ranges where the most significant emission lines fall in F115W, F150W or F200W grism spectrum. The redshift distribution of multiple images is shown in Fig. 7. Most new multiple images with NIRISS redshifts are found around redshifts 2 and 3 where at least two of the [O II], [O III] and H_α emission lines fall in the most sensitive range of the grism spectrum. In that range, NIRISS also successfully recovers most redshifts measured by MUSE.

Fig. 6. Redshift and F200W magnitude of multiple images in our gold category with known spectroscopic redshifts. Grey and black points represent known images with MUSE redshifts. For black squares, we also obtained NIRISS redshifts. Gold squares represent newly identified images for which we could obtain NIRISS redshift. Blue diamonds represent new images without NIRISS but with NIRSpec redshift. Images for which the redshift could not be reliably measured but with measured system redshift from other images with NIRISS or NIRSpec are marked with gold and blue circles, respectively. Each multiple image is considered once (i.e. different clumps are not considered separately). The coloured bands at the bottom of the figure represent redshift ranges where some of the most prominent emission lines fall in the wavelength range, covered by each of the three NIRISS grism filters, used in this work. The darker shade of the grism range represents a more sensitive region (see Sect. 2.2).

Fig. 7. Redshift distribution of multiple images included in our gold catalogue with known spectroscopic redshift. Images for which the spectroscopic redshift was obtained only by MUSE are shown in grey. Images for which we provided spectroscopic redshifts with JWST are shown in orange. Striped regions represent images with both JWST and MUSE redshift measurements. In the histogram, each multiple image is counted once (i.e. different multiply lensed clumps are not considered separately). The red solid line represents the expected redshift distribution of multiple-image systems, considering the redshift distribution of sources in the field and multiply imaged area in the source plane, which increases with redshift.

Furthermore, we roughly estimate the expected redshift distribution of multiple-image systems (red line in Fig. 7). This was done by first calculating the median distribution of photometric redshifts in the NIRCam flanking fields of all 5 CANUCS clusters. We considered all galaxies with reliable photometric redshifts regardless of their magnitudes. The redshift distribution of the sources peaks between redshift 1 and 2 due to the maximum angular diameter distance (and thus the largest volume covered by the same field-of-view). This unlensed distribution was then multiplied by the area of the multiply imaged region in the source plane at each redshift. The area is equal to the area bounded by the tangential caustics, which could be estimated with magnification and deflection maps produced by Lenstool. The area increases rapidly with redshift. The interplay between increasing lensed area, changing angular diameter distance and the decline in the number of detected sources with redshift results in a distribution that peaks at the redshift at which we expect to find the largest number of strongly lensed systems.

In Fig. 7, we show that in MACS0416, we expect to find most multiply lensed sources around redshift 2. At this redshift, [O II], [O III] and H_α emission lines all fall in the range of our F115W, F150W and F200W grism spectra, making NIRISS spectroscopy ideal for finding multiple images up to faint magnitudes (even fainter than magnitude 27 in F200W filter; see Fig. 6). MUSE spectral range is, on the other hand, less ideal for this redshift range, as it cannot detect any of the above-mentioned emission lines above z ∼ 1.5 and Lyman α emission below z ∼ 3. Given sufficiently deep observation, such as in this cluster, MUSE is still effective at finding sources based on fainter emission lines (e.g. C III] at 1907 Å; see the emission line catalogues from Richard et al. 2021). However, it is evident that NIRISS grism spectroscopy, such as that used in this work, complements even the deepest MUSE observations.

NIRISS grism spectroscopy also provides emission line maps of bright emission lines (e.g. H_α) with high spatial resolution, which can, in combination with photometry, help identify multiply imaged clumps in extended multiple images and giant arcs. High magnifications in the giant arcs make them especially important for probing faint small-scale systems at high redshift. On the other hand, high magnification gradients near the critical curves make magnifications highly uncertain and any additional constraints invaluable. We demonstrated the use of emission line maps with system K8 (Fig. 2), where H_α flux complemented colour information, obtained with broadband NIRCam imaging. It helped us rule out some of the clumps with high H_α flux in the spiral galaxy as potential counter images of the brightest clump in the giant arc. We also point out that a significant fraction of systems, including some of the giant arcs (see Di24), are found in the redshift range where at least two of the [O II], [O III] and H_α emission lines fall in our grism range (Fig. 6). Thus, this method could be extended from using individual emission line maps to using quantities such as emission line ratios, which are independent of magnification (provided they originate from the same source area) and could become a useful tool in future lensing studies.

In addition to NIRISS spectroscopy, this work also demonstrated the value of NIRSpec MSA prism spectroscopy for strong lens modelling. Its wide wavelength coverage enables redshift measurements from multiple emission lines in a wide redshift range (e.g. see the spectra in Fig. C.1). For instance, it allowed us to measure the redshift of system K100, z_nirspec = 7.24 ± 0.02, which became the highest spectroscopic redshift of a multiple-image system in this cluster. However, unlike NIRISS and MUSE observations, NIRSpec MSA observations target preselected objects, which requires a reliable multiple-image candidate identification based on prior imaging data and photometric redshifts. Furthermore, if other follow-up targets are prioritised based on various science goals, not all multiple images can be observed due to the limited possible configurations of MSA masks. In this work, we obtained only one NIRSpec redshift per multiple-image system (see Table A.1). Although this is enough for determining the system redshift z_sys, it is, on its own, insufficient to confirm each multiply imaged candidate. To assess the reliability of our gold multiple-image systems, it was thus crucial to use the NIRSpec data in combination with our imaging data and NIRISS spectroscopy.

6. Conclusion

In this work, we have presented an updated strong lensing model of galaxy cluster MACS J0416.1-2403 (MACS0416) that leverages the new JWST imaging and spectroscopic data from CANUCS. We expanded the existing catalogues of multiple images by searching for new multiple images in our NIRCam imaging data. To this end, we first predicted multiple images from sources in CANUCS photometric catalogues and then visually examined the region around each prediction in several HST and JWST/NIRCam RGB combinations. Our search independently discovered most of the already known candidates while also finding 15 previously unknown systems. The catalogue was further augmented by our spectroscopic redshift measurements. We examined the spectra obtained with NIRISS wide-field slitless spectroscopy of all our candidates and observed some of them with the NIRSpec multi-object prism spectroscopy. We measured spectroscopic redshift for 17 multiple-image systems with previously unknown spectroscopic redshifts. We furthermore corrected the redshift of two known systems and confirmed two new counter images. We also re-evaluated known sources with MUSE redshift (e.g. Be23, Richard et al. 2021, Bergamini et al. 2021, Caminha et al. 2017) and candidates found by previous studies that used NIRCam imaging data (Di24). Our complete catalogue of all candidates comprises 415 multiple images from 150 systems, 124 of which are known to have spectroscopic redshift. Of these, we selected 303 reliable multiple images from 111 systems with spectroscopic redshift that we used for constraining the lens model (i.e. the gold catalogue). This makes our catalogue the largest catalogue of spectroscopic multiple images in a galaxy cluster field to date.

Our lens model was constrained using the parametric lens modelling code Lenstool with the parameterisation presented in Bergamini et al. (2021) and Be23. The new model is able to reproduce the multiple-image position with ${\mathrm{\Delta }}_{\mathrm{r}\mathrm{m}\mathrm{s}}=0\stackrel{″}{.}52$ and is in good agreement with the Be23 lens model, which uses the same Lenstool parameterisation and pre-JWST data. The differences between our and the Be23 lens model are larger in poorly constrained southern regions of the cluster. The median relative κ difference between our models is 1.6%, and the radial profiles of the total mass in a circular aperture agree within a few percent. We find larger disagreements between our and the Be23 lens model in the mass profiles of the galaxy sub-halo mass component, resulting from different scaling relation normalisation parameter values. This is indicative of a known degeneracy between cluster members and the smooth mass components despite well-constrained total mass distribution (e.g. Limousin et al. 2016). We also compared magnifications of multiple images and found that the relative μ differences given by our model differ from those obtained with the Be23 model by 10% for images with our μ < 10 (and 30 % if μ > 10). The differences are larger when compared with Di24, which uses a different lens modelling method (as well as a different dataset), amounting to 41% for μ < 10 (85% for μ > 10). Large μ differences between different modelling techniques are consistent with the findings from pre-JWST data (Priewe et al. 2017; Raney et al. 2020).

In this work, we also demonstrated the utility of NIRISS wide-field slitless spectroscopy in F115W, F150W, and F200W filters for strong lensing studies. We show that such spectroscopy is the most effective at measuring redshifts between redshifts 2 and 3, where at least two of the three significant emission lines (H_α, [O III] λλ4959, 5007 and [O II] λ3727) can be detected. With a simple estimate, considering the redshift distribution of galaxies in the field and the multiply lensed area as a function of redshift, we showed that this redshift range is particularly useful for strong lensing studies, as it contains the most multiple images. With system K8, we further demonstrated how spatially resolved NIRISS spectroscopy and emission line maps can aid in understanding of the extended multiple images and giant arcs. Since MACS0416 has been observed by some of the deepest HST imaging and MUSE spectroscopy programs, our study thus serves as a comparison between the results attainable with JWST and results currently attainable with the best ground- or space-based facilities using the largest known dataset in a single cluster.

Data availability

Appendices A, B, and C with tables summarizing our updates to the multiple-image catalogue, the spectra, and the RGB cutouts are available at https://zenodo.org/records/14850937. We provide our catalogue of multiple images used to constrain the model with their magnifications and the full catalogue of multiple-image candidates on CANUCS website (https://niriss.github.io/) together with the lens model presented in this work and a tool for computing magnifications of arbitrary sources. Following the first CANUCS data release, the lens model will also be available on STScI Mikulski Archive for Space Telescopes (MAST) with DOI 10.17909/ph4n-6n76.

Acknowledgments

We thank the anonymous referee whose comments greatly enhanced the quality of this paper. GR, MB, AH, VM and RT acknowledge support from the ERC Grant FIRSTLIGHT and from the Slovenian national research agency ARRS through grants N1-0238, P1-0188. MB acknowledges support from the program HST-GO-16667, provided through a grant from the STScI under NASA contract NAS5-26555. This research was enabled by grant 18JWST-GTO1 from the Canadian Space Agency and funding from the Natural Sciences and Engineering Research Council of Canada. This research used the Canadian Advanced Network For Astronomy Research (CANFAR) operated in partnership by the Canadian Astronomy Data Centre and The Digital Research Alliance of Canada with support from the National Research Council of Canada the Canadian Space Agency, CANARIE and the Canadian Foundation for Innovation.

References

All Tables

All Figures

Fig. 1. RGB image of the central region of MACS0416 with multiple images. For the RGB image, we used several JWST/NIRCam and HST/WCS filters. F277W, F356W, F410M, and F444W are shown in red; F115W, F150W, and F200W are in green; and F814W, F606W, F435W, and F090W are in blue. Newly identified systems or multiply imaged clumps belonging to the gold, silver, or bronze category are shown in green, cyan, or magenta, respectively. Known systems and clumps from Di24 and Be23 are depicted in gold, silver, or bronze, with a colour corresponding to their category in this work. The bronze category also includes quartz images. Squares and diamonds represent images with new spectroscopic redshift, obtained with NIRISS and NIRSpec, respectively. Empty circles represent images for which we could not measure spectroscopic redshift directly but have new spectroscopic redshifts inferred from other images of the system. Stars represent newly identified candidate systems for which we could not obtain spectroscopic redshift. The red solid line represents the tangential (outermost) critical curve from our lensing model for redshift z = 9.

In the text

	Fig. 2. Multiply imaged clumps a, b, and c of system K8. The lower panels show the extended arc, with two counter images of the spiral galaxy shown in the upper panels. The two panels on the left show two RGB combinations of NIRCam images. In the right panel, we show Hα emission line flux obtained from the F200W grism spectra (see Appendix B), which was used to identify image K8a.1.
In the text

Fig. 3. Upper panel: RGB cutouts (with the same filter combinations as used in Fig. 1) of cluster members 8971 and 5068 with surrounding multiple-image systems K33, K77, and K111. The inset cutout shows the counterimage K33.4. Gold, silver or bronze colours indicate the category (reliability) of multiple images. Dotted and dashed lines with arrows represent the orientation of 150C and 150R grism spectrum, respectively. They are also indicated on the bottom two panels. The dotted line indicating the 150C spectrum orientation extends to galaxy 8971, positioned 6″ to the northeast of galaxy 5068. Second panel: RGB cutouts with removed cluster galaxies, enabling a better colour comparison. Third panel: F200W 150C grism spectrum of galaxy 5068 showing potential [O III] emission lines observed at 2.1 μm from K33.6 and K33.7, and emission from K33.2, K77.1, and K33.5 on the left. The wavelength calibration matches the spectrum of galaxy 5068. Fourth panel: F200W 150R grism spectrum of galaxy 5068, indicating the [O III] emission of K33.6 and K33.7.

In the text

Fig. 4. Relative difference between the convergence (κ) from the updated strong lensing model with JWST data and the lens model from Be23, which uses the same parameterisation. The positions of four main cluster halos (H1–H4) and the two galaxies not included in the scaling relations (galaxy 8971 and the foreground galaxy) are marked with red pluses. Positions of the northern and southern BCGs are indicated with blue circles. Green crosses represent new multiple images and old multiple images with new spectroscopic redshifts, not included in the Be23 nor the Di24 catalogues. Black solid and dashed lines represent the tangential (outermost) critical curve for sources at redshift 9, derived from the new and Be23 lensing models, respectively.

In the text

Fig. 5. Top: Projected mass Mr in the circular aperture with radius r centred at the northern (red) or the southern BCG (yellow) computed from the best-fit model (solid line) with indicated 68% confidence intervals. The profiles from the best-fit Be23 models are shown in green and blue. The dashed lines represent the cumulative mass Mr contained in the cluster member halos (modelled with the scaling relations). Small vertical bars represent the distances of multiple images from the nothern or southern BCG. Middle: Difference between the total cumulative mass Mr from our model and the Be23 best-fit model relative to our Mr. The confidence intervals are computed from our Mr uncertainty. Bottom: Relative difference between the cumulative mass contained in the cluster member halos from our model and that of Be23.

In the text

Fig. 6. Redshift and F200W magnitude of multiple images in our gold category with known spectroscopic redshifts. Grey and black points represent known images with MUSE redshifts. For black squares, we also obtained NIRISS redshifts. Gold squares represent newly identified images for which we could obtain NIRISS redshift. Blue diamonds represent new images without NIRISS but with NIRSpec redshift. Images for which the redshift could not be reliably measured but with measured system redshift from other images with NIRISS or NIRSpec are marked with gold and blue circles, respectively. Each multiple image is considered once (i.e. different clumps are not considered separately). The coloured bands at the bottom of the figure represent redshift ranges where some of the most prominent emission lines fall in the wavelength range, covered by each of the three NIRISS grism filters, used in this work. The darker shade of the grism range represents a more sensitive region (see Sect. 2.2).

In the text

Fig. 7. Redshift distribution of multiple images included in our gold catalogue with known spectroscopic redshift. Images for which the spectroscopic redshift was obtained only by MUSE are shown in grey. Images for which we provided spectroscopic redshifts with JWST are shown in orange. Striped regions represent images with both JWST and MUSE redshift measurements. In the histogram, each multiple image is counted once (i.e. different multiply lensed clumps are not considered separately). The red solid line represents the expected redshift distribution of multiple-image systems, considering the redshift distribution of sources in the field and multiply imaged area in the source plane, which increases with redshift.

In the text