© The Authors 2023

1. Introduction

ACT-CL J0102-4915, known as El Gordo and arguably the most famous galaxy cluster at redshift z > 0.8, has been observed with the James Webb Space Telescope (JWST) and its Near Infrared Camera (NIRCam) as part of the Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS) Guaranteed Time Observing (GTO) program (Windhorst et al. 2023). At z = 0.870, El Gordo is a merging cluster with a double-peaked galaxy distribution (Williamson et al. 2011; Menanteau et al. 2012). In X-rays, El Gordo exhibits a cometary structure clearly visible in Chandra data. This suggests that El Gordo is being observed right after a collision of two subgroups (Molnar & Broadhurst 2015; Zhang et al. 2015), similar to the iconic Bullet cluster. The presence of two radio relics ahead of and behind the X-ray cometary structure supports this interpretation (Molnar & Broadhurst 2015; Lindner et al. 2014). Ng et al. (2015) argue that El Gordo is in a return phase after first core passage. That is, the cluster is being observed after the phase of maximum separation, and the two groups are moving back toward each other. This interpretation is, however, challenged by lens models based on strong lensing, which place most of the mass in the southeastern group (Zitrin et al. 2013; Cerny et al. 2018; Diego et al. 2020). The massive nature of El Gordo is already demonstrated by its strong Sunyaev-Zeldovich (SZ) effect. Despite its relatively high redshift, El Gordo shows up in Planck maps as a massive cluster with a signal-to-noise ratio of ≈13 (Planck Collaboration XXVII 2016). Some of the SZ signal may be due to increased pressure from the ongoing merger.

El Gordo’s mass has been estimated through a variety of techniques, as shown in Table 1. These methods are in agreement and indicate that El Gordo is the most massive cluster in its redshift range with an estimated mass above 10¹⁵ M_⊙. The most recent measurement (Caminha et al. 2022) is notable because its strong-lensing model was based on a set of 23 spectroscopically confirmed systems. Prior to that work, all lens models had relied on photometric redshifts.

At the upper end, the derived cluster masses are in tension with the standard ΛCDM model. Jee et al. (2014), for example, found a total mass ≈2.8 × 10¹⁵ M_⊙, exceeding the maximum mass ≲2 × 10¹⁵ M_⊙ (Harrison & Coles 2012) expected at this redshift. A similar conclusion was reached from large N-body simulations. Using the very large 630 Gpc³N-body simulation Jubilee (based on a standard ΛCDM model), Watson et al. (2014, their Fig. 5) found that the most massive cluster at z = 0.9 is expected to have M_178 m ≈ 1.7 × 10¹⁵ M_⊙. (Symbols for M are defined in a footnote to Table 1.) The observed masses are mostly above this maximum mass, but they can become smaller if one considers the Eddington bias (Waizmann et al. 2012). Despite that, El Gordo is still uncomfortably close to the expected ΛCDM limit. El Gordo was found in a survey – the original Atacama Cosmology Telescope (ACT) survey – that covered less than 2% of the area of the sky. The low likelihood of finding such an extreme object in such a small area raises the tension with the ΛCDM model. Asencio et al. (2021) used the same Jubilee simulation to argue that, for El Gordo, the tension is high. On the other hand, after correcting for redshift and the Eddington bias, Waizmann et al. (2012) found El Gordo not to be in tension with ΛCDM. However, their analysis assumed the survey area to be 3.7 times higher than the area where El Gordo was originally found. Any tension can be reduced if previous mass estimates are too high. Additional mass estimates, based on alternative methods, can explore the uncertainty as to the mass of the cluster.

The X-ray emission exhibits an interesting offset between the peak of the X-ray emission and the position of the brightest cluster galaxy (BCG). Contrary to what happens in the Bullet cluster, the X-ray peak seems to be ahead of the BCG. However, in the interpretation of Ng et al. (2015), the BCG would be moving toward the second group, so the X-ray peak would be trailing the BCG. The returning-phase interpretation of El Gordo is challenged by dedicated N-body and hydrodynamical simulations reproducing most of the observations of El Gordo (Molnar & Broadhurst 2015; Zhang et al. 2015). Molnar & Broadhurst (2018) demonstrate that the speed of the outgoing shocks can be very large (4000–5000 km s⁻¹) in a massive, merging cluster, such as El Gordo, therefore leaving the system before the first turnaround.

El Gordo is also unique in that it is a powerful lens at a relatively high redshift. One of the features that makes El Gordo an attractive target for lensing studies is the fact that for sources at a high redshift, critical curves form at relatively large distances from the member galaxies. This is particularly true in the gap between the two clusters, where the critical curves are relatively undisturbed by nearby member galaxies. Having undisturbed critical curves is relevant to observe caustic crossing events of distant stars (Kelly et al. 2018; Diego et al. 2018) because, in this case, the maximum magnification can be larger than in situations where critical curves are affected by microlenses in member galaxies or the intracluster medium. Caustic crossing events have been proposed as a technique useful to study Population III stars and stellar-mass black-hole accretion disks with JWST (Windhorst et al. 2018).

Because El Gordo is the highest-known-redshift cluster with potentially significant transverse motion—based on the X-ray morphology and the two lensing-mass centers discussed in this paper–it is an ideal target for JWST follow-up to search for caustic transits at z ≫ 1 and possibly for caustic transits at z > 7. For this reason, El Gordo was selected as a JWST GTO target. This paper uses the new JWST data to derive the mass distribution. We used our free-form lensing reconstruction code WSLAP+ (Diego et al. 2005, 2007, 2016, 2022; Sendra et al. 2014), which does not rely on assumptions about the distribution of dark matter. Our results offer an important cross-check with previous results because any disagreement between our free-form method and the results obtained by previous parametric methods could signal potential systematic problems in one (or both) types of modeling.

This paper is organized as follows. Section 2 describes the data and simulations used in this work. Section 3 describes the lensing constraints, and Sect. 4 explains the lens-modeling method and gives results for three models. The integrated mass and estimated virial masses for all three models are discussed in Sect. 5 together with a comparison with earlier models in the literature. Section 6 discusses small-mass perturbers found near a giant arc nicknamed “La Flaca.” Section 7 presents candidates to be compact, luminous objects near caustics and introduces “Quyllur”, likely a red supergiant star that is being magnified by a factor of several thousand. Section 8 discusses the key results, and Sect. 9 summarizes our conclusions. We adopt a standard flat ΛCDM cosmological model with Ω_m = 0.3 and h = 0.7. At the redshift of the lens, this cosmological model implies that 1″ corresponds to 7.8 kpc.

2. JWST observations and ancillary data

2.1. JWST data

