El Gordo as seen by JWST
Diego et al.

JWST's PEARLS: a new lens model for ACT-CL J0102-4915, “”, and the first red supergiant star at cosmological distances discovered by JWST

Jose M. Diego<EMAIL_ADDRESS>
0000-0001-9065-3926]Jose M. Diego
Instituto de Física de Cantabria (CSIC-UC). Avda. Los Castros s/n. 39005 Santander, Spain

0000-0002-7876-4321]Ashish K. Meena
Physics Department, Ben-Gurion University of the Negev, P.O. Box 653, Beer-Sheva 8410501, Israel

0000-0003-4875-6272]Nathan J. Adams
Jodrell Bank Centre for Astrophysics, Alan Turing Building,
University of Manchester, Oxford Road, Manchester M13 9PL, UK

0000-0002-5807-4411]Tom Broadhurst
Department of Theoretical Physics, University of the Basque
Country UPV-EHU, 48040 Bilbao, Spain
Donostia International Physics Center (DIPC), 20018 Donostia, The
Basque Country
IKERBASQUE, Basque Foundation for Science, Alameda Urquijo, 36-5
48008 Bilbao, Spain

0000-0003-2091-8946]Liang Dai
Department of Physics, 366 Physics North MC 7300, University of California, Berkeley, CA 94720, USA

0000-0001-7410-7669]Dan Coe
AURA for the European Space Agency (ESA), Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21210, USA

0000-0003-1625-8009]Brenda Frye
Steward Observatory, University of Arizona, 933 N Cherry Ave,
Tucson, AZ, 85721-0009, USA

0000-0003-3142-997X]Patrick Kelly
School of Physics and Astronomy, University of Minnesota, 116 Church Street SE, Minneapolis, MN 55455, USA

0000-0002-6610-2048]Anton M. Koekemoer
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21210, USA

0000-0002-2282-8795]Massimo Pascale
Department of Astronomy, University of California, 501 Campbell
Hall #3411, Berkeley, CA 94720, USA

0000-0002-9895-5758]S. P. Willner
Center for Astrophysics  Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

0000-0003-1096-2636]Erik Zackrisson
Observational Astrophysics, Department of Physics and Astronomy, Uppsala University, Box 516, SE-751 20 Uppsala, Sweden

0000-0002-0350-4488]Adi Zitrin
Physics Department, Ben-Gurion University of the Negev, P.O. Box
653, Beer-Sheva 8410501, Israel

0000-0002-7908-9284]Larry D. Bradley
Space Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA

0000-0001-6052-3274]Gabriel Caminha
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany   (ORCID: 0000-0001-6052-3274)

0000-0001-8183-1460]Karina Caputi
Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV Groningen, The Netherlands
Cosmic Dawn Center (DAWN), Copenhagen, Denmark

The first JWST data on the massive colliding cluster El Gordo confirm 23 known families of multiply lensed images and identify 8 new members of these families. Based on these families, which have been confirmed spectroscopically by MUSE, we derived an initial lens model. This model guided the identification of 37 additional families of multiply lensed galaxies in the  JWST data, among which 28 of them are entirely new system candidates and 9 are previously known. The initial lens model also helped determine geometric redshifts for the 37 additional systems. The geometric redshifts agree reasonably well with spectroscopic or photometric redshifts when those are available. The geometric redshifts enable two additional models that include all 60 families of multiply lensed galaxies spanning a redshift range $2 \lesssim z \lesssim 6$.  The derived distribution of dark matter confirms the double-peak configuration of the mass distribution of El Gordo found by earlier work with the southern and  northern clumps having similar masses. We confirm that El Gordo is the most massive known cluster at $z>0.8$ and has an estimated  virial mass close the maximum mass allowed by standard cosmological models.
The  JWST images also reveal the presence of small-mass perturbers  that produce small lensing distortions. The smallest of these  is consistent with being a dwarf galaxy at $z=0.87$ and has an estimated mass of $3.8\times10^9$ , making it the smallest substructure found at $z>0.5$. The JWST images also show several candidate caustic-crossing events. One of them is detected at high significance at the expected position of the critical curve and is likely a red supergiant star at  $z=2.1878$. This would be the first red supergiant found at cosmological distances. The cluster lensing should magnify background objects at $z\ga6$, making more of them visible than in blank fields of similar size, but there appears to be a deficiency of such objects.

§ 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 JWST/NIRCam as part of the PEARLS GTO project  Windhorst et al., 2022.
At $z=0.870$,  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] argued 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 towards each other. This interpretation is, however, challenged by lens models based on strong lensing,  which place most of the mass in the SE 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,  shows up in Planck maps as a massive cluster with a signal-to-noise ratio of $\approx 13$ [Planck Collaboration et al., 2016]. Some of the SZ signal may be due to increased pressure from the ongoing merger.

El Gordo Mass Estimates

Mass ($10^{15}$ M$_{\odot}$)

SZ     M_200m     1.89±0.15     [Williamson et al., 2011]

SZ     M_200m     2.16±0.32
[Menanteau et al., 2012]

SL     M_200m     2.3     [Zitrin et al., 2013]

WL     M_200c     3.13±0.56     [Jee et al., 2014]

SZ     M_500c     1.07±0.05     [Planck Collaboration et al., 2016]

WL     M_200c     1.11     [Schrabback et al., 2018]

SL     M(500 kpc)     1.1     [Cerny et al., 2018]

SL     M_200c     1.08±0.70     [Diego et al., 2020]

WL     M_200c     2.13     [Kim et al., 2021]

SL     M(1 Mpc)     1.84     [Caminha et al., 2022]

SL     M(500 kpc)     0.80–0.86     This paper

SL     M(R_vir)     2.09–2.24     This paper

aSZ = Sunyaev-Zeldovich, SL = strong lensing, WL = weak lensing.
b$M_{200m}$ is the mass within a sphere of radius $r_{200m}$, which is the radius above which the mass density drops below 200$\times$ the mean mass density of the Universe.  $M_{200c}$ is similar but measured within the radius where mass density drops below 200$\times$ the critical density, and $M_{500c}$ is the mass enclosed within the radius above which the density falls below 500 times the critical density. Other masses are measured within the specified physical radii. All $M_{200m}$ and $M_{200c}$ SL masses are obtained by extrapolating the profile obtained in the constrained region.

's mass has been estimated  through a variety of techniques as shown in Table <ref>.

