© The Authors 2024

1 Introduction

The James Webb Space Telescope (JWST) provides unique insight into the distant Universe thanks to its superior capabilities (spatial resolution and sensitivity) in the infrared. The power of JWST increases even more when pointing it toward a massive galaxy cluster, acting as an additional lens in front of the telescope. Images from objects behind the lens can get amplified by factors O(10) if they are on a galaxy scale, to factors O(10 000) if they are the size of stars. Through these natural pinholes to the distant Universe we can see objects that are a few to ≈10 magnitudes fainter than it would have been observed without the gravitational lens. That is, one can observe objects that (without the aid of lensing) would have had an apparent magnitude AB ≈40 in JWST’s IR bands with an integration time of just a few hours. To put this in perspective, for the most extreme magnification factors expected (~10₄, similar to the magnification estimated for Earendel, the farthest star detected through gravitational lensing and at z ≈ 6, and other distant stars at z > 4 Welch et al. 2022; Meena et al. 2023; Furtak et al. 2024) the combined effect of JWST plus the gravitational lens is equivalent to a JWST-like telescope but with a 600 meter diameter mirror. Precise modeling of these gravitational lenses is needed in order to identify and characterize these regions of maximum magnification, which in turn is needed in order to interpret the observations.

One of these lenses is the galaxy cluster MACS J0416.1-2403 (MACS0416 hereafter; Ebeling et al. 2001; Mann & Ebeling 2012) at redshift z = 0.396. This cluster is one of the six clusters studied in detail by the Hubble Space Telescope (HST) as part of the Hubble Frontier Field (HFF) program (Lotz et al. 2017), and selected based on the presence of several strongly lensed galaxies (Zitrin et al. 2013). As part of this program, the central ~10 arcmin² region was observed in wavelengths ranging from 0.45 µm to 1.6 µm and to a depth of ~28.5 mag in the visible and IR bands. The HFF data from MACS0416 have been extensively used by lens modelers (Jauzac et al. 2014, 2015; Johnson et al. 2014; Diego et al. 2015, 2024b; Kawamata et al. 2016; Hoag et al. 2016; Caminha et al. 2017; Richard et al. 2021; Bergamini et al. 2023) in the search for high-redshift galaxies (McLeod et al. 2015; Oesch et al. 2015; Kawamata et al. 2016; Ishigaki et al. 2018), and for the discovery of some of the first lensed stars at z > 0.9 (Rodney et al. 2018; Chen et al. 2019; Kaurov et al. 2019). The area observed by HST around MACS0416 was doubled thanks to the Beyond Ultra-deep Frontier Fields and Legacy Observations (BUFFALO) program (Steinhardt et al. 2020), although with shallower observations than in the HFF program.

More recently, JWST has observed MACS0416 as part of the Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS, Windhorst et al. 2023), extending the wavelength coverage to λ ≈ 5 µm, and with a superior resolution and sensitivity in the IR range than previously observed by HST. The PEARLS observations of MACS0416 consisted of three epochs, which, together with an additional fourth epoch from the CAnadian NIRISS Unbiased Cluster Survey (CANUCS, Willott et al. 2017), allowed for several transient events to be detected, out of which ≈10 are suspected to be extremely magnified objects (EMO), in particular EMO stars, undergoing microlensing events near the regions of maximum magnification (Yan et al. 2023). These events are mostly due to the existence of two galaxies (dubbed Warhol and Spock) at z ≈ 1 that are crossing the cluster caustic. One of these events was studied in detail by Diego et al. (2023c) and corresponds to a double star being magnified by a millilens with a mass smaller than 2.5 × 10⁶ M_⊙, making it the smallest structure whose mass has been estimated through gravitational lensing above z = 0.3, and with implications for dark matter (DM) models. Earlier HST observations have revealed at least nine additional transient events in the Spock and Warhol galaxies (Rodney et al. 2018; Chen et al. 2019; Kaurov et al. 2019; Kelly et al. 2022), all of them consistent with being microlensing events, bringing the total number to a record 19 alleged microlensing events in a single cluster field. Because of this wealth of flickering events, MACS0416 is referred to as the Christmas tree galaxy cluster and is so far the most remarkable cluster in terms of the number of transient events, allegedly due to microlensing.

The low redshift of the two aforementioned caustic crossing galaxies, Warhol and Spock, correlates with the large number of microlensing transient events of distant stars detected in these galaxies. Without lensing, stars in these galaxies would have apparent magnitudes in the range 36–40. With magnification factors between 10³ and 10⁴ , some of these stars could be detected by JWST (Miralda-Escude 1991). At lower redshifts, lensing can more easily promote fainter but more abundant stars above the detection limit (typically AB ~29, see also Golubchik et al. 2023). At z ≳ 2, detecting such faint stars would require more extreme and less likely magnification factors to compensate for the greater luminosity distance. Irrespective of the redshift, the stars that one expects to detect are always very luminous (L > 104 L_ø, see for instance Kelly et al. 2018). Among these, hot blue stars can have luminosities exceeding 10⁶ L_⊙. Red (cooler) stars are typically less luminous but are more often variable, facilitating their detection as transients.

These transient events are all found near the critical curve (CC) of the cluster and are the result of a temporary alignment between stars responsible for the intracluster light (ICL), and stars in the background galaxy, which increase their apparent brightness by a few magnitudes during the microlensing event. Other small structures in the cluster, such as globular clusters (GCs), can produce longer-lasting millilensing magnification when aligned with background stars. One such millilensed star, dubbed Mothra, was discovered in MACS0416 (Diego et al. 2023c). JWST’s superior sensitivity and resolution allow the best view of both the ICL and the GCs in MACS0416, and hence valuable insight into the observed micro- and millilensing events can be gained, as well as into the expected correlation between the spatial distribution of these events and the distribution of ICL and GCs.

Earlier results based on observations and N-body simulations have shown a tight correlation between the number of GCs in a halo and its virial mass (Blakeslee et al. 1997; Harris et al. 2017; Burkert & Forbes 2020; Valenzuela et al. 2021; Berkheimer et al. 2024). This relation is satisfied also when relating the total virial mass of a halo and the total mass contained in the GCs (Dornan & Harris 2023). This relation has been verified multiple times at low redshifts (z ≲ 0.2) with observations, although in most cases, the masses are derived from kinematics, X-ray or weak lensing. On galaxy cluster scales and at higher redshifts, the first results from JWST on SMACS0723 (z = 0.39) have shown the strong correlation between the distributions of mass, GCs, and ICL (Lee et al. 2022; Faisst et al. 2022; Montes & Trujillo 2022; Diego et al. 2023b). This is observed also in highly irregular galaxy clusters such as A2744 at z = 0.308 Harris & Reina-Campos (2023). MACS0416, with its wealth of strong lensing constraints and relatively deep observations, allows us to explore the correlation between the distribution of mass and the spatial distribution of GCs in greater detail.

Properly characterizing the abundance and distribution of GCs in gravitational lenses is also important in order to correctly interpret past and future transient events near the region of high magnification, since GCs can introduce significant distortions into the magnification pattern.

The study of the topics described above (high-z stars, magnification and time delay near CCs, or the relation between total mass and distribution of GCs) is only possible if accurate lens models are available. Lens models are also needed to correctly interpret these events, for instance by setting a constraint on their physical size that can later be used to identify some objects near CCs as stars.

This paper presents a new lens model derived using new lensing constraints from the latest observations of MACS0416 with JWST. The paper focuses on the lensing properties (magnification and time delay) of the region of maximum magnification around galaxies that are crossing caustics in order to provide the proper context for the observed transients and future ones that will continue to be discovered in this cluster. We also pay special attention to the population of GCs found in this cluster. Their distribution can be used as a mass proxy in regions of the lens plane where lensing constraints are rare or absent. Also, GCs can act as millilenses, so it is important to understand their distribution and be able to estimate the probability of the alignment between GCs and the background luminous stars that can be magnified by these GCs. The paper is organized as follows. Section 2 describes the observations carried out during the PEARLS program. The lensing constraints (including the new ones discovered by JWST) used to derive the lens model are discussed in Section 3. Section 4 presents the lens model, and Section 5 focuses on the lensing implications for galaxies crossing caustics and the possibility of detecting stars in lensed arcs. Section 6 discuses GCs as mass tracers and lensing perturbers. Section 7 considers in more detail the complexity of the lens model around the Spock arc, one of the transients in Yan et al. (2023), and the lensing distortions introduced by GCs. Finally, Section 8 summarizes our conclusions.

We adopt a standard flat cosmological model with Ω_M = 0.3 and h = 0.7. At the redshift of the lens (z = 0.396), and for this cosmological model, one arcsecond corresponds to 5.39 kpc. Unless otherwise noted, magnitudes are given in the AB system (Oke & Gunn 1983). Common definitions used throughout the paper are the following. “Macromodel” refers to the global lens model derived in Section 4. The main caustic and main CCs are from this model. “Micromodel” means perturbations introduced by microlenses in the lens plane (stars and remnants), and “microcaustics” means regions in the source plane of divergent magnification associated with microlenses. The term millilens is used for objects in the mass range 10⁴–10⁶ M_⊙ (typically GCs and dwarf galaxies in the galaxy cluster).

2 New JWST observations of MACS0416

MACS0416 was observed in three different epochs by JWST during October 2022 and February 2023 as part of the PEARLS program. Observations were done in 8 filters, F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W. Integration times in each filter were between 2920 seconds and 3779 seconds reaching 5σ limiting magnitudes between 28.45 and 29.9 per epoch and for point sources. Details of the observations for this cluster can be found in Windhorst et al. (2023). The data processing for this cluster is described in detail by Yan et al. (2023). The only special steps involved removing artifacts (1/ f noise, pedestal, and wisps, that were removed with ProFound Robotham et al. 2017, 2018) and aligning to Gaia DR3. As part of the process, the available data from the HFF and BUFFALO programs with HST (including data taken with the relevant bluer filters, F435W, F606W, and F814W) were realigned to the same Gaia data. In the central region of the cluster, HST data have similar depth to JWST thanks to the considerably longer exposures, but it proves very valuable because it extends the wavelength coverage down to 0.4 µm, providing useful and complementary photometric measurements at these wavelengths.