El Gordo was observed by JWST on 2022 July 29 as part of the PEARLS Cycle 1 GTO program (pid #1176, P.I. R. Windhorst). Data were obtained in the IR filters F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W with effective exposure times of 1889 s for F150W and F356W, 2104 s for F200W and F277W, and 2490 s for F090W, F115W, F410M, and F444W. The 5σ magnitude limit is ≈29.9 in the long-wavelength filters. The best spatial resolution was obtained in the F150W filter with a measured Full-Width-Half-Maximum (FWHM) of 0063. Further details are in the PEARLS overview paper (Windhorst et al. 2023). Figure 1 shows an overview of the field.

Fig. 1. Composite image of El Gordo obtained after combining HST and JWST bands (filters spanning ≈0.5 μm to ≈5 μm). Top panel: White circles mark the systems confirmed spectroscopically by MUSE. Bottom panel: Different color rendering that highlights fainter objects. Yellow circles mark new system candidates identified with the new JWST data. In each label, numbers identify multiple images of the same background source, and letters identify the individual images. The image orientation and scale are indicated in the lower left.

JWST reveals details of this cluster with unprecedented clarity. Galaxies barely detected or not detected at all in Hubble Space Telescope (HST) data are easily detected by JWST, and many are much brighter at wavelengths longer than 2 μm than at the longest HST wavelength (1.6 μm). Figure 2 shows one such galaxy, which is heavily distorted by cluster-galaxy lensing into a fishhook shape. We refer to this galaxy as “El Anzuelo” (Spanish for fishhook). Pairs of multiply lensed features are easily identified in the image. Image pairs constrain the position of critical curves, which are expected to pass between the images. The red nucleus of the extended background galaxy is lensed into at least three images clearly visible in the reddest JWST bands. The three images of the nucleus coincide with the position of Atacama Large Millimeter/submillimeter Array (ALMA) sources EG-SMG 2 and EG-SMG 4 (Cheng et al. 2023), the latter containing two images barely resolved by ALMA. Previous HST data did not show El Anzuelo as multiply lensed because this galaxy is red and therefore faint at wavelengths shorter than 1.6 μm. HST data do, however, identify a multiply lensed blue feature.

Fig. 2. El Anzuelo galaxy in the northwestern part of the cluster. The galaxy has a photo-z of 3.58 (Cheng et al. 2023). The color-coded circles mark image pairs of three lensed sources. Each pair identifies a critical point between the images.

Another example illustrating the JWST data quality is shown in Fig. 3, which shows the giant arc already observed in pre-JWST data. We refer to this galaxy as “La Flac” (“the thin one” in Spanish). Unlike El Anzuelo, La Flaca is easily visible in HST data. The new JWST data show with greater detail the multiply lensed features within this arc. In particular, the nucleus of the arc appears three times (2.1a, 2.2a, and 2.4a) instead of the expected two times. The extra image (2.4a, already seen in HST data) is in fact a double image (unresolved) that appears as a consequence of a small perturber marked with a white arrow in the figure. Another perturber splits counterimage 2.2b into three images (2.2b, 2.4b, and 2.5b). Section 6 discusses these perturbers in greater detail. Unlike El Anzuelo, La Flaca does not show any clear pair of counterimages that can be used to constrain the position of the critical curve. However, it is still possible to get an approximate location based on simple arguments. One can use the ratio of distances between knots 2.1a–2.1e and 2.2a–2.2e to derive the relative magnification between counterimages (Fig. 3). This ratio is 1.66. If we assume that this magnification ratio is maintained between knots 2.1e and 2.2e, the critical curve is approximately southwest from knot 2.2e as shown in Fig. 3. (The distance between knots 2.1e and 2.2e is 505.) This is an approximation because the magnification ratio may change slightly as one approaches the critical curve, but it should be a good approximation. A pair of features near the expected position of the critical curve could be counterimages of each other bracketing the critical curve.

Fig. 3. Perturbers in the giant arc La Flaca. Scale and orientation are shown on the image. Color-coded circles mark multiple images of three background sources. La Flaca is extremely thin and one of the longest known arcs. It is located between the two main mass clumps in El Gordo. Small-mass perturbers create additional multiple images. These perturbers are marked with arrows. The approximate position of the critical curve can be estimated by extrapolating the ratio of distances between knots 2.1a–2.1e and knots 2.2a–2.2e. The point where one expects the critical curve is marked with the cyan circle having 1″ radius. The two cyan ellipses within this circle mark a possible double counterimage. If these two objects are counterimages of each other, the critical curve would be very close to the middle point between the two. The small cyan circle labeled A marks a possible caustic-crossing object.

2.2. HST data

In addition to JWST data, we used public HST imaging data from programs GO 12755 (P.I. J. Hughes), GO 12477 (P.I. F. High), and GO 14096 (P.I. D. Coe). These ACS and WFC3/IR observations include data in 10 filters spanning wavelengths ∼0.4–1.6 μm. The Reionization Lensing Cluster Survey (RELICS; Coe et al. 2019) delivered reduced images combining data from all of these HST programs, including their own (14096). The RELICS data release also includes galaxy catalogs with photometry and photometric redshifts.

The HST data add valuable information, especially at wavelengths shorter than 0.8 μm. The HST data also cover a wider field of view than JWST. This is important for lens modeling because some member galaxies are outside the field of view of JWST but within the HST field. An especially important one is a massive member galaxy at RA = 15.752995, Dec = −49.281048 and spectroscopically confirmed at z = 0.8749 (Caminha et al. 2022) that falls on the edge of the JWST image.

The color images in this paper include both HST and JWST data. A particular combination, which we call RGB6, offers a good compromise between low noise levels, sensitivity, and color range. RGB colors are defined by

(1)

(2)

(3)

Unless otherwise noted, all color images shown in this work are based on RGB6.

2.3. X-ray data

Chandra data from the ACIS instrument acquired in 2011–2012 (ObsID 12258, 14022, and 14023, PI.. J. Hughes) total ≈350 ks. The X-ray data were smoothed using the code ASMOOTH (Ebeling et al. 2006). Figure 4 shows the BCG galaxy with overlaid contours of the smoothed X-ray data. Gas falling into the BCG forms a filamentary structure that coincides with the peak X-ray emission, suggesting a powerful cooling flow. The infalling gas has strong emission in [O II] and correlates spatially with the peak of the X-ray emission. MUSE [O II] data reveal a velocity structure in the infalling gas. The intense X-ray emission combined with the infalling velocity structure suggest this filament is the optical counterpart of a powerful cooling flow, where large amounts of material fall toward the BCG. The strong [O II] emission in the filament can be understood as the plasma being chemically enriched by feedback process from the BCG itself or from material being stripped away from other member galaxies. Lindner et al. (2014) found an unresolved radio source “U7” at one of the extremes of the filament. No obvious counterpart for U7 can be identified in JWST images, but within the beam size of the radio data, U7 coincides with the peak of the X-ray emission and the infalling [O II]-emitting gas. U7’s radio emission may be synchrotron emission arising as the result of rapid deceleration of electrons where the cooling plasma encounters denser material surrounding the BCG.

Fig. 4. JWST image of the central BCG with smoothed X-ray contours shown in green. Scale bar shows physical distance corresponding to 256. The white ellipse marks the position and beam shape (ATCA at 2.1 GHz) of the unresolved radio source U7 (Lindner et al. 2014). Blue knots in and above the circle agree in position with [O II] emission seen by MUSE.

3. Lensing constraints

In the strong lensing regime, lens models are primarily constrained by the observed positions of multiply lensed galaxies. The first step is recognizing which observed objects constitute an “image family” of a single background source. Once a suite of image families from sources at different redshifts are available, the lens model is optimized by searching for a solution that focuses each family of images into a single position at the corresponding redshift. The redshift information is ideally obtained from spectroscopy. When this is not available, photometric redshifts can be used instead, or the redshift can be treated as a free parameter.

With the new data from JWST, some candidate image families found in earlier work (Zitrin et al. 2013; Cerny et al. 2018; Diego et al. 2020) are now confirmed, thanks to the improved resolution, depth, and color information. The improvement in spatial resolution and depth also unveils new candidates. We visually inspected the JWST images and compiled a new list of candidate families. Spectral information from MUSE and redshifts derived by Caminha et al. (2022) confirmed some of these candidates as real multiple-image systems.

Using MUSE data, Caminha et al. (2022) confirmed 23 families of lensed galaxies, 12 of which were previously known and 11 of which were new systems. System 8 of Caminha et al. (2022) is a redefinition of System 6 of Diego et al. (2020) but with two counterimages incorrectly identified. Of the 23 Caminha et al. (2022) systems, 12 had only two images visible. The third image is presumably too faint to be identified in the HST images or show up in the MUSE spectra. Thanks to JWST’s increased depth and spectral range, we found candidates for the third image for 8 out of these 12 systems. For the remaining 4, the counterimage is expected to be too faint to be visible even in the current JWST data. The list of all available constraints with spectroscopic redshift confirmation is given in Table A.1.

In addition to the classic constraints given by the observed positions of strongly lensed images, we also used the position of critical points that act as anchors for the critical curve. Critical points are positions that a critical curve is known to pass through. These can be identified as symmetry points in arcs that are formed by the merging of two images. Several of these merging arcs can be found in the JWST data. One needs to know the redshift of the arc in order to use the critical point as a valid constraint, and therefore only critical points of galaxies that have spectroscopic redshift can be used. There are three such critical points in systems 1, 2, and 23. Their positions are listed at the end of Table A.1.

4. Free-form lens modeling of El Gordo

4.1. Modeling method

The mass reconstruction is based on our method WSLAP+. Previous papers (Diego et al. 2005, 2007, 2016; Sendra et al. 2014) give details. As a brief summary, the lens equation is

(4)

where θ is the observed position of the lensed source, α is the deflection angle, Σ(θ) is the surface mass-density of the lens (El Gordo cluster in our case) at position θ, and β is the true position of the background source. Both the strong-lensing and weak-lensing data can be expressed in terms of derivatives of the lensing potential ψ:

(5)

where D_l, D_s, and D_ls are the angular-diameter distances to the lens, to the source, and from the lens to the source, respectively. The unknowns of the lensing problem are the surface mass density of the lens and the positions of the background sources in the source plane.

As shown by Diego et al. (2005, 2007), the lensing constraints can be expressed as a system of linear equations

(6)

where the observables (positions of strongly lensed galaxies, positions of critical points, and weak lensing if available) are contained in the array θ, and the unknown surface mass density, rescaling factors for the different layers discussed below, and true source positions are in the array X. The matrix Γ is known and given by the position of the grid points and positions of the constraints. In our case, θ contains the positions listed in Table A.1 (lensed galaxies and critical points). (No weak-lensing constraints were used in this work.) Details on how the critical points are added to the system of linear equations are given by Diego et al. (2022). A solution, or lens model, is found by minimizing a quadratic function derived from the system of linear equation with the constraint X > 0. That is, all masses and source positions (referred to a corner of the field of view) must be positive.

The model surface mass-density for El Gordo was described by a combination of two components: (i) a soft (or diffuse) component and (ii) a compact component that accounts for the mass associated with the individual halos (galaxies) in the cluster. The diffuse component was computed as a set of mass centers at predefined positions. Each mass center was assumed to represent a surface mass-distribution defined by a Gaussian with a predefined full width at half maximum (FWHM). The algorithm then optimized the mass of each mass center to best fit the constraints.

Mass distributions other than Gaussians could have been used, but Gaussian functions provide a good compromise between the desired compactness and smoothness of the basis function. A Gaussian basis offers several advantages, including a fast analytical computation of the integrated mass for a given radius, a smooth and nearly constant amplitude between overlapping Gaussians (with equal amplitudes) located at the right distances, and orthogonality between relatively distant Gaussians that help reduce unwanted correlations. (Diego et al. 2007 discussed alternative basis functions including polynomial, isothermal, and power laws.) The spatial distribution of the grid can be uniform or adaptive. The uniform distribution assumes a flat prior for the mas while the adaptive grid is based on a previous solution (for instance derived with a regular grid) and increases the resolution of the uniform component in areas where more mass is expected. This work explores both configurations, and Fig. 5 shows the adaptive grid used here.

Fig. 5. Multiresolution adaptive grid. The number density of grid points (shown as 891 red crosses) traces a smooth version of the solution in the spectroscopic lens model, increasing the resolution in the densest regions. Ten additional grid points are shown as small black crosses. Two of them mark the perturbers in the La Flaca arc with scales of 012 and 024, and eight more are in the inner region of the Anzuelo arc. Each one of these 8 grid points has a scale of 18.

For the compact component, we adopted directly the light distribution around the brightest member galaxies in the cluster. For each location, we assigned a mass proportional to the galaxy surface-brightness at that location with the mass-to-light ratio (M/L) being solved for as part of the optimization process. We used the surface brightness in the F160W HST image because its wide field of view includes some galaxies not covered during by JWST. Directly adopting the light distribution for these galaxies avoids additional modeling parameters. The compact component was divided into four independent layers with independent M/L because M/L can differ for different galaxy types. The first layer contains only the main BCG. The second layer contains the three galaxies acting on El Anzuelo. Layer 3 contains all remaining members in the cluster, and Layer 4 contains a small group of galaxies in the foreground at z = 0.63 (Caminha et al. 2022). As discussed by Caminha et al. (2022), although this is a small group, it has a non-negligible contribution to the lensing. For simplicity we assumed a single lens plane so this second group was assumed to be at the same redshift as the cluster. This results in a small bias for the true mass of this group, but this work needs only the effective lensing mass at the redshift of the cluster.

4.2. Spectroscopic lens model

We derived a first model (the “spectroscopic model”) using only the 23 families that have spectroscopic redshifts¹ These are identified as rank A in Table A.1. The soft-component mass-centers were defined on a uniform 32×32 grid. This is equivalent to adopting a flat prior for the mass distribution. The grid spacing was 46875, and the FWHM of each Gaussian mass center was taken to be 24.

The solution for the spectroscopic model is shown in Fig. 6. The distribution of the soft component (mostly dark matter) correlates well with that of the member galaxies. A clear peak in the soft component appears at the position of the BCG in the southeast part of the cluster. In contrast, the northwest clump has a less-concentrated peak in the soft component, in agreement with previous results.

Fig. 6. Spectroscopic lens model. Contours represent the mass distribution of the smooth component of the lens model, that is, dark matter, diffuse baryons such as stars from the ICL, and X-ray-emitting plasma. The galaxies used to describe the compact contribution to the lens model are shown in gray. North is up and east is left, and the image is 2.5 arcmin across.

Figure 7 compares the soft mass component with the X-ray emission from Chandra. The mass in the southeast clump peaks close to the position of the apparent cooling flow (Sect. 2). On larger scales, the distribution of X-rays shows a cometary structure similar to the Bullet cluster. The X-ray morphology suggests that the southern clump is moving in the southwest direction, and therefore that the cluster has already had one encounter. The northern clump is more irregular and may have been affected by the dynamics of the collision. There is an excess of mass east of the local X-ray peak. The south clump is much more massive than the northern one but also more elongated along the apparent direction of motion.

Fig. 7. Diffuse mass component from the spectroscopic model vs. X-ray emission from Chandra. Color shows the convergence at z = 3 as indicated by the color bar. Black contours show the X-ray emission from Chandra. The image is oriented north up, east to the left, and a scale bar is shown at top left.

4.3. Geometric redshifts from the spectroscopic lens model

The spectroscopic lens model relies only on confirmed systems, and therefore it gives a robust solution unaffected by redshift uncertainties of the background sources. This model can be used to confirm new system candidates and predict their redshifts. These redshifts are known as geometric redshifts (or geo-z) because they are the redshifts at which the multiple images from a given family focus at a single point. That is, the redshift corresponds to the focal plane of the lens for that particular family of images.

The spectroscopic model gives geometric redshifts for all system candidates listed in Table A.1. The probability of a system to be at redshift z is

(7)

where V(z) is the variance between the multiple positions of a given system projected on the source plane at redshift z. The projection was calculated from the deflection field of the spectroscopic model (computed at z = 3) which was rescaled for each system to the corresponding redshift. The dispersion, σ, in the expression above was fixed to 018 (about 3 pixels in the LW channels). This is a reasonable choice for well-constrained systems, resulting in relatively narrow distributions for the redshift. Systems that are well reproduced by the spectroscopic model give small V(z) near the optimal redshift, which in turn results in maximum values of P(z) close to 1. Systems that are poorly reproduced by the spectroscopic model have larger values of V(z) at its minimum, which reduces the maximum value of P(z). A low maximum probability for P(z) does not necessarily mean that the system is a bad candidate. This can simply be the result of the spectroscopic model not being well constrained in that part of the lens plane. Systems at high redshift tend to have broader probabilities because the deflection field varies slowly with redshift for z ≳ 4.

As a sanity check, we used the spectroscopic model to derive geometric redshifts of the systems used to build the model. The result is shown in Fig. 8. As expected, most systems were recovered with a redshift close to the true redshift and with small uncertainties. Systems 22 and 23 lie at the edge of the region that has fewer lensing constraints. Their uncertainties are larger, and the model may not be reliable in those areas. The only other clear outlier is system 16, which has only two counterimages (as does System 14, which agrees reasonably well). System 17 has a large uncertainty because its counterimages are close to each other and near member galaxies that are not individually constrained.

Fig. 8. Geometric redshifts compared with the input spectroscopic redshifts. This plot shows a sanity test where the input redshifts are properly recovered by the lens model. Families with discrepant redshifts are labeled. Black points show results of the spectroscopic model, and red points show results of a bootstrap analysis where a new model was derived without that particular system and that model used to predict the system’s redshift.

A better way to check the model is to bootstrap the lensed systems. We removed one system at a time and derived a new lens model, then used that model to predict the redshift of the system that was removed. In general the spectroscopic model still recovered the redshift of these systems with good accuracy. The exception are Systems 22 and 23, which are at the outskirts of the constrained region. Therefore removing one of these systems results in an accentuated degradation of the spectroscopic model. Where several spectroscopic systems exist, redshifts are well recovered. In earlier work, a similar bootstrap confirmed the excellent performance of geo-z estimates of well-calibrated lens models (Chan et al. 2017, 2020).

Geometric redshifts are a useful alternative to photometric redshifts, especially for high-redshift candidates for which photometry may be poor or nonexistent (for example in the case of dropouts). Also, photometric redshifts may be unreliable for dusty galaxies at z ≈ 4 because their Balmer break may be misinterpreted as Ly-α, placing them at much higher redshifts (Naidu et al. 2022; Zavala et al. 2023; Harikane et al. 2023).

To provide photometric redshifts, we used LePhare (Arnouts et al. 1999; Arnouts & Ilbert 2011). The setup was the same as used with NIRCam data by Adams et al. (2023). In brief, sources were selected in the F444W image using SExtractor (Bertin & Arnouts 1996). Its dual-image mode was then used to derive photometry (MAG_AUTO) in all other JWST filters. Galaxy templates were the BC03 set (Bruzual & Charlot 2003) with 57 ages, a Chabrier IMF (Chabrier et al. 2000), constant or exponential star formation histories, solar or 20% solar metallicity, dust extinction in the range of 0 < E(B − V) < 1.5 (Calzetti et al. 2000), the IGM treatment from Madau (1995), and 0 ≤ z < 25. Figure 9 shows the results. In general, photo-z agrees with spec-z and geo-z. The outlier points in the horizontal branch are mostly biased low photo-z estimates although three of the seven labeled outliers in Fig. 9 (Systems 4, 8, and 21) have secondary solutions much closer to the spectroscopic redshift. In a wider sample of galaxies with spectroscopic redshifts (not necessarily multiply lensed), ≈86% have LePhare redshifts with an error δz/(1 + z) < 0.15.

Fig. 9. Comparison of JWST photometric redshifts with other redshifts. Black points represent spectroscopic redshifts, and red points represent geometric redshifts, both of which are shown on the x axis. Dashed line shows equality, and discrepant systems are labeled. Only systems having consistent photo-z estimates for at least two counterimages are shown.

A general problem with photometric redshifts based on NIRCam observations alone is that they can give multiple redshift solutions. Some of this degeneracy can be rectified by the inclusion of other data such as from HST. Also, residual systematics in zero-point calibrations may still be affecting the photo-z estimates with JWST (Boyer et al. 2022). RELICS photometric redshifts (based on BPZ: Benítez 2000), given in Table A.1), offer examples of the improvement that can be obtained. In particular, adding HST data gives fewer small-phot-z outliers. For example, for System 8 with z_spec = 4.3175, RELICS predicts z_phot ≈ 4.3 while JWST photometry alone predicts z_phot ≈ 0.7. System 10 at z_spec = 4.3275 is consistent with the RELICS z_phot ≈ 4.2–4.6 and in better agreement than JWST-based z_phot ≈ 0.7. This is almost identical to System 12 (z_spec = 4.7042), where HST z_phot ≈ 4.6–5.3 while JWST z_phot ≈ 0.7. System 21 at z_spec = 5.5811 is predicted to be at z_phot ≈ 6 according to HST while JWST predicts z_phot ≈ 1.3. Frye et al. (in prep.) will provide additional photometric redshifts that combine HST and JWST photometry. The important lesson from Figs. 8 and 9 is that geo-z can be reliable estimates of the redshift.

