PASSAGES: The Wide-ranging, Extreme Intrinsic Properties of Planck-selected, Lensed Dusty Star-forming Galaxies

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-resolution 1.1 mm Atacama Large Millimeter/submillimeter Array and 0.″4 5 cm Jansky Very Large Array continuum imaging, in tandem with 1.6 μm Hubble and optical imaging with Gemini-S. Magnifications range from μ = 2 to 28 (median μ = 7), yielding intrinsic infrared luminosities of L IR = 0.2–5.9 × 1013 L ⊙ (median 1.4 × 1013 L ⊙) and inferred star formation rates of 170–6300 M ⊙ yr−1 (median 1500 M ⊙ 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 e = 1.7–4.3 kpc; median 3.0 kpc) but consistent with L IR–R 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.

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 , reflecting star formation rates > 100 − 1000 M 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 Eddingtonlike 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 starforming 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. 2006Narayanan et al. , 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(Negrello et al. , 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 − 100pc scales where starforming 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 Kcorrection in the submillimeter (e.g.Blain & Longair 1993), owing to steep Rayleigh-Jeans tails in the restframe 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 shorterwavelength (and therefore brighter) parts of the SED.For example, at 850µm in the observer-frame, these effects roughly equilibrate, and the flux for an object of fixed luminosity appears uniform from z ≈ 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 DS-FGs.
The Planck All-Sky Survey to Analyze Gravitationallylensed Extreme Starbursts (PASSAGES) sample was introduced in Harrington et al. 2016 andBerman 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 (µL IR = 0.1 − 3.1 × 10 14 L , median 1.2 × 10 14 L ; also Harrington et al. 2021).Even assuming a fiducial magnification factor of µ ≈ 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 IR .Such objects with IR luminosity in excess of 10 13 L have typically been referred to RGB color panels show imaging with Gemini-S r and z filters (where available) with HST H-band (where H is shown as red, z as green, and r 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).
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 spectroscopicallyconfirmed 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.
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-arcsecondresolution millimeter imaging-will be described in future works.Some models have already been presented publicly for PASSAGES, including: PJ020941.3(9io9; Geach et al. 2015Geach et al. , 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.

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 followup 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 ≈ 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).

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 (λ effective = 1.54µm,H-band) was used to image each member of our sample for one orbit each (with an additional orbit in F110W, 1.15µm, 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σ sen-  Berman et al. 2022;Harrington et al. 2021;(3) Geach et al. 2015;(4) Harrington et al. 2016;(5) Cañameras et al. 2015, 2017aand Frye et al. 2019.sitivity of m AB ≈ 28.7 was reached for an unresolved point source.Due to the sub-pixel dithering setup (0.572 spacing), the native pixel scale was improved by a factor of 2 to 0.065 (with a point spread function FWHM = 0.2 ).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-ofview (141 × 125 ), we can apply a 3-dimensional astrometric solution.In fields with fewer matched Gaia detections, we apply only 2-dimensional shifts in RA/Dec.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 arcs1 .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.

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 3hour tracks, to improve uv-coverage.The most extended A-configuration was used to provide optimal resolution ( θ min θ maj ≈ 0.3 − 0.7 , median 0.4 ).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 (Mc-Mullin 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σ threshold.Achieved sensitivities ranged from 1σ = 2.6 − 11.7 µJy (median 2.9 µ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 PAS-SAGES 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 6GHz ) is reported in Table 4.
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 filter for 1500s to achieve S/N = 10 for an object with r = 25.0, and in the z filter for 2100s to achieve S/N = 10 for a z = 23.0 object.These limiting magnitudes are each chosen to be one magnitude fainter than a typical earlytype L cluster member at z ∼ 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 × 5.5 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 (dragons2 ) software, as will be described in Cooper et al., in prep.

GRAVITATIONAL LENS MODELING WITH LENSTOOL
All gravitational lens modeling in this work was performed with the publicly-available software, lenstool 3  (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, where β, θ, and α(ξ) 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 ξ.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: where G is the Newtonian gravitational constant and c is the speed of light.Σ is surface mass density evaluated at the location of the impact parameter ξ, after projecting the volumetric density distribution on to the lens plane (Schneider et al. 1992).Thus, provided the locations θ (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 PAS-SAGES 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 PAS-SAGES sample will be formulated using both parametric modeling (e.g.lenstool) and light-traces-mass models, in an approach similar to Zitrin et al. (2009Zitrin et al. ( , 2015) ) and Frye et al. (2019) 4 , 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.

Visual identification of multiple image families
Identifying multiple image families with longwavelength, 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 µ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 2dimensional 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. 2012Bussmann et al. , 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σ 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 optimization5 .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 galaxyor 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.

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 χ 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 PAS-SAGES objects in this work are primarily galaxy-scale lenses, fewer free parameters are involved than clusterscale 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 κ, takes the form where f ≡ b/a is the axis ratio (for semi-major and semi-minor axes a and b; 0 < f ≤ 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 ρ ∼ 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. 2006Koopmans et al. , 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), andSpilker et al. (2016).The profile can be fully parameterized by position (α and δ, or x and y), projected 2dimensional ellipticity e = (a 2 − b 2 )/(a 2 + b 2 ), position angle (PA or θ; measured counterclockwise from east), and velocity dispersion, σ 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 σ ≈ 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 PAS-SAGES objects (e.g., Geach et al. 2015Geach et al. , 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 parameters6 .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 ≈ 95% of galaxy-scale, 10 12 M halos at z ∼ 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 ∼ 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 θ, and velocity dispersion σ 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 MCMCsampled 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 χ 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 χ 2 value for all median solutions7 is χ2 med = 1.9, and the median for the best solutions is χ2 best = 0.3.This suggests that the goodness-of-fit of our lens models are reasonable, in aggregate, although the individual χ 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.

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, nonuniform 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. 2018Sharma et al. , 2021)).These individual reconstructions are compared visually to ensure that the total reconstruction of all image-plane light is appropriate and self-consistent.

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 wavelengthdependent, 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 Table 2. Summary of strong lensing properties derived from the models described in Table 6, including weighted total magnification µtot for different wavelengths and the equivalent Einstein radius θ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 µ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). −11.9 3.61 +0.36 1.17 +0.11 13.6 +15.7 −6.7 6.52  1), so approximate bounds of enclosed masses are also given for a set of redshifts representative of the range of the full sample.* * 870 µm (Band 7) measurement (see §2.4).The 1 mm magnification is µ ≈ 11.8 ± 2.1 for the 0.07 image, but the 870 µm value is adopted for µ 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;slightly different magnification, measured as the area above a 3σ 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 misidentified 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 (∼ 100 pc) or cut radius (∼ 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.

Measuring Einstein radii
The Einstein radius θ E is a characteristic scale relating directly to the total enclosed mass of the foreground lens inside a projected circle in the sky M (< θ E ), with some additional dependence on the redshift geometry of the lens and source planes: 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).θ 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 θ E from lens models.Here, we opt for the equivalent Einstein radius, θ Ein,eq , defined as the radius enclosing an average surface mass density equal to the critical surface density (or equivalently, where the average convergence κ 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(Meneghetti et al. , 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 θ Ein in Table 2, as the median value and asymmetric 1σ uncertainties.The Einstein radii are computed analogously to the magnification, where 100 samples are drawn from the modelingderived 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.

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. 2012Bussmann et al. , 2013Bussmann et al. , 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 size8 .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 FWHM = 2×R e .
In our method, we therefore use the geometric mean of the semi-major and semi-minor axes as the effective radius9 .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.

PASSAGES sources have modest magnifications
Eleven out of the 15 objects in this work for which we present new lensing models were previously undiscovered10 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 resolution) and JVLA (6 GHz, 0.3 − 0.7 resolution) observations for all members.For 11 out of 15, we also used ALMA imaging (1 mm, 0.4 − 0.8 resolution), and for 8 out of 15, we used Gemini r and z 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.(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.
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.
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 dis-tribution: 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 ), 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 ∼ 0.4 resolution comfortably resolves the lensed images in most cases, but does not offer much detail of substructure inside each of the images.Higher resolution (∼ 0.1 ) 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 wellseparated 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-thanexpected 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µm 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 (θ ≈ 0.5 ) was insufficient to spatially resolve emission within individual lensed images, leading to halflight radius estimates skewed to larger values, in turn giving smaller magnification factors.Indeed, later works showed that higher-resolution imaging (≈ 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).

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 PAS-SAGES (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.
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 µ 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 µ 6GHz factors can thus be used to roughly estimate intrinsic L 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 yr −1 (median 1500 M 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 be-tween SFR and M , 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 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 [SFR/M yr for galaxies with log 10 (M /M ) 10.5.For this work, the inner 68% confidence interval in star-formation surface density is Σ SFR ≈ 20 − 160 (see §4.4), suggesting distances to the main sequence of ∆log 10 (sSFR) 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 Σ 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 opportu-nity 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.

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 e = 1.7−4.3kpc (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.

Possible source size bias?
The intrinsic sizes of the DSFGs that make up the PAS-SAGES 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.µ > 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 lens ≈ 0.37, Robertson et al. (2020) found source-plane cross-section solid areas at z s = 2 satisfying |µ| > 10 to be σ S lens ≈ 0.008 and 0.08 arcsec 2 for halo masses of 10 12 and 10 13 M , respectively.For large magnification thresholds (µ 0 1), this cross-section drops off as σ S lens ∝ µ −2 0 (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 (µ ∼ 10) could preferentially magnify diffuse (or extended) components in the source plane, contrary to the common intuition (e.g. ).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 Re > 0.3 have magnifications µ ≤ 10.
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 (µ < 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) deter- mined that the cross section for µ > 10 was optimally maximized when the source-plane object had an angular size of ≈ 30% the Einstein radius (in contrast to ≈ 15% for the threshold µ > 20).We make further comparisons to these predictions in §4.5.
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µm magnifications for SPT lenses, along with Herschel lenses from Bussmann et al. 2013 andBussmann 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 (µ 870µ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 µm SCUBA-2 Cosmology Legacy Survey (Simpson et al. 2015a,b; median FWHM 0.30±0.04 ) and 1.1-

IR
. 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 LIR − λ 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 yr −1 kpc −2 are shown as dashed lines.The sample from this work is located at higher luminosities (LIR 10 13 L ), but effective radii are in general correspondingly larger, and thus consistent with the previous trends.
mm AzTEC sources (Ikarashi et al. 2015(Ikarashi et al. , 2017; median FWHM 0.20 ± 0.04 and 0.31 ± 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 0.2 , a region of the parameter space for which our sample contributes significantly.Objects with sizes larger than 0.3 all have lower magnifications, µ ≤ 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 µL IR ∼ 10 14 L 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 µ and R e .In tandem with the positive slope of L 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 µ and L 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.

Intrinsic far-IR source size
With this work, we seek to place the PASSAGES objects modeled in this work in the context of other submmselected 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-005343or SDP.130, R e = 4.14 ± 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 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 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 andSpilker 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 deg −2 , in comparison with the Herschel samples' density of 0.02−0.1 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 ∼ 3.5 − 5 primarily).While the angular size does not evolve significantly for a given physical size between z ∼ 4 and z ∼ 2, Fujimoto et al. ( 2017) find a detectable evolution in physical size for a given L IR (L IR ∼ 10 12 − 10 13 L ), 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. (2013Bussmann et al. ( , 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(Swinbank et al. , 2011;;Danielson et al. 2011;Hodge et al. 2012Hodge et al. , 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 counterevidence 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 (∼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 (∼ 0.5 ) 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 Table 3. Lensing-corrected source-plane effective radii at 1.1 mm and intrinsic total IR luminosities LIR, using µ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, Σ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
Re,1100µm Re,1100µm θmaj θmin 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 ± 0.2 kpc vs. R e = 1.6 ± 0.3 for star-formation dominated or R e = 1.5±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. for L FIR 10 12 L , 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 ∝ L 0.37±0.03FIR (in a much smaller sample of 18 DSFGs, median R e = 0.32 ± 0.09 kpc).However, these trends have yet to be tested for a large ensemble of objects at L IR ∼ 10 13 L , 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 DS-FGs 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: , with which the PAS-SAGES objects concur (Fig. 7).In keeping with the Stefan-Boltzmann relation, this would necessitate that T d ∼ R 0.2−0.3 e , 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 ] ≈ −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 thin13 (which fortuitously allows for the direct estimation of ISM mass from flux; Scoville et al. 2016Scoville et al. , 2017)).

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 endproduct 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 FIR is computed as where we determine L 1.4 GHz from 6 GHz flux as for luminosity distance D L in meters and spectral index α (assumed here to be 0.8; Condon 1992) such that S ν ∝ ν −α .These luminosities (corrected for lensing magnification µ radio ) are included in Table 4, along with the corresponding values of q FIR .We find values ranging from q FIR = 2.0 − 3.1, with a median of 2.3 and 1σ dispersion of 0.3, in agreement with the mean found by Yun Table 4. Properties of JVLA 6 GHz continuum imaging, including natural-weighting flux (S6GHz, uncorrected for lensing), inferred 1.4 GHz luminosity, and the derived radio-FIR correlation parameter, qFIR.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
S6GHz Note- † Lensed images of PJ231356 are blended with a foreground double-lobed radio jet at 6 GHz, so the 6 GHz flux and effective radius of the lensed DSFG cannot be determined.The reported flux includes significant contribution from the bright foreground object.
et al. ( 2001) for local IR-selected galaxies, q 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σ dispersion of q FIR ≈ 1.9 − 2.9.Unsurprisingly, none of the star-formation-dominated PAS-SAGES galaxies fall below the q 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 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 ∼ 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 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 starforming galaxies could result in a steepening of the evolution with redshift.Alternatively, it is possible the calculation of q FIR for PASSAGES galaxies is weakened by the assumption of a single radio spectral index α, 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 ν ∝ ν −0.1 ) may contribute non-trivially at restframe ν 20 GHz (Condon 1992).
Under a different assumption of a flatter spectral index α = 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 FIR greater than those reported in Table 4 by a small amount, q FIR,α=0.5 − q FIR,α=0.8 = 0.3 − 0.4.This change in q FIR is not insignificant, but it is smaller than the dispersion we find in q FIR for our sample, and not much larger than the uncertainty on q FIR for any individual object, so we conclude that the possible effect on our interpretation is minimal.
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  Helou 1998;Murphy et al. 2006b) might be responsible, relative to the smaller ∼ 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 IR and Σ 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 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 ∼ 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 Σ 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 Σ 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, Σ 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 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 halflight sizes.
4.3.4.Redshift evolution of galaxy size  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≤ log 10 M < 12.0), mid (10.5 ≤ log 10 M < 11.0), and low (10.0 ≤ log 10 M < 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).Fujimoto et al. (2017) found modest evolution in the R e − L FIR relationship between z = 1 − 2 and z = 2 − 4, and a more significant evolution (by a factor of ∼ 2) to z = 4 − 6.This decrease in effective radius with redshift for a given luminosity is consistent with restframe UV/optical samples.For example, Shibuya et al. (2015) found that the UV/optical effective radius of UVluminous star-forming galaxies evolves with redshift as R e,UV /kpc ∝ (1 + z) −0.84±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,radio /kpc ∝ (1 + z) −0.3±0.3 from z ≈ 0.3 − 3, while Lindroos et al. 2018 found steeper evolution at 1.4 GHz of R e,radio /kpc ∝ (1 + z) −1.7 out to z ≈ 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 den-sity 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 Herschelselected 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≤ log 10 M < 12.0), mid (10.5 ≤ log 10 M < 11.0), and low (10.0 ≤ log 10 M < 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σ 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 ∼ 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.

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 ∼ 500 L /M 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 radi-ation 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).
Yet, a vexing result has been the discovery that highredshift DSFGs almost universally do not exceed the Eddington limit even at ∼500pc scales, as accessed by very high-resolution ALMA imaging and/or strong lensing (e.g.Bussmann et al. 2012Bussmann et al. , 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 ± 20 M yr −1 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 IR /M ISM = 500L /M , but most of these are coincidentally excluded from our measurement, owing to their lack of either ALMA observations or fullydeveloped lens models.
In line with Hodge et al. (2015Hodge et al. ( , 2019)), we also estimate the peak values of Σ 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σ confidence interval of pixel values) as a lower bounds on Σ 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 Σ 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 Σ SFR .Still, it is not likely that this effect can fully account for the order-of-magnitude discrepancy in Σ SFR below the Eddington limit.Σ SFR scales linearly with starformation 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.

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 ∼ 1.5 , in contrast with the smaller median of ∼ 0.8 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 × 10 10 to ∼ 3 × 10 13 M , 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 θ Ein ∝ √ 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 µ ∝ √ M 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 ≡ R e /θ 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 highmagnification region of the source-plane when positioned at random.Ideally, one can describe a balance between these two competing effects.
For a general SIS lens, de Freitas et al. ( 2018) derived analytic, perturbative solutions for the source-plane crosssections σ µ -or the solid angle of the region with magnifications exceeding some threshold µ thr -as a piecewise function of R 0 : where the joining point R J = 2.15/µ thr marks a transitional point to the Einstein ring regime.As σ µ is relative to the Einstein radius, it can be multiplied by R 2 0 to obtain an absolute solid angle (e.g. in arcsec 2 ).This treatment implies that the peak area for µ thr = 5 is approximately R 0 = 0.6 (i.e., source size 60% of the Einstein radius).
These analytic curves are shown in the bottom panel of Fig. 12.For given magnification µ, the cross-sectional area σ µ essentially describes the probability of a given R 0 ratio.We find that the peaks for the selected magnification thresholds (µ thr = 5, 10, and 20) are approximately concordant with the negative correlation seen in the upper panel of Fig. 12 for the PASSAGES sample, as indicated by the conjoined vertical/horizontal dotted lines.Specifically, the local maxima of σ µ in Equation 8are found to occur at µ(R 0 ) ≈ 2.9522/R 0 , shown as a solid black line in Fig. 12.There is perhaps greater deviation from this relation for smaller values of R 0 , but there are also fewer objects occupying this regime, so the deviation is not a robust conclusion.
The relative probability of lower-magnification lenses is higher, as σ µ ∼ µ −2 thr , which would suggest that they would be more numerous.However, this will not be fully reflected for the PASSAGES sample, as the selection is not uniformly sensitive to all magnifications.This is because selection by apparent flux serves to maximize the multiplicative product of magnification and intrinsic luminosity (the latter of which is itself correlated with source-plane size, as discussed in §4.3.2).For this reason, the recovered distributions of intrinsic luminosities, magnifications, Einstein radii, and source-plane sizes are all interconnected for PASSAGES and other lensing samples defined by flux.
In future work (Englert at al., in prep.), we intend to place the PASSAGES sample in context with other studies that have used magnification and Einstein radius distributions as constraints on the halo mass function and other cosmological parameters (e.g., Eales 2015;Amvrosiadis et al. 2018).This is an analysis which benefits immensely from a large sample with consistent selection effects. .This ratio can serve as a proxy for the compactness of the background object relative to the foreground lens.There is an apparent negative correlation between this ratio and magnification, suggesting possibly that smaller background objects can align more precisely with the highest-magnification regions in the source plane (which increase in size with Einstein radius).The redshift uncertainty of the foreground for 3 objects does not strongly influence these values, as θEin is computed in angular units.Bottom panel: The predicted source-plane cross sections σµ as a function of the same ratio, for magnification thresholds of µ thr = 5, 10, 20, using the analytic perturbative solutions calculated by de Freitas et al. (2018).The relative solid angle σµ can be converted to absolute angular units by multiplying by θ 2 Ein .Dotted vertical lines indicate the optimal source size for given magnification thresholds (shown as the conjoined horizontal lines in the upper panel) to maximize the crosssectional area.These maxima are located at µ ≈ 2.9522/R0, which is shown as a black solid line in the top panel.The color gradient in the top panel shows the distribution of σµ by computing these curves as a function of R0 for a range of continuous µ values.This indicates that smaller magnifications have larger cross sections, as expected, but there is also a thin diagonal band of elevated cross sections around the black solid line.

CONCLUSIONS
We have presented detailed lens modeling for 15 members of the PASSAGES (Planck All-Sky Survey to Analyze Gravitationally-lensed Extreme Starbursts) sample of dusty star-forming galaxies (z = 1.1 − 3.3) that were introduced in Harrington et al. (2016) and Berman et al. (2022).To do this, we gather complementary information from an extensive set of recent multi-wavelength sub-arcsecond continuum imaging with ALMA, JVLA, HST, and Gemini-S.We implement parametric lens modeling with lenstool, where priors on the nature of the foreground mass distribution are established by optical/near-IR imaging, and multiple image positions of background objects (from imaging at all wavelengths) provide modeling constraints.In analyzing the results of this modeling and the resultant source-plane reconstructions, we find: • The 1 mm and 6 GHz magnification factors of the DSFGs modeled in this sample range from µ = 2 − 28, with more than half at lower magnification, µ < 10.Correcting the apparent infrared luminosities from Harrington et al. (2016) and Berman et al. (2022) with these magnification factors, we find intrinsic values of L IR = 0.2−5.9×10 13L (median 1.4×10 13 L ), placing them solidly within the regime of hyperluminous infrared galaxies, or HyLIRGs (Fig. 3).The corresponding inferred star formation rates (also corrected for lensing) span 200 − 6300 M yr −1 (median 1500 M yr −1 ).Their extreme, rare properties are likely to be a direct result of the all-sky selection method for PASSAGES.While initially assumed to be starbursts based on their large star formation rates, evidence from their molecular gas depletion times and star formation rate surface densities indicates that-while nearly all above the center of the star-forming main sequence-they largely are not predicted to surpass a ∆MS > 0.6 dex threshold typically imposed to define starbursts.Future work to attempt to constrain stellar masses is thereby warranted in order to directly establish their relation to the main sequence.
• In characterizing their de-lensed, spatial extents in the source plane at rest-frame far-IR and radio, we find that the PASSAGES DSFGs are generally more extended than typical DSFGs (both lensed and unlensed).Since they appear to be consistent with size-luminosity scaling relations for unlensed submm-bright objects (Fujimoto et al. 2017), we interpret their larger sizes to be a consequence of their higher star formation rates (or vice versa); see Figs. 6 and 7.They are consistent with a framework of some DSFGs being the result of active growing of and accretion of gas onto a starforming disk, with a wider distribution of star formation than low-z ULIRGs, which form stars in compact nuclear regions.
• We compute the radio-FIR correlation parameter q FIR by assuming a single radio spectral index; the range of values for q FIR is in line with those from Giulietti et al. (2022) for Herschel-selected lenses, but we find no clear evolution with redshift, as predicted by some recent works.We also find that PASSAGES objects do not show signs of significant contributions from AGN, which might be indicated by a radio excess, q FIR < 1.8.
• In comparing their radio vs. sub-mm sizes, we find some weak dependence on L IR , Σ SFR , and the radio-FIR correlation parameter q FIR , with the more luminous and more densely star-forming galaxies having preferentially more compact 1 mm sizes (or more extended 6 GHz halos).Relative to q FIR , we find that galaxies closer to being AGNpowered (i.e.showing a radio excess) had generally smaller radio sizes, which might be indicating that deeply-buried AGN are leading to more compact half-light radii.At present, our conclusions are limited by the small number of 10 objects with both radio and FIR sizes measured.
• The larger sizes of PASSAGES galaxies (relative to the broader population of DSFGs) may also owe to their concentration around Cosmic Noon, z < 3.5, in contrast to higher-z DSFGs like Spilker et al. (2016).Nonetheless, we observe no clear trend in size as a function of redshift for the narrow range covered by these PASSAGES objects (Fig. 10).This agrees with recent results by Wang et al. (2022) for the FIR size of stacked main-sequence galaxies and Jiménez-Andrade et al. ( 2019) for the radio size of star-forming galaxies, which evolves only shallowly with redshift.
• There may be some size bias present in the sample (Fig. 5), by which more compact objects are capable of producing higher lensing magnifications (Hezaveh et al. 2012;Spilker et al. 2016).We parameterize this effect by comparing magnifications with the ratio of the angular size of background objects to the Einstein radius of the foreground, and find a modest downward trend as expected (Fig. 12), which is in line with predictions for isothermal lenses (e.g. de Freitas et al. 2018).
That the PASSAGES sample is not fully coincident with the distribution of lensing cross-sections σ µ as a function of the size ratio and magnification is an indication of the selection function at play, as we are not uniformly sensitive to all magnifications.Given the selection of these objects by flux, they are likely of an optimal size that is simultaneously sufficiently large to host high intrinsic luminosities and sufficiently small to align with higher-magnification regions of the source plane.
• Since the spatial extents of their star-forming regions are proportionally larger to match their extreme luminosities, the star-formation surface densities for PASSAGES are significantly sub-Eddington on global scales (in agreement with most DSFGs for which this measurement has been made).There is still certainly a possibility that their rapid stellar mass assembly is substantially regulated by radiation pressure, but this effect must be on sub-kpc scales not captured by these current 0.4 (image-plane) resolution ALMA images, perhaps over short timescales that are not contemporaneous between sites of star formation.For these observations, ALMA does not significantly resolve the far-IR continuum, and so the peak values of Σ SFR in each object are not much larger than the respective globally-averaged values.In the future, we hope to repeat this analysis with higher angular resolution from ALMA to properly assess the densities of star formation, which we expect may be clumpy and extended beyond just nuclear regions.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(Tody , 1993 Harrington et al. (2016) and Appendix B of Berman et al. (2022) provide notes on each member of the larger PASSAGES sample (in addition to Table 1 of Harrington et al. 2021).In this section, we review the information necessary for our lens models, including spectroscopic redshifts of the background DSFGs, photometric/spectroscopic redshifts of foreground objects, and multiple image morphologies, summarized in Table 1.
PJ011646.This galaxy-scale lens system appears to consist of multiple background components lensed by an elliptical galaxy with a spectroscopic redshift of z spec = 0.555 determined with VLT MUSE observations (Kamieneski et al., in prep.).There is a smaller, secondary foreground galaxy of unknown redshift to the northeast, but it has minimal impact on the lensing morphology and so we exclude it from our model.The foreground elliptical is modeled with our standard approach: an SIE profile with centroid, ellipticity, orientation, and velocity dispersion all kept as free parameters.The DSFG, with a CO(3-2) detection at z = 2.125 by LMT/RSR (Berman et al. 2022), is detected as two primary arcs in ALMA and JVLA imaging and appears to spatially coincide with a bright double image detected by HST (labeled 2ab in Figure 1).Another image family (labeled 1abcd) is a clear quad image, nearly spatially coincident with family 2ab but noticeably bluer in the RGB image with Gemini r and z (Fig. 1).It is not entirely clear if this image family is located at the same redshift as 2ab, but since the two families have very similar Einstein radii and appear to be lensed by virtually the same foreground mass profile, then this should dictate that the background objects lie at similar redshifts.Within the Einstein radius uncertainty for this object (see discussion in Section 3.5), θ Ein,eq = 2.34 − 2.41 , and holding other quantities fixed, this corresponds to a redshift uncertainty of z ≈ 2.0 − 2.4 (assuming z lens = 0.555).For the purposes of deriving our lens models, it is sufficient to assume that they lie at the same redshift.Given the proximity of this quad component and the doubly-imaged DSFG component in the source plane ( 1 , 8.5 kpc at the CO-derived redshift), it seems most likely that these two components may be interacting (or else a very serendipitous alignment of two objects around z ∼ 2).We also note that these multiple components are aligned with the foreground in such a way as to cover most of the caustic curve region so that there is a very complex set of concentric, partial Einstein ring structures in the image plane, offering excellent constraints on the mass of the foreground elliptical, (1.54 ± 0.03) × 10 12 M within the Einstein radius at z lens = 0.555 (see Table 2).
PJ014341.This compact galaxy-scale lens is one of the lowest-redshift members of our sample, with a CO(2-1) detection from the background DSFG at z = 1.096.The lensing foreground galaxy is detected spectroscopically at z spec = 0.594 from SDSS (Berman et al. 2022).Its small Einstein radius of ≈ 0.5 (on the lower end of our sample, see Table 2) is due primarily to the small separation in redshift space between the lens and source planes.The enclosed foreground mass is (1.3 ± 0.3) × 10 11 M .The redshift geometry also has the consequence of making lensed image identification difficult, as they are likely heavily blended with foreground light in the optical (Fig. 1), and only marginally resolved beyond the beam size in the interferometric ALMA/JVLA images.Nonetheless, the image morphology is consistent with a classic 4-image fold-caustic geometry, with two isolated images and an additional merging pair (Fig. 2).The radio morphology is similar (Fig. 2) but more difficult to interpret due to its lower signal-to-noise, so it is not included as a constraint in this iteration of the model.Given the compactness of the lensing configuration, the source-plane area subtended by the caustic curve is likely smaller than or similar in size to the DSFG, so slight offsets between radio and far-IR emitting regions can result in very different image morphologies.There appear to be other foreground objects nearby, but given the compactness of the lensing configuration, we only parameterize the primary deflector with our standard SIE model.The best-fit model finds a high ellipticity, e = 0.77 +0.15 −0.21 , which suggests these other foreground deflectors may contribute non-trivially.A refined model with more available constraints in the future might allow for a more complex mass distribution.
PJ020941.Also known as 9io9 or the "Red Radio Ring", PJ020941 has been studied at length by Geach et al. (2015Geach et al. ( , 2018)), Harrington et al. (2016Harrington et al. ( , 2019)), Rivera et al. (2019), Doherty et al. (2020), andLiu et al. (2022).A comparison of our lens model with existing models is provided in Appendix D. The target was first revealed by the gravitational lens citizen science project Space Warps (Geach et al. 2015;Marshall et al. 2016;More et al. 2016) before it was identified independently as a strong sub-millimeter source in surveys by Planck (Harrington et al. 2016), the Atacama Cosmology Telescope (Su et al. 2017;Gralla et al. 2020), andHerschel (Viero et al. 2014).Near-infrared i-band and J/K s -band imaging from the VISTA-CFHT Stripe 82 survey (Geach et al. 2017) and CFHT-MegaCam Stripe 82 survey (Moraes et al. 2014) revealed a nearly complete, red Einstein ring (Geach et al. 2015), confirmed with high-resolution H-band HST imaging (Geach et al. 2018, Lowenthal et al. in prep.).The 1.4 GHz eMERLIN and 5 GHz JVLA radio images (Geach et al. 2015) show a partial Einstein ring that largely coincides spatially with the near-IR emission, which is also borne out by ALMA continuum and spectral line images (Geach et al. 2018;Doherty et al. 2020;Berman et al. 2022;and this work).
PJ020941 has a CO(3-2) detection by LMT/RSR at z = 2.553 (Harrington et al. 2016(Harrington et al. , 2021)).It is strongly lensed by an elliptical galaxy (SDSS J020941.27+001558.5) at z fg = 0.202, and by a secondary galaxy assumed to lie at a similar redshift (which is expected to have a non-negligible influence on the lens, Rivera et al. 2019).Both are included in our model, with a standard SIE approach for the primary elliptical, and an SIS model (with position held fixed to the centroid of the luminous component) for the secondary object.
The overall lens morphology consists of a double image, with an extended arc to the west and a less-magnified counterimage to the northeast.Recent Cycle 7 ALMA θ ∼ 0.15 Band 6 (1.1 mm) continuum imaging (2019.1.01197.S, PI: P. Kamieneski; to be presented in forthcoming work by Kamieneski et al. in prep.)provides further information.Embedded inside the extended arc are 3 distinct quad-imaged regions, suggestive of a cusp-caustic configuration, which can be identified by matching clumps and using parity information.This is a commonly-observed lensing configuration where some portion of a galaxy crosses into the inner region of the caustic, which is readily apparent from the morphology in Fig. 2.This quadruply-imaged portion is also evident from subtle surface brightness peaks in the HST image (Fig. 1) that coincide with the 1 mm peaks.
PJ022633.This partial Einstein ring consists of a background object lensed primarily by a group of three foreground objects, assumed to lie at the same redshift.The DSFG CO(3-2) detection is at z = 3.120 (Berman et al. 2022), with the brightest lensing galaxy at z = 0.41 (Wen & Han 2015).Obvious lensing features include a bright, extended partial Einstein ring to the southeast (images 1bc) and a counter-image to the northeast (image 1a), with faint emission connecting the two (more apparent in the optical image).There is some evidence that the lensed image expected to be located on the western side of the primary deflector (opposite the semi-ring) is split into two images by the secondary foreground galaxy to the northwest, visible faintly in H-band and at 6 GHz (images 1de, Figs. 1 and 2).The radio images are nearly spatially coincident with an extended double-lobed radio jet emanating from the primary foreground elliptical, and the HST image is affected by a bright star nearby, so confirmation of this proposed morphology is still pending.However, this interpretation is motivated by the image predictions of a simple zeroth-order lens model including only the mass of the primary and secondary ellipticals.The primary is modeled as an SIE profile, with the secondary as an SIS model, with its position allowed to vary to a small amount.
PJ030510.This object is a galaxy-scale (or small group-scale) lens, with a CO(3-2) detection of the DSFG at z = 2.263 (Harrington et al. 2021;Berman et al. 2022).A preliminary red cluster sequence analysis (Gladders & Yee 2000) yields a foreground redshift of z phot = 0.4 − 0.5.Near-IR imaging by HST and Gemini reveal ring-like features surrounding two foreground galaxies in a dense environment.Our model parameterizes the combined mass contribution from both as an SIE model.ALMA 1 mm continuum imaging has a slightly larger beam size of 0.85 × 0.59 , which leaves the DSFG only marginally resolved.It appears to consist of two lensed images, one of which is notably fainter, on opposite sides of the foreground pair.It is possible that higher-resolution imaging might reveal more of a cusp-like quadruple-image morphology, but at present, a double-image configuration is the best explanation.
PJ105353.This object was discovered independently by Harrington et al. 2016 andCañameras et al. 2015 (also referred to as "the Ruby" or PLCK G244.8+54.9),and has also been studied by Cañameras et al. (2017aCañameras et al. ( ,b, 2021)).We compare our modeling results with theirs in Appendix D. The lensing galaxy is at a higher redshift than most other PASSAGES objects, z = 1.525, with a source redshift of z = 3.005.It is modeled as an SIE profile, with position, ellipticity, orientation, and velocity dispersion as free parameters.As we discuss in Appendix D, CO(4-3) imaging presented in Cañameras et al. (2017a) reveals a complicated image-plane structure, dominated by an apparent quadruply-imaged component with no clear counterpart in rest-frame optical (see also Cañameras et al. 2017a,b), likely due to blending with the foreground as suggested by Frye et al. (2019).High-resolution (0.07 ) 1 mm imaging with ALMA (PID 2015.1.01518.S, PI: N. Nesvadba), which was not directly included in the model by Cañameras et al. (2017a), leads us to a slightly different conclusion on the quad-image configuration.For this work, we utilize the Band 7 (870µm) continuum image to derive the dust magnification, as the observing setup (resolution 0.15 ) is closer to that of the other objects.
PJ112713.This object has not been imaged by AzTEC or ALMA, but a CO(2-1) line was detected by RSR at z = 1.303 (Harrington et al. 2021;Berman et al. 2022).A photometric redshift for the lens was derived to be z phot = 0.42 (Lowenthal et al., in prep.) using SDSS.JVLA imaging indicates two arcs (1ab) separated by ∼ 1 , assumed to be the lensed DSFG, as there are no other bright radio sources within 1 .There is also indication of an incomplete Einstein ring in HST imaging, surrounding a single foreground object, with an Einstein radius that approximately matches the 6 GHz arcs, θ Ein = 0.6 .The foreground galaxy is modeled with the standard SIE approach of this work.Intriguingly, there are two bright peaks in the optical Einstein ring (denoted 2ab) that roughly correspond to minima in the partial radio ring.We interpret this as possibly the result of radio continuum tracing sites of heavily dust-obscured active star formation, whereas the optical peaks arise from less-extincted sightlines.Both 1ab and 2ab are utilized as constraints for our lens model.
PJ113805.A single CO(2-1) line was detected by RSR and interpreted to lie at z = 2.019 (Harrington et al. 2021;Berman et al. 2022).ALMA imaging reveals a compact galaxy-scale lens that is not resolved into multiple components (due in part to a lower-resolution synthesized beam, 1.03 × 0.66 ).The system is however resolved into two primary components by JVLA 6 GHz imaging (Fig. 2).These double images appear to nearly coincide with faint arcs detected by HST (1ab), although the latter are heavily contaminated by light from the foreground lensing elliptical, which has a photometric redshift from SDSS of z phot = 0.52.Since the radio peaks are slightly offset from the optical images, we add them as independent constraints (2ab), although they are likely to provide only minimal new information on the foreground mass.Given the sparse number of model constraints, this object is modeled only as an SIS profile, with position and velocity dispersion as free parameters.
PJ113921.This object was discovered independently by Cañameras et al. (2015), identified as PLCK G231.3+72.2, and later reported in Berman et al. (2022).Our model-derived magnification is compared with that of (Cañameras et al. 2018a) in Appendix D. The RSR CO(3-2) line is interpreted from photometric support to be at z = 2.858, in agreement with additional spectroscopic coverage (Cañameras et al. 2018a;Nesvadba et al. 2019;Harrington et al. 2021).ALMA partially resolves what appears to be a quad-image in 260 GHz continuum (1abcd), centered around a foreground elliptical with a photometric redshift of z = 0.57 (Berman et al. 2022), but likely perturbed by at least one other nearby galaxy.Our model here includes the primary lens as an SIE, and the deflecting galaxy to the west as a simple SIS profile.A faint red feature in HST is possibly affiliated with the background DSFG, as it is close to the ALMA continuum emission, or it may be associated with an apparent foreground spiral galaxy to the southeast.In this iteration of the model, we do not include this feature in constraining our lens model.PJ132630.This object was discovered independently through our sample (Berman et al. 2022) and in the H-ATLAS survey (Bussmann et al. 2013;Yang et al. 2017), known also as NAv1.195.A single CO(3-2) line was detected with RSR at z = 2.951 (confirmed with spectroscopic follow-up by Yang et al. 2017 andHarrington et al. 2021).Bussmann et al. 2013 published a lens model for 340 GHz imaging by the SMA (µ 340GHz = 4.1 ± 0.3), with a lens redshift of z fg = 0.786.We compare their results with our model in Appendix D. ALMA and JVLA imaging reveals that the DSFG is doubly-imaged, with images separated by 3.7 in the image plane (3ab).An optical counterpart appears to separate into two components that surround the long-wavelength emission; we label these image families 1ab and 2ab and include them separately in the model.There is evidence for additional background sources lensed by the foreground elliptical, which we model using a standard SIE, but they are at an unknown redshift and are thus not considered as constraints.
PJ133634.Two CO lines (J = 3−2 and J = 4−3) were detected by RSR at z = 3.254 (Harrington et al. 2021;Berman et al. 2022).A pair of primary lensing galaxies separated by 0.9 have an SDSS photometric redshift of z = 0.26.The 6 GHz continuum reveals a clumpy, partial Einstein ring, centered ≈ 0.2 east of the center of the foreground elliptical (likely due to a small perturbation by the secondary).This object is too northern to be visible to ALMA, but the 1.1 mm continuum from AzTEC (8.5 beam) coincides with the radio continuum, which we thus interpret to originate from the Planck DSFG.The radio ring consists of multiple clumps larger than the synthesized beam, and we identify 3 doubly-imaged families in the eastern counter-image and western semi-ring to be used in the model (1ab, 2ab, and 3ab), making use of the expected parity flip between the counter-image and more extended arc.The foreground mass is modeled as a single SIE profile, as the effect of the secondary perturber is not large enough to break the degeneracy with the primary lens.
PJ144653.This lower-redshift member of our sample has a CO(2-1) line detected by RSR at z = 1.084 (Harrington et al. 2021;Berman et al. 2022).The lensing galaxy has an SDSS spectroscopic redshift of z = 0.493, and has the appearance of a face-on spiral galaxy.However, some of the spiral-like structure may actually be the result of the multiply-imaged background DSFG.At present, it is not feasible to make this distinction without imaging in an additional filter with comparable quality to the HST image.ALMA and JVLA images are consistent with each other and reveal the DSFG to follow a fold-caustic lensing configuration, with two isolated images and an additional merging image pair towards the southeast (labeled 1abc).The images form a partial Einstein ring with a radius of ≈ 0.9 .The foreground is modeled with a single SIE profile.The optimized model has a high ellipticity, e = 0.65 +0.18  −0.25 , possibly because of the small number of available constraints currently available, which necessitates a perhaps too-simplistic model.
PJ144958.This cluster lens has a CO(3-2) line detection at z = 2.153 (Harrington et al. 2021;Berman et al. 2022) with a preliminary foreground cluster redshift measured at z phot = 0.4.Optical emission probed by HST and Gemini appears to be more extended than at 1 mm and 6 GHz, which reveal a set of more point-source-like images.However, there appears at first to be four distinct images all on one side of a foreground galaxy group, with no obvious counterimage.We consider it likely that the middle pair of images (1bc) is in fact a merging pair on either side of the critical curve (see Fig. 2).The northernmost image 1d also appears to be split into four images by another intervening cluster member (4abc), visible only with HST 's higher resolution and unresolved by ALMA and JVLA.As constraints on our model, we include images 1abcd (JVLA/ALMA locations), 2abc for the peaks of the HST images, and 3abc as the fainter tail visible only in optical (adjacent to 2abc).Lastly, we include images 4abc for the point-like sources surrounding the perturbing foreground.The lensing distribution is parameterized by a cluster-scale NFW profile, 3 smaller SIS profiles describing the group of elliptical galaxies interior to the arc, and a fourth SIS profile for the perturbing galaxy at the northern end of the arc.
PJ160722.This member of our sample has a CO(2-1) line observed by RSR at z = 1.482, confirmed with photometric information and subsequent spectroscopic follow-up (Harrington et al. 2016(Harrington et al. , 2021)).HST imaging reveals two apparent foreground lensing objects with a photometric redshift of z = 0.65, determined from the optical to mid-IR SED in (Harrington et al. 2016).A bright radio AGN is detected in the northern foreground (Fig. 2).Both near-IR and radio reveal extended ring-like arcs tightly surrounding the foreground.Four peaks embedded with the surrounding ring are apparent in the radio and, to a lesser extent, in the optical image, which we interpret as a quadruply-imaged region of the source plane.The lensing mass distribution is represented by a single SIE model, given the proximity of the two foregrounds (and therefore strong degeneracies in optimal parameters between the two).
PJ231356.This object has a CO(3-2) line observed by RSR at z = 2.215, confirmed by spectroscopic follow-up (Harrington et al. 2021;Berman et al. 2022).HST /WFC3 1.6µm imaging appears to show a primary lensing galaxy, with a secondary galaxy to the north that appears to have a substantial impact on the lensing morphology, leading to a surprisingly straight arc towards the north.Numerous background arcs are visible in HST, some of which coincide with the arcs in ALMA 260 GHz continuum.Most of these multiple images have not been identified definitively, but the quadruple-image family 1abcd is a sufficient constraint for the model at present.The lensing mass distribution is parameterized by the primary elliptical (SIE) and the secondary perturber to the north (SIS).While the x-position of the primary galaxy is left as a free parameter, we hold fixed its y-position to avoid degeneracy with the secondary lens.

B. LENSED MULTIPLE IMAGE SYSTEMS USED AS CONSTRAINTS
In §3.1, we discuss our approach to identifying multiply-imaged features of the background DSFGs using HST, Gemini, JVLA, and ALMA observations.In Table 5, we summarize these features that were used as constraints for our models, which are also marked on Fig. 1.

C. BEST-FIT LENS MODEL PARAMETERS
In Table 6, we summarize the best-fit lens models for each object.In most cases, we list "median,""best," and "mode" solutions, which represent the set of parameters consisting (respectively) of the median of the posterior distribution, the highest-likelihood solution recovered from the MCMC iterations, and the mode of the posterior distribution.The latter solution is usually employed when the posterior distribution is bimodal, especially in terms of position angle where there may be solutions offset by 90 • .) given are reference positions (to which ∆α and ∆δ are relative), and are typically chosen as the centroid of the brightest foreground.The positional uncertainty σposition is supplied for each object, and is equal to the size of the synthesized beam for ALMA or JVLA (whichever is lowest resolution).RMSim is the RMS deviation between observed and modeled multiple image positions for each field, calculated in the image plane for the lowest-χ 2 solution.The reduced χ 2 ν is reported for the given number of degrees of freedom ν. a In some cases, the mode solution of the optimization is used instead of the median, which is less useful when the posterior distribution is bimodal or if one parameter's optimization is sensitive to the choice of parameter bounds (e.g., when ellipticity is poorly constrained by the multiple images, the posterior tends to follow more of a uniform distribution, so the median tends towards the midpoint of the parameter bounds).For these solutions, the uncertainty in the median is taken to be representative of the uncertainty in the mode.

D. COMPARISON WITH EXISTING LENS MODELS
Several of the objects included in this work have already been the subject of gravitational lens modeling, but we independently derive models here to both incorporate information from new observations and to ensure consistency between models for this subsample.In particular, we examine discrepancies in derived magnifications as a benchmark.
First, PJ020941 (or 9io9) is perhaps the most well-studied DSFG included, with independent modeling efforts by Geach et al. (2015Geach et al. ( , 2018)); Rivera et al. (2019), andLiu et al. (2022).The initial model by Geach et al. (2015) was based on 4 different radio observations and near-IR ground-based imaging (seeing ≈ 0.8 ), and was parameterized by an isothermal ellipsoid and a satellite isothermal spheroid for the secondary foreground, with added external shear.With the gravlens software (Keeton 2011), they modeled the near-IR as a single Sérsic profile and the radio as two Gaussian components, finding a magnification of µ ≈ 10, in agreement with what we find, µ 1mm = 10.5 ± 3.7 and µ 6GHz = 12.0 ± 3.5 (Table 2).Revising this model using CO(4-3) ALMA spectro-imaging and a semi-linear inversion approach (Warren & Dye 2003;Dye et al. 2018), Geach et al. (2018) found a slightly higher magnification of µ CO4−3 = 14.7 ± 0.3.Rivera et al. (2019) incorporated CO(3-2) imaging with NOEMA, with 32 independent velocity channels as constraints, finding a smooth velocity gradient consistent with the structure found by Geach et al. (2018).Also using a 2-component lens model with external shear, they found that the magnification varied strongly with velocity, ranging from µ ≈ 7 − 22, but with a luminosity-weighted mean of µ ≈ 13.Most recently, Liu et al. (2022) used newly-collected adaptive optics H and K s -band imaging, alongside HST /WFC3 F125W, SMA 870 µm, and multi-J CO observations with ALMA in order to build their own lens model.As with our work, Liu et al. used lenstool to fit the primary and secondary foreground with an SIE and SIS profile, constrained by the F125W image.This approach is the one most similar to our own, so an equivalent comparison is possible.The authors found optimized parameters that agreed with all in our model (Table 6) within uncertainties, except for only a small disagreement for position angle (PA = −8 ± 2 • vs. our −14 ± 1 • ), and their source-plane reconstruction is visually very similar to ours (Fig. 2).As a result, Liu et al. found a dust continuum magnification of µ dust = 12.8 ± 0.3 and a stellar component magnification of µ star = 13.6 ± 0.4, consistent with our value of µ dust ≈ 10.5, especially given the different interferometric configurations of ALMA vs. SMA.PJ105353 (or G244.8+54.9, the "Ruby") was previously discussed and modeled by Cañameras et al. (2017a,b); Frye et al. (2019).Cañameras et al. used lenstool to model the z = 1.525 foreground with a pseudo-isothermal elliptical mass density (PIEMD) profile (with the poorly-constrained core and cut radii held fixed at 0.15 kpc and 100 kpc), based on image locations and parities from CO(4-3) observations (angular resolution 0.1 ).However, the authors acknowledged that these image identifications were not completely unambiguous, but settled ultimately on a pair of source-plane objects imaged 2 and 4 times, respectively.An additional 0.07 -resolution Band 6 (1 mm) observation with ALMA (Program 2015.1.01518.S, PI: N. Nesvadba) was not included as a constraint, but our interpretation of the image configuration differs slightly.Whereas Cañameras et al. identify a cusp-like configuration, the 1 mm image appears to instead show more of an Einstein cross, with 4 well-separated images (nearly aligned north, south, east, and west).However, the second, doubly-imaged system is in accord with Cañameras et al.. Unfortunately, the non-detection of the background object with HST (as it is blended with the foreground; Frye et al. 2019) makes confirmation of this morphology difficult.Cañameras et al. find a dust magnification of µ 3mm = 21.8 ± 0.6 (reported in Cañameras et al. 2018a), with magnifications for individual clumps ranging from µ ≈ 7 − 38.In our case, we find a significantly different total value of µ 870µm = 7.6 ± 0.7 (or µ 1mm ≈ 12).Given the difficulty in deriving robust constraints on the model, a discrepancy is not particularly unexpected, especially as both interpretations agree that the source-plane galaxy lies near to a caustic, where the magnification gradient is steep and highly sensitive to the model.On the other hand, values like the Einstein radius are not as sensitive, and our value is in close agreement with that of Cañameras et al. (2017a).
While magnifications are generally robust between various lens modeling approaches, the non-trivial differences between this work and others should not be discounted.Spilker et al. (2016) computes magnification as flux-weighted average over the elliptical model components describing the source plane, whereas Bussmann et al. (2013Bussmann et al. ( , 2015) ) use the ratio of image-plane to reconstructed source-plane flux.On the other hand, Dye et al. (2018) explore the sensitivity of magnification to imaging resolution and depth by calculating magnification as a function of fraction of source-plane flux (i.e.surface brightness threshold).The authors found that, while variation can be small, in the vicinity of a caustic in particular, magnification can vary by up to 25% as a function of interferometric configuration.Our work does not make explicit assumptions of source-plane structure, but also does not take advantage of the entirety of information present in the interferometric images like the visibility-modeling approach of these other works.Future work to apply a visibility-based modeling approach to these objects will allow for a more direct comparison.An additional promising method that we hope to explore in the future involves Regularized Semi-linear Inversion (Warren & Dye 2003;Nightingale & Dye 2015;Enia et al. 2018), by which the source plane is tessellated but not assumed to follow a specific parametric form.
In summary, the general agreement for most objects with prior models (derived from independent constraints) helps to lend credence to their robustness.We interpret any tensions in results-in particular, for PJ105353-to be more influenced by different observing setups than by different lens modeling approaches.

Figure 1 .
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 and z filters (where available) with HST H-band (where H is shown as red, z as green, and r 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).

Figure 2 .
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σ to the 99th percentile of the image (unless this value is less than 5σ).The ALMA image shown for PJ105353 is at Band 7 to ensure a resolution more consistent with that of the other fields.

Figure 3 .
Figure3.Apparent IR luminosity (µIRLIR) vs. magnification µ at 1100 µm (or 870µm if not available) and at 6 GHz (see Table2).Dashed lines indicate intrinsic, magnificationcorrected star formation rates from the Kennicutt 1998 calibration under assumption of a Kroupa initial mass function (IMF), SFR = LIR/(9.4× 10 9 )M 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.

Figure 4 .
Figure 4. Histogram of magnifications for the PASSAGES sample (derived from lens models presented in this work and in the references included in Table2), compared with Herschel lenses fromBussmann et al. (2013Bussmann et al. ( , 2015)),Dye et al. (2015, 2018), and Massardi et al. (2018); and with SPT lenses fromSpilker 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.

Figure 5 .
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 µm for the Herschel (Bussmann et al. 2015) and SPT (Spilker et al. 2016) samples.All measurements are made consistently for ≈ 0.5 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 Re and µ 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 bySpilker 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 Re > 0.3 have magnifications µ ≤ 10.

Figure 6 .
Figure 6.Relation between far-IR effective radius Re (measured at 1 mm) and the magnification-corrected, intrinsic IR luminosity (see Table3) for the 15 PASSAGES objects studied in this work.The gray shaded area shows the Eddington-limited region (1000 M yr −1 kpc −2 ) fromAndrews & Thompson (2011).For comparison, we also plot data points from similar samples of lensed DSFGs studied byBussmann et al. (2013) and Spilker et al. (2016), with corrected LIR values from Reuter et al. 2020 for the latter.For additional comparison, see Figure 5 of Enia et al. 2018 (not shown).

Figure 8 .
Figure 8.The far-IR-radio correlation parameter qFIR as a function of redshift for the PASSAGES sample.The 2σ range of qFIR 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 IRselected 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σ uncertainties.The gray hatch-shaded region at qFIR is the proposed threshold by Condon et al. (2002) for radio-loud galaxies powered by AGN, rather than star formation.

Figure 9 .
Figure9.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 Σ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 qFIR.The radio excess (AGN-powered) regime from Fig.8is again shown as a hatch-shaded region.

Figure 10 .
Figure 10.Effective radius vs. redshift for the PAS-SAGES sample (Tables3 and 4) and the other lensed DSFGs also shown inFig.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 byWang et al. (2022), in stellar mass bins of high (11.0≤ log 10 M < 12.0), mid (10.5 ≤ log 10 M < 11.0), and low (10.0 ≤ log 10 M < 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).

Figure 11 .
Figure11.Magnification µ 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 Table2, 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 µ ∝ √ M lens , a rough expected scaling relation.

Figure 12 .
Figure12.Top panel: Magnifications vs. the ratio Re/θEin, or source-plane effective radius (in arcsec) to Einstein radius (also in arcsec).This ratio can serve as a proxy for the compactness of the background object relative to the foreground lens.There is an apparent negative correlation between this ratio and magnification, suggesting possibly that smaller background objects can align more precisely with the highest-magnification regions in the source plane (which increase in size with Einstein radius).The redshift uncertainty of the foreground for 3 objects does not strongly influence these values, as θEin is computed in angular units.Bottom panel: The predicted source-plane cross sections σµ as a function of the same ratio, for magnification thresholds of µ thr = 5, 10, 20, using the analytic perturbative solutions calculated by de Freitas et al. (2018).The relative solid angle σµ can be converted to absolute angular units by multiplying by θ 2 Ein .Dotted vertical lines indicate the optimal source size for given magnification thresholds (shown as the conjoined horizontal lines in the upper panel) to maximize the crosssectional area.These maxima are located at µ ≈ 2.9522/R0, which is shown as a black solid line in the top panel.The color gradient in the top panel shows the distribution of σµ by computing these curves as a function of R0 for a range of continuous µ values.This indicates that smaller magnifications have larger cross sections, as expected, but there is also a thin diagonal band of elevated cross sections around the black solid line.
Fixed parameters are denoted with an asterisk.Asymmetric error bars indicate 1σ confidence levels and are determined from the resulting MCMC posterior distribution.Ellipticity is defined as e = (a 2 − b 2 )/(a 2 + b 2 ), for semi-major and semi-minor axes a and b.The position angle θ is measured counterclockwise from east.Velocity dispersion σ for singular isothermal ellipsoids are functionally equivalent to physical velocity dispersions.The Right ascension (RA) and Declination (Dec.

Table 1 .
Summary of PASSAGES objects included in this lens modeling study.A brief description of each object is provided in Appendix A. 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) Harrington et al. 2021;Berman et al. 2022)WISE cross-matching.bSourceredshiftsaredeterminedfromCO detections, reported in the listed references.cForegroundlensredshifts(photometricandspectroscopic) are taken from the supplied references.dSpectroscopicforegroundredshiftdetermined with the Multi Unit Spectroscopic Explorer (MUSE) on the Very Large Telescope (VLT).*Preliminaryphotometricredshiftestimate(Cooperetal. in  prep.).For the lensing mass estimates from Einstein radii in §3.5, we assume a fiducial uncertainty in the lens redshift of σz = ±0.2(given the inner 68% range of lens redshifts for the PASSAGES sample, ≈ 0.3 − 0.7;Harrington et al. 2021;Berman et al. 2022); see Table2.
Rodighiero et al. (2011)ses 11 of 10 10.5 M , 10 10.75 M , and 10 11 M 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 τ ∼ 150 − 350 Myr-we obtain ∆MS = −0.02−0.69 for M = 10 10.5 M , ∆MS = 0.02 − 0.78 for M = 10 10.75 M , and ∆MS = 0.08 − 0.88 for M = 10 11 M .This suggests that a clear majority of the PASSAGES objects should have specific star formation rates (sSFR ≡ SFR/M ) 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 = 10 10.75 M ) yields ∆MS = 0.32, which falls below the ∆MS > 0.6 dex threshold for starbursts proposed byRodighiero et al. (2011).

Table 5 .
Multiple image positions used as constraints in the lensing models (presented in Table6), as per §3.1.In some rare cases, the location and orientation of caustic curves as inferred by lensing morphologies are used as constraints.

Table 5 .
(continued) a Assumed to be located at the same redshift as the DSFG.

Table 6 .
Optimized gravitational lens model parameters.Unless otherwise noted, all potentials are parameterized with singular isothermal ellipsoid (or spheroid) models.For NFW potentials, scale radius in arcsec (rs) replaces velocity dispersion as a parameter.While the median of the posterior distribution for each parameter is typically the preferable statistical estimator, multimodal distributions and highly non-Gaussian distributions may be better summarized by the mode of the posterior or the highest-likelihood ("best") sample.The first solution set listed (median, best, or mode) is taken to be the preferred model, but others are listed for completeness.