Where are the Eddington-limited starbursts? Gravitational lensing provides a way forward for sub-kiloparsec views of star formation

Abstract In the past decade, submillimeter surveys have been employed to define samples of gravitationally-lensed dusty star-forming galaxies (DSFGs) at z ∼ 1 − 4. These extreme objects () appear to form stars prodigiously at rates of . Using all-sky Planck and WISE surveys, and wide-area Herschel surveys, we have identified the PASSAGES sample, with some of the rarest hyper-luminous IR galaxies ever discovered. We have found that their globally-averaged star formation surface densities are always sub-Eddington, typically by an order of magnitude. This may suggest that our understanding of how radiation pressure from massive stars disrupts the collapse of molecular clouds (thereby quenching star formation) is flawed—or simply that smaller physical resolutions are necessary. With the aid of lensing, we can now capture the source-plane distribution of star formation at ∼ 100pc scales, letting us identify isolated super-Eddington regions where quenching is occurring.


Introduction
Stellar feedback-the deposition of momentum and energy into the interstellar medium (ISM) by the formation of stars-has come to be accepted as a fundamental process governing the evolution of a galaxy.In particular, this feedback appears to be responsible for preventing gas-rich galaxies from rapidly converting gas into stars at a rate consistent with the free-fall time (see review McKee and Ostriker 2007).This feedback has impacts on global scales of a galaxy, but is ultimately a more local process.The efficiency of star formation is low on galaxy-wide scales; defined as ε ff ≡ τ ff • SFR/M gas -with star-formation rate SFR, total gas mass M gas , and free-fall time τ ff -typical values are of order 1% (e.g.Kennicutt 1998; Krumholz and McKee 2005;Krumholz et al. 2012;Padoan and Nordlund 2011;Federrath and Klessen 2012;Utomo et al. 2018).Yet, the efficiency for giant molecular clouds (GMCs) themselves can be much higher, even approaching unity in some instances (e.g., Murray 2011).These observational constraints to test models of star formation have been difficult to derive, so it remains unclear what physical processes dominate in regulating the collapse of molecular gas.Scoville (2003), Murray et al. (2005), and Thompson et al. (2005) first developed theoretical models by which radiation pressure in dusty environments could become the predominant source of stellar feedback.In particular, where the ISM is optically thick to IR emission (resulting from the reprocessing of UV photons from massive stars), the momentum of the photons against dust grains coupled to the gas becomes an efficient support against the self-gravity of molecular clouds.Such an environment is expected for objects like ultra-luminous infrared galaxies (ULIRGs; Lonsdale et al. 2006)-or more generally the population of dusty star-forming galaxies (DSFGs)-which host extreme dust-enshrouded starbursts.This Eddington-like limit at which radiation pressure exceeds self-gravity is estimated to occur at star formation rate surface densities of around Σ SFR ≡ SFR/(2πR 2 eff ) ∼ 1000 M yr −1 kpc −2 (Andrews and Thompson 2011;Simpson et al. 2015).However, as this limit is for globally-averaged Σ SFR , galaxies forming stars at overall sub-Eddington densities can still host super-Eddington clumps of star formation.Moreover, the measurement of effective radii R eff for high-z (often compact) galaxies requires finer angular resolution than local ULIRGs.For these reasons, the application of gravitational lensing as a "cosmic telescope" is required to make progress in testing models of stellar feedback.

Sample Definition and Data
For this analysis, we employ the lensed DSFGs selected with Planck as part of the PASSAGES sample (Planck All-Sky Survey to Analyze Gravitationally-lensed Extreme Starbursts; Harrington et al. 2016Harrington et al. , 2021;;Berman et al. 2022; lens models from Kamieneski et al. 2023).We also make some comparisons to similar samples defined with the South Pole Telescope (SPT; e.g., Weiss et al. 2013; lens models from Spilker et al. 2016) and Herschel (e.g., Negrello et al. 2010Negrello et al. , 2017;;Wardlow et al. 2013; lens models from Bussmann et al. 2013Bussmann et al. , 2015)).As the Planck-selected objects occupy the upper echelon of apparent infrared luminosities for their redshift range of z ∼ 1 − 3.5 (see also objects identified independently by Cañameras et al. 2015), they offer perhaps the ideal view of extreme starburst events.Moreover, radiative transfer modeling of multi-transition molecular gas lines for the PASSAGES sample by Harrington et al. (2021) suggested that some objects would approach this limit, based on the ratio of their L IR to ISM mass (relative to the threshold of 500 L /M ; see Scoville 2003).Since the total flux-weighted magnifications for extended objects (i.e., starbursts rather than quasars) generally tend not to exceed μ ∼ 10 − 20 (e.g., Hezaveh et al. 2012), their large apparent luminosities were perhaps more easily explained by larger intrinsic luminosities rather than larger magnifications.As such, they are among the best candidates for maximal starbursts forming stars at surface densities close to (or above) the Eddington limit.
Apparent star formation rates were derived from their total infrared luminosities (Harrington et al. 2016;Berman et al. 2022), which captures the UV emission from massive stars that is reprocessed into the IR by dust.As this covers a wide wavelength regime, for which the spatial distribution of continuum can vary, it is not entirely trivial to characterize the size of the dust-emitting region.For our purpose, we use 1.1mm continuum imaging (rest-frame ∼ 250 − 500μm) on the Rayleigh-Jeans tail.Band 6 (1.1 − 1.4 mm) observations were taken with the Atacama Large Millimeter/submillimeter Array as part of Cycle 5 program 2017.1.01214.S (PI: M. Yun), reaching synthesized beam resolutions of 0.4 − 0.8 and sensitivities of 0.07 − 0.33 mJy.The PASSAGES sub-sample that has been observed with ALMA and has a credible strong lens model consists of 13 objects.Further details on the observations are provided in Berman et al. (2022) and Kamieneski et al. (2023).