4.4. Lens model with geometric redshifts

One of the factors determining the uncertainty of a lens model is the possible error in the redshifts. This is particularly problematic for high-redshift galaxies that are visible only in the reddest JWST bands, mimicking galaxies at even higher redshifts as recently studied by Naidu et al. (2022), Zavala et al. (2023), Harikane et al. (2023). The spectroscopic model is a robust solution because it relies on systems with spectroscopic redshifts, but the model is limited by the redshifts available.

In order to improve the resolution of the spectroscopic model, we need to include systems without spectroscopic redshifts. To do so without adding bias from inaccurate photometric redshifts, we relied on the geometric redshifts. This adds 37 additional system candidates labeled in Fig. 1 (bottom). The “Full-r model” maintains the 32 × 32 regular grid used in the spectroscopic redshift model and simply increases the number of constraints by using all 60 systems with rank A or B listed in Table A.1.

The Full-r model solution is shown in Fig. 10. As in the spectroscopic model, there is a good correlation between the smooth and compact components. The new constraints reveal more detail in the distribution of the smooth component. The clump in the northwest shows a more irregular distribution than in the spectroscopic model. In between the two clumps, we see excess mass going in the east–west direction. This excess mass correlates well with the wings of the X-ray emission. Simulations of the cluster collision (Molnar & Broadhurst 2015) show similar wings forming in the direction perpendicular to the axis of the collision provided the impact parameter is ≤100 kpc and infalling velocities are ≳2500 km s⁻¹. The smooth component traces mostly the dark matter and the X-ray-emitting plasma, and the excess mass in the direction perpendicular to the collision axis could be partially due to the plasma. We also see an elongation of the smooth component in both clumps, especially the southern clump. The peak of the smooth component correlates well with the peak in the X-ray emission (Fig. 11) although there is a small offset. This is partially due to the possible cooling flow that is the brightest feature in X-rays. Finally, the amplitudes of both mass peaks are slightly smaller in the Full-r model than in the spectroscopic model.