These methods agree that
El Gordo is the most massive cluster in its redshift range with an estimated mass
above $10^{15}$ .

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 $\Lambda$CDM model <cit.>, which finds a total mass exceeding  the maximum mass expected at this redshift at ${\lesssim} 2 \times 10^{15}~M_{\odot}$ [Harrison & Coles, 2012]. A similar conclusion is reached from large N-body simulations. Using the very large 630 Gpc$^3$ N-body simulation Jubilee (based on a standard $\Lambda$CDM model), <cit.> found that the most massive cluster at $z=0.9$ is expected to have $M_{178_m} \approx 1.7\times 10^{15}~M_{\odot}$. ([Watson et al., 2014] defined the masses as $M_{178_m}$ rather than $M_{200_m}$ or $M_{200c}$. For an NFW profile, $M_{200_m} \approx 1.2M_{200c}$, and $M_{178_m}\approx 4\%$ times higher than $M_{200_m}$, hence this mass would correspond to  $M_{200_c} \approx 1.4\times 10^{15}M_{\odot}$.)
The observed masses are mostly above this maximum mass but can become smaller if one considers  Eddington bias [Waizmann et al., 2012].

Despite that, El Gordo is still uncomfortably close to the expected  $\Lambda$CDM limit.  El Gordo was found in a survey (the original ACT survey) that covered less than 2% the area of the sky. The low likelihood of finding such an extreme object in such a small area raises the  tension with the $\Lambda$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 Eddington bias, [Waizmann et al., 2012] found El Gordo not to be in tension with $\Lambda$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 in 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 towards 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/hydrodynamical simulations reproducing most of the observations of El Gordo [Molnar & Broadhurst, 2015, Zhang et al., 2015].
[Molnar & Broadhurst, 2018] demonstrated that the speed of the outgoing shocks can be very large (4000–5000 km s$^{-1}$) in a massive, merging cluster like El Gordo, therefore leaving the system before the first turnaround.

El Gordo is also unique in that it is a powerful lens at relatively high redshift.  One of the features that makes El Gordo an attractive target for lensing studies is the fact that for sources at 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-redshift known 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\gg1$ 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 use our free-form lensing reconstruction code WSLAP+ [Diego et al., 2005, Diego et al., 2007, Sendra et al., 2014, Diego et al., 2016, Diego et al., 2022], 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 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 <ref> describes the data and simulations used in this work. Section <ref> describes the lensing constraints, and Section <ref> 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 Section <ref> together with a comparison with earlier models in the literature.
Section <ref> discusses small-mass perturbers found near a giant arc nicknamed “La Flaca.”
Section <ref> 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 <ref> discusses the key results, and Section <ref> summarizes our conclusions.

We adopt a standard flat $\Lambda$CDM cosmological model with $\Omega_m=0.3$ and $h=0.7$. At the redshift of the lens, this cosmological model implies that 1 corresponds to 7.8 kpc.

§ JWST OBSERVATIONS AND ANCILLARY DATA

Composite image of El Gordo obtained after combining HST and JWST bands (filters spanning from $\approx 0.5$ micron to $\approx 5$ micron). Top panel: white circles mark the systems confirmed spectroscopically by MUSE.  Bottom panel: 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. Image orientation and scale are indicated on the lower left.

§.§ 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$\sigma$ magnitude limit is $\approx$29.9 in the long-wavelength filters. The best spatial resolution was obtained in the F150W filter with a measured FWHM of 0063. Further details are  in the PEARLS overview paper [Windhorst et al., 2022]. Figure <ref> shows an overview of the field.

El Anzuelo galaxy in the northwestern part of the cluster. The galaxy has a photo-$z$ of $3.54$ [Cheng, 2022]. The color-coded circles mark image pairs of three lensed sources.  Each pair identifies a critical point between the images.

JWST reveals details of this cluster with unprecedented clarity.  Galaxies barely detected or not detected at all in HST data are easily detected by JWST, and many are much brighter at wavelengths longer than 2  than at the longest HST wavelength (1.6 ).  Figure <ref> shows one such galaxy, which is heavily distorted by cluster-galaxy lensing into a fishhook shape. We will 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 ALMA sources EG-SMG 2 and EG-SMG 4 [Cheng, 2022], 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 $\mu$m. HST data do, however, identify a multiply-lensed blue feature.

Another  example illustrating the JWST data quality is shown in Figure <ref>, which shows the giant arc already observed in pre-JWST data. We refer to this galaxy as "La Flaca" ("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 <ref> 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 (Figure <ref>). 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 $5\farcs05/(1+1.66)=1\farcs9$ southwest from knot 2.2e as shown in Figure <ref>. (The distance between knots 2.1e and 2.2e (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.

§.§ 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 $\sim$0.4–1.6$\mu$m. The Reionization Lensing Cluster Survey <cit.> 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 . 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 $\rm RA=15.752995$, $\rm Dec=-49.281048$ and spectroscopically confirmed at $z=0.8749$ by [Caminha et al., 2022]) that falls on the edge of the JWST image.

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 with in 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.

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

\begin{align}
\nonumber
R&=\rm 0.3\times F277W({JWST})
+0.35\times F356W({JWST})\\
&\rm +0.15\times F410M({JWST})
+0.2\times F444W({JWST}),\\
G&=\rm 0.45\times F150W({JWST})
+0.55\times F200W({JWST}), \\
\nonumber
B&=\rm 0.1\times F606W({HST})
+0.1\times F625W({HST})\\
&\rm +0.3\times F090({JWST})
+0.4\times F115({JWST})\quad.
\end{align}

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

§.§ X-ray data

Chandra data from the ACIS instrument acquired in 2011–2012 (ObsID 12258, 14022, and 14023, PI.. J. Hughes) total $\approx$350 ks. The X-ray data were smoothed using the code asmooth [Ebeling et al., 2006].
Figure <ref> 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 [O2] and correlates spatially with the peak of the X-ray emission. [O2] line in MUSE data reveals 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 towards the BCG. The strong [O2] 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 a radio point source “U7” at one of the extremes of the filament. No obvious counterpart for U7 can be identified in JWST images, but when one accounts for the beam of the radio data, U7 coincides with the peak of the X-ray emission and the infalling [O2]-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.

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 radio point source U7  [Lindner et al., 2014].  Blue knots in and above the circle agree in position with [O2] emission seen by MUSE.

§ 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 <ref>.

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 <ref>.

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.

§ FREE-FORM LENS MODELING OF EL GORDO

§.§ Modeling method

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

\begin{equation}
\beta = \theta - \alpha(\theta,\Sigma(\theta)),
\label{eq_lens}
\end{equation}

where $\theta$ is the observed position of the lensed source, $\alpha$ is the deflection angle, $\Sigma(\theta)$ is the surface mass-density of the lens (El Gordo cluster in our case) at  position $\theta$, and $\beta$ 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 $\psi$:[Note however, that through observations one measures the reduced shear, $\gamma_r=\gamma/(1-\kappa)$ where $\kappa$ is the convergence.]

\begin{equation}
\psi(\theta) = \frac{4 G D_{l}D_{ls}}{c^2 D_{s}} \int  \Sigma(\theta')\ln(|\theta - \theta'|)\,d^2\theta' \quad,
\label{eq_psi}
\end{equation}

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 and the positions of the background sources in the source plane.

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