The three epochs from the PEARLS program are taken at different position angles and cover an area of ≈2.8 × 2.8 arcminute² in the cluster core field, with additional area being covered further away by the second module. This work focuses only in the central cluster core field which is the only one containing lensing constraints. A portion of the area observed by JWST in this cluster core region is shown in Fig. 1. The area shown in this figure contains all multiple images identified in the data and marked with circles. This image (and others in this paper) is a color composition made after combining previous HST data from the HFF program with the new JWST data. The blue channel is a combination of HST’s F435W, F606W, F814W, and JWST’s F090W. The green channel combines HST’s F105W, and F160W with JWST’s F115W, F150W, and F200W. The red channel combines JWST’s F277W, F356W, F410M, and F444W.

3 Lensing constraints

Existing strong-lensing constraints comprise 77 multiply lensed galaxies that produce 223 individual images. Diego et al. (2024b) (their Table A.1) gave a complete list compiled from the literature (Jauzac et al. 2014; Johnson et al. 2014; Diego et al. 2015; Kawamata et al. 2016; Caminha et al. 2017; Richard et al. 2021; Bergamini et al. 2023). All galaxies in that list have spectroscopic redshifts. Additional lensed system candidates that lack spectroscopic redshift were identified in HST images in the references above. This set is enlarged with more system candidates identified with the new JWST data. We refer to this second set of constraints as the non-spectroscopic sample. All available lensing constraints are listed in Table B.1. These include some of the missing counterimages in the spectroscopic sample, identified in the new JWST data, and marked with † in the table.

Identification of new lensed candidates in JWST images is done in the usual manner. Families of multiple images with compatible photo-z (when available), color, morphology, parity, and consistency with the previous lens model by Diego et al. (2024b) are identified in the new JWST images. Based on the large mismatch between the geometric and photometric redshifts in some systems, we expect a few of these new candidates (≈4) to be false positives (see discussion at the end of Section 4.3). These low confidence systems are not used as lensing constraints. In Table B.1 the confidence of all systems is ranked with A for the spectroscopically confirmed systems, B for systems without spectroscopy but reliable according to the criteria listed above, C for systems that are less reliable (for instance, lacking sufficient morphological information), and D for the least reliable systems (for instance, the lens model is not accurate enough around these systems). Images with score C and D are not used as lensing constraints but are added here as potentially interesting targets for spectroscopic follow-up.

All arcs listed in Table B.1 are shown in Fig. 1. The systems with spectroscopic confirmation are marked with white circles, while system candidates without spectroscopic confirmation are marked with yellow circles.

Fig. 1 Color image of MACS0416. The scale and orientation are indicated. White circles mark previously known counterimages with spectroscopic redshifts from Richard et al. (2021); Bergamini et al. (2023). Yellow circles mark reliable (rank B) HST- and JWST-identified counterimage candidates that lack spectroscopic confirmation. Magenta and red circles mark the least reliable counterimage candidates with rank C and D respectively.

4 Lens model

4.1 Lens model method

The lens model is based on the code WSLAP+ (Diego et al. 2005, 2007; Sendra et al. 2014; Diego et al. 2016a). This is a hybrid type of model that combines a free-form decomposition of the smooth, large-scale cluster component with a small-scale contribution from cluster galaxies. The code combines weak (when available) and strong lensing in a natural way with changes in the large-scale and small-scale components equally affecting the strong- and weak-lensing observables. For this work, only strong-lensing constraints are used.

As a brief description, WSLAP+ solves a system of linear equations that can be represented in compact form as $$\Theta =\text{Γ}X,$$(1)

where Θ is a vector containing the observed sky-plane positions of the multiple images, X is a vector containing the unknown source-plane positions of the background sources and the model mass distribution in the lens-plane, and Γ is a matrix representing the ray-trace calculations that predict observed positions, Θ′, given X for a fiducial mass. The solution X is found by iteratively minimizing a function of Θ′ – Θ with the constraint X > 0 (that is, masses and positions are positive, where the origin for the positions is the bottom-left corner of the area being considered, hence by definition all coordinates within this area are positive; Diego et al. 2005). In practice, the lens-plane mass is divided into N_c grid points or cells. These cells are associated with a smooth matter distribution representing the sum of DM, intra-cluster medium (gas and stars), and any other smooth component. These cells have preselected positions and sizes chosen by the modeler. Each of the N_s source galaxies is represented by two variables in X, and WSLAP+ allows “layers” representing independent mass normalizations for different groups of member galaxies. The number of layers, N₁ , can vary between 1 when all member galaxies are included in one single layer, to the total number of member galaxies. In practise, the number of layers is usually between 1 and 5. With these components, vector X has dimension $$N_{\text{X}}=N_{\text{c}}+2N_{\text{s}}+N_{\text{l}}.$$(2)

If there are N_s1 strong-lensing observables, the vector Θ has dimension N_Θ = 2N_sl because each multiple image contributes two position constraints, θ_x and θy, and distances (or redshifts) are considered known. Hence, the matrix Γ has dimension N_X × 2N_sl and is also known.

Our MACS0416 models used different sets of constraints and different numbers and distributions of cells as described below. All models used three layers in addition to the dark-matter layer. Layer 1 contained the two brightest cluster galaxies (BCGs); Layer 2 contained the remaining most-prominent cluster galaxies (444 galaxies brighter than magnitude 23.5) with the bulk of them identified as members of the red sequence (414 galaxies, out of which 298 are also spectroscopically confirmed by MUSE) and an additional 30 galaxies confirmed spectroscopically by MUSE (Caminha et al. 2017) but not included in our red sequence sample. For the spectroscopic galaxies we considered the range 0.39 < z < 0.41 in order to account for the large velocity dispersion (and likely elongation along the line- of-sight) of this cluster; and Layer 3 contained a few z = 0.12 foreground galaxies. This third layer contains only 3 galaxies but is included here because these foreground galaxies are near elongated arcs and play a role in their distortion. We expect a negligible contribution from the smooth component associated with these three galaxies so the smooth component at z = 0.12 is not included in our model. We also expect a negligible contribution from small member galaxies not included in Layer 2. We assigned each galaxy’s mass distribution to be proportional to the galaxy’s surface brightness in JWST’s F356W filter. This filter has the best sensitivity in the long wavelength channel of JWST. The proportionality constant could differ for the three layers and was readjusted as part of the optimization process. The large scale mass around the member galaxies is modeled by the Gaussians distributed in the grid. Diego et al. (2005) gave a complete description of the iterative solution process, and Sendra et al. (2014) described the model’s convergence and performance based on simulated data. For a direct comparison with other lens modeling techniques see Meneghetti et al. (2017).

4.2 Lens model with spectroscopic systems only

First, we derive a solution using only image systems with spectroscopic redshifts as constraints. This initial trial used a regular 20 × 20 grid for the smooth mass component. At this resolution, a cell has a dimension of 45.3 kpc on a side. Adopting a regular grid is equivalent to assuming no prior information for the mass distribution because all grid points contribute similarly to the system of equations in Eq. (1). The next solution used an adaptive grid with more grid points around the densest regions of the first solution, to increase the resolution in the central region. The cell size varies ranging from ≈15 kpc near the two BCGs to ≈100 kpc in the outer regions of the field. This applies a prior to the mass distribution in which more mass is expected in regions where the distribution of grid points is denser. The prior helps to stabilize the solution and reduces the scatter from the choice of the initial grid. Because the largest source of variability in the derived models is the choice of grid, we obtained multiple solutions with different grid configurations. Solutions obtained with an adaptive grid having a number of grid points smaller than the regular grid number of grid points (N < 400) do not reproduce some arcs as well as solutions obtained with an adaptive grid having a number of grid points N > 600 result in artifacts (spurious large scale mass structures) in the regions of the cluster with no lensing constraints. We find that solutions obtained with a multiresolution grid of 495 grid points offers robust and more stable solutions, that do no introduce artifacts in the unconstrained boundaries of the field of view. Small variations (δN ≈ 20) in the number of grid points around N ≈ 500 result in almost indistinguishable solutions.

4.3 Geometric redshifts

Once a lens model is created from the set of spectroscopic systems, it predicts “geometric redshifts” (geo-z or z_geo) for lensed systems with and without spectroscopic redshifts. We used the solution obtained with the multiresolution grid described above to derive geo-z values for all arcs following the same procedure as Diego et al. (2023b,a). For systems with spectroscopic redshifts, the lens model predicts z_geo close to z_spec as shown in Table B.1. Therefore z_geo is likely to be accurate for candidate image systems that lack spectroscopic confirmation and that appear near regions where a spectroscopic system is acting as a constraint.

As a consistency test, we compare geo-z with photometric redshifts for systems that lack z_spec, and results are shown in Fig. 2. The photometric redshifts were obtained from or following the methods of Adams et al. (2024). In brief, the photo-z’s were calculated from JWST photometry, extracted in 0716-radius apertures, and aperture corrected based on PSF models. The photometry was then run through the EAZY phot-z code (Brammer et al. 2008) using the template combination from Larson et al. (2023). In MACS0416, these photo-z’s recovered 235/261 (90%) of known MUSE spectroscopic redshifts (Caminha et al. 2017). Approximately half of our lensed sources are in the Adams et al. (2024) catalogs, while the rest are blended or contaminated with foreground cluster galaxies or the intra-cluster light. These required manual aperture placement to obtain z_phot estimates.