Fig. 10. Diffuse mass distribution and galaxy locations for the Full-r model. Contours represent the smooth component (dark matter and diffuse baryons such as stars from the ICL and X-ray-emitting plasma). The galaxies used to describe the compact contribution to the lens model are shown in gray.

Fig. 11. Diffuse mass component from the Full-r model vs. X-ray emission from Chandra. Color shows the convergence at z = 3 as indicated by the color bar. Black contours show the X-ray emission from Chandra. The image is oriented north up, east to the left, and a scale bar is shown at top left.

4.5. Lens model with adaptive grid

To refine the mass distribution, we built a third model called “Full-a.” It has the same constraints as the Full-r model but placing the soft-component mass-centers on an irregular grid. This provides increased resolution in regions with more mass. The grid was built from a Monte Carlo realization of a smoothed version of the mass distribution of the spectroscopic model, and instead of fixed sizes of 24, the 891 main mass centers have widths proportional to the average separation from their neighbors. This grid is shown in Fig. 5. It can be interpreted as using the prior that the mass distribution is the one indicated by the spectroscopic model, but of course the model optimization allows the derived mass distribution to deviate from the prior. The grid also includes two Gaussians with widths 012 and 024 at the position of two perturbers near the arc nicknamed La Flaca (Fig. 12) and eight additional Gaussians with FWHM 18 around the group lensing the galaxy El Anzuelo. The perturbers in La Flaca have to have small FWHMs because the perturbers cannot constrain masses larger than the image separation.

Fig. 12. Small-scale deflectors. The giant, thin arc nicknamed La Flaca runs horizontally across the image. The heavy red, green, and blue curves are the critical curves of El Gordo at the redshift of the arc (zspec = 2.8254) for the spectroscopic, Full-r and Full-a models, respectively. The Full-a model (but not the other models) includes small deflectors at the two positions marked by the white and yellow arrows. Dotted blue curves show the critical curves of the deflectors. The deflector marked with a yellow arrow has a MUSE spectrum with z = 0.7842. With a Gaussian FWHM of 024, its mass is 2.7 × 1010 M⊙. The deflector marked with a white arrow is barely detected, and we have assumed it is at El Gordo’s redshift. With , its mass is 3.8 × 109 M⊙.

The Full-a model reveals even more details as shown in Figs. 13 and 14. The mass distribution has a more concentrated peak around the BCG in the southeastern clump, while the northwestern clump is less centrally concentrated though about the same mass. Neither the Full-a model nor any of the others shows excess mass at the position of the [O II] filament. This suggests that the contribution of the dense filament to the projected mass must be relatively small. The critical curves for the three lens models considered in this work are shown in Fig. 15, assuming a source at redshift z_s = 2.8254.

Fig. 13. Mass vs. galaxies for the Full-a model. Contours represent the mass distribution of smooth component (dark matter and diffuse baryons such as stars from the ICL and X-ray emitting plasma). The galaxies used to describe the compact contribution to the lens model are shown in gray.

Fig. 14. Mass vs. X-ray for the model Full-a. The color plot shows the convergence (at z = 3) for the spectroscopic model. Only the dark matter component of the lens model is shown. X-ray emission from Chandra are shown as black contours. It is important to note the excess of mass near the X-ray local peak at the east. The south clump is now clearly more massive than the northern one but also more elongated along the direction of motion.

Fig. 15. Comparison of the z = 2.8254 critical curves for the three models. The red curve corresponds to the spectroscopic model, the green curve is for model Full-r, and the blue curve is for model Full-a. The color image is made from the combination of HST+JWST filters RGB6 discussed in Sect. 2.

5. El Gordo mass estimate

All mass models show a clear double-peaked distribution. The southern group is centered near the BCG, and the northern group is centered near a luminous galaxy at , . Masses of each clump for the various models are given in Table 2. The table includes masses from the parametric model of Caminha et al. (2022) for comparison. The Caminha et al. (2022) masses are lower because they exclude the mass of member galaxies, while we quote total masses. Otherwise all masses agree reasonably well.

The virial mass in El Gordo has been used in earlier work to study possible tension with the standard cosmological model, which predicts the most massive cluster at this redshift should have a virial mass M_200c ≲ 2 × 10¹⁵ M_⊙ (Harrison & Coles 2012; Watson et al. 2014). Given the similarity in mass of the two clumps, the cluster center of mass should be near the midpoint between the two groups, and we adopted a galaxy at RA = 15.7313639, . Figure 16 shows the integrated masses as a function of radius for the models derived here and some others. At 500 kpc, the masses for the models described here are between 8.0 and 8.6 × 10¹⁴ M_⊙. These are higher than the Diego et al. (2020) model based on RELICS HST data (M_500 kpc = 6.6 × 10¹⁴ M_⊙) but in reasonable agreement given the different data sets. The Caminha et al. (2022) model used exactly the same constraints as our spectroscopic model but a completely different approach. The agreement within the constrained region is excellent, and M_500 kpc = 7.98 × 10¹⁴ M_⊙, also in good agreement with our estimates. Finally, the N-body simulation of Molnar & Broadhurst (2018) agrees very well with the lensing models up to R ≈ 500 kpc, the limit of the strong-lensing constraints, and M_500 kpc = 7.20 × 10¹⁴ M_⊙, also in good agreement. Above 500 kpc, one expects El Gordo’s mass profile to continue as in the simulation. The virial radius is 1.75 Mpc, and the simulation gives M(R_vir) = 1.88 × 10¹⁵ M_⊙. At 500 kpc, the lens models give masses 11–19% higher than the N-body model, and one might correct the N-body mass higher by the same factor. This is at least close to, and probably above, the cosmological limit. A more robust estimate of the virial mass could be obtained by combining the strong lensing constraints with weak lensing measurements extending beyond the virial radius.

Fig. 16. Total El Gordo mass as a function of radius. The red curve shows values for the spectroscopic model, the green curve for the Full-r model, and the cyan curve for the Full-a model. The purple line shows the mass obtained derived from the RELICS data using the same algorithm but different constraints (Diego et al. 2020). The orange dot-dashed line shows the recent model of Caminha et al. (2022). The black dashed line shows values from a simulation (Molnar & Broadhurst 2015). All mass curves are based on a common center, and the masses are cylindrical, that is, they correspond to the mass within a cylinder of radius R projected along the line of sight. The light turquoise line shows the mass enclosed in a sphere with constant density equal to 200 times the critical density. While the spherical mass is not identical to the cylinder mass, for an NFW with concentration c = 7, the 3D mass computed up to the virial radius is only ≈4% smaller than the 2D mass. The vertical dotted line marks the maximum radius (70″) at which lensing constraints exist. The integrated masses will be biased low beyond this point.

6. Small-scale substructures in La Flaca

As discussed in Sect. 2 and shown in Fig. 3, the giant arc La Flaca contains two small-scale perturbers. We can use the observed lensing distortions to infer their masses as WSLAP+ optimizes the mass in all grid points. The model Full-a contains two small Gaussian at the perturbers’ positions. This model finds masses of 2.7 × 10¹⁰ M_⊙ and 3.8 × 10⁹ M_⊙ for the larger and smaller perturbers, respectively. Profiles other than Gaussian should change the masses within the Einstein radius of the perturber by no more than a factor of a few tens of percent. Figure 12 shows how the critical curve of the larger deflector winds around knots 2.2b, 2.5b, and 2.4b to reproduce the observed triple image. The curve of the smaller deflector passes between knots 2.2a and 2.4a. We expect 2.4a to be a double image, but if so, it is unresolved in JWST images. The critical curve is close to the middle point in 2.4a but is still off by a small fraction of an arcsecond, suggesting there is still room for a small improvement in modeling this double knot.

The masses inferred by the lens model are consistent with small galaxies. In particular, the smaller 2.4a perturber has mass consistent with being a dwarf galaxy if it is at El Gordo’s redshift. As far as we know, this would be the smallest mass measured at z > 0.5. The effective mass of a perturber scales as its mass times the macromodel magnification, and new JWST images near critical curves perturbed by small halos will allow for masses for even smaller substructures to be measured, potentially constraining models of dark matter, as Diego et al. (2022) recently did.

Assuming M/L ∼ 100 for the smaller perturber, the stellar component of such a galaxy would be barely detected, especially if the light is spread over hundreds of parsecs. This is consistent with the faint detection in the JWST data. A precise measurement of the luminosity of the smaller perturber is challenging because its light overlaps with a bright feature in the giant arc.

7. Caustic-crossing candidates

Among the greatest achievements by HST in the last few years was its discovery of extremely magnified stars at high redshift. Such objects were first predicted by Miralda-Escude (1991). The discovery of Icarus (Kelly et al. 2018) marks the starting point of this new field, which allows study not only of individual stars at z > 1 but also compact star clusters (Dai 2021). The discovery also opens the door to novel studies of dark-matter structures (Diego et al. 2018, 2022; Oguri et al. 2018; Venumadhav et al. 2017; Dai et al. 2018, 2020; Dai & Miralda-Escudé 2020; Meena et al. 2022). After the discovery of Icarus, other examples followed (Chen et al. 2019, 2022; Kaurov et al. 2019; Diego et al. 2022), culminating with the recent discovery of Earendel at z ≈ 6.2 (Welch et al. 2022a,b).