\begin{equation}
\vectTheta = \matrGamma \vectX\quad,
\label{eq_lens_system}
\end{equation}

where the observables (positions of strongly lensed galaxies, positions of critical points,  and weak lensing if available) are contained in the array $\vectTheta$, and the unknown surface mass density, re-scaling factors for the different layers discussed above, and true source positions are in the array $\vectX$. The matrix $\matrGamma$ is known and given by the position of the grid points and positions of the constraints. In our case, $\vectTheta$ contains the positions listed in Table <ref>  (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  $\vectX>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 pre-defined positions. Each mass center was assumed to represent a surface mass-distribution defined by a Gaussian with a pre-defined 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.)

For the compact component, we adopted  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.
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$ <cit.>. As discussed by [Caminha et al., 2022], although this is a small group, it has a non-negligible contribution to the lensing.

Spectroscopic lens model. Contours represent the mass distribution of the smooth component of the lens model, i.e., 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. The image is 2.5 arcminutes across.

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.

§.§ Spectroscopic lens model

We derived a first model (the “spectroscopic model”) using only 23 families

that have spectroscopic redshifts. These are identified as rank A  in Table <ref>.
The soft-component mass-centers were defined on a uniform 32$\times$32 grid.  This is equivalent to adopting a flat prior for the mass distribution.
The grid spacing was 4.6875, and the FWHM of each Gaussian mass center was taken to be 24.

The solution for the spectroscopic model is shown in Figure <ref>. 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.

Figure <ref>  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 (Section <ref>).

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.

§.§ 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 <ref>. The probability of a system to be at redshift $z$ is

\begin{equation}
P(z) = \exp(-V(z)/[2\sigma^2)]\quad,
\label{eq_Pz}
\end{equation}

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 re-scaled for each system to the corresponding redshift. The dispersion, $\sigma$, 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 result 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\ga4$.

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.

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 Figure <ref>. 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 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.

Comparison of 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 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, Chan et al., 2020].

Geometric redshifts are a useful alternative to photometric redshifts, especially for high-redshift candidates for which photometry may be poor or nonexistent (e.g., dropouts). Also, photometric redshifts may be unreliable for dusty galaxies at $z\approx4$ because their Balmer break may be misinterpreted as Ly-$\alpha$, placing them at much higher redshifts [Naidu et al., 2022, Zavala et al., 2022, Harikane et al., 2022].

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., 2022]. 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 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\le z< 25$.   Figure <ref> 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 Figure <ref> (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), $\approx$86% have LePhare redshifts with an error $\delta z/(1+z)<0.15$.

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.

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 <ref>), 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_{\rm spec}=4.3175$,  RELICS predicts $z_{\rm phot}\approx 4.3$ while JWST photometry alone predicts $z_{\rm phot}\approx 0.7$. System 10 at $z_{\rm spec}=4.3275$ is consistent with the RELICS $z_{\rm phot}\approx$ 4.2–4.6, and in better agreement than JWST-based $z_{\rm phot}\approx 0.7$. This is almost identical to System 12 ($z_{\rm spec}=4.7042$) where HST  $z_{\rm phot}\approx$ 4.6–5.3 while JWST  $z_{\rm phot}\approx 0.7$.
System 21 at $z_{\rm spec}=5.5811$ is predicted to be at $z_{\rm phot}\approx 6$ according to HST while JWST predicts $z_{\rm phot}\approx1.3$.
[Frye, 2022] will provide additional photometric redshifts that combine HST and JWST photometry.
The important lesson from Figures <ref> and <ref> is that geo-$z$ can be reliable estimates of the redshift.

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.

§.§ 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., 2022, Harikane et al., 2022].
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 Figure <ref> (bottom).  The “Full-r model” maintains the $32\times32$ 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 <ref>.

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 .

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. 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.

The Full-r model  solution  is shown in Figure <ref>. 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 $\le$100 kpc and infalling velocities are $\ga$2500 km s$^{-1}$). 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 (Figure <ref>) 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.

Comparison of the critical curves for the three models. The critical curves are computed at $z=2.8254$. The red curve corresponds to the spectroscopic model, the green curve is for model Full-r, and the blue curve is fr model Full-a. The color image is made from the combination of HST+JWST filters RGB6 discussed in section <ref>.

§.§ 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 Figure <ref>.
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 (Figure <ref>) and  eight additional Gaussians with FWHM 18 around the group lensing the galaxy El Anzuelo.

The Full-a model reveals even more details as shown in Figures <ref> and <ref>. The mass distribution has a more concentrated peak around the BCG in the SE clump, while the NW 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 [O2] filament. This suggests that the contribution of the dense filament to the projected mass must be relatively small.

§ EL GORDO MASS ESTIMATE

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.  The light turquoise curve shows the total mass enclosed in a sphere with constant density equal to the critical density. The vertical dotted line marks the maximum radius (70) at which lensing constraints exist. The  integrated masses will be biased low beyond this point.

El Gordo Clump Mass Estimates

$M(\rm 300~kpc)$

SE clump:

[Caminha et al., 2022]    2.29$^a$

Spectroscopic     3.46

Full-r     3.43

Full-a     3.29

NW clump:

[Caminha et al., 2022]    2.19$^a$

Spectroscopic     3.85

Full-r     3.79

Full-a     3.67

aMasses in [Caminha et al., 2022] exclude the contribution from member galaxies

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 $\rm RA=15\fdg7210573$, $\rm Dec=-49\fdg2528437$. Masses of each clump for the various models are given in Table <ref>. 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.

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 ($z_{\rm spec}=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\times10^{10}$ . The deflector marked with a white arrow is barely detected, and we have assumed it is at El Gordo's redshift. With $\rm FWHM = 0\farcs12$, its mass is  $3.8\times10^{9}$ .

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}\lesssim 2\times 10^{15}$  [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  $\rm RA=15.7313639$, $\rm Dec=-49\fdg2605609$.
Figure <ref> 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\times10^{14}$ .

These are higher than the [Diego et al., 2020] model based on RELICS HST data (M(500 kpc)$=6.6\times10^{14}$ ) 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~\rm kpc})=7.98\times 10^{14}$ , 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}\approx 500$ kpc, the limit of the strong-lensing constraints, and    $M(500~{\rm kpc})=7.20\times 10^{14}$ , 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_{\rm vir})=1.88\times 10^{15}$ .  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.

§ SMALL SCALE SUBSTRUCTURES IN LA FLACA

As discussed in Section <ref> and shown in Figure <ref>, 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.
Figure <ref> 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 measuring masses for even smaller substructures, potentially constraining models of dark matter, as  [Diego et al., 2022] recently did.