To first order, there is a good correspondence between all three redshifts, but there are several examples of significant mismatches. (All results for all systems are given in Table B.1.) Many of the outliers in the photometric redshift concentrate below z_phot ≈ 1. This may be simply because these sources are fainter or contaminated by intracluster light as suggested by the relatively large concentration of outliers near z_phot ≈ 0 4 (the redshift of the cluster). The number and distribution of mismatches between geo-z and photo-z are comparable to the mismatch between photo-z and the spectroscopic redshift. Focusing on the biggest outliers, there are ≈4 red points (geo-z) in Fig. 2 that are more offset (when compared with the z_phot–spect-z outliers or blue points outside the perfect correlation). These could be false positives in our sample (that is, incorrectly identified lensed systems). Two of these systems are already ranked D in Table B.1 as unreliable. These are systems 113 (z_geo = 0.69, z_phot = 2.02) and 114 (z_geo = 16.67, z_phot = 0.51) in Table B.1. A third suspicious system is 116 (z_geo = 16.69, z_phot = 3.11) which is ranked C (and hence not used as a constraint) and with very discrepant geo-z and photo-z. System 81 (z_geo = 11.3, z_phot = 2.56), in this case two counterimages form a pair of symmetric features near the expected position of the CC and are ranked B and used as constraints in the full lens model. The third counterimage candidate is very faint, ranked C and hence not used as a constraint. Moreover, when using only the two images ranked B, we derive $z_{geo}={2.8}_{-1.1}^{+4.4}$, in better agreement with z_phot, thus confirming our low confidence on the third image of system 81. Even if the adopted redshift for System 81 is wrong, the fact that we only use the relatively close pair 81.a and 81.b as constraints results in very little weight in the final model, since this pair of images is already very close in the source plane for a wide range of redshifts, with a separation of just ≈1.5 arcsec in the image plane. Despite the suspicion that the high geo-z of System 81 may be biased (when considering the third counterimage candidate), we opt for including this system with the high-z estimate of the redshift, in order to maintain consistency with the other rank B images. We conclude from this comparison that for systems ranked B, z_geo is a reliable estimate of the redshift of the system. Errors for the geo-z estimations are listed in Table B.1. A more detailed study of the errors in geo-z from estimates obtained with WSLAP+ can be found in Diego et al. (2023a) where it is demonstrated, with a bootstrap analysis, how the errors are representative of the true errors (see their figure 8), except for systems in the outskirts of the lens where the number density of lensing constraints is low.

Since z_geo are by construction correlated with a previously derived lens model (and may be affected by biases in that model), using high-fidelity z_phot for the lensed galaxies is in general a preferable option than using z_geo for the same galaxies when deriving lens models. However, for our particular case, we are interested in having the most reliable redshift that allow us to predict magnification and time delays for counterimages of sources at those redshifts. Given the large number of spectroscopic systems available in MACS0416, the uniform and dense spatial distribution of the lensed images, and the relatively poor photometry (most counterimage candidates are faint) affecting the z_phot estimates, z_geo estimates are more reliable for this particular case and are used in the next sections to derive the final lens model and the derived lensing properties of all arcs. This option is also the safest, since an inaccurate z_phot can introduce a large bias in the derived model, while z_geo estimates are expected to introduce smaller biases since they produce lens models which are consistent with the model derived using spectroscopic systems only.

As a final note, during the last phases of the refereeing process, Rihtaršič et al. (2024) publish new spectroscopic redshifts for 13 of the systems candidates in Table B.1. For these 13 systems, we find that our previously derived geometric redshifts are accurate estimates of the true redshift (the comparison is shown in Figure B.1) with a dispersion in the error of 6.8%, but with a median in the distribution of errors of less than 1% error.

Fig. 2 Comparison of photometric redshifts with spectroscopic (blue symbols) or geometric (red symbols) redshifts for the lensed systems in Table B.1. For the photometric redshifts, we chose the value (from the multiple counterimages of the same system) closest to the spect-z or geo-z. The dotted line shows equality. The comparison between spect-z and geo-z is not shown but it follows almost perfectly the dotted line (see Table B.1).

4.4 Lens model with all systems

A final lens model comes from adopting z_geo for candidate image systems that lack zspec and with rank A or B. We refer to this sample of lensing constraints as the full sample. The model, based on the 495 multiresolution grid and the full sample, does not differ much from the previous one because the model was already well constrained with the spectroscopic sample and the new systems have, by construction, redshifts consistent with the spectroscopic model. However, the new constraints help to reduce uncertainties, particularly in regions where the density of spectroscopic systems is low. Regarding the precision in the predicted positions of the counterimages, this model has a root mean square (RMS) between the observed and predicted positions in the image plane of 1.07″. When considering only the images with spectroscopic redshift, the RMS is 0.91″, while images with geo-z have an RMS of 1.36″. These values are typical of lens models derived with WSLAP+.

Fig. 3 shows the resulting mass profiles centered on each of the two groups (BCGs). The biggest differences between models come from the choice of grid configuration rather than from the set of constraints. Nevertheless, differences are minuscule in the region constrained by lensing.