JWST is expected to see farther than HST and possibly detect individual stars close to the beginning of cosmic reionization (Windhorst et al. 2018). The first public data from JWST already show candidate lensed stars at cosmological distances (Pascale et al. 2022). A search for similar candidates in El Gordo found a clear example in a strongly lensed arc shown in Fig. 17. Given its proximity to the midpoint between the pair of knots 23.4a and 23.4b (yellow circles in Fig. 17), this source (marked with a white circle in Fig. 17) is likely right on the critical curve and therefore highly magnified. Because images near a critical curve always come in pairs, the lack of counterimage and unresolved nature of the source implies that the image pair must have a counterimage separation no more than approximately 1 pixel or 30 mas. This situation is similar to the cases of Godzilla (Diego et al. 2022) and Earendel (Welch et al. 2022a,b), where those sources must also form unresolved pairs of counterimages.

Fig. 17. Possible caustic crossing event in system 23 at z = 2.1878. Top: Color image combining F115W, F200W, and F356W as blue, green, and red, respectively. The outer, yellow circles mark the positions of 23.4a and 23.4b, two counterimages of a source that brackets the position of the critical curve. The white, central circle marks a bright source (called “Quyllur”) that lies approximately at the midpoint. Bottom: color image with the RGB6 combination. Red, green, and blue curves show the critical curves for the spectroscopic, Full-r, and Full-a models respectively. The two panels have the same scale and orientation as indicated at the left.

Several alternative possibilities are discussed below, but the most likely explanation for the source is a single, highly magnified, red supergiant star. We nickname it “Quyllur,” which is the Quechua term for star².

Figure 18 shows Quyllur’s SED. The SED shows a rapid decline below 3 μm consistent with a surface temperature T ≈ 3500 K at z = 2.1878, the redshift of 23.1a. Assuming 30 mas as the maximum separation between counterimages, and accounting for both images (a factor 2 in flux), the magnification must be μ > 4000, a boost of at least 9 magnitudes when considering the flux from both images. If the separation between images is smaller than 30 mas, then the magnification can be even larger. The observed 3.6 μm magnitude is 25.5, which corresponds to > 34.5 unmagnified. At a distance modulus of ∼45, and depending on the magnification, the absolute magnitude at 1.1 μm rest wavelength would be ∼−10.5, consistent with an M2–M4 supergiant.

Fig. 18. Quyllur compared to SEDs of single stellar populations at El Gordo’s redshift z = 0.87. Blue lines show the unreddened models as labeled, and green lines show the models reddened with E(B − V) = 3.5 mag. Black points with error bars show the Quyllur photometry. For a comparison with the SED of individual stars see Fig. 21.

The required magnification limits the source size and therefore the possible source types. For μ = 4000, the distance to the caustic must be less than 7.5 microarcsec or ≈0.06 parsec. The only known sources that can satisfy the size and luminosity constraints are supergiant stars or accretion discs around very massive back holes. Some young open clusters with an age around 10–20 Myr show remarkably rich populations of red supergiants (Davies et al. 2008, 2007; Alexander et al. 2009; Froebrich & Scholz 2013). If the source is not a single red supergiant star, it might be a group of red supergiant stars in an open cluster. This scenario is however disfavored because the constraint on the source size ≲0.06 pc is significantly smaller than the ≳several pc typical sizes of open clusters. Moreover, red supergiant stars are short-lived, and it is likely that a group of them would be accompanied by B or even O stars, whose bluer light would be detected unless dust extinction is severe.

The absence of emission shortward of 1.5 μm (0.47 μm rest) argues against an accretion disk unless there is extinction corresponding to E(B − V) > 4. Another argument against an accretion disk is that source-plane reconstruction through the lens model reveals that the source is not located near the nucleus of the host galaxy. Photometry performed between Quyllur’s position and 23.4a and 23.4b yields a color index E(G − J)≈2.5. (This color index corresponds to the difference between the F410M and F150W filters.) Quyllur has a much redder index E(G − J)≈5. This rules out significant reddening affecting the portion of the galaxy surrounding Quyllur and also rules out multiple images because only Quyllur shows such extreme color.

The possibility that Quyllur is a transient event showing (at the time of observation) only one of the two lensed images is unlikely based on the expected time delay. Based on the time separation between neighboring knots 23.4a and 23.4b (see time delay column in Table A.1), the time delay between Quyllur and an unseen counterimage must be significantly shorter than 0.28 year because it is closer to the critical curve. (Time delays are inversely proportional to magnification.) For Quyllur’s position, the lens model predicts a magnification μ ≳ 30″/d, where d is the distance to the critical curve in arcseconds measured along the stretched arc. In order to observe a transient event such as a supernova and not its counterimage, the JWST observation must have taken place within the first days of the event, an unlikely coincidence.

Lacking spectroscopic confirmation, we consider the possibility that Quyllur is an interloper galaxy, most likely in the El Gordo cluster. However, Quyllur’s SED is remarkably red. A quiescent stellar population even as old as the age of the Universe (6.3 Gyr at z = 0.87) would appear bluer unless dust extinction is also present. As shown in Fig. 18, Quyllur can be consistent with either a 6 Gyr stellar population with stellar mass M_* ∼ 4 × 10⁹ M_⊙ or a 100 Myr stellar population with M_* ∼ 2 × 10⁸ M_⊙. Either would be of substantial mass but would have to be more compact than ∼0.5 kpc, that is, a dwarf spheroidal (dSph) galaxy. However, either scenario requires dust reddening E(B − V) = 3.5 (assuming Milky-Way-like reddening with R_V = 3.1). Without observational evidence for such a galaxy population inside clusters, the scenario seems less credible than a caustic event. If Quyllur is an interloper with z < 0.87 or moderately higher, it would still be unusually red and would require dramatic dust reddening following the same argument. If the redshift is even higher, it would be strongly lensed by El Gordo and show counter-lensed images, which are not found.

A second candidate to be a compact source crossing a caustic is within the merging arc of system 1 (Fig. 19). As before, the two knots 1.3a and 1.3b are very close to the critical curve and determine the position of the critical curve. A very faint feature is northwest of this pair and is another candidate for an extremely magnified source no more than a fraction of a parsec from the caustic. The pair of knots 1.3a and 1.3b themselves have magnification factors > 100 and are interesting objects to study in more detail.

Fig. 19. Possible caustic crossing in system 1. Top: Color image combining F090W, F200W, and F356W as blue, green, and red, respectively. Bottom: Color image in the RGB6 pallette with critical curves from the spectroscopic, Full-r, and Full-a models shown as red, green, and blue, respectively. In both images, the two yellow circles mark images 1.3a and 1.3b, a pair of counterimages close to the expected position of the critical curve. The white circle contains a faint source that is barely resolved.

As discussed in Sect. 2, a possible additional candidate is barely detected in the giant arc La Flaca. However, similar to the candidate in system 1, it is too faint to extract any useful information. Finally, as shown in Fig. 2, small, unresolved image pairs can be found in the Anzuelo galaxy. These are not as magnified as Quyllur and therefore are likely bigger and brighter sources such as compact star-forming regions or globular clusters.

8. Discussion

In contrast to the lensing cluster WHL J013725.2+140341 (Welch et al. 2022a,b), where few new families of lensed galaxies have been identified, El Gordo acts as a magnificent lens amplifying the flux of tens of distant galaxies. This is partially due to the existence of an overdensity at z ≈ 4.3 (including some of the multiple images already identified by Zitrin et al. 2013), first identified by Caputi et al. (2021) and later confirmed by Caminha et al. (2022) thanks to ALMA and MUSE data respectively, that intersects the caustic region of El Gordo at this redshift. The dual mass peaks and resulting stretched caustics in El Gordo also help increase its cross section for large magnification factors. In contrast, earlier lens models of WHL suggest a rounder shape (Welch et al. 2022a), resulting in a more concentrated caustic region. These characteristics make El Gordo a good target for lensing studies.

8.1. Possible tension with ΛCDM

The highest mass estimate for El Gordo, M_vir ≈ 2.8 × 10¹⁵ M_⊙ (Jee et al. 2014, their Table 2 based on weak lensing) has significant tension with ΛCDM. Our M(vir)≈2 × 10¹⁵ is easier to accommodate, but extrapolating the mass from the region constrained by strong lensing out to the virial radius is less than ideal. Recent work (Asencio et al. 2021; María Ezquiaga et al. 2022) has used weak-lensing limits to illustrate the tension of this cluster with ΛCDM. At the core of this debate is the true mass of El Gordo. More recent work (Table 1) has tended to give lower masses, but the mass of El Gordo remains uncertain within a factor of two. Extrapolated masses are less reliable than this, and we cannot reach firm conclusions about the virial mass. Better weak-lensing data are needed to settle the debate. Ideally, weak-lensing measurements should be combined with strong-lensing constraints in a self-consistent manner, and the strong-lensing mass presented in this work can be used to break the mass-sheet degeneracy in future studies. Strong lensing anchors the mass within the Einstein radius, and weak lensing can extend the mass estimation to beyond the virial radius. The addition of new spectroscopically confirmed strong lensing systems, especially at the highest redshifts, can also help reduce the uncertainty in the mass. Additionally, galaxy–galaxy lensing in the outskirts of El Gordo can be used to constrain the cluster potential, and hence mass, at large radii, because the effective lensing mass of member galaxies scales as their true mass times the large-scale magnification. Constraints on the magnification at a larger radius would constrain the mass within that radius. The limited spatial coverage of our JWST data does not allow a weak-lensing analysis beyond the virial radius, but additional pointings that extend the coverage would better constrain the virial mass.

8.2. Flux ratios

Comparing the observed flux ratios to the predicted magnification ratios can identify regions where lensing models are inadequate. Differences between the ratios can arise from substructures missing from the lens model or from transient events affecting some but not all images in a family. Figure 20 shows that most images fall within a factor two of expectation, comparable to the expected uncertainty in the lens models. (Uncertainties are typically larger for large magnification factors (Zitrin et al. 2015; Meneghetti et al. 2017)).