Results and Discussion
Fortunately, the measurement of Σ SFR is relatively insensitive to statistical and systematic uncertainties in the lens model, as it is connected to surface brightness, which is conserved by gravitational lensing.For example, if a lens model erroneously overestimates a galaxy size by a factor of k relative to the "true" value, then the magnification measured from the ratio of image-plane to source-plane area will be underestimated by a factor of ∼ k 2 ; this then leads to factor of k 2 overestimate of the intrinsic luminosity L IR .As a result of this, both the numerator and denominator of Σ IR = L IR /(2πR 2 eff ) are overestimated by k 2 , and the effect essentially cancels out.Considered slightly differently, Σ IR and Σ SFR could even be measured directly from the image-plane, lensing-uncorrected data, were it not for the challenge that shear distorts the shape of the galaxy and makes the determination of effective radius more complicated.As a more direct comparison with unlensed DSFGs is desirable, we carry out lens modeling in order to estimate the global surface densities of star formation.
Full details of our lens modeling approach for the PASSAGES sample are presented in Kamieneski et al. (2023), but we summarize it in brief here.Only the foreground lensing mass is parameterized (typically with singular isothermal ellipsoid profiles, Kormann et al. 1994), with constraints set by only the locations of multiply-imaged components.These image families are determined using multi-wavelength information from Hubble Space Telescope (HST) at 1.6 μm, Gemini-S Observatory at r and z bands, ALMA at 1.1mm, and the Karl G. Jansky Very Large Array at 6 GHz.Models are optimized through Markov Chain Monte Carlo with lenstool (Kneib et al. 1996;Jullo et al. 2007;Jullo and Kneib 2009), effectively minimizing the RMS of image-plane separations between observed vs. predicted image locations.Magnifications (and associated uncertainties) are estimated by randomly sampling from the multi-dimensional posterior distribution, reconstructing in the source plane by ray-tracing, and computing the ratio of image-plane to source-plane areas.Objects are reconstructed in the source plane without parameterization or regularization, as these extreme starbursts may be clumpy and highly irregular.A similar approach is taken for measuring intrinsic effective radii and uncertainties (see also Section 3 of Kamieneski et al. 2023).
Figure 1 shows the globally-averaged Σ SFR vs. R eff of the dust continuum for three samples of lensed DSFGs, with intrinsic values derived with lens models: PASSAGES (this work, and see Table 3 from Kamieneski et al. 2023), SPT (Spilker et al. 2016), and Herschel (Bussmann et al. 2013).Additionally, unlensed DSFGs from the ALESS sample with resolved size measurements from Hodge et al. (2019) are shown.We find that very few of the objects shown approach the theoretical Eddington limit (with none surpassing it), and the median value of the four samples (Σ SFR ∼ 58 M yr −1 kpc −2 , R eff ∼ 1.6 kpc) is more than an order of magnitude below the threshold.This is in apparent contrast to local ultra-luminous infrared galaxies (ULIRGs; L IR = 10 12−13 L ), which have been found to frequently exceed the Eddington limit, even by over an order of magnitude in some cases (e.g., Barcos-Muñoz et al. 2017;Pereira-Santaella et al. 2021).The canonical and well-studied ULIRG Arp 220, for example, has regions forming stars at surface densities of Σ SFR ∼ 10 4.1±0.1 M yr −1 kpc −2 (Barcos- Muñoz et al. 2015).On the other hand, Song et al. (2022) found predominantly sub-Eddington star formation in a broader sample of local luminous infrared galaxies (LIRGs; L IR = 10 11−12 L ) and ULIRGs.If the extreme hyper-luminous IR galaxies (HyLIRGs; L IR > 10 13 L ) seen at high-z like the PASSAGES objects are in fact sub-Eddington, it might suggest simply that their orderof-magnitude increase in IR luminosity over ULIRGs is accompanied by a comparable ∼ 0.5 dex increase in radius.Indeed, for a very large sample of 1000+ DSFGs covering LIRGs, ULIRGs, and a handful of HyLIRGs from z = 0 − 6, Fujimoto et al. (2017) found a highly significant correlation between R eff and L IR , fit approximately by R eff ∝ L 0.28±0.07IR (albeit with large scatter).The authors speculate that the power-law slope could be set by processes of disk formation, given the similarity in slope for UV-detected starforming galaxies (e.g., α UV = 0.27 ± 0.01, Shibuya et al. 2015).While perhaps a too simplistic picture for complicated high-z systems, it is useful to consider that the Stefan-Boltzmann law suggests that an optically-thick spherical blackbody should radiate as However, as noted by Simpson et al. (2015), "clumpy" distributions of star formation with isolated sub-galactic regions-e.g., star-forming complexes and individual GMCs-forming stars above the Eddington limit could still appear sub-Eddington in its galaxy-wide star formation rate surface density.For this reason, in future work, we intend to directly map Σ SFR at ∼ 100 − 500pc scales with newly-obtained θ ∼ 0.1 resolution ALMA observations of 16 PASSAGES objects at 870μm and 1.1mm.Even with very high-resolution θ ∼ 0.07 ALMA imaging of six ALESS DSFGs, reaching 500pc physical scales, Hodge et al. (2019) found that both the peak and the globally-averaged Σ SFR values remained sub-Eddington.However, these values were actually consistent with those of local ULIRGs from Barcos-Muñoz et al. (2017) with comparable physical sizes, and only ULIRGs with smaller effective radii ( < ∼ 0.2 kpc) were found to be super-Eddington.Hodge et al. suggested this as a sign that even smaller physical scales need to be probed to find super-Eddington star formation in the distant Universe.Indeed, reaching physical scales of ∼ 60 − 160 pc with lensing, Cañameras et al. (2017) found local densities of ∼ 1000 − 2000 M yr −1 kpc −2 .At present, these resolutions are only feasibly accessible with the longest-baseline configurations of facilities like ALMA, and with the aid of gravitational lensing.However, the move to 100 pc scales also brings additional challenges.For example, measurement of the properties of small star-forming clumps may be more affected by gravitational microlensing (e.g., Fian et al. 2021), or by millilensing from dark, low-mass substructures in the lensing galaxy or anywhere along the line-of-sight (e.g., Mao and Schneider 1998;Vegetti et al. 2012Vegetti et al. , 2014;;Inoue 2016;Nierenberg et al. 2017Nierenberg et al. , 2020)).Such effects are negligible at coarser resolutions, but they should be incorporated into uncertainties when approaching GMC scales.Fortunately, the conservation of surface brightness from lensing may again work to our advantage when mapping Σ SFR , whereas the measurement of absolute properties like the mass and luminosity of clumps is made more difficult.
Novel selection criteria and new facilities will continue to increase this sample of lensed submillimeter-bright objects in the coming decade.As one notable example, the TolTEC camera of the Large Millimeter Telescope will have unprecedented sensitive mapping capabilities at 1.1mm, 1.4mm, and 2.1mm simultaneously (Wilson et al. 2020).This longer-wavelength regime could particularly help in accessing higher-redshift (z > ∼ 4) lensed DSFGs (see, e.g., Casey et al. 2021).With this growing sample of lensed, intrinsically extreme starbursts-alongside the anticipated large ensemble of more than 10 5 optically-selected lenses with new mapping facilities like Euclid (e.g.Laureijs et al. 2011;Serjeant 2014)-new frontiers are being opened up for testing theoretical frameworks for star formation and the dynamics of the interstellar medium.Strong gravitational lensing will continue to be a crucial tool in helping to bridge the gap in the physical resolutions that we can reach for the local vs.high-redshift Universe.

Figure 1 .
Figure 1.The galaxy-averaged star formation rate surface density ΣSFR vs. dust continuum effective radius R eff for members of the PASSAGES sample (Kamieneski et al. 2023), as measured in the model-reconstructed source plane.Also shown are the source-plane values for samples of lensed DSFGs identified with SPT (Spilker et al. 2016) and Herschel (Bussmann et al. 2013), alongside unlensed DSFGs from the ALESS sample (Hodge et al. 2019).The Eddington limit is typically computed at values of ΣSFR ∼ 650 − 1000 M yr −1 kpc −2 (depending on choice of gas mass fraction and dust-to-gas ratio; Andrews and Thompson 2011; Simpson et al. 2015; Hodge et al. 2019), which is shown as a red shaded region.Diagonal dotted, dashed, and dashdotted lines indicate the loci for IR luminosities of 10 11 L (LIRGs), 10 12 L (ULIRGs), and 10 13 L (hyper-luminous infrared galaxies or HyLIRGs).