#  PASSAGES: the wide-ranging, extreme intrinsic properties of
$\boldsymbol{\it Planck}$-selected, lensed dusty star-forming galaxies

Patrick S. Kamieneski Department of Astronomy, University of Massachusetts,
Amherst, MA 01003, USA  School of Earth and Space Exploration, Arizona State
University, Tempe, AZ 85287-6004, USA Patrick S. Kamieneski<EMAIL_ADDRESS>Min S. Yun Department of Astronomy, University of Massachusetts, Amherst, MA
01003, USA  Kevin C. Harrington European Southern Observatory, Alonso de
Córdova 3107, Viticura, Casilla 19001, Santiago de Chile, Chile  James D.
Lowenthal Smith College, Northampton, MA 01063, USA  Q. Daniel Wang Department
of Astronomy, University of Massachusetts, Amherst, MA 01003, USA  Brenda L.
Frye Department of Astronomy/Steward Observatory, 933 North Cherry Avenue,
University of Arizona, Tucson, AZ 85721, USA  Eric F. Jiménez-Andrade
Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de
México, Antigua Carretera a Pátzcuaro # 8701, Ex-Hda. San José de la Huerta,
Morelia, Michoacán, México C.P. 58089  Amit Vishwas Department of Astronomy,
Cornell University, Space Sciences Building, Ithaca, NY 14853, USA  Olivia
Cooper Department of Astronomy, The University of Texas at Austin, 2515
Speedway Boulevard Stop C1400, Austin, TX 78712, USA  Massimo Pascale
Department of Astronomy, University of California, 501 Campbell Hall #3411,
Berkeley, CA 94720, USA  Nicholas Foo Department of Astronomy/Steward
Observatory, 933 North Cherry Avenue, University of Arizona, Tucson, AZ 85721,
USA  Derek Berman Department of Earth and Atmospheric Sciences, Cornell
University, Ithaca, NY 14853, USA  Anthony Englert Department of Physics,
Brown University, Providence, RI 02912, USA  Carlos Garcia Diaz Department of
Astronomy, University of Massachusetts, Amherst, MA 01003, USA

###### Abstract

The PASSAGES (Planck All-Sky Survey to Analyze Gravitationally-lensed Extreme
Starbursts) collaboration has recently defined a sample of 30 gravitationally-
lensed dusty star-forming galaxies (DSFGs). These rare, submillimeter-selected
objects enable high-resolution views of the most extreme sites of star
formation in galaxies at Cosmic Noon. Here, we present the first major
compilation of strong lensing analyses using lenstool for PASSAGES, including
15 objects spanning $z=1.1-3.3$, using complementary information from
$0.6^{\prime\prime}$-resolution 1 mm Atacama Large Millimeter/submillimeter
Array (ALMA) and $0.4^{\prime\prime}$ 5 cm Jansky Very Large Array continuum
imaging, in tandem with 1.6$\micron$ Hubble and optical imaging with Gemini-S.
Magnifications range from $\mu={2-28}$ (median $\mu=7$), yielding intrinsic
infrared luminosities of $L_{\rm IR}={0.2-5.9}\times 10^{13}~{}L_{\odot}$
(median ${1.4}\times 10^{13}~{}L_{\odot}$) and inferred star formation rates
of ${170-6300}~{}M_{\odot}~{}{\rm yr}^{-1}$ (median ${1500}~{}M_{\odot}~{}{\rm
yr}^{-1}$). These results suggest that the PASSAGES objects comprise some of
the most extreme known starbursts, rivaling the luminosities of even the
brightest unlensed objects, further amplified by lensing. The intrinsic sizes
of far-infrared continuum regions are large ($R_{\rm e}={1.7-4.3}$ kpc; median
$3.0$ kpc) but consistent with $L_{\rm IR}-R_{\rm e}$ scaling relations for
$z>1$ DSFGs, suggesting a widespread spatial distribution of star formation.
With modestly-high angular resolution, we explore if these objects might be
maximal starbursts. Instead of approaching Eddington-limited surface
densities, above which radiation pressure will disrupt further star formation,
they are safely sub-Eddington—at least on global, galaxy-integrated scales.

Gravitational lensing (670) — Ultraluminous infrared galaxies (1735) —
Starburst galaxies (1570)

††journal: ApJ††facilities: HST(WFC3), LMT, ALMA, JVLA, Gemini††software:
astropy (Astropy Collaboration et al. 2013, 2018), lenstool(Kneib et al. 1993,
1996; Jullo et al. 2007; Jullo & Kneib 2009), APLpy (Robitaille & Bressert
2012), IRAF (Tody 1986, 1993; Science Software Branch at STScI 2012), CASA
(McMullin et al. 2007), Ned Wright’s Cosmology Calculator (Wright 2006)

PASSAGES: Lens modeling of Planck DSFGsPASSAGES: Lens modeling of Planck DSFGs

##  1 Introduction

Dusty star-forming galaxies (DSFGs) host some of the most extreme infrared
luminosities in the known Universe, often in excess of
$10^{12-13}~{}L_{\odot}$, reflecting star formation rates $>100-1000$ $\rm
M_{\odot}~{}yr^{-1}$. Observed mainly during the peak of the cosmic star
formation history at $z=1-4$, they are fundamental to our understanding of
stellar mass assembly in the last 10 billion years (see reviews by Blain et
al. 2002; Casey et al. 2014). The current galaxy evolution paradigm favors
that DSFGs are progenitors of massive, quiescent elliptical galaxies (e.g.,
Lilly et al. 1999; Brodwin et al. 2008; Tacconi et al. 2008; Daddi et al.
2009; Toft et al. 2014), in part because of their association with galaxy
over-densities and mergers. However, numerous questions remain regarding the
rapid conversion of cold gas reservoirs into new stars, often at star
formation surface densities approaching theoretical maxima. Beyond this
Eddington-like limit, stellar feedback should hinder or even quench subsequent
star formation (Scoville et al. 2001; Scoville 2003; Murray et al. 2005;
Thompson et al. 2005; Andrews & Thompson 2011; Hodge et al. 2019). Galaxies
that are sub-Eddington when integrated over star-forming regions may still
have Eddington-limited clumps of star formation (e.g., Simpson et al. 2015a;
Barcos-Muñoz et al. 2017). Moreover, the predominant triggering mechanism for
fueling this star formation has lingering uncertainties. There is now
plentiful evidence (observational and theoretical) that major galaxy mergers
play an important role (Ivison et al. 2002; Narayanan et al. 2006, 2010;
Tacconi et al. 2008, among many others), but cold-mode accretion, in which
star formation in massive galaxies is fed by smooth infall of gas, has also
been put forward as a viable explanation (Kereš et al. 2005; Dekel et al.
2009b, a). Answering these open questions absolutely necessitates resolutions
corresponding to sub-galactic physical scales.

The past decade has ushered in the discovery of a substantial number of DSFGs
that are amplified in submillimeter flux by strong gravitational lensing. Due
to the steep drop-off in submillimeter number counts for high-$z$ galaxies
(Blain 1996; Negrello et al. 2007), lensed DSFGs can be efficiently identified
using a simple flux threshold. This was first demonstrated with large-area
surveys undertaken with Herschel (Negrello et al. 2010, 2017; Conley et al.
2011; Wardlow et al. 2013; Nayyeri et al. 2016), the South Pole Telescope
(SPT; Vieira et al. 2010; Weiss et al. 2013; Mocanu et al. 2013; Vieira et al.
2013), the Atacama Cosmology Telescope (ACT; Marsden et al. 2014), and more
recently with Planck (Planck Collaboration et al. 2015; Cañameras et al. 2015;
Harrington et al. 2016; Berman et al. 2022). Such lens samples are invaluable
for the study of star formation at sub-kpc scales in the early Universe,
thanks to the magnification in angular size afforded by strong lensing.
Without lensing, the angular resolution of Hubble and JWST is only sufficient
to resolve physical sizes at $z>1$ of 1 kpc or greater, so the critical
$10-100$pc scales where star-forming and feedback processes are most relevant
remain inaccessible. These samples offer unique opportunities for statistical
studies of the foreground lensing population (e.g., Eales 2015; Amvrosiadis et
al. 2018)—with implications for the halo mass function and mass density
profiles for individual halos—as the selection function depends only on
intrinsic flux and magnification. This stands in contrast to the complex
biases of lens detection methods based on arc morphology and foreground halo
masses.

DSFGs also benefit from a strongly negative $K$-correction in the
submillimeter (e.g. Blain & Longair 1993), owing to steep Rayleigh-Jeans tails
in the rest-frame far-IR regime of the spectral energy distribution (SED) from
dust emission. Increasing the distance to a galaxy of fixed luminosity would
dim its observable flux at a given wavelength, were it not also the case that
the increased redshift results in capturing shorter-wavelength (and therefore
brighter) parts of the SED. For example, at 850$\micron$ in the observer-
frame, these effects roughly equilibrate, and the flux for an object of fixed
luminosity appears uniform from $z\approx 1-8$ (e.g., Casey et al. 2014). It
is for this reason that DSFGs can be efficiently identified in the
submillimeter regime up to high redshifts (modulo their evolving number
density), and this partially explains the preponderance of submillimeter-
bright galaxies (SMGs; Smail et al. 1997; Hughes et al. 1998; Barger et al.
1998; review by Blain et al. 2002), a term now largely synonymous with DSFGs.

Figure 1: Postage stamp of near-IR images for the 15 PASSAGES objects modeled
in this work. RGB color panels show imaging with Gemini-S $r^{\prime}$ and
$z^{\prime}$ filters (where available) with HST $H$-band (where $H$ is shown
as red, $z^{\prime}$ as green, and $r^{\prime}$ as blue). One exception is
PJ105353, which instead includes an HST F110W image as the green/blue filters
(§2.3). Grayscale panels are those where only $H$-band is available. A white
scalebar in the upper right corner of each panel represents 2′′, and the
target ID is included in the upper left corner. Each panel is north-aligned
and re-centered to best showcase visible lensing features, with the WISE
centroids (or target coordinates for the HST observations) indicated as green
stars (to facilitate easy comparison with Fig. 2). Labeled cyan crosses
indicate the locations of image families, summarized in Table 5 (see also
discussion in Section 3.1).

The Planck All-Sky Survey to Analyze Gravitationally-lensed Extreme Starbursts
(PASSAGES) sample was introduced in Harrington et al. 2016 and Berman et al.
2022, with the scientific aims of exploring the gas fueling and induced
starbursting phases in some of the most apparently IR-luminous objects yet
identified ($\mu L_{\rm IR}=0.1-3.1\times 10^{14}~{}L_{\odot}$, median
$1.2\times 10^{14}~{}L_{\odot}$; also Harrington et al. 2021). Even assuming a
fiducial magnification factor of $\mu\approx 10$, typical for galaxy-scale
strong lenses (e.g. Negrello et al. 2017), they remain intrinsically hyper-
luminous, an order of magnitude higher than local ultra-luminous infrared
galaxies (ULIRGs, Sanders & Mirabel 1996; Lonsdale et al. 2006), which appear
to be the closest low-$z$ analog for DSFGs in terms of observed $L_{\rm IR}$.
Such objects with IR luminosity in excess of $10^{13}L_{\odot}$ have typically
been referred to as hyper-luminous infrared galaxies, or HyLIRGs (Cutri et al.
1994). As suggested by Berman et al. 2022, the selection from an all-sky
facility like Planck can identify the rarest, most extreme star-forming
galaxies. On one hand, strong gravitational lensing may boost the flux of
fainter populations that would otherwise be undetected by an observation (the
well-known magnification bias, e.g., Turner 1980; Turner et al. 1984; Scranton
et al. 2005; Hildebrandt et al. 2009), but it may also be the case that
selection from Planck reveals objects that would be detected within flux
limits even without the modulation by lensing.

A detailed understanding of the PASSAGES objects requires refined and
meticulously-constructed gravitational lens models. As the null geodesics that
light rays from distant objects follow on their path to Earth depend on all
matter—indiscriminate of luminous vs. dark matter—it is impossible to
perfectly describe the exact deformation caused by lensing. While luminous,
baryonic matter can offer great insight into the total underlying mass
distribution, most gravitational lens modeling techniques hinge on
characterizing the distortion of the background and determining what
foreground mass profile is required to induce this distortion. Fortunately,
lensing is achromatic, as photons of all wavelengths follow the same null
geodesics. Images of different wavelengths (which presumably originate in
slightly disparate spatial distributions in the source plane) offer
complementary information by piercing different locations in the foreground
lens. In this context, robust lens models of strongly-lensed sources enable
investigations of $z>1$ galaxies at a (source-plane) spatial resolution that
is currently possible only for local, nearby objects. Therefore, in this
paper, we derive our first multi-wavelength-based lens models for 15 (out of
30) members of PASSAGES (all with spectroscopically-confirmed background
redshifts). In so doing, we use optical and near-infrared information from
Hubble and Gemini, (sub-)millimeter information from the Atacama Large
Millimeter/submillimeter Array, and radio information from the Karl G. Jansky
Very Large Array.

This paper is organized as follows: in Section 2, we introduce the subset of
the Planck lens sample discussed in this work and the relevant near-IR,
submillimeter, and radio observations from which our lens models are derived.
In Section 3, we review our approach to gravitational lens modeling with
lenstool, and discuss the identified lensing features and constraints for each
object. In Section 4, we examine some of the overall properties of this
PASSAGES subsample, including de-lensed luminosity and intrinsic source sizes.
Finally, we summarize our conclusions in Section 5.

In this work, we adopt a flat cosmological model with $H_{0}=70$ km s-1 Mpc-1,
$\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$. At the median source redshift of
this work, $z\approx 2.2$, $1^{\prime\prime}$ corresponds to a physical scale
of 8.27 kpc; at the median lens redshift, $z\approx 0.5$,
$1^{\prime\prime}=6.10$ kpc (Wright 2006).

##  2 Data

###  2.1 Sample selection

PASSAGES consists of 30 high-$z$ objects that were identified from the Planck
Catalogue of Compact Sources (Planck Collaboration et al. 2014) and cross-
identified with Herschel and WISE, the details of which are provided in
Harrington et al. (2016) and Berman et al. (2022). Here, we select a subset of
15 members of the larger current sample, restricting our focus to primarily
those objects with both a simpler lensing morphology, arising from galaxy-
scale or small group-scale lens potentials (Fig. 1), and sufficient multi-
wavelength information. These lower-mass deflectors account for more than 75%
of the total collection of lensed objects from PASSAGES. All of the objects in
this subsample have been imaged with high angular resolution with
optical/near-IR, radio, and (in most cases) submillimeter telescopes. These
data were sufficient to identify the families of lensed arcs and images that
were responsible for the large submillimeter fluxes in these fields, in
addition to the foreground lensing galaxies (which are predominately very
faint at longer wavelengths). We summarize this subsample of objects in Table
1, and provide further details on each source in Appendix A. Lens models for
the remaining objects—most of which are cluster-scale lenses or fields lacking
sub-arcsecond-resolution millimeter imaging—will be described in future works.
Some models have already been presented publicly for PASSAGES, including:
PJ020941.3 (9io9; Geach et al. 2015, 2018; Rivera et al. 2019; Liu et al.
2022; and this work), PJ105322.6 (G145.2+50.9; Frye et al. 2019), PJ105353.0
(G244.8+54.9; Cañameras et al. 2017a, b; Frye et al. 2019; and this work),
PJ112714.5 (G165.7+67.0; Frye et al. 2019; Pascale et al. 2022), PJ132630.3
(NAv1.195; Bussmann et al. 2013; and this work), PJ132934.2/PJ132935.3 (the
“Cosmic Eyebrow;" Díaz-Sánchez et al. 2017), and PJ142823.9 (HBootes03; Borys
et al. 2006). For the 3 previously-modeled targets studied independently in
this work (in addition to PJ113921.7, for which a magnification is given in
Cañameras et al. 2018a), we will examine any inconsistencies in Appendix D.

In deriving the lens models in this work, we take advantage of multi-
wavelength observations from the Large Millimeter Telescope, the Hubble Space
Telescope, Gemini South Observatory, the Karl G. Jansky Very Large Array, and
the Atacama Large Millimeter/submillimeter Array, which we outline in this
section. These observations serve to spectroscopically identify source
redshifts, pinpoint lensed arcs associated with the Planck DSFGs, provide
priors on the foreground lensing mass structure, and offer multiple sightlines
through the deflecting foreground for the purpose of constraining lens mass
models. Expanded details on this workflow are given in Section 3.