Fig. 20. Model magnification ratio versus measured F200W flux ratio for image pairs with accurate photometry. The flux ratio is defined as the counterimage with the largest flux divided by counterimage with the smallest flux. Solid lines mark equality, and dashed lines mark a factor of two deviation, as expected from uncertainty in the model magnifications. Points are color-coded for each model as shown in the legend, and two outliers are marked.

Only one image family with spectroscopic redshift, System 18, has a flux ratio clearly deviating from the expectations of all three lens models. The models predict similar magnification for images 18b and 18a, yet the observed flux of 18b is ≈7 times larger than that of 18a. Image 18a has negative parity and is close to a nearby member galaxy, and its flux may be demagnified by the galaxy. This would require the galaxy to be more massive than predicted by the lens model. (This particular galaxy was not optimized individually but as part of the collection of member galaxies in layer 3 of the lens model.) System 30 is also an outlier, in particular for the spectroscopic model that does not use this system as a constraint. The Full-a model optimizes the local mass density, and the model ratio falls more in line with the observations.

8.3. Quyllur

Very luminous stars are rarely found isolated and normally live near stellar groups. The multiply lensed pair 23.4a–23.4b is at a projected distance of only 10–15 pc from Quyllur, depending on the magnification. This bluer source could be a star-forming region with Quyllur in its outskirts. Assuming Quyllur is a red supergiant, the magnification is likely significantly higher than our lower limit ≈4000. Betelgeuse (d = 197 pc with luminosity 1.26 × 10⁵ L_⊙) has a flux density of ≈3 × 10⁻⁹ ergs s⁻¹ cm⁻² Å⁻¹ at 6500 Å (Levesque & Massey 2020). This wavelength would redshift into the F200W JWST filter, for which we find M_AB ≈ 28, (see Fig. 21). At the distance of Quyllur (d_L = 1.76 × 10¹⁰ pc) and without magnification, Betelgeuse would have an apparent magnitude of ≈39.75 − 2.5 × log₁₀(1 + z) = 38.5 in F200W (ignoring precise color- or K-correction and extinction). This requires magnification ≈15 000 to match Quyllur’s apparent magnitude in F200W. Red hypergiants (for example UY Scuti and Stephenson 2-18) are rare but have luminosities that can exceed that of Betelgeuse by factors of 3–5 (Humphreys & Davidson 1979). At these luminosities, the required magnification drops to ≈4000–7000 (Fig. 19).

Fig. 21. SED of Quyllur compared to stellar models. The black dots show the observed aperture-corrected magnitudes. The blue, green, and orange curves show three models from Coelho (2014) redshifted to z = 2.1878 and magnified by a factor of 20 000. Model temperatures are shown in the legend. Squares show the expected flux in the JWST bands for a star with T ≈ 3500 K and luminosity 105.1 L⊙, similar to Betelgeuse. A star with a similar temperature but being a few times brighter, such as UY Scuti, would require magnification to be less than 10 000.

At magnification factors of several thousand, microlenses are expected to introduce temporary fluctuations in the flux even for a relatively low surface mass-density of microlenses (Σ ≲ 10 M_⊙ pc⁻²). This is a typical value for the surface mass-density of stars at distances of ≈100 kpc from the BCG (Montes & Trujillo 2022), the approximate distance from Quyllur to the BCG. The total projected mass around critical curves is usually about 2500 M_⊙ pc⁻², and 1% of this in the form of stars is well within the range for microlensing to be common. At the redshift of El Gordo, its intracluster light and stellar surface mass-density are difficult to measure, but we can approximately estimate the surface mass-density of microlenses if we conservatively assume the stellar mass contributes ≈0.2%–1% to the total mass. At the redshift of El Gordo, the ICL is still forming, so we should expect no more than 1% contribution to the projected mass from stars at 100 kpc from the BCG (Montes & Trujillo 2022). At the position of Quyllur, the convergence from our lens models κ = 0.63–0.66, and the critical surface mass-density for the redshifts of El Gordo and Quyllur is 2300 M_⊙ pc⁻². This would make the surface mass-density of microlenses between 3 and 15 M_⊙ pc⁻². Combined with the model magnification of a few thousand, this implies an effective surface mass density, Σ_eff = μ × Σ, much larger than the critical surface mass density. In this regime, microlensing events are common (Venumadhav et al. 2017; Diego et al. 2018), so we should expect changes in the flux of Quyllur by approximately a factor ≈2 when observing at intervals separated by several months (Welch et al. 2022a). It is possible that we are observing Quyllur during one of these episodes of increased magnification due to microlenses. Alternatively, red supergiants can exhibit quasi-regular changes in their flux, and we could be observing Quyllur during an episode of increased luminosity.

The ubiquitous presence of microlenses sets a natural limit to the maximum magnification and hence minimum luminosity of detectable stars. Magnifications of order 10⁵ are possible only very close to the cluster caustic and only for very low surface mass-densities of microlenses Σ_* < 1 M_⊙ pc⁻² (Diego et al. 2018; Venumadhav et al. 2017). Typical values of Σ_* are above this, limiting the maximum theoretical magnification to < 10⁵, as for Earendel (Welch et al. 2022a). Moreover, at these short separations from the cluster caustic, the stars’ relative motion means they can maintain such large magnifications for only a few years, making these extreme magnification events very unlikely. Moderate magnification factors of order 10⁴ or less are two to three orders of magnitude more likely and are sufficient to explain Quyllur as a red supergiant sar.

If Quyllur is confirmed to be a red supergiant, it would be the first example of many to come. So far, all previous stars discovered near a caustic and extremely magnified have been hot and blue stars. This is partially an observational bias because cool stars, such as red supergiants, at high redshift emit most of their light at wavelengths that are beyond the reach of HST. JWST is most sensitive at these wavelengths (Dai et al. 2018), but red supergiants are intrinsically rarer than blue supergiants because they correspond to a shorter-lived stage of massive-star evolution.

At larger magnification factors, the prospect of detecting the intrinsically fainter red giants at cosmological distances opens the interesting prospect of using them as standard candles (TRGB or “tip of the red giant branch”). The TRGB is commonly used in our local Universe as part of the distance ladder. Extending this ladder to z ≳ 1 would provide an alternative to type Ia supernovae with different (and sometimes easier to manage) systematic effects. The uncertain magnification may be the main obstacle to use these newly found stars in practice although due to the lower luminosity of red giants compared to red supergiants, only those at the most extreme magnification values are expected to be observed. Monitoring for lensing events (that scale as the unknown magnification) will constrain the magnification and open the door to using these stars as standard candles up to z ≈ 1 and possibly beyond.

Other exotic objects show a red SED and could in principle mimic the observed photometry of Quyllur. One example is luminous transients that can temporarily reach ∼10⁶ L_⊙ and can be linked to proto-planetary nebulae (Prieto et al. 2009). Another possibility is a compact luminous object such as an accretion disk around a black hole. We have considered two hypothetical scenarios of obscured AGN that could reproduce the observed colors and magnitudes at redshifts 2.54 < z < 4.4. However, at these redshifts we would expect counterimages with similar colors and fluxes which are not observed.

Regarding the probability of observing a star similar to Quyllur, one interesting aspect to consider is the large volume accessible through lensing, V ≈ 1.5 × 10¹¹ Mpc³ between z = 2 and z = 2.3 compared with V ≈ 3.3 × 10⁴ Mpc³ up to ≈20 Mpc, which is the largest distance to which we have observed the brightest stars in our local Universe. In other words, through lensing we can probe a volume ≈5 × 10⁶ times larger at 2 < z < 2.3 than locally. With a larger volume and the lower metallicity expected at z ≈ 2, one should expect to see even brighter stars, and μ ≈ 1000 may suffice to observe them with JWST. Only a small fraction of the stars at z > 2 can be magnified above μ ≈ 1000 because P(μ > 1000)≈10⁻⁸ (Diego 2019). Therefore we would naively expect to see O(1000) stars between 2 < z < 2.3 in the entire sky, assuming there are O(1000) stars in our local r < 20 Mpc volume with L > 10⁵ L_⊙ and that these are ∼10 times more abundant at z > 2. In the next years, both JWST and HST will continue increasing the number of extremely magnified stars, significantly improving the statistics.

8.4. Dearth of high-redshift galaxies behind El Gordo

One of the promises of JWST is the observation of the most distant galaxies. In the case of El Gordo, the existence of overdensities at z > 4 that fall between the caustics at those redshifts offers a unique opportunity to study these groups in greater detail (Frye et al., in prep.). Our new system candidates include many galaxies at high redshift, but the lens models do not place any of them at z > 6. Because gravitational lenses such as El Gordo amplify faint objects, the probability of observing a z > 6 galaxy is non-negligible. Average magnifications μ ≳ 10 are possible near the critical curve. Galaxies found in these areas should be stretched in a well-predicted direction. Simply scanning the predicted high-z critical curve, looking for red sources stretched in the expected direction, can easily identify candidates to be high-redshift galaxies. A search near the z = 10 critical curve revealed some candidates, but all of them are very faint and are consistent with being stretched arcs at lower redshifts. The clearest candidate to be a strongly lensed galaxy at high z is at RA = 15.730584,  Dec = −49.271130. It is stretched to a length of ∼13 in the expected direction, and its position intersects the expected high-redshift critical curves. No other clear candidates were found. For a more detailed search in this cluster using JWST data see Bhatawdekar et al. (in prep.).

If we assume that there are no highly magnified galaxies at z > 6 behind El Gordo, we can estimate the likelihood of this lack of galaxies given the lens model and the expected luminosity function at high redshift. The area with total magnification > 30 is easily computed from the lens model by inverse ray tracing. (Total magnification means after adding all multiple images.) In a typical lensing configuration involving large magnification factors, three images form with two of them carrying most of the magnification, and the third one being significantly less magnified and in some cases falling below the detection limit. (This is the case for some of the systems in Table A.1.) We adopt a ratio of 2.5 between the total magnification and the magnification of the brightest counterimage. As shown by Vega-Ferrero et al. (2019), this is a common ratio found in the Hubble Frontier Fields clusters for systems with large magnification factors. This translates into magnification factors of at least 12 (or 2.7 mag) for the brightest counterimage. Our three lens models give areas at z = 6 with μ > 30 of ≈0.015 arcmin². This is comparable to the area of the six Hubble Frontier Fields clusters, placing El Gordo as the second most efficient lens at high redshift, behind only the train wreck MACS J0717.5+3745 (e.g., Vega-Ferrero et al. 2019).