§ 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, Oguri et al., 2018, Venumadhav et al., 2017, Dai et al., 2018, Dai & Miralda-Escudé, 2020, Dai et al., 2020, Diego et al., 2022, Meena et al., 2022]. After the discovery of Icarus, other examples followed [Chen et al., 2019, Kaurov et al., 2019, Diego et al., 2022, Chen et al., 2022], culminating with the recent discovery of Earendel at $z\approx 6.2$ [Welch et al., 2022, Welch et al., 2022].

JWST is expected to see further 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 Figure <ref>.  Given its proximity to the midpoint between the pair of knots

23.4a and 23.4b (yellow circles in Figure <ref>), this source (marked with a  white circle in Figure <ref>) 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., 2022, Welch et al., 2022], where those sources must also form unresolved pairs of counterimages.

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.[Pronounced Koijur in English. Quechua was spoken by the Incas before the arrival of  Europeans, and it is still spoken by millions of people in different parts of South America. Similar to the disappearance of the Incan empire, Quyllur must have vanished eons ago. But similarly to the Quechua language that is still alive, Quyllur's light still traverses the Universe].

Figure <ref> shows Quyllur's SED. The SED shows a rapid decline below 3  consistent with a surface temperature $T\approx 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 $\mu>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  magnitude is 25.5, which corresponds to $>$34.5 unmagnified.   At a distance modulus of $\sim$45, and depending on the magnification, the absolute magnitude at 1.1  rest wavelength would be ${\sim}-$10.5, consistent with an M2–M4 supergiant.

The required magnification limits the source size and therefore the possible source types.
For $\mu=4000$, the  distance to the caustic must be less than 7.5 microarcseconds or $\approx$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, Davies et al., 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 $\lesssim$0.06 pc is significantly smaller than the $\ga$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 shortwards of 1.5  (0.47  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 render a color index E(G-J)$\approx 2.5$ (this color index corresponds to the difference between the F410M and F150W filters respectively), while for Quyllur we find a much redder index E(G-J)$\approx 5$. This rules out significant reddening affecting the portion of the galaxy surrounding Quyllur, and being multiply imaged, since only Quyllur shows such extreme color variation.

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.

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 <ref>), 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  $\mu \ga 30\arcsec/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 SNe, and not its counterimage, the JWST observation must have taken place  within the first days of the event, which is very unlikely.

Quyllur compared to SED fits to a single stellar population 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 Figure <ref>.

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 Figure <ref>, Quyllur can be consistent with either a 6 Gyr stellar population with stellar mass $M_* \sim 4\times 10^9$   or a 100 Myr stellar population with $M_*\sim 2 \times 10^8$. Either would be of substantial mass but would have to be more compact than $\sim$0.5 kpc, i.e., a dwarf spheroidal (dSph) galaxy. However, either scenario requires dust reddening $E(B-V)=3.5$ (assuming Milky Way like reddening curve with $R_V=3.1$). Without observational evidences 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 (Figure <ref>). 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.

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

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.

§ DISCUSSION

In contrast to the lensing cluster WHL J013725.2+140341 [Welch et al., 2022, Welch et al., 2022], 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\approx 4.3$ <cit.>, 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., 2022], resulting in a  more concentrated caustic region.  These characteristics make El Gordo a good target for lensing studies.

§.§ Possible tension with $\Lambda$CDM

The highest mass estimate for El Gordo,
$M_{\rm vir}\approx 3\times 10^{15}$   <cit.> has significant tension with $\Lambda$CDM. Our $M({\rm vir})\approx 2\times10^{15}$ 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, Ezquiaga et al., 2022] has used weak-lensing limits to illustrate the tension of this cluster with $\Lambda$CDM. At the core of this debate is the true mass of El Gordo.

More recent work (Table <ref>) has tended to give lower masses, but the mass of El Gordo remains uncertain within a factor two. 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. 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.

§.§ 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 <ref> 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.)

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 $\approx$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.

Model magnification ratio versus measured F200W flux ratio for image pairs with more 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.

§.§ Quyllur
Very luminous stars are rarely found isolated and normally live near stellar groups.  The nearby, multiply-lensed pair 23.4a–23.4b is estimated to be at a distance of 10–15 pc from Quyllur, depending on the magnification. This bluer source could be a star-forming region, and Quyllur would be in its outskirts.
Assuming Quyllur is a red supergiant, the magnification is likely significantly higher than our lower limit $\approx$4000.
Betelgeuse ($d=197$ pc with luminosity  $1.26\times10^5$ )
has a flux density of
$\approx 3\times 10^{-9}$ ergs s$^{-1}$ cm$^{-2}$ Å$^{-1}$ at 6500 Å [Levesque & Massey, 2020].
This wavelength would redshift into the F200W JWST filter, for which we find $M_{AB}\approx 28$.  At the distance of Quyllur ($d=1.76\times10^{10}$ pc) and without magnification, Betelgeuse would have an apparent magnitude of $\approx 39.75-2.5*{\rm log}_{10}(1+z)=38.5$ in  F200W (ignoring precise color- or K-correction and extinction).
This implies the magnification $\approx$15 000 in order to match Quyllur's apparent magnitude in F200W.

Red hypergiants are rare but have luminosities that can exceed the luminosity of Betelgeuse by a factor of $\approx$ 3–5 (for example UY Scuti and Stephenson 2-18) but still below below the observational limit  [Humphreys & Davidson, 1979]. At these luminosities, the required magnification drops to $\approx$4000–7000 (Figure <ref>).

At magnification factors of several thousand, microlenses are expected to introduce temporary distortions in the flux even for relatively low surface mass density of microlenses ($\Sigma \la 10$  pc$^{-2}$). Given 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 assume the stellar mass contributes $\approx$0.2%–1% to the total mass. At the position of Quyllur, the convergence from our lens models  $\kappa=0.63$–0.66, and the critical surface mass-density for the redshifts of El Gordo and Quyllur is 2300  pc$^{-2}$. This would make the surface mass-density of microlenses  between 3 and  15  pc$^{-2}$. Combined with the model magnification of a few thousand, this implies an effective surface mass density, $\Sigma_{\rm eff}=\mu\times\Sigma$, 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 $\approx 2$ when observing at intervals separated by several months [Welch et al., 2022].  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.

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 colder 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 tip of the red giant branch is commonly used in our local Universe as part of the distance ladder.

Extending this ladder up to or beyond $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 allow to constrain the magnification and open the door to using these stars as standard candles up to $z\approx 1$ and possibly beyond.

Although not discussed in detail is worth mentioning other exotic objects that show a red SED and could in principle mimic the observed photometry of Quyllur,  including luminous transients that can temporarily reach ${\sim}10^6$  and can be linked to proto-planetary nebulae [Prieto et al., 2009].

Finally we can not rule out compact luminous objects such as accretion disks around black holes. 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.

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}\approx 3500$ K and luminosity $10^{5.1}$ ,  similar to Betelgeuse. A star with similar temperature but a few times brighter, such as UY Scuti, would require magnification less than 10 000.