Table 1: Summary of PASSAGES objects included in this lens modeling study. A brief description of each object is provided in Appendix A.  ID | RAa | Dec.a | $z_{\rm s}^{b}$ | $z_{\rm lens}^{c}$ | Ref.
---|---|---|---|---|---
PJ011646 | 01:16:46.77 | -24:37:01.9 | 2.125 | 0.555d | 1
PJ014341 | 01:43:41.20 | -01:47:26.0 | 1.096 | 0.594 | 2
PJ020941 | 02:09:41.30 | 00:15:59.0 | 2.554 | 0.202 | 3, 4
PJ022633 | 02:26:33.98 | 23:45:28.3 | 3.120 | 0.41 | 2
PJ030510 | 03:05:10.62 | -30:36:30.3 | 2.263 | 0.5∗ | 2
PJ105353 | 10:53:53.15 | 05:56:18.8 | 3.005 | 1.525 | 4, 5
PJ112713 | 11:27:13.44 | 46:09:24.1 | 1.303 | 0.42 | 2
PJ113805 | 11:38:05.53 | 32:57:56.9 | 2.019 | 0.52 | 2
PJ113921 | 11:39:21.74 | 20:24:50.9 | 2.858 | 0.57 | 2
PJ132630 | 13:26:30.25 | 33:44:07.4 | 2.951 | 0.786 | 2
PJ133634 | 13:36:34.94 | 49:13:13.6 | 3.254 | 0.26 | 2
PJ144653 | 14:46:53.20 | 17:52:33.3 | 1.084 | 0.493 | 2
PJ144958 | 14:49:58.59 | 22:38:36.8 | 2.153 | 0.4∗ | 2
PJ160722 | 16:07:22.77 | 73:47:02.2 | 1.482 | 0.65 | 4
PJ231356 | 23:13:56.64 | 01:09:17.7 | 2.217 | 0.560 | 2

###  2.2 Spectroscopic source redshifts

Spectroscopic redshifts of background objects are critical for deriving
accurate lens models and inferring intrinsic properties in the source plane.
An initial follow-up of candidate lensed objects was performed with 1.1 mm
AzTEC imaging and 3 mm Redshift Search Receiver (RSR) spectroscopy with the
Large Millimeter Telescope (LMT), described in Harrington et al. (2016) and
Berman et al. (2022). The wide bandwidth of RSR ($73-111$ GHz) captures at
least two CO transitions at $z>3.15$, and also in a very narrow window around
$z\approx 2.2$ (Yun et al. 2015), which accounts for 3 members of this
subsample (PJ011646, PJ133634, and PJ144958). The photometrically-supported
redshifts of all other fields were later confirmed spectroscopically to high
precision with the detection of more than 160 redshifted CO and 35 [CI]
emission lines by single-dish facilities—including the Green Bank Telescope
(GBT), the IRAM 30 m telescope, and the Atacama Pathfinder Experiment
(APEX)—as summarized by Harrington et al. (2021).

###  2.3 HST WFC3/IR

All objects in this work were observed with the IR channel of the Wide Field
Camera 3 (WFC3) on the Hubble Space Telescope (HST) during Cycle 24 (Program
GO-14653, PI: J. Lowenthal), except for PJ105353, which was observed in Cycle
23 (Program GO-14223, PI: B. Frye; Frye et al. 2019). Full details of the
observation and reduction of the HST images are given in Lowenthal et al., in
prep., but we provide a brief summary here. The F160W wide-band filter
($\lambda_{\rm effective}=1.54\micron$, $H$-band) was used to image each
member of our sample for one orbit each (with an additional orbit in F110W,
$1.15\micron$, for PJ105353). Each orbit was composed of a 5-point dithering
pattern for 500 seconds each, which were combined during data reduction with
astrodrizzle (final_pixfrac of 0.9). A fiducial $5\sigma$ sensitivity of
$m_{AB}\approx 28.7$ was reached for an unresolved point source. Due to the
sub-pixel dithering setup ($0.572^{\prime\prime}$ spacing), the native pixel
scale was improved by a factor of 2 to $0.065^{\prime\prime}$ (with a point
spread function FWHM = $0.2^{\prime\prime}$). Before combining, the individual
exposures were pipeline-reduced, flat-fielded, and calibrated with standard
Space Telescope Science Institute routines (e.g., Deustua 2016). The absolute
astrometry of the drizzle-combined exposures was subsequently calibrated using
the position of stars within the image frame from the Gaia Data Release 1
catalog (Gaia Collaboration et al. 2016). For fields containing at least 3
Gaia objects detected within the HST/WFC3 field-of-view
($141^{\prime\prime}\times 125^{\prime\prime}$), we can apply a 3-dimensional
astrometric solution. In fields with fewer matched Gaia detections, we apply
only 2-dimensional shifts in RA/Dec.

###  2.4 ALMA 1.1 mm continuum

We obtained Band 6 ($1.1-1.4$ mm, $211-275$ GHz) observations with the Atacama
Large Millimeter/submillimeter Array (ALMA) during Cycle 5 (Program
2017.1.01214.S, PI: M. Yun) for eleven of these 15 targets, executed between
March 22 and May 5, 2018 (with the exception of one executed on August 24,
2018). Additional details of the program are presented in Berman et al.
(2022), but we summarize relevant details here. A default continuum setup was
employed to target rest-frame FIR/sub-mm dust emission in the DSFGs, with a
total bandwidth of 8 GHz ($250-254$ GHz and $266-270$ GHz). A target angular
resolution of $0.4-1.0^{\prime\prime}$ and largest resolvable angular scale of
6′′ was requested, with a target sensitivity of 0.1 mJy over the full
bandwidth. The minimum baseline utilized was $15$ meters, and the maximum
baselines ranged from $484-784$ meters, enabling natural-weighted synthesized
beam sizes of $\sqrt{\theta_{\rm min}\theta_{\rm maj}}\approx
0.4-0.8^{\prime\prime}$. Total on-source integration times ranged from 17 - 21
minutes, resulting in achieved sensitivities of $1\sigma\approx 0.07-0.33$
mJy. Superb weather conditions were achieved for the observations, with the
mean precipitable water vapor ranging between $0.7-1.7$ mm.

As one exception, PJ105353 was observed separately with Band 6 during Cycle 3
(Program 2015.1.01518.S, PI: N. Nesvadba) and with Band 7 ($0.8-1.1$ mm,
$275-373$ GHz) during Cycle 7 (Program 2019.1.01636.S, PI: M. Yun). As the
Band 7 observation had an angular resolution closer to that of our Cycle 5
program than the Band 6 observations (0.35′′ vs. 0.07′′; Cañameras et al.
2017a, b), we opt to use the Band 7 image for our analysis in this work so as
to mitigate the effects of highly disparate beams on derived magnifications.
However, we refer to the Band 6 observation when deriving a lens model; see
Appendix A.

The 1 mm fluxes from ALMA recovered those measured with the LMT/AzTEC
bolometer camera, suggesting that the higher-resolution interferometric images
do not resolve out significant flux at large angular scales (see Figure 8 and
relevant discussion in Berman et al. 2022), ensuring a complete assessment of
lensed arcs111Cluster-scale lenses, most of which are not included in this
analysis, are more likely to be affected by this.. This is important as,
throughout this work, we use lensing magnifications at 1 mm as a proxy for the
total infrared magnification, e.g. in tabulating intrinsic luminosities in
Table 3. In theory, source-plane structure (and in turn magnification) can
vary continuously with infrared wavelength, but this is clearly not feasible
to measure at high spectral resolution in practice.

###  2.5 JVLA 6 GHz

All members of this sub-sample were observed with the Karl G. Jansky Very
Large Array (JVLA) (Program 18A-399, PI: P. Kamieneski). The observations were
carried out between March 29 and April 30, 2018 in 13 execution blocks
totaling 38.9 hours, each targeting two PASSAGES objects. All objects were
observed at 5 cm with C-band ($4-8$ GHz) and full polarization for 1.5 hours
each. With the WIDAR correlator configured to 3-bit sampling, the effective
bandwidth of the two basebands is 4 GHz, centered at 6 GHz. The targets were
scheduled to be observed in pairs at similar right ascensions, alternating
between the two during 3-hour tracks, to improve $uv$-coverage. The most
extended A-configuration was used to provide optimal resolution
($\sqrt{\theta_{\rm min}\theta_{\rm maj}}\approx 0.3-0.7^{\prime\prime}$,
median $0.4^{\prime\prime}$). With maximum and minimum baselines of 36.4 km
and 0.68 km, the largest recoverable angular scale was 8.9′′, which we assume
to be larger than most PASSAGES objects. The 6 GHz continuum data were reduced
using the Common Astronomy Software Applications, casa (McMullin et al. 2007).
Basic flagging and calibration were performed with the VLA Calibration
pipeline (version 2018.1). Each object was first imaged using natural
weighting (maximal sensitivity at the expense of slightly degraded
resolution). After creating a ‘dirty’ image with no deconvolution, the sky
noise level is estimated, before using the CLEAN algorithm (Högbom 1974) to
deconvolve down to a $2\sigma$ threshold. Achieved sensitivities ranged from
$1\sigma=2.6-11.7~{}\mu$Jy (median 2.9 $\mu$Jy). For fields where higher
resolution is desired, we also created images with Briggs weighting (typically
robust $=0.5$).

Photometry for natural-weighted 6 GHz imaging of the background DSFGs is
performed using the flood-filling source extraction software blobcat (Hales et
al. 2012). Gravitational lensing renders the majority of the PASSAGES sample
to be resolved, extended sources, inconsistent with simple 2-dimensional
Gaussian profiles. For this reason, we use the uncorrected integrated surface
brightness when performing continuum photometry (see discussion in §3.3 of
Hales et al.). The total observed 6 GHz flux ($S_{\rm 6GHz}$) is reported in
Table 4.

###  2.6 Gemini $r^{\prime}$ & $z^{\prime}$

We have additionally obtained Gemini Multi-Object Spectrograph (GMOS)
$r^{\prime}$ and $z^{\prime}$ imaging for eight (out of 15) of the members of
this sample as part of Gemini-South programs GS-2018A-Q-216 (PI: J. Lowenthal;
PJ011646, PJ014341, PJ144653, PJ231356), GS-2018B-Q-123 (PI: J. Lowenthal;
PJ020941, PJ030510), and GS-2020A-Q-217 (PI: J. Lowenthal, PJ113921,
PJ144958). In tandem with our $H$-band images from HST, this three-color
imaging (shown in Fig. 1) allows us to identify multiple images, separate
foreground and background objects, and identify lensing groups and clusters
with techniques such as the Red Cluster Sequence, which is effective to at
least $z=1.4$ (Gladders & Yee 2000, 2005). Each target was observed in the
$r^{\prime}$ filter for 1500s to achieve $S/N=10$ for an object with
$r^{\prime}=25.0$, and in the $z^{\prime}$ filter for 2100s to achieve
$S/N=10$ for a $z^{\prime}=23.0$ object. These limiting magnitudes are each
chosen to be one magnitude fainter than a typical early-type L⋆ cluster member
at $z\sim 0.5$ (Gladders & Yee 2005), in order to capture the foreground
objects most likely to be responsible for lensing. Each field was observed
with dithered 300-second exposures to cover chip gaps within the GMOS
$5.5^{\prime}\times 5.5^{\prime}$ field-of-view. The observations have a pixel
scale of 0.080′′. Observations were taken on July 13-14, 2018 for semester
2018A, on October 7, 2018 for semester 2018B, and on February 20, 2020 for
semester 2020A. Absolute image astrometry was corrected in the same manner as
the HST images using Gaia catalog DR1 before creating RGB images (see Fig. 1).
Images were bias-subtracted and flat-fielded using the Data Reduction for
Astronomy from Gemini Observatory North and South
(dragons222https://github.com/GeminiDRSoftware/DRAGONS) software, as will be
described in Cooper et al., in prep.

##  3 Gravitational lens modeling with LENSTOOL

All gravitational lens modeling in this work was performed with the publicly-
available software, lenstool333https://projets.lam.fr/projects/lenstool/wiki
(Kneib et al. 1993, 1996; Jullo et al. 2007; Jullo & Kneib 2009). lenstool
uses parametric forms to model the foreground lens mass distribution that
contributes sufficiently to deflect light from background objects. This
deflection is described by the simple lens equation,

$\boldsymbol{\beta}=\boldsymbol{\theta}-\frac{D_{LS}}{D_{S}}\boldsymbol{\hat{\alpha}(\xi)},$
(1)

where $\boldsymbol{\beta}$, $\boldsymbol{\theta}$, and
$\boldsymbol{\hat{\alpha}(\xi)}$ are the vector quantities describing
(respectively) the intrinsic angular position of the source, the observed
angular position after deflection, and the deflection angle at impact
parameter $\boldsymbol{\xi}$. $D_{S}$ and $D_{LS}$ are the angular diameter
distances to the source plane and from the lens plane to the source plane,
respectively. The deflection angle itself is the integral of surface mass
density in the lens plane:

$\boldsymbol{\hat{\alpha}(\xi)}=\frac{4G}{c^{2}}\int_{\mathds{R}^{2}}d^{2}\xi^{\prime}~{}\Sigma(\boldsymbol{\xi^{\prime}})\boldsymbol{\frac{\xi-\xi^{\prime}}{|\xi-\xi^{\prime
2}|}}$ (2)

where $G$ is the Newtonian gravitational constant and $c$ is the speed of
light. $\Sigma$ is surface mass density evaluated at the location of the
impact parameter $\xi$, after projecting the volumetric density distribution
on to the lens plane (Schneider et al. 1992). Thus, provided the locations
$\boldsymbol{\theta}$ (and redshifts) of multiple images identified visually
by the user, one can use this formalism to model the distribution of mass in
the lens plane. The process of isolating image families is usually iterative,
where the model may predict unidentified image locations. When confirmed,
these added images can be incorporated into the model for further refinement.
However, for many of the lensed systems in this sample, which are
predominantly galaxy-scale lenses, identifying image families can be
straightforward. We describe this step in more detail in the following
section, and label the image systems in Fig. 1.

The only information included in the derived lens models for lenstool are the
multiple image catalogs supplied by the user, and the parametric mass models
are likewise defined by the user. The input image catalogs we constructed are
provided in Table 5. For some lenses, this can grossly oversimplify the
information provided by imaging itself, such as the extended ring-like
features seen in ALMA and HST observations for many PASSAGES objects (e.g.,
Fig. 1). When these arcs and rings are clumpy, it may be possible to identify
sub-images, or smaller-scale features that appear inside the broader multiple
images. For interferometric imaging, this may be challenging, as certain
spatial scales are filtered out, and correlated noise can introduce spurious
point-like sources. As lens models for the PASSAGES sample are continually
refined in the future, we will ideally incorporate information about extended
rings, which often closely trace the critical curve.

Lens models for most cluster-scale members of the PASSAGES sample will be
formulated using both parametric modeling (e.g. lenstool) and light-traces-
mass models, in an approach similar to Zitrin et al. (2009, 2015) and Frye et
al. (2019)444See also Broadhurst et al. 2005, Acebron et al. 2018, and Cibirka
et al. 2018. , the latter of which covers several members of the full PASSAGES
sample that were identified independently by Cañameras et al. 2015. Here, we
narrow our focus primarily to the galaxy-scale and group-scale lenses, which
are more likely to be adequately described by simpler parametric mass models
than larger, cluster-scale lenses.

|
---|---
|
|
|
|
|
Figure 2: Gravitational lens models for ALMA 1.1 mm (left column) and JVLA 6 GHz (right column) continuum imaging. Each row shows the observed image for different PASSAGES members (left), the model-derived image plane structure as a result of ray-tracing the observed image to the source plane and then back to the image plane, with gold contours showing the observed data for comparison (center), and the source-plane reconstruction (for all multiple images combined), zoomed-in to show more detail (right). All images are shown with the same colorscale limits. Red stars in each panel indicate the location of the WISE centroid, with each set of panels centered on the phase center of the interferometric pointing. The synthesized beam is shown in cyan in the lower left of the left panels. The model-derived caustic and critical curves are shown in the center panel in cyan and purple, respectively. The contours in the middle panel are typically chosen to range evenly from $3\sigma$ to the 99th percentile of the image (unless this value is less than $5\sigma$). The ALMA image shown for PJ105353 is at Band 7 to ensure a resolution more consistent with that of the other fields.  |
---|---
|
|
|
|
|
|
Figure 2: (Continued.)  |
---|---
|
Figure 2: (Continued.) The 6 GHz image for PJ231356 is not modeled, as the
emission is dominated by an apparent AGN and jet in the foreground lensing
elliptical, and any background lensing structure is not easily recovered.

###  3.1 Visual identification of multiple image families

Identifying multiple image families with long-wavelength, high-resolution
radio/millimeter imaging can be easier than at shorter wavelengths, as the
ubiquity of dust in the background DSFGs leads to strong obscuration in the
optical regime and a much brighter rest-frame far-IR arising from the peak of
the dust SED. Moreover, foreground galaxies are largely filtered out at longer
wavelengths. There is a notable exception of radio AGN and jets propagating
from the massive elliptical galaxies in the lens plane (e.g., Fanaroff & Riley
1974, Smith et al. 1986), which can contaminate image classification. For this
reason, our method of identifying image families relies on the combined
information of optical/near-IR imaging from HST, millimeter imaging with ALMA,
and radio imaging from the JVLA. This coverage spanning 0.6 $\mu$m to 6 cm
gives us a large number of unique sightlines through foreground mass
distributions, which we use to constrain our models.

The interferometric observations by ALMA and the JVLA do not directly probe
the sky brightness pattern, but rather the complex visibility, which is a
2-dimensional Fourier transform of the brightness pattern. Translating the
visibilities into an image introduces correlated noise, which is a function of
the sampling pattern of the interferometer. For this reason, a number of works
have performed gravitational lens modeling in the visibility $uv$-plane
instead of the surface brightness plane (including Bussmann et al. 2012, 2015;
Hezaveh et al. 2013; Spilker et al. 2016). In this work, as we are merging
information from traditional optical telescopes and radio interferometers, we
apply lens modeling to the transformed image plane. This has the additional
benefit of improved computational efficiency. As a proof of concept, Dye et
al. (2018) used both image-plane and visibility-plane lens modeling towards
ALMA imaging of six Herschel-detected lenses for direct comparison, and found
minimal difference in the derived parameters and source reconstructions, and
derived total magnifications were in agreement within 1$\sigma$ uncertainties.
Dye et al. note that there might be greater disparities in cases where the
visibility sampling is sparse, but we do not expect this to be a major concern
in our case. Both ALMA and JVLA observations make use of a large number of
antennas and long integration times to maximize $uv$-coverage (and as
mentioned previously, JVLA observations of fields were carried out in pairs
during 3-hour tracks to ensure more distributed $uv$-sampling).

Each lensed image family of multiplicity $n$ provides $2(n-1)$ constraints to
be used in the optimization555The right ascension and declination of each
image provide 2 constraints for each image, but one member does not add
constraints because the actual source position of all images is unknown.. If
relative fluxes (or equivalently, solid angle) of each family member are
included, this adds an extra $n-1$ constraints (Blandford & Narayan 1992). For
most models herein, we use only image positions, which helps to minimize the
systematic uncertainty introduced in the model. However, some lensed systems
on galaxy- or group-scales with few multiple image constraints are poorly
suited for strong lens modeling (see e.g. Limousin et al. 2009; Verdugo et al.
2014). In such cases, it may be necessary to include additional information in
the model, which we discuss in the following section where applicable.

###  3.2 Lens model mass profile optimization

lenstool employs Bayesian Markov chain Monte Carlo (MCMC) methods to
thoroughly explore the posterior distribution for each supplied parameter in
order to find a best-fit solution. While this method is more impervious to
local $\chi^{2}$ minima, it is more computationally expensive than other
minimization techniques when the parameter space is complex. Because of the
possible parameter degeneracies and non-Gaussian distributions that can arise
in the underlying posterior, it is usually highly beneficial to employ MCMC.
Performing this optimization in the source plane (i.e. minimizing separation
between image family members once ray-traced to the source plane, rather than
making the comparison in the image plane) can reduce the requisite
computational time. However, this also results in lower precision and poorer
estimation of model uncertainties. As the PASSAGES objects in this work are
primarily galaxy-scale lenses, fewer free parameters are involved than
cluster-scale systems, so we perform all optimization in the image plane
unless otherwise noted. Another quantification of goodness-of-fit can be the
root-mean-square of the image-plane observed vs. model-reconstructed image
locations, as this is independent of the positional uncertainty, therefore
facilitating an easier comparison of different lens models (e.g. Caminha et
al. 2016).

One of the priors provided to lenstool is the functional form of each mass
potential, including the relevant free parameters (each of which has their own
prior distribution, specified by the user). For these models, we use a
singular isothermal ellipsoid (SIE) potential as a basis for all mass profiles
(Kormann et al. 1994). With this model, the normalized surface mass density,
or convergence $\kappa$, takes the form

$\kappa(x_{1},x_{2})=\frac{\sqrt{f}}{2\sqrt{x_{1}^{2}+f^{2}x_{2}^{2}},}$ (3)

where $f\equiv b/a$ is the axis ratio (for semi-major and semi-minor axes $a$
and $b$; $0<f\leq 1$) and $x_{1}$ and $x_{2}$ are normalized Cartesian
coordinates. As noted by Treu (2010), the SIE—a generalization of the singular
isothermal sphere (SIS) where 3-dimensional density follows $\rho\sim
r^{-2}$—appears to be the simplest profile that effectively describes galaxy-
scale strong lensing configurations (see also work by Treu & Koopmans 2004;
Koopmans et al. 2006, 2009; Barnabè et al. 2009). A number of studies of
lensed DSFGs have had success using this profile, including Fu et al. (2012),
Bussmann et al. (2012, 2015), Hezaveh et al. (2013), Calanog et al. (2014),
and Spilker et al. (2016). The profile can be fully parameterized by position
($\alpha$ and $\delta$, or $x$ and $y$), projected 2-dimensional ellipticity
$e=(a^{2}-b^{2})/(a^{2}+b^{2})$, position angle (PA or $\theta$; measured
counterclockwise from east), and velocity dispersion, $\sigma_{\rm SIE}$. This
velocity dispersion is generally in good agreement with observable stellar
velocity dispersions for elliptical galaxy lenses (e.g., Bolton et al. 2008
for $\sigma\approx 175-400$ km s-1). For the majority of the galaxy-scale
lenses in the PASSAGES sample, the arcs and rings are visible in the
intermediate range between core and cut radii (e.g., Table 2), and so these
parameters of pseudo-isothermal elliptical mass distributions (PIEMD; Kassiola
& Kovner 1993) and dual pseudo-isothermal elliptical distributions (dPIE;
Elíasdóttir et al. 2007) would be poorly constrained.

As galaxy-scale systems have comparatively fewer constraints than clusters
that are rich with multiply-imaged arcs, the number of available free
parameters is usually quite limited. In some cases, this necessitates a
simple, single SIE potential, with free parameters including the location of
the profile center, the velocity dispersion, and the ellipticity (and
associated position angle). Where possible (and appropriate), a secondary
singular isothermal ellipsoidal/spheroidal potential can be added to account
for the effect of any secondary deflectors. This 2nd-order correction is
similar to including an external shear field induced by underlying large scale
structure (see review of weak lensing by Bartelmann & Schneider 2001), which
is a common approach to modeling galaxy-scale lenses, including for a number
of PASSAGES objects (e.g., Geach et al. 2015, 2018; Rivera et al. 2019). In
some galaxy-scale cases, the location of the profile center is fixed to the
centroid of the visible baryonic component revealed by HST, which reduces the
number of free parameters666For cluster-scale lenses, the lens position is
kept as a free parameter, although often the coordinates of the BCG may be
assumed to represent the cluster center.. An upper limit is generally imposed
on the ellipticity of $e<0.75$, following Acebron et al. (2017). Theoretical
predictions by Despali et al. (2017) indicate that $\approx 95\%$ of galaxy-
scale, $10^{12}~{}M_{\odot}$ halos at $z\sim 0.5$ will have a projected
ellipticity less than $e=0.5$. In our sample, there are several cases where
the optimized lens models have ellipticity values larger than expected (up to
$e\sim 0.7$). This suggests that these models may benefit from additional
complexity and/or substructure (e.g. Johnson et al. 2014), but this is often
limited by the small number of constraints presently available.

The results of our MCMC optimization are summarized in Table 6, including the
positional center of each mass profile (relative to reference coordinates),
and the ellipticity $e$, position angle $\theta$, and velocity dispersion
$\sigma$ in km s-1. For each lensing field, we provide at least two optimized
solutions, labeled median and best. The median solution comprises the median
values of the MCMC-sampled posterior for each parameter, with associated
asymmetric uncertainties from the inner 68% confidence interval of the
distribution, whereas the best solution is only those MCMC realizations with
the highest likelihood (or lowest $\chi^{2}$). In some cases, the mode
solution is also provided, comprising the peak of the probability density for
each parameter. This latter solution may be preferable in situations where the
posterior distribution is multimodal, for which the median may lie in troughs
between multiple peaks in parameter space. Additionally, the mode solution can
be useful when a parameter is only weakly constrained in the model, such that
the choice of parameter bounds can have a strong impact on the distribution.
Here, the median might be sensitive to these bounds, while the mode might be
more resistant.

In general, we select the median of the posterior distribution, which is less
prone to local likelihood maxima, where small perturbations to any parameter
can drastically reduce the likelihood. However, in the ideal case, the median,
mode, and best solutions would be in agreement within error, but this requires
a parameter space that is sufficiently simple. Altogether, the median
$\chi^{2}$ value for all median solutions777For fields where the mode solution
is also provide, we use the higher $\chi^{2}$ value of the two. is
$\tilde{\chi}_{\rm med}^{2}={1.9}$, and the median for the best solutions is
$\tilde{\chi}_{\rm best}^{2}={0.3}$. This suggests that the goodness-of-fit of
our lens models are reasonable, in aggregate, although the individual
$\chi^{2}$ values are influenced by the complexity of the foreground
environment and the amount of lensing evidence available from current
observations. The lensing models we select for analysis are shown in Fig. 2,
which depicts the observed image-plane structure relative to the model-
reconstructed image-plane and source-plane structure, the latter of which we
discuss in the next section. Qualitatively, the image-plane reconstructions
are in agreement with the observations, mainly in that they reproduce the
correct number of multiple images.

###  3.3 Source-plane reconstruction

With an optimized potential model in hand, we can invert the observed arcs to
approximate the unlensed brightness distribution in the source plane using the
cleanlens function of lenstool. Fundamentally, this inversion can be quite
trivial, in that it consists only of ray-tracing a (usually sub-sampled) image
plane distribution to the lens plane, computing the deflection vectors based
on the lens mass potential, and continuing until the source-plane redshift is
reached. Typically, sub-pixel sampling is also used for the source plane. The
direction and magnitude of the deflection are determined by a surface integral
of the surface mass density, weighted by the distance from where each light
ray intersects the lens potential (see e.g., Schneider et al. 1992, Refsdal &
Surdej 1994).

While the optimization and goodness-of-fit measurement with lenstool is done
in the image plane, one of the tests that we use to evaluate the goodness-of-
fit of a lens model is evaluating the spatial coincidence of the source plane
reconstruction of individual multiple images. Every observation has an
intrinsic point spread function (PSF) with which the image-plane structure is
convolved, which is carried through as a distorted, non-uniform PSF in the
source plane. This will be slightly different for each individual multiple
image (or even dramatically different, in the case of a source near the
lensing caustics; see e.g. Sharma et al. 2018, 2021). These individual
reconstructions are compared visually to ensure that the total reconstruction
of all image-plane light is appropriate and self-consistent.

Table 2: Summary of strong lensing properties derived from the models
described in Table 6, including weighted total magnification $\mu_{\rm tot}$
for different wavelengths and the equivalent Einstein radius $\theta_{\rm
Ein,eq}$. Under the Literature heading, these values are listed for other
PASSAGES objects not included in this study that have published lens models.
Here, magnifications are derived for 870 $\mu$m, and Einstein radii are
defined as the radius of the circle with equivalent area to that inside the
critical curve (see discussion in Section 3.5).

Field | $\mu_{\rm 1mm}$ | $\mu_{\rm 6GHz}$ | $\theta_{\rm Ein}$ | $M_{<\theta_{\rm E}}$
---|---|---|---|---
|  |  | [′′] | [$10^{11}M_{\odot}$]
PJ011646 | $7.0^{+0.3}_{-0.3}$ | $7.0^{+0.4}_{-0.3}$ | $2.37^{+0.04}_{-0.03}$ | $15.4\pm 0.3$
PJ014341 | $7.1^{+2.2}_{-2.2}$ | $5.9^{+2.0}_{-1.4}$ | $0.53^{+0.07}_{-0.09}$ | $1.3\pm 0.3$
PJ020941 | $10.5^{+4.2}_{-3.2}$ | $12.0^{+3.7}_{-3.3}$ | $2.55^{+0.14}_{-0.20}$ | $6.6\pm 0.6$
PJ022633 | — | $28.2^{+10.7}_{-11.9}$ | $3.61^{+0.36}_{-0.30}$ | $24.7\pm 3.2$
PJ030510 | $2.2^{+0.6}_{-0.5}$ | $2.6^{+1.0}_{-1.0}$ | $0.58^{+0.12}_{-0.15}$ | $0.8\pm 0.3^{*}$
|  |  |  | $[0.5-1.1]$
PJ105353 | ${7.6^{+0.5}_{-0.8}}^{**}$ | $7.9^{+0.6}_{-0.8}$ | $0.71^{+0.03}_{-0.05}$ | $3.6\pm 0.3$
PJ112713 | — | $4.5^{+0.5}_{-0.6}$ | $0.58^{+0.07}_{-0.04}$ | $0.8\pm 0.1$
PJ113805 | $2.8^{+0.4}_{-0.3}$ | $4.4^{+1.2}_{-1.0}$ | $0.40^{+0.06}_{-0.04}$ | $0.4\pm 0.1$
PJ113921 | $4.8^{+1.3}_{-0.9}$ | $8.6^{+2.8}_{-2.6}$ | $0.71^{+0.47}_{-0.39}$ | $1.3\pm 1.1$
PJ132630 | $4.3^{+0.7}_{-0.5}$ | $11.4^{+4.5}_{-2.2}$ | $1.78^{+0.21}_{-0.14}$ | $11.0\pm 1.5$
PJ133634 | — | $8.3^{+4.2}_{-1.9}$ | $1.17^{+0.11}_{-0.08}$ | $1.7\pm 0.2$
PJ144653 | $3.9^{+0.6}_{-0.7}$ | $7.7^{+1.5}_{-2.0}$ | $0.82^{+0.05}_{-0.04}$ | $1.6\pm 0.1$
PJ144958 | $10.1^{+5.7}_{-4.7}$ | $13.6^{+15.7}_{-6.7}$ | $6.52^{+1.31}_{-1.16}$ | $84.0\pm 22.5^{*}$
|  |  |  | $[43.3-126.0]$
PJ160722 | — | $6.8^{+2.1}_{-1.3}$ | $1.00^{+0.11}_{-0.08}$ | $4.0\pm 0.5$
PJ231356 | $6.1^{+1.2}_{-1.2}$ | —† | $2.05^{+0.26}_{-0.23}$ | $11.5\pm 1.9$
Literature | $\mu_{870\mu{\rm m}}$ |  | $\theta_{\rm Ein}$ | $M_{<\theta_{\rm E}}$
|  |  | [′′] | [$10^{11}M_{\odot}$]
PJ105322a | $7.6\pm 0.5$ |  | $5.9$ | 120.3
PJ112714b | $29.4\pm 5.9$ |  | $13$ | 283.7
PJ132302c | $11.2\pm 0.7$ |  | $\lesssim 13$ | $\lesssim 380.1$
PJ132934d | $11\pm 2$ |  | $11.0\pm 0.4$ | $267.9\pm 13.8$
PJ142823e | $3.0\pm 1.5$ |  | $0.10\pm 0.03$ | $0.13\pm 0.06$
PJ154432f | $14.7\pm 0.8$ |  | $7$ | 152.6
PJ160918g | $15.4\pm 1.0$ |  | — | —

Note. — a J2000 object positions are derived from WISE cross-matching. b
Source redshifts are determined from CO detections, reported in the listed
references. c Foreground lens redshifts (photometric and spectroscopic) are
taken from the supplied references. d Spectroscopic foreground redshift
determined with the Multi Unit Spectroscopic Explorer (MUSE) on the Very Large
Telescope (VLT). ∗ Preliminary photometric redshift estimate (Cooper et al. in
prep.). For the lensing mass estimates from Einstein radii in §3.5, we assume
a fiducial uncertainty in the lens redshift of $\sigma_{z}=\pm 0.2$ (given the
inner 68% range of lens redshifts for the PASSAGES sample, $\approx 0.3-0.7$;
Harrington et al. 2021; Berman et al. 2022); see Table 2. Properties such as
magnification are not impacted by this uncertainty, as parameters
encapsulating the mass of the lens are largely degenerate with lens redshift.
References: (1) Kamieneski et al., in prep.; (2) Berman et al. 2022;
Harrington et al. 2021; (3) Geach et al. 2015; (4) Harrington et al. 2016; (5)
Cañameras et al. 2015, 2017a and Frye et al. 2019.

Note. — ∗ Only a preliminary photometric redshift measurement (see Table 1),
so approximate bounds of enclosed masses are also given for a set of redshifts
representative of the range of the full sample. ∗∗ 870 $\mu$m (Band 7)
measurement (see §2.4). The 1 mm magnification is $\mu\approx 11.8\pm 2.1$ for
the $0.07^{\prime\prime}$ image, but the 870 $\mu$m value is adopted for
$\mu_{\rm dust}$ to more closely match the imaging limitations for the other
targets. † 6 GHz image of PJ231356 is not resolved from a possible radio jet,
so it is not currently possible to estimate magnification of the background
DSFG. a Also known as PLCK_G145.2+50.9; Cañameras et al. (2018a); Frye et al.
(2019). b Also known as PLCK_G165.7+67.0; Cañameras et al. (2018a); Frye et
al. (2019). c Also known as PLCK_G113.7+61.0; Cañameras et al. (2018a).
Einstein radius upper limit estimated from VLA FIRST survey image (Becker et
al. 1995). d Also known as the Cosmic Eyebrow; Díaz-Sánchez et al. (2017);
Dannerbauer et al. (2019). e Also known as HBoötes3; Bussmann et al. (2013);
Borys et al. (2006). f Also known as PLCK_G080.2+49.8; Cañameras et al.
(2018a); Frye et al. (2019). g Also known as PLCK_G092.5+42.9; Cañameras et
al. (2018a); Frye et al. (2019).

###  3.4 Total magnification factors

Total magnification of each lens is determined as the ratio of the total solid
angle subtended in the image plane to that of the source plane (as
reconstructed by our lens model). These measurements are wavelength-dependent,
as the source-plane structure varies with wavelength. To estimate model-
dependent magnification uncertainties, we follow an approach similar to Sharon
et al. 2012. 300 samples are drawn from the full MCMC-sampled posterior
distribution for each lens mass profile parameter optimized with lenstool, and
the observed image-plane arcs are ray-traced into the source plane for each
realization, each resulting in a slightly different magnification, measured as
the area above a $3\sigma$ signal-to-noise threshold (for both the source
plane and image plane). The median of the distribution (and uncertainty in the
form of 16th and 84th percentiles) are reported in Table 2 for each object.
However, we caution that these statistical values can underestimate larger
systematic uncertainties—such as mis-identified multiple images and
inappropriate choices of parameterized mass distributions—which are difficult
to properly take into account. For example, Limousin et al. (2016) considered
the effect of using cored vs. non-cored cluster mass models—i.e., where the
isothermal mass profile of cluster members is modified within an inner core
radius. The authors found that systematic uncertainties in the magnification
could be nearly an order of magnitude greater than statistical model-specific
uncertainties. Here, this is unlikely to have a strong impact on the
predominantly galaxy-scale lenses, which are influenced more by mass at
intermediate radii than the typical core radius ($\sim 100$ pc) or cut radius
($\sim 30-100$ kpc), as noted by Cañameras et al. (2017a). Regardless, this
provides some motivation for model-independent characterizations, such as
Wagner & Bartelmann (2016) and Wagner & Tessore (2018), as magnifications can
vary across different works. We reserve the discussion of the distribution in
magnification factors in the PASSAGES sample for Section 4.1.

###  3.5 Measuring Einstein radii

The Einstein radius $\theta_{E}$ is a characteristic scale relating directly
to the total enclosed mass of the foreground lens inside a projected circle in
the sky $M(<\theta_{E})$, with some additional dependence on the redshift
geometry of the lens and source planes:

$\theta_{E}=\sqrt{\frac{4GM(<\theta_{E})}{c^{2}}\frac{D_{LS}}{D_{L}D_{S}}},$
(4)

where $D_{L}$, $D_{S}$, and $D_{LS}$ are the angular diameter distances to the
lens plane, to the source plane, and from the lens to source plane,
respectively (Narayan & Bartelmann 1996). $\theta_{E}$ is usually quite
consistent across different lens models (e.g., Kochanek 1991; Wambsganss &
Paczynski 1994), making it a valuable measurement. Observationally, the
Einstein radius is comparable to half the typical separation between multiple
images.

There are a number of methods for deriving $\theta_{E}$ from lens models.
Here, we opt for the equivalent Einstein radius, $\theta_{Ein,eq}$, defined as
the radius enclosing an average surface mass density equal to the critical
surface density (or equivalently, where the average convergence $\kappa$ is
equal to 1; see Richard et al. 2010; Zitrin et al. 2011). Other definitions
include the effective Einstein radius (Redlich et al. 2012), or the radius
enclosing an area equal to that interior to the tangential critical curve, and
the median Einstein radius, defined as the median distance from points on the
tangential critical to the lens centroid (Meneghetti et al. 2011, 2013). Since
all definitions are more or less connected to the area interior to the
critical curves, they are also related to the lensing cross section, which is
the source-plane area inside the caustics (or sometimes defined as the area of
the region where the magnification is greater than a certain threshold).

We report the Einstein radii $\theta_{\rm Ein}$ in Table 2, as the median
value and asymmetric $1\sigma$ uncertainties. The Einstein radii are computed
analogously to the magnification, where 100 samples are drawn from the
modeling-derived posterior distribution, with a new convergence map created
for each iteration. Centered on the peak pixel of each convergence map (the
point of greatest surface mass density), the radius is increased linearly
until the average convergence crosses below 1. The distribution of Einstein
radii for this sample is discussed in §4.5.

###  3.6 Source-plane size

As our lens-modeling algorithm does not assume a parametric form of the
source-plane structure, the intrinsic size of the far-IR continuum-emitting
region is measured from the source-plane reconstruction using casa imfit. As a
caveat, this measurement is perhaps best made in the visibility-plane in order
to deconvolve the effect of the synthesized beam. Since our source-plane
reconstruction retains a distorted version of the beam that varies across the
source, removing this is not trivial without a $uv$-based forward-modeling
approach (e.g. visilens by Hezaveh et al. 2013; Spilker et al. 2016 or
uvmcmcfit by Bussmann et al. 2012, 2013, 2015). Moreover, since this lens
modeling approach does not impose a restriction on the nature of the source
plane emission, small but non-negligible offsets might be expected between the
source-plane positions of reconstructions from different lensed images (§3.2).
We measure the effective radius from the reconstruction of all lensed images
together, so these offsets can artificially change the apparent size888It is
not advisable here to measure sizes based on the reconstruction of individual
lensed images, as some images may not include the full source-plane structure.
As an example, for galaxies where only a small region crosses into a
quadruply-imaged region, two of the four images will contain only this region
alone (while the remaining two images will show the full galaxy).. To capture
some of this uncertainty from lens modeling, we randomly select 10 of the MCMC
lens modeling iterations to create different source-plane reconstructions,
measure effective radii independently, and take the mean and standard
deviation of this distribution. Lastly, we convert these angular measurements
to physical sizes in kpc given the redshift of the source plane. Our observing
setup with ALMA is not markedly different from that of Spilker et al. (2016)
and Bussmann et al. (2015), for example, and the relative uncertainty in
effective radius for this work is likewise similar. For this reason, we
consider a direct comparison to be feasible.

However, as an added complication in comparing to the broader literature,
there are different approaches to quantifying source size. At optical
wavelengths, it is standard to use the effective radius or half-light radius
within which half the object’s light is contained (under an assumption of
circular symmetry). This in turn is determined through fitting a Sérsic model
to the radial profile (Sérsic 1963; e.g., with GALFIT, Peng et al. 2002), for
which the degree of central concentration is parameterized by a Sérsic index
$n$. For radio imaging, it’s more common to assume a two-dimensional Gaussian
functional form (equivalent to a Sérsic model with $n=0.5$). The full-width at
half-maximum of the Gaussian is related to the effective radius as ${\rm
FWHM}=2\times R_{\rm e}$. In our method, we therefore use the geometric mean
of the semi-major and semi-minor axes as the effective radius 999 Murphy et
al. (2017) additionally demonstrate that the effective radius of marginally
resolved sources can be connected to the deconvolved major axis of an inclined
disk as $\theta_{\rm maj}\approx 2.430\times R_{\rm e}$. . While the lensing-
reconstructed PSF will vary in dimensions for the source plane, we do not
account for this here for these global size measurements, as we expect this
effect to be smaller than the uncertainties introduced from modeling.

Figure 3: Apparent IR luminosity ($\mu_{\rm IR}L_{\rm IR}$) vs. magnification
$\mu$ at 1100 $\micron$ (or 870$\micron$ if not available) and at 6 GHz (see
Table 2). Dashed lines indicate intrinsic, magnification-corrected star
formation rates from the Kennicutt 1998 calibration under assumption of a
Kroupa initial mass function (IMF), ${\rm SFR}=L_{\rm IR}/(9.4\times
10^{9})M_{\odot}{\rm yr}^{-1}$. Purple connecting lines are drawn between sub-
mm and radio magnifications where both are measured. While the far-IR
magnifications are perhaps more representative of the factor modifying the
total IR luminosity than those for the radio continuum, the two are largely
consistent, and so the radio magnifications can still offer a representative
correction to estimate intrinsic star formation rates.

##  4 Discussion

###  4.1 PASSAGES sources have modest magnifications

Eleven out of the 15 objects in this work for which we present new lensing
models were previously undiscovered101010We note that some were since
identified independently by Trombetti et al. (2021). before PASSAGES
(Harrington et al. 2016; Berman et al. 2022). These new lens models therefore
provide many of the first measurements of magnification at any wavelength for
PASSAGES. In our analysis, we utilized HST (1.6µm $0.15^{\prime\prime}$
resolution) and JVLA (6 GHz, $0.3-0.7^{\prime\prime}$ resolution) observations
for all members. For 11 out of 15, we also used ALMA imaging (1 mm,
$0.4-0.8^{\prime\prime}$ resolution), and for 8 out of 15, we used Gemini
$r^{\prime}$ and $z^{\prime}$ observations. In future work, we intend to
refine and apply these lens models towards multi-J CO and [CI] high angular
resolution datasets from ongoing campaigns with ALMA, the Submillimeter Array
(SMA), and the Northern Extended Millimeter Array (NOEMA). These will be used
to interpret the kinematic structure of different gas components and build a
systematic set of differential magnification corrections. After accounting for
this work, only 4 high-$z$ members of the PASSAGES sample lack published lens
models: PJ074851.7 (Garcia Diaz et al., in prep.), PJ084650.1 (Foo et al., in
prep.), and PJ132217.5 (clusters with background DSFGs at $z=2.76$, $z=2.66$,
and $z=2.07$), and the presumed galaxy-galaxy lens PJ141230.5 ($z=3.31$),
which currently lacks high-resolution mm/radio follow-up and does not have
discernible lensing features in the near-IR.

The PASSAGES objects are detected in near-IR imaging with Hubble (Lowenthal et
al. in prep.), in contrast to the often optical faintness of DSFGs, especially
at $z>3$ (Franco et al. 2018; Wang et al. 2019; Barrufet et al. 2022; Nelson
et al. 2022; Zavala et al. 2022). This latter point is due in part to the
fraction of obscured star formation increasing with SFR (Whitaker et al. 2017)
and to the more dust-attenuated shorter wavelengths probed at higher-$z$.
Given the availability of lensing evidence in the rest-frame optical, in
complement to the rest-frame far-IR, our approach in using lenstool is less
vulnerable to the pitfalls of modeling using solely interferometric imaging,
which would typically necessitate direct modeling of $uv$-plane visibilities.

Figure 4: Histogram of magnifications for the PASSAGES sample (derived from
lens models presented in this work and in the references included in Table 2),
compared with Herschel lenses from Bussmann et al. (2013, 2015), Dye et al.
(2015, 2018), and Massardi et al. (2018); and with SPT lenses from Spilker et
al. (2016). “Barcode" lines at the bottom indicate the actual distribution
values, for which kernel density estimations are shown as solid curves. The
smoothing kernel is chosen according to Scott’s Rule (Scott 1992). We caution
that, crucially, lens models have not been developed for all PASSAGES members,
and that the magnification distribution may be skewed towards lower values due
to several more complex, cluster-scale halos being excluded from this analysis
at this juncture. Nonetheless, from the information currently available, the
PASSAGES selection does not appear to be subject to drastically different
lensing systematics from similar works, as discussed in §4.1.

The distribution of magnifications for (sub-)mm continuum in this work is
largely consistent with those from similar lensed surveys defined from
Herschel and SPT, as shown in Figure 4. This may be surprising given the
larger apparent luminosities of PASSAGES. However, we caution that there are
still several objects in the PASSAGES sample that lack detailed lens modeling,
so this comparison is not robust at this point. In addition, there are
external factors at play in this distribution: in this work, we have largely
reserved the time-intensive modeling of cluster lenses for subsequent
analyses. Likewise, the lenses identified with Herschel preferentially exclude
cluster-scale deflectors, as noted by Bussmann et al. (2015), due to their
deconvolution of the Herschel Spectral and Photometric Imaging REceiver
(SPIRE; Griffin et al. 2010) PSF (beam FWHM $18-36^{\prime\prime}$), which
removes extended objects close to the size of the beam. The SPT sample,
however, contained four cluster-lensed objects, but these were not analyzed by
Spilker et al. (2016) as not all images were captured within the ALMA field-
of-view. Both this work and others are limited to objects with sufficient
high-resolution, high-sensitivity data that are suitable for modeling.

The resolution of data can have some effect on the derived magnifications
(e.g., Hezaveh et al. 2012). For example, the current $\sim
0.4^{\prime\prime}$ resolution comfortably resolves the lensed images in most
cases, but does not offer much detail of substructure inside each of the
images. Higher resolution ($\sim 0.1^{\prime\prime}$) might reveal chance
alignments of compact, bright star-forming clumps with high-magnification
areas near the caustic curve. Lower resolutions that are blind to these
details could bias the magnification factors. As an added example, two well-
separated clumps lying on opposite sides of the critical curve might result in
an over-estimated magnification factor if they are blended in lower-resolution
imaging, thus appearing to cross the critical curve contiguously. Such effects
can not be mitigated completely, but future higher-resolution observations
will be key to ruling out these biases.

Bussmann et al. (2013) remarked on lower-than-expected magnifications for
Herschel-identified lenses, given the objects’ apparent luminosities.
Theoretical predictions from Wardlow et al. (2013) based on submillimeter
galaxy number counts to estimate the average magnification factor as a
function of 500$\micron$ flux density were consistently higher than model-
derived magnifications in all but three cases. In their case, an assumption of
a single Sérsic profile to represent the source plane might systematically
underestimate magnifications. Bussmann et al. (2015) offered an explanation
that the resolution of their SMA and ALMA observations ($\theta\approx
0.5^{\prime\prime}$) was insufficient to spatially resolve emission within
individual lensed images, leading to half-light radius estimates skewed to
larger values, in turn giving smaller magnification factors. Indeed, later
works showed that higher-resolution imaging ($\approx 100$ milliarcsecond) of
the same objects (Dye et al. 2015; Rybak et al. 2015; Tamura et al. 2015)
could lead to a factor of $1.5-2$ increase in magnification (Bussmann et al.
2015).

###  4.2 Intrinsic IR luminosities of highly-magnified DSFGs and their
relation to the main sequence

Apart from a different modeling approach, it is also possible that the all-sky
Planck selection used in PASSAGES (contrasted with the survey area-limited
Herschel and SPT counterparts) is simply more advantageous for selecting the
rarest objects: intrinsically extreme starbursts that are also serendipitously
subject to strong lensing. While it is clear that identifying lensing
candidates by setting simple submillimeter flux thresholds is very efficient,
there is not an obvious relation between flux and magnification for the
resultant sample. With a larger ensemble, it will be logical to compare both
the distribution functions of magnifications and the luminosity functions of
the de-lensed objects with unlensed populations, but that is beyond the scope
of this current work. In the next section, we discuss the apparent balance for
a source size that is large enough to sustain intense star formation but small
enough to avoid the dilution of magnification due to size bias.

Using the lens model-derived magnifications at 1.1 mm, we can derive
intrinsic, de-lensed IR luminosities (given in Table 3 and shown in Figure 3),
which range from $L_{\rm IR}=0.2-5.9\times 10^{13}~{}L_{\odot}$ (median
$1.4\times 10^{13}~{}L_{\odot}$). Since we expect the ALMA 1.1 mm dust
continuum to capture the structure of the IR-emitting regions of the source
plane, we use this magnification factor as a proxy in place of $\mu_{\rm IR}$.
Deriving the intrinsic IR SED would ideally involve de-lensing each of the IR
photometric measurements, but the resolution of Planck, WISE, Herschel, and
AzTEC in the mid-IR and far-IR is dramatically coarser than the high-
resolution ALMA 1.1 mm imaging. Using the magnification from a representative
high-resolution rest-frame far-IR image remains the best approach. However,
for the objects modeled in this work and in others summarized in Table 2, four
are lacking comparable ALMA imaging (due primarily to being inaccessible
northern sources). In Figure 3, we include the magnifications measured from 6
GHz radio continuum imaging at similar angular resolution, which are
consistent within 40% relative difference for the majority of cases. These
$\mu_{\rm 6GHz}$ factors can thus be used to roughly estimate intrinsic
$L_{\rm IR}$, as indicated in Table 3. However, because of these possible
differences in magnification, we do not include the objects lacking far-IR
magnifications in the subsequent discussion on source sizes and star formation
rate surface densities.

Despite their high inferred SFRs of $170-6300~{}M_{\odot}~{}{\rm yr}^{-1}$
(median $1500~{}M_{\odot}~{}{\rm yr}^{-1}$; Table 3), it is not clear where
exactly the PASSAGES members fall in relation to the star-formation main
sequence, a tight correlation between SFR and $M_{\star}$, the normalization
of which evolves with redshift (e.g. Brinchmann et al. 2004; Noeske et al.
2007; Elbaz et al. 2007; Daddi et al. 2007; Whitaker et al. 2012; Speagle et
al. 2014). The peak and Rayleigh-Jeans tail of their dust SEDs are well
sampled by Planck, Herschel, AzTEC, and ALMA photometry, leading to secure SFR
estimates; extensive multi-line/continuum modeling by Harrington et al. (2021)
results in robust molecular gas estimates. Yet, the dearth of stellar mass
estimates makes it impossible to reliably place these other properties in
context. This is not a novel problem for DSFGs: it has not been settled
whether they largely represent abnormal outliers forming stars at elevated
rates compared to other objects of similar mass, or if their stellar masses
are proportionally higher to match their extreme SFRs (Casey et al. 2014).
Typically poor constraints on rest-frame optical luminosities due to heavy
dust attenuation make the assessment of $M_{\star}$ difficult for this
population. This is additionally complicated by both the higher sensitivity to
the initial mass function (IMF), as shorter-wavelength optical and UV fluxes
predominantly trace the contribution of the most massive stars, and the less-
constrained contribution of AGN to the near-IR continuum (Hainline et al.
2011). Recent work with other submillimeter-selected samples—primarily with
$\log_{10}[{\rm SFR}/M_{\odot}~{}{\rm yr}^{-1}]\sim 2.5-3.5$, in line with
PASSAGES—by Miettinen et al. 2017 (for $z\sim 1-7$ DSFGs in the COSMOS field;
Scoville et al. 2007) and Barrufet et al. 2020 (for $z\sim 1-5$ DSFGs in the
North Ecliptic Pole field) suggested that around 60% of DSFGs might be
consistent with the main sequence at their respective redshifts, whereas
Ikarashi et al. (2017) suggested 72% were safely above the main sequence
($z\approx 1.4-2.5$).

Harrington et al. (2021) and Berman et al. (2022) highlighted the PASSAGES
sample’s short gas depletion times ($\tau_{\rm dep}=100-400$ Myr), but these
values are not particularly meaningful without an estimate of the average
depletion time of a main sequence galaxy at the same redshift and of similar
stellar mass (e.g. Scoville et al. 2017; Tacconi et al. 2018). Following Liu
et al. (2019a), Wang et al. (2022) recently derived a best-fit relationship
for molecular gas depletion time $\tau_{\rm mol}\equiv M_{\rm mol}/{\rm SFR}$
as a function of stellar mass, redshift/cosmic age, and distance to the main
sequence, $\Delta{\rm MS}\equiv{\rm log}_{10}({\rm SFR}/{\rm SFR}_{\rm MS})$
(see their equation 9). For fiducial stellar masses111111These values are
chosen in line with the characteristic mass (and $1\sigma$ confidence
interval) of Leslie et al. 2020 that contributes most to the cosmic SFR
density at $z=2.5$. Moreover, abundance matching by Behroozi et al. (2013)
suggests an average upper stellar mass limit of $5\times 10^{11}~{}M_{\odot}$
within halos of mass $10^{12}-10^{13}~{}M_{\odot}$. of $10^{10.5}M_{\odot}$,
$10^{10.75}M_{\odot}$, and $10^{11}M_{\odot}$ at a representative redshift of
$z=2.5$, this yields main-sequence depletion times of 340, 360, and 380 Myr,
respectively. Conversely, for our observed range of depletion times—with an
inner 68% confidence interval of $\tau\sim 150-350$ Myr—we obtain $\Delta{\rm
MS}=-0.02-0.69$ for $M_{\star}=10^{10.5}~{}M_{\odot}$, $\Delta{\rm
MS}=0.02-0.78$ for $M_{\star}=10^{10.75}~{}M_{\odot}$, and $\Delta{\rm
MS}=0.08-0.88$ for $M_{\star}=10^{11}~{}M_{\odot}$. This suggests that a clear
majority of the PASSAGES objects should have specific star formation rates
(sSFR $\equiv{\rm SFR}/M_{\star}$) that lie at least above the center of the
main sequence, but raises doubts that they could all be classified uniformly
as traditional starbursts. Specifically, a median depletion time and stellar
mass (250 Myr and $M_{\star}=10^{10.75}~{}M_{\odot}$) yields $\Delta{\rm
MS}=0.32$, which falls below the $\Delta{\rm MS}>0.6$ dex threshold for
starbursts proposed by Rodighiero et al. (2011).

Another independent approach is to use a relation between star-formation rate
surface density $\Sigma_{\rm SFR}$ and $\Delta{\rm log}_{10}({\rm sSFR})_{\rm
MS}\equiv{\rm log}[{\rm sSFR}/{\rm sSFR}_{\rm MS}(M_{\star},z)](=\Delta{\rm
MS})$, as proposed by Jiménez-Andrade et al. (2019). For the highest-redshift
bin121212We note that the next highest redshift bin, $1.3<z\leq 1.75$, follows
effectively the same relation, but where the vertical offset is 1.1 instead of
1.2., $1.75<z\leq 2.25$, Jiménez-Andrade et al. find a best-fit power-law
relation of ${\rm log}(\Sigma_{\rm SFR})=[1.5\pm 0.2]\times\Delta{\rm
log}({\rm sSFR})_{\rm MS}+[1.2\pm 0.1]$ for galaxies with ${\rm
log}_{10}(M_{\star}/M_{\odot})\gtrsim 10.5$. For this work, the inner 68%
confidence interval in star-formation surface density is $\Sigma_{\rm
SFR}\approx 20-160$ (see §4.4), suggesting distances to the main sequence of
$\Delta{\rm log}_{10}({\rm sSFR})_{\rm MS}=0.07-0.67$, in rather close
agreement with the values estimated from depletion times. However, this
relation seems to evolve somewhat from $z=0.5-2$, where the slope of the
relation becomes shallower and the offset increases with redshift, meaning
that surface densities are higher in the Universe, with smaller discrepancies
in $\Sigma_{\rm SFR}$ between starbursts and main sequence galaxies, so its
predictive power may be weakened for higher-$z$ PASSAGES members.

Forthcoming dynamical measurements of the gravitational potentials for
PASSAGES will offer helpful context to this discussion, but at present, their
actual specific star formation rates remain an open question. Sensitive
telescopes like JWST will also offer the opportunity to better constrain the
rest-frame UV and optical SED of such DSFGs, rendering stellar mass
measurements slightly less uncertain.

To reiterate, based on the gravitational lens modeling in this work, we have
found that the PASSAGES sample overwhelmingly consists of objects that are
intrinsically luminous and not extraordinarily amplified by lensing. This
discovery prompts a more thorough exploration of their intrinsic properties,
which we discuss in the next section.

###  4.3 Source-plane reconstructed galaxy sizes

In §3.6, we discuss our computation of galaxy size or effective radius as the
geometric mean of the semi-major and semi-minor axes of a best-fit Gaussian.
At 1 mm, this results in dust continuum sizes of $R_{\rm e}=1.7-4.3$ kpc
(median 3.0 kpc; Table 3). As we discuss in this section, the PASSAGES
galaxies are larger at this wavelength than most DSFGs studied thus far, but
to a degree that is consistent with their larger luminosities. However, as
these are lensed objects selected by their image-plane fluxes, selection
effects may bias the intrinsic sizes to which we are sensitive.

####  4.3.1 Possible source size bias?

Figure 5:  Our model-derived 1.1 mm magnification factors vs. (angular)
effective radii for the PASSAGES members of this work (closely following Fig.
5 of Spilker et al. 2016), in the context of matching measurements at 870
$\mu$m for the Herschel (Bussmann et al. 2015) and SPT (Spilker et al. 2016)
samples. All measurements are made consistently for $\approx
0.5^{\prime\prime}$ ALMA imaging, but Bussmann et al. (2013, 2015) and Spilker
et al. (2016) use a parameterized source-plane in their lens modeling
approaches. Black markers indicate binned values of $R_{e}$ and $\mu$ for the
three samples together (with error bars showing the 16th and 84th percentiles
for each bin). An apparent trend where larger magnifications are correlated
with smaller intrinsic source sizes would indicate the presence of a possible
size bias. While there is not a clear such correlation at all scales, it does
appear to be the case that the highest-magnification objects are
preferentially more compact, as noted by Spilker et al. (2016). With the
addition of our sample, it is also evident that the most extended objects tend
to have lower magnifications; all objects with $R_{e}>0.3^{\prime\prime}$ have
magnifications $\mu\leq 10$.

The intrinsic sizes of the DSFGs that make up the PASSAGES sample have
implications for both the maximum lensing magnifications and the maximum star
formation surface densities that can be observed. For the former, there is a
small region of the source plane capable of producing high magnifications
(e.g. $\mu>10$), especially in the case of galaxy-scale lenses where the
caustic network is comparable in size to that of typical background objects.
For EAGLE-simulated lenses at $z_{\rm lens}\approx 0.37$, Robertson et al.
(2020) found source-plane cross-section solid areas at $z_{s}=2$ satisfying
$|\mu|>10$ to be $\sigma^{\rm S}_{\rm lens}\approx 0.008$ and 0.08 arcsec2 for
halo masses of $10^{12}$ and $10^{13}~{}M_{\odot}$, respectively. For large
magnification thresholds ($\mu_{0}\gg 1$), this cross-section drops off as
$\sigma^{\rm S}_{\rm lens}\propto\mu_{0}^{-2}$ (see also de Freitas et al.
2018).

In the vicinity of a caustic curve (where magnification diverges), increasing
source size has the effect of diluting the overall magnification, as the outer
extent of the source reaches into lower magnification regions (such as the
local minimum inside the caustic curves, or the asymptotically-decreasing
magnifications exterior to the caustics). However, in simulations to study the
predicted size bias of flux-selected lens samples, Hezaveh et al. (2012) found
the intriguing result that intermediate total magnifications ($\mu\sim 10$)
could preferentially magnify diffuse (or extended) components in the source
plane, contrary to the common intuition (e.g. Serjeant 2012) that more compact
components are amplified more strongly (which Hezaveh et al. did indeed find
to be the case for smaller and larger magnifications). This latter effect may
be due to the nuanced consideration that a larger source still has a greater
chance of being modestly magnified ($\mu<10$). In other words, as described by
Robertson et al. (2020), the maximum possible magnification (optimized as a
function of source position) decreases for larger source sizes (as a
decreasing fraction of the source will be some arbitrary distance from the
caustic), but the probability of some part of the source being near to the
caustic network increases for larger sources. In the case of a singular
isothermal sphere acting as a lens, de Freitas et al. (2018) determined that
the cross section for $\mu>10$ was optimally maximized when the source-plane
object had an angular size of $\approx 30\%$ the Einstein radius (in contrast
to $\approx 15\%$ for the threshold $\mu>20$). We make further comparisons to
these predictions in §4.5.

Figure 6:  Relation between far-IR effective radius $R_{e}$ (measured at 1 mm)
and the magnification-corrected, intrinsic IR luminosity (see Table 3) for the
15 PASSAGES objects studied in this work. The gray shaded area shows the
Eddington-limited region (1000 $M_{\odot}$ yr-1 kpc-2) from Andrews & Thompson
(2011). For comparison, we also plot data points from similar samples of
lensed DSFGs studied by Bussmann et al. (2013) and Spilker et al. (2016), with
corrected $L_{\rm IR}$ values from Reuter et al. 2020 for the latter. For
additional comparison, see Figure 5 of Enia et al. 2018 (not shown).  Figure
7:  As with Fig. 6, but with unlensed samples shown for comparison. Objects
are drawn from the ALESS sample (Hodge et al. 2016; Swinbank et al. 2014) and
from Fujimoto et al. (2018), comprising 33 stacked ASAGAO objects (ALMA
twenty-Six Arcmin2 survey of GOODS-S at One-millimeter) at $L_{\rm IR}\approx
10^{12}~{}L_{\odot}$ and the median of 12 other bright ALMA sources ($L_{\rm
IR}\approx 10^{12.6}~{}L_{\odot}$). The shaded red region illustrates the
best-fitting power-law relation found by Fujimoto et al. (2017) in a large
sample of $\approx 500$ archival 1 mm maps from ALMA, $R_{e}\propto L_{\rm
IR}^{0.28\pm 0.07}$. Actual dispersion in the sample is much larger than the
shaded region, and is consistent with that of the PASSAGES objects. For a
secondary comparison, the green shaded region found by Burnham et al. (2021)
(using the $L_{\rm IR}-\lambda_{\rm peak}$ trend from Casey et al. 2018) is
consistent but slightly steeper. The hypothesized Eddington limit from Andrews
& Thompson (2011) is again shown as a gray forbidden region, and star
formation surface densities of 10, 100, and 1000 $M_{\odot}$ yr-1 kpc-2 are
shown as dashed lines. The sample from this work is located at higher
luminosities ($L_{\rm IR}\gtrsim 10^{13}L_{\odot}$), but effective radii are
in general correspondingly larger, and thus consistent with the previous
trends.

This lensing size bias effect has been tested empirically by many works,
including Serjeant (2012), and Spilker et al. (2016). In the case of the
latter, Spilker et al. (2016) showed that 870$\micron$ magnifications for SPT
lenses, along with Herschel lenses from Bussmann et al. 2013 and Bussmann et
al. 2015, did not appear to correlate strongly with source size across the
entire range of observed sizes or magnifications. In fact, the largest sources
covered a wide range of magnifications consistent with the remaining objects.
However, the authors did note that the highest-magnification objects
($\mu_{870\mu m}>10$) were preferentially more compact. Notwithstanding this,
Spilker et al. found that the distribution of source sizes in the SPT and
Herschel lens samples was statistically consistent with that of unlensed
DSFGs, such as the 850 $\micron$ SCUBA-2 Cosmology Legacy Survey (Simpson et
al. 2015a, b; median FWHM $0.30\pm 0.04^{\prime\prime}$) and 1.1-mm AzTEC
sources (Ikarashi et al. 2015, 2017; median FWHM $0.20\pm 0.04^{\prime\prime}$
and $0.31\pm 0.03$). Together, these tests suggested that the hypothesized
size bias was not observable with the then-available sample size.

In Fig. 5, we similarly test for the presence of any size bias for the
PASSAGES sample (extending Fig. 5 of Spilker et al. 2016). In our case, there
is some evidence of anti-correlation between magnification and effective
radius, especially for $R_{e}\gtrsim 0.2^{\prime\prime}$, a region of the
parameter space for which our sample contributes significantly. Objects with
sizes larger than $0.3^{\prime\prime}$ all have lower magnifications, $\mu\leq
10$, whereas smaller objects exhibit a large range of magnifications. These
findings raise the question of the nature of selection biases in the
definition of PASSAGES (Harrington et al. 2016; Berman et al. 2022), as
objects capable of producing observed luminosities $\mu L_{\rm IR}\sim
10^{14}L_{\odot}$ must have some combination of high total lensing
magnification and high intrinsic luminosity. If galaxy-scale lenses were to
preferentially magnify more compact objects, then there would be a negative
power-law slope between $\mu$ and $R_{e}$. In tandem with the positive slope
of $L_{\rm IR}$ vs. $R_{e}$, it might be the case that the selection function
of this sample identifies objects of an optimal size that maximizes the
product of $\mu$ and $L_{\rm IR}$ (e.g. Lutz 2014). Of course, the
distribution of intrinsic luminosities is also subject to the luminosity
function at high-redshift, as there is an effective cutoff on the high-
luminosity end.

Table 3:  Lensing-corrected source-plane effective radii at 1.1 mm and
intrinsic total IR luminosities $L_{\rm IR}$, using $\mu_{\rm 1.1mm}$
magnification factors as representative of dust emission magnification. The
continuum synthesized beam parameters for each field are shown for context.
The star formation rate surface density, $\Sigma_{\rm SFR}$ is given as both a
global average value (total SFR divided by the area of the disk with given
effective radius), and as a peak value given in brackets (see §4.4).

Field | $R_{{\rm e},{1100\mu{\rm m}}}$ | $R_{{\rm e},{1100\mu{\rm m}}}$ | $\theta_{\rm maj}$ | $\theta_{\rm min}$ | PA | ${L_{\rm IR}}^{a}$ | ${{\rm SFR}_{\rm IR}}^{b}$ |  $\Sigma_{\rm SFR}$ [$\Sigma_{\rm SFR,peak}$]
---|---|---|---|---|---|---|---|---
| [′′] | [kpc] | [′′] | [′′] | [∘] | [$10^{13}L_{\odot}$] | [$M_{\odot}~{}{\rm yr}^{-1}$] | [$M_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$]
PJ011646 | $0.45\pm 0.02$ | $3.7\pm 0.2$ | 0.43 | 0.38 | -83.2 | $1.7\pm 0.4$ | $1820\pm 460$ |  $42\pm 11$ [41]
PJ014341 | $0.31\pm 0.02$ | $2.5\pm 0.2$ | 0.49 | 0.40 | 67.7 | $0.16\pm 0.05$ | $170\pm 50$ |  $9\pm 3$ [16]
PJ020941 | $0.29\pm 0.07$ | $2.3\pm 0.6$ | 0.54 | 0.40 | 63.8 | $1.9\pm 0.8$ | $2030\pm 820$ |  $119\pm 56$ [136]
PJ022633 | — | — | — | — | — | $1.0\pm 0.4^{\dagger}$ | $1020\pm 470^{\dagger}$ | —
PJ030510 | $0.44\pm 0.05$ | $3.6\pm 0.4$ | 0.85 | 0.59 | 80.9 | $5.9\pm 1.7$ | $6290\pm 1850$ |  $153\pm 48$ [187]
PJ105353 | $0.22\pm 0.01^{\dagger\dagger}$ | $1.7\pm 0.0$ | 0.45 | 0.27 | 63.1 | $2.5\pm 0.6$ | $2610\pm 590$ |  $290\pm 66$ [264]
PJ112713 | — | — | — | — | — | $0.6\pm 0.1^{\dagger}$ | $590\pm 120^{\dagger}$ | —
PJ113805 | $0.36\pm 0.02$ | $3.0\pm 0.2$ | 1.03 | 0.66 | 3.5 | $0.9\pm 0.2$ | $910\pm 190$ |  $32\pm 7$ [37]
PJ113921 | $0.40\pm 0.07$ | $3.1\pm 0.5$ | 0.82 | 0.69 | -14.0 | $2.5\pm 0.7$ | $2660\pm 750$ |  $87\pm 29$ [79]
PJ132630 | $0.38\pm 0.03$ | $2.9\pm 0.2$ | 0.99 | 0.64 | -2.2 | $1.7\pm 0.4$ | $1860\pm 470$ |  $68\pm 18$ [61]
PJ133634 | — | — | — | — | — | $2.2\pm 0.9^{\dagger}$ | $2310\pm 990^{\dagger}$ | —
PJ144653 | $0.38\pm 0.01$ | $3.1\pm 0.1$ | 0.76 | 0.68 | 11.7 | $0.7\pm 0.2$ | $710\pm 180$ |  $24\pm 6$ [24]
PJ144958 | $0.33\pm 0.10$ | $2.7\pm 0.8$ | 0.81 | 0.67 | 6.3 | $1.0\pm 0.6$ | $1050\pm 580$ |  $45\pm 28$ [78]
PJ160722 | — | — | — | — | — | $0.2\pm 0.1^{\dagger}$ | $220\pm 110$ | —
PJ231356 | $0.52\pm 0.05$ | $4.3\pm 0.4$ | 0.78 | 0.61 | 68.0 | $1.0\pm 0.3$ | $1060\pm 310$ |  $18\pm 6$ [29]
Literature | $R_{{\rm e},{870\mu{\rm m}}}$ | $R_{{\rm e},{870\mu{\rm m}}}$ | $\theta_{\rm maj}$ | $\theta_{\rm min}$ | PA | ${L_{\rm IR}}$ | ${{\rm SFR}_{\rm IR}}^{c}$ | $\Sigma_{\rm SFR}$
| [′′] | [kpc] | [′′] | [′′] | [∘] | [$10^{13}L_{\odot}$] | [$M_{\odot}~{}{\rm yr}^{-1}$] | [$M_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$]
PJ105322d | — | — | — | — | — | $2.9\pm 0.2$ | $3090\pm 210$ | —
PJ112714e | $0.17\pm 0.04$ | $1.4\pm 0.3$ | 0.90 | 0.75 | 39 | $0.4\pm 0.1$ | $430\pm 110$ | $70\pm 20$
PJ132302d | — | — | — | — | — | $0.9\pm 0.1$ | $960\pm 110$ | —
PJ132934f | — | — | — | — | — | $1.3\pm 0.1$ | $1380\pm 110$ | —
PJ142823g | $0.08\pm 0.05$ | $0.71\pm 0.43$ | $\approx 0.7$ | $\approx 0.7$ | — | $0.6\pm 0.3$ | $640\pm 320$ | $400\pm 320$
PJ154432d | — | — | — | — | — | $0.31\pm 0.02$ | $330\pm 20$ | —
PJ160918d | — | — | — | — | — | $1.6\pm 0.1$ | $1710\pm 130$ | —

Note. — a Lensing-corrected intrinsic luminosities derived by dividing
apparent far-infrared luminosities ($8-1000\micron$) from Table 6 of Berman et
al. (2022) (or Table 4 of Harrington et al. 2016) by magnification factors in
Table 2. b Lensing-corrected IR-inferred star formation rates derived by
dividing Starburst SED star formation rates from Table 6 of Berman et al. 2022
(or Table 4 of Harrington et al. 2016) by magnification factors in Table 2. c
Star formation rates derived from Kennicutt 1998 calibration, ${\rm
SFR}=L_{\rm IR}/(9.4\times 10^{9})M_{\odot}{\rm yr}^{-1}$. d De-lensed $L_{\rm
IR}$ from Cañameras et al. (2015, 2018a). e De-lensed $L_{\rm IR}$ from
Cañameras et al. (2015, 2018a, 2018b), using latest magnification factor from
Pascale et al. (2022). f De-lensed $L_{\rm IR}$ from Díaz-Sánchez et al.
(2017). g De-lensed $L_{\rm IR}$ from Bussmann et al. (2013). † Approximate
intrinsic IR luminosity estimated using 6 GHz magnification factor (as high-
resolution 1 mm/870 $\mu$m image is not available, see discussion in §4.2). ††
Low-resolution ALMA Band 7 ($870\mu$m) measurement from program 2019.1.01636.S
(PI: M. Yun).

####  4.3.2 Intrinsic far-IR source size

With this work, we seek to place the PASSAGES objects modeled in this work in
the context of other submm-selected lensed systems, as seen in Figures 5 and
6. In comparison to the Herschel and SPT lensed DSFGs, the PASSAGES DSFGs
encompass larger intrinsic source sizes in physical units. However, the
largest PASSAGES object is comparable in size to the largest member of the
Bussmann et al. (2013) sample (J091305.0-005343 or SDP.130, $R_{e}=4.14\pm
0.72$). Yet, despite occupying the upper end of luminosities covered by the
other samples, PASSAGES objects show notable evidence for lower star-formation
rate surface densities, which are primarily between $10-100~{}M_{\odot}$ yr-1
kpc-2, consistent with the vast majority of over 1000 dusty starbursts studied
by Fujimoto et al. (2017).

To further elucidate this discrepancy, we extrapolate our analysis to include
unlensed DSFGs. Figure 7 reveals that PASSAGES is in accordance with a scaling
relation of $L_{\rm IR}-R_{e}$ found by Fujimoto et al. (2017), but both
larger and more luminous than the physical quantities derived by Bussmann et
al. 2013 and Spilker et al. 2016 (in members for which redshifts are
measured). Berman et al. (2022) found that the PASSAGES DSFGs account for an
average number density of about $10^{-2}~{}{\rm deg}^{-2}$, in comparison with
the Herschel samples’ density of $0.02-0.1~{}{\rm deg}^{-2}$ (Vieira et al.
2013; Wardlow et al. 2013; Weiss et al. 2013), suggesting that PASSAGES
includes some of the very rarest all-sky objects. It is thus not entirely
surprising that they would have higher intrinsic luminosities (when
considering comparable lensing magnifications) and larger source sizes than
other DSFG samples, which probe deeper into the regime of less extreme DSFGs
but are restricted to a smaller survey area.

It might also be the case that other factors are at play. In the case of
Spilker et al. (2016) in particular, the more compact sizes might be due to
their higher redshift ($z\sim 3.5-5$ primarily). While the angular size does
not evolve significantly for a given physical size between $z\sim 4$ and
$z\sim 2$, Fujimoto et al. (2017) find a detectable evolution in physical size
for a given $L_{\rm IR}$ ($L_{\rm IR}\sim 10^{12}-10^{13}~{}L_{\odot}$),
albeit with a small sample size ($<10$ galaxies at $z>4$). On the other hand,
Enia et al. (2018) do not observe a clear redshift evolution, and note that
some apparent evolution could be due to selection effects biased against the
decrease in surface brightness (for fixed luminosity and size) at higher
redshifts. Additionally, the lens modeling approach for PASSAGES in this work
is fundamentally different from Bussmann et al. (2013, 2015), and Spilker et
al. (2016), which all optimize model parameters directly over interferometric
visibilities. This requires an adoption of a parametric form of the source-
plane distribution (in contrast to our method, described in Section 3.6). Enia
et al. (2018) note that estimating source size from parametric modeling
consistently resulted in systematically lower values, in comparison with an
approach of measuring the area of the source above a certain signal-to-noise
threshold. This is especially true for clumpy morphologies, where a Gaussian
or Sérsic profile might underestimate the true spatial extent.

An ample number of recent studies have revealed the presence of significant
sub-kpc inhomogeneous, clumpy substructure that can dominate the sub-mm flux
of DSFGs (Swinbank et al. 2010, 2011; Danielson et al. 2011; Hodge et al.
2012, 2019; Alaghband-Zadeh et al. 2012; Menéndez-Delmestre et al. 2013;
Hatsukade et al. 2015; Iono et al. 2016; Cañameras et al. 2018b; Tadaki et al.
2018; Dessauges-Zavadsky et al. 2019; Ramasawmy et al. 2021; Spilker et al.
2022), although some counter-evidence to the ubiquity of clumps in DSFGs does
exist (e.g. Hodge et al. 2016; Ivison et al. 2020). These clumpy sites of star
formation appear to be somewhat akin to larger ($\sim$100pc), more luminous
analogs to giant molecular clouds (Swinbank et al. 2010; Dessauges-Zavadsky et
al. 2019). Similar substructure is also apparent in rest-frame UV/optical
imaging of star-forming galaxies (e.g. Förster Schreiber et al. 2011; Iani et
al. 2021, and references therein). While the angular resolution of this
current ALMA imaging ($\sim 0.5^{\prime\prime}$) is consistent with that of
the Herschel and SPT lenses, it is challenging at present to assess the degree
of clumpiness at sub-kpc scales (Fig. 2). In the future, we intend to compare
these results with higher-resolution ALMA observations, which enables physical
resolutions closer to 100 pc with lensing, at which point we can adequately
map the spatial variation of $\Sigma_{\rm SFR}$ (e.g. Hodge et al. 2019).

The intrinsic far-IR continuum sizes in this study span a range of $1.7-4.3$
kpc, consistent with the typical predicted size of $\sim 3$ kpc from radiative
transfer modeling by Harrington et al. (2021), especially when adopting a
revised fiducial magnification of $\mu\sim 10-20$. For optically-thick
regions, the Stefan-Boltzmann Law applied to typical ULIRG luminosities and
temperatures gives $R_{e}=1.09\cdot[L_{\rm
IR}/10^{12}~{}L_{\odot}]^{1/2}\cdot[T_{\rm dust}/35~{}{\rm K}]^{-2}$ kpc
(Scoville 2013), so the luminosities for this sample ($L_{\rm
IR}\approx(2-60)\times 10^{12}L_{\odot}$) translate to expected sizes of
$R_{e}\approx 1.5-8.5$ kpc. The far-IR physical sizes we find are
comparatively a bit larger than those found by recent studies of
submillimeter-bright galaxies, including Simpson et al. 2015a (median
$R_{e}=1.2\pm 0.1$ kpc, interquartile range $1.8-3.2$ kpc), Tadaki et al. 2015
($R_{e}=0.7\pm 0.1$ kpc), Barro et al. 2016a and Oteo et al. 2016 ($R_{e}\sim
1$ kpc), Hodge et al. 2016 (median $R_{e}=1.8\pm 0.2$ kpc), Rujopakarn et al.
2016 (median $R_{e}=2.1\pm 0.9$ kpc), Ikarashi et al. 2017 ($R_{e}=1.6\pm 0.3$
kpc), Tadaki et al. 2017 ($R_{e}\sim 1.5$ kpc), and Gómez-Guijarro et al. 2022
(median $R_{e}=0.73\pm 0.29$ kpc).

As emphasized by radiative transfer modeling of FIRE-2 (Hopkins et al. 2014,
2018) zoom-in cosmological simulations by Cochrane et al. (2019), the spatial
extent of single-band dust continuum emission is determined by a combination
of the spatial extent of dust mass and of recent star formation, given the
mechanism of dust heating. Cochrane et al. confirmed that spatially-resolved
dust temperature gradients (measured via spectral indices of multiple sub-mm
bands) is an effective method to disentangle these combined effects (e.g. as
done with local galaxies by Galametz et al. 2012).

A consistent picture has arisen that the far-IR emission of high-$z$ star-
forming galaxies originates from a substantially more compact spatial region
(by a typical factor of $2-4$) than that traced by rest-frame UV and optical
imaging (e.g. Calanog et al. 2014; Simpson et al. 2015a; Chen et al. 2015;
Ikarashi et al. 2015; Barro et al. 2016b; Talia et al. 2018; Tadaki et al.
2020; Pantoni et al. 2021). This is likewise generally the case for radio
sizes as well, which are comparable to far-IR sizes (e.g. Bondi et al. 2018;
Jiménez-Andrade et al. 2019). This compact phase of star formation is
consistent with the evolutionary schemes for progenitors of massive
elliptical/quiescent galaxies (e.g. Barro et al. 2013; Toft et al. 2014; Lapi
et al. 2018). Intriguingly, Ikarashi et al. (2017) found that DSFGs with
evidence for composite star-forming/AGN growth were more compact than star
formation-dominated or AGN-dominated objects ($R_{e}=1.0\pm 0.2$ kpc vs.
$R_{e}=1.6\pm 0.3$ for star-formation dominated or $R_{e}=1.5\pm 0.6$ for AGN-
dominated), which the authors suggest might be due to supermassive black hole
growth during a compact star-forming phase of merger coalescence (Springel et
al. 2005). In contrast, extended FIR sizes of star-formation-dominated DSFGs
may arise from an intermediary merger stage, while extended morphology of AGN-
dominated DSFGs may be the result of positive feedback inducing star formation
at larger radii (Ishibashi & Fabian 2012). At present, though, there is not
overwhelming evidence to support or rule out this theory. As an added caveat,
the extended rest-frame UV/optical sizes could be the result of a radial
dependence of dust attenuation (e.g. Nelson et al. 2016), where a strongly
centrally-concentrated distribution could result in artificially-extended
half-light radii larger than half-mass radii (e.g., Tacchella et al. 2018;
Bondi et al. 2018; Suess et al. 2019). Early studies with JWST have begun to
test this hypothesis, and initial results (e.g. Suess et al. 2022) indicate
the important role that spatial variations in dust attenuation and stellar
population ages can play.

The elevated dust-emitting sizes of PASSAGES starbursts may be explained
primarily by the larger intrinsic IR luminosities. Fujimoto et al. (2017)
uncovered a power-law scaling relation between FIR effective radius and
luminosity, $R_{e}\propto L_{\rm FIR}^{0.28\pm 0.07}$ for $L_{\rm FIR}\gtrsim
10^{12}L_{\odot}$, in studying a large sample of several hundred 1 mm ALMA
detections, with effective radii spanning physical scales of $0.2-5.0$ kpc.
Burnham et al. (2021) offered some confirmation with a comparable power-law
slope of $R_{e}\propto L_{\rm FIR}^{0.37\pm 0.03}$ (in a much smaller sample
of 18 DSFGs, median $R_{e}=0.32\pm 0.09$ kpc). However, these trends have yet
to be tested for a large ensemble of objects at $L_{\rm IR}\sim
10^{13}L_{\odot}$, i.e. HyLIRGs. As we demonstrate in Fig. 7, the PASSAGES
objects in this work are consistent with the previously-discovered trends.
Intriguingly, in comparison with the Herschel and SPT lensed DSFGs of Bussmann
et al. (2013) and Spilker et al. 2016, as shown in Fig. 6, the PASSAGES
objects appear to be larger on average for a given IR luminosity. This may be
a result of the lower median redshift of PASSAGES, as discussed in the next
section.

Some works have applied the Stefan-Boltzmann law at galaxy scales, where the
IR luminosity can be related to effective galaxy size and dust temperature:
$L_{\rm IR}\approx 4\pi R_{e}^{2}\sigma T_{d}^{4}$, for constant
$\sigma=5.670\times 10^{-5}$ erg s-1 cm-2 K-4. Yan & Ma (2016) suggested the
use of a generalized, temperature-dependent $\tilde{\sigma}(T)=\sigma\times
10^{-3}(-3.03T^{1.5}+45.55T-127.53)$ in the case of a modified blackbody
spectrum, whereas Ma & Yan (2015) used a simpler power-law relation to arrive
at $L_{\rm IR}=4\pi R_{e}^{2}\sigma T_{d}^{4.32}$. Here we adopt
$\tilde{\sigma}(T)\propto T^{0.8}$, which closely approximates the Yan & Ma
relation as a single power-law for the primary range of interest,
$T_{d}=20-50$ K, such that $L_{\rm IR}\propto R_{e}^{2}T_{d}^{4.8}$.
Empirically, Fujimoto et al. (2017) and Burnham et al. (2021) found a relation
of $L_{\rm IR}\propto R_{e}^{2.7-3.6}$, with which the PASSAGES objects concur
(Fig. 7). In keeping with the Stefan-Boltzmann relation, this would
necessitate that $T_{d}\sim R_{e}^{0.2-0.3}$, indicating a weak dependence,
but this is contrary to some expectations that dust temperature decreases with
size. For example, for 16 objects studied by Hodge et al. 2016, there is an
approximate correlation of $\log[T_{d}]/\log[R_{e}]\approx-2$ to $-1$
(although there is also an anti-correlation of temperature with redshift,
which may play a larger role). Regardless, it may be the case that this simple
framework does not adequately describe the infrared emission of dusty star-
forming objects, as it assumes that they radiate as optically-thick, spherical
blackbodies. For this work in particular, sizes are measured at observed-frame
1 mm, in the Rayleigh-Jeans regime often assumed to be optically thin131313
Harrington et al. (2021) tested this assumption with masses derived from large
velocity gradient (LVG; Goldreich & Kwan 1974; Scoville & Solomon 1974)
radiative transfer models, finding that the single-band estimation can
systematically over-predict mass (see also Förster Schreiber et al. 2018; Liu
et al. 2019b). (which fortuitously allows for the direct estimation of ISM
mass from flux; Scoville et al. 2016, 2017).

Table 4:  Properties of JVLA 6 GHz continuum imaging, including natural-weighting flux ($S_{6{\rm GHz}}$, uncorrected for lensing), inferred 1.4 GHz luminosity, and the derived radio-FIR correlation parameter, $q_{\rm FIR}$. Also included are the lensing-corrected source-plane effective radii (as with Table 3), and image-plane synthesized beam parameters for the image weightings used in the size calculation (and as presented in Fig. 2).  Field | $S_{6{\rm GHz}}$ | $\log[L_{\rm 1.4GHz}]$ | $q_{\rm FIR}$ | $R_{{\rm e},{6GHz}}$ | $R_{{\rm e},{6GHz}}$ | $\theta_{\rm maj}$ | $\theta_{\rm min}$ | PA
---|---|---|---|---|---|---|---|---
| [$\mu$Jy] | $[{\rm W~{}Hz}^{-1}]$ |  | [′′] | [kpc] | [′′] | [′′] | [∘]
PJ011646 | $664\pm 37$ | 24.9 $\pm$ 0.1 | 2.3 $\pm$ 0.2 | $0.67\pm 0.01$ | $5.5\pm 0.1$ | 0.73 | 0.36 | 15.9
PJ014341 | $76\pm 7$ | 23.4 $\pm$ 0.3 | 2.8 $\pm$ 0.4 | $0.32\pm 0.05$ | $2.6\pm 0.4$ | 0.38 | 0.29 | 29.9
PJ020941 | $909\pm 37$ | 25.0 $\pm$ 0.3 | 2.3 $\pm$ 0.5 | $0.23\pm 0.05$ | $1.9\pm 0.4$ | 0.51 | 0.40 | 25.0
PJ022633 | $1295\pm 59$ | 25.0 $\pm$ 0.4 | 2.0 $\pm$ 0.6 | $0.57\pm 0.39$ | $4.3\pm 3.0$ | 0.48 | 0.37 | 64.3
PJ030510 | $221\pm 22$ | 24.9 $\pm$ 0.4 | 2.9 $\pm$ 0.5 | $0.79\pm 0.03$ | $6.5\pm 0.2$ | 1.21 | 0.43 | 27.2
PJ105353 | $874\pm 46$ | 25.3 $\pm$ 0.1 | 2.1 $\pm$ 0.3 | $0.25\pm 0.01$ | $1.9\pm 0.1$ | 0.42 | 0.36 | 24.6
PJ112713 | $678\pm 35$ | 24.6 $\pm$ 0.1 | 2.2 $\pm$ 0.2 | $0.27\pm 0.01$ | $2.3\pm 0.1$ | 0.61 | 0.37 | 74.6
PJ113805 | $525\pm 27$ | 25.0 $\pm$ 0.3 | 2.0 $\pm$ 0.3 | $0.21\pm 0.03$ | $1.7\pm 0.2$ | 0.33 | 0.30 | -67.7
PJ113921 | $249\pm 14$ | 24.7 $\pm$ 0.3 | 2.7 $\pm$ 0.4 | $0.53\pm 0.05$ | $4.1\pm 0.4$ | 0.44 | 0.39 | 47.2
PJ132630 | $92\pm 6$ | 24.2 $\pm$ 0.3 | 3.1 $\pm$ 0.4 | $0.34\pm 0.11$ | $2.7\pm 0.9$ | 0.43 | 0.38 | 80.9
PJ133634 | $403\pm 24$ | 25.0 $\pm$ 0.4 | 2.3 $\pm$ 0.6 | $0.44\pm 0.21$ | $3.3\pm 1.5$ | 0.33 | 0.28 | -69.2
PJ144653 | $901\pm 49$ | 24.3 $\pm$ 0.2 | 2.5 $\pm$ 0.4 | $0.27\pm 0.03$ | $2.2\pm 0.2$ | 0.41 | 0.37 | 49.4
PJ144958 | $400\pm 16$ | 24.4 $\pm$ 0.8 | 2.6 $\pm$ 1.0 | $0.30\pm 0.18$ | $2.5\pm 1.5$ | 0.41 | 0.38 | 61.7
PJ160722 | $162\pm 9$ | 23.9 $\pm$ 0.3 | 2.4 $\pm$ 0.6 | $0.34\pm 0.05$ | $2.9\pm 0.4$ | 0.47 | 0.37 | -1.1
PJ231356 | $3244\pm 163^{\dagger}$ | —† | —† | —† | —† | 0.43 | 0.29 | 44.2
Figure 8:  The far-IR-radio correlation parameter $q_{\rm FIR}$ as a function
of redshift for the PASSAGES sample. The 2$\sigma$ range of $q_{\rm FIR}$ for
Herschel-selected lensed DSFGs (i.e., also selected in the IR) from Giulietti
et al. (2022) are shown as a shaded solid gray region. The median value from
the IR-selected local sample from Yun et al. (2001) and from the radio-
selected sample of Ivison et al. (2010) are shown as a blue dashed and purple
dotted line, respectively. Best-fit redshift evolutions derived from joint
radio- and IR-selected samples from Magnelli et al. (2015) and Delhaize et al.
(2017) are shown as red and green lines, respectively, with shaded regions
showing 1$\sigma$ uncertainties. The gray hatch-shaded region at $q_{\rm FIR}$
is the proposed threshold by Condon et al. (2002) for radio-loud galaxies
powered by AGN, rather than star formation.

####  4.3.3 Radio-FIR correlation and intrinsic radio source size

The well-studied correlation between radio and FIR flux (e.g., Helou et al.
1985; Condon 1992; Yun et al. 2001; Bell 2003; Murphy et al. 2006a, b) is
interpreted to be the result of the far-IR probing heated dust surrounding
star-forming regions and the radio capturing synchrotron emission from
relativistic electrons originating from supernova remnants. As the latter are
the end-product of massive star formation, these two proxies are tightly
correlated when averaged over galactic scales. The correlation even holds
seemingly at sub-galactic scales (e.g. Tabatabaei et al. 2007; Dumas et al.
2011), which may be explained by the more efficient propagation of
relativistic cosmic rays in the dense ISM, with its higher magnetic field
density. The correlation parameter $q_{\rm FIR}$ is computed as

$\displaystyle q_{\rm FIR}=$ $\displaystyle\log[L_{\rm FIR}/(3.75\times
10^{12}~{}{\rm W})]$ $\displaystyle-\log[L_{\rm 1.4~{}GHz}/(1~{}{\rm
W~{}Hz}^{-1})]$ (5)

where we determine $L_{\rm 1.4~{}GHz}$ from 6 GHz flux as

$\displaystyle\frac{L_{\rm 1.4~{}GHz}}{{\rm W~{}Hz}^{-1}}=$
$\displaystyle\frac{1}{\mu_{\rm
radio}}\frac{4\pi}{(1+z)^{1-\alpha}}\bigg{(}\frac{D_{L}(z)}{1~{}{\rm
m}}\bigg{)}^{2}$ (6) $\displaystyle\times\bigg{(}\frac{S_{\rm
6GHz}}{10^{32}~{}\mu{\rm Jy}}\bigg{)}\bigg{(}\frac{\rm 1.4~{}GHz}{\rm
6~{}GHz}\bigg{)}^{-\alpha}$ (7)

for luminosity distance $D_{L}$ in meters and spectral index $\alpha$ (assumed
here to be 0.8; Condon 1992) such that $S_{\nu}\propto\nu^{-\alpha}$. These
luminosities (corrected for lensing magnification $\mu_{\rm radio}$) are
included in Table 4, along with the corresponding values of $q_{\rm FIR}$. We
find values ranging from $q_{\rm FIR}=2.0-3.1$, with a median of 2.3 and
$1\sigma$ dispersion of 0.3, in agreement with the mean found by Yun et al.
(2001) for local IR-selected galaxies, $q_{\rm FIR}=2.34$. This is also in
good agreement with more recent results for Herschel-selected lenses from
Giulietti et al. (2022), who found a $2\sigma$ dispersion of $q_{\rm
FIR}\approx 1.9-2.9$. Unsurprisingly, none of the star-formation-dominated
PASSAGES galaxies fall below the $q_{\rm FIR}<1.8$ threshold for AGN-powered
radio sources proposed by Condon et al. (2002).

As shown in Fig. 8, there is not an evident decline in the $q_{\rm FIR}$
parameter with redshift for PASSAGES galaxies. As discussed extensively by
Sargent et al. (2010a), any observed trend with redshift is heavily dependent
on selection biases from objects selected only in the IR or radio. Accounting
for this affect, Sargent et al. found no evolution in the correlation out to
at least $z\sim 1.4$ (see also Sargent et al. 2010b). However, Magnelli et al.
(2015) and Delhaize et al. (2017), both using a joint radio/FIR selection,
observed similar redshift evolution in $q_{\rm FIR}$ to smaller values.
Moreover, for Herschel-selected lensed DSFGs, Giulietti et al. (2022) remarked
upon a weak but detectable trend in line with previous results. Delhaize et
al. (2017) consider the possibility that contributions of AGN in the radio
regime alone for star-forming galaxies could result in a steepening of the
evolution with redshift. Alternatively, it is possible the calculation of
$q_{\rm FIR}$ for PASSAGES galaxies is weakened by the assumption of a single
radio spectral index $\alpha$, rather than direct measurement for each galaxy.
As also noted by Delhaize et al. (2017), for these higher-frequency
measurements at 6 GHz, free-free emission (following instead
$S_{\nu}\propto\nu^{-0.1}$) may contribute non-trivially at rest-frame
$\nu\gtrsim 20$ GHz (Condon 1992).

Under a different assumption of a flatter spectral index $\alpha=0.5$ (similar
to that of extreme local starbursts, e.g., Condon et al. 1991; Clemens et al.
2008), we would find values of $q_{\rm FIR}$ greater than those reported in
Table 4 by a small amount, $q_{{\rm FIR},\alpha=0.5}-q_{{\rm
FIR},\alpha=0.8}=0.3-0.4$. This change in $q_{\rm FIR}$ is not insignificant,
but it is smaller than the dispersion we find in $q_{\rm FIR}$ for our sample,
and not much larger than the uncertainty on $q_{\rm FIR}$ for any individual
object, so we conclude that the possible effect on our interpretation is
minimal.

Figure 9:  Top: The ratio of 1 mm to 6 GHz continuum size, for the 10 objects
in this work with both measurements available, vs. magnification-corrected
infrared luminosity (blue circles, bottom axis) and vs. star formation rate
surface density $\Sigma_{\rm SFR}$ (red diamonds, top axis). Vertical
errorbars are not shown for the latter for legibility, but they are equal to
those for the blue points. The gray shaded region shows delineates a factor of
$>2$ difference in size. A possible trend where the more luminous, more
densely star-forming objects have a preferentially more compact dust emission
region is detectable, but small sample size makes interpretation difficult.
Bottom: The same ratio of sizes is shown relative to the far-IR-radio
correlation parameter $q_{\rm FIR}$. The radio excess (AGN-powered) regime
from Fig. 8 is again shown as a hatch-shaded region.

In Figs. 9 and 10, we examine the distribution of effective radii for our
sample in both the rest-frame far-IR and radio, in comparison with far-IR
sizes from Bussmann et al. (2013) and Spilker et al. (2016). As both are
expected to delineate the extent of ongoing massive star formation, the radio
and far-IR sizes are reasonably consistent (where both are measurable), with
discrepancies possibly explained partly by different observing resolutions.
While not for the same set of galaxies, Murphy et al. 2017 found 10 GHz radio
sizes comparable to the median dust-emission sizes of DSFGs from Ikarashi et
al. (2015) and Simpson et al. (2015a); see also Jiménez-Andrade et al. (2019)
at 3 GHz. As a caution, however, the 1.4 GHz radio sizes for the same DSFGs
from Simpson et al. 2015a (presented in Biggs & Ivison 2008) are typically
twice as large as the 870 $\mu$m sizes. This comparison of continuum size at
different wavelengths could instead be the result of different optical depths,
or even different surface brightness sensitivities between the two sets of
observations. Simpson et al. (2015a) suggested that the larger diffusion scale
for cosmic rays (of order $\sim 1-2$ kpc; Bicay & Helou 1990; Marsh & Helou
1998; Murphy et al. 2006b) might be responsible, relative to the smaller $\sim
100$ pc diffusion length for the far-IR photons (which also explains why the
radio-FIR correlation begins to decouple at resolved sub-kpc scales). Fig. 9
shows the ratio of 1 mm to 6 GHz effective radii vs. $L_{\rm IR}$ and
$\Sigma_{\rm SFR}$ for the 10 objects with both radii measured, and there is a
discernible trend where the more luminous, more densely star-forming galaxies
have preferentially more compact dust emission regions (or conversely, perhaps
more extended radio emission). When excluding the lower-luminosity object
PJ014341, this negative correlation with luminosity is significant at
$p\lesssim 0.01$. Murphy et al. (2006b) suggested that galaxies with larger IR
surface densities had undergone a more recent episode of star formation, such
that young cosmic rays would only have had time to travel $\sim 100$ pc,
leading to smaller radio scale lengths. Curiously, this appears not to be the
case for our sample, which instead shows larger radio scales with larger
$\Sigma_{\rm SFR}$, but it’s not currently practical to draw conclusions based
on this sample. However, in local edge-on galaxies, Wiegert et al. (2015) also
found preliminary evidence for a correlation of radio halo size with
$\Sigma_{\rm SFR}$, such that a compact star formation distribution was
advantageous to creating radio halos (also Dahlem et al. 2006). Yet, Heesen et
al. (2018) did not find a strong correlation between radio scale height and
SFR, $\Sigma_{\rm SFR}$, or galaxy rotation speed, for both diffusion- vs.
advection-dominated (or wind-driven) modes of cosmic ray transport. Lastly, in
Fig. 9, we examine the same size ratio as a function of the radio-FIR
correlation parameter, $q_{\rm FIR}$. There is again some weak indication of a
negative correlation, but not with sufficient statistical significance given
uncertainties. Such a correlation might reveal the influence of an AGN in
driving towards more compact radio half-light sizes.

Figure 10:  Effective radius vs. redshift for the PASSAGES sample (Tables 3
and 4) and the other lensed DSFGs also shown in Fig. 6 (Herschel, Bussmann et
al. 2013; SPT, Spilker et al. 2016). There is some evidence for a downwards
trend in galaxy size with increasing redshift, but this is illustrated mostly
by the higher-$z$ SPT objects (likely owing to their selection at a longer
wavelength of 1.4mm; Weiss et al. 2013). We compare with the FIR sizes of
stacked main sequence galaxies measured by Wang et al. (2022), in stellar mass
bins of high ($11.0\leq{\rm log}_{10}M_{\star}<12.0$), mid ($10.5\leq{\rm
log}_{10}M_{\star}<11.0$), and low ($10.0\leq{\rm log}_{10}M_{\star}<10.5$).
For $0.4<z<3.6$, Wang et al. find minimal size evolution with redshift, in
accordance with the lack of a strong evolution of radio continuum sizes found
by Jiménez-Andrade et al. (2019).

####  4.3.4 Redshift evolution of galaxy size

Fujimoto et al. (2017) found modest evolution in the $R_{e}-L_{\rm FIR}$
relationship between $z=1-2$ and $z=2-4$, and a more significant evolution (by
a factor of $\sim 2$) to $z=4-6$. This decrease in effective radius with
redshift for a given luminosity is consistent with rest-frame UV/optical
samples. For example, Shibuya et al. (2015) found that the UV/optical
effective radius of UV-luminous star-forming galaxies evolves with redshift as
$R_{e,{\rm UV}}/{\rm kpc}\propto(1+z)^{-0.84\pm 0.11}$. While the slope
between $R_{e}$ and UV luminosity is consistent with redshift, the
normalization decreases for the range $z=0-8$. At longer wavelengths, Jiménez-
Andrade et al. (2021) found at 3 GHz that $R_{e,{\rm radio}}/{\rm
kpc}\propto(1+z)^{-0.3\pm 0.3}$ from $z\approx 0.3-3$, while Lindroos et al.
2018 found steeper evolution at 1.4 GHz of $R_{e,{\rm radio}}/{\rm
kpc}\propto(1+z)^{-1.7}$ out to $z\approx 3$. These results suggest that the
correlation between luminosity and effective radius holds across cosmic time,
but the size of an object with fixed luminosity decreases with lookback time.
In turn, the surface density of star formation appears to increase with
lookback time.

In Fig. 10, we do not observe a clear indication of redshift evolution for the
PASSAGES sample alone, covering only $z=1-3.5$, as is also the case for the
Herschel-selected objects at a similar redshift range (Bussmann et al. 2013).
While the average far-IR size is elevated for PASSAGES relative to the
Bussmann et al. sample, there is significant overlap in the observed size
range. This minimal evolution in far-IR size matches the results found by Wang
et al. (2022) for stacked main-sequence galaxies from $z=0.4-3.6$, in three
stellar mass bins: high ($11.0\leq{\rm log}_{10}M_{\star}<12.0$), mid
($10.5\leq{\rm log}_{10}M_{\star}<11.0$), and low ($10.0\leq{\rm
log}_{10}M_{\star}<10.5$). We find also that these main-sequence sizes are
typically smaller than those found for the DSFGs in this work, but not much
further than the level of $1\sigma$ uncertainties of each bin. On the other
hand, for the higher-redshift SPT sample ($z=3-6$), there is hardly any
overlap in effective radius with PASSAGES galaxies, the latter of which are
systematically larger in size. This systematic discrepancy in sizes above vs.
below $z\sim 3$ be a signpost of redshift evolution in DSFGs, but this direct
comparison is complicated slightly by the different selection effects at play
for each sample.

###  4.4 Eddington-limited star-formation surface densities?

As the DSFGs in this sample are among the most strongly star-forming systems
presently known (Berman et al. 2022), they offer key insight into effective
maximum thresholds of star formation (e.g. Elmegreen 1999; Tacconi et al.
2006). In particular, stellar mass growth is understood to be a self-
regulating process: as gas collapses under self-gravity to form stars, the
short-lived, massive stars inject a sizable radiation pressure in opposition
to dust grains, which are coupled to the remaining gas (Thompson et al. 2005).
The concept that radiation pressure can have significant influence on star
formation is not novel (e.g. Elmegreen 1983; Scoville et al. 2001; Scoville
2003); it can even induce star formation in other regions of a cloud (e.g.
Elmegreen & Lada 1977). Scoville et al. (2001) derived that a luminosity-to-
mass ratio of a star cluster of $\sim 500~{}L_{\odot}/M_{\odot}$ would provide
a radiation pressure in excess of its self-gravity and halt accretion to the
cloud core—analogous to the so-called Eddington limit—and by extension, on
galaxy-wide scales as well. This is a conservative theoretical limit, assuming
only free-fall collapse and no other turbulence- or magnetic-driven obstacles
to molecular core accretion. It also does not consider the likely possibility
of asymmetric cloud geometries reducing the importance of radiation pressure.
Additionally, this limit does not account for mechanical feedback contributed
by stellar winds, which may well be an important factor (e.g. Tan & McKee
2001; Harper-Clark & Murray 2009; Rogers & Pittard 2013).

The vast majority of DSFGs and luminous IR galaxies appear to form stars at
sub-Eddington rates (Tacconi et al. 2006; Hodge et al. 2015, 2019) with the
exception of some local LIRGs and ULIRGs (e.g. Barcos-Muñoz et al. 2017). One
example is the prototypical Arp 220 (Scoville et al. 1997; Downes & Solomon
1998), which has star formation surface densities—defined as $\Sigma_{\rm
SFR}={\rm SFR}/(2\pi{R_{\rm eff,dust}}^{2})$—on the order of
$10^{4}~{}{M}_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$, likely the highest-
known value (Barcos-Muñoz et al. 2015). Here, $\Sigma_{\rm SFR}\approx
1000~{}{M}_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$, or $\Sigma_{\rm
IR}\approx 10^{13}~{}L_{\odot}~{}{\rm kpc}^{-2}$ (Thompson et al. 2005;
Andrews & Thompson 2011), is the typical Eddington limit for a star-forming
system in equilibrium. There is some evidence for near-Eddington star
formation for a small number of high-$z$ objects, including the submillimeter-
bright galaxies GN20 and AzTEC-1 (Younger et al. 2008; Daddi et al. 2009);
AzTEC2-A (Jiménez-Andrade et al. 2020); AzTEC-3 (Riechers et al. 2014); HFLS3
(Riechers et al. 2013); the Cosmic Eyelash (Thomson et al. 2015); SGP38326
(Oteo et al. 2016); SPT0346-52 (Ma et al. 2016); ALMACAL-1 and ALMACAL-2 (Oteo
et al. 2017); ADFS-27 (Riechers et al. 2017); and, more recently, the Sunrise
Arc (Vanzella et al. 2019; Welch et al. 2022). Some quasar hosts have also
shown evidence for such compact, maximal starbursts, at $z\sim 2$ (Stacey et
al. 2021) and $z>6$ (Walter et al. 2009). However, angular resolution plays a
critical role in this measurement, as it may be possible for objects to be
sub-Eddington on global, galaxy-integrated scales but super-Eddington in local
regions of dense star formation, as suggested by Simpson et al. (2015a) and
Barcos-Muñoz et al. (2017). For example, for local ULIRGs, Song et al. (2022)
found radio continuum clumps smaller than 100 pc with large surface densities
up to $\Sigma_{\rm SFR}\approx 1600~{}{M}_{\odot}~{}{\rm yr}^{-1}~{}{\rm
kpc}^{-2}$.

Yet, a vexing result has been the discovery that high-redshift DSFGs almost
universally do not exceed the Eddington limit even at $\sim$500pc scales, as
accessed by very high-resolution ALMA imaging and/or strong lensing (e.g.
Bussmann et al. 2012, 2013; Rybak et al. 2015; Enia et al. 2018; Hodge et al.
2019; Dudzevičiūtė et al. 2020). In particular, in the case of Rybak et al.
(2015), the star formation surface density of SDP.81 is mapped at sub-50pc
scales (unprecedented at high-$z$), giving a maximum of $190\pm
20~{}{M}_{\odot}~{}{\rm yr}^{-1}~{}{\rm kpc}^{-2}$, significantly under the
theoretical limit. This also appears to be the case for the PASSAGES sample
studied so far, as shown in Fig. 6 and 7 and Table 3. In contrast, turbulence-
based ISM modeling of the full sample by Harrington et al. (2021) revealed
several objects approaching $L_{\rm IR}/M_{\rm ISM}=500L_{\odot}/M_{\odot}$,
but most of these are coincidentally excluded from our measurement, owing to
their lack of either ALMA observations or fully-developed lens models.

In line with Hodge et al. (2015, 2019), we also estimate the peak values of
$\Sigma_{\rm SFR}$ in Table 3 by normalizing the 1 mm flux density to equal
the SFR, effectively assuming that the spatial distribution at sub-mm
wavelengths directly traces the distribution of star formation (also Hatsukade
et al. 2015; Tadaki et al. 2018; Sharon et al. 2019). Unlike Hodge et al.,
however, this calculation applied towards lensed galaxies is complicated by
the varying source-plane PSF. To mitigate this, we compute the 99.7th
percentile (i.e. the 3$\sigma$ confidence interval of pixel values) as a lower
bounds on $\Sigma_{\rm SFR,peak}$. This helps to mitigate any spurious non-
physical artifacts introduced by the lensing reconstruction, and also accounts
partially for the correlation of adjacent pixels within a beam size. In future
work, with higher angular resolution observations, it will be worthwhile to
construct smoothed (i.e. uniform PSF) maps of $\Sigma_{\rm SFR}$, but the
current ALMA data are only marginally able to resolve this distribution. For
this reason, the peak values given in Table 3 are usually consistent with
globally-averaged values, which is what one would expect for source sizes
comparable in extent to the resolving beam.

While it is clear that 100pc-scale imaging is essential to properly assess
these maximum surface densities, there may be other factors responsible. One
such factor, as pointed out by Hodge et al. (2019), might be the assumption of
a single dust temperature, leading to a possible underestimation of the
overall SFR (Berman et al. 2022). For example, Calistro Rivera et al. (2018)
find evidence for linear temperature gradients in four DSFGs, decreasing from
the center to outer parts of the disk, which would directly affect the peak
values of $\Sigma_{\rm SFR}$. Still, it is not likely that this effect can
fully account for the order-of-magnitude discrepancy in $\Sigma_{\rm SFR}$
below the Eddington limit. $\Sigma_{\rm SFR}$ scales linearly with star-
formation rate, but scales with inverse-square dependence on galaxy size, so
the latter is a more facile explanation. Additionally, Andrews & Thompson
(2011) and Murray et al. (2010) point out that intermittency and non-
uniformity of star formation over the extent of a galaxy can lead to a global
star-formation surface density that appears significantly more sub-Eddington
than if higher spatial and temporal resolution were achieved. Ultimately, a
larger sample of sub-100pc dust continuum observations of DSFGs is required to
make more substantive claims, and we intend to pursue this in the future with
a more thorough, higher-resolution analysis of the PASSAGES sample.

Figure 11:  Magnification $\mu$ vs. mass enclosed within the Einstein radius
(as derived from our lensing models) for the PASSAGES sample (wavelength
denoted by color). Open symbols denote 3 objects for which the foreground
redshift is preliminary; uncertainties include the wider ranges given in Table
2, which encapsulate some of the effect this has on enclosed mass
determination. While shown in this plot, they are not included in our analysis
in §4.5. Purple lines connect radio and far-IR measurements where both exist
for the same object. A fiducial line segment shows the slope of
$\mu\propto\sqrt{M_{\rm lens}}$, a rough expected scaling relation.

###  4.5 Connection of lensing halo mass to magnification

The distribution of Einstein radii (and masses enclosed within these radii) is
summarized in Table 2 (also Lowenthal et al., in prep.). As remarked by Frye
et al. (2019) (see their Figure 1), the objects selected via Planck-Herschel
or Planck-WISE have a tendency for larger far-IR flux and Einstein radii,
perhaps resulting from the larger-area footprint of Planck. The median
Einstein radius is found to be $\sim 1.5^{\prime\prime}$, in contrast with the
smaller median of $\sim 0.8^{\prime\prime}$ for the Herschel lenses of Eales
(2015) and Amvrosiadis et al. (2018). While the clear majority are galaxy-
scale lenses, the total halo masses enclosed within the Einstein radius for
PASSAGES ranges from $1\times 10^{10}$ to $\sim 3\times 10^{13}~{}M_{\odot}$,
covering galaxy, group, and low-mass cluster scales. There is not a systematic
effect where the most luminous objects in the sky are all high-magnification
cluster-lensed DSFGs; galaxy-scale lenses are just as capable of producing
lensed DSFGs with large fluxes.

The source-plane area inside the caustic curves that is capable of producing
high magnifications (or lensing cross-section; Turner et al. 1984; Meneghetti
et al. 2003) has been found to correlate tightly with the Einstein radius
(regardless of definition; Meneghetti et al. 2011). By extension, as
$\theta_{\rm Ein}\propto\sqrt{M}$, this cross-section also should correlate
with lensing mass. This effect might be at play in Fig. 11, which shows a
positive correlation between magnification and the mass enclosed within the
Einstein radius. As a general rule-of-thumb, magnification scales
approximately as $\mu\propto\sqrt{M_{\rm lens}}$ for a circularly symmetric
lens (assuming fixed source position; see e.g. Narayan & Bartelmann 1996; Frye
et al. 2019). Additional factors play a large role in determining
magnification—including perhaps most notably the alignment of background
relative to the foreground mass profile—so there should be large scatter in
this relation. We thus do not consider it instructive to parameterize it, but
we remark that the distribution of PASSAGES objects in this work is broadly
consistent with this trend (shown in Fig. 11).

Building on this, since we recall that more extended source-plane objects may
also be subject to lower overall magnification factors than compact ones (e.g.
Hezaveh et al. 2012), then a useful parameter might be the ratio of galaxy
angular size (effective radius) to the angular Einstein radius of the lens,
$R_{0}\equiv R_{\rm e}/\theta_{\rm Ein}$. If this ratio is close to unity, the
size of the background object would be comparable to the areal extent of the
caustic network. Larger values would necessitate (in general) that the source
extends further into lower-magnification regions, resulting in lower, more
diluted magnifications. Smaller values, on the other hand, would suggest that
the background object can carefully align inside the caustics. As remarked in
§4.3.1, however, smaller source-plane objects have a lower probability of
being near to the high-magnification region of the source-plane when