The UV luminosity function (Ishigaki et al. 2018) at z = 6, derived from galaxies behind the Hubble Frontier Fields, reaches M_UV ≈ −16. A galaxy at z = 6 (distance modulus 48.85) with absolute magnitude −16 and one of its counterimages having magnification 12 would appear at apparent magnitude ≈30, that is, within reach of JWST.

The volume between redshifts 5.5 < z < 6.5 is V ≈ 360 Gpc³ or 2400 Mpc³ arcmin⁻². The volumetric number density of all 5.5 < z < 6.5 galaxies with M_UV < −16 is N ≈ 0.1 galaxies per Mpc³. Hence we expect a surface number density, N(z ≈ 6) = 240 galaxies arcmin⁻² in the same redshift interval. Most of these galaxies would be out of reach for JWST without gravitational lensing, but as stated above we expect ≈ 0.015 arcmin² behind El Gordo to be magnified by at least a factor 30, with one counterimage having magnification at least 12. That is, we expect ≈3.6 galaxies at z ≈ 6 strongly magnified and potentially detectable in the JWST images. At higher redshifts the magnification needs to be larger in order to compensate for the increase in distance modulus. The area above a given magnification μ scales as μ⁻², reducing the number of detectable galaxies. Also, the galaxy number density is smaller at higher redshift, so we expect many fewer galaxies at z ≫ 6 than at z ≈ 6.

A deficit of strongly magnified galaxies z > 6 could be the result of a patchy reionization scheme, where the UV and visible emission from these galaxies is still being absorbed by a neutral intergalactic medium. Alternative models based on fuzzy dark matter also predict a lack of galaxies at high redshift due to the existence of a minimum scale, which is determined by the mass of the dark matter particle (Leung et al. 2018). More clusters are needed to improve the statistical significance.

9. Conclusions

The El Gordo cluster is a unique object and of great relevance not only to understand galaxy evolution but also for cosmology. With its large mass and redshift, it has been at the center of a debate regarding a possible tension with ΛCDM. In the new JWST data, we have identified 37 lensed-system candidates (with 28 of them being new ones), which can be added to the 23 spectroscopically confirmed systems from MUSE (Caminha et al. 2022) to create three new free-form lens models. The models show correlations between the derived mass distribution and the X-ray emission. Extrapolating from a radius of 500 kpc, where the lens models constrain the mass, to the virial radius of ∼2 Mpc gives a mass of about 2.1 × 10¹⁵ M_⊙, which is close to the limit for standard ΛCDM cosmology. Future analyses combining weak and strong lensing will provide observational mass constraints out to the virial radius.

Interesting individual objects found include a candidate red supergiant star (nicknamed Quyllur) at z = 2.1878 and with an estimated magnification of at least 4000. The color of Quyllur is consistent with a red supergiant star, in contrast to other distant, magnified stars, which all appear as blue supergiants. This is in line with the prediction that JWST should greatly outperform HST in its ability to detect highly magnified stars with cool surface temperatures (Dai et al. 2018). Another object of interest is the smallest substructure measured to date at z = 0.87 with an estimated mass of 3.8 × 10⁹ M_⊙. This object creates additional images in portions of “La Flaca,” a lensed galaxy 22″ long. Finally, “El Anzuelo” is the image of a submillimeter galaxy that can be studied in great detail because of the lensing, but it is not modeled in detail here because it still lacks spectroscopic confirmation.

One anomaly is an apparent deficiency of z ≳ 6 galaxies seen by virtue of El Gordo’s magnification. More systematic searches and quantitative limits are needed. Spectroscopic confirmation by JWST of faint sources magnified by El Gordo may soon confirm the existence of more z ≳ 6 galaxies in this field.

Acknowledgments