Regarding the probability of observing a star like Quyllur, one interesting aspect to consider is the large volume accessible through lensing, $V\approx 1.5E11$ Mpc$^3$ between $z=2$ and $z=2.3$ compared with a volume $V\approx 3.3E4$ up to $\approx 20$ Mpc, which is the largest distance up to which we have observed the brightest stars in our local universe. This implies through lensing we can probe a volume that is a factor ${\approx} 5 \times 10^6$ times larger at $2<z<2.3$. A larger volume and at the lower metallicity expected at $z\approx 2$, one should expect to see even brighter stars,  and $\mu \approx 1000$ may suffice to observe them with JWST. Only a small fraction of the stars at $z>2$ can be magnified above $\mu\approx 1000$ because $P(\mu>1000) \approx 10^{-8}$ [Diego, 2018].

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^5$ and that these are $\sim$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.

§.§ 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, 2022]. 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 like El Gordo amplify faint objects, the probability of observing a  $z>6$ galaxy is non-negligible.
Average magnifications $\mu\ga10$ 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 $\rm RA=15.730584,\  Dec=-49.271130$. It is  stretched to a length of $\sim$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, 2022].

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 <ref>.) 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 $\mu> 30$ of $\approx$0.015 arcmin$^2$. 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 <cit.>.

The UV luminosity function [Ishigaki et al., 2018] at $z=6$, derived from galaxies behind the Hubble Frontier Fields, reaches  $M_{\rm UV}\approx -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 $\approx 30$, that is, within reach of JWST.

The volume between redshifts $5.5<z<6.5$ is $V\approx360$ Gpc$^3$ or 2400 Mpc$3$ arcmin${-2}$.
The volumetric number density of all $5.5<z<6.5$ galaxies with $M_{UV}<-16$ is $N \approx 0.1$ galaxies per Mpc$^{3}$. Hence we expect a surface number density, $N(z=6)=240\, {\rm gal}\, {\rm arcmin}^{-2}$ 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 $\approx$ 0.015 arcmin$^2$ behind El Gordo to  be magnified by at least a factor 30, with one counterimage having magnification at least 12. That is, we expect an average of $N \approx 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 $\mu$ scales as $\mu^{-2}$, reducing the number of detectable galaxies. Also, the galaxy number density is smaller at higher redshift, so we expect many fewer galaxies at $z\gg6$ than at $z\approx6$.

A deficit of strongly magnified galaxies above $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.

§ 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 $\Lambda$CDM.

In the new JWST data we have identified 37 lensed-system candidates (28 of them 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 $\sim$2 Mpc gives a mass of about $2.1\times10^{15}$ ,
close to the limit for standard $\Lambda$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 the 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\times10^9$ , and that 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 modelled in detail here because it still lacks spectroscopic confirmation.

One anomaly is an apparent deficiency of $z\ga6$ galaxies seen by virtue of El Gordo's magnification.  More systematic searches and quantitative limits are needed.

§ ACKNOWLEDGEMENTS

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.
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 acknowledge 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.

Astropy: <http://www.astropy.org> [Astropy Collaboration et al., 2013, Astropy Collaboration et al., 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, Robotham et al., 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.

Hubble and James Webb Space Telescope Mikulski Archive

§ COMPILATION OF ARC POSITIONS

Table <ref> 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 in [Diego et al., 2020].
Finally, systems 7 and 8 were identified and used as constraints in [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 the spectroscopic 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, 2022]. 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 is modeled together with other members galaxies (ellipticals) and assuming the same light-to-mass ratio which is probably not fully correct. We fix 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 neighbouring 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 as $>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 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 constrints 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 OII emission at z=1.0714 which would rule out this candidate as strongly lensed. Another faint feature is observed in the same arc and corresponding to H$\delta$ at the same redshift. System 55 is composed of two counterimages separated by only 0.75". Only one of these counterimages shows spectral features in MUSE. 55a is also bluer than 55b. We can not rule out projection effects for 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 very limited constraining power so even if its finally proven not to be strongly lensed, it has a small impact on our lens models.

Fluxes in table A1 are 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 change in the calibration. The correction factors are smaller at longer wavelengths, with 0.9893,  1.0108,  1.0790, and  1.1128 for the F277W, F356W, F410M, and F444W filters respectively. Since these are the filters in which we detect Quyllur, its corrected  spectrum would be a bit redder ($\approx 0.11$ magnitudes brighter in F444W and $\approx 0.17$ magnitudes fainter in F200W) and would correspond to a slightly cooler star.

[Adams et al., 2022]
Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2022, arXiv
e-prints, arXiv:2207.11217.

[Alexander et al., 2009]
Alexander, M. J., Kobulnicky, H. A., Clemens, D. P., et al. 2009, ,
137, 4824, 10.1088/0004-6256/137/6/4824

[Arnouts et al., 1999]
Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, , 310,
540, 10.1046/j.1365-8711.1999.02978.x

[Arnouts & Ilbert, 2011]
Arnouts, S., & Ilbert, O. 2011, LePHARE: Photometric Analysis for
Redshift Estimate, Astrophysics Source Code Library, record ascl:1108.009.

[Asencio et al., 2021]
Asencio, E., Banik, I., & Kroupa, P. 2021, , 500, 5249,

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
, 558, A33, 10.1051/0004-6361/201322068

[Astropy Collaboration et al., 2018]
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al.
2018, , 156, 123, 10.3847/1538-3881/aabc4f

[Benítez, 2000]
Benítez, N. 2000, , 536, 571, 10.1086/308947

[Bertin & Arnouts, 1996]
Bertin, E., & Arnouts, S. 1996, , 117, 393,

[Bhatawdekar, 2022]
Bhatawdekar, R. A. . 2022, in preparation

[Boyer et al., 2022]
Boyer, M. L., Anderson, J., Gennaro, M., et al. 2022, arXiv e-prints,

[Bradley et al., 2020]
Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils:
1.0.0, 1.0.0,  Zenodo, 10.5281/zenodo.4044744

[Bruzual & Charlot, 2003]
Bruzual, G., & Charlot, S. 2003, , 344, 1000,

[Calzetti et al., 2000]
Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, , 533, 682,

[Caminha et al., 2022]
Caminha, G. B., Grillo, C., Rosati, P., et al. 2022, arXiv e-prints,

[Caputi et al., 2021]
Caputi, K. I., Caminha, G. B., Fujimoto, S., et al. 2021, , 908,
146, 10.3847/1538-4357/abd4d0

[Cerny et al., 2018]
Cerny, C., Sharon, K., Andrade-Santos, F., et al. 2018, , 859, 159,

[Chabrier et al., 2000]
Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ,
542, 464, 10.1086/309513

[Chan et al., 2017]
Chan, B. M. Y., Broadhurst, T., Lim, J., et al. 2017, , 835, 44,

[Chan et al., 2020]
—. 2020, , 888, 35, 10.3847/1538-4357/ab44a4

[Chen et al., 2019]
Chen, W., Kelly, P. L., Diego, J. M., et al. 2019, , 881, 8,

[Chen et al., 2022]
Chen, W., Kelly, P. L., Treu, T., et al. 2022, arXiv e-prints,

[Cheng, 2022]
Cheng, C. . 2022, in preparation

[Coe et al., 2019]
Coe, D., Salmon, B., Bradac, M., et al. 2019, arXiv e-prints.

[Coelho, 2014]
Coelho, P. R. T. 2014, , 440, 1027, 10.1093/mnras/stu365

[Dai, 2021]
Dai, L. 2021, , 501, 5538, 10.1093/mnras/stab017

[Dai & Miralda-Escudé, 2020]
Dai, L., & Miralda-Escudé, J. 2020, , 159, 49,

[Dai et al., 2018]
Dai, L., Venumadhav, T., Kaurov, A. A., & Miralda-Escud, J. 2018,
, 867, 24, 10.3847/1538-4357/aae478

[Dai et al., 2020]
Dai, L., Kaurov, A. A., Sharon, K., et al. 2020, , 495, 3192,

[Davies et al., 2007]
Davies, B., Figer, D. F., Kudritzki, R.-P., et al. 2007, , 671,
781, 10.1086/522224

[Davies et al., 2008]
Davies, B., Figer, D. F., Law, C. J., et al. 2008, , 676, 1016,

[Diego, 2018]
Diego, J. M. 2018, arXiv e-prints.

[Diego et al., 2022]
Diego, J. M., Pascale, M., Kavanagh, B. J., et al. 2022, arXiv
e-prints, arXiv:2203.08158.

[Diego et al., 2005]
Diego, J. M., Protopapas, P., Sandvik, H. B., & Tegmark, M. 2005a,
, 360, 477, 10.1111/j.1365-2966.2005.09021.x

[Diego et al., 2007]
Diego, J. M., Tegmark, M., Protopapas, P., & Sandvik, H. B. 2007,
, 375, 958, 10.1111/j.1365-2966.2007.11380.x

[Diego et al., 2016]
Diego, J. M., Broadhurst, T., Chen, C., et al. 2016, , 456, 356,

[Diego et al., 2018]
Diego, J. M., Kaiser, N., Broadhurst, T., et al. 2018, , 857, 25,

[Diego et al., 2020]
Diego, J. M., Molnar, S. M., Cerny, C., et al. 2020, , 904, 106,

[Ebeling et al., 2006]
Ebeling, H., White, D. A., & Rangarajan, F. V. N. 2006, , 368, 65,

[Ezquiaga et al., 2022]
Ezquiaga, J. M., García-Bellido, J., & Vennin, V. 2022, arXiv
e-prints, arXiv:2207.06317.

[Froebrich & Scholz, 2013]
Froebrich, D., & Scholz, A. 2013, , 436, 1116,

[Frye, 2022]
Frye, B. . 2022, in preparation

[Harikane et al., 2022]
Harikane, Y., Ouchi, M., Oguri, M., et al. 2022, arXiv e-prints,

[Harrison & Coles, 2012]
Harrison, I., & Coles, P. 2012, , 421, L19,

[Humphreys & Davidson, 1979]
Humphreys, R. M., & Davidson, K. 1979, , 232, 409,

[Ishigaki et al., 2018]
Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, , 854, 73,

[Jee et al., 2014]
Jee, M. J., Hughes, J. P., Menanteau, F., et al. 2014, , 785, 20,

[Kaurov et al., 2019]
Kaurov, A. A., Dai, L., Venumadhav, T., Miralda-Escudé, J., &
Frye, B. 2019, , 880, 58, 10.3847/1538-4357/ab2888

[Kelly et al., 2018]
Kelly, P. L., Diego, J. M., Rodney, S., et al. 2018, Nature Astronomy,
2, 334, 10.1038/s41550-018-0430-3

[Kim et al., 2021]
Kim, J., Jee, M. J., Hughes, J. P., et al. 2021, , 923, 101,

[Landsman, 1993]
Landsman, W. B. 1993, in Astronomical Society of the Pacific Conference
Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed.
R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 246

[Leung et al., 2018]
Leung, E., Broadhurst, T., Lim, J., et al. 2018, , 862, 156,

[Levesque & Massey, 2020]
Levesque, E. M., & Massey, P. 2020, , 891, L37,

[Lindner et al., 2014]
Lindner, R. R., Baker, A. J., Hughes, J. P., et al. 2014, , 786,
49, 10.1088/0004-637X/786/1/49

[Madau, 1995]
Madau, P. 1995, , 441, 18, 10.1086/175332

[Meena et al., 2022]
Meena, A. K., Arad, O., & Zitrin, A. 2022, , 514, 2545,

[Menanteau et al., 2012]
Menanteau, F., Hughes, J. P., Sifón, C., et al. 2012, , 748, 7,

[Meneghetti et al., 2017]
Meneghetti, M., Natarajan, P., Coe, D., et al. 2017, , 472, 3177,

[Miralda-Escude, 1991]
Miralda-Escude, J. 1991, , 379, 94, 10.1086/170486

[Molnar & Broadhurst, 2015]
Molnar, S. M., & Broadhurst, T. 2015, , 800, 37,

[Molnar & Broadhurst, 2018]
—. 2018, , 862, 112, 10.3847/1538-4357/aad04c

[Naidu et al., 2022]
Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022, arXiv e-prints,

[Ng et al., 2015]
Ng, K. Y., Dawson, W. A., Wittman, D., et al. 2015, , 453, 1531,

[Oguri et al., 2018]
Oguri, M., Diego, J. M., Kaiser, N., Kelly, P. L., & Broadhurst, T.
2018, , 97, 023518, 10.1103/PhysRevD.97.023518

[Pascale et al., 2022]
Pascale, M., Frye, B. L., Diego, J., et al. 2022, arXiv e-prints,

[Planck Collaboration et al., 2016]
Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, ,
594, A27, 10.1051/0004-6361/201525823

[Prieto et al., 2009]
Prieto, J. L., Sellgren, K., Thompson, T. A., & Kochanek, C. S. 2009,
, 705, 1425, 10.1088/0004-637X/705/2/1425

[Robotham et al., 2018]
Robotham, A. S. G., Davies, L. J. M., Driver, S. P., et al. 2018,
, 476, 3137, 10.1093/mnras/sty440

[Robotham et al., 2017]
Robotham, A. S. G., Taranu, D. S., Tobar, R., Moffett, A., & Driver,
S. P. 2017, , 466, 1513, 10.1093/mnras/stw3039

[Schrabback et al., 2018]
Schrabback, T., Applegate, D., Dietrich, J. P., et al. 2018, ,
474, 2635, 10.1093/mnras/stx2666

[Sendra et al., 2014]
Sendra, I., Diego, J. M., Broadhurst, T., & Lazkoz, R. 2014, ,
437, 2642, 10.1093/mnras/stt2076

[Vega-Ferrero et al., 2019]
Vega-Ferrero, J., Diego, J. M., & Bernstein, G. M. 2019, ,

[Venumadhav et al., 2017]
Venumadhav, T., Dai, L., & Miralda-Escudé, J. 2017, , 850, 49,

[Waizmann et al., 2012]
Waizmann, J. C., Ettori, S., & Moscardini, L. 2012, , 420, 1754,

[Watson et al., 2014]
Watson, W. A., Iliev, I. T., Diego, J. M., et al. 2014, , 437,
3776, 10.1093/mnras/stt2173

[Welch et al., 2022]
Welch, B., Coe, D., Diego, J. M., et al. 2022a, , 603,
815, 10.1038/s41586-022-04449-y

[Welch et al., 2022]
Welch, B., Coe, D., Zackrisson, E., et al. 2022b, arXiv
e-prints, arXiv:2208.09007.

[Williamson et al., 2011]
Williamson, R., Benson, B. A., High, F. W., et al. 2011, , 738,
139, 10.1088/0004-637X/738/2/139

[Windhorst et al., 2018]
Windhorst, R. A., Timmes, F. X., Wyithe, J. S. B., et al. 2018, ,
234, 41, 10.3847/1538-4365/aaa760

[Windhorst et al., 2022]
Windhorst, R. A., Cohen, S. H., Jansen, R. A., et al. 2022, arXiv
e-prints, arXiv:2209.04119.

[Zavala et al., 2022]
Zavala, J. A., Buat, V., Casey, C. M., et al. 2022, arXiv e-prints,

[Zhang et al., 2015]
Zhang, C., Yu, Q., & Lu, Y. 2015, , 813, 129,

[Zitrin et al., 2013]
Zitrin, A., Menanteau, F., Hughes, J. P., et al. 2013, , 770, L15,

[Zitrin et al., 2015]
Zitrin, A., Fabris, A., Merten, J., et al. 2015, , 801, 44,

[Adams et al., 2022]
Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2022, arXiv
e-prints, arXiv:2207.11217.

[Alexander et al., 2009]
Alexander, M. J., Kobulnicky, H. A., Clemens, D. P., et al. 2009, ,
137, 4824, 10.1088/0004-6256/137/6/4824

[Arnouts et al., 1999]
Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, , 310,
540, 10.1046/j.1365-8711.1999.02978.x

[Arnouts & Ilbert, 2011]
Arnouts, S., & Ilbert, O. 2011, LePHARE: Photometric Analysis for
Redshift Estimate, Astrophysics Source Code Library, record ascl:1108.009.

[Asencio et al., 2021]
Asencio, E., Banik, I., & Kroupa, P. 2021, , 500, 5249,

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
, 558, A33, 10.1051/0004-6361/201322068

[Astropy Collaboration et al., 2018]
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al.
2018, , 156, 123, 10.3847/1538-3881/aabc4f

[Benítez, 2000]
Benítez, N. 2000, , 536, 571, 10.1086/308947

[Bertin & Arnouts, 1996]
Bertin, E., & Arnouts, S. 1996, , 117, 393,

[Bhatawdekar, 2022]
Bhatawdekar, R. A. . 2022, in preparation

[Boyer et al., 2022]
Boyer, M. L., Anderson, J., Gennaro, M., et al. 2022, arXiv e-prints,

[Bradley et al., 2020]
Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils:
1.0.0, 1.0.0,  Zenodo, 10.5281/zenodo.4044744

[Bruzual & Charlot, 2003]
Bruzual, G., & Charlot, S. 2003, , 344, 1000,

[Calzetti et al., 2000]
Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, , 533, 682,

[Caminha et al., 2022]
Caminha, G. B., Grillo, C., Rosati, P., et al. 2022, arXiv e-prints,

[Caputi et al., 2021]
Caputi, K. I., Caminha, G. B., Fujimoto, S., et al. 2021, , 908,
146, 10.3847/1538-4357/abd4d0

[Cerny et al., 2018]
Cerny, C., Sharon, K., Andrade-Santos, F., et al. 2018, , 859, 159,

[Chabrier et al., 2000]
Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ,
542, 464, 10.1086/309513

[Chan et al., 2017]
Chan, B. M. Y., Broadhurst, T., Lim, J., et al. 2017, , 835, 44,

[Chan et al., 2020]
—. 2020, , 888, 35, 10.3847/1538-4357/ab44a4

[Chen et al., 2019]
Chen, W., Kelly, P. L., Diego, J. M., et al. 2019, , 881, 8,

[Chen et al., 2022]
Chen, W., Kelly, P. L., Treu, T., et al. 2022, arXiv e-prints,

[Cheng, 2022]
Cheng, C. . 2022, in preparation

[Coe et al., 2019]
Coe, D., Salmon, B., Bradac, M., et al. 2019, arXiv e-prints.

[Coelho, 2014]
Coelho, P. R. T. 2014, , 440, 1027, 10.1093/mnras/stu365

[Dai, 2021]
Dai, L. 2021, , 501, 5538, 10.1093/mnras/stab017

[Dai & Miralda-Escudé, 2020]
Dai, L., & Miralda-Escudé, J. 2020, , 159, 49,

[Dai et al., 2018]
Dai, L., Venumadhav, T., Kaurov, A. A., & Miralda-Escud, J. 2018,
, 867, 24, 10.3847/1538-4357/aae478

[Dai et al., 2020]
Dai, L., Kaurov, A. A., Sharon, K., et al. 2020, , 495, 3192,

[Davies et al., 2007]
Davies, B., Figer, D. F., Kudritzki, R.-P., et al. 2007, , 671,
781, 10.1086/522224

[Davies et al., 2008]
Davies, B., Figer, D. F., Law, C. J., et al. 2008, , 676, 1016,

[Diego, 2018]
Diego, J. M. 2018, arXiv e-prints.

[Diego et al., 2022]
Diego, J. M., Pascale, M., Kavanagh, B. J., et al. 2022, arXiv
e-prints, arXiv:2203.08158.

[Diego et al., 2005]
Diego, J. M., Protopapas, P., Sandvik, H. B., & Tegmark, M. 2005a,
, 360, 477, 10.1111/j.1365-2966.2005.09021.x

[Diego et al., 2007]
Diego, J. M., Tegmark, M., Protopapas, P., & Sandvik, H. B. 2007,
, 375, 958, 10.1111/j.1365-2966.2007.11380.x

[Diego et al., 2016]
Diego, J. M., Broadhurst, T., Chen, C., et al. 2016, , 456, 356,

[Diego et al., 2018]
Diego, J. M., Kaiser, N., Broadhurst, T., et al. 2018, , 857, 25,

[Diego et al., 2020]
Diego, J. M., Molnar, S. M., Cerny, C., et al. 2020, , 904, 106,

[Ebeling et al., 2006]
Ebeling, H., White, D. A., & Rangarajan, F. V. N. 2006, , 368, 65,

[Ezquiaga et al., 2022]
Ezquiaga, J. M., García-Bellido, J., & Vennin, V. 2022, arXiv
e-prints, arXiv:2207.06317.

[Froebrich & Scholz, 2013]
Froebrich, D., & Scholz, A. 2013, , 436, 1116,

[Frye, 2022]
Frye, B. . 2022, in preparation

[Harikane et al., 2022]
Harikane, Y., Ouchi, M., Oguri, M., et al. 2022, arXiv e-prints,

[Harrison & Coles, 2012]
Harrison, I., & Coles, P. 2012, , 421, L19,

[Humphreys & Davidson, 1979]
Humphreys, R. M., & Davidson, K. 1979, , 232, 409,

[Ishigaki et al., 2018]
Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, , 854, 73,

[Jee et al., 2014]
Jee, M. J., Hughes, J. P., Menanteau, F., et al. 2014, , 785, 20,

[Kaurov et al., 2019]
Kaurov, A. A., Dai, L., Venumadhav, T., Miralda-Escudé, J., &
Frye, B. 2019, , 880, 58, 10.3847/1538-4357/ab2888

[Kelly et al., 2018]
Kelly, P. L., Diego, J. M., Rodney, S., et al. 2018, Nature Astronomy,
2, 334, 10.1038/s41550-018-0430-3

[Kim et al., 2021]
Kim, J., Jee, M. J., Hughes, J. P., et al. 2021, , 923, 101,

[Landsman, 1993]
Landsman, W. B. 1993, in Astronomical Society of the Pacific Conference
Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed.
R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 246

[Leung et al., 2018]
Leung, E., Broadhurst, T., Lim, J., et al. 2018, , 862, 156,

[Levesque & Massey, 2020]
Levesque, E. M., & Massey, P. 2020, , 891, L37,

[Lindner et al., 2014]
Lindner, R. R., Baker, A. J., Hughes, J. P., et al. 2014, , 786,
49, 10.1088/0004-637X/786/1/49

[Madau, 1995]
Madau, P. 1995, , 441, 18, 10.1086/175332

[Meena et al., 2022]
Meena, A. K., Arad, O., & Zitrin, A. 2022, , 514, 2545,

[Menanteau et al., 2012]
Menanteau, F., Hughes, J. P., Sifón, C., et al. 2012, , 748, 7,

[Meneghetti et al., 2017]
Meneghetti, M., Natarajan, P., Coe, D., et al. 2017, , 472, 3177,

[Miralda-Escude, 1991]
Miralda-Escude, J. 1991, , 379, 94, 10.1086/170486

[Molnar & Broadhurst, 2015]
Molnar, S. M., & Broadhurst, T. 2015, , 800, 37,

[Molnar & Broadhurst, 2018]
—. 2018, , 862, 112, 10.3847/1538-4357/aad04c

[Naidu et al., 2022]
Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022, arXiv e-prints,

[Ng et al., 2015]
Ng, K. Y., Dawson, W. A., Wittman, D., et al. 2015, , 453, 1531,

[Oguri et al., 2018]
Oguri, M., Diego, J. M., Kaiser, N., Kelly, P. L., & Broadhurst, T.
2018, , 97, 023518, 10.1103/PhysRevD.97.023518

[Pascale et al., 2022]
Pascale, M., Frye, B. L., Diego, J., et al. 2022, arXiv e-prints,

[Planck Collaboration et al., 2016]
Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, ,
594, A27, 10.1051/0004-6361/201525823

[Prieto et al., 2009]
Prieto, J. L., Sellgren, K., Thompson, T. A., & Kochanek, C. S. 2009,
, 705, 1425, 10.1088/0004-637X/705/2/1425

[Robotham et al., 2018]
Robotham, A. S. G., Davies, L. J. M., Driver, S. P., et al. 2018,
, 476, 3137, 10.1093/mnras/sty440

[Robotham et al., 2017]
Robotham, A. S. G., Taranu, D. S., Tobar, R., Moffett, A., & Driver,
S. P. 2017, , 466, 1513, 10.1093/mnras/stw3039

[Schrabback et al., 2018]
Schrabback, T., Applegate, D., Dietrich, J. P., et al. 2018, ,
474, 2635, 10.1093/mnras/stx2666

[Sendra et al., 2014]
Sendra, I., Diego, J. M., Broadhurst, T., & Lazkoz, R. 2014, ,
437, 2642, 10.1093/mnras/stt2076

[Vega-Ferrero et al., 2019]
Vega-Ferrero, J., Diego, J. M., & Bernstein, G. M. 2019, ,

[Venumadhav et al., 2017]
Venumadhav, T., Dai, L., & Miralda-Escudé, J. 2017, , 850, 49,

[Waizmann et al., 2012]
Waizmann, J. C., Ettori, S., & Moscardini, L. 2012, , 420, 1754,

[Watson et al., 2014]
Watson, W. A., Iliev, I. T., Diego, J. M., et al. 2014, , 437,
3776, 10.1093/mnras/stt2173

[Welch et al., 2022]
Welch, B., Coe, D., Diego, J. M., et al. 2022a, , 603,
815, 10.1038/s41586-022-04449-y

[Welch et al., 2022]
Welch, B., Coe, D., Zackrisson, E., et al. 2022b, arXiv
e-prints, arXiv:2208.09007.

[Williamson et al., 2011]
Williamson, R., Benson, B. A., High, F. W., et al. 2011, , 738,
139, 10.1088/0004-637X/738/2/139

[Windhorst et al., 2018]
Windhorst, R. A., Timmes, F. X., Wyithe, J. S. B., et al. 2018, ,
234, 41, 10.3847/1538-4365/aaa760

[Windhorst et al., 2022]
Windhorst, R. A., Cohen, S. H., Jansen, R. A., et al. 2022, arXiv
e-prints, arXiv:2209.04119.

[Zavala et al., 2022]
Zavala, J. A., Buat, V., Casey, C. M., et al. 2022, arXiv e-prints,

[Zhang et al., 2015]
Zhang, C., Yu, Q., & Lu, Y. 2015, , 813, 129,

[Zitrin et al., 2013]
Zitrin, A., Menanteau, F., Hughes, J. P., et al. 2013, , 770, L15,

[Zitrin et al., 2015]
Zitrin, A., Fabris, A., Merten, J., et al. 2015, , 801, 44,