All models show a striking similarity between the NE and SW halos, consistent with other recent models Bergamini et al. (2023); Diego et al. (2024b). The profiles can be well reproduced by the same generic Navarro–Frenk–White profile (Navarro et al. 1997; Hernquist 1990; Zhao 1996; Wyithe et al. 2001; Nagai et al. 2007): $$\rho (r)=\frac{{\rho }_o}{{\left(r/r_s\right)}^{\gamma }[{\left(1+{\left(r/r_s\right)}^{\alpha }\right]}^{(\beta -\gamma )/\alpha }}$$(3)

where γ, α, and β are the inner, intermediate, and outer slopes, and rs is the scale radius. In all models, the inner 10 kpc is dominated by the baryonic component of the BCGs surrounded by a shallow convergence profile out to 200 kpc. This type of shallow profile extending beyond 100 kpc has been observed in other merging systems similar to MACS0416, for example A370 (Lagattuta et al. 2017; Diego et al. 2018b).

In our case, the maximum radius constrained by lensing is 200 kpc, and the mass interior to that M(<200 kpc) = 1.72 × 10¹⁴ M_⊙ for the NE halo and M(<200kpc) = 1.77 × 10¹⁴ M_⊙ for the SW halo. In a larger volume, Bergamini et al. (2023) found masses (from their best-fitting model) M(<400 kpc)≈3.4 × 10¹⁴ M_⊙ for both halos. At this distance, we find M(< 400 kpc) = 3.81 × 10¹⁴M_⊙ and M(< 400kpc) = 3.93 × 10¹⁴ M_⊙ for the NE and SW halos, respectively. Earlier results derived using the same WSLAP+ code, BUFFALO data from HST, the spectroscopic systems only, and including weak-lensing measurements, obtained a mass M(<400 kpc) = 3.6 × 10¹⁴ M_⊙ (Diego et al. 2024b) for their best model, also in reasonable agreement with our results.

For the compact components we find a mass of 2.28 × 10¹² M_⊙ for Layer 1 (BCG1+BCG2), 1.03 × 10¹³ M_⊙ for Layer 2 (red-sequence and spectroscopic member galaxies) and 3.31 × 10⁹ M_⊙ for Layer 3 (foreground galaxies).

We also find very similar magnification patterns to earlier models. The CCs for three different source redshifts are shown in Fig. 4. In the same image, some of the arcs that are crossing caustics at their respective redshifts are shown in the zoomed rectangular regions. The white (rank A) and thin yellow (rank B) circles mark the position of the lensing constraints from these arcs used to derive the lens model. The thick yellow circles mark the positions of unresolved sources close to the CC. The source labeled D21_W1 is a peculiar transient presented in Yan et al. (2023) and discussed in more detail later. Since the magnification predicted by the lens model at these positions is very large (typically hundreds), the unresolved nature of these sources at these magnification factors imply that these objects are good candidates to be EMO stars. Some of the objects marked with a thick yellow circle have indeed shown variability, which is a fundamental prediction for EMO stars moving across a web of microcaustics. JWST is expected to detect many more microlensing events in this cluster, in particular from the low redshift galaxies Warhol and Spock Diego et al. (2024b); Yan et al. (2023). Section 5 discusses in more detail the lensing properties of the highly magnified arcs shown in panels A–G of Fig. 4.

We conducted a search for dropout galaxies intersecting the CC at high redshift but found no convincing candidate although a few very faint and stretched features are found, which are high- redshift galaxy candidates intersecting the cluster caustic. One such candidate is shown in panel F of Fig. 4 (system candidate 116 and barely visible in JWST images). The CC in that panel is at z = 6.15, suggesting that this arc must be at z > 6

The solution presented above relies on geo-z estimates. Earlier geo-z predictions based on WSLAP+ models have been confirmed with subsequent spectroscopic follow up. For instance System 25 in Diego et al. (2015) predicted at z = 1.05 and later measured at z = 0.94 (System 12 in Caminha et al. 2017), System 10 in Diego et al. (2016b) at a predicted redshift of z = 0.78 and later measured at z = 0.73 (B4 in Table 3 of Caminha et al. 2016), systems 7 and 19, proposed in Diego et al. (2018b) as being part of the same system with a predicted redshift of z = 2.8 and later confirmed by Lagattuta et al. (2017) as being indeed the same system and at z = 2.7512, or System 11 in Diego et al. (2020) predicted at 2.2 ≲ z ≲ 3.1 (see their Fig. 5) and measured at z = 2.1887 in Caminha et al. (2023) (their System 23). When reliable photo-z are not available because of insufficient photometric information or the sources are too faint, geo-z are in general a safer option since they are derived from a model that relies on systems having spectroscopic redshifts. Because of this, solutions obtained with spectroscopic plus geo-z systems are often very similar to solutions obtained with spectroscopic systems only. The similarity between the profiles discussed earlier demonstrates this (Figure 3). Earlier results based on WSLAP+ and JWST data from PEARLS also show a high degree of similarity between both type of solutions (Diego et al. 2023a). In particular, in Diego et al. (2023a) it is shown with a bootstrap method how the errors in the geo-z obtained by WSLAP+ are representative of the true errors, except in regions where the spectroscopic constraints are the most scarce, such as in the outskirts of the region constrained by strong lensing. Given the very high number density of spectroscopic systems surrounding the non-spectroscopic systems (see Figure 1) we expect the geo-z predictions for MACS0416 to be better than in earlier work, and hence very reliable.

Fig. 3 Mass profile of the lens models. Blue curves correspond to profiles centered on the NE BCG while red curves are for profiles centered at the SW BCG. Solid lines represent the models derived using only spectroscopically confirmed (rank A in Table B.1) galaxies as constraints and an adaptive grid with 604 grid points for the smooth DM. Dotted lines show the corresponding solution for an alternative model with the same constraints (rank A) but with only 495 grid points. Dashed lines represent the same 495-point grid but with the full sample of lensed systems ranked A and B (geometric redshifts) as constraints. The mass profiles are given in dimensionless κ units assuming a source at z = 3. The long-dashed line shows a gNFW profile with parameters (Eq. (3)) γ = 0.9, α = 2, β = 3, and scale radius rs = 640 kpc. The two vertical black segments indicate the distance interval where the profiles are constrained by lensing. The inset shows the integrated mass as a function of radius in units of 1014 M⊙ for the same three profiles shown in the main panel. The black long-dashed line is the integrated profile when the center is chosen as the midpoint between the two BCGs (RA = 4hr16m08s.428, Dec = –24°04′21″.0).

5 Magnification and time delay for known caustic-crossing galaxies

Yan et al. (2023) found eleven transients in three highly magnified galaxies (Warhol, Spock, and Mothra) using the new JWST PEARLS+CANUCS data. Most of these transients are consistent with being microlensed stars. Previously, nine transient events were identified in Warhol and Spock in HST observations Rodney et al. (2018); Chen et al. (2019); Kaurov et al. (2019); Kelly et al. (2022). The lensed star in the Mothra arc was also seen in previous HST data, but no flux variations were noticed in the HST data so prior to JWST observations it was not recognized as a lensed star. All these transients are expected to be very luminous stars (L > 10⁴ L_⊙) that temporarily increase their observed flux as a microlens from the galaxy cluster intersects the line of sight. Interpreting these events requires knowing their magnifications in the absence of microlensing, that is, the macromodel magnifications at the stars’ positions. The lens model predicts the magnifications, but the accuracy depends on whether sufficient local constraints on the model are available. A feature of highly magnified stars is that they usually produce two counterimages, one on each side of the CC. Given the difference in paths and lensing potential, the two counterimages suffer different time delays, which are predicted by the lens model. The delays are important when a star exhibits intrinsic variability (a SN being the most extreme example) to predict the time when the trailing counterimage will appear.

Rather than presenting lens-model-predicted magnifications and time delays for specific events, it is more useful to present predictions for a range of distances from the CC because new transients can appear anywhere along these arcs. This needs to be done for each lensed galaxy individually because differences in redshift and lensing potential change the model predictions. Also, this approach reduces the uncertainties due to imperfect lens modeling that often places the CC at positions a fraction of an arcsecond away from the true (unknown) position. Even though these offsets are small, they are sufficiently large for accurately predicting the magnification of unresolved sources very close to the CC. In these situations, a fit to the expected scaling laws combined with a model independent estimation of the distance to the CC (for instance from symmetry arguments between a pair of counterimages) results in more accurate predictions for the magnification and time delays. In this section, we present models for the magnification and time delay along 12 prominent arcs that are crossing the cluster caustic at their corresponding redshifts.

The magnification along an arc and as a function of distance to the CC can be fitted with the canonical law for CCs µ_th = µ₀ /d_cc. Here µ₀ is a constant (expressed in arcseconds) that depends on the redshift of the source and position along the CC, and dcc is the distance to the CC expressed in arcseconds. The time delay has a similar scaling with distance, ∆T_th = T₀ * d_cc, where T₀ is another constant expressed in years per arcsecond. For a given source redshift, the constant T₀ varies also from position to position along the CC at that redshift. The time delay ∆T is measured between two counterimages of the same source position on each side of the CC and separated by a distance d = 2d_cc. These scalings in µ and ΔT are expected for a lens with an isothermal-sphere density distribution, and they normally reproduce well the observed magnifications and time delays in real lenses.

Departures from the linear scalings with d_cc are seen when substructures perturb the deflection field near the CCs. To account for these possible perturbations, we computed the deviation of the observed magnifications from these simple scalings. We did so at three fixed distances from the CC, 0.″6, 1″, and 2″ and on either side of the CC, that is, for images with positive and negative parity. The deviation from the scaling is computed as: $${\sigma }_{\mu }={\displaystyle \sum _i{\left(\frac{\left|{\mu }_{\text{pred}}-{\mu }_{\text{th}}\right|}{{\mu }_{\text{th}}}\right)}_i}$$(4)

for the magnification, and similarly for the time delay $${\sigma }_{\text{Δ}T}={\displaystyle \sum _i{\left(\frac{\left|\text{Δ}{T}_{\text{pred}}-\text{Δ}{T}_{\text{th}}\right|}{\text{Δ}{T}_{\text{th}}}\right)}_i}.$$(5)

The index i runs over three separations (i = 1, 2, 3 for d_cc = 0″.6, 1″, and 2″, respectively). In the equations above, the subscript "pred" refers to the predicted magnification while the subscript "th" refers to the expected scaling laws. For convenience, we took µ₀ as the magnification measured at 1″, and similarly we took T₀ as the time delay measured when d_cc is 1″. With this definition, the terms where i = 2 (that is, where d_cc = 1″) in the equations above are identically zero. Because the magnification can scale differently depending on the parity of the image, for each arc we derived two normalization factors, ${\mu }_0^{+}$ and ${\mu }_0^{-}$ for the sides with positive and negative parity, respectively (both following laws $\mu ={\mu }_0^{+}/d_{\text{cc}}$ and $\mu ={\mu }_0^{-}/d_{\text{cc}}$). The time delay is simpler because the image with positive parity always arrives first, so the relative time delay is given by the delay of the image with negative parity with respect to the image with positive parity. Hence only one normalization factor, T0, is needed for each arc.

Using the three fixed relative separations to the CC at the position of each of the 12 arcs, we computed ${\mu }_0^{+},{\mu }_0^{-},T_0,{\sigma }_{{\mu }^{+}},{\sigma }_{{\mu }^{-}}$, and σ_ΔT for the 12 caustic-crossing arcs. The results are summarized in Table 1, and the 12 arcs are shown in panels A–G of Fig. 4. The values listed in this table are derived from measured distances to the lens model CC and not from observed sources along these arcs. Hence, they are not affected by the small offsets (typically a fraction of an arcsecond) between the predicted positions of the CC or lensed arcs and the true (unknown in the case of the CC) positions. Very small and local modifications in the lens model (too small to be captured by the grid of Gaussians in our model) can place the predicted CC or arcs at the true position while leaving the values listed in Table 1 virtually unchanged. When computing distances, these are measured in a direction perpendicular to the CC. The arcs studied in this section are all at directions which are almost orthogonal to the CC, so distances computed this way should be accurate enough. A simple estimation of the uncertainty in the magnification and time delay can be obtained from the σ values listed in Table 1. Relative uncertainties for the magnification are between a few percent to ≈20%. For time delays, the uncertainty in σ is larger, specially for arcs with small time delay amplitudes, T₀, for which the denominator in Eq. (5), can be very small. These uncertainties are for one lens model only and does not account for uncertainty in the lens model itself which typically is larger and can exceed 50% very close to the CC (see for instance Meneghetti et al. 2017). Since the new candidate systems lack spectroscopic confirmation we avoid the lengthy exploration of uncertainties in the lens model, which has been done in the past for the spectroscopic sample.

Below we briefly discuss some of the relevant arcs in Table 1 and shown in Fig. 4:

• System 1 corresponds to the Warhol galaxy (panel C in Fig. 4). This galaxy is extreme on several levels. It has the lowest redshift among all the lensed galaxies and is the closest to the BCG. The BCG’s proximity gives System 1 the most (alleged) microlensing events: five seen by HST plus seven seen by JWST. Microlensing events are relatively short-lived, with durations between a few days to a few weeks, depending on several factors such as relative velocity, radius of the background star, mass of the microlens, or redshift of the background star (Kelly et al. 2018; Diego et al. 2018a, 2024a). It also has the largest time-delay factor with T₀ = 3.14 years (Table 1, again mostly due to proximity to the BCG, which results in a large gradient in the lensing potential and hence a large contribution from the Shapiro term of the time delay, see Eq. (6) below). This system contained the brightest event found by Yan et al. (2023), denoted D21_W1, which reached AB mag 27.1 in F277W during PEARLS’s epoch 2. This event was located at d_cc = 0″.43 in the negative-parity region, where the lens model predicts macromodel magnification µ = 28. The unusual brightness of the event at such a small macromodel magnification suggests the microlensing event may not have been a typical one. Section 7 discusses some possibilities.

• System 2 corresponds to the Spock galaxy (panel D in Fig. 4). Rodney et al. (2018); Diego et al. (2024b) described the lensing properties of this arc and discussed the difficulty of modeling the CC around the arc’s position. Near this arc, there are several cluster galaxies with slightly different red-shifts (a bimodal distribution), there is a z = 0.12 foreground galaxy to the NW, and there are no other constraints for z = 1 . These properties make modeling this portion of the lens plane challenging and uncertain. Rodney et al. (2018) described how different models predict CCs that intersect the elongated Spock arc in different positions. There may be more than one caustic that crosses this arc, and if so, the distribution of magnifications and time delays would depart from the assumed scaling laws. (Section 7 has further discussion.) This departure would be more evident on the side with negative parity, where member galaxies contribute more to the deflection, as shown by the large value of ${\sigma }_{{\mu }_0}^{-}$ in Table 1. The time delay would also be affected, resulting in a poor fit to the ΔT = T₀ * d_cc scaling law.

• System 6 is the most prominent arc in MACS0416 and one of the longest (panel B in Fig. 4). No evidence of individual lensed stars has been found so far, but this arc is a prime target for future observations. Its higher redshift (z = 1.895) makes it more challenging to identify lensed stars than in the Spock or Warhol galaxies because larger magnification is required to make any given star detectable. The macromodel magnification is also relatively modest, ${\mu }_0^{-}=16$ and ${\mu }_0^{+}=27$, partially explaining the lack of success at finding EMO stars in this arc. Southwest of System 6 is System 11 (panel B in Fig. 4), which is close to the caustic and shows a clear asymmetry between the two counter-images, with the west one being partially demagnified. This is possible if a small perturber is near that counterimage. A faint galaxy is barely visible at that location (and is marked with a white arrow in Fig. 4). This galaxy is likely the one responsible for the distortion in the magnification.

• System 8 is a colorful and very red arc containing two candidate EMO stars (panel E in Fig. 4). The two candidates are very close to the expected position of the CC. An alternative is that these two objects could be GCs in the cluster. Further evidence, for instance flux variability, is needed to confirm these candidates: GCs are not expected to vary, but EMO stars do so because of both the ubiquitous microlenses and intrinsic variability. System 8 has a redshift similar to that of System 6, but it is crossing a more powerful caustic. This makes it a better candidate to search for EMO stars. Next to System 8 is another caustic-crossing galaxy, System 101, at lower redshift but without any obvious EMO stars.

• System 12 is the Mothra galaxy. This galaxy (panel C of Fig. 4) is crossing one of the most powerful portions of the caustic, making it a great candidate for detecting EMO stars. So far, though, only one EMO star (Mothra) has been detected. Mothra’s counterimage remains elusive, although it may be marginally detected in JWST data. In their detailed study, Diego et al. (2023c) found that Mothra is most likely to be a double star millilensed by an undetected halo of mass ~10⁴–10⁶ M_⊙. Diego et al. estimated Mothra to be only d_cc = 0″.07 from the CC (this conclusion is reached in a model independent way as detailed in that work). At that distance, our lens model gives macro-model µ = 972 and 716 for the counterimages with negative (the observed image) and positive parity, respectively, before accounting for milli- or micro-lensing. Both Diego et al. (2023c) and Yan et al. (2023) found that the red component of the Mothra binary is intrinsically variable, a common feature of red supergiant stars. Deeper observations should reveal the positive counterimage, which should exhibit matching time variations in the red component with a delay of just 13.5 days. (Photons from the observed counterimage arrive first.)

• System 19 is crossing a moderately powerful caustic. It shows two resolved, blue counterimages meeting at the CC (panel B of Fig. 4). The resolved nature of these counterimages and their length, L ≈ 0″.72, indicate that the source must have a size R ≈ 8 pc. It is most likely a star-forming region or a young GC. This source probably harbors a multitude of luminous stars, making this arc in principle a good candidate to detect individual micro-lensing events. However, as noted by Dai (2021), microlensing events in regions of high stellar density are harder to identify because they have smaller flux changes than the same micro-lensing of an isolated star.

• System 34 (panel A of Fig. 4) is in principle a challenging target for EMO stars because of its high redshift (Table 1).

However, this system is crossing the most powerful portion of the caustic among all galaxies listed in Table 1, making it an interesting target. On the negative side, only the outer portion of the galaxy is touching the caustic, thus reducing the probability of seeing an EMO star in this arc. The time delay for this arc is small (although with large uncertainty): an intrinsically variable source located in the middle of one of the two blue images near the CC would have its counterimage arrive just 11 days later (ignoring micro-lensing, which could be different in each counterimage).

• System 75 (panel A of Fig. 4) was barely detected in the JWST images. It was originally identified with MUSE (Vanzella et al. 2021; Bergamini et al. 2023), thanks to emission lines. An EMO star in this arc would be the most distant one ever observed, surpassing the current record holder, Earendel at z ≈ 6 (Welch et al. 2022). No EMO-candidate was found near the CC position, but the relatively large values of ${\mu }_0^{-}$ and ${\mu }_0^{+}$ offer hope for a future detection of an EMO star in this arc during a microlensing boost in its magnitude.

• System 82 (panel G of Fig. 4) has a good chance of revealing some of its brightest stars but has not yet done so. Two pairs of symmetric, unresolved counterimages are easily identified near the CC. (One close pair is inside the yellow circle in panel G.) The outer pair is at d_cc ≈ 0″.3, resulting in predicted macromodel magnifications of 89 and 99 for the counterimages with negative and positive parity, respectively. The time delay between these counterimages is 38 days. No variability was seen in the PEARLS+CANUCS epochs, so these images are probably not microlensed and therefore are magnified by factors close to the macromodel values. These magnifications are insufficient to make a red supergiant star detectable at z = 2, even if caught during an outburst. The source in this case is probably a compact star-forming region with R < 250 pc. The inner pair is at d_cc ≈ 0″.1, resulting in macromodel magnifications three times larger than for the outer pair (270 and 300) and time delays three times shorter (≈13 days). The constraint on the source size is likewise three times stronger, R < 83 pc.

• System 94 (panel E of Fig. 4) crosses the caustic at one of the weakest points. That means EMO stars can be found only very close to the CC even though the system has a low red-shift. There is a red, unresolved source halfway between the arc’s counterimages, making the source a good candidate to be a pair of unresolved counterimages from a star located <1 pc from the cluster caustic. In this case, the net magnification (from both counterimages) must be µ > 1500 and the time delay between the two counterimages (assuming intrinsic variability) would be very short, ΔT < 4 days, but these fluctuations would be observed at the same red dot but separated by regular intervals of ΔT days.

• System 115 (panel A of Fig. 4) shows two pairs of resolved feature on either side of the CC. (As in System 82, one pair is inside the yellow circle in panel A.) The situation is very similar to System 82 discussed above, but the System 115 caustic is significantly less powerful. Therefore the upper limits on the sizes are larger than the ones in System 82.

• System 116 (panel F of Fig. 4) was barely detected in JWST images. A long, elongated arc runs almost perpendicular to the CC, as expected for a strongly lensed galaxy in this position of the image plane. The geometric redshift is poorly constrained, but the photometric redshift is consistent with z ≈ 8. Like System 75, this arc is an excellent candidate to search for very high redshift EMO stars, although it is even more challenging than System 75 given the higher redshift and smaller ${\mu }_0^{+}$ and ${\mu }_0^{-}$. Detecting an EMO star in this arc will require a considerable boost from micro-lensing. No candidates were found in the current observations.

The model parameters discussed above imply a relation between the normalization constants T₀ and µ₀. In general, the time delay $$\text{Δ}T=\frac{D_s{D}_d}{D_{\text{ds}}}\frac{\left(1+z_s\right)}{c}\left(\frac{1}{2}{(\theta -\beta )}^2-\psi (\theta )\right),$$(6)

where the first and second terms are known as the geometric and Shapiro contributions to the time delay. For the popular singular isothermal sphere (SIS) model, the time delay is simply $\text{Δ}T\propto \left({\theta }_1^2-{\theta }_2^2\right)$, where θ₁ and θ₂ are the observed positions of the two counterimages. A simple rotation to make the deflection field parallel to the line connecting the two counterimages gives θ₁ = C + d_cc and θ₂ = C − d_cc, where C is the position of the CC. That makes ΔT ∝ d_cc. Because near a CC $\mu \propto d_{\text{cc}}^{-1}$, the SIS model gives $T_0\propto {\mu }_0^{-1}$.

The observed relation between T₀ and µ₀, shown in Fig. 5, follows the simple expectation but with some dispersion. The simple scaling T_o = 10/ < µ₀ > reproduces to first order the relation between the amplitudes T₀ and µ₀. This comes from the irregular and dynamically perturbed nature of MACS0416. There is an interesting trend with redshift. Arcs with lower red-shif tend to have larger time delays and smaller magnification factors. This is a consequence of low-z images appearing closer to the inner regions of the cluster. In these regions, the lensing potential changes more rapidly with distance. Because the magnification near CCs is inversely proportional to the gradient of the lensing potential, caustic-crossing arcs near the center of the cluster are generally less magnified. On the other hand, the larger potential near the center increases the Shapiro term in the time delay, resulting in counterimages from these galaxies having relatively longer time delays. The largest normalization factor T₀ is observed for System 1 at z = 0.9397. This system is very close to the main BCG of the cluster, where the Shapiro delay is largest.

Fig. 4 Critical curves from the lens model superposed on the MACS0416 color image. The image scale and orientation are indicated. The red curve is for sources at z = 6.145 (systems 73 and 74), the green curve is for z = 2.091 (System 12), and the magenta curve is for z = 0.9397 (System 1). The image colors have been stretched to better show the distribution of ICL. Panels A through G show enlarged images of areas indicated by white rectangles in the main image. These areas contain some of the most prominent arcs crossing a caustic, where candidate EMOs (marked by thick yellow circles) can be found in JWST images. The CCs in panels A–G are the same as in the larger plot. The numbers in each panel indicate the system ID, and white (rank A) or yellow (rank B) thin circles indicate system counterimages. The lensing properties of these systems are described in more detail in Section 5. Galaxies Warhol (System 1), Spock (System 2), and Mothra (System 12) are shown in panels C and D.

Fig. 5 Magnification versus time delay normalization factors for the 12 caustic-crossing galaxies in MACS0416. Each point represents one of the 12 galaxies color-coded by their redshift. Low-redshift arcs normally appear near the center of the cluster, which gives larger Shapiro delays and steeper potentials. That largely explains the larger delays for smaller magnification factors. The x-axis shows the average of the two normalization factors $\langle {\mu }_0\rangle =\left({\mu }_0^{+}+{\mu }_0^{-}\right)/2$ for each arc. The dashed line shows the expected scaling for an isothermal model.

6 The relation between mass and globular clusters

A scaling between the number of GCs and the mass of the halo hosting them has been established both from observations at lower redshifts and from N-body simulations (Blakeslee et al. 1997; Harris et al. 2017; Burkert & Forbes 2020; Valenzuela et al. 2021; Dornan & Harris 2023). This relation is a powerful proxy for halo masses when lensing constraints or spectroscopic measurements (dynamical mass) are not available. Another use is to set prior masses for cluster galaxies and DM halos before starting to optimize a lens model. The relation can also be used as a sanity check of lens models because the derived lensing mass for each halo is expected to correlate with the observed number of GCs within the halo.

The JWST data have the resolution and sensitivity needed to see many GCs in MACS0416. The smallest ones (in terms of mass and hence luminosity) will remain undetected, but there should still be a scaling between the observed number of GCs and the halo masses given by our lens model. Appendix A gives details about the definition of the GC sample, and Fig. 6 shows the GCs detected in the region constrained by multiple lensed galaxies. This region contains N_GC ≈ 7000 GCs. As expected, the GCs concentrate around the massive halos. Fig. 6 also shows a good correspondence between the distribution of GCs, the distribution of projected mass (contours), and the ICL. Fig. 7 shows, as one example, GCs detected around Halo 1, which is particularly challenging for GC identification because it corresponds to a massive and luminous galaxy (with increased photon noise at its center). The image also shows what appears to be a binary supermassive black hole (SMBH) in the center of Halo 1. The binary becomes evident only after high-pass filtering. The compact core surrounding the binary nucleus is intriguing. Free-floating stars in the central region are expected to be scattered away by the scouring action of the putative binary SMBH Postman et al. (2012). Additional examples of GCs detected in JWST images around other massive halos are shown in Figs. A.1 and A.2.

At a more quantitative level, Fig. 8 shows the integrated profiles for seven well-constrained halos (marked by big yellow circles in Fig. 6). At distances ≲30 kpc, the number density of GCs scales as the projected mass with a slope in the number (or mass) surface density α ≲ −0.3. At larger radii, the surface number density of GCs falls faster, with α ≈ −1.3, that is ≈1 dex faster than the mass surface density, in good agreement with results from simulations (Schaller et al. 2015; Pillepich et al. 2018; Alonso Asensio et al. 2020) as well as earlier observational results from JWST in galaxy-cluster lenses (Diego et al. 2023b). The number of GCs found within each of the seven correlates with the halo projected mass within the same radius. The halo mass is reasonably well approximated by the law M/M_⊙ ≈ 1.2 × 10¹⁰N_GC. (The circular regions where we computed the halo mass and N_GC are shown in Fig. 6.) We took smaller circles for smaller halos, but the circle size makes little difference provided the radius is in the range where the slopes of the solid-line and dashed-line profiles in Fig. 8 are similar.

7 Discussion

Despite the large number of available lensing constraints, the lens model still struggles to offer a satisfactory answer for some lensed galaxies. A remarkable example is the Spock arc (System 2), where the lens model is unable to place the CC at a viable symmetry point. In part this is because only two positions (marked by two blue ellipses in Fig. 9) are available as constraints. These ellipses mark the only easily recognizable pair of counterimages in this arc. JWST images show several unresolved features along the arc, but the correspondence (in terms of pairs of counterimages) between them is unclear. Applying a high-pass filter to the data reveals a possible symmetry plane close to the midpoint of the arc (long dashed line in Fig. 9). However, our lens model (solid white curve) does not place the CC close enough to that point. The same plot shows two alternative CCs derived by changing the grid configuration (one of the free variables in our reconstruction algorithm). The CC (white, red, and blue curves in Fig. 9) intersects the arc just once, but a loop-like structure forms at ≈0″.5 from the arc and south of the possible symmetry point.

A small modification in the lens model near the Spock arc could bring its loop-like CC to an intersection point with the arc at the presumed symmetry point. This would require a slight reduction in the mass in the diffuse component (parameterized by the grid of Gaussian points). This may be possible if the number of constraints nearby, or in the arc itself, is increased to force the CC to pass through the intersection between the larger dashed line and the arc. However, these constraints cannot be confirmed valid at the moment, so they are not used. Reducing the mass of the smooth component would also close the loop-like CC around the cluster galaxy at the bottom of Fig. 9. This naturally produces a CC that grazes the SE portion of the arc and runs between the triplet marked by three green circles. This triplet’s image system would then include the counterimage in the fourth green circle in the NW portion of the arc. What would have been a single counterimage in the SE gets split into three counterimages because of the cluster galaxy to the south. In this interpretation, where the true CCs pass through the dashed lines, only two of the three white circles could be counterimages of a single source, and the third would be a different source seen in projection. The object marked with a blue circle and a question mark has no obvious counterpart. It could also be a projection effect (one of the many GCs) or a small structure (possibly an EMO star) being magnified locally by a microlens or millilens. Additional lens models based on the available and future spectroscopic confirmation of some of the lensed galaxy candidates will provide further evidence in support of the single versus multiple caustic crossings of this arc.

In conjunction with the lens model, the number of transient events reported by Yan et al. (2023) constrains the massive-star population in the Spock galaxy. The observed number agrees with the low-end expectation of Diego et al. (2024b), corresponding to a situation where SG stars above the Humphreys-Davidson (HD, Humphreys & Davidson 1979) limit do not exist. The HD limit is an observational upper limit for the luminosity of red SG stars derived in the local Universe. In other words, the results from Yan et al. (2023) are consistent with a scenario where the HD limit is also valid at z = 1. This limit may be a consequence of wind-driven mass loss in red SG stars (Vink & Sabhahit 2023), a possibility that can be tested with observations of low-metallicity EMO stars at high-z.

The amplitudes listed in Table 1 provide context to past transient events observed in arcs such as Spock, Warhol, and Mothra as well as for future transient events in these arcs. JWST observes an average of ≈3 events per epoch at a depth of AB 29 magnitude in F200W (or ≈1 hour of integration time, Yan et al. 2023). One of the most intriguing events found by Yan et al. (2023) was D21_W1 in the Warhol galaxy (panel C in Fig. 4). Both D21_W1 and its counterimage, D21_W1′, were detected in previous HST images from the HFF program as well as in the new JWST images. Given the moderate macromodel magnification factor predicted from our lens model, µ⁻ = 28 and µ⁺ = 46 for D21_W1 and D21_W1′, respectively, the source’s intrinsic luminosity must be significantly larger than the luminosity of SG stars and more consistent with a GC or small star-forming region. If this source is indeed a GC at z ≈ 1 (or similar object containing multiple stars), that would explain its brightness in previous HST images and during the JWST observations (in epoch 1 and before the transient was observed). However, this would make it difficult to explain the flux increase of ≈0.5 magnitude in F200W between epoch 1 and epoch 2 of PEARLS, separated by 80 days).

A single star in a GC being microlensed results in a relatively small change in the observed flux because the increase in flux from the microlensed star still represents a relatively small fraction of the total flux from the entire GC (Dai 2021). This is a similar argument to the one that applies to larger sources, with radii much larger than the Einstein radius of a microlens (typically much smaller than 1 parsec) and for which microlensing effects are very small. Thus any source that shows significant change due to microlensing must be much smaller than 1 parsec. If the source is a single star, it must have an extraordinary brightness. At the redshift of Warhol, the distance modulus is 43.96. At macromodel magnification µ = 30, the source must have absolute magnitude ≈−12 in order to appear as AB ≈ 28 magnitude in JWST F200W. This would make D21_W1 comparable to the brightest stars known in our local Universe, such as V4650 Sagittarii, a luminous blue variable (LBV). However, D21_W1 exhibits a SED more compatible with a cooler star (typically fainter), making it difficult to explain its luminosity. A long-lasting (years) outburst episode in a cooler SG star could explain the high luminosity. This would be similar to the case of Godzilla, another EMO star at z = 2.38, that is suspected to be undergoing a major outburst lasting several years, similar to the Great Eruption of Eta Carinae in the 19^th century (Diego et al. 2022).

Considering both luminosity and light curve, the transient event D21_W1 (Yan et al. 2023) can be explained in two alternative ways: 1) a secondary minor outburst within a major stellar-outburst event, or 2) a microlensing event of a single luminous star within a group of stars. In the first scenario, the outburst should be observed in both counterimages because it is intrinsic to the source. The counterimage, D21_W1′, showed no variability in JWST data. However, this is consistent with our model prediction of ΔT = 1.35 years (d_cc = 0″.743) while the JWST observations span approximately 5 months. According to our model, if an outburst took place D21_W1 during the second epoch of JWST observations, it would have appeared in D21_W1′ ≈ 1.15 years before JWST first pointed to this cluster. In this case, the fact that no variability was observed in the counterimage during the JWST observations implies that the minor outburst lasted less than ≈0.5 years in the rest frame of the source. The explanation of the second scenario above is much simpler. If a microlens temporarily boosted the magnification of D21_W1, we expect no correlation between flux changes in the two images, and D21_W1′ may maintain constant flux for long periods of time. Future observations spanning a period longer than ΔT are needed to establish correlations between the two counterimages’ variability. Also, as discussed earlier, if the microlensing interpretation is correct, the source must be a small group of stars with the star undergoing microlensing being one of the dominant stars in that group, so its relative change in flux due to microlensing is significant when compared to the total flux of the group. The asymmetry in the light curves (Yan et al. 2023) favors the outburst scenario because microlensing events tend to produce symmetric light curves. However, the evidence in favor of the outburst scenario is not conclusive, and more data in future campaigns are needed.

Another interesting issue is the validity of the predicted correlation between the number of GCs and the lensing mass. Before JWST, GCs in galaxy clusters could be observed only in clusters at relatively low redshifts (z ≲ 0.2). At these redshifts, the optical depth for lensing is also very small, so reliable lensing masses are not available for such low-z halos. The notable exception is Abell 1689, for which HST was able to observe both GCs in the cluster and a relatively large number (30) of lensed galaxies Broadhurst et al. (2005); Alamo-Martínez et al. (2013). JWST can extend this type of study to larger redshifts and simultaneously detect hundreds of lensed galaxies behind massive clusters and thousands of GCs in the lensing cluster. As described in Section 6, there is a linear correlation between halo mass and number of GCs. This relation is similar to the one obtained from observations, M_vir = 5 × 10⁹ × N_GC M_⊙ (Burkert & Forbes 2020) but with an amplitude a factor 2.4 times higher. The larger amplitude can be interpreted as a result of the fraction of fainter GCs that are missed in our sample due to the relatively high redshift of the cluster and the limited depth of the observations. This would imply that the total number of GCs in this cluster is ≈2.4 times higher than the observed number. This estimate is relevant for discussions of the lensing effect of GCs near the CC region because one can estimate the true number density of GCs near a certain region by simply counting the number of GCs in that region and correcting for a factor 2.4. The actual correction factor depends on the redshift of the lens and depth of the observations. Thus one can assess the probability that a GC is aligned with a particular feature in a lensed arc. For instance, Mothra was interpreted (Diego et al. 2023c) as a star at z = 2.091 being magnified by a millilens in the cluster. The most likely candidate for the millilens is one of the multiple GCs around the NE BCG. The probability of a chance alignment of a GC with the background star was estimated as 8% based on the observed local number density of GCs near the Mothra position (≈1.6 GCs per arcsecond²). Correcting for the factor of 2.4, this probability increases to ≈20%, making a GC an even more likely interpretation for the source millilensing Mothra.

The logic for the abundance of GCs can be applied to other arcs. For instance, the Spock arc has a local number density of GCs similar to that of Mothra, so the probability of finding a GC millilens in a given position in that arc should also be ≈20%. For Warhol, the probability should be a bit higher because the local number density of observed GCs is also higher because Warhol is closer to the BCG.

The linear scaling between the number of GCs and halo mass is not maintained at larger radii. This is because the number density of GCs falls off more sharply beyond ≈30 kpc, and their density cannot keep up with the DM density, which falls more slowly with distance. A similar trend was observed in Abell 1689 Alamo-Martínez et al. (2013). Inside the regions shown in Fig. 6, M = 2.7 × 10¹⁴ M_⊙, and N_GC = 6500, while the M–N_GC law above would predict M = 7.8 × 10¹³ M_⊙, that is ≈3.4 times less mass. The scaling between M and N_GC is thus a good proxy for halo mass only up to a few tens of kpc from the center of the halo. As found previously, the DM halo for MACS0416 is shallower than the profile of the distribution of GCs. This is fully consistent with the CDM model, where DM is composed of much lighter particles than the baryonic structures they interact with (gravitationally) such as stars, GCs, or even member galaxies, hence acquiring larger velocities in the constant 3-body interactions. Another way to look at the same process is that dynamical friction causes the most massive objects, such as GCs, to migrate toward the center of the halo.

Fig. 6 Globular clusters in MACS0416 (green dots) superposed on the JWST image. The seven yellow circles mark the positions of massive halos for which the lens model places good constraints on the halo masses. The size of each circle indicates the radius within which the mass was computed. The magenta contours are the convergence map, κ, from the lens model computed at z = 3 and in intervals of 0.1, starting at κ = 0.5.

Fig. 7 Central region of Halo 1 after high-pass spatial filtering. The scale and orientation are indicated. Green circles mark sources identified as GCs, and yellow circles mark the components of an apparent binary AGN. The two candidate SMBHs are separated by less than 1 kpc in projection. The image was created by applying a 3-pixel Gaussian high-pass filter to three JWST images as marked at the top of the figure before combining them into the color composite.

Fig. 8 Profiles of the lensing mass, number of GCs, and light around seven halos. The solid lines show the projected mass profiles from the lens model as a function of radius, the dashed lines show the corresponding number-density profile of GCs, and the dotted lines show the F200W surface brightness. The profiles are centered on each of the seven halos identified in Fig. 6 and plotted in different colors as shown in the legend. The inset shows the total mass versus number of GCs for each halo integrated up to the radius of the halo’s yellow circle in Fig. 6. The solid line in the inset shows the relation M = 1.2 × 1010NGC.

Fig. 9 High-pass filtered image of the Spock arc with CCs and the possible interpretation of multiply imaged sources. The white curve is the CC predicted by the full model (systems with rank A+B and adaptive grid with 495 grid points) at z = 1.0054. The image shows multiple unresolved features along the Spock arc. Sources marked with the same color circle are assumed to be the same source that is multiply imaged, and the dashed lines mark the positions of the CC based on these sources. The white circles mark a lensed galaxy candidate at a slightly different redshift. If the mass of the smooth component is reduced slightly in this portion of the lens plane, the main CC would move toward the NW up to the position of the long dashed line (≈0″.5 away). The critical line around the member galaxy in the south of the image would then close (smaller dashed line), naturally producing the triplet marked with three green circles. The red and blue CCs are from two alternative models where the grid configuration is changed allowing a slightly different distribution for the smooth mass component. Because we did not use multiple knots in this arc as constraints, the lensing constraints are insufficient to properly reproduce the multiple crossing positions of the CC.

8 Conclusions

The new data from JWST reveal a wealth of new lensed galaxy candidates behind MACS0416. The data bring the total number of lensed background galaxies to 119 galaxies (including candidates, a few of which may turn out to be false positives) producing a record number of 343 multiple images. With such a large number of constraints, free-form algorithms can produce reliable lens models while making a minimum number of assumptions regarding the distribution of the invisible dark matter. The free-form WSLAP+ model presented here was based on all the spectroscopic systems and the full set of constraints that includes also system candidates without spectroscopic confirmation. The redshifts for the systems without spectroscopic redshift were derived from the lens model (geo-z), which shows a good correlation with independent photo-z estimates. The model also provides magnification and time-delay estimates for all lensed images. The derived lens model consists of two large halos with nearly identical mass profiles, in agreement with earlier results derived using parametric models and a similar set of spectroscopic systems.

MACS0416 is exceptional in its number of caustic-crossing galaxies. We modeled the magnification and time delay for 12 of them, paying special attention to possible EMO stars in these arcs. The amplitudes of the time-delay and magnification scaling relations are anticorrelated, as expected for smooth analytical models, but there are departures from perfect anti-correlation. These can be attributed to galactic-scale contributions to the deflection angle and the non-relaxed nature of the cluster.

For the Spock arc, the large-scale lens model fails to produce a CC at a viable symmetry point. Small changes to the lens model can make the CC intersect the arc at three points, thus explaining unresolved features seen in the new JWST data around these points. The candidate EMO star, D21_W1, is an intriguing and relatively bright transient found by Yan et al. (2023) in the Warhol arc at z = 0.9397. Based on the model’s time-delay prediction between this source and its counterimage D21_W1′, the source may be better described by an intrinsically variable source than by a microlensing event. Spectral information rules out a hot LBV. A cooler variable star that temporarily increased in brightness may offer a better explanation, but further monitoring is needed to rule out micro-lensing.

For the most massive halos in MACS0416, there is a simple linear correlation between the lensing mass and number of GCs. This is the same as found at lower redshifts although the normalization is higher (i.e., more mass for a given number of GCs). This is likely related to the incompleteness of our sample of GCs. The correlation is valid up to a radius of a few tens of kpc but not beyond.

Data availability

Table B.1 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/690/A114

Acknowledgements

We dedicate this paper to the memory of our dear PEARLS colleague Mario Nonino, who was a gifted and dedicated scientist, and a generous person. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with JWST programs 1176 and 2738. We thank the anonymous referee for useful comments and suggestions. JMD acknowledges the support of project PID2022-138896NB-C51 (MCIU/AEI/MINECO/FEDER, UE) Ministerio de Ciencia, Investigación y Universidades. RAW, SHC, and RAJ acknowledge support from NASA JWST Interdisciplinary Scientist grants NAG5-12460, NNX14AN10G and 80NSSC18K0200 from GSFC. Work by CJC acknowledges support from the European Research Council (ERC) Advanced Investigator Grant EPOCHS (788113). BLF thanks the Berkeley Center for Theoretical Physics for their hospitality during the writing of this paper. MAM acknowledges the support of a National Research Council of Canada Plaskett Fellowship, and the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE17010001. CNAW acknowledges funding from the JWST/NIRCam contract NASS-0215 to the University of Arizona. We thank the CANUCS team for generously providing early access to their proprietary data of MACS0416. 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); by the Ministry of Science & Technology, Israel; and by the Israel Science Foundation Grant No. 864/23. We also acknowledge the indigenous peoples of Arizona, including the Akimel O’odham (Pima) and Pee Posh (Maricopa) Indian Communities, whose care and keeping of the land has enabled us to be at ASU’s Tempe campus in the Salt River Valley, where much of our work was conducted. SOFTWARE: Astropy: http://www.astropy.org (Astropy Collaboration 2013; Astropy Collaboration 2018); IDL Astronomy Library: https://idlastro.gsfc.nasa.gov (Landsman 1993); Photutils: https://photutils.readthedocs.io/en/stable/ (Bradley et al. 2020); ProFound: https://github.com/asgr/ProFound (Robotham et al. 2018); ProFit: https://github.com/ICRAR/ProFit (Robotham et al. 2017); SExtractor: SourceExtractor: https://www.astromatic.net/software/sextractor/orhttps://sextractor.readthedocs.io/en/latest/ (Bertin & Arnouts 1996). FACILITIES: Hubble and James Webb Space Telescope Mikulski Archive https://archive.stsci.edu

Appendix A Globular clusters

Fig. A.1 Detail of globular clusters found in MACS0416 (marked with green circles) detected around the NE BCG. The galaxy Warhol is shown in the top right. Unresolved features in Warhol and other resolved galaxies are not included in the sample of GCs. The image is a composition made from two filtered images (using a Gaussian Kernel with 3 and 15 pixels) and the raw data. The red channel is the F150W filter, the blue channel is the F115W. The most sensitive F200W is used for the green channel in order to maximize the effectiveness of the visual inspection.

To search for GCs, we first high-pass filtered the three bands F115W, F150W, and F200W. Bands in the LW channel lack the resolution for GCs while the F090W is less sensitive than F115W. The high-pass filter is defined as the difference between the direct image (PSF convolved) and a filtered image with a Gaussian kernel with 3 pixels FWHM. This results in a filter very similar to an A Trous wavelet. After filtering the three bands, we visually search for unresolved sources in the combination color image. Since visual inspection is most sensitive to green-yellow colors, we set the filtered F200W as the green channel, since it offers the best combination between resolution and sensitivity to detected GCs. Blue and red channels are set to F115W and F150W, respectively. Most GCs in our sample are seen in all three bands, but some are seen in just two. We reject candidates that are only seen in one band. Some of these single-band candidates can be noisy fluctuations but also faint GCs or very small high-redshift dropout galaxies (see Berkheimer et al. 2024, for a discussion of contaminants in JWST images).

Detection of GCs is more challenging near crowded areas. For instance, around the BCGs, the larger fluctuations from photon noise make it harder to detect point sources. Detection of the same source in at least two bands is used to reduce false positives. The presence of background lensed galaxies adds another complication since these galaxies may contain unresolved features on their own that can be confused with GCs. This is shown in Fig. A.1, where the lensed galaxy Warhol overlaps with a region containing a relatively high number density of GCs near the NE BCG. We avoid these regions and as a consequence we detect fewer GCs near background galaxies.

A different example is shown in Fig. A.2, this time around the compact halo (labeled CH in Fig. 6 and with steep profiles for the ICL and GCs, as shown in Fig. 8). This halo is not affected by foreground lensed galaxies and due to its compact nature we can detect GCs much closer to the center of the halo.

Fig. A.2 Detail of the rich GC structure found around the compact halo at RA=4:16:08.84, DEC=-24:04:20.12 (CH in Fig. 6). The image has been processed in a similar way to Fig. A.1.

Appendix B Lensing constraints

The lensing constraints used in this work are summarized in Table B.1.¹ The first 77 systems in our sample are the spectroscopic galaxies and are the same as the systems used in Bergamini et al. (2023); Diego et al. (2023b). System 77 has the same redshift as System 33 forming a pair of lensed galaxies. Systems with ID 78 and above are candidate systems. Most of them are found in the new JWST data but some where also proposed as system candidates from previous HST data (Jauzac et al. 2014, 2015; Johnson et al. 2014; Diego et al. 2015; Kawamata et al. 2016). Arcs in the spectroscopic sample that are marked with † are new counterimages found in JWST data. Columns four, five, and six show the spectroscopic redshift, geometric red-shift, and photometric redshift, respectively for each system or arc. Images with negative photometric redshift are too faint to produce meaningful redshift estimates. The geometric redshift for systems 113, 114, 116, and 117 are all above 10 and unreliable. These systems are either poorly fitted by the spectroscopic lens model (systems 113 and 114) or the geometric redshift prediction is also poor because the system has only two counterimages, with a relatively small angular separation between counterimages (systems 116 and 117). However, for System 117, both the photo-z and geo-z predict this galaxy to be at high redshift. Some of the photo-z estimates may be affected by contamination from the member galaxies (z ≈ 0.4) resulting in photo-z estimates close to the cluster redshift. Column seven shows the magnification predicted by the full model and for each arc at the RA and DEC positions listed in columns 2 and 3, while column eight shows the time delay (at the same RA and DEC positions) relative to the first image to arrive from the same model. Finally, the last column indicates the rank of the images. Rank A is for images that are spectroscopically confirmed. Rank B are for images which are good candidates but lack spectroscopic confirmation (color match, morphology and parity, and good consistency with the spectroscopic lens model). Rank C are similar to Rank B images but usually lack morphological confirmation because they are unresolved. Finally images with rank D (two systems) are the least reliable because they are in tension with the spectroscopic lens model. System 113 (Rank D) is an interesting candidate for follow up since it could consist of a background galaxy being lensed by an AGN or QSO at RA=4:16:08.5462, DEC=−24:05:22.004. It is not used as a lensing constraint but is included here as a possible interesting candidate. System 114 (D) appears lensed by a member galaxy but in a region of the cluster that is poorly constrained. The small separation between the counterim-ages of System 114 make it a relatively poor lensing constraint. Other interesting candidates, not included in the list are a small Einstein ring of a potentially high-redshift background galaxy lensed by a galaxy at RA=4:16:11.1176, DEC=-24:05:24.916. Several galaxy-galaxy lenses are also easily identifiable in the new JWST data and beyond the region constrained by the multiply lensed galaxies. Errors for geo-z are derived from the lens model that relies in rank A systems only. Errors in the time delay and magnification are not included because they are derived from the lens model including rank B systems, which lack spectro-scopic confirmation. Hence, these magnification and time delays should be taken only as indicative values.

Fig. B.1 A posteriori confirmation of zɡeo with the new MUSE redshifts in Rihtaršič et al. (2024). Uncertainties have been multiplied by a factor 4 to better show some of the error bars. The bottom panel shows (zgeo − Zspec)/zspec. The two dotted lines in the bottom panel show the 5% relative change.

During the last phases of the refereeing process, new spectroscopic redshifts for 13 of our system candidates were obtained by MUSE in Rihtaršič et al. (2024). The new redshifts confirm these systems but more importantly demonstrates the accuracy of the geometric redshifts obtained from the well calibrated lens model. The comparison between our z_geo and the new MUSE z_spec is shown in Figure B.1. The 13 systems are 82, 87, 89, 90, 91, 95, 97, 98, 99, 100, 106, 109, and 115 in Table B.1 below. The bottom panel of the figure shows the relative difference between z_geo and z_spec. Most systems fall within 5% and five of them are within 2% of the correct value. The worst outlier is 15% above the correct value. The two high-z systems (97 and 100) have z_geo within 12% and 8% of the correct value, respectively.

References

All Tables

All Figures

Fig. 1 Color image of MACS0416. The scale and orientation are indicated. White circles mark previously known counterimages with spectroscopic redshifts from Richard et al. (2021); Bergamini et al. (2023). Yellow circles mark reliable (rank B) HST- and JWST-identified counterimage candidates that lack spectroscopic confirmation. Magenta and red circles mark the least reliable counterimage candidates with rank C and D respectively.

In the text

Fig. 2 Comparison of photometric redshifts with spectroscopic (blue symbols) or geometric (red symbols) redshifts for the lensed systems in Table B.1. For the photometric redshifts, we chose the value (from the multiple counterimages of the same system) closest to the spect-z or geo-z. The dotted line shows equality. The comparison between spect-z and geo-z is not shown but it follows almost perfectly the dotted line (see Table B.1).

In the text

Fig. 3 Mass profile of the lens models. Blue curves correspond to profiles centered on the NE BCG while red curves are for profiles centered at the SW BCG. Solid lines represent the models derived using only spectroscopically confirmed (rank A in Table B.1) galaxies as constraints and an adaptive grid with 604 grid points for the smooth DM. Dotted lines show the corresponding solution for an alternative model with the same constraints (rank A) but with only 495 grid points. Dashed lines represent the same 495-point grid but with the full sample of lensed systems ranked A and B (geometric redshifts) as constraints. The mass profiles are given in dimensionless κ units assuming a source at z = 3. The long-dashed line shows a gNFW profile with parameters (Eq. (3)) γ = 0.9, α = 2, β = 3, and scale radius rs = 640 kpc. The two vertical black segments indicate the distance interval where the profiles are constrained by lensing. The inset shows the integrated mass as a function of radius in units of 1014 M⊙ for the same three profiles shown in the main panel. The black long-dashed line is the integrated profile when the center is chosen as the midpoint between the two BCGs (RA = 4hr16m08s.428, Dec = –24°04′21″.0).

In the text

Fig. 4 Critical curves from the lens model superposed on the MACS0416 color image. The image scale and orientation are indicated. The red curve is for sources at z = 6.145 (systems 73 and 74), the green curve is for z = 2.091 (System 12), and the magenta curve is for z = 0.9397 (System 1). The image colors have been stretched to better show the distribution of ICL. Panels A through G show enlarged images of areas indicated by white rectangles in the main image. These areas contain some of the most prominent arcs crossing a caustic, where candidate EMOs (marked by thick yellow circles) can be found in JWST images. The CCs in panels A–G are the same as in the larger plot. The numbers in each panel indicate the system ID, and white (rank A) or yellow (rank B) thin circles indicate system counterimages. The lensing properties of these systems are described in more detail in Section 5. Galaxies Warhol (System 1), Spock (System 2), and Mothra (System 12) are shown in panels C and D.

In the text

Fig. 5 Magnification versus time delay normalization factors for the 12 caustic-crossing galaxies in MACS0416. Each point represents one of the 12 galaxies color-coded by their redshift. Low-redshift arcs normally appear near the center of the cluster, which gives larger Shapiro delays and steeper potentials. That largely explains the larger delays for smaller magnification factors. The x-axis shows the average of the two normalization factors $\langle {\mu }_0\rangle =\left({\mu }_0^{+}+{\mu }_0^{-}\right)/2$ for each arc. The dashed line shows the expected scaling for an isothermal model.

In the text

Fig. 6 Globular clusters in MACS0416 (green dots) superposed on the JWST image. The seven yellow circles mark the positions of massive halos for which the lens model places good constraints on the halo masses. The size of each circle indicates the radius within which the mass was computed. The magenta contours are the convergence map, κ, from the lens model computed at z = 3 and in intervals of 0.1, starting at κ = 0.5.

In the text

Fig. 7 Central region of Halo 1 after high-pass spatial filtering. The scale and orientation are indicated. Green circles mark sources identified as GCs, and yellow circles mark the components of an apparent binary AGN. The two candidate SMBHs are separated by less than 1 kpc in projection. The image was created by applying a 3-pixel Gaussian high-pass filter to three JWST images as marked at the top of the figure before combining them into the color composite.

In the text

Fig. 8 Profiles of the lensing mass, number of GCs, and light around seven halos. The solid lines show the projected mass profiles from the lens model as a function of radius, the dashed lines show the corresponding number-density profile of GCs, and the dotted lines show the F200W surface brightness. The profiles are centered on each of the seven halos identified in Fig. 6 and plotted in different colors as shown in the legend. The inset shows the total mass versus number of GCs for each halo integrated up to the radius of the halo’s yellow circle in Fig. 6. The solid line in the inset shows the relation M = 1.2 × 1010NGC.

In the text

Fig. 9 High-pass filtered image of the Spock arc with CCs and the possible interpretation of multiply imaged sources. The white curve is the CC predicted by the full model (systems with rank A+B and adaptive grid with 495 grid points) at z = 1.0054. The image shows multiple unresolved features along the Spock arc. Sources marked with the same color circle are assumed to be the same source that is multiply imaged, and the dashed lines mark the positions of the CC based on these sources. The white circles mark a lensed galaxy candidate at a slightly different redshift. If the mass of the smooth component is reduced slightly in this portion of the lens plane, the main CC would move toward the NW up to the position of the long dashed line (≈0″.5 away). The critical line around the member galaxy in the south of the image would then close (smaller dashed line), naturally producing the triplet marked with three green circles. The red and blue CCs are from two alternative models where the grid configuration is changed allowing a slightly different distribution for the smooth mass component. Because we did not use multiple knots in this arc as constraints, the lensing constraints are insufficient to properly reproduce the multiple crossing positions of the CC.

In the text

Fig. A.1 Detail of globular clusters found in MACS0416 (marked with green circles) detected around the NE BCG. The galaxy Warhol is shown in the top right. Unresolved features in Warhol and other resolved galaxies are not included in the sample of GCs. The image is a composition made from two filtered images (using a Gaussian Kernel with 3 and 15 pixels) and the raw data. The red channel is the F150W filter, the blue channel is the F115W. The most sensitive F200W is used for the green channel in order to maximize the effectiveness of the visual inspection.

In the text

	Fig. A.2 Detail of the rich GC structure found around the compact halo at RA=4:16:08.84, DEC=-24:04:20.12 (CH in Fig. 6). The image has been processed in a similar way to Fig. A.1.
In the text

	Fig. B.1 A posteriori confirmation of zɡeo with the new MUSE redshifts in Rihtaršič et al. (2024). Uncertainties have been multiplied by a factor 4 to better show some of the error bars. The bottom panel shows (zgeo − Zspec)/zspec. The two dotted lines in the bottom panel show the 5% relative change.
In the text