The authors would like to thank Ignacio Gonzalez-Serrano, and Ignacio Negueruela for useful comments and discussion, and Wilfredo Villca for assisting with the Quechua translation and pronunciation. J.M.D. acknowledges the support of project PGC2018-101814-B-100 (MCIU/AEI/MINECO/FEDER, UE) Ministerio de Ciencia, Investigación y Universidades. This project was funded by the Agencia Estatal de Investigación, Unidad de Excelencia María de Maeztu, ref. MDM-2017-0765. AZ and AKM acknowledge support by Grant No. 2020750 from the United States-Israel Binational Science Foundation (BSF) and Grant No. 2109066 from the United States National Science Foundation (NSF), and by the Ministry of Science & Technology, Israel. RAW, SHC, and RAJ acknowledge support from NASA JWST Interdisciplinary Scientist grants NAG5-12460, NNX14AN10G and 80NSSC18K0200 from GSFC. EZ acknowledges funding from the Swedish National Space Agency. Work by CJC acknowledges support from the European Research Council (ERC) Advanced Investigator Grant EPOCHS (788113). LD acknowledges the research grant support from the Alfred P. Sloan Foundation (Award Number FG-2021-16495). 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. GBC acknowledges the Max Planck Society for financial support through the Max Planck Research Group for S. H. Suyu and the academic support from the German Centre for Cosmological Lensing. The authors would like to thank the RELICS team for making the reduced data set available to the community. The scientific results reported in this article are based in part on data obtained from the Chandra Data Archive (ivo://ADS/Sa.CXO#obs/12258; ivo://ADS/Sa.CXO#obs/14022; ivo://ADS/Sa.CXO#obs/14023). 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 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, 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. 2017, 2018); ProFit: https://github.com/ICRAR/ProFit (Robotham et al. 2018); SExtractor: SourceExtractor: https://www.astromatic.net/software/sextractor/ or https://sextractor.readthedocs.io/en/latest/ (Bertin & Arnouts 1996). We would also like to thank Harald Ebeling for making the code ASMOOTH (Ebeling et al. 2006) available. Facilities: Hubble and James Webb Space Telescope Mikulski Archive (https://archive.stsci.edu).

References

Appendix A: Compilation of arc positions

Table A.1 lists the complete sample of images used as constraints.³ Systems 1, 2, 5, 10, 23 (images a and b), 28, 29, and 30 were originally used in the lens model of Zitrin et al. (2013). Systems 25, 32, and 35 were first added by Cerny et al. (2018) to their lens model. Systems 11, 23 (image c), 33, 34, and 36 were first identified and used as additional constraints by Diego et al. (2020). Finally, systems 7 and 8 were identified and used as constraints by Caputi et al. (2021). Systems 1 to 23 are taken directly from Caminha et al. (2022) and use the same ID numbers, which are in column 1. Those marked with symbol † indicate new counterimage candidates found in JWST images in systems that were missing a third counterimage. RA and Dec positions (FK5) of the counterimages are given in columns 2 and 3 respectively. When available, spectroscopic redshifts from MUSE are shown in column 4. Column 5 lists the geometric redshifts predicted by the lensing model based on spectroscopic redshifts. Geo-z marked with ‡ (systems 24 and 53) are unconstrained or poorly constrained by this lens model. System 24 (El Anzuelo) is not included in the set of constraints for the spectroscopic model, and the member galaxies constraining this model are completely unconstrained making impossible any redshift prediction for this system. The assigned redshift corresponds to the photometric redshift derived by Cheng et al. (2023). The redshift of system 53 is unconstrained by the lens model, indicating a possible issue with this system, which is close to a spiral member galaxy. This galaxy was modeled assuming the same mass-to-light ratio as other member galaxies, most of which are ellipticals. This is probably not correct. We fixed the redshift to z = 4 based on the color and location of multiple images of this system. Column 6 lists two photometric redshifts derived respectively from the RELICS program (HST) and from the JWST photometry. Images with no photo-z are marked with −1. Column 7 lists the flux in the F200W band. A −1 in this column indicates that there is no detection in this band. Fluxes followed by the symbol ¶ indicate possible systematic errors in the flux estimate due to neighboring member galaxies or diffraction spikes from bright stars. Column 8 lists the magnification predicted by the spectroscopic, Full-r, and Full-a models respectively. The magnification is computed at the corresponding redshift of the system. Magnification values shown as 99.99 are calculated to be > 100, but such high values are unreliable. The magnification is at the exact position given in columns 2 and 3. Changes of a fraction of an arcsecond can result in big changes in the magnification. Column 9 lists the time delays predicted by the spectroscopic, Full-r, and Full-a models respectively. The time delay is expressed in years and is relative to the image that arrives first, so by construction is always positive. Finally, the column labeled Rank shows the quality of the system. Systems marked with rank A are the most reliable and were used to derive the spectroscopic model. Systems marked with rank B were used to derive (together with systems having rank A) the Full-r and Full-a models. Counterimages marked with C are less reliable because they cannot be confirmed based on morphology arguments, but they are still consistent with the spectroscopic lens model. Counterimages ranked C were not used to constrain any of the models. The last three rows with IDs CP1, CP2, and CP3 show the positions used as critical point constraints with their corresponding redshifts.

We have checked the MUSE data cube for spectroscopic confirmation of the new system candidates. Only counterimage 55a shows a hint of a line, possibly [O II] emission at z = 1.0714, which would rule out this candidate as strongly lensed. Another faint feature is observed in the same arc and corresponds to Hδ at the same redshift. System 55 is composed of two counterimages separated by only 075. Only one of these counterimages shows spectral features in MUSE. 55a is also bluer than 55b. We cannot rule out projection effects being responsible for the spectral features in 55a and bluer color, and as mentioned earlier, this arc could simply not be strongly lensed. Given the small separation between counterimages 55a and 55b, this arc has little constraining power, so even if it is proven not to be strongly lensed, it has small impact on our lens models.

Fluxes in table A1 were obtained from data reduced with Pipeline version 1.6.2 in early August 2022 using the context file jwst_0942.pmap_filters. At the time of finalizing this paper a new calibration became available that affects the photometric measurements. Since we relied on geometric redshifts for models Full-r and Full-a, our lens models are insensitive to this change, but fluxes listed in column 8 need to be multiplied by factors, 0.8885, 0.8285, 0.8986, 0.7998 for images falling in modules B1, B2, B3, and B4 respectively in order to reflect the new calibration. The correction factors are smaller at longer wavelengths: 0.9893, 1.0108, 1.0790, and 1.1128 for the F277W, F356W, F410M, and F444W filters respectively. These are the filters in which we detected Quyllur, and its corrected spectrum is ≈0.11 magnitudes brighter in F444W and ≈0.17 magnitudes fainter in F200W and would correspond to a slightly cooler star.

All Tables

All Figures

Fig. 1. Composite image of El Gordo obtained after combining HST and JWST bands (filters spanning ≈0.5 μm to ≈5 μm). Top panel: White circles mark the systems confirmed spectroscopically by MUSE. Bottom panel: Different color rendering that highlights fainter objects. Yellow circles mark new system candidates identified with the new JWST data. In each label, numbers identify multiple images of the same background source, and letters identify the individual images. The image orientation and scale are indicated in the lower left.

In the text

	Fig. 2. El Anzuelo galaxy in the northwestern part of the cluster. The galaxy has a photo-z of 3.58 (Cheng et al. 2023). The color-coded circles mark image pairs of three lensed sources. Each pair identifies a critical point between the images.
In the text

Fig. 3. Perturbers in the giant arc La Flaca. Scale and orientation are shown on the image. Color-coded circles mark multiple images of three background sources. La Flaca is extremely thin and one of the longest known arcs. It is located between the two main mass clumps in El Gordo. Small-mass perturbers create additional multiple images. These perturbers are marked with arrows. The approximate position of the critical curve can be estimated by extrapolating the ratio of distances between knots 2.1a–2.1e and knots 2.2a–2.2e. The point where one expects the critical curve is marked with the cyan circle having 1″ radius. The two cyan ellipses within this circle mark a possible double counterimage. If these two objects are counterimages of each other, the critical curve would be very close to the middle point between the two. The small cyan circle labeled A marks a possible caustic-crossing object.

In the text

	Fig. 4. JWST image of the central BCG with smoothed X-ray contours shown in green. Scale bar shows physical distance corresponding to 256. The white ellipse marks the position and beam shape (ATCA at 2.1 GHz) of the unresolved radio source U7 (Lindner et al. 2014). Blue knots in and above the circle agree in position with [O II] emission seen by MUSE.
In the text

Fig. 5. Multiresolution adaptive grid. The number density of grid points (shown as 891 red crosses) traces a smooth version of the solution in the spectroscopic lens model, increasing the resolution in the densest regions. Ten additional grid points are shown as small black crosses. Two of them mark the perturbers in the La Flaca arc with scales of 012 and 024, and eight more are in the inner region of the Anzuelo arc. Each one of these 8 grid points has a scale of 18.

In the text

	Fig. 6. Spectroscopic lens model. Contours represent the mass distribution of the smooth component of the lens model, that is, dark matter, diffuse baryons such as stars from the ICL, and X-ray-emitting plasma. The galaxies used to describe the compact contribution to the lens model are shown in gray. North is up and east is left, and the image is 2.5 arcmin across.
In the text

	Fig. 7. Diffuse mass component from the spectroscopic model vs. X-ray emission from Chandra. Color shows the convergence at z = 3 as indicated by the color bar. Black contours show the X-ray emission from Chandra. The image is oriented north up, east to the left, and a scale bar is shown at top left.
In the text

Fig. 8. Geometric redshifts compared with the input spectroscopic redshifts. This plot shows a sanity test where the input redshifts are properly recovered by the lens model. Families with discrepant redshifts are labeled. Black points show results of the spectroscopic model, and red points show results of a bootstrap analysis where a new model was derived without that particular system and that model used to predict the system’s redshift.

In the text

	Fig. 9. Comparison of JWST photometric redshifts with other redshifts. Black points represent spectroscopic redshifts, and red points represent geometric redshifts, both of which are shown on the x axis. Dashed line shows equality, and discrepant systems are labeled. Only systems having consistent photo-z estimates for at least two counterimages are shown.
In the text

	Fig. 10. Diffuse mass distribution and galaxy locations for the Full-r model. Contours represent the smooth component (dark matter and diffuse baryons such as stars from the ICL and X-ray-emitting plasma). The galaxies used to describe the compact contribution to the lens model are shown in gray.
In the text

	Fig. 11. Diffuse mass component from the Full-r model vs. X-ray emission from Chandra. Color shows the convergence at z = 3 as indicated by the color bar. Black contours show the X-ray emission from Chandra. The image is oriented north up, east to the left, and a scale bar is shown at top left.
In the text

Fig. 12. Small-scale deflectors. The giant, thin arc nicknamed La Flaca runs horizontally across the image. The heavy red, green, and blue curves are the critical curves of El Gordo at the redshift of the arc (zspec = 2.8254) for the spectroscopic, Full-r and Full-a models, respectively. The Full-a model (but not the other models) includes small deflectors at the two positions marked by the white and yellow arrows. Dotted blue curves show the critical curves of the deflectors. The deflector marked with a yellow arrow has a MUSE spectrum with z = 0.7842. With a Gaussian FWHM of 024, its mass is 2.7 × 1010 M⊙. The deflector marked with a white arrow is barely detected, and we have assumed it is at El Gordo’s redshift. With , its mass is 3.8 × 109 M⊙.

In the text

	Fig. 13. Mass vs. galaxies for the Full-a model. Contours represent the mass distribution of smooth component (dark matter and diffuse baryons such as stars from the ICL and X-ray emitting plasma). The galaxies used to describe the compact contribution to the lens model are shown in gray.
In the text

Fig. 14. Mass vs. X-ray for the model Full-a. The color plot shows the convergence (at z = 3) for the spectroscopic model. Only the dark matter component of the lens model is shown. X-ray emission from Chandra are shown as black contours. It is important to note the excess of mass near the X-ray local peak at the east. The south clump is now clearly more massive than the northern one but also more elongated along the direction of motion.

In the text

	Fig. 15. Comparison of the z = 2.8254 critical curves for the three models. The red curve corresponds to the spectroscopic model, the green curve is for model Full-r, and the blue curve is for model Full-a. The color image is made from the combination of HST+JWST filters RGB6 discussed in Sect. 2.
In the text

Fig. 16. Total El Gordo mass as a function of radius. The red curve shows values for the spectroscopic model, the green curve for the Full-r model, and the cyan curve for the Full-a model. The purple line shows the mass obtained derived from the RELICS data using the same algorithm but different constraints (Diego et al. 2020). The orange dot-dashed line shows the recent model of Caminha et al. (2022). The black dashed line shows values from a simulation (Molnar & Broadhurst 2015). All mass curves are based on a common center, and the masses are cylindrical, that is, they correspond to the mass within a cylinder of radius R projected along the line of sight. The light turquoise line shows the mass enclosed in a sphere with constant density equal to 200 times the critical density. While the spherical mass is not identical to the cylinder mass, for an NFW with concentration c = 7, the 3D mass computed up to the virial radius is only ≈4% smaller than the 2D mass. The vertical dotted line marks the maximum radius (70″) at which lensing constraints exist. The integrated masses will be biased low beyond this point.

In the text

Fig. 17. Possible caustic crossing event in system 23 at z = 2.1878. Top: Color image combining F115W, F200W, and F356W as blue, green, and red, respectively. The outer, yellow circles mark the positions of 23.4a and 23.4b, two counterimages of a source that brackets the position of the critical curve. The white, central circle marks a bright source (called “Quyllur”) that lies approximately at the midpoint. Bottom: color image with the RGB6 combination. Red, green, and blue curves show the critical curves for the spectroscopic, Full-r, and Full-a models respectively. The two panels have the same scale and orientation as indicated at the left.

In the text

	Fig. 18. Quyllur compared to SEDs of single stellar populations at El Gordo’s redshift z = 0.87. Blue lines show the unreddened models as labeled, and green lines show the models reddened with E(B − V) = 3.5 mag. Black points with error bars show the Quyllur photometry. For a comparison with the SED of individual stars see Fig. 21.
In the text

Fig. 19. Possible caustic crossing in system 1. Top: Color image combining F090W, F200W, and F356W as blue, green, and red, respectively. Bottom: Color image in the RGB6 pallette with critical curves from the spectroscopic, Full-r, and Full-a models shown as red, green, and blue, respectively. In both images, the two yellow circles mark images 1.3a and 1.3b, a pair of counterimages close to the expected position of the critical curve. The white circle contains a faint source that is barely resolved.

In the text

Fig. 20. Model magnification ratio versus measured F200W flux ratio for image pairs with accurate photometry. The flux ratio is defined as the counterimage with the largest flux divided by counterimage with the smallest flux. Solid lines mark equality, and dashed lines mark a factor of two deviation, as expected from uncertainty in the model magnifications. Points are color-coded for each model as shown in the legend, and two outliers are marked.

In the text

Fig. 21. SED of Quyllur compared to stellar models. The black dots show the observed aperture-corrected magnitudes. The blue, green, and orange curves show three models from Coelho (2014) redshifted to z = 2.1878 and magnified by a factor of 20 000. Model temperatures are shown in the legend. Squares show the expected flux in the JWST bands for a star with T ≈ 3500 K and luminosity 105.1 L⊙, similar to Betelgeuse. A star with a similar temperature but being a few times brighter, such as UY Scuti, would require magnification to be less than 10 000.

In the text