# Are JWST/NIRCam color gradients in the lensed $z=2.3$ dusty star-forming
galaxy El Anzuelo due to central dust attenuation or inside-out galaxy growth?

Patrick S. Kamieneski School of Earth and Space Exploration, Arizona State
University, PO Box 871404, Tempe, AZ 85287-1404, USA  Brenda L. Frye Steward
Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ, 85721-0009
Massimo Pascale Department of Astronomy, University of California, 501
Campbell Hall #3411, Berkeley, CA 94720, USA  Seth H. Cohen School of Earth
and Space Exploration, Arizona State University, PO Box 871404, Tempe, AZ
85287-1404, USA  Rogier A. Windhorst School of Earth and Space Exploration,
Arizona State University, PO Box 871404, Tempe, AZ 85287-1404, USA  Rolf A.
Jansen School of Earth and Space Exploration, Arizona State University, PO Box
871404, Tempe, AZ 85287-1404, USA  Min S. Yun Department of Astronomy,
University of Massachusetts, Amherst, MA 01003, USA  Cheng Cheng Chinese
Academy of Sciences South America Center for Astronomy, National Astronomical
Observatories, CAS, Beijing 100101, China  Jake S. Summers School of Earth and
Space Exploration, Arizona State University, PO Box 871404, Tempe, AZ
85287-1404, USA  Timothy Carleton School of Earth and Space Exploration,
Arizona State University, PO Box 871404, Tempe, AZ 85287-1404, USA  Kevin C.
Harrington European Southern Observatory, Alonso de Córdova 3107, Vitacura,
Casilla 19001, Santiago de Chile, Chile  Jose M. Diego Instituto de Física de
Cantabria (CSIC-UC). Avda. Los Castros s/n. 39005 Santander, Spain  Haojing
Yan Department of Physics and Astronomy, University of Missouri, Columbia, MO
65211, USA  Anton M. Koekemoer Space Telescope Science Institute, 3700 San
Martin Dr., Baltimore, MD 21218, USA  Christopher N. A. Willmer Steward
Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ, 85721-0009
Andreea Petric Space Telescope Science Institute, 3700 San Martin Dr.,
Baltimore, MD 21218, USA  Lukas J. Furtak Physics Department, Ben-Gurion
University of the Negev, P.O. Box 653, Be’er-Sheva 84105, Israel  Nicholas Foo
University of Arizona, Department of Astronomy/Steward Observatory, 933 N
Cherry Ave, Tucson, AZ 85721, USA  Christopher J. Conselice Jodrell Bank
Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13
9PL, UK  Dan Coe AURA for the European Space Agency (ESA)  Space Telescope
Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  Simon P.
Driver International Centre for Radio Astronomy Research (ICRAR) and the
International Space Centre (ISC), The University of Western Australia, M468,
35 Stirling Highway, Crawley, WA 6009, Australia  Norman A. Grogin Space
Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA
Madeline A. Marshall National Research Council of Canada, Herzberg Astronomy &
Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7,
Canada  ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions
(ASTRO 3D), Australia  Nor Pirzkal Space Telescope Science Institute, 3700 San
Martin Dr., Baltimore, MD 21218, USA  Aaron S. G. Robotham International
Centre for Radio Astronomy Research (ICRAR) and the International Space Centre
(ISC), The University of Western Australia, M468, 35 Stirling Highway,
Crawley, WA 6009, Australia  Russell E. Ryan, Jr. Space Telescope Science
Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  Scott Tompkins
School of Earth and Space Exploration, Arizona State University, PO Box
871404, Tempe, AZ 85287-1404, USA

###### Abstract

Gradients in the mass-to-light ratio of distant galaxies impede our ability to
characterize their size and compactness. The long-wavelength filters of JWST’s
NIRCam offer a significant step forward. For galaxies at Cosmic Noon ($z\sim
2$), this regime corresponds to the rest-frame near-infrared, which is less
biased towards young stars and captures emission from the bulk of a galaxy’s
stellar population. We present an initial analysis of an extraordinary lensed
dusty star-forming galaxy (DSFG) at $z=2.3$ behind the El Gordo cluster
($z=0.87$), named El Anzuelo (“The Fishhook") after its partial Einstein-ring
morphology. The FUV-NIR SED suggests an intrinsic star formation rate of
$81^{+7}_{-2}~{}M_{\odot}~{}{\rm yr}^{-1}$ and dust attenuation $A_{V}\approx
1.6$, in line with other DSFGs on the star-forming main sequence. We develop a
parametric lens model to reconstruct the source-plane structure of dust imaged
by the Atacama Large Millimeter/submillimeter Array, far-UV to optical light
from Hubble, and near-IR imaging with 8 filters of JWST/NIRCam, as part of the
Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS)
program. The source-plane half-light radius is remarkably consistent from
$\sim 1-4.5~{}\mu$m, despite a clear color gradient where the inferred galaxy
center is redder than the outskirts. We interpret this to be the result of
both a radially-decreasing gradient in attenuation and substantial spatial
offsets between UV- and IR-emitting components. A spatial decomposition of the
SED reveals modestly suppressed star formation in the inner kiloparsec, which
suggests that we are witnessing the early stages of inside-out quenching.

Strong gravitational lensing(1643) — Starburst galaxies(1570) — James Webb
Space Telescope(2291)

††facilities: JWST(NIRCam), HST(ACS,WFC3), ALMA††software: APLpy (Robitaille &
Bressert 2012; Robitaille 2019) astropy (Astropy Collaboration et al. 2013),
blobcat (Hales et al. 2012), casa (McMullin et al. 2007), eazy (Brammer et al.
2008; Brammer 2021), galfit(Peng et al. 2002), lenstool (Kneib et al. 1993,
1996; Jullo et al. 2007; Jullo & Kneib 2009), Ned Wright’s Cosmology
Calculator (Wright 2006), photutils (Bradley et al. 2022) SExtractor (Bertin &
Arnouts 1996)

##  1 Introduction

The recent advent of JWST has already offered substantial new insight into the
Universe at the epoch of Cosmic Noon, $1\lesssim z\lesssim 3$, during which
the cosmic star-formation history was at its zenith (Madau & Dickinson 2014).
Dusty star-forming galaxies (DSFGs)—also referred to historically as
submillimeter galaxies, or SMGs (see reviews by Blain et al. 2002; Casey et
al. 2014a)—have been studied extensively in the past few decades as tracers of
the rapid stellar mass assembly of the Universe that peaked 10 billion years
ago. Their substantial shrouds of dust are efficient at reprocessing the
intense ultraviolet (UV) and visible light from sites of rapid, active star
formation into infrared (IR) radiation, resulting in large observed-frame
(sub)-millimeter fluxes. Yet, their dust content also has the effect of
strongly attenuating and reddening their rest-frame UV/optical emission, such
that it has been a challenge to study them with even the most sensitive
optical telescopes, including the Hubble Space Telescope (HST). When detected,
the peak of UV emission in DSFGs is often significantly offset from that of
the far-IR (e.g., Goldader et al. 2002; Chapman et al. 2004; Hodge et al.
2015; Chen et al. 2017; Calistro Rivera et al. 2018; Cochrane et al. 2021).
Yet, a subset of these objects have gone undetected in deep $H$-band imaging
(to limiting $5\sigma$ depths of $H>27$ mag), leading to their classification
as “HST-dark" or “OIR-dark" galaxies (e.g., Huang et al. 2011; Simpson et al.
2014; Franco et al. 2018; Wang et al. 2019b, and references therein). JWST has
proven capable of shedding new light on these elusive galaxy populations,
which will offer substantial refinements to theories of galaxy evolution
(e.g., Barrufet et al. 2023; Ferreira et al. 2022).

Galaxy size is a key observable property that has been used to advance our
understanding of the assembly history of galaxies throughout cosmic time (Mo
et al. 1998). In particular, the size-mass relation ($M_{\star}-R_{e}$; e.g.,
Shen et al. 2003; Trujillo et al. 2004; van der Wel et al. 2014a; Mowla et al.
2019) indicates that more massive galaxies are typically larger in size, with
markedly steeper slopes for quiescent (late-type) galaxies than for star-
forming (early-type) galaxies. Ryan et al. (2012) found that the evolution of
the size of passive galaxies with redshift is itself stellar mass-dependent
(such that the most massive galaxies experienced the greatest evolution from
$z$$\sim$2 to the present). However, it is not trivial to assign the effective
radius of a galaxy, as the distribution of light is known to vary strongly
with wavelength, such that $K$-band sizes can be as small as half those at
$r$-band (e.g., Evans 1994; Cunow 2001; La Barbera et al. 2002, 2010; Kelvin
et al. 2012; Häußler et al. 2013; Vulcani et al. 2014; Kennedy et al. 2015,
2016). This often results in radial color gradients, or variations in the
mass-to-light ratios, which are usually explained by some combination
of111Observation-specific effects, such as unmatched PSFs at different
wavelengths, can also be responsible.:

  * •

centrally-concentrated dust, leading to greater attenuation in the inner
regions of a galaxy (e.g., Jansen et al. 1994; Peletier & Balcells 1996;
Möllenhoff et al. 2006; Graham & Worley 2008; Pastrav et al. 2013; Pantoni et
al. 2021)

  * •

inside-out galaxy growth and quenching (Carrasco et al. 2010; van Dokkum et
al. 2010; Hopkins et al. 2010), by which stellar populations in the center are
predominantly older, especially when a spheroidal/bulge component is present
(e.g., Driver et al. 2007a, b; La Barbera et al. 2010; Marian et al. 2018)

  * •

gradients in metallicity, often driven by inflow of primordial gas to
galaxies’ gravitational potential (e.g., Cresci et al. 2010; Jones et al.
2010; Yuan et al. 2011; Pilkington et al. 2012)

Typically, galaxies with color gradients are redder in the center, but this is
especially true for more luminous objects (Jansen et al. 2000) and for late-
type disk-dominated galaxies (de Jong 1996; Gadotti & dos Anjos 2001).

In effect, the implication of these color gradients is that optical filters
are not always effective proxies for the underlying stellar mass distribution,
and the mass-to-light ratio can vary with radius in even a single galaxy
(e.g., Suess et al. 2019). More worrisome is the finding by Suess et al. that
the amplitude of color gradients can correlate with properties like stellar
mass or surface mass density, such that our interpretation of scaling
relations (like $M_{\star}-R_{e}$) that depend on light-weighted sizes might
be impacted significantly. Rest-frame near-infrared light, on the other hand,
arises primarily from stellar populations of cooler, low-mass stars that make
up most of a galaxy’s stellar mass. Additionally, near-IR is comparatively
much less affected by the dust attenuation that plagues the rest-frame UV and
bluer optical regime (by approximately an order of magnitude; e.g., Li &
Draine 2001). Until recently, however, observations of the high-$z$ Universe
at this wavelength have been limited to unresolved, galaxy-integrated studies.

Now, JWST has enabled resolved measurements of the light distribution at
wavelengths beyond $\sim$2 $\mu$m in $z>1$ galaxies for the very first time.
In an examination of $\sim 1000$ galaxies from the Cosmic Evolution Early
Release Science program (CEERS; see Finkelstein et al. 2023), Suess et al.
(2022) found that galaxies were systematically smaller at 4.4 $\mu$m
($\lambda_{\rm rest}\sim 1.6\mu$m) than at 1.5 $\mu$m ($\lambda_{\rm rest}\sim
550$ nm) by about 9% on average (but reaching $\sim 30\%$ for high-mass
galaxies, $M_{\star}\sim 10^{11}~{}M_{\odot}$). It is now feasible to directly
probe the distribution of stellar mass in the extragalactic Universe with the
long-wavelength filters of JWST’s NIRCam instrument, without needing to model
the galaxy spectral energy distribution (SED) or infer any gradients in mass-
to-light ratio (e.g., Suess et al. 2019).

To fully take advantage of this new capability, the best approach is to use
JWST to observe objects affected by strong gravitational lensing. The
amplification of flux helps to facilitate the mapping of the stellar continuum
out to even larger radii, where the light profile would otherwise fall below
the detection limit. This amplification also allows one to probe less-luminous
galaxies that are more representative of the luminosity function at $z>1$. Of
paramount importance to this work is the areal magnification provided by
lensing, which helps push angular resolution beyond even the JWST diffraction
limit, albeit with a complex pattern of distortion that must be accounted for.
A prime opportunity for using JWST to map the assembly of stellar mass at
high-redshift is with the Prime Extragalactic Areas for Reionization and
Lensing Science (PEARLS) program, which targets seven massive lensing clusters
(Windhorst et al. 2023). One such cluster, ACT-CL J0102-4915, or “El Gordo"
(Menanteau et al. 2012), has been identified as a foreground for a rich galaxy
group at $z$=4.3 (Caputi et al. 2021; Frye et al. 2023), including DSFGs (see
also Cheng et al. 2023). In this work, we focus on an extraordinary, recently-
discovered lensed DSFG J010249-491507 (shown in Fig. 1), which we find to be
near the peak of Cosmic Noon, at $z$=2.291. The two primary lensed images were
denoted as EG-ALMA2a and EG-ALMA2b by Cheng et al. (2023). Diego et al. (2023)
first gave this object the moniker El Anzuelo (Spanish for “Fishhook," given
the shape of the lensing morphology), which we will use for the remainder of
this work. This lensed system has been revealed in much greater detail than
was possible in previous optical imaging with Hubble. Thus, we consider it an
excellent prospect for an in-depth analysis of mapping the light vs. mass
profile in an extreme star-forming object—all with superb, sub-diffraction-
limit angular resolution thanks to lensing.

Figure 1:  RGB image of El Anzuelo, which is aligned with north, $13.5\arcsec$
on each side, and centered on $(\alpha,\delta)=$
(01h02m49$\fs$452,$-49^{d}$15m07$\fs$05). For this JWST-only image, the red,
green, and blue channels are constructed by weighting filters as follows:
$R=[(0.3\times{\rm F356W})+(0.6\times{\rm F410M})+(1.0\times{\rm F444W})]$;
$G=[(0.8\times{\rm F200W})+(1.0\times{\rm F277W})+(0.5\times{\rm F356W})]$;
$B=[(1.0\times{\rm F090W})+(1.0\times{\rm F115W})+(1.0\times{\rm
F150W})+(1.0\times{\rm F200W})+(1.0\times{\rm F277W})+(0.8\times{\rm
F356W})+(0.5\times{\rm F410M})+(0.5\times{\rm F444W})]$.

This paper is organized as follows: in Section 2, we describe the JWST, Hubble
Space Telescope (HST), and Atacama Large Millimeter/submillimeter (ALMA) data
used in this work. In Section 3, we detail our analysis, including lens
modeling, subtraction of foreground light, photometric redshift estimation,
characterization of source-plane size and light profile, and determination of
the (wavelength-dependent and location-dependent) magnification. In Section 4,
we discuss this differential magnification, and weigh the evidence of various
effects that might be responsible for the observed radial color gradients. We
discuss a decomposition of the spectral energy distribution for the inner-
versus outer-disk, and what this reveals about possible inside-out galaxy
quenching. Our conclusions are summarized in Section 5.

Throughout this work, we use a concordance $\Lambda$CDM cosmological model,
with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}=70~{}{\rm
km~{}s}^{-1}~{}{\rm Mpc}^{-1}$. At the redshift of El Anzuelo ($z=2.291$), the
angular-to-physical size conversion with this cosmology is $1\arcsec=8.209$
kpc, and at the redshift of the El Gordo cluster ($z=0.87$), the conversion is
$1\arcsec=7.713$ kpc (Wright 2006). The magnitudes in this paper are given
according to the AB system. All images are aligned with north up and east to
the left, unless otherwise noted.

##  2 Data

###  2.1 Object selection

El Anzuelo was identified as a candidate for an in-depth analysis because of
the relative simplicity of the foreground lensing environment, in contrast
with the large-scale cluster mass profile. As there are fewer cluster galaxies
contributing significantly to the lensing deflection, we can derive a robust
model with few parameters. Moreover, we can fine-tune the model to minimize
the separations between observed and reconstructed multiple image locations
for this lensed object alone, rather than a large collection of image systems
(e.g. Diego et al. 2023; Frye et al. 2023). As we discuss in §3.5, we adopt a
source redshift of $z=2.291$ based on SED-fitting and a CO line detection.

###  2.2 JWST NIRCam

The El Gordo cluster was selected to be included in JWST’s PEARLS program
(Windhorst et al. 2023) because of its extremely large halo mass, on the order
of $M_{h}$ $\sim$ 2$\times 10^{15}~{}\rm M_{\odot}$ (Menanteau et al. 2012). A
more detailed description of this data, collected on 2022 July 29, and its
subsequent reduction are provided in Windhorst et al. (2023), Diego et al.
(2023), and Frye et al. (2023), but we summarize relevant details here. El
Gordo was observed in the wide-band filters F090W, F115W, F150W, F200W, F277W,
F356W, F444W, and the medium-band filter F410M. Exposure times were 2491
seconds for F090W, F115W, F410M, and F444W; 2104 seconds for F200W and F277W;
and 1890 seconds for F150W and F356W. The data were processed with pipeline
version 1.6.2 and the context file jwst_0942.pmap. The $5\sigma$ limiting
magnitudes for all filters are in the range $m_{\rm AB}\simeq 28-29$ mag.

###  2.3 HST

Existing ancillary imaging with the Hubble Space Telescope (HST) provides
helpful constraints for the rest-frame UV regime in El Anzuelo. The filters
used in this analysis include Advanced Camera for Surveys (ACS) Wide Field
Channel (WFC) F606W, F625W, F775W, F814W, F850LP; and Wide Field Camera 3
(WFC3) IR channel F105W, F125W, F140W, and F160W. These data were collected as
part of the Reionization Lensing Cluster Survey222https://relics.stsci.edu/
(RELICS), through HST programs GO 12755 (PI: J. Hughes), GO 12477 (PI: F.
High), and GO 14096 (PI: D. Coe). The observations were taken in 2012
September, 2012 October, and 2016 July/August, respectively. Data were
retrieved from the Mikulski Archive for Space
Telescopes333https://archive.stsci.edu/ (MAST). We refer the reader to Coe et
al. (2019) for a more complete description of these observations.

###  2.4 ALMA Bands 3, 6, and 7

El Anzuelo in the El Gordo field was observed with the Atacama Large
Millimeter/submillimeter Array (ALMA) in programs 2015.1.01187.S (PI: K. Basu)
and 2017.1.01621.S (PI: K. Basu) in Band 3 (3 mm or 100 GHz), 2013.1.01358.S
(PI: A. Baker) and 2018.1.00035.L (PI: K. Kohno) in Band 6 (1 mm or 270 GHz),
and 2013.1.01051.S (PI: P. Aguirre) in Band 7 (870 $\mu$m or 340 GHz). All
data are public and were retrieved from the ALMA
archive444https://almascience.nrao.edu/aq/.

The reduced 3 mm data are presented by Basu et al. (2016). Coarse-resolution
imaging with a synthesized beam size of $3.6\arcsec\times 2.7\arcsec$ (PA
$=88^{\circ}$) was taken in 2015 December with the most compact antenna
configuration, with a noise rms of 6 $\mu$Jy beam-1. While designed as a
continuum experiment, these data were used for spectroscopic redshift
confirmation of the background DSFG (see §3.5). The angular resolution is not
sufficient to perform a reasonable source-plane reconstruction.

The two 1 mm measurement sets were combined in the complex visibility (i.e.,
Fourier) space and imaged together by Cheng et al. (2023). The resulting
mosaic has a beam size of $1.28\arcsec\times 0.93\arcsec$ (PA $=88^{\circ}$)
and a rms sensitivity of $\sim$0.065 mJy beam-1. Likewise, the 870 $\mu$m map
has a beam size of $0.33\arcsec\times 0.40\arcsec$ (PA $=-35^{\circ}$), with
noise rms 0.25 mJy beam-1. For both the 870 $\mu$m and 1mm data, images were
created in casa v.6.2.1 (McMullin et al. 2007) with Briggs weighting
(robust$=2$) and cleaned with tclean down to the $3\sigma$ level.

##  3 Analysis

For this case study of El Anzuelo, we develop a gravitational lens model,
measure filter-dependent magnifications, reconstruct the observed data in the
source plane, and characterize the spatial extent as a function of wavelength.

Table 1: Lens modeling results for El Anzuelo. All potentials are held fixed
at the cluster redshift of $z=0.87$.

Potential ID & type | RA | Dec | $e$ | PA |  $\sigma$ (or $\gamma$)
---|---|---|---|---|---
| [h:m:s] | [d:m:s] |  | [∘] | [km s-1]
L1 (SIE) | 1:02:49.347† | -49:15:06.60† | 0.16 $\pm$ 0.13 | 27 $\pm$ 39 | 302.6 $\pm$ 4.6
L2 (SIE) | 1:02:49.269† | -49:15:05.33† | 0.44 $\pm$ 0.19‡ | 71 $\pm$ 29 | 219.7 $\pm$ 7.4
L3 (SIE) | 1:02:49.665† | -49:15:08.22† | 0.75† | -50† | 27 $\pm$ 68‡
L4 (ex. shear) | — | — | — | 99 $\pm$ 6 | 0.36 $\pm$ 0.08

^†^†footnotetext: Parameter value held fixed. For lens L3, the ellipticity and
orientation are set at those derived with galfit, as the lens model appears to
be quite insensitive to the mass at this projected distance.

^‡^‡footnotetext: Parameters for which there is a sizable disagreement between
the highest-likelihood and median solution include the ellipticity of L2
(which has a median $e=0.23$) and the velocity dispersion of L3 (which has a
median $\sigma=98$ km s-1).

Note. — Here we report the parameter values for the highest-likelihood MCMC
iteration. These are consistent within 1$\sigma$ uncertainties with the median
of the posterior distribution except for the two cases mentioned above. Where
possible in this work, we account for the uncertainty inherent to our lens
model by randomly selecting from the full posterior, as discussed in §3.1.3.

###  3.1 Lens modeling

Figure 2: Observed and modeled image plane light distribution for three short-
wavelength NIRCam filters in an RGB image (red: F200W, green: F150W, blue:
F115W). The caustic curve is shown in cyan, and the corresponding image-plane
critical curve is shown in violet. Left: Observed image, after removing
foreground lens galaxy light with galfit. Any residuals nearest the foreground
galaxies (i.e., for F115W) are not of great concern, as they are masked out
during the ray-tracing. Identified image families used as constraints in the
lens model are labeled, including $1abcde$ (identified primarily with the
long-wavelength NIRCam filters) $2ab$, and $3ab$. Center: Model-reconstructed
image-plane distribution of the same filters, made by ray-tracing the observed
data to the source-plane with the best-fit model, then ray-tracing back to the
image plane (and finally, convolving with a smoothing kernel to remove ray-
tracing discontinuity artifacts). Right: Model-reconstructed source-plane
distribution, zoomed-in to show more spatial detail. Gray contours show ALMA 1
mm continuum, ray-traced to the source-plane.

Figure 3:  As with Fig. 2, but for 3 long-wavelength NIRCam filters: F444W
(red), F356W (green), F277W (blue). In the observed image (left), we also show
the image-plane ALMA 1 mm contours in gray. Green crosses identify the
locations of apertures used for photometry, discussed in §3.3.

The lensing distortion of El Anzuelo is dominated primarily by a small group
of cluster galaxies, and not by the cluster halo itself, as evidenced by its
small Einstein radius of $\theta_{E}=1\farcs 70\pm 0\farcs 06$. For this
reason, we construct a separate lens model (in parallel to those developed
recently by Diego et al. 2023 and Frye et al. 2023) that is focused only on El
Anzuelo, so as to tailor the optimization to a small number of constraints and
parameters.

Parametric gravitational lens modeling was performed with
lenstool555https://projets.lam.fr/projects/lenstool/wiki (Kneib et al. 1993,
1996; Jullo et al. 2007; Jullo & Kneib 2009). lenstool solves for the simple
lens equation (see e.g., Schneider et al. 1992; Narayan & Bartelmann 1996),
which relates the source-plane distribution of light with the observed image-
plane distribution via the deflection angle, which is derived from an integral
of the surface mass density in the lens plane. With a set of observed
locations for a multiply-imaged family, one can constrain the foreground mass
distribution (where each separate set of $n$ images provides $2n-2$
constraints). This is done completely independently of the light profile of
the foreground mass, although this is used for setting priors on each
parameter, and for establishing the number of foreground lenses to include
(and sometimes the nature of each profile).

lenstool uses a Bayesian Markov chain Monte Carlo algorithm to sample the
multi-dimensional posterior distribution for the set of free parameters. The
likelihood function connects the image-plane positions of multiple images used
as constraints, and the resulting positions of these locations obtained by
ray-tracing to the source-plane and back to the image-plane with each model
iteration (as shown in Figs. 2 and 3). We use 3000 iterations to extensively
sample the posterior space. The best-fit parameters are provided in Table 1,
which result in an RMS deviation between observed and modeled image positions
of $0.06\arcsec$. This value closely matches the FWHM of the PSF for F444W,
indicating a high goodness-of-fit. While one could choose the median or mode
of each parameter’s posterior distribution as the solution, we opt for the
highest-likelihood (minimum $\chi^{2}$) set of parameters, after verifying
that they are consistent with median and mode values within uncertainties.

The results of our lens models are depicted in Figs. 2 and 3, where the panels
show RGB images of the observed light distribution (after subtracting
foreground light with galfit, as described in Section 3.1.2) for the short-
wavelength and long-wavelength filters separately. The middle panels show a
model-reconstructed image-plane, created by ray-tracing the original observed
image into the source plane, before then ray-tracing this result back to the
image plane (and smoothing by a kernel equal to half the PSF FWHM to reduce
artifacts). As lensing is not a one-to-one transformation—e.g., ray-tracing
one multiply-imaged feature to the source-plane and back to the image-plane
should recover all members of the image family—this translation to the source-
plane and back is likely to introduce unwanted features in the case of poor
goodness-of-fit. In the right panels, we show the source-plane reconstruction
itself. In Fig. 5, we show the same results for the ALMA 870 $\mu$m and 1 mm
continuum images.

Table 2: Locations of secure multiple image family members that were used in the lens model, which are also labeled in Figs. 2 and 4.  Image ID | RA | Dec
---|---|---
| [h:m:s] | [d:m:s]
1a | 1:02:49.188 | -49:15:08.77
1b | 1:02:49.298 | -49:15:04.58
1c | 1:02:49.431 | -49:15:06.12
1d | 1:02:49.332 | -49:15:04.92
1e | 1:02:49.320 | -49:15:05.49
2a | 1:02:49.122 | -49:15:06.04
2b | 1:02:49.130 | -49:15:05.50
3a | 1:02:49.411 | -49:15:05.80
3b | 1:02:49.375 | -49:15:05.45

####  3.1.1 Lens model parameters

Here, we have used a common approach where each foreground elliptical is
parameterized as a singular isothermal ellipsoid (SIE) or spheroid (SIS),
which is efficient for describing galaxy- or group-scale lensing (e.g.,
Kormann et al. 1994; Treu & Koopmans 2004). We also include an external shear
term to account for any weak lensing effect from the El Gordo cluster.
However, recent work by Etherington et al. (2023) suggests that the vast
majority of models with external shear are not actually affected by weak
lensing, and this added factor only compensates for insufficient model
complexity with the chosen parameterization (albeit, this result is primarily
based on isolated field galaxies). We note that the best-fit orientation of
the external shear term is PA$\approx$99∘ (counterclockwise from east), which
is not closely aligned with the major or minor axis of the most massive
foreground galaxy (PA$\approx$27∘), so we preserve this as a free parameter,
with the acknowledgment that this is not necessarily representing actual
external shear.

Figure 4:  Image families used as constraints for the lens model. The center
panel shows an RGB image of all 8 NIRCam filters, with the lens model-derived
caustic and critical curves shown in cyan and purple. Insets show the
identified image families $1abcde$ (most prominent in F444W) and $2ab$ and
$3ab$ (prominent in F200W). Pink contours in the $1abcde$ cutouts show ALMA
1mm emission at levels of $[5,10,15]\sigma$.

We identify 3 image families (containing multiplicities of 2, 2, and 5 images,
respectively) labeled in Fig. 4. With $2(n-1)$ constraints provided by each
family of $n$ images, this results in 8 constraints from image set $1abcde$
and 2 constraints each from $2ab$ and $3ab$. The locations of each image are
provided in Table 2. The observed dust continuum structure suggests an overall
double-image morphology with a southwest and northeast component, but precise
image locations are chosen carefully from NIRCam imaging. In particular, the
structure in F444W has a prominent peak in the southwest image ($1a$), which
appears as multiple peaks with a similar color in the northeast side of the
arc, which we identify as multiple images of the same region. An initial,
oversimplified mass model of the two primary foreground galaxies—with
centroids, ellipticities, and orientations held fixed to the best-fit galfit
results (and only velocity dispersions as free parameters)—predicted a similar
five-image morphology, helping to confirm our identification.

For image families $2ab$ and $3ab$, F200W in particular reveals an observable
symmetry that arises from images in close proximity on opposite sides of the
critical curve. The pairs of prominent clumps of similar size and luminosity
are circled as $2ab$ and $3ab$ in Fig. 4. These are incorporated into the
model iteratively, by first using the initial mass model to determine that
these pairs of clumps could logically lie in the vicinity of the critical
curve. Their eventual inclusion in the model helps to refine and further
constrain the mass distribution beyond just the sightlines probed by images
$1abcde$, but adding them does not significantly alter the best-fit
parameters. It may be feasible to identify additional multiply-imaged
features, but this set of secure constraints is conservative and more than
sufficient given the simplicity of our lens model parameterization, in
contrast with the models of the full El Gordo cluster. In the future, resolved
spectral line observations would add further confirmation of these image
identifications (which should share similar line profiles). However, the
fidelity of the image-plane reconstructions in Figs. 2 and 3 (even in regions
not directly constrained as an image family, such as the highly-magnified
structure between images $1a$ and $2a$) lends further credence that our image
identifications are appropriate.

The two primary foreground elliptical galaxies are modeled as singular
isothermal ellipsoids, as is the third foreground at greater projected
distance (but with ellipticity and position angle held fixed to galfit fitting
results). The $x$ and $y$ (RA and Dec.) centers of each lens are held fixed to
the median light centroid of the JWST filters, as otherwise large degeneracies
between the different lensing galaxies can be introduced. These galaxies are
labeled in Fig. 1, all of which are spectroscopically-confirmed cluster
members at $z=0.87$ (Caminha et al. 2022). A perturbative shear factor is
added to account for any missing complexity from our model, as we don’t
account for the halo potential of the El Gordo cluster. In total, this
combines to 9 free parameters (Table 1), leaving 3 degrees of freedom in the
model. While additional free parameters could be introduced, they are not
overwhelmingly physically-motivated, and so we opt for the simplest
parameterization feasible.

Figure 5:  As with Fig. 2, but for the ALMA 1mm continuum image (with the
lower signal-to-noise 870 $\mu$m image shown as red contours). A different
field-of-view is chosen to show more of the noise variations. The synthesized
beam is shown in cyan (1 mm) and filled red (870 $\mu$m) in the lower-left
corner. Noise in the 870 $\mu$m image increases towards the southwest, which
is closer to the edge of the primary beam. In the middle panel, contours of
the observed image are plotted over the model reconstruction for ease of
comparison. Here, the red contours now show the modeled 870 $\mu$m image-plane
emission. In the last panel showing the source plane (870 $\mu$m as red
contours), the beams are shown as dashed ellipses. Lensing allows one to over-
resolve the structure (and the actual source-plane PSF varies due to the
magnification gradient). It is unclear if the dust continuum actually extends
further into the center of caustic curve relative to the stellar continuum
(cf. Figs. 2 and 3), which might suggest an outflow from the plane of the
disk. However, given the much larger PSF for the ALMA 1 mm image relative to
JWST ($\theta$ $\sim$1.1$\arcsec$ vs. $\theta$$\sim$0.03-0.15$\arcsec$), it’s
not possible to answer this at present.

The more disk-like lensing galaxy a little further to the southeast, L3, is
somewhat poorly constrained by our model, and the best-fit velocity dispersion
is an un-physical value of $\sigma\approx 27$ km s-1. The uncertainty on this
value is large and highly asymmetric, but we still include it in our model
with the intention of examining more closely the amplitude of the external
shear, which we find to be $\gamma=0.36\pm 0.08$. This value is large, but not
entirely unexpected for a complex mass distribution like a group of galaxies
embedded within a much larger, extremely massive cluster. Moreover, the unique
morphology of the lensed background galaxy suggests a more complicated set of
lensing caustics than are likely to be provided by two foreground ellipticals.

To test this further, we re-run the optimization and hold the shear amplitude
fixed at $\gamma=0.2$, to be approximately consistent with the contribution
from the cluster potential as estimated by models from Diego et al. (2023) and
Frye et al. (2023). This results in a consistent set of parameters—i.e.,
within 1$\sigma$ from the highest-likelihood values—except for the ellipticity
of L1 and the velocity dispersion of L3. These shift to $e_{1}\approx 0.32$
and $\sigma_{3}\approx 169$ km s-1, but both of these are well within the
respective 2$\sigma$ bounds from the best-fit solution, and are actually
within $1\sigma$ of the median. This offers further evidence that the
amplitude of external shear included in the best-fit model may be driven by
insufficient flexibility in the parameterization of the foreground galaxies.

On the other hand, the orientation of the shear component is consistent with
what would be expected for that contributed by the cluster. The center of mass
adopted by Diego et al. (2023) between the two peaks of the bimodal mass
distribution is oriented at PA=$27^{\circ}$ (CCW from east) from El Anzuelo,
while the nearest, northern mass component of the cluster is centered on a
luminous galaxy that is oriented at PA=$6^{\circ}$. This central galaxy also
coincides with the surface mass density peak of the northwest component from
the independent lens model of Frye et al. (2023). External shear at PA=$99\pm
6^{\circ}$ is thus very close to perpendicular to the direction of the
cluster’s mass center, and consistent within uncertainty with the direction of
the northern component’s mass center. With this context, it seems likely that
the external shear we find may be due in part to weak lensing from the
cluster, while the large amplitude may be explained by some level of
additional shear that is required to best describe the mass distribution of
the galaxy group lensing the background DSFG.

####  3.1.2 galfit foreground subtraction

Before reconstructing the observed data in the source-plane, we subtract a
Sérsic model of the foreground light of the massive lensing ellipticals using
galfit (Peng et al. 2002). We note that a detailed assessment of the light
profiles is not directly related to the aims of this paper, and the central
regions of each foreground lens (where there might be large residuals from
galfit subtraction) are masked out before ray-tracing to the source plane, and
thus don’t significantly impact our analysis.

Two Sérsic profiles were used to model the light from the two foreground
elliptical cluster galaxies (L1 and L2) closest to the light from El Anzuelo
(see Fig. 4) The position, flux, effective radius, Sérsic index, axis ratio,
and position angle were all kept as free parameters. Each filter (for both
JWST and HST) was optimized independently of each other, but checked for
consistency. For the F150W and F200W filters, it was necessary for the
centroids of each Sérsic profile to be held fixed to match the average
centroid of the six other JWST filters in order for the models to converge.
The weight map was constructed by galfit, applied to a $5\arcsec\times
7\arcsec$ cutout as shown in Figs. 2 and 3.

There is non-negligible residual at the far eastern end of El Anzuelo (the
“barb" of the fishhook) in some of the shorter-wavelength filters (visible
predominantly in the left panel of Fig. 2). This is likely due to significant
blending of emission from the background galaxy with the core of the
foreground elliptical. In the middle panel of Fig. 2, we create an image-plane
reconstruction by masking out this portion of the arc and ray-tracing to the
source-plane and back to the image-plane. As a result, the lens model predicts
that there is significant highly-magnified structure from a small portion of
the background galaxy in this region, which means that simply adding a second
Sérsic component to the foreground model is unlikely to improve the residual
meaningfully. For our analysis in this work, we are able to mask out this
region and exclude it from our measurements, as the southwestern image ($1a$)
that we focus on contains the entire background object and is significantly
less affected by foreground contamination. However, future work investigating
this highly-magnified structure would warrant a more careful removal of the
foreground light to reveal the lensed structure.

####  3.1.3 Magnification

The total magnification at each optical/near-IR wavelength is determined by
computing the flux-weighted average of the magnification map. To measure
uncertainties, we randomly draw 100 different models from each parameter’s
full posterior distribution from MCMC optimization. For each iteration, we
reconstruct the source-plane structure, and a map of source-plane
magnification at the same pixel scale. Since no assumption is made about the
distribution of each parameter’s posterior distribution, the resultant
distribution of magnifications should accurately capture random uncertainties
(although there may still be systematic uncertainties). In most cases, the
magnifications are normally-distributed, so we cite $1\sigma$ errors on the
median magnification in Table 3.

For the 1mm ALMA images, we take a slightly different approach, as the flux
distribution between pixels is highly correlated within the beam size, which
complicates a calculation of the flux-weighted pixel-based average of the
magnification map. Instead, we estimate magnification to be $\mu_{1{\rm
mm}}=9.5\pm 0.9$ by using the ratio of image-plane to source-plane areas above
an equivalent noise threshold. For the 870 $\mu$m image, on the other hand,
the DSFG is only marginally detected (Cheng et al. 2023), so we estimate the
magnification by using the range of values in a small aperture centered on the
brightest pixel (see Fig. 5), resulting in $\mu_{870\mu{\rm m}}\approx 6-9$.

###  3.2 Source-plane size and light profile

To reconstruct the source-plane light distribution, we ray-trace all lensed
images (e.g., $1abcde$) to the source-plane simultaneously. As the PSF varies
over the extent of the source due to lensing distortion, the choice of a
single kernel for galfit is nuanced. Since the PSF is broadest and varies the
least for image $1a$, it is the best choice for describing the galaxy’s
overall spatial extent. For our purpose, we place a simulated PSF from WebbPSF
v.1.1.0 (Perrin et al. 2012, 2014) at 10 random locations within a
$2\arcsec\times 2\arcsec$ box centered at the peak of the F444W emission in
the counterimage ($1a$), then reconstruct in the source-plane with our lensing
model. This creates a set of representative PSFs in the vicinity of the region
of greatest interest. Next, we convolve the model-reconstructed source-plane
structure with this source-plane kernel in order to mitigate small-scale
substructure within the galaxy (primarily arising from the higher-
magnification images, $1b-1d$) and to ensure a more uniform spatial
resolution. We use a single-component Sérsic profile and a sky background
component in our models for each filter. We note that the masking of
foreground objects when reconstructing in the source plane leads to only a
small area of blank sky being included. As the Sérsic indices we derive are
generally very small, this sky component likely has little effect on the
estimation of the effective radius. This is more of a concern for objects with
large best-fit $n$ indices, which have significant signal out to large radii.

Of primary interest are the effective radius along the semi-major axis $R_{e}$
(encapsulating half the total flux) and the Sérsic index (describing the
radial concentration of light). The results for each filter are provided in
Table 3. We find minimal scatter in the axis ratio $b/a=0.49\pm 0.05$, which
corresponds to an estimated inclination666We use an intrinsic axial ratio of
$q=0.2$ (e.g., Unterborn & Ryden 2008), and estimate inclination from
$\cos^{2}{i}=[(b/a)^{2}-q^{2}]/(1-q^{2})$. of $i=63^{\circ}\pm 7^{\circ}$. The
uncertainties reported from galfit $\chi^{2}$ optimization are likely severe
underestimates of the actual, empirical dispersion in $R_{e}$ and $n$. For
this reason, we instead adopt the relative uncertainties calculated for the
lensing magnification in each filter (ranging from $\Delta\mu/\mu\sim
10-30\%$; §3.1.3), which are likely to dominate the uncertainty
budget777Lensing magnifies each dimension by approximately $\mu^{0.5}$, so the
relative uncertainty in size from lensing would actually propagate to be half
the relative uncertainty in magnification. However, given the substantial
deviations from a smooth light profile, we maintain the full $10-30\%$
fractional uncertainties to be conservative.. Importantly, the scatter in
measured properties from the set of representative source-plane PSFs is also
smaller than this adopted uncertainty.

In some cases, a single-component Sérsic profile does not adequately capture
the various spatial frequencies of the galaxy’s light. In these scenarios, a
second component is added (and often held fixed at the center of the galaxy’s
nucleus when necessary), as in the case of a bulge $+$ disk model (e.g. Gao &
Ho 2017). Our criterion for requiring this extra component is a $>5\sigma$
peak in the residual near the galaxy’s center. An equivalent Sérsic index and
effective radius in these cases are determined by first performing
optimization with two components plus a sky background, followed by a second
iteration reducing to a single Sérsic component but holding the sky parameters
fixed to the best-fit results. These effective values are shown separately in
Table 3. In all cases, the effective radius remains constant, but there is
some very slight modification in Sérsic indices to smaller values (but well
within the uncertainties of the original, one-component values).

We calculate the effective radius of the ALMA 1 mm continuum image by
reconstructing in the source plane for a number of model iterations (as
described in §3.1.3 for estimating magnification. We use the casa (v6.2.1)
imfit algorithm, finding $R_{e}=0.53\pm 0.03\arcsec$. This size is
substantially smaller than the image-plane beam size $1.28\arcsec\times
0.93\arcsec$, due to the lensing magnification. A more careful estimation of
the dust continuum size may be warranted in the future with higher-resolution
submillimeter data (i.e. with unconvolved, forward modeling of the $uv$-plane
visibilities). However, at present, the source-plane effective radius is
smaller than the beam to a degree that such an analysis is not likely to yield
much improved results.

Table 3: Total magnification, spatial extent, and photometry for each filter.

Filter | $\lambda_{\rm eff,rest}$ | Magnification $\mu$ | $R_{e}$ | $R_{e}$ | Sérsic $n$ | $m_{\rm AB}-2.5\log{\mu_{1a}}$
---|---|---|---|---|---|---
| [nm] |  | [″] | [kpc] |  | [mag]
HST/F606W | 180 | 2.3 $\pm$ 0.3 | 0.24 $\pm$ 0.03 | 1.9 $\pm$ 0.2 | 0.23 $\pm$ 0.03 | 24.49 $\pm$ 0.04
HST/F625W | 192 | 2.5 $\pm$ 0.4 | 0.24 $\pm$ 0.04 | 2.0 $\pm$ 0.3 | 0.23 $\pm$ 0.04 | —
HST/F775W | 234 | 2.6 $\pm$ 0.4 | 0.30 $\pm$ 0.05 | 2.5 $\pm$ 0.4 | 0.31 $\pm$ 0.05 | —
HST/F814WaaFor the F814W filter, a detector gap obscures the primary lensing galaxy in the image. While the lensed structure is largely still visible, we exclude this from our source size analysis, as a comparable initial galfit subtraction of the foreground is not possible, and the F775W/F850LP filters can effectively be used as proxy. | 244 | 3.1 $\pm$ 0.9 | — | — | — | 23.87 $\pm$ 0.04
HST/F850LP | 274 | 2.7 $\pm$ 0.3 | 0.35 $\pm$ 0.04 | 2.9 $\pm$ 0.3 | 0.23 $\pm$ 0.03 | —
JWST/F090W | 274 | 3.8 $\pm$ 0.8 | 0.29 $\pm$ 0.06 | 2.4 $\pm$ 0.5 | 0.30 $\pm$ 0.06 | 23.66 $\pm$ 0.01
HST/F105W | 321 | 8.4 $\pm$ 1.5 | — | — | — | —
JWST/F115W | 351 | 6.9 $\pm$ 1.6 | 0.37 $\pm$ 0.09 | 3.1 $\pm$ 0.7 | 0.16 $\pm$ 0.04 | 22.85 $\pm$ 0.01
HST/F125W | 379 | 5.8 $\pm$ 1.5 | — | — | — | —
HST/F140W | 423 | 7.8 $\pm$ 1.5 | — | — | — | —
JWST/F150W | 456 | 7.6 $\pm$ 1.7 | 0.38 $\pm$ 0.08 | 3.1 $\pm$ 0.7 | 0.23 $\pm$ 0.05 | 21.886 $\pm$ 0.005
HST/F160W | 467 | 8.4 $\pm$ 2.3 | — | — | — | —
JWST/F200W | 604 | 7.8 $\pm$ 2.3 | 0.41 $\pm$ 0.12 | 3.4 $\pm$ 1.0 | 0.25 $\pm$ 0.08 | 21.227 $\pm$ 0.004
JWST/F277W | 839 | 8.3 $\pm$ 2.0 | 0.41 $\pm$ 0.10 | 3.4 $\pm$ 0.8 | 0.22 $\pm$ 0.06 | 20.807 $\pm$ 0.004
JWST/F356W | 1084 | 10.9 $\pm$ 2.5 | 0.40 $\pm$ 0.09 | 3.3 $\pm$ 0.7 | 0.42 $\pm$ 0.10 | 20.373 $\pm$ 0.004
(2-comp.$\dagger$$\dagger$footnotemark: ) | — | — | (0.40 $\pm$ 0.09) | (3.3 $\pm$ 0.7) | (0.42 $\pm$ 0.10) | —
JWST/F410M | 1240 | 8.3 $\pm$ 2.0 | 0.39 $\pm$ 0.09 | 3.2 $\pm$ 0.8 | 0.32 $\pm$ 0.08 | 20.202 $\pm$ 0.004
(2-comp.$\dagger$$\dagger$footnotemark: ) | — | — | (0.39 $\pm$ 0.09) | (3.2 $\pm$ 0.8) | (0.32 $\pm$ 0.08) | —
JWST/F444W | 1349 | 9.8 $\pm$ 2.2 | 0.39 $\pm$ 0.09 | 3.2 $\pm$ 0.7 | 0.45 $\pm$ 0.11 | 20.146 $\pm$ 0.004
(2-comp.$\dagger$$\dagger$footnotemark: ) | — | — | (0.39 $\pm$ 0.09) | (3.2 $\pm$ 0.7) | (0.44 $\pm$ 0.10) | —
| $\lambda_{\rm eff,rest}$ | — | — | —— | — | Flux $\mu S_{\nu}$
| [$\mu$m] |  |  |  |  | [mJy]
ALMA/870$\mu$mbbThe signal-to-noise ratio of the Briggs-weighted (robust$=2$, close to natural weighting) 870 $\mu$m image is low, so we do not attempt to derive an effective radius or Sérsic index. | 264 | $\approx 6-9$ | — | — | — | 1.1 $\pm$ 0.7
ALMA/1100$\mu$m | 334 | 9.5 $\pm$ 0.9 | 0.53 $\pm$ 0.03 | 4.4 $\pm$ 0.2 | 0.5 (fixedccThe ALMA data is fit with a 2-dimensional Gaussian, which is nearly equivalent to a $n=0.5$ Sérsic profile.) | 3.8 $\pm$ 0.2

††footnotetext: For the longest-wavelength filters (F356W, F410M, F444W), an
additional Sérsic component is added to account for large residuals in the
galaxy center. As discussed in §3.2, we also calculate effective values for
the combined components. In all cases, the effective radius remains the same,
but there is some slight change in $n$.

Note. — Reported image-plane photometry values are only for image $1a$, and
are not corrected for magnification $\mu_{1a}$, the values for which are
reported in Table 5. For this image, differential magnification is minimized
and approximately consistent from rest-frame far-UV to near-IR. We also
tabulate the effective (or pivot) wavelength of each filter converted to the
rest-frame at $z=2.291$. For the WFC3/IR filters of HST, we do not measure
effective radii or Sérsic index, as superior imaging is available from JWST at
comparable wavelengths.

###  3.3 Photometry of JWST, HST, and ALMA data

In this study, we require photometry of the rest-frame near-UV, optical, and
near-IR regimes in order to confirm the redshift of El Anzuelo, and to measure
some initial basic properties. Image-plane photometry for El Anzuelo was
performed for all JWST/NIRCam filters, in addition to some of the bluer HST
filters, F606W and F814W, which provide constraints on the rest-frame UV not
possible with JWST. We refer the reader to Frye et al. (2023) for more details
on the careful photometric measurements for the entire El Gordo field, which
we make use of here. Briefly, WebbPSF is used to generate PSF models for the
short- and long-wavelength filters, and modeled ACS filter PSFs are used for
the HST images. All images are registered astrometrically onto the Gaia Data
Release 3 reference frame (Gaia Collaboration et al. 2016, 2021). The PSF-
convolved matched-aperture photometry (see Pascale et al. 2022) was collected
using SExtractor (Bertin & Arnouts 1996), with F200W serving as the detection
band. Automatic aperture fluxes from SExtractor are measured using Kron-like
elliptical apertures (Kron 1980). All measurements are PSF-matched to the
reference F200W filter: for longer-wavelength filters with broader PSFs, the
F200W filter is convolved to each respective PSF to determine how much flux
would be lost within the aperture888See
https://www.stsci.edu/~dcoe/ColorPro/color for details.. This flux correction
factor is then applied to each filter. Magnitudes are reported in Table 3,
without correction for lensing magnification, for apertures located at
$(\alpha,\delta)=$ (01h02m49 $\fs$169, $-49^{d}$15m07 $\fs$99)999For
comparison, we also include photometry from Frye et al. (2023) at an
independent aperture located at (01h02m49 $\fs$202, $-49^{h}$15m04$\fs$32),
part of the highly-distorted northern portion of the arc, as shown in the
right panel of Fig. 7. This is labeled Aperture 2 in Fig. 3., shown as
Aperture 1 in Fig. 3. We also measure the flux of the 870 $\mu$m and 1 mm
continuum with blobcat (Hales et al. 2012), a flood-filling detection
algorithm designed for radio images.

For this initial study, we opt to use only image-plane photometry for fitting
the SED, as source-plane photometry is plagued by a number of issues. For
example, the magnification gradient of El Anzuelo is rather steep in the
region of the background DSFG, which means that the source-plane PSF varies
significantly over the extent of the source (although this primarily impacts
the highest-magnification multiple images towards the northeast). This is
because the regular, consistent PSF for the observed data becomes distorted
and magnified inconsistently in the source plane. This makes it very difficult
(or impossible) to perform PSF-convolved aperture photometry directly on the
demagnified source-plane emission. Moreover, it is much more difficult to
incorporate any local sky subtraction in the source plane, as foreground
objects must be masked before ray-tracing with the lens model (and doing this
for too large an area becomes computationally prohibitive).

###  3.4 Inner vs. outer disk decomposition

To analyze the SED of El Anzuelo for the inner bulge (or central kiloparsec)
relative to the outer disk, we simultaneously de-magnify the image-plane
photometry and estimate the bulge-to-total fraction of light ($B/T$) in the
source-plane using photutils (Bradley et al. 2022). Before performing this
analysis, we first convolve all filters to match the broadest PSF of the F444W
filter, to ensure that our spatially-resolved color analysis is not influenced
by the PSF. By correcting the image-plane photometry rather than performing
the measurements directly on the source-plane, we can avoid the effects of a
varying PSF and poor sky subtraction.

To determine a sensible threshold for the inner bulge versus outer disk, we
follow the method employed by Cutler et al. (2023) to determine appropriate
sizes for a fixed aperture and annulus. Cutler et al. used $z_{\rm F850LP}$
and $H_{\rm F160W}$ in order to span the 4000Å-break at $z\sim 2.3$ as an age
indicator. Here we instead use F277W and F444W (spanning the wavelength at
which point the light profile begins to steepen for El Anzuelo, which is
suggestive of a bulge/spheroidal component being assembled). The [2.8 - 4.4]
color is measured in apertures of varying radius from 0.03$\arcsec$ to
0.6$\arcsec$, in over-sampled increments of 4 milliarcseconds. The smallest
aperture that encloses the reddest emission (i.e. maximizes [2.8 - 4.4]) is
chosen, under the condition that the next local minimum at a larger radius is
less than 1% of the local maximum’s color (to avoid secondary peaks). In our
case, there is only one local maximum, which is at 0.13$\arcsec$, or
approximately 1.1 kpc. While this is comparable to the scale of the NIRCam PSF
FWHM in F444W (Rigby et al. 2023), lensing magnification makes it so that this
inner region is safely resolved in all filters. The fraction of light inside
this aperture to the total flux of the disk (or bulge-to-total fraction,
$B/T$) is given for each filter in Table 5. These fractions are then applied
to the de-magnified photometry, and eazy is run for the inner, outer, and
total disk. Coincidentally, this inner vs. outer threshold of $\sim$1 kpc is
consistent with that used often to measure the central surface mass density,
as an alternative to a Sérsic index (e.g., Cheung et al. 2012; Ji & Giavalisco
2022).

The available photometric coverage helps to break the degeneracy between
reddening from stellar population age vs. from dust attenuation, in particular
through the well-studied $UVJ$ diagram (e.g., Labbé et al. 2005; Wuyts et al.
2007; Patel et al. 2012). The F115W, F150W, and F200W filters together offer
strong constraints on the rest-frame $U$ and $V$ fluxes, while F356W now
offers comparable spatial-resolution coverage of rest-frame $J$ flux for the
first time (Miller et al. 2022). As discussed by Leja et al. (2019), however,
SED-fitting is often more precise and accurate than color-color diagrams in
identifying properties. This may be especially true for spatially-resolved
$UVJ$ diagrams, which have only recently become feasible at $z\sim 2$ and have
not yet been fully calibrated.

Regarding the demagnification itself, we carry this out in two different ways.
For the first approach, we correct image-plane fluxes by the magnification
that we derive on a filter-by-filter basis, using the values given in Table 5.
However, given the expected correlation in magnifications for adjacent
filters, this may not be the most physically-motivated approach. Instead, it
may make more sense to adopt a regime-specific magnification (i.e. separate
values for UV, for optical, and for near-IR). Yet, for image $1a$, we observe
very little variation in magnification for all filters, so we instead adopt a
median magnification that we apply to all filters for our second approach.
Moreover, we consider that the uncertainties for each filter might be overly-
conservative estimates, as the empirical scatter in magnification from far-UV
to near-IR is notably smaller than any of the quoted $1\sigma$ uncertainties.

For this approach, we adopt the dispersion in $\mu_{1a}$ as the $1\sigma$
uncertainty in magnification. While this may be an under-estimate of any
systematic uncertainty, we consider it likely that these would essentially
affect all filters equally, which would modify only the amplitude of the SED
and not its shape. Lastly, we note that the derived properties in Table 4 are
consistent between the two methods within uncertainties, so the difference has
minimal impact on our interpretations.

Figure 6:  The assumed CO(3-2) line detection for El Anzuelo (top) from ALMA
program 2015.1.01187.S (Basu et al. 2016) at observed frequency 105.07 GHz
($z=2.291$). While the target is near the edge of the primary beam of this
pointing, an independent detection of the same line from a separate pointing
(program 2017.1.01621.S, bottom) confirms that it is very likely not a
spurious artifact (although this pointing is even further from the primary
beam center).

Figure 7:  SED-fitting results with eazy for image $1a$ (left) and the
northernmost part of the arc where there is significant distortion (right),
which we use as independent verification of the photometric redshift estimate.
In the left part of each panel, the SED templates at the best-fit and
spectroscopically-confirmed redshifts (in blue and red, respectively) are
shown: $z\approx 2.34$ (left) and $z\approx 2.28$ (right) vs. spectroscopic
$z=2.291$. The measured photometry points are shown as black squares, while
the red circles show the expected fluxes (including continuum and line
emission) given the filter throughput. In the right part of each panel, the
probability distribution function of redshift given the SED fitting is shown,
with the spectroscopic redshift shown as a red vertical line. In both cases,
the redshift solution is extremely narrow (with uncertainties more typical of
spectroscopic redshifts), thanks to the careful sampling of the rest-frame UV
through near-IR spectrum. In particular, the 4000Å break and/or Balmer break
(at $\sim$3600Å) is well-constrained with the F115W and F150W filters, as the
break falls right in the small gap between these wideband filters.

###  3.5 Source redshift confirmation and measurement of basic properties with
eazy

Cheng et al. (2023) included El Anzuelo in their preliminary analysis of
objects with joint ALMA/JWST detections. Their preliminary photometry with
only JWST NIRCam filters resulted in photometric redshift solutions of
$z=3.54\pm 0.15$ and $z=3.03\pm 0.06$ for the primary southwest and northeast
lensed components. However, this photometry did not specifically account for
the contribution of light from the foreground elliptical galaxies. Instead,
Cheng et al. employed small 0.6″ apertures centered on the peak pixel in F444W
in order to minimize the contamination. As there is a spatial offset in the
continuum structure between the short- and long-wavelength filters of NIRCam,
this had the effect of underestimating the flux below 2 $\mu$m, thereby
biasing to higher-redshift solutions.

For the SED-fitting and photometric-redshift estimation in this work, we use
eazy (Brammer et al. 2008; Brammer 2021) to optimize the rest-frame UV through
near-IR SED ($0.18-1.35\mu$m; observed-frame $0.59-4.44\mu$m). We note,
however, the need for future work using software with a complex treatment of
star-formation histories, such as prospector (Johnson et al. 2021). We used
the tweak_fsps_QSF_12_v3.param eazy templates, with a uniform prior on the
redshift from $z=0.01-6$, which encompass a range of realistic star formation
histories and levels of dust attenuation (up to $A_{V}\sim 3$). We test
several other template sets, including the eazy_v1.x.param files (which
include a dusty galaxy template), but find minimal observable differences in
the results.

The photometric redshift for El Anzuelo is remarkably precise ($z=2.34\pm
0.01$), although the redshift uncertainty is very likely underestimated, since
it disagrees with the spectroscopic-redshift, $z=2.291$, as we discuss in the
next paragraph. We interpret the reason for this precision to be the
fortuitous alignment of the Balmer break at $\lambda=1.2\mu$m and the 4000Å
break at $\lambda=1.3\mu$m, such that F115W exclusively samples continuum
short of the 4000Å break, and F150W samples only continuum longward of the
break.

With this precise prediction for the redshift, we re-examined ALMA archival
data that might provide coverage of bright lines at this redshift. The 3 mm
continuum data (program 2015.1.01187.S, PI: K. Basu) has low spectral
resolution and coarse angular resolution; moreover, the location of El Anzuelo
falls near the edge of the primary beam, where the sensitivity is degraded.
However, we are still able to detect what we infer to be the CO(3-2)
transition at 105.07 GHz with high significance at the center of the bandpass
(see Fig. 6). Additionally, a separate pointing (also at 3 mm; program
2017.1.01621.S, PI: K. Basu) has a different phase center, but the line is
again detected (despite falling once again near the very edge of the usable
primary beam response). We consider this to be a secure detection from these 3
independent sets of observations, but we will pursue follow-up observations in
the future to offer further confirmation. It may also be possible that the
detected line is CO(4-3), which would indicate a redshift solution of
$z=3.388$, although this is in worse agreement with the photometric evidence.
For this work, we adopt $z=2.291$ as the most likely solution, and simply note
that the angular-to-physical size ratio changes only modestly from 8.21 kpc
arcsec-1 at $z=$2.291 to 7.41 kpc arcsec-1 at $z$=3.388. Moreover, most
results discussed in this work pertain to relative sizes as a function of
wavelength, so this has little impact.

Unfortunately, it is not possible to derive geometric redshifts based on the
lens model (see, e.g., §4.3 of Diego et al. 2023). This is because the lens
model is only constrained by El Anzuelo itself, and the background redshift is
more or less degenerate with the mass of the foreground lens. One could
hypothetically use Fundamental Plane relations (e.g. Dressler et al. 1987) and
the spectroscopic redshift of the cluster members to estimate the mass
enclosed within the Einstein radius of the lens (thereby allowing for a
determination of the background redshift, as this radius depends only on the
foreground mass and the redshift geometry of lens and source planes). However,
given that the lens consists of a group of galaxies and not a single massive
elliptical, the dark matter distribution is more uncertain, and this method
may be too tenuous.

We adopt henceforth a redshift of $z$=2.291, and re-run eazy with the redshift
held fixed, in order to estimate more accurate uncertainties on the derived
properties, which we discuss in §4 and present in Table 4. In short, the star-
formation rate (SFR) and stellar mass ($M_{\star}$) are inferred using only
the far-UV-to-near-IR photometry, as the angular resolution of the far-IR
observations is much coarser, and thus insufficient for a spatially-resolved
SED decomposition. In the future, we intend to incorporate far-IR information
at a more comparable angular resolution with a more flexible approach to SED
modeling. For this present analysis, properties derived from the SED are used
primarily to support inferences from the surface brightness distributions.

##  4 Discussion

Recently-collected JWST/NIRCam imaging reveals El Anzuelo to be an excellent
candidate for follow-up studies, especially with targeted, higher-resolution
interferometric radio/millimeter imaging and spectroscopy. This will allow for
a more complete analysis of the interstellar medium driving the active star
formation. For this initial case study, we focus primarily on the distribution
of light in the rest-frame UV, optical, near-infrared, and far-infrared, in
order to understand the distributions of new star formation and existing
stellar mass from older populations. As part of this, we must take into
account the fact that the effect of gravitational lensing is not nearly
spatially uniform across El Anzuelo.

###  4.1 Large differential magnification between rest-frame UV, near-IR, and
far-IR

Figure 8:  Derived magnification $\mu_{\rm total}$ is shown for each HST and
JWST filter available for El Anzuelo. Corresponding filter throughput curves
are shown at the bottom. A dramatic rise in magnification with wavelength is
due to a steepening in the light profile, as shown also with the RGB inset
images showing rest-frame UV, optical, and near-IR images of the source-plane
reconstruction. There is a notable discontinuity in magnification around
$1\mu$m, which is due primarily to the lack of signal in the central core,
where magnification is greatest. This is likely due to some combination of a
dust attenuation and an age gradient, especially given the location of the
4000Å break, redshifted to $1.3\mu$m for $z=2.291$.

Figure 9:  Variation in the reconstructed source-plane light profile for the
HST/ACS and JWST/NIRCam filters. Top: Best-fit Sérsic index $n$ determined
with galfit for each filter. For the three longest-wavelength filters, excess
model residuals near the galaxy nucleus suggests the need for a 2-component
(bulge and disk) model. Here, an effective Sérsic index is shown, as discussed
in §3.2. Filters are denoted by color as in Fig. 8. As the uncertainties
reported by galfit are likely underestimations, we scale them up to match the
fractional uncertainties in magnification for each filter (Table 3). Bottom:
As above, but showing effective radii for the Sérsic profiles for each filter,
in angular and physical units (at the redshift of the source, $z=2.291$).

Figure 10:  Major axis 1-dimensional surface brightness profiles of the
source-plane reconstruction for each filter. The large inferred inclination,
clear asymmetry, and lack of an obvious center for all filters make this a
more robust approach than a radial profile, as tested by Devour & Bell (2017,
2019). Left: Filter-by-filter surface brightness profiles along the major axis
(as determined by the median position angle from galfit models for all
filters). Each bin along the profile, which shows the sum of pixels above a
$2\sigma$-clipped threshold, is extracted from a 1.5$\arcsec$-wide box
perpendicular to the major axis. The profile for each filter (each convolved
to match the F444W PSF) is normalized between 0 and 1 and has arbitrary flux
units, and they are shifted vertically for the purpose of legibility. Bins are
spaced at intervals of $\approx 0.045\arcsec$ along the $2\arcsec$-long slit,
comparable to the pixel size. The thickness of the curve indicates signal-to-
noise, which is determined from the signal-to-noise ratio of the 95th
percentile value of pixels in each bin. Differing pixel scales and depths
likely play a role in discrepancies between adjacent filters. Vertical dashed
lines show the smallest effective radius (for F606W). Center: An RGB image
showing all 8 NIRCam filters (red: F356W + F410M + F444W, green: F200W +
F277W, blue: F090W + F115W + F150W). The white box shows the slit used to
compute the surface brightness profiles, with the middle line showing the
major axis orientation. The beams in the lower left show the FWHM and
orientation of representative source-plane PSFs for F115W (cyan) and F444W
(red). These are created by ray-tracing a model PSF placed at the location of
the galaxy centroid in the less-magnified counter-image ($1a$) into the
source-plane, which has the least spatial variation in PSF. While the more
magnified images can offer better resolution on their own, this source-plane
reconstruction combines all images together, effectively smearing the
structure to the resolution of the counter-image. Right: Inferred surface
brightness profiles in the same manner as the left panel, but applied to our
galfit-modeled light profiles. Dashed curves for the three reddest NIRCAM
filters show the 2-component models, which shift the profiles slightly.

While gravitational lensing is nominally achromatic (i.e., all wavelengths
originating from the same location of the source-plane are deflected in the
same manner), differential magnification can arise due to different source-
plane distributions of light at different wavelengths. This effect is easily
observable for El Anzuelo, as the redder, central core of the galaxy (which is
significantly more prominent beyond 2 $\mu$m) is precisely aligned with some
of the highest-magnification regions of the source-plane near the caustic
curves (as we discuss in more detail in §4.3). As the total magnification for
each filter is a flux-weighted average of the magnification map, the long-
wavelength NIRCam filters with this prominent core are amplified much more
strongly.

In Fig. 8, we illustrate this trend towards larger magnifications at longer
wavelengths between 0.5 and 5 $\mu$m. Shorter than 1 $\mu$m (observed-frame),
the emission is confined to an outer region of the disk, with a substantially
different morphology from that of the long-wavelength bands of NIRCam. Given
this, it is not surprising that the magnification for rest-frame UV differs
from that for rest-frame near-IR. The magnitude of this difference is notable:
the near-IR stellar continuum is amplified by a magnification $\mu\sim 8-12$,
while the rest-frame UV is amplified by only $\mu\sim 2-3$. This factor of
$\sim 4$ difference in magnification appears to be driven primarily by the
alignment of the red galaxy core with two caustic cusps in very close
proximity (see Fig. 3). As total magnification is weighted by flux, the lack
of signal in this region for shorter-wavelength filters dilutes their
magnification and weights more highly the lower magnification ($\mu\sim 2$)
from the outer regions of the galaxy. One can readily observe this in Figs. 2
and 3, where the northeast image is dominated by emission at wavelengths
longer than $\sim$1.5$\mu$m. This is effectively a manifestation of the well-
known lensing size bias, in which more compact (or more highly concentrated)
objects are often magnified more strongly (e.g., Hezaveh et al. 2012;
Robertson et al. 2020).

At present, there are relatively few robust measurements of differential
magnification to which we can compare. In simulations of lensed DSFG samples,
Serjeant (2012) predicted differential magnifications between older to younger
stellar populations of log[$\mu_{\star}$/$\mu_{\rm SFR}$] = $0\pm 0.26$ in the
case of $\mu>10$, or log[$\mu_{\star}$/$\mu_{\rm SFR}$] = $0\pm 0.14$ for
modest magnifications $2<\mu<5$. For observed lensed DSFGs, Calanog et al.
(2014) found a typical differential magnification between near-IR and far-IR
(880 $\mu$m) of $\mu_{\rm FIR}$/$\mu_{\rm NIR}\sim 1.5$ (i.e. +0.18 dex).
Thus, the $\sim$0.6 dex difference in our case clearly represents an extreme
case of differential magnification. By assuming a median magnification of
$\mu\sim 7$ for all filters, for example, this would have the effect of
underestimating the UV flux by a factor of $\sim 2-3$, and overestimating the
IR flux by a factor of $\sim 1.5$, the combined effect of which could severely
impact the inferred galaxy properties. However, this is only the case for the
full galaxy-integrated photometry. For the aperture we use at the location of
image $1a$, as we discuss in §3.4, the median magnification for all HST/JWST
filters is $\mu_{1a,{\rm med}}=2.7\pm 0.2$. This smaller uncertainty reflects
the empirical scatter as a function of wavelength, but may underestimate
systematic uncertainties that affect all filters.

###  4.2 A dust-obscured nucleus or inside-out growth—or both?

As shown in Figs. 9 and 10, there is actually not a high degree of variation
in the light profile from $\sim 1.5-4.4\mu$m, all of which is longward of the
4000Å-break (observed-frame $1.3~{}\mu$m). At first, this consistency in light
profile from rest-frame near-UV to near-IR might appear to rule out large
gradients in dust attenuation, which would push the half-light radii of
shorter wavelengths to larger values101010This is because the expected central
peak in surface brightness becomes increasingly suppressed for bluer filters
when there is a radially-decreasing profile of dust attenuation, such that an
increasingly large aperture is required to contain half of the galaxy’s
observed flux. (Nelson et al. 2016a; Bondi et al. 2018; Tacchella et al. 2018;
Suess et al. 2019; Popping et al. 2022). This runs somewhat contrary to the
expectation that IR-luminous galaxies would have high concentrations of dust
(and high IR luminosity surface densities) within the central kiloparsec (as
for ultra-luminous infrared galaxies, or LIRGs and ULIRGs, e.g., Sanders &
Mirabel 1996; Lonsdale et al. 2006; Armus et al. 2009; Rujopakarn et al.
2011). As we point out later in this section, this could also be the result of
inclination effects. Recently, Miller et al. (2022) found that 70% of star-
forming galaxies at $z\sim 2$ demonstrated $U-V$/$V-J$ color gradients caused
primarily by gradients in dust attenuation (also Liu et al. 2017; Wang et al.
2017). For the same sample, 23% of the objects showed color gradients that
were in line with gradients in stellar population properties (such as specific
star formation rate).

At wavelengths shorter than the Balmer and 4000$\AA$ breaks ($<1.2\mu$m
observed-frame), we find that the effective radius decreases gradually to
$0.24\arcsec$, about 40% smaller than the sizes in the rest-frame near-IR. In
an analysis of the resolved dust continuum of star-forming galaxies simulated
with FIRE (Feedback in Realistic Environments; Hopkins et al. 2014, 2018),
Cochrane et al. (2019) showed that rest-frame UV was significantly correlated
spatially with both cool, dense gas (i.e. regions of active star formation)
and holes in the dust distribution. In both cases, this leads to a very clumpy
spatial distribution, often offset from the near-IR structure. Similarly, Ma
et al. (2018) found that these UV-bright clumps (which can dominate the
emission) are likely poor, biased tracers of stellar mass. On the other hand,
while only for $z>7$ galaxies, Yang et al. (2022) found only a modest decrease
in size from rest-frame optical to UV. Shibuya et al. (2015) likewise found
close agreement in median UV and optical sizes for $z\sim 1-2$ and stellar
mass $\log(M_{\star}/M_{\odot})=9-11$, although with large scatter. However,
Ma et al. (2018) also point out that parametric size measurements for very
clumpy, irregular galaxies will break down and will primarily capture only the
separations between individual bright clumps. Given the patchy, nearly non-
contiguous structure of El Anzuelo in the UV, we consider this to be a likely
explanation for the suppressed UV sizes.

In Table 4, we report basic (magnification-corrected) properties derived from
fitting the SED of El Anzuelo as a whole, while also repeating this for the
inner and outer disk regions separately. In correcting the photometry for
lensing magnification, we compare two options: $a)$ applying filter-by-filter
magnification corrections, versus $b)$ applying a median magnification factor,
for which uncertainty is given by the empirical scatter for the set of
filters, as discussed in §3.4. For our purposes, we primarily use the median
magnification, which is more physically motivated in that magnification is
unlikely to change significantly between adjacent filters for image $1a$
(given the extensive photometric sampling available, and given the minimal
magnification gradient in this region of the image plane).

A similar analysis of an un-lensed galaxy would likely be significantly
impeded by the PSF, but the magnification of El Anzuelo means that this
central region is resolved by NIRCam (such that the observed color gradient is
not driven only by variation in PSF with wavelength). The characterization of
PSFs for NIRCam is still ongoing; so far it has been found that the PSFs can
vary spatially across the field-of-view (e.g. Zhuang & Shen 2023) and
temporally over timescales of a month (e.g. Nardiello et al. 2022), both of
order $\sim 5\%$. For this reason, it is desirable to examine and calibrate
color gradients in objects with angular sizes much larger than the
instrumental resolution, such as lensed galaxies. While gravitational lens
modeling introduces its own systematic uncertainties, some quantities (such as
color variations) can actually be interpreted in the image plane, requiring no
assumptions about the intervening mass.

When we perform SED fitting on the inner vs. outer disk in the prescribed
manner, we find that the implied attenuation is larger in the galaxy center
than in the outskirts ($A_{V}=2.0$ vs. $A_{V}=1.4$). We also observe that the
half-light radius for the dust continuum is $R_{e}\approx 0.5\arcsec$, which
is reasonably consistent with the UV/optical sizes. This may be due in part to
the coarser resolution of the 1 mm ALMA continuum imaging relative to JWST,
but there is not an obvious factor of $2-4\times$ difference in effective
radius for the UV vs. far-IR, which has been observed in many other studies of
DSFGs (including Barro et al. 2013, 2016; Calanog et al. 2014; Chen et al.
2015; Ikarashi et al. 2015; Simpson et al. 2015; Tadaki et al. 2015; Hodge et
al. 2016; Lutz et al. 2016; Massardi et al. 2018; Talia et al. 2018; Gullberg
et al. 2019; Tadaki et al. 2020; Pantoni et al. 2021, and references therein,
for a sampling). Moreover, Lang et al. (2019) found that the rest-frame far-IR
continuum emission was not just more compact than the rest-frame optical
light, but also more compact than the stellar mass distribution itself
(described by the half-mass radius), indicating a higher central concentration
of star formation relative to existing stellar mass.

The central light profile for $\lambda=3-5~{}\mu$m (Sérsic $n\sim 0.4$) is
also steeper than shorter wavelengths (Fig. 9), implying a greater
concentration in stellar mass than in the distribution of more recent star
formation (probed by UV continuum). We caution still that these values might
be subject to systematic effects arising from applying galfit to the source
plane. Given the highly asymmetric nature of the disk, with the red central
bulge appearing offset to one end, we also examine the surface brightness
profiles in a non-parametric fashion in Fig. 10. This approach is similar to
the inclination-independent linear surface brightness profile developed by
(Devour & Bell 2017, 2019). We extract the 1-dimensional surface brightness
profiles along a slit $2\arcsec$ in length for each filter (oriented at the
median position angle from galfit and centered at the F444W peak, with a width
of $1.5\arcsec$ perpendicular to the major axis). Each image was first
convolved in the image-plane with a kernel to match the PSF of the F444W
filter, before reconstructing in the source plane. Pixels in each bin are
sigma-clipped at a conservative $2\sigma$-level and then summed. The curve
thickness denotes signal-to-noise, determined by taking the 95th-percentile
value of the distribution of S/N pixel values in each bin (which indirectly
incorporates the number of pixels above the noise threshold). In the right
panel of Fig. 10, we repeat this process exactly, but for the galfit modeling
results (including the 2-component models for F356W, F410M, and F444W) to show
the inferred, modeled profiles. We can visualize the increase in galaxy size
from $\sim 0.6-1.2~{}\mu$m (observed), but relative consistency at longer
wavelengths, although the peaks are not co-spatial.

As discussed previously, a number of issues plague the measurement of
photometry in the source plane. These include the need for a proper accounting
of the sky background (which becomes distorted and limited in scope after ray-
tracing the lensed structure to the source plane). In capturing one-
dimensional surface brightness profiles along a slit in Fig. 10, we can
examine consistency between filters while re-normalizing the profiles so that
the sky background does not need to be removed. Moreover, while the source-
plane PSF may vary across the object, this effect may be mitigated by
averaging within each bin along the length of the slit. Variation in the PSF
from filter to filter can be observed in the smoothness of each profile, but
the structure is resolved for all filters and so the overall profile shape is
not hindered by the PSF.

Figure 11:  SED-fitting results for image $1a$, as with Fig. 7, but fluxes
from Table 4 are de-magnified using the median magnification $\mu_{\rm
1a}=2.7\pm 0.2$ (see Table 5). The SED is spatially decomposed into the total
disk (top), inner $\sim$1 kpc (middle), and outer disk (bottom) according to
the method described in §3.4.

The spatial offsets we observe might be due to a patchy and uneven dust
distribution (e.g., Liu et al. 2013; Casey et al. 2014b; Popping et al. 2017;
Narayanan et al. 2018; Ma et al. 2019). Likewise, the western portion of the
source appearing bluer than the eastern side (see e.g., middle panel of Fig.
10) could be the result of a nearly edge-on inclination ($i=63^{\circ}\pm
7^{\circ}$; §3.2) coupled with a dust lane. This interpretation would indicate
that the western portion is the far side of the disk, for which less stellar
light is attenuated (e.g., Hubble 1943; de Vaucouleurs 1958). In this
scenario, it is possible that the size measured in the bluer filters is
dominated by the light profile of the far side of the disk. However, given
that the inclination-independent 1-dimensional surface brightness profiles in
Fig. 10 are consistent in width beyond $\sim$1.2 $\mu$m (rest-frame $\sim$350
$\mu$m), this explanation appears to be unlikely.

Another interpretation might be that the spatial offset is the result of dust
driven off the plane of the disk by being entrained within outflows and strong
stellar winds, as seen locally in the starburst galaxy M82 (Sanders & Balamore
1971; Alton et al. 1999; Roussel et al. 2010; Hutton et al. 2014). Local LIRGs
and ULIRGs also contain a significant fraction of active galactic nuclei
(AGN), of order 20% (Goulding et al. 2010; Petric et al. 2011) and 40-60%
(Armus et al. 2007; Desai et al. 2007; Veilleux et al. 2009), respectively.
Similarly, this could also be reminiscent of star formation triggered by jet-
driven shocks compressing the interstellar medium, as in the case of Centaurus
A (Crockett et al. 2012). Additionally, such color asymmetries and spatial
offsets can be seen in the local Universe in (ultra-)luminous infrared
galaxies (LIRGs and ULIRGs). One particularly illustrative example is the mid-
stage merger, VV 114 (e.g., Knop et al. 1994; Yun et al. 1994; Iono et al.
2004; Saito et al. 2015; Evans et al. 2022), for which the gas kinematics are
dominated by the remnants of orbiting disks. In this case, a dust lane from
one of the disk remnants reddens the other component, which may be quite
similar to what we are seeing with El Anzuelo.

####  4.2.1 SED-derived SFR and $M_{\star}$

In order to have an effective radius consistent with that of near-IR, but with
a shallower core from dust attenuation, the rest-frame UV light would need to
truncate more steeply at large radii. This itself would suggest an age
gradient with massive, more short-lived and more UV-emitting stellar
populations in the center than in the outskirts. A radially-increasing age
gradient would also manifest as a truncation in the extent of UV emission,
from the lack of hot, young stars in the galaxy outskirts and from metal
absorption in the cooler stars that are present.

Yet, this presumed age gradient is not reflected in the SED-derived specific
star formation rate (sSFR $\equiv{\rm SFR}/M_{\star}$), which is greater (i.e.
closer to a starburst) for the outer disk relative to the inner disk by a
factor of 3 ($1.5\pm 0.9~{}{\rm Gyr}^{-1}$ versus $0.51\pm 0.04~{}{\rm
Gyr}^{-1}$). The respective SED fits are shown in Fig. 11; see also Tacchella
et al. (2018), for example. At the redshift of El Anzuelo, the specific SFR at
the main sequence is ${\rm sSFR}_{\rm MS}\sim 2~{}{\rm Gyr}^{-1}$ (Speagle et
al. 2014). We find the total sSFR for the galaxy to be consistent with the
main sequence at $1.7\pm 0.4~{}{\rm Gyr}^{-1}$, but we caution that SFR
derived from the UV-to-NIR SED can be underestimated without incorporating the
far-IR that probes the dust SED (e.g., Elbaz et al. 2018). Likewise, stellar
masses from the UV-to-NIR SED alone can also be underestimated by up to
$\sim$0.3 dex for massive DSFGs ($M_{\star}\gtrsim 10^{10}~{}M_{\odot}$;
Battisti et al. 2019). Unfortunately, the expected differential magnification
and the low angular resolution of existing millimeter imaging (or low
sensitivity, for 870 $\mu$m) makes a more robust measurement unfeasible at
present. Still, it is clear that this object does not match the hyper-luminous
regime ($L_{IR}>10^{13}~{}L_{\odot}$) seen for the most extreme DSFGs and the
AGN-driven, hot Dust-Obscured Galaxies (DOGs; Dey et al. 2008; Pope et al.
2008; Eisenhardt et al. 2012; Wu et al. 2012, 2014; Tsai et al. 2015; Penney
et al. 2020) identified with the Wide-field IR Survey Explorer (WISE; Wright
et al. 2010) or the Spitzer Space Telescope (Werner et al. 2004).

To test the robustness of the rising specific star-formation rate that we find
in the source plane, especially given the influence of the lensing-distorted
PSF, we repeat this analysis in the image plane. We use a fixed elliptical
aperture ($0.3\arcsec\times 0.2\arcsec$, PA$=50^{\circ}$) centered on the
F444W peak for the inner component, and an elliptical annulus (outer extent
$1.5\arcsec\times 1.0\arcsec$, PA$=50^{\circ}$) for the outer component, and
perform PSF-matched photometry on images convolved to the F444W PSF. These
elliptical apertures roughly account for the distortion by lensing of the
circular apertures used in the source plane when mapped to the image plane.
Since specific star-formation rate is the ratio of SFR to $M_{\star}$, both of
which are approximately equally scaled by lensing magnification $\mu$, one
would expect the sSFR to remain consistent in the image plane. With this
approach, we find an inner sSFR $=0.44\pm 0.03$ Gyr-1 and outer sSFR $=1.3\pm
0.2$ Gyr-1, both slightly lower than the values measured in the source-plane,
but in agreement within the $1\sigma$ uncertainties. The lower sSFR estimates
could also simply be the result of the elliptical image-plane apertures not
matching the exact source-plane regions covered by the circular apertures.
Thus, since the rise in sSFR towards outer radii that we observe in the source
plane is also recovered through PSF-matched image-plane photometry, we have
reason to believe this trend is real.

Similarly, we find relative consistency in the surface density of star
formation111111$\Sigma_{\rm SFR}\equiv{\rm SFR}/A_{d}$, where $A_{d}$ is the
disk surface area, $2\pi(R_{\rm out}^{2}-R_{\rm in}^{2})$ for the inner and
outer radial bounds of each region, $R_{\rm in}$ and $R_{\rm out}$.
($\Sigma_{\rm SFR}$) between the inner 1 kpc and outer disk, at $\Sigma_{\rm
SFR}=0.88\pm 0.03~{}M_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$, versus
$0.7\pm 0.3~{}M_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$, respectively. These
are also in line with dust-corrected SFR surface densities found by Tacchella
et al. (2018) for $z$$\sim$2 galaxies ($10<\log[M_{\star}/M_{\odot}]<11$) on
the star-forming main sequence For context, Heckman (2001) suggest that
$0.1~{}M_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$ is the threshold beyond
which starburst-driven galactic winds become ubiquitous (see also Sharma et
al. 2016).

In the context of the inside-out growth framework of stellar mass assembly (an
extensive set of literature, including Carrasco et al. 2010; van Dokkum et al.
2010, 2014; Oser et al. 2010; La Barbera et al. 2012; Nelson et al. 2012,
2016b; van der Wel et al. 2014b; Wuyts et al. 2013; Dimauro et al. 2022),
these findings support a scenario where El Anzuelo has begun quenching its
star formation in an inner bulge. With a mass $M_{\star}=[1.3\pm 0.1]\times
10^{10}~{}M_{\odot}$, this is roughly a quarter of the total stellar mass,
$M_{\star}=[4.9\pm 1.4]\times 10^{10}~{}M_{\odot}$. This is in line with the
observation that galaxies above $M_{\star}=10^{10}~{}M_{\odot}$ seem to have
an increasing percentage of active star-forming regions with increasing radius
(e.g., Rowlands et al. 2018), interpreted as a signpost of inside-out
quenching. It may also be the case that feedback from a (likely buried) AGN is
contributing to quenching in the galaxy center, but this is challenging to
confirm directly (e.g., review by Fabian 2012; also Smethurst et al. 2016,
2017; Luo et al. 2022). Recent findings have also shown that the most strongly
star-forming systems at high-$z$ can host active star formation widely
distributed on kiloparsec scales throughout the disk (e.g., Bussmann et al.
2012, 2013; Hodge et al. 2019; Kamieneski et al. 2023). These often-clumpy
sites of star formation may be driven by gravitational instabilities in the
outer gas disk, and may eventually migrate radially inwards and add to the
growth of a central bulge (e.g., Ceverino et al. 2010, 2012; Wuyts et al.
2012). This results in morphological quenching (Martig et al. 2009), by which
the bulge stabilizes gas in the galaxy from gravitational collapse.

Still, central starbursts are not uncommon (e.g., Ellison et al. 2018, 2020;
Wang et al. 2019a), and galaxies above the main sequence are more likely to
have amplified star formation rate surface densities in their central region
(out to $\sim$0.5 $R_{e}$). Likewise, the ratio of inner-to-outer sSFR is
expected to have a positive correlation with distance above the main sequence
(Tacchella et al. 2016a, b, 2018; Zolotov et al. 2015). In a decomposition of
the bulges vs. disks of main sequence galaxies (also at $z$$\sim$2.3), Cutler
et al. (2023) found that central regions had star formation histories (SFH)
dominated by recent ($<100$ Myr) bursts, whereas outer portions of the disk
showed much steadier SFHs. They suggested that this could be driven by
dissipative gas compaction events driven by disk instabilities (e.g., Dekel &
Burkert 2014) or accretion of star-forming clumps (Dekel et al. 2009; van
Dokkum et al. 2013). In the case of El Anzuelo, major or minor mergers might
have been responsible for tidal disruption of gas towards the center (Bournaud
et al. 2011), as is the case with ULIRGs in the local Universe. However,
unlike local ULIRGs, strongly star-forming systems at Cosmic Noon are a much
more heterogeneous sample beyond only major mergers (e.g., Hayward et al.
2011). We presume that El Anzuelo is not yet a dynamically-settled disk,
although follow-up kinematic observations would provide helpful context as to
its level of rotational support.

We now revisit the observation that rest-frame near-UV sizes are similar to
those at near-IR, which seems to stand in contrast to the simple picture of
widely-distributed star formation and a central concentration of greater dust
attenuation (which would actually broaden the UV sizes relative to longer
wavelengths). Perhaps the simplest observation comes from the shifting peak
along the major axis from far-UV to near-IR, seen in Fig. 10, with offsets of
order $\sim$0.3$\arcsec$, or $\sim$2.5 kpc in physical units. This separation
is substantial, as it is comparable to (or even exceeds) the effective radius
in the far-UV (Table 3). This suggests that the far-UV and near-IR are
primarily tracing spatially-distinct components—similar to the common offsets
in UV versus far-IR—as might be the case for the late stages of a galaxy
merger like VV 114. We conclude that it’s likely that the complex distribution
of dust and ongoing star formation in El Anzuelo is not easily characterized
by a 1-dimensional surface brightness profile on its own, and a full
spatially-resolved mapping of dust, gas, and star formation (e.g., through
integral field spectroscopy) is a better-suited observational approach.

Table 4: Basic galaxy properties derived with eazy SED-fitting, including an
inner vs. outer decomposition.

Component | SFR | $M_{\star}$ | sSFR $\equiv$ SFR / $M_{\star}$ | $\Sigma_{\rm SFR}$ | $\Sigma_{M_{\star}}$ | $A_{V}$
---|---|---|---|---|---|---
| [$M_{\odot}~{}{\rm yr}^{-1}$] | [$10^{10}~{}M_{\odot}$] | [Gyr-1] | [$M_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$] | [$10^{8}~{}M_{\odot}~{}{\rm kpc}^{-2}$] |
Median$a$$a$footnotemark:  |  |  |  |  |  |
Total | $81_{-2}^{+7}$ | 4.9 $\pm$ 1.4 | 1.7 $\pm$ 0.4 | 1.3 $\pm$ 0.1 | 7.6 $\pm$ 2.2 | $1.6\pm 0.3$
Inner disk | $6.7_{-0.2}^{+0.2}$ | 1.3 $\pm$ 0.1 | 0.51 $\pm$ 0.04 | 0.88 $\pm$ 0.03 | 17 $\pm$ 1 | $2.04\pm 0.03$
Outer disk | $71_{-33}^{+16}$ | 4.9 $\pm$ 2.1 | 1.5 $\pm$ 0.9 | 0.7 $\pm$ 0.3 | 5.0 $\pm$ 2.1 | $1.4\pm 0.3$
Filter-by-filter$b$$b$footnotemark:  |  |  |  |  |  |
Total | $80_{-4}^{+9}$ | 5.2 $\pm$ 2.0 | 1.6 $\pm$ 0.5 | 1.2 $\pm$ 0.1 | 8.1 $\pm$ 3.1 | $1.8\pm 0.3$
Inner disk | $6.6_{-0.5}^{+0.5}$ | 1.3 $\pm$ 0.2 | 0.5 $\pm$ 0.1 | 0.87 $\pm$ 0.07 | 17 $\pm$ 3 | $2.0\pm 0.1$
Outer disk | $73_{-3}^{+10}$ | 4.7 $\pm$ 1.6 | 1.7 $\pm$ 0.5 | 0.74 $\pm$ 0.07 | 5 $\pm$ 2 | $1.6\pm 0.4$

aafootnotetext: Median magnification for Aperture 1 (see Fig. 3): $\mu_{\rm
med}=2.7\pm 0.2$.

bbfootnotetext: See magnifications in Table 5.

Note. — Inner vs. outer decomposition is performed at projected radius
$0.13\arcsec\approx 1.1$ kpc. The outer bound for the outer disk is taken to
be at 4.1 kpc, consistent with the effective radius of the dust continuum. As
the SED fit for the inner vs. outer disk is completely independent, summed
values may not agree with total values. All values are corrected for lensing
magnification, but using different methods—i.e. corrected with a single median
value, corrected by regime (UV, optical, or near-IR), or corrected filter-by-
filter. We don’t use the filter-based method for the decomposition due to the
large uncertainties.

Table 5: Magnification and bulge-to-total fraction ($B/T$) for image $1a$
only, used to correct the photometry in Table 3.

Filter | Magnification $\mu_{1a}$ | $B/T$
---|---|---
HST/F606W | 2.7 $\pm$ 0.4 | 0.04
HST/F625W | 2.7 $\pm$ 0.4 | 0.10
HST/F775W | 2.7 $\pm$ 0.4 | 0.10
HST/F814W | 3.0 $\pm$ 0.6 | —$a$$a$footnotemark:
HST/F850LP | 2.9 $\pm$ 0.5 | 0.19
JWST/F090W | 2.6 $\pm$ 0.4 | 0.11
HST/F105W | 3.0 $\pm$ 0.6 | 0.15
JWST/F115W | 2.6 $\pm$ 0.4 | 0.09
HST/F125W | 2.8 $\pm$ 0.5 | 0.07
HST/F140W | 3.1 $\pm$ 0.7 | 0.06
JWST/F150W | 2.8 $\pm$ 0.5 | 0.07
HST/F160W | 2.9 $\pm$ 0.5 | 0.10
JWST/F200W | 2.8 $\pm$ 0.5 | 0.09
JWST/F277W | 2.7 $\pm$ 0.4 | 0.13
JWST/F356W | 2.7 $\pm$ 0.4 | 0.20
JWST/F410M | 2.6 $\pm$ 0.4 | 0.19
JWST/F444W | 2.6 $\pm$ 0.4 | 0.22

aafootnotetext: F814W is excluded from the SED decomposition due to its
coincidence with a detector gap.

Note. — In using an aperture centered on this lensed counter-image,
differential magnification is minimized, while still effectively capturing
flux from the entire background object.

###  4.3 Five-image morphology from a bimodal mass distribution

One of the more remarkable aspects of El Anzuelo’s lensing morphology is the
complex structure apparent primarily in the longest wavelength filters of
NIRCam ($2.7-4.4~{}\mu$m), which reveal that the bright central, reddish core
of the background galaxy seen in the southwest image (1a) is multiply imaged
an additional 4 times towards the northeast ($1b$ through $1e$). This was
identified in the F444W image, which showed 5 distinct peaks. Using
predictions from our lens model, this appears to be due to a “beak-to-beak"
caustic metamorphosis, where the usual diamond-shaped caustics of each
foreground mass merges together into a single contiguous caustic network. The
position of El Anzuelo’s core at a short fold region between two closely-
separated cusps (see Figs. 2 and 3) closely mirrors Fig. 6 of Schneider &
Weiss (1986). There is a point-like source (most prominent in the F200W
filter) that lies barely on the interior of the southernmost caustic
(something akin to an intermediate location of what is shown in Figs. 6b and
6d of Schneider & Weiss). This alignment results in 5 multiple images, but it
is not obvious to what degree the core image is demagnified (in fact, it is
likely magnified). This analysis was extended from the binary point masses of
Schneider & Weiss (1986) to be applied to binary isothermal potentials by Shin
& Evans (2008), who found that the image configurations could include 3, 5, or
even 7 multiple images. Here, the 5-image morphology is closest to what Shin &
Evans label as configuration 5B-1 (their Fig. 11).

Based on our current lens model, we estimate the magnification near the core
image (labeled $1e$) to be $\mu\approx$1.5–8. This is seemingly confirmed by
the reasonably comparable flux observed in the demagnified image vs. the other
magnified images. However, there are likely some parts of the background DSFG
that are lensed to have a demagnified image, in particular near the core of
the secondary lensing galaxy. Pursuing this would benefit greatly from
spectroscopic confirmation with other tracers (e.g. molecular gas), or using
dust continuum at high-resolution. For now, we consider this interpretation of
a magnified 5th image to be the most likely, given the 5-image morphology that
we observe and the bimodal (nearly equal) mass distribution of the foreground.

However, there is still a chance that the core image is demagnified (although
the demagnification would need to be very weak). While central demagnified
images have long been predicted (Dyer & Roeder 1980; Burke 1981), they are
quite rarely observed, for three primary reasons: $a)$ they are typically
expected near the bright central cores of galaxies, where very high dynamic
range would be required to identify the small contrast of the light of the
lensed object (except at long wavelengths where foreground contamination is
low); $b)$ the images are demagnified, sometimes very strongly, so that they
often fall below the limiting magnitudes of observations; and $c)$ they
require a central surface mass density slope in the deflector that is
shallower than $\Sigma\propto r^{-1}$ (Rusin & Ma 2001). Few have been
discovered to date (e.g., Colley et al. 1996; Winn et al. 2004; Inada et al.
2005; Frye et al. 2007; Sharon et al. 2012; Collett et al. 2017; Ostrovski et
al. 2018; Muller et al. 2020), but they are uniquely valuable tools for
probing the inner density profiles of galaxies (Wallington & Narayan 1993;
Rusin & Ma 2001; Keeton 2003; Tamura et al. 2015; Wong et al. 2015). Whether
the 5th central image is demagnified or not, the discovery of this exotic
morphology is exciting and will provide observational tests of rare lensing
geometries, which in turn help constrain the mass density profiles of the
lensing (foreground) galaxies (e.g., Bozza & Melchiorre 2016; Bozza et al.
2020).

The lensing morphology also lends some clues to what we might expect for the
resolved distribution of dust continuum. The northeast component arises from
only the small inner portion of the disk interior to the lensing caustic
($\sim 0.1\arcsec$, or $\sim 800$pc, in radius), yet it is actually $\sim 2$
times brighter than the southwest component at 1mm (see Cheng et al. 2023).
This implies that a significant fraction of the total 1mm flux arises from
this high-magnification region. It is not possible to constrain this fraction
at present, however, as the northeast component consists of multiple images
blended together (see e.g. the 870$\mu$m structure in Fig. 5), each of varying
magnifications. Higher-resolution follow-up observations could test this and
directly map the extent of dust-obscured star formation in El Anzuelo, which
will have important implications for its current evolutionary stage and
relation to the main sequence.

Finally, we also remark that the best-fit caustic network is asymmetric, in
that there are two cusps in close proximity near the core of El Anzuelo, but
only a single cusp on the opposite end, towards the northwest. This is
possibly due to the external shear included in the model, or to the
orientation of the two SIE lenses relative to each other. In future work, it
may be worthwhile to employ the latest lensing models of the entire El Gordo
cluster (e.g. Diego et al. 2023; Frye et al. 2023) to make predictions for the
external shear that the cluster potential contributes at this location, and to
compare with the best-fit external shear amplitude of this work
($\gamma=0.36\pm 0.08$), which is quite substantial. We will refrain from
drawing definitive conclusions in this regard, as recent results from
Etherington et al. (2023) show that this added shear might only be to
compensate for insufficient model complexity. Yet, as discussed in §3.1.1, the
orientation of the best-fit shear is what would be expected for weak lensing
from the El Gordo cluster. We reserve a more in-depth analysis of these finer
details for a future work, as additional constraints will help to refine and
further test our lens model.

##  5 Conclusions

El Anzuelo is a $z$=2.291 ALMA-detected dusty star-forming galaxy that is
strongly lensed by the massive $z=0.87$ galaxy cluster, El Gordo. New
JWST/NIRCam imaging reveals an extraordinary lensing morphology, with some
portions of the background source being multiply-imaged five times. After
developing a new parametric lens model specifically tailored to this object
using lenstool, we take advantage of the magnification and amplification to
examine the surface brightness distribution of the rest-frame UV through near-
IR ($0.18-1.3~{}\mu$m) in this pilot study. The near-IR is expected to be a
better proxy for stellar mass, as it is more sensitive to light from redder,
lower-mass stars that account for most of a galaxy’s stellar mass. This is in
contrast to rest-frame UV or blue optical imaging, which is both more affected
by dust obscuration and more easily biased by the light from short-lived
massive (and more luminous) stars.

We find the peculiar result that the source-plane surface brightness
distribution is largely consistent between 1.2 and $4.4~{}\mu$m (or rest-frame
$0.35-1.3~{}\mu$m).We interpret the ostensible match of the near-UV light and
the stellar distribution as the combined effect of a central concentration of
dust and widespread recent star formation. Dust attenuation may truncate the
peak of UV light that would otherwise be seen at relatively smaller radii than
longer wavelengths, which should have the effect of broadening the UV sizes.
Since this is not the case, we interpret that the UV and near-IR are tracing
different components (suggested also by the $\sim$2.5 kpc offset in their
peaks), perhaps from an ongoing late-stage coalescence of a galaxy merger.

In a decomposition of the SED for an inner, bulge-like component for the
central kiloparsec versus the outskirts of the disk, we find elevated specific
star formation rates at larger radii: $1.5\pm 0.9~{}{\rm Gyr}^{-1}$ in the
outer disk and $0.51\pm 0.04~{}{\rm Gyr}^{-1}$ for the inner disk. We
interpret the suppressed sSFR in the galaxy center to be a possible indication
of the early stages of inside-out galaxy quenching, possibly driven by the
inward radial migration of star-forming clumps distributed throughout the
disk. The overall (magnification-corrected) intrinsic SFR is estimated to be
$81^{+7}_{-2}~{}M_{\odot}~{}{\rm yr}^{-1}$, consistent with the star-forming
main sequence at $z=2.3$ (and suggesting an inferred intrinsic IR luminosity
of $L_{IR}\approx 5\times 10^{11}~{}L_{\odot}$). This may be an underestimate,
however, as it is derived only from the UV, optical, and near-IR SED. Higher-
resolution far-IR measurements are needed to test if the DSFG might lie above
the main sequence, as the peak of archival low-resolution imaging coincides
with the region of greatest magnification, such that the intrinsic (de-
magnified) far-IR SED is poorly constrained.

We note the large differential magnification indicated by our lens model, by
which the near-IR ($\mu_{\rm tot}\sim$ 8–12) is magnified approximately four
times more strongly than the UV ($\mu_{\rm tot}\sim$ 2–3). This is directly
related to the surface brightness profiles we observe, as dust in the galaxy
center (which is coincidentally aligned with the high-magnification region of
the source plane) preferentially attenuates shorter wavelengths. Finally, we
briefly discuss the implications of a 5-image lensing morphology arising from
a bimodal foreground mass distribution. These include the possibility of
detecting central demagnified images, which are a powerful probe of the
central mass densities of lensing galaxies and groups, but which have so far
proven quite elusive to discover. The serendipitous alignment of the galaxy
nucleus with a part of the source-plane that is both highly-magnified and
multiply-imaged 5 times makes El Anzuelo a very attractive target for future
study (e.g., in examining possible variability in an AGN).

We thank the anonymous referee for their suggestions to improve this
manuscript. PSK would like to thank Ian Smail, Allison Noble, and Alex
Pigarelli for helpful discussion in framing this work. 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 program # 1176 (PEARLS;
PI: R. Windhorst). RAW, SHC, and RAJ acknowledge support from NASA JWST
Interdisciplinary Scientist grants NAG5-12460, NNX14AN10G and 80NSSC18K0200
from GSFC. This research is also based on observations made with the NASA/ESA
Hubble Space Telescope obtained from the Space Telescope Science Institute,
which is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS 5—26555. These observations are
associated with programs GO-12477, GO-12755, and GO-14096. The specific
observations from JWST and HST are available at:
https://doi.org/10.17909/7cnj-5y94 (catalog doi:10.17909/7cnj-5y94). This
paper makes use of the following ALMA data: ADS/JAO.ALMA# 2013.1.01051.S,
2013.1.01358.S, 2015.1.01187.S, 2017.1.01621.S, and 2018.1.00035.L. ALMA is a
partnership of ESO (representing its member states), NSF (USA) and NINS
(Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio
Astronomy Observatory is a facility of the National Science Foundation
operated under cooperative agreement by Associated Universities, Inc. This
research has made extensive use of NASA’s Astrophysics Data System.

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.

## References

  * Alton et al. (1999) Alton, P. B., Davies, J. I., & Bianchi, S. 1999, A&A, 343, 51
  * Armus et al. (2007) Armus, L., Charmandaris, V., Bernard-Salas, J., et al. 2007, ApJ, 656, 148, doi: 10.1086/510107
  * Armus et al. (2009) Armus, L., Mazzarella, J. M., Evans, A. S., et al. 2009, PASP, 121, 559, doi: 10.1086/600092
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Barro et al. (2013) Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2013, ApJ, 765, 104, doi: 10.1088/0004-637X/765/2/104
  * Barro et al. (2016) Barro, G., Kriek, M., Pérez-González, P. G., et al. 2016, ApJ, 827, L32, doi: 10.3847/2041-8205/827/2/L32
  * Barrufet et al. (2023) Barrufet, L., Oesch, P. A., Weibel, A., et al. 2023, MNRAS, 522, 449, doi: 10.1093/mnras/stad947
  * Basu et al. (2016) Basu, K., Sommer, M., Erler, J., et al. 2016, ApJ, 829, L23, doi: 10.3847/2041-8205/829/2/L23
  * Battisti et al. (2019) Battisti, A. J., da Cunha, E., Grasha, K., et al. 2019, ApJ, 882, 61, doi: 10.3847/1538-4357/ab345d
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Blain et al. (2002) Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002, Phys. Rep., 369, 111, doi: 10.1016/S0370-1573(02)00134-5
  * Bondi et al. (2018) Bondi, M., Zamorani, G., Ciliegi, P., et al. 2018, A&A, 618, L8, doi: 10.1051/0004-6361/201834243
  * Bournaud et al. (2011) Bournaud, F., Chapon, D., Teyssier, R., et al. 2011, ApJ, 730, 4, doi: 10.1088/0004-637X/730/1/4
  * Bozza & Melchiorre (2016) Bozza, V., & Melchiorre, C. 2016, J. Cosmology Astropart. Phys, 2016, 040, doi: 10.1088/1475-7516/2016/03/040
  * Bozza et al. (2020) Bozza, V., Pietroni, S., & Melchiorre, C. 2020, Universe, 6, 106, doi: 10.3390/universe6080106
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.6.0, 1.6.0, Zenodo, Zenodo, doi: 10.5281/zenodo.7419741
  * Brammer (2021) Brammer, G. 2021, gbrammer/eazy-py: Tagged release 2021, 0.5.2, Zenodo, Zenodo, doi: 10.5281/zenodo.5012705
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Burke (1981) Burke, W. L. 1981, ApJ, 244, L1, doi: 10.1086/183466
  * Bussmann et al. (2012) Bussmann, R. S., Gurwell, M. A., Fu, H., et al. 2012, ApJ, 756, 134, doi: 10.1088/0004-637X/756/2/134
  * Bussmann et al. (2013) Bussmann, R. S., Pérez-Fournon, I., Amber, S., et al. 2013, ApJ, 779, 25, doi: 10.1088/0004-637X/779/1/25
  * Calanog et al. (2014) Calanog, J. A., Fu, H., Cooray, A., et al. 2014, ApJ, 797, 138, doi: 10.1088/0004-637X/797/2/138
  * Calistro Rivera et al. (2018) Calistro Rivera, G., Hodge, J. A., Smail, I., et al. 2018, ApJ, 863, 56, doi: 10.3847/1538-4357/aacffa
  * Caminha et al. (2022) Caminha, G. B., Grillo, C., Rosati, P., et al. 2022, arXiv e-prints, arXiv:2209.02718.  https://arxiv.org/abs/2209.02718
  * Caputi et al. (2021) Caputi, K. I., Caminha, G. B., Fujimoto, S., et al. 2021, ApJ, 908, 146, doi: 10.3847/1538-4357/abd4d0
  * Carrasco et al. (2010) Carrasco, E. R., Conselice, C. J., & Trujillo, I. 2010, MNRAS, 405, 2253, doi: 10.1111/j.1365-2966.2010.16645.x
  * Casey et al. (2014a) Casey, C. M., Narayanan, D., & Cooray, A. 2014a, Phys. Rep., 541, 45, doi: 10.1016/j.physrep.2014.02.009
  * Casey et al. (2014b) Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014b, ApJ, 796, 95, doi: 10.1088/0004-637X/796/2/95
  * Ceverino et al. (2010) Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151, doi: 10.1111/j.1365-2966.2010.16433.x
  * Ceverino et al. (2012) Ceverino, D., Dekel, A., Mandelker, N., et al. 2012, MNRAS, 420, 3490, doi: 10.1111/j.1365-2966.2011.20296.x
  * Chapman et al. (2004) Chapman, S. C., Smail, I., Windhorst, R., Muxlow, T., & Ivison, R. J. 2004, ApJ, 611, 732, doi: 10.1086/422383
  * Chen et al. (2015) Chen, C.-C., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 799, 194, doi: 10.1088/0004-637X/799/2/194
  * Chen et al. (2017) Chen, C.-C., Hodge, J. A., Smail, I., et al. 2017, ApJ, 846, 108, doi: 10.3847/1538-4357/aa863a
  * Cheng et al. (2023) Cheng, C., Huang, J.-S., Smail, I., et al. 2023, ApJ, 942, L19, doi: 10.3847/2041-8213/aca9d0
  * Cheung et al. (2012) Cheung, E., Faber, S. M., Koo, D. C., et al. 2012, ApJ, 760, 131, doi: 10.1088/0004-637X/760/2/131
  * Cochrane et al. (2019) Cochrane, R. K., Hayward, C. C., Anglés-Alcázar, D., et al. 2019, MNRAS, 488, 1779, doi: 10.1093/mnras/stz1736
  * Cochrane et al. (2021) Cochrane, R. K., Best, P. N., Smail, I., et al. 2021, MNRAS, 503, 2622, doi: 10.1093/mnras/stab467
  * Coe et al. (2019) Coe, D., Salmon, B., Bradač, M., et al. 2019, ApJ, 884, 85, doi: 10.3847/1538-4357/ab412b
  * Collett et al. (2017) Collett, T. E., Buckley-Geer, E., Lin, H., et al. 2017, ApJ, 843, 148, doi: 10.3847/1538-4357/aa76e6
  * Colley et al. (1996) Colley, W. N., Tyson, J. A., & Turner, E. L. 1996, ApJ, 461, L83, doi: 10.1086/310015
  * Cresci et al. (2010) Cresci, G., Mannucci, F., Maiolino, R., et al. 2010, Nature, 467, 811, doi: 10.1038/nature09451
  * Crockett et al. (2012) Crockett, R. M., Shabala, S. S., Kaviraj, S., et al. 2012, MNRAS, 421, 1603, doi: 10.1111/j.1365-2966.2012.20418.x
  * Cunow (2001) Cunow, B. 2001, MNRAS, 323, 130, doi: 10.1046/j.1365-8711.2001.04140.x
  * Cutler et al. (2023) Cutler, S. E., Giavalisco, M., Ji, Z., & Cheng, Y. 2023, ApJ, 945, 97, doi: 10.3847/1538-4357/acb5e9
  * de Jong (1996) de Jong, R. S. 1996, A&A, 313, 377, doi: 10.48550/arXiv.astro-ph/9604010
  * de Vaucouleurs (1958) de Vaucouleurs, G. 1958, ApJ, 127, 487, doi: 10.1086/146476
  * Dekel & Burkert (2014) Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870, doi: 10.1093/mnras/stt2331
  * Dekel et al. (2009) Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451, doi: 10.1038/nature07648
  * Desai et al. (2007) Desai, V., Armus, L., Spoon, H. W. W., et al. 2007, ApJ, 669, 810, doi: 10.1086/522104
  * Devour & Bell (2017) Devour, B. M., & Bell, E. F. 2017, MNRAS, 468, L31, doi: 10.1093/mnrasl/slx021
  * Devour & Bell (2019) —. 2019, ApJS, 244, 3, doi: 10.3847/1538-4365/ab339c
  * Dey et al. (2008) Dey, A., Soifer, B. T., Desai, V., et al. 2008, ApJ, 677, 943, doi: 10.1086/529516
  * Diego et al. (2023) Diego, J. M., Meena, A. K., Adams, N. J., et al. 2023, A&A, 672, A3, doi: 10.1051/0004-6361/202245238
  * Dimauro et al. (2022) Dimauro, P., Daddi, E., Shankar, F., et al. 2022, MNRAS, 513, 256, doi: 10.1093/mnras/stac884
  * Dressler et al. (1987) Dressler, A., Lynden-Bell, D., Burstein, D., et al. 1987, ApJ, 313, 42, doi: 10.1086/164947
  * Driver et al. (2007a) Driver, S. P., Allen, P. D., Liske, J., & Graham, A. W. 2007a, ApJ, 657, L85, doi: 10.1086/513106
  * Driver et al. (2007b) Driver, S. P., Popescu, C. C., Tuffs, R. J., et al. 2007b, MNRAS, 379, 1022, doi: 10.1111/j.1365-2966.2007.11862.x
  * Dyer & Roeder (1980) Dyer, C. C., & Roeder, R. C. 1980, ApJ, 238, L67, doi: 10.1086/183260
  * Eisenhardt et al. (2012) Eisenhardt, P. R. M., Wu, J., Tsai, C.-W., et al. 2012, ApJ, 755, 173, doi: 10.1088/0004-637X/755/2/173
  * Elbaz et al. (2018) Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110, doi: 10.1051/0004-6361/201732370
  * Ellison et al. (2018) Ellison, S. L., Sánchez, S. F., Ibarra-Medel, H., et al. 2018, MNRAS, 474, 2039, doi: 10.1093/mnras/stx2882
  * Ellison et al. (2020) Ellison, S. L., Thorp, M. D., Pan, H.-A., et al. 2020, MNRAS, 492, 6027, doi: 10.1093/mnras/staa001
  * Etherington et al. (2023) Etherington, A., Nightingale, J. W., Massey, R., et al. 2023, arXiv e-prints, arXiv:2301.05244, doi: 10.48550/arXiv.2301.05244
  * Evans et al. (2022) Evans, A. S., Frayer, D. T., Charmandaris, V., et al. 2022, ApJ, 940, L8, doi: 10.3847/2041-8213/ac9971
  * Evans (1994) Evans, R. 1994, MNRAS, 266, 511, doi: 10.1093/mnras/266.2.511
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455, doi: 10.1146/annurev-astro-081811-125521
  * Ferreira et al. (2022) Ferreira, L., Adams, N., Conselice, C. J., et al. 2022, ApJ, 938, L2, doi: 10.3847/2041-8213/ac947c
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13, doi: 10.3847/2041-8213/acade4
  * Franco et al. (2018) Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152, doi: 10.1051/0004-6361/201832928
  * Frye et al. (2007) Frye, B. L., Coe, D., Bowen, D. V., et al. 2007, ApJ, 665, 921, doi: 10.1086/519244
  * Frye et al. (2023) Frye, B. L., Pascale, M., Foo, N., et al. 2023, ApJ, 952, 81, doi: 10.3847/1538-4357/acd929
  * Gadotti & dos Anjos (2001) Gadotti, D. A., & dos Anjos, S. 2001, AJ, 122, 1298, doi: 10.1086/322126
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2, doi: 10.1051/0004-6361/201629512
  * Gaia Collaboration et al. (2021) —. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gao & Ho (2017) Gao, H., & Ho, L. C. 2017, ApJ, 845, 114, doi: 10.3847/1538-4357/aa7da4
  * Goldader et al. (2002) Goldader, J. D., Meurer, G., Heckman, T. M., et al. 2002, ApJ, 568, 651, doi: 10.1086/339165
  * Goulding et al. (2010) Goulding, A. D., Alexander, D. M., Lehmer, B. D., & Mullaney, J. R. 2010, MNRAS, 406, 597, doi: 10.1111/j.1365-2966.2010.16700.x
  * Graham & Worley (2008) Graham, A. W., & Worley, C. C. 2008, MNRAS, 388, 1708, doi: 10.1111/j.1365-2966.2008.13506.x
  * Gullberg et al. (2019) Gullberg, B., Smail, I., Swinbank, A. M., et al. 2019, MNRAS, 490, 4956, doi: 10.1093/mnras/stz2835
  * Hales et al. (2012) Hales, C. A., Murphy, T., Curran, J. R., et al. 2012, MNRAS, 425, 979, doi: 10.1111/j.1365-2966.2012.21373.x
  * Häußler et al. (2013) Häußler, B., Bamford, S. P., Vika, M., et al. 2013, MNRAS, 430, 330, doi: 10.1093/mnras/sts633
  * Hayward et al. (2011) Hayward, C. C., Kereš, D., Jonsson, P., et al. 2011, ApJ, 743, 159, doi: 10.1088/0004-637X/743/2/159
  * Heckman (2001) Heckman, T. M. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 240, Gas and Galaxy Evolution, ed. J. E. Hibbard, M. Rupen, & J. H. van Gorkom, 345, doi: 10.48550/arXiv.astro-ph/0009075
  * Hezaveh et al. (2012) Hezaveh, Y. D., Marrone, D. P., & Holder, G. P. 2012, ApJ, 761, 20, doi: 10.1088/0004-637X/761/1/20
  * Hodge et al. (2015) Hodge, J. A., Riechers, D., Decarli, R., et al. 2015, ApJ, 798, L18, doi: 10.1088/2041-8205/798/1/L18
  * Hodge et al. (2016) Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103, doi: 10.3847/1538-4357/833/1/103
  * Hodge et al. (2019) Hodge, J. A., Smail, I., Walter, F., et al. 2019, ApJ, 876, 130, doi: 10.3847/1538-4357/ab1846
  * Hopkins et al. (2010) Hopkins, P. F., Bundy, K., Hernquist, L., Wuyts, S., & Cox, T. J. 2010, MNRAS, 401, 1099, doi: 10.1111/j.1365-2966.2009.15699.x
  * Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581, doi: 10.1093/mnras/stu1738
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Huang et al. (2011) Huang, J. S., Zheng, X. Z., Rigopoulou, D., et al. 2011, ApJ, 742, L13, doi: 10.1088/2041-8205/742/1/L13
  * Hubble (1943) Hubble, E. 1943, ApJ, 97, 112, doi: 10.1086/144504
  * Hutton et al. (2014) Hutton, S., Ferreras, I., Wu, K., et al. 2014, MNRAS, 440, 150, doi: 10.1093/mnras/stu185
  * Ikarashi et al. (2015) Ikarashi, S., Ivison, R. J., Caputi, K. I., et al. 2015, ApJ, 810, 133, doi: 10.1088/0004-637X/810/2/133
  * Inada et al. (2005) Inada, N., Oguri, M., Keeton, C. R., et al. 2005, PASJ, 57, L7, doi: 10.1093/pasj/57.3.L7
  * Iono et al. (2004) Iono, D., Ho, P. T. P., Yun, M. S., et al. 2004, ApJ, 616, L63, doi: 10.1086/420784
  * Jansen et al. (2000) Jansen, R. A., Franx, M., Fabricant, D., & Caldwell, N. 2000, ApJS, 126, 271, doi: 10.1086/313303
  * Jansen et al. (1994) Jansen, R. A., Knapen, J. H., Beckman, J. E., Peletier, R. F., & Hes, R. 1994, MNRAS, 270, 373, doi: 10.1093/mnras/270.2.373
  * Ji & Giavalisco (2022) Ji, Z., & Giavalisco, M. 2022, ApJ, 935, 120, doi: 10.3847/1538-4357/ac7f43
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22, doi: 10.3847/1538-4365/abef67
  * Jones et al. (2010) Jones, T., Ellis, R., Jullo, E., & Richard, J. 2010, ApJ, 725, L176, doi: 10.1088/2041-8205/725/2/L176
  * Jullo & Kneib (2009) Jullo, E., & Kneib, J.-P. 2009, MNRAS, 395, 1319, doi: 10.1111/j.1365-2966.2009.14654.x
  * Jullo et al. (2007) Jullo, E., Kneib, J.-P., Limousin, M., et al. 2007, New Journal of Physics, 9, 447, doi: 10.1088/1367-2630/9/12/447
  * Kamieneski et al. (2023) Kamieneski, P. S., Yun, M. S., Harrington, K. C., et al. 2023, arXiv e-prints, arXiv:2301.09746.  https://arxiv.org/abs/2301.09746
  * Keeton (2003) Keeton, C. R. 2003, ApJ, 582, 17, doi: 10.1086/344539
  * Kelvin et al. (2012) Kelvin, L. S., Driver, S. P., Robotham, A. S. G., et al. 2012, MNRAS, 421, 1007, doi: 10.1111/j.1365-2966.2012.20355.x
  * Kennedy et al. (2015) Kennedy, R., Bamford, S. P., Baldry, I., et al. 2015, MNRAS, 454, 806, doi: 10.1093/mnras/stv2032
  * Kennedy et al. (2016) Kennedy, R., Bamford, S. P., Häußler, B., et al. 2016, MNRAS, 460, 3458, doi: 10.1093/mnras/stw1176
  * Kneib et al. (1996) Kneib, J.-P., Ellis, R. S., Smail, I., Couch, W. J., & Sharples, R. M. 1996, ApJ, 471, 643, doi: 10.1086/177995
  * Kneib et al. (1993) Kneib, J. P., Mellier, Y., Fort, B., & Mathez, G. 1993, A&A, 273, 367
  * Knop et al. (1994) Knop, R. A., Soifer, B. T., Graham, J. R., et al. 1994, AJ, 107, 920, doi: 10.1086/116906
  * Kormann et al. (1994) Kormann, R., Schneider, P., & Bartelmann, M. 1994, A&A, 284, 285
  * Kron (1980) Kron, R. G. 1980, ApJS, 43, 305, doi: 10.1086/190669
  * La Barbera et al. (2002) La Barbera, F., Busarello, G., Merluzzi, P., Massarotti, M., & Capaccioli, M. 2002, ApJ, 571, 790, doi: 10.1086/340021
  * La Barbera et al. (2010) La Barbera, F., de Carvalho, R. R., de La Rosa, I. G., et al. 2010, MNRAS, 408, 1313, doi: 10.1111/j.1365-2966.2010.16850.x
  * La Barbera et al. (2012) La Barbera, F., Ferreras, I., de Carvalho, R. R., et al. 2012, MNRAS, 426, 2300, doi: 10.1111/j.1365-2966.2012.21848.x
  * Labbé et al. (2005) Labbé, I., Huang, J., Franx, M., et al. 2005, ApJ, 624, L81, doi: 10.1086/430700
  * Lang et al. (2019) Lang, P., Schinnerer, E., Smail, I., et al. 2019, ApJ, 879, 54, doi: 10.3847/1538-4357/ab1f77
  * Leja et al. (2019) Leja, J., Tacchella, S., & Conroy, C. 2019, ApJ, 880, L9, doi: 10.3847/2041-8213/ab2f8c
  * Li & Draine (2001) Li, A., & Draine, B. T. 2001, ApJ, 554, 778, doi: 10.1086/323147
  * Liu et al. (2017) Liu, F. S., Jiang, D., Faber, S. M., et al. 2017, ApJ, 844, L2, doi: 10.3847/2041-8213/aa7cf5
  * Liu et al. (2013) Liu, G., Calzetti, D., Hong, S., et al. 2013, ApJ, 778, L41, doi: 10.1088/2041-8205/778/2/L41
  * Lonsdale et al. (2006) Lonsdale, C. J., Farrah, D., & Smith, H. E. 2006, Ultraluminous Infrared Galaxies (Springer-Verlag), 285, doi: 10.1007/3-540-30313-8_9
  * Luo et al. (2022) Luo, Y., Rowlands, K., Alatalo, K., et al. 2022, ApJ, 938, 63, doi: 10.3847/1538-4357/ac8b7d
  * Lutz et al. (2016) Lutz, D., Berta, S., Contursi, A., et al. 2016, A&A, 591, A136, doi: 10.1051/0004-6361/201527706
  * Ma et al. (2018) Ma, X., Hopkins, P. F., Boylan-Kolchin, M., et al. 2018, MNRAS, 477, 219, doi: 10.1093/mnras/sty684
  * Ma et al. (2019) Ma, X., Hayward, C. C., Casey, C. M., et al. 2019, MNRAS, 487, 1844, doi: 10.1093/mnras/stz1324
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Marian et al. (2018) Marian, V., Ziegler, B., Kuchner, U., & Verdugo, M. 2018, A&A, 617, A34, doi: 10.1051/0004-6361/201832750
  * Martig et al. (2009) Martig, M., Bournaud, F., Teyssier, R., & Dekel, A. 2009, ApJ, 707, 250, doi: 10.1088/0004-637X/707/1/250
  * Massardi et al. (2018) Massardi, M., Enia, A. F. M., Negrello, M., et al. 2018, A&A, 610, A53, doi: 10.1051/0004-6361/201731751
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Menanteau et al. (2012) Menanteau, F., Hughes, J. P., Sifón, C., et al. 2012, ApJ, 748, 7, doi: 10.1088/0004-637X/748/1/7
  * Miller et al. (2022) Miller, T. B., Whitaker, K. E., Nelson, E. J., et al. 2022, ApJ, 941, L37, doi: 10.3847/2041-8213/aca675
  * Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319, doi: 10.1046/j.1365-8711.1998.01227.x
  * Möllenhoff et al. (2006) Möllenhoff, C., Popescu, C. C., & Tuffs, R. J. 2006, A&A, 456, 941, doi: 10.1051/0004-6361:20054727
  * Mowla et al. (2019) Mowla, L. A., van Dokkum, P., Brammer, G. B., et al. 2019, ApJ, 880, 57, doi: 10.3847/1538-4357/ab290a
  * Muller et al. (2020) Muller, S., Jaswanth, S., Horellou, C., & Martí-Vidal, I. 2020, A&A, 641, L2, doi: 10.1051/0004-6361/202038978
  * Narayan & Bartelmann (1996) Narayan, R., & Bartelmann, M. 1996, ArXiv Astrophysics e-prints
  * Narayanan et al. (2018) Narayanan, D., Davé, R., Johnson, B. D., et al. 2018, MNRAS, 474, 1718, doi: 10.1093/mnras/stx2860
  * Nardiello et al. (2022) Nardiello, D., Bedin, L. R., Burgasser, A., et al. 2022, MNRAS, 517, 484, doi: 10.1093/mnras/stac2659
  * Nelson et al. (2012) Nelson, E. J., van Dokkum, P. G., Brammer, G., et al. 2012, ApJ, 747, L28, doi: 10.1088/2041-8205/747/2/L28
  * Nelson et al. (2016a) Nelson, E. J., van Dokkum, P. G., Momcheva, I. G., et al. 2016a, ApJ, 817, L9, doi: 10.3847/2041-8205/817/1/L9
  * Nelson et al. (2016b) Nelson, E. J., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2016b, ApJ, 828, 27, doi: 10.3847/0004-637X/828/1/27
  * Oser et al. (2010) Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312, doi: 10.1088/0004-637X/725/2/2312
  * Ostrovski et al. (2018) Ostrovski, F., Lemon, C. A., Auger, M. W., et al. 2018, MNRAS, 473, L116, doi: 10.1093/mnrasl/slx173
  * Pantoni et al. (2021) Pantoni, L., Massardi, M., Lapi, A., et al. 2021, MNRAS, 507, 3998, doi: 10.1093/mnras/stab2346
  * Pascale et al. (2022) Pascale, M., Frye, B. L., Dai, L., et al. 2022, ApJ, 932, 85, doi: 10.3847/1538-4357/ac6ce9
  * Pastrav et al. (2013) Pastrav, B. A., Popescu, C. C., Tuffs, R. J., & Sansom, A. E. 2013, A&A, 553, A80, doi: 10.1051/0004-6361/201220962
  * Patel et al. (2012) Patel, S. G., Holden, B. P., Kelson, D. D., et al. 2012, ApJ, 748, L27, doi: 10.1088/2041-8205/748/2/L27
  * Peletier & Balcells (1996) Peletier, R. F., & Balcells, M. 1996, AJ, 111, 2238, doi: 10.1086/117958
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266, doi: 10.1086/340952
  * Penney et al. (2020) Penney, J. I., Blain, A. W., Assef, R. J., et al. 2020, MNRAS, 496, 1565, doi: 10.1093/mnras/staa1582
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 91433X, doi: 10.1117/12.2056689
  * Perrin et al. (2012) Perrin, M. D., Soummer, R., Elliott, E. M., Lallo, M. D., & Sivaramakrishnan, A. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, ed. M. C. Clampin, G. G. Fazio, H. A. MacEwen, & J. Oschmann, Jacobus M., 84423D, doi: 10.1117/12.925230
  * Petric et al. (2011) Petric, A. O., Armus, L., Howell, J., et al. 2011, ApJ, 730, 28, doi: 10.1088/0004-637X/730/1/28
  * Pilkington et al. (2012) Pilkington, K., Few, C. G., Gibson, B. K., et al. 2012, A&A, 540, A56, doi: 10.1051/0004-6361/201117466
  * Pope et al. (2008) Pope, A., Bussmann, R. S., Dey, A., et al. 2008, ApJ, 689, 127, doi: 10.1086/592739
  * Popping et al. (2017) Popping, G., Puglisi, A., & Norman, C. A. 2017, MNRAS, 472, 2315, doi: 10.1093/mnras/stx2202
  * Popping et al. (2022) Popping, G., Pillepich, A., Calistro Rivera, G., et al. 2022, MNRAS, 510, 3321, doi: 10.1093/mnras/stab3312
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001, doi: 10.1088/1538-3873/acb293
  * Robertson et al. (2020) Robertson, A., Smith, G. P., Massey, R., et al. 2020, MNRAS, 495, 3727, doi: 10.1093/mnras/staa1429