Hidden Giants in JWST's PEARLS: An Ultramassive z = 4.26 Submillimeter Galaxy that Is Invisible to HST

We present a multiwavelength analysis using the Submillimeter Array (SMA), James Clerk Maxwell Telescope, NOEMA, JWST, the Hubble Space Telescope (HST), and the Spitzer Space Telescope of two dusty strongly star-forming galaxies, 850.1 and 850.2, seen through the massive cluster lens A 1489. These SMA-located sources both lie at z = 4.26 and have bright dust continuum emission, but 850.2 is a UV-detected Lyman-break galaxy, while 850.1 is undetected at ≲ 2 μm, even with deep JWST/NIRCam observations. We investigate their stellar, interstellar medium, and dynamical properties, including a pixel-level spectral energy distribution analysis to derive subkiloparsec-resolution stellar-mass and A V maps. We find that 850.1 is one of the most massive and highly obscured, A V ∼ 5, galaxies known at z > 4 with M * ∼1011.8 M ⊙ (likely forming at z > 6), and 850.2 is one of the least massive and least obscured, A V ∼ 1, members of the z > 4 dusty star-forming population. The diversity of these two dust-mass-selected galaxies illustrates the incompleteness of galaxy surveys at z ≳ 3–4 based on imaging at ≲ 2 μm, the longest wavelengths feasible from HST or the ground. The resolved mass map of 850.1 shows a compact stellar-mass distribution, Remass ∼1 kpc, but its expected evolution means that it matches both the properties of massive, quiescent galaxies at z ∼ 1.5 and ultramassive early-type galaxies at z ∼ 0. We suggest that 850.1 is the central galaxy of a group in which 850.2 is a satellite that will likely merge in the near future. The stellar morphology of 850.1 shows arms and a linear bar feature that we link to the active dynamical environment it resides within.


Introduction
Obscuration by dust has been known to be a complicating factor in understanding the properties of high-redshift galaxies for the last ∼30 yr.The first studies of z  3 star-forming galaxies employed samples selected through their rest-frame ultraviolet (UV) emission (Steidel & Hamilton 1993).However, the subsequent identification in the submillimeter wave band of highly dust-obscured galaxies at similar redshifts (e.g., Ivison et al. 1998) demonstrated a much wider diversity in the properties of galaxies at high redshifts than just the subset detectable in the rest-frame UV (Ivison et al. 2000).The potential for entire classes of dusty, massive galaxies (detectable in the (sub)millimeter or the mid-infrared with Spitzer/IRAC and variously termed optically/near-infrared faint/dark or less precisely "HST-dark") to be missed from UV, visible-and near-infrared-selected samples, has been a concern for studies attempting to construct stellar-mass-limited galaxy surveys at z > 1 (e.g., Dey et al. 1999;Caputi et al. 2012;Simpson et al. 2014;Wang et al. 2019).With the launch of JWST, there has been a renewed appreciation of the significant role of dust in defining the observed visible to nearinfrared properties of high-redshift galaxy populations (e.g., Barrufet et al. 2023;Bisigello et al. 2023;Kokorev et al. 2023;Magnelli et al. 2023;Pérez-González et al. 2023).
The potential significance of missing optically undetected populations was underlined by the analysis of large samples of dust-obscured galaxies.These populations are generally optically faint, but the majority of star formation in the Universe at z  4-5 was estimated to occur in these obscured systems (e.g., Bouwens et al. 2020;Dudzevičiūtė et al. 2020).It was the rapid evolution in this dust-obscured activity that created the peak in the star formation rate (SFR) density at z ∼ 2 (e.g., Magnelli et al. 2013).Discovering the physical processes driving this dust-obscured activity is thus crucial for understanding galaxy evolution.The most massive and obscured examples of this high-redshift, dust-rich population are known as submillimeter galaxies (SMGs).They are found in surveys at 0.85-1.25 mm that select galaxies based on cool dust mass out to z ∼ 6 (e.g., Barger et al. 1998;Eales et al. 1999;Smail et al. 2002;Borys et al. 2003;Webb et al. 2003;Greve et al. 2004;Laurent et al. 2005;Coppin et al. 2006;Bertoldi et al. 2007;Scott et al. 2008;Weiß et al. 2009;Aretxaga et al. 2011;Hatsukade et al. 2011;Lindner et al. 2011;Casey et al. 2013;Umehata et al. 2014;Hsu et al. 2016;Cowie et al. 2017).These galaxies are proposed to evolve into the cores of massive galaxies at z ∼ 0 (e.g., Eales et al. 1999;Dudzevičiūtė et al. 2020;Amvrosiadis et al. 2023).
Current studies of the properties of dust-obscured galaxies are hampered by both their complex dust obscuration and high redshift.In particular, their rest-frame UV to visible morphologies (available from HST) may not provide true indicators of their structures, necessary to determine the role of interactions and mergers in driving their evolution (e.g., Gillman et al. 2023;Huang et al. 2023).JWST promises to revolutionize this field by providing high-spatial-resolution imaging in the restframe near-infrared wave band, necessary to map the internal variations in the dust reddening and star formation, as well as isolating any contributions from reddened active galactic nucleus (AGN) emission, which may have confounded previous attempts to estimate the stellar masses of some of these systems (e.g., Hainline et al. 2011).The analysis of resolved JWST rest-frame near-infrared imaging will thus hopefully deliver much more reliable stellar masses (Sorba & Sawicki 2018).In particular, ∼0 1 FWHM NIRCam imaging will yield more robust stellar morphologies on kiloparsec scales that may indicate what drives the growth of the population of massive, and gas-and metal-rich galaxies responsible for the increasing obscuration of star formation at z  4-5.
To exploit the unique JWST imaging from the Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS) Guaranteed Time program (Windhorst et al. 2023), we collated the submillimeter imaging from the SCUBA-2 camera on the James Clerk Maxwell Telescope (JCMT) of all of the PEARLS fields visible from Maunakea, using the JCMT archive at CADC.The analysis of a 30-minute SCUBA-2 exposure of the cluster A 1489 (z = 0.35) revealed two bright and previously unknown 850 μm sources, SMM J121223.0 +273351 and SMM J121220.4+273410(hereafter, 850.1 and 850.2) with S 850μm ∼60 and 20 mJy, respectively.This snapshot observation was so shallow that it would not normally be expected to yield any significant detections, so the presence of two bright sources was a surprise.The total area of the archival SCUBA-2 data that was searched meant there was only a ∼0.1%-1% chance of finding a pair of such bright sources (e.g., Garratt et al. 2023), where the lower probability assumes the two sources are not physically associated.Initial analysis of these sources suggested that they were likely to be behind the cluster and hence high-redshift, dusty star-forming galaxies gravitationally magnified by the massive foreground cluster (Zitrin et al. 2020).
Gravitational lensing has been used as a tool to aid the study of dusty galaxies at high redshifts since their discovery.The first example of a high-redshift dust-obscured galaxy, FSC 10214 at z = 2.286 (Rowan-Robinson et al. 1991), was uncovered by virtue of the gravitational magnification caused by a serendipitously positioned foreground galaxy.Lensing by clusters and individual galaxies has been exploited ever since, to study high-redshift dusty star-forming galaxies (e.g., Smail et al. 1997;Swinbank et al. 2010;Ciesla et al. 2020;Fujimoto et al. 2023), as lensing both boosts the apparent fluxes and magnifies sources across all wavelengths, making it easier to study these otherwise very faint galaxies.
Follow-up of the two A 1489 submillimeter sources with the Submillimeter Array (SMA), NOEMA, and JCMT gave precise locations and redshifts for the two bright sources and constrained their far-infrared emission.These observations showed that the two sources were indeed very dusty, starforming galaxies lying ∼100 kpc apart in a structure at z ∼ 4.26.One counterpart was undetected in the HST and Spitzer 3.6/4.5 μm imaging of this field, while the other corresponds to a UV-bright Lyman-break galaxy.They are thus examples of the two populations that dominate the SFR density in the transition from primarily unobscured at z  5 to predominantly obscured at z  4-5, with the benefit that the lens magnification will help JWST resolve their internal structures.
The structure of this paper is as follows: Section 2 describes our new and archival multiwavelength observations, their reduction and analysis, including spatially resolved spectral energy distribution (SED) fitting.Section 3 presents our results and discusses these, and our conclusions are given in Section 4. We assume a cosmology with Ω M = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 .In this cosmology, the luminosity distance to z = 4.26 is 39.4 Gpc, and 1″ corresponds to 6.9 kpc before correcting for lens magnification.All quoted magnitudes are on the AB system, and uncertainties on median values were derived from bootstrap resampling.The quoted stellar properties assume a Chabrier (2003) initial mass function (IMF).

JCMT
The initial 30-minute SCUBA-2 observation of A 1489 came from program M15AI29, taken in very good weather, τ 225GHz = 0.03, on 2015 May 6.The 850 μm map showed two significantly detected point sources, one of which was exceptionally bright.The presence of these unusual sources in this JWST survey field motivated us to obtain higherresolution continuum observations at 880 μm with SMA as well as deeper continuum observations at 450 and 850 μm with SCUBA-2, and spectroscopic observations covering the 3 mm band from NOEMA.
To improve the depth of the submillimeter maps, the cluster was observed for program M23AP006 with SCUBA-2 simultaneously at 450 and 850 μm (at resolutions of 7 5 and 14 5 FWHM, respectively) in very good weather conditions τ 225GHz = 0.03-0.04 on the nights of 2023 February 15, March 21, and May 14.The total exposure time was 3.5 hr using a standard constant-velocity daisy mapping pattern, and the 30 minutes of exposure taken in similar conditions from M15AI29 was added to give a total exposure time of 4 hr in each band.
Individual maps for each 30-minute observation were reduced using the Dynamic Interactive Map-Maker (DIMM) tool of the Sub-Millimetre User Reduction Facility (SMURF; Chapin et al. 2013) with the blank field configuration in order to detect point sources within the maps.The bright source 850.1 was well detected in each 30-minute exposure, so a shiftand-add method (Ivison 2006) was applied to the data, coaligning the positions at 850 μm to the median position from the M23AP006 observations.The same offsets were applied to the 450 μm maps, and this modestly improved the signal-to-noise for sources in both bands.
The individual maps were calibrated using a flux conversion factor of FCF 450μm = 491 Jy beam −1 pW −1 and FCF 850μm = 537 Jy beam −1 pW −1 (Dempsey et al. 2013) and then combined using inverse-variance weighting to create a final map at each wavelength.To improve point-source detection, the resulting 450 μm and 850μm maps were match-filtered with their respective beams: 7 5 and 14 5 FWHM Gaussians.The maps were cropped to radii of 4 0, where the noise was low and more uniform, providing coverage of all of the primary JWST field.The 1σ sensitivity in the region of 850.1 and 850.2 was 1.5 mJy at 850 μm and 21 mJy at 450 μm after matched filtering.The resulting 450 and 850 μm maps are shown in Figure 1.
The new maps identified a third bright submillimeter source to the northwest of the cluster, SMM J121213.4+273517(850.3), with S 850μm = 16.5 mJy, but undetected at 450 μm (Table 1).Unfortunately this source is ∼2′ from 850.2 and lies outside the coverage of the existing HST, JWST, SMA, and NOEMA observations, and there was no obvious counterpart in the Spitzer imaging.This source is not discussed further, beyond noting that it is likely to be intrinsically luminous as the lensing magnification from the cluster will be modest (1.5 ×) and the 450/850 μm flux ratio suggests this source may lie at z > 3; thus, it could be another member of the z = 4.26 structure.
Positions and flux densities for all three 850 μm sources were derived by fitting Gaussian beams to the 850 and 450 μm maps.The resulting fluxes are given in Table 1, and these were deboosted at 450 μm following Geach et al. (2017) to reflect the lower signal to noise.

SMA
SMA observed the JCMT-determined positions for 850.1 and 850.2 during a single transit between UT 6.1 and 16.3 on 2022 March 14 for the DDT Project 2021B-S072 (PI: S. Willner).At the time, the SMA was operating with all eight antennas in the EXT configuration (maximum baselines of ∼220 m).The array receiving system was tuned to an LO frequency of 346 GHz (giving an effective band center of ∼880 μm), providing LSB coverage from 330-342 GHz and USB coverage from 350-362 GHz in each of two orthogonally polarized receivers, resulting in 24 GHz of dual-polarization continuum bandwidth.The weather was good, with the precipitable water vapor column ranging from 1.5-2 mm and stable, low atmospheric phase variation ("seeing").
Observations were made in a short repeated cycle, measuring the nearby phase and amplitude gain calibrator J 1159+292 for 1.5 minutes, then the positions of 850.1 and 850.2 for 4 minutes each.This resulted in high-quality gain calibration transfer to each target.In total, 850.1 was observed for 3.75 hr and 850.2 for 3.85 hr.Passband calibration was determined from observations of J 0854+201, J 1159+292, and BL Lac, and the absolute flux scale was derived from observations of Ceres and MWC 349A, taken immediately before and after the targets.
For each source, the resolution of the synthesized beam was 0 85× 0 63-around 10 × better than from SCUBA-2-and for sources with lensing magnification of μ ∼ 4-6, the circularized beam corresponds to 2.0-2.5 kpc FWHM in the source plane.The rms of the maps was σ 880μm ∼0.8 mJy beam −1 , and both 850.1 and 850.2 were well detected in these subarcsecond resolution maps, accurately locating the submillimeter sources.Table 1 reports their precise positions.
Note that 850.1 was detected with a peak flux density around 19 mJy.However, this source was clearly resolved in both dimensions (Figure 2).Similarly, 850.2 had a peak flux density of ∼9 mJy, although it also appeared resolved in the east−west direction.The vector-average amplitude versus (u, v) distance plots for both sources showed increasing amplitudes on short baselines.This demonstrated that both sources were resolved.Fitting the visibility data and extrapolating these to zero spacing indicated total flux densities of S 880μm = 42.3 ± 1.6 mJy for 850.1 and S 880μm = 22.5 ± 1.6 mJy for 850.2.Table 1 reports these values, along with the other longwavelength measurements.The estimated flux density for 850.2 was in good agreement with that measured at 850 μm by SCUBA-2, but the SMA flux for 850.1 was ∼20% lower, a ∼6.5σ difference.Half of this may be accounted for the difference in the effective band centers of the two instruments, but the remainder suggests that there may be further extended continuum, or [CII] emission (Smail et al. 2011), around 850.1 that was resolved out with the subarcsecond SMA beam.
At the redshifts of 850.1 and 850.2 determined from their NOEMA 3 mm spectra (see Section 2.3), the [CII] 1900.54GHz line falls in the spectral coverage of the SMA observations.The 1σ sensitivity of these data after binning to ∼120 km s −1 (143 MHz) channels was ∼25 mJy.However, this sensitivity would likely be insufficient to detect [CII], assuming typical [CII] to far-infrared luminosity ratios for the sources.Indeed, a search for the line in the sources did not show any significant line emission associated with the continuum sources at the relevant frequencies.

NOEMA
NOEMA observations were obtained as part of the S22CX program in several partial tracks during 2022 July.An earlier Gemini-N GMOS program (GN-2020A-Q-903, PI: Zitrin) had yielded an optical spectroscopic redshift for the SMA-located HST counterpart to 850.2 from template-fitting of z = 4.267 (Nonino & Zitrin, private communication).This allowed the LO tunings to be chosen to cover the expected frequencies of the CO(4-3), CO(5-4), and [C I] lines of this source, in the expectation that 850.1 was likely to lie at high redshift and potentially the same redshift as 850.2.Hence, two setups were used in Band 1 (3 mm) with LO frequencies of 93.5 GHz and 101.244GHz.The wide bandwidth of the PolyFiX correlator means that these two setups covered a contiguous frequency range of ∼31 GHz from 81.88-112.86GHz.All tracks were observed in the compact D-configuration (maximum baseline length ∼175 m) with 10 antennas in the array for five partial tracks, nine antennas for three partial tracks, and eight antennas for one partial track.
Calibration was performed using the Grenoble Image and Line Data Analysis Software (GILDAS) package.Either MWC 349 or LkHa 101 was observed for use as the flux calibrator; in cases where one of these sources was either not observed or the data were of poor quality, the gain calibrator 1156+295 was used as the flux calibrator, bootstrapping its flux from other tracks with reliable observations of MWC 349 or LkHa101.Bandpass calibration was obtained from observations of 3C 273, 3C 279, and BL Lac.After flagging for bad visibilities, the effective on-source time for a 10-antenna array for each of 850.1 and 850.2 was 1.6 hr from the lowerfrequency setup and 1.8 hr from the higher-frequency setup.Finally, dirty cubes were produced with MAPPING (part of the GILDAS package) using natural weighting from uv-tables resampled to 20 MHz.The synthesized beam sizes for the cubes ranged between 3 8 and 6 1 FWHM. is shown as blue, F200W+F277W is green, and F356W+F444W is red).Signal-to-noise contours are shown, representing the SCUBA-2 850 μm emission (starting at 5σ, in 5σ increments, 14 5 FWHM beam), the SCUBA-2 450 μm emission (3σ and 4σ, 7 5 FWHM beam), and the NOEMA 3 mm continuum emission (10σ, ∼5″ FWHM beam).Note that 850.1 and 850.2 are the two bright sources visible at 450 μm, 850 μm, and 3 mm to the northeast of the cluster core, which is identified by several bright elliptical galaxies, as well as associated strong lensing features.The counterparts to these sources are unambiguously identified using the higherresolution SMA 880 μm observations shown in Figure 2.
The reduction of the NOEMA observations showed 3 mm continuum detections of both sources with positions that agree within the relative astrometric uncertainty with those from SMA (Figure 1).Spectra extracted at the positions of the two galaxies yielded three lines in each spectrum (Figure 3).The frequencies of the brighter two lines in each source were 87.65 GHz and 109.56GHz (Figure 3, Table 2) corresponding to CO(4-3) and CO(5-4), at 461.041 GHz and 576.268GHz in the rest frame, respectively.The third, fainter line in each source (at ∼93.6 GHz) was [C I](1-0) at 492.161 GHz.These low-resolution observations in NOEMA's compact configuration showed no significant evidence for extended emission in either the continuum or lines.
Gaussian fits to all three emission lines in each source yielded redshifts of z = 4.2599 ± 0.0001 for 850.1 and z = 4.2593 ± 0.0003 for 850.2 (Figure 3).The latter was in reasonable agreement with the optical measurement for 850.2.The Gaussian fits were then used to derive line centers and widths, which were used to integrate the emission to derive the intensity (after continuum subtraction) within a window of ± FWHM around the peak.For the weak [C I] line in 850.2 the FWHM of the CO lines was used to determine the integration window.From these, the line intensity, center, and width (expressed as an equivalent Gaussian FWHM), were measured along with their uncertainties.The second moment of the line yields FWHM-equivalent estimates that were ∼20% smaller than the Gaussian fitted values.Note that 850.2 shows evidence of a double-peak profile in both the CO(4-3) and CO(5-4) lines (Figure 3), and so for that source, the velocity difference between the two peaks was instead used as the best estimate of the kinematics of the source.The line and continuum properties of both sources are given in Tables 1,  2, and 3.The CO line luminosities were used to estimate line luminosity ratios between transitions from J up = i to J up = j, r ij , following Solomon & Vanden Bout (2005).

HST, JWST, and SST
The details of the observations, reduction, and analysis of the HST imaging of A 1489 (PID: 15959) were reported in Zitrin et al. (2020).The observations were obtained in 2020 March and comprised a total of five orbits, including one orbit (∼1.9 ks) in each of the F435W, F606W, and F814W filters with the Advanced Camera for Surveys (ACS) and shorter exposures with WFC3 in the F105W, F125W, F140W, and F160W filters.The WFC3 data were not used in our analysis as they were shallower than the equivalent JWST imaging and also do not cover both targets.As described by Zitrin et al. (2020), mosaicked images were produced for all of the HST filters using the calibrated exposures at a pixel scale of 0 03 pixel −1 (matching that used for NIRCam) on an astrometric grid aligned to Gaia DR2.See Zitrin et al. (2020) for more details of the reduction and analysis.
JWST observed A 1489 as part of GTO program 1176 (PEARLS; Windhorst et al. 2023) targeting moderate-redshift clusters and blank fields to study galaxy evolution.A 1489 was imaged in 2023 January with NIRCam filters F090W, F150W, F200W, F277W, F356W, and F444W.Observations of F090W and F444W were taken with the MEDIUM8 readout with six groups per integration and four integrations for a total exposure time of 2.5 ks.F150W and F356W were observed with an SHALLOW4 readout with nine groups per integration and four integrations for a total exposure of 1.9 ks, while F200W and F277W were observed with an SHALLOW4 readout with 10 groups per integration and four integrations for a total exposure of 2.1 ks.
The NIRCam images were initially processed using the STScI JWST data reduction pipeline version 1.8.4 (Bushouse et al. 2022) with calibration reference files corresponding to the JWST_1039.PMAP context.Further processing was then applied to remove artifacts in the images.These steps were described in detail by Windhorst et al. (2023).In brief, "snowballs" (resulting from energetic cosmic-ray hits) were identified and masked; then "wisps" caused by stray light from the secondarymirror support were corrected using filter-specific templates.Next, striping due to amplifier differences was corrected by matching the background level in adjacent columns using the PROFOUND code (Robotham et al. 2018).Finally, detectorlevel offsets in the short-wavelength camera were corrected by constructing a "super-sky" with PROFOUND and scaling each detector to that sky level.The final stacked images were created at a pixel scale 0 03 and were aligned with the existing HST images on the same scale.
Figure 1 shows a JWST and HST view of the wider environment of 850.1 and 850.2.More detailed zooms are shown in Figure 2, including the SMA 880 μm maps that precisely locate the counterparts to the submillimeter sources.Figure 4 illustrates the variation in the brightness of the two counterparts with wavelength across the nine JWST and HST bands.Sources were identified in a sequence starting with "A" for distinct sources, "B" for potential subcomponents, and "C" for likely emission-line features.These are labeled on the F090W images.
To extend the wavelength coverage of these dusty, highredshift galaxies, the pipeline reductions of the Spitzer IRAC and MIPS 24 μm images of A 1489 were retrieved from the SST (Spitzer Space Telescope) archive.This was undertaken prior to both the launch of JWST and the SMA observations, finding potential 8 μm bright counterparts within the SCUBA-2 error circles for both 850.1 and 850.2, that turned out to be the correct identifications.These data were included in our SED analysis to extend the wavelength coverage in the gap between JWST and SCUBA-2.Neither source was detected at 24 μm with MIPS.This was unsurprising given their high redshifts and the modest depth of the data, and so a 3σ limit was adopted at that wavelength (Table 1).The emission in the IRAC bands was blended with the nearby galaxies for both 850.1 and 850.2, and therefore the flux estimates were based on small 2 4-3 6 diameter apertures.However, the SMGs were increasingly dominant at the longer wavelengths, and 850.2 was effectively uncontaminated by 8 μm.For both sources, excluding the IRAC and MIPS data points, or adopting limits instead, did not alter the solutions for the best-fit SEDs obtained below.
To measure integrated photometry for 850.1, 850.2, and the other sources in close proximity to them in the JWST and HST data, large circular apertures were used (except for 850.1) with sizes listed in Table 4.These provided close-to-total fluxes for most of the sources.These measurements used the photometric zero-points from Windhorst et al. (2023) for NIRCam and from +F277W as green, and F356W+F444W as red.Panels (c) and (e) show ¢¢5 5″ regions using F277W, F356W, and F444W as blue, green, and red and the same filter combination is used in panels (d) and (f), which show views of 850.1 and 850.2 "de-sheared" to correct for the lensing magnification (the scale bar indicates 1 kpc).Panels (a) and (b) also show the SMA 880 μm contours (starting from 2σ in 7σ increments), these unambiguously identify the two submillimeter sources.Note that 850.1 corresponds to an extremely red, highly structured source lying close to a bright, foreground disk galaxy with several fainter sources in close proximity.In contrast, 850.2 is more isolated and appears to be a bluer, disk-like system.The various sources and features discussed in the text are identified on each panel: distinct sources have labels starting with "A," "B" for potential subcomponents, and "C" for likely emission-line features.The colors and photometric redshifts suggest that "knot" B1 near 850.1 is likely to lie at z < 4.26, while B2 is probably a close/interacting companion seen in projection to 850.1 or a less-obscured component within the galaxy.The photometric-redshift analysis also suggests that A1, A4, and A5 are likely to be cluster members at z ∼ 0.35, while A2 and A3 are probably behind the cluster, but foreground to 850.1 and 850.2.
the ACS pipeline (as given in the FITS file headers) for the ACS imaging.Point-spread function (PSF) matching employed model PSFs created with WEBBPSF v1.1.0or TINYTIM for HST/ACS, utilizing the ODT files closest to the observation date, and convolved to match the F444W PSF.Due to the large aperture sizes, PSF corrections for 850.1 and 850.2 were modest in all cases.Sky corrections were derived from large annuli around the targets.An additional correction was applied for the contribution from the disk galaxy, A4, to those sources around 850.1.This was modeled in each filter assuming rotational symmetry and rotating the images by 180°around the center of A4, then smoothing slightly, before subtracting this image from the original.The resulting photometry is listed in Table 4 with 3σ limits for nondetections.
To determine the sizes and profiles of the two sources, the GALFIT code was applied to the different JWST wave bands assuming single Sérsic profiles and correcting for the varying resolution between the bands using an empirical PSF taken from stars in the images.This process was challenging in the crowded environment of 850.1, especially given the fine structures visible in the surrounding galaxies as a result of the superlative resolution and depth of the data.However, aided by the fact that the bulk of 850.1 was undetected in the F090W image, a model was constructed in this band of the surrounding galaxies and used as the basis for constrained fitting to the redder passbands, where a further component was included to represent 850.1.The structural parameters of both sources from this analysis are reported in Tables 2 and 5.

Integrated SED Modeling
The SED modeling of the two submillimeter sources was iterated with the aim of understanding the nature of the various subcomponents in these systems and so derive robust physical properties for the two galaxies.The MAGPHYS code (da Cunha et al. 2015;Battisti et al. 2019) was adopted for this analysis as it uses an approximate energy-balance calculation to ensure self-consistent modeling of the dust absorption in the UV-tooptical and the reprocessed far-infrared-to-submillimeter emission that is particularly important in highly obscured systems.Applying MAGPHYS to the sources in this study also removes systematic uncertainties when comparing to the large samples of SMGs that have been previously analyzed with this code (e.g., da Cunha et al. 2015;Smolčić et al. 2015;Miettinen et al. 2017;Simpson et al. 2020;Dudzevičiūtė et al. 2021).MAGPHYS has also been extensively validated using tests on simulated galaxies (Dudzevičiūtė et al. 2020).
The analysis here employed both the photometric-redshift version of the MAGPHYS code (Battisti et al. 2019) as well as the high-redshift version (da Cunha et al. 2015) to derive parameters at the known CO redshifts, both using the standard priors (including a star formation history that rises at early times, before declining exponentially, with superimposed random bursts).The codes were applied to the integrated photometry of the sources (including potentially unrelated subcomponents), to the integrated photometry of the individual subcomponents, and also to the resolved "pixel-level" photometry (described in Section 2.6).In principle, the latter resolved SED fitting should be most robust, but the modest resolution of some of the longer-wavelength observations (Spitzer IRAC, SCUBA-2, and NOEMA) means that there were fewer constraints available on the form of the SEDs.For this reason, the resolved SED fitting is presented as a proof of concept, and for most of the quantitative analysis, the fitting to the integrated photometry was used to derive the physical properties.
The photometric analysis of 850.2 was straightforward owing to the relative isolation of the galaxy from other sources in the field, as well as the relatively good agreement between the redshifts derived from the CO and [C I] emission lines and the archival rest-frame UV spectroscopy.The latter confirmed that the emission seen in both submillimeter and the rest-frame UV-optical all arises from the same physical system, although not necessarily precisely the same components within this system.
In contrast, 850.1 has a complex environment, and the analysis of this system required iteration.Based on the source morphology and its variation with wavelength shown in Figure 4, as well as the results from the resolved SED modeling discussed in Section 2.6, the emission around 850.1 was divided into three components: two blue "knots" B1 and B2, and the bulk of the remainder of the source, 850.1, which shows much redder colors and corresponds to the peak of the SMA map (Figure 2).The analysis started by estimating likely redshifts for these three components as well as the other sources seen close to 850.1 and 850.2 that lacked spectroscopic redshifts (those marked in Figure 2 and listed in Table 4 and 5), with the photometric-redshift version of MAGPHYS (da Cunha  et al. 2015; Battisti et al. 2019).These photometric redshifts are reported in Table 5.These were then used to assess the level of agreement of the spectroscopic redshift measurements for 850.1 and 850.2 with their photometric properties.
For 850.2, the photometric redshift, z ph = 4.28 - + 0.10 0.08 , driven by the strong Lyman break in the broadband SED (Figures 4  and 5) was in excellent agreement with the spectroscopic  1 (left) and 850.2 (right), showing the CO(4-3), [C I](1-0), and CO(5-4) emission lines.The panels are centered on the expected frequencies of the lines at the adopted redshifts for each source from Table 2 and a model fit comprising a Gaussian and a constant continuum level is shown for each line.Both sources have relatively weak [C I] emission compared to their CO(4-3) and CO(5-4).A single Gaussian provided an adequate description of the emission for the lines in 850.1, but all three lines in 850.2 display a dip at the systemic velocity, suggesting a double-peaked profile, with the higher-frequency CO peak close to the redshift derived from the rest-frame UV spectroscopic features.The spectra for 850.2 have been rebinned to 40 MHz per channel, except for the [C I] line, which has been rebinned to 100 MHz per channel (and as a result is shown over a slightly wider frequency window to illustrate the off-line noise).

Table 2
Kinematic Properties of Sources

Notes.
a Axis ratio measured from the F444W isophote shape and corrected for lensing shear.b To estimate the dynamical mass, we used the velocity difference between the two components of the CO(4-3) and CO(5-4) lines, 330 ± 70 km s −1 , rather than the FWHM.redshift z = 4.26 derived from the multiple lines seen in the NOEMA spectrum (Figure 3), and also the Gemini-N/GMOS spectrum (Nonino & Zitrin, private communication).The latter spectrum showed a pronounced decline in flux blueward of ∼1215 Å, no Lyman-α emission but a series of strong absorption features including SiII λ1260, OI+SiII λ1303 (blend), CII λ1334, SiII λ1526 and CIV λ1548, FeII λ1608 and AlII λ1670.The SiII λ1260 and CIV λ1548 lines displayed P-Cygni-like profiles.These low-ionization interstellar absorption features are commonly seen in Lyman-break galaxies (Shapley et al. 2003).In 850.2, their measured line centers displayed a median blueshift of ∼ −190 km s −1 relative to the CO-derived systemic redshift.This velocity offset may be a signature of a wind, or it indicates that the UV-bright element of the system corresponds to the higher-frequency peak in the CO spectrum of the source (Figure 3), with the redshifted CO component then potentially arising in the redder part of the system.
A similar photometric analysis for 850.1 used a photometric aperture that includes the emission from the two knots B1 and B2 (Figures 2 and 4, Table 4).This resulted in a photometric redshift that was significantly lower than the CO redshift: z ph = 2.63 - + 0.23 0.46 (Table 5).Figure 4 shows clearly that B1 and B2 are considerably bluer than the bulk of 850.1, both dominating the emission shortward of ∼2 μm, and B1 is visible all the way down to F435W, where the wavelength coverage of the filter falls almost completely blueward of 912 Å in the rest frame at z = 4.26.Photometric redshifts measured for the three components individually yielded estimates of z ph = 1.12 - + 0.12 0.11 and z ph = 4.17 - + 3.28 0.16 for B1 and B2, respectively, while the very red 850.1 gave z ph = 4.41 - + 1.95 0.77 in reasonable agreement with the CO redshift (Table 5).
The photometric-redshift estimate for B1 indicated that the source was unlikely to be associated with either 850.1 or the foreground cluster.Thus, while the alignment of B1 with features in 850.1 is very suggestive that they are related, B1ʼs detection in the F435W filter is hard to reconcile with it lying at z = 4.26, and it is more likely to be a superimposed lowerredshift source.The situation is more complicated for B2, where the estimated photometric redshift suggests it could be associated with 850.1 at z = 4.26, but the uncertainty allows for much lower redshift solutions as well.We therefore retain B2 as a potential UV-bright component within 850.1 or a close/interacting companion and discuss how this influences the inferred properties of the source.
Having confirmed that the measured photometry for 850.1, B2, and 850.2 were all consistent with the CO-derived redshifts, more representative estimates of their physical properties were derived by running the high-redshift version of MAGPHYS (da Cunha et al. 2015) with the redshifts fixed at z = 4.26 for all three components.The main parameters from these fits are reported in Table 5, and the best-fit SEDs are plotted in Figure 5.For 850.1 and B2, the two best-fit SEDs were summed and labeled "850.1+B2fits," and this is plotted compared to the integrated photometry of the system, although 850.1 dominates the emission longward of ∼3 μm.In addition, Figure 5(c) show the results of fitting solely to 850.1, illustrating the very red colors of this source.
The fit to the combined photometry of 850.1 and B2 using a single MAGPHYS model and labeled "850.1+B2unresolved" is also shown in Figure 5(c).This was much less successful at reproducing the optical/near-infrared SED of the combined source.This resulted in a statistically much poorer fit to the observed photometry that predicted a much lower mass, SFR, A V , and age (see Table 5) but a comparably high dust mass to the preferred sum of the individual fits to 850.1 and B2.These parameters would imply a galaxy more similar to 850.2 in its physical properties, but we believe that this poor fit was erroneous due to the distinct low and high extinction regions that were being modeled with a single dust extinction term.Similar model biases may also be responsible for some apparently low-mass SMGs identified from SED fitting of systems with a wide range in A V (e.g., Pope et al. 2023).This erroneous fit also produces an apparent power-law-like excess in the observed fluxes longward of 3.6 μm compared to the model.This is reminiscent of similar features seen in the SEDs of ∼10% of submillimeter sources analyzed by Hainline et al. (2011), who attributed it to hot dust emission from dustobscured AGNs.If the origin of this excess is similar in those sources to 850.1, then an obscured AGN may not be the cause as there is no such bright compact, red point source seen in 850.1.Instead the red mid-infrared excess may be a result of fitting a single SED model to a source with a mixture of internal dust obscuration, including both weakly and highly reddened regions.Finally, a non-energy-balance code, FAST++ (Kriek et al. 2009;Schreiber et al. 2018), was used to fit the SEDs of 850.1, B2, and 850.2 covered by HST, JWST, and IRAC (so excluding the mid-, far-infrared, and submillimeter).This modeling adopted the Bruzual & Charlot (2003) stellar tracks, a Chabrier (2003) IMF, Calzetti et al. (2000) dust law, and exponentially declining star formation histories.As the sources are expected to be dusty, the A V was allowed to range between 0 and 6, and metallicities of 0.008, 0.02, and 0.05 were used for the stellar tracks.We fixed the redshifts to z = 4.26 for all three sources.For 850.1, the resulting best-fit parameters were: μM * = 5. 20 M e yr −1 , and 1 .The sum of the FAST++ fits to 850.1+B2 yielded a stellar mass that is dominated by 850.1 but is ∼4 × lower than that from MAGPHYS (Table 5) with a lower A V and a ∼7 × lower SFR that is driven by B2 with 850.1 effectively passive.This fit to the rest-frame UV to near-infrared underestimates the far-infrared luminosity of the system by around an order of magnitude and so is not reliable.These biases underline the need to use energybalance codes, such as MAGPHYS, when modeling the SEDs of galaxy populations that are suspected to have significant dust obscuration.
For the less-obscured 850.2, the best-fit FAST++ model yielded: 600 M e yr −1 , and . These parameters are closer to those from MAG-PHYS.The A V , stellar mass, and SFRs all agree within ∼1σ of the corresponding MAGPHYS measurements, given the significant uncertainties on the FAST++ estimates.This is as expected, as the obscured activity within 850.2 is less significant than in 850.1 and so the absence of energy balance in the SED fitting has less of an impact on the estimated properties.

Resolved SED Modeling
The high resolution and depth of the JWST and HST imaging of these sources can be used to go beyond the analysis just described to undertake resolved SED modeling of the sources to investigate their internal structures.Specifically, to map the distribution of stellar mass and dust reddening within the two systems (see also, Duncan et al. 2023;Kokorev et al. 2023).To perform this analysis, 9″ × 9″ thumbnails were extracted from the JWST and HST imaging centered on 850.1 and 850.2.The JWST thumbnails were weighted and coadded, the combined image smoothed with a Gaussian filter with an FWHM of 0 07, and then SEXTRACTOR, which Bertin & Arnouts (1996) used to construct segmentation maps.The segmentation maps were employed to identify regions of blank sky that were used to estimate the pixel-to-pixel noise and also to define areas associated with the various sources.To aid this for 850.1, the large-scale emission from the bright galaxy, A4, was first removed by rotating the images by 180°centered on A4 and subtracting these from the original frames.
To ensure that the final spatially resolved maps were not oversampled, the JWST and HST images were PSF-matched to the resolution of the F444W observations, ∼0 15.These images (and the associated segmentation masks) were then rebinned from the default 0 03 pixel −1 sampling to 0 18 pixel −1 , to provide nearly independent pixels.From this binned data, photometry was extracted in all of the JWST and HST bands for each 0 18 pixel, including the uncertainties measured from the pixel-to-pixel variations in the background regions.
To include information about the longer-wavelength emission, the SMA 880 μm map was used to provide a constraint on the long-wavelength pixel-to-pixel photometry.For each source, a model was constructed from the deconvolved 880 μm source properties measured by SMA and sampled at 0 18.For the 880 μm photometry, liberal uncertainties of 2 mJy per pixel were adopted to avoid over-constraining the MAGPHYS fits and also to reflect the fact that the resolution of the SMA map was coarser than the JWST and HST imaging.
For the MAGPHYS analysis of the pixel-based photometry, at least three detections were required at >3σ in the visible or near-infrared wave bands to include the fitted parameters in the analysis.Fits to pixels with too few photometric constraints  resulted in poorly constrained properties.However, the results did not change significantly if 2σ limits in three bands or 3σ in two bands were used instead, although the uncertainties became more significant.To derive the physical properties of each source, the high-redshift version of MAGPHYS was applied to the photometry for each pixel within the segmentation mask that defines each of the sources.A fixed redshift of z = 4.26 was used, but all other parameters were allowed to vary according to the standard priors used in the previous section, and also adopted in previous studies (e.g., Dudzevičiūtė et al. 2020).For the analysis of 850.1, the pixels that were associated with B1 were included, but then removed as this appears to be an unrelated foreground source.This conclusion was also supported by the χ 2 of the MAGPHYS fits for B1 at z = 4.26 being significantly worse than those in the remainder of the source (? 20, compared to ∼1).
Two further requirements were applied in the resolved analysis for a fit to a pixel to be included in the final maps.First, the MAGPHYS fit had to provide a reasonable description of the photometric observations, yielding a χ 2 < 10.This removed seven pixels in 850.2 and none in 850.1.Second, the predicted 850 μm flux density from the best-fit MAGPHYS SED for the pixel must not exceed 10 mJy.This caught situations where the predicted rest-frame far-infrared SED for the pixel, which was only weakly constrained by the SMA observations, was strongly at odds with the expected flux distribution for the source.This affected only four pixels from 116 in 850.1 and 12 from 415 in 850.2.Typically, these were in the lowest surface brightness regions on the periphery of the sources, and these pixels were replaced by the local median fit.With these cuts applied, the total submillimeter flux densities for the sources from summing the predicted MAGPHYS 850 μm fluxes were 52 mJy for 850.1, compared to 42.3 ± 1.6 mJy and 54.1 ± 1.5 mJy observed by SMA and SCUBA-2, respectively, and 14 mJy for 850.2, compared to 22.5 ± 1.6 mJy and 21.2 ± 1.5 mJy for SMA and SCUBA-2, respectively.Figure 5 shows the predicted full SEDs for both sources, constructed by summing the individual best-fit pixel SEDs.These provide very good fits to the visible to near-infrared, including the redder Spitzer bands that were not included in the resolved analysis.At longer wavelengths in the rest frame (far-infrared to submillimeter), the agreement was poorer, but  1 & 4) and the best-fit MAGPHYS model at the redshifts determined from the NOEMA CO observations, z = 4.26.The lower two panels, (c) and (d), show an expanded view of the restframe UV-near-infrared region of the SEDs.Two SEDs are shown for each source, one derived from summing the resolved-pixel-level SED fits (Σ pix , including B2) and one from fitting to the integrated source photometry.For the latter in 850.1, the individual fits to B2 and 850.1 were summed, although the emission longward of ∼3 μm is dominated by 850.1.Wave bands used in the resolved-pixel fits are shown as solid points.For the source-integrated fits, these bands were used as well as the unresolved photometry in the wave bands marked with open symbols.Also shown for comparison are the composite SED from the MAGPHYS analysis of z 3 SMGs from the AS2UDS survey by Dudzevičiūtė et al. (2020) and the SED of Arp 220, both normalized to the SMA 880 μm flux density.These illustrate that the combined SED of 850.1+B2 is similar to the AS2UDS composite between rest-frame B band (observed ∼2μm) and the near-infrared, but due to the contribution from B2 it becomes much bluer into the rest-frame UV, better matching Arp 220.850.1 alone is redder than the AS2UDS composite shortward of 2 μm.While 850.2 has an SED comparable to Arp 220 longward of ∼3 μm (rest frame ∼6000 Å) but bluer shortward of this and considerably bluer than the AS2UDS composite.In panel (c) the results of fitting to the photometry and limits for just 850.1 are shown (gray open squares), as well as a single MAGPHYS fit to the combined ("unresolved") photometry from 850.1 and B2, illustrating that this provides a poorer description of the observations, especially longward of 3 μm.In all panels, the resolved-pixel modeling, Σ pix , had only a single constraint beyond 4.4 μm from the SMA at 880 μm, yet it does at least as well out to 8 μm as the fits to the integrated source photometry (that included the IRAC 5.8 and 8.0 μm fluxes in the fit), although the absence of spatially resolved constraints at 450 μm resulted in poorer agreement in the rest-frame far-infrared region.
given the limited information available at these wavelengths for a spatially resolved analysis at JWST resolution, we view the fits as reasonably successful.
This resolved analysis allowed the internal variation in some of the key physical properties of the sources to be mapped on ∼0 18 scales, equivalent to 0.5-0.6 kpc in the source plane after accounting for lensing magnification.The discussion below focuses on the distribution of stellar mass and dust reddening, A V , as these two quantities are most closely tied to the photometric properties in the JWST and HST wave bands, where the observations have the highest resolution and signalto-noise.The absolute normalization of these maps will depend on the SED modeling code used (e.g., Carnall et al. 2019;Pacifici et al. 2023), but for comparison with the earlier results and published samples, we use MAGPHYS.The resolved 2D maps of the estimated stellar-mass density and dust reddening for 850.1 and 850.2 are shown in Figure 6.

Results and Discussion
Figure 2 illustrates the JWST/HST view of the surroundings of the two sources as well as their morphologies after correcting for the image shear due to lensing, while the variation in their brightness across the nine JWST and HST bands is shown in Figure 4. Strikingly, 850.1 is undetected in all of the HST ACS bands and the bluer JWST bands (the bulk of the galaxy only being significantly detected at 2 μm or above), while 850.2 displays emission down to F606W but then a strong drop between F606W and F435W (see Figures 4  and 5).
There are a number of sources detected with JWST and HST lying close to 850.1 and 850.2, with photometric redshifts reported in Table 5.These redshifts indicated that two of the galaxies near 850.1, including the bright disk galaxy A4, are probably members of the cluster A 1489 at z ∼ 0.35 and that similarly, the disk galaxy A5, near 850.2, is also likely to be a cluster member.This is unsurprising given that the two submillimeter sources are seen through the dense central regions of the foreground cluster.This analysis also suggested that none of the other galaxies visible in close proximity to 850.1 and 850.2 are likely to be physically associated with the submillimeter sources.

Gravitational Lens Model
A first strong-lens model for A 1489 was presented by Zitrin et al. (2020) based on multiple images identified in their HST observations.The analysis here uses a new lens model based on the same set of multiple images but employing a revised and updated version of the Zitrin et al. (2015) parametric code (Pascale et al. 2022;Furtak et al. 2023).In this model, cluster galaxies were modeled as double pseudo-isothermal ellipticals, and dark matter halos were modeled using pseudo-isothermal elliptical mass distributions.Two such halos were used for A 1489, centered on the two main mass clumps identified by Zitrin et al. (2020), and mass components were included to represent the cluster galaxies near 850.1 and 850.2, including A4 and A5.
The updated lens model predicts similar lensing magnification for 850.1 and 850.2 (Table 5) of μ = 4.0 - + 2.2 1.0 and μ = 5.6 - + 3.3 1.0 , respectively, where the uncertainties were judged from the range of acceptable models and are expected to be conservative.These values were employed when comparing the properties of these sources to other populations.Table 6 lists intrinsic, lensing-corrected properties for both sources.The NIRCam morphologies of the sources after correcting for the lensing magnification are shown in Figures 2(d) and (f).The stellar-mass surface density is per 0 18× 0 18 pixel, corresponding to ∼0.6 × 0.6 kpc in the source plane for 850.1 and ∼0.5 × 0.5 kpc for 850.2.B1 was excluded from the map for 850.1 due to the poor quality of the SED fits at z = 4.26, supporting its likely lower redshift.B2 has much lower A V and stellar-mass surface density compared to the bulk of the massive and dusty remainder of 850.1.These very distinct properties motivate the decision to fit B2 and the bulk of 850.1 separately and then sum the results.In contrast, 850.2 exhibits a smaller variation in A V and stellar-mass surface density across the source.
Accounting for lensing magnification, 850.1 has an intrinsic 850 μm flux density of S 850μm = 13.5 ± 0.4 mJy, while 850.2 has S 850μm = 3.8 ± 0.3 mJy.This shows that 850.1 is an example of a rare population of very luminous submillimeter sources with a surface density of ∼10 ± 5 degree −2 (Geach et al. 2017;Simpson et al. 2019;Garratt et al. 2023), while 850.2 comes from a much more numerous population around the knee in the 850 μm counts (Stach et al. 2018) with a surface density ∼50 × higher (∼ 500 degree −2 ).
Corrected for lens magnification, Δ ∼ 1.5 mag, 850.1 would have H ∼ 29 and K ∼ 27 for an extended source (not a point source as frequently adopted to define survey limits).This demonstrates that neither previous HST nor ground-based Kband surveys would have been able to detect this class of source, let alone if there are less-massive or higher-redshift examples.Thus, 850.1 represents an example of the K-faint subset of the dusty star-forming galaxy population (e.g., Simpson et al. 2014;Smail et al. 2021;Manning et al. 2022).These galaxies are massive but are missed by optical and nearinfrared galaxy surveys that seek to construct "stellar-mass" selected samples out to z ∼ 3-4.
After correcting their sky positions for lensing, the relative separation of 850.1 and 850.2 was estimated to be ∼100 kpc (∼ 15″) in the source plane.Thus, without the lensing magnification of the foreground cluster, these two galaxies would appear as a single blended bright (S 850μm ∼ 17 mJy) source in a single-dish submillimeter survey, similar to the associated multiple z ∼ 4.62 components found in two SCUBA-2 sources by Mitsuhashi et al. (2021).This small spatial separation, along with the line-of-sight velocity difference between 850.1 and 850.2 of only Δv ∼ 30 ± 20 km s −1 , indicates that these two galaxies are part of a small group at z ∼ 4.26 lying behind the core of the foreground cluster.Indeed, given the stellar mass for 850.1 estimated below, unless 850.2 has a very high transverse velocity; then it is very likely on a bound orbit and will merge at some point in the future.
This group behind A 1489 adds to a growing list of similar systems comprising one or more SMGs and other companion galaxies in groups at z ∼ 2-4, lensed by massive foreground clusters (e.g., Ledlow et al. 2002;Borys et al. 2004;Kneib et al. 2005;Caputi et al. 2021;Frye et al. 2023).In view of the large volumes accessible in surveys for submillimeter sources due to the negative K-correction, the frequency of identification of such systems may reflect the similarity between the spatial size of group-sized halos with mass comparable to those of SMGs (∼ 10 13 M e ; e.g., Stach et al. 2021) and the source-plane extent of the high-magnification region of the most massive clusters at z  1 (see Frye et al. 2023).

SEDs and Stellar Properties
The SEDs of the combined 850.1+B2 and 850.2 shown in Figure 5 are strikingly different.Shortward of 2 μm (rest-frame U band), the only emission from 850.1+B2 comes from the component/companion B2, with 850.1 itself undetected (see Figure 5(c)).In contrast, both the disk and the high surface brightness feature in 850.2 display a blue continuum that is detected down to the rest-frame far-UV (Figure 4).
The MAGPHYS analysis provided a lensing-corrected stellar mass from the sum of the fits to 850.1 and B2 M * = 5.5 - 400 M e yr −1 (Tables 5 and 6).These estimates are in reasonable agreement with the sum of the fits to the individual resolved pixels M * = 4.0 - 300 M e yr −1 , even though the latter fit lacked several of the mid-infrared and (sub)millimeter photometric constraints used in the former.If B2 is taken to be a separate galaxy from 850.1, the estimates for the system do not change appreciably, as the latter contains the majority of the mass and SFR.For 850.2, the estimates were: M * = 0.18 ± 0.07 × 10 11 M e and an SFR = 210 ± 20 M e yr −1 , with the resolved-pixel fits again in reasonable agreement: M * = 0.40 - 160 M e yr −1 .AGN contamination is not biasing the high stellar mass derived for 850.1 (or indeed 850.2), as the F444W morphology indicates only a small fraction of the light in the rest frame near-infrared could be from a central point source (see Figures 2(d) and 2(f)).Similarly, removing the IRAC 5.8 and 8.0 μm photometry from the MAGPHYS analysis did not change the best-fitting SEDs.
After correcting for lensing magnification, 850.1 is fainter than 850.2 at all wavelengths shortward of 8 μm (rest-frame H band), yet the estimated stellar mass for 850.1 is 30 × higher, while its SFR and dust mass are only ∼3-4 × higher than those derived for 850.2.The primary driver of this is the differences in ages for the best-fit SEDs.Note that 850.1 has a best-fit mass-weighted SED with an age of ∼450 Myr, corresponding to a mass-to-light ratio for the stellar population in the restframe H band (observed ∼8.5 μm) M/L H ∼4.5.In contrast, 850.2 has a much younger best-fit SED with an age of ∼50 Myr and M/L H ∼0.3.This difference in age accounts for a factor of ∼15 in the predicted masses, with the remainder due to the difference in extinction between the two galaxies, A V ∼ 5 mag and ∼1 mag, equivalent to ΔA H ∼1. The time needed to form the stellar masses of 850.1 and 850.2, assuming the current SFRs (averaged over a 100 Myr timescale), was ∼600 Myr and ∼90 Myr, respectively.These are comparable to the MAGPHYS mass-weighted ages for both galaxies, indicating that it is possible that the bulk of the stellar mass in both systems was formed in the ongoing star formation events, with 850.1 potentially forming at z form ∼ 6.

Morphologies and Structures
The extended emission of the two galaxies is well detected in the redder JWST bands and shows that 850.1 is morphologically more complex than 850.2 (Figure 2).The F444W 1.9 1.0 a Dual error bars give first the random and then the systematic uncertainties due to the lensing magnification.
b Sum of the parameter estimates for the fit to 850.1 and the fit to B2.
morphology of 850.1 in Figure 4 suggests a nearly face-on system with two arm-like features extending north and south from the ends of a linear structure of clumps.This bar-like feature has a source plane radius of ∼2-2.5 kpc, and spans a slightly brighter clump at the center of the galaxy.The 880 μm emission detected by SMA in 850.1 peaks on the central clump (Figure 2), and the resolved SED analysis shown in Figure 6 indicates that the stellar-mass surface density and the dust reddening also both peak around this region.The map shows a maximum stellar-mass density of 14 × 10 10 M e kpc −2 and a mean density within R e mass of (6.4 ± 1.4)× 10 10 M e kpc −2 .As surface densities, these quantities are magnification invariant.The resolved SED maps also indicate that B2, seen in projection against the northern part of the disk of 850.1, has a lower stellar mass surface density, ∼10 9 M e kpc −2 , than the bulk of 850.1 and also a much lower reddening, A V ∼ 0.3-1.
The variation as a function of scale of the orientation of the main structural features in 850.1 is illustrated in Figure 4. Fitted ellipses are overlaid at three surface brightness levels onto the F356W frame, showing a change in orientation and ellipticity in the inner regions, compared to the outskirts.This variation is similar to the characteristic change in orientation and ellipticity seen from a bar.More qualitatively, the morphology of this source (seen in the rest-frame I band) is reminiscent of the bar and arms in NGC 1365, the dominant galaxy of the Fornax group, as well as showing features similar to the arms and bars seen in dust continuum in high-redshift SMGs (Hodge et al. 2019).
Although the internal colors of 850.1 appear uniformly red in Figure 2(a), this is misleading, as the source is only well detected in the three filters longward of 2.8 μm, two of which were combined into the red channel.Figure 2 (panels (c) and (d)) uses the three reddest JWST bands and shows variation in color within the system, with the reddest region coinciding with the central clump.The map of A V from the resolved-pixel-level SED fits in Figure 6 also indicates a reddening gradient in the 850.1 that ranges from A V = 6.8 in the center to A V ∼ 3 in the outskirts.However, these reddening estimates are for light detected in F444W, and will be lower than what would be measured if longer-wavelength resolved photometry was included.This may explain why the median F444W-weighted A V from the pixel-level SEDs was A V = 3.6 - + 1.2 1.5 compared to A V ∼ 5 measured from the integrated photometry including longer wavelength constraints.
While the map of A V in 850.1 indicates that the central regions are highly obscured (probably even more obscured than shown in the map, as this corresponds to light that was still detectable at rest-frame I band; see Simpson et al. 2017), the outer parts of the galaxy are also surprisingly highly and uniformly reddened out to radii of ∼3-4 kpc.The simplest explanation for the widespread high A V in 850.1 may be that this galaxy has a very dusty disk (even if the gas fraction for the system is low; see Section 3.5), as well as an even more obscured nuclear region.Kokorev et al. (2023) presented a similar spatially resolved JWST analysis of a z = 2.58 Atacama Large Millimeter/submillimeter Array (ALMA)-detected galaxy that shows a comparably high dust extinction over a very wide area, although their system appears to be much more highly inclined than 850.1.Such extended dust disks have also been reported in similarly dust-rich galaxies (e.g., Hodge et al. 2016;Gullberg et al. 2019), in conjunction with more compact dust continuum structures that have been interpreted as bars and rings (Hodge et al. 2019).
Another, perhaps more speculative, explanation for this widespread (and relatively uniform) dust reddening in this (modestly inclined, in contrast to Nelson et al. 2023; see also Kokorev et al. 2023) galaxy is a dusty wind, as is seen on similar scales in M 82 (e.g., Yamagishi et al. 2012;Beirão et al. 2015), Mg II absorbers at higher redshifts (Ménard & Fukugita 2012), and also as theoretically expected (e.g., Krumholz & Thompson 2013).To illustrate the plausibility of this explanation, we assume a Milky Way gas-to-dust ratio, for which an A V ∼ 3 corresponds to a hydrogen column density of ∼10 22 cm −2 (Bouchet et al. 1985).Integrating over the extent of 850.1 gives a gas mass of ∼10 10 M e .If this gas (and associated dust) is being expelled in a wind with an outflow rate comparable to the current SFR, then this mass represents just ∼10 Myr of activity, during which time the material would cover the extent of 850.1 if it was traveling at the escape velocity.Thus, it is possible that the relatively uniform red colors of 850.1 may be due to a foreground screen of dust expelled by a starburst-driven wind.
Turning now to 850.2, this has an asymmetric structure with a bright elongated region (sheared by lensing) and a western extension.The overall appearance is of a dusty, disk-like galaxy (Figure 2).In contrast to the very red colors of the bulk of 850.1, both the high-surface brightness region of 850.2 and the western extension are detected in F606W, corresponding to rest frame ∼1100 Å, with neither region detected in F435W, rest frame ∼830 Å, indicating the presence of a strong Lyman break in the SED (Figure 4).However, in addition to these UVbright components, there is also faint red emission to the north that shows up as a band of higher-A V in Figure 6.This northern region includes a compact red component, seen in F277W, F356W, and F444W (Figures 2 and 4), that appears as a distinct entity in the source plane (Figure 2(f)).The peak of the SMA 880 μm emission lies between this red clump and the UV-bright region, suggesting that this may be an interacting system with both obscured and unobscured components, and with the unobscured part corresponding to the higher-frequency peak seen in the CO spectrum of this source.
The pixel-level SED fits for 850.2 give A V ∼ 0.3-1.0 over most of the UV-detected region, but rising to A V ∼ 3 at the red clump on the northern side.This behavior contrasts with that seen in local galaxies, where extinction typically rises toward the galaxy center (e.g., Holwerda et al. 2005;Keel et al. 2023).From the stellar mass map, the median stellar-mass density within the effective radius in 850.2 was (0.10 ± 0.01) × 10 10 M e kpc −2 with a peak of 0.7 × 10 10 M e kpc −2 .Thus, 850.2 has a more uniform distribution of stellar mass density and A V than 850.1.
Although 850.2 is very blue compared to 850.1, the UV light traces only ∼10% of the star formation in this system.In addition, there is a slight excess of emission in the Spitzer 5.8 and 8 μm filters above the best-fit model SED in Figure 5 that may indicate a steeply rising SED for the dustier red northern region (but this is far from dominant in F444W).Taken together with the strong far-infrared/submillimeter emission, indicating the presence of highly obscured star formation within the system and the velocity offsets seen between the optical and CO spectroscopy, this suggests that there could be two spatially distinct components within this source.However, in contrast to B2 and 850.1 in 850.2, the low-obscuration component dominates from the optical out to at least 4.4 μm, with the fainter and more obscured component only appearing at > 2 μm.
The GALFIT analysis quantified the stellar continuum structure of the two sources.The analysis of 850.1 in the F444W passband yielded a Sérsic parameter n = 0.9 ± 0.1 (with lower, but increasingly uncertain, values in bluer bands).A fixed n = 1 was therefore adopted to simplify comparison of the derived sizes between the bands.This showed a constant effective radius with wavelength of R e F444W ∼3.8 ± 0.4 kpc after correcting for lensing (see also Chen et al. 2022;Kamieneski et al. 2023).The similar analysis of 850.2 yielded n ∼ 1.9 ± 0.1 in F444W, with higher values in the bluer bands, although a combination of a de Vauncoleurs (Sérsic n = 4) and an n = 1 Sérsic did a better job of reproducing the asymmetric, extended emission in the source.However, for comparison to 850.1, n = 1 was adopted, which gave sizes that modestly decrease in bluer bands, but with low significance, and R e F444W ∼2.5 ± 0.3 kpc in F444W (also corrected for lensing).These sizes are comparable to those derived in the F444W wave band for SMGs with flux densities of S 850μm 1 mJy in the literature: R e F444W = 3.9 ± 0.8 kpc (Chen et al. 2022;Cheng et al. 2022Cheng et al. , 2023)).There was no evidence for a trend in rest-frame near-infrared size with either dust mass or redshift for the combined samples.
The dust continuum radii, R e dust , derived from the SMA observations (corrected for lensing magnification) were 2.8 ± 0.7 kpc for 850.1 and 1.9 ± 1.0 kpc for 850.2.Compared to the F444W sizes, these gave identical ratios of dust to F444W effective radii of 0.7 ± 0.2 and 0.7 ± 0.4, respectively, similar to previous estimates (e.g., Gullberg et al. 2019;Lang et al. 2019;Chen et al. 2022).These indicate the presence of compact dust-enshrouded regions in the dense cores of both galaxies, linked to the high A V regions (Simpson et al. 2017).
Finally, GALFIT was applied to the stellar mass distributions from the resolved, pixel-level SED modeling to attempt to derive stellar-mass radii with a Sérsic profile with n = 1.However, this proved very unstable, and so instead, simple numerical integration was used to derive the half-mass radii for the two sources: R e mass = 1.0 ± 0.2 kpc for 850.1 and 2.2 ± 0.4 kpc for 850.2, both corrected for lensing magnification.While the stellar mass-radius determined for 850.2 was in agreement with that measured at F444W, for 850.1 the stellarmass size was around ∼4× smaller than seen at F444W.This demonstrates that a significant fraction of the stellar mass in this massive system lies in a compact and very highly obscured nuclear region with the ratios of dust-to-stellar-mass sizes: 2.8 ± 0.9 and 0.9 ± 0.5 for 850.1 and 850.2, respectively.

Comparisons to Dust-mass-selected Populations
Figure 7 compares the derived properties of 850.1 and 850.2 with a representative sample of dust-mass-selected galaxies derived from ALMA interferometric identifications of ∼1000 sources found in ∼3 degrees 2 of wide-field 850 μm SCUBA-2 surveys in the UDS and COSMOS fields: AS2UDS (Stach et al. 2019) and AS2COSMOS (Simpson et al. 2020).These samples have been analyzed in an identical manner using MAGPHYS, allowing us to perform a simple comparison of their derived properties (Dudzevičiūtė et al. 2020;Simpson et al. 2020;S. Ikarashi et al. 2023, in preparation).After accounting for lensing magnification, Figure 7(a) shows that 850.1 and 850.2 lie at either end of the mass range seen in 850 μmselected galaxies at z > 4: 850.1 being one of the more massive galaxies found (from dust-selected samples or other methods) at these redshifts (see also, Casey et al. 2019), while 850.2 is one of the least massive in dust-selected samples (see also, Pope et al. 2023).Looking at their specific SFR, 850.1 lies below and 850.2 lies just inside the claimed "starburst" regime at this redshift.The strong linear trend in Figure 7(a) likely reflects the fact that these samples were selected primarily on their cold dust mass, the formation of which requires an equivalent characteristic stellar mass of ∼10 11 M e , with an observed dispersion of only ∼0.4 dex.
Figure 7(b) shows the variation of A V with stellar mass for 850.1 and 850.2 compared with the ∼40 ALMA-identified z > 4 dusty star-forming galaxies selected from SCUBA-2 mapping in the UDS and COSMOS fields as well as less-active K-selected star-forming galaxies at z > 4.This figure demonstrates one of the more striking differences between 850.1 and 850.2: they lie at diametrically opposite ends of the trend of A V with stellar mass for dusty galaxies at z ∼ 4. Thus, these two closely associated sources span the full range in reddening seen in the high-redshift submillimeter population, underlining the diversity in the restframe UV-to-visible properties for dusty star-forming galaxies, which was already present in the population, even in the same environment, at z  4. The two sources broadly follow the trend seen in the z > 4 K-selected galaxy population, which has the form A log V 10 ( ) = 0.42 + M log 10 10 11 ( ) * 0.23 with a median absolute deviation of ± 0.32, although this has a steeper slope than seen at z ∼ 0 (Garn & Best 2010; but see also Shapley et al. 2022;Gómez-Guijarro et al. 2023).
To test if there were other observed properties that would provide a better predictor of A V for this sample of z > 4 SMGs, the correlation of various properties (M * , M d , L FIR , SFR, T d , age, and Σ d ) with A V was analyzed using maximal informationbased nonparametric exploration (MINE; Reshef et al. 2011).This showed that the most significant correlation with A V for the SMG sample was indeed with stellar mass, with a 0.3% chance of random correlation using a jackknife test on the sample and calibrating the MINE likelihoods with a Monte Carlo simulation.There was only a minor reduction in scatter when additional corrections were included for either L FIR or M d (Gómez-Guijarro et al. 2023), and so it appeared that stellar mass was the most likely driver of the differences in dust reddening between 850.1 and 850.2 (although L FIR or SFR did nearly as well).
The final panel in Figure 7 showing purely SED properties is Figure 7(c), which compares the dust temperature and farinfrared luminosities from the MAGPHYS analysis of 850.1 and 850.2 to the wider dusty z > 4 population in UDS and COSMOS, where T d and M d were derived in very similar ways with MAGPHYS.The two galaxies lie within the scatter of the T d -L FIR trend at z  4, with 850.1 having a higher farinfrared luminosity and a marginally higher dust temperature (although with considerable uncertainties, even with the inclusion of the 450 μm SCUBA-2 photometry).The similarity in the dust temperatures is due to the similar ratios of farinfrared luminosity to dust mass-or equivalently interstellar medium (ISM) mass, as 850.1 and 850.2 have very similar dust-to-gas ratios (see below)-for the two sources:  and (7 ± 1)× 10 3 L e -M 1  for 850.1 and 850.2, respectively.This suggests that their far-infrared luminosities predominantly arise in similarly optically thick starburst regions (Scoville 2013).

Gas and Dynamical Masses
The lines detected in the NOEMA 3 mm spectra of the sources, along with their dust continuum properties from SMA and SCUBA-2, provide a wealth of information about the properties of the galaxies including both their gas dynamics and the characteristics of their ISM.Looking first at their ISM properties, the CO line luminosity ratios for the two galaxies are very similar: r 54 = 0.54 ± 0.01 and [C I]/CO(4-3) = 0.29 ± 0.01 for 850.1; r 54 = 0.61 ± 0.02 and [C I]/CO (4-3) = 0.29 ± 0.04 for 850.2.These compare to r 54 = 0.69 ± 0.06 and [C I]/CO(4-3) = 0.82 ± 0.06 for the well-studied z = 2.3 SMG the Cosmic Eyelash (Swinbank et al. 2010;Danielson et al. 2011), indicating that 850.1 and 850.2 have marginally less excited CO spectral line energy distributions (SLEDs) compared to the Eyelash (i.e., suggesting a flattening/turnover of the SLED around J up = 4-5).
The NOEMA spectra also provide kinematic information about the motion of cool gas in the galaxies.The Gaussian FWHM of the CO lines is the best measure of the kinematics of the gas in 850.1.For 850.2, both CO transitions show evidence for double-peaked lines, and so the mean velocity difference between the peaks, 330 ± 70 km s −1 , was used as the measure of the projected circular velocity.Inclination was estimated from the apparent axial ratio of the source measured in the F444W band and corrected for lensing shear (Table 2).To calculate the dynamical masses, a fixed radius of 10 kpc was chosen, as this should be sufficient to contain the majority of Figure 7.The stellar and interstellar medium (ISM) properties of 850.1 and 850.2.For comparison, several panels show samples of 870 μm selected ALMA-identified galaxies from the AS2UDS (Stach et al. 2019;Dudzevičiūtė et al. 2020) and AS2COSMOS (Simpson et al. 2019;S. Ikarashi et al. 2023, in preparation) surveys and the CO properties of a sample of submillimeter-continuum selected galaxies from Birkin et al. (2021).(a) The trend of specific SFR with stellar mass (SFR/M * -M * ) where the dotted line shows the starburst boundary at z = 4 from Speagle et al. (2014).(b) A V vs. M * , where the dashed and dotted lines show the best-fit relation and 1σ scatter for a K-selected z > 4 galaxies sample with M * > 10 10 M e drawn from the UKIDSS UDS survey (Lawrence et al. 2007) analyzed by Dudzevičiūtė et al. (2020).(c) The dust temperature vs. far-infrared luminosity (L IR -T d ) relation.(d) CO-derived gas fraction, M g /(M g + M * ), as a function of redshift, including lowredshift SMG-analogs from Oteo et al. (2017).(e) The ratio of the [C I]/CO(4-3) line luminosities vs. the [C I] to far-infrared luminosity ratio, which also include sources and photodissociation region (PDR) model tracks from Alaghband-Zadeh et al. (2013) and Valentino et al. (2020).The literature sources are flagged if they are estimated to have 20% contributions from an AGN to their far-infrared luminosities (Valentino et al. 2020).(f) Dust-to-stellar mass ratio, M d /M * , as a function of gas consumption timescale, M g /SFR, including two model tracks from Calura et al. (2017) for halos with baryonic masses of 10 11 M e and 10 12 M e roughly matching 850.2 and 850.1, respectively.The tracks were rescaled from their assumed Larson IMF with m c = 1.2 for comparison to the MAGPHYS derived values.
the stellar and cold-gas mass in these systems. 25Using r = 10 kpc and 10 11 M e , which both exceed the dynamical estimates, although they are still formally in agreement due to the large uncertainties on the latter from the inclination corrections and adopted sizes.
Using the estimated gas masses assuming α CO = 1, Figure 7(d) illustrates the variation in gas fraction of galaxy populations as a function of redshift for dust-mass-selected samples.The gas fraction of 0.16 ± 0.03 for 850.1 is at the low end of the distribution for dust-selected systems at z ∼ 4, while the estimate for 850.2 of 0.62 ± 0.07 is at the upper end and so more consistent with the rough trend for higher gas fractions in dusty galaxies at higher redshifts.

ISM Properties
Both 850.1 and 850.2 show lower [C I] luminosities relative to either their CO(4-3) or far-infrared luminosities compared to many of the high-redshift dust-mass selected galaxies or other [C I]-detected star-forming galaxies shown in Figure 7 ( ) ∼3.5-4.3 (in Habing units), depending upon the adopted tracks (see the discussion in Valentino et al. 2020).These ISM properties are within the ranges previously estimated for dusty star-forming galaxies and local ULIRGs (see Alaghband-Zadeh et al. 2013).
Assuming a simple uniform density sphere for the geometry of the ISM in these sources, then the ISM densities and gas masses indicate a radius of ∼0.4 kpc for 850.1 and ∼0.25 kpc for 850.2 (see Yan & Ma 2016).Adopting a more plausible thin-disk geometry would result in a larger radius by a factor of 1.9-2.5 for disks with thickness of ∼0.1-0.05their diameter.This would correspond to diameters of ∼1.8 kpc and ∼1.1 kpc, a factor of ∼1.6 smaller than the measured dust continuum sizes in both sources (Table 6).However, a more realistic model would involve a multiphase ISM, with much of the cool gas in a component that is potentially distinct from the starforming material, thus arguing for using a lower gas mass (relevant for the star-forming phase) in these calculations and hence a smaller volume.Even without this extra complication, it is clear that the dense ISM traced by the moderate-J CO lines in these galaxies is limited to a relatively small volume, either distributed as smaller clouds spread throughout the disks or more likely mostly restricted to a compact region in the galactic centers.
A similar calculation can be undertaken using the ISM radiation fields.For solar metallicities, a Salpeter IMF, and constant-SFR bursts with ages of 50-500 Myr, ∼10% of the bolometric luminosity arises in ionizing radiation with 6-13.6 eV.Assuming a centrally illuminated spherical shell geometry, the estimated radiation field from the PDR models implies shell radii r ∼ 1.1-2.9kpc for 850.1 and r ∼ 0.6-1.4kpc for 850.2.These sizes were smaller than those measured at F444W (Table 6) but are closer to those derived for the total stellar mass (and dust continuum) in Section 3.2 as well as by previous dust continuum studies (Simpson et al. 2015;Gullberg et al. 2019), although these assume different geometries.

Gas Masses and Fractions
Figure 7(f) shows the variation of dust-to-stellar mass ratio with the gas consumption timescale at the current SFR for the two galaxies.The trend in Figure 7(f) reflects the fact that galaxy-integrated dust mass provides a crude proxy for gas mass; hence M d ∝ M g , with similar dust-to-gas ratios of 69 ± 16 and 61 ± 8 for 850.1 and 850.2, slightly lower than the median of 120 ± 50 for the z > 4 subset from Birkin et al. (2021), with 60 ± 37 for their whole sample.The gas consumption timescales of both galaxies are ∼100 Myr, so they will consume their current gas reservoirs by z ∼ 4, unless further material is accreted.While the dust-to-stellar mass ratio for 850.2 is an order of magnitude higher than for 850.1, at M d /M * = 0.026 ± 0.004 compared to 0.0027 ± 0.0009, reflecting the very different stellar masses of the two galaxies.To interpret this difference, Figure 7(f) shows two tracks for proto-spheroid galaxy models (top-hat collapse) from Calura et al. (2017).These show the variation of dust-to-stellar mass ratio with gas consumption timescale for halos with baryonic masses of 10 12 M e and 10 11 M e , which roughly bracket the estimates for 850.1 and 850.2.These evolutionary tracks start at high M g /SFR and initially show rapidly increasing M d /M * until a peak at ∼150-200 Myr before declining, the tracks then terminate at the point that an outflow is expected to remove the remainder of the cold gas in the system.The higher mass model peaks at a lower dust-to-stellar mass ratio owing to the shorter star formation timescale leading to increased dust destruction, before the galactic wind terminates all star formation.These theoretical models demonstrate that the observed dust-to-stellar ratios for 850.1 and 850.2 can be achieved by models at typical ages of ∼100-300 Myr, with their relative dust-to-stellar mass ratios being consistent with 850.1 being older than 850.2.
The ISM properties of 850.1 and 850.2 are surprisingly similar in terms of dust temperature, ionization field, and density, given the very different stellar masses and rest-frame UV-visible obscuration of the stellar populations in the two galaxies.This similarity has its basis in the comparable ratios of dust mass to far-infrared luminosity (or SFR) for the two systems, which indicates similar relative amounts of ISM mass available to absorb and then reradiate the luminosity arising from the ongoing star formation (Scoville 2013).In contrast, the dust-to-stellar mass ratios of the two galaxies are very different, reflecting their different evolutionary stages.Note that 850.1 is a much more massive system and is relatively gas poor, suggesting that it may be in the process of quenching its activity, while 850.2 is a young and gas-rich system that is rapidly growing its stellar and dust mass.

Evolutionary Links
850.1 is one of the most massive galaxies known at z > 4 in terms of stellar or total baryonic mass at ∼10 11.8 M e (e.g., Caputi et al. 2011;Stefanon et al. 2015;Marsan et al. 2022), with a space density of ∼3 × 10 −7 Mpc −3 estimated from the surface density of comparably bright 850 μm sources and assuming a redshift range of z = 2-5 for this population (Simpson et al. 2020;S. Ikarashi et al. 2023, in preparation).This estimate is consistent with expected limits for such massive systems from Λ cold dark matter (CDM; < 5 × 10 −6 Mpc −3 ; Behroozi & Silk 2018).Unlike most previous studies of high-redshift massive galaxies, the high-resolution JWST imaging used here tightly constrains the rest-frame optical and infrared SED, and also excludes significant contamination from an obscured AGN in our estimates of the stellar mass.Likewise, the submillimeter imaging from SMA provided both the identification of this massive (but optically undetected) high-redshift galaxy and also constrained the influence of dust on 850.1ʼsSED necessary to derive a robust stellar mass (see also Ikarashi et al. 2022).
This analysis also benefited from the high-spatial resolution and long-wavelength coverage of JWST to go beyond what has been previously possible in resolved studies of the stellar content of dusty galaxies at high redshifts (e.g., Lang et al. 2019).The resolved stellar-mass maps from the pixel-level SED fitting indicated stellar mass surface density estimates within the central regions of 850.1 and 850.2 that are comparable to the stellar density of early-type galaxies at the present day, long believed to be the potential descendants of these systems.Moreover, the small stellar-mass sizes for the two sources, 1-2 kpc, are comparable to those claimed for the most compact, massive "quiescent/evolved" galaxies at z ∼ 1-3, ∼1-5 kpc (e.g., Waddington et al. 2002;Trujillo et al. 2006;Longhetti et al. 2007;Buitrago et al. 2008;Damjanov et al. 2009;Mowla et al. 2019;Lustig et al. 2021).The stellar-mass sizes for 850.1 and 850.2 are also consistent with the expectations from hydrodynamic simulations (e.g., Wellons et al. 2015), which indicate sizes of ∼1 kpc for systems with masses of ∼10 11 M e .
The mean stellar-mass density within R e mass for 850.1 is estimated to be (6.4 ± 1.4)× 10 10 M e kpc −2 , which is in excellent agreement with the value of 7 × 10 10 M e kpc −2 predicted by extrapolating the trend for quiescent galaxies (S e Q ) in (Barro et al. 2017).The formation of such a dense core is thought to be a key requirement for the subsequent quenching of star formation in the progenitors of quiescent galaxies at high redshift (e.g., van Dokkum et al. 2014;Barro et al. 2017).Theoretical simulations suggest that the core formation process involves compaction driven by strongly dissipative processes (Dekel et al. 2009;Tacchella et al. 2016), including interactiondriven instabilities that likely form bars as we see in 850.1, followed by quenching, which may link to the relatively low gas fraction we infer for 850.1.
Figure 8 shows the intrinsic sizes and stellar-mass densities for 850.1 and 850.2 compared to their lensing-corrected stellar masses.These are plotted alongside estimates from integrated SED fitting and JWST NIRCam F444W sizes for z = 1-3 SMGs with S 850μm ∼2-4 mJy from Chen et al. (2022).The F444W-based sizes may overestimate the true stellar-mass sizes for these systems, particularly at higher redshifts due to the centrally concentrated dust, which is still optically thick at the observed wavelengths (Lang et al. 2019).Although, the examples shown from Chen et al. (2022) lie at z  3, where F444W is sampling the rest-frame near-infrared, and so the influence of dust extinction is expected to be less.In addition, Figure 8 includes a sample of high-mass spectroscopically classified quiescent galaxies at somewhat lower redshifts, z ∼ 1.2-1.7,with estimated stellar masses above 3 × 10 11 M e , and sizes derived from HST F160W imaging (Longhetti et al. 2007).These comparisons show that 850.2 has a stellar-mass size and stellar-mass surface density consistent with those seen in early-type galaxies at the present day, as do the other lowerredshift SMGs with comparable 850 μm flux densities from Chen et al. (2022).However, 850.1 is an outlier in both panels, with a size that is a factor of ∼10 smaller than comparably massive early-type galaxies at z ∼ 0. Nevertheless, the peak stellar-mass surface density for 850.1, 1.4 × 10 11 M e kpc −2 , is similar to the characteristic maximum stellar-mass density estimated in the central regions local early-type galaxies (e.g., Hopkins et al. 2009), suggesting a potential link.
Empirical estimates of the size evolution of massive, earlytype galaxies out to z ∼ 3-4 indicate that these increase at lower redshifts as (1 + z) −1.39±0.13(e.g., Trujillo et al. 2007;Cimatti et al. 2012;Holwerda et al. 2020;Ito et al. 2023).This implies a size increase of ∼4 × for 850.1 from z = 4.26 to z = 1.5 and ∼10 × to the present day.The evolution in the size of massive galaxies is driven by minor mergers (such as the potential accretion of B2 and 850.2 onto 850.1) that increase the total stellar mass of the system modestly but boost the effective radius considerably (Hopkins et al. 2009;Naab et al. 2009;Oser et al. 2012).The size evolution for lower-mass galaxies is less extreme (Genel et al. 2018) and may be moot for 850.2 if it subsequently merges with 850.1.As Figure 8 shows, if 850.1 undergoes size evolution as suggested by simulations and observations at z < 3, then at z ∼ 1.5 it will end up in the region of the R e -M * and Σ * -M * planes that are populated by the high-mass quiescent galaxies from Longhetti et al. (2007), before evolving to match the properties of the most massive early-type galaxies at z ∼ 0. These massive galaxies at z ∼ 0 are also those that appear to have formed the bulk of their stars at z > 4, again in good agreement with the observations of 850.1.Overall, this appears to be strong circumstantial evidence that 850.1 represents the formation phase of an ultramassive early-type galaxy, which will subsequently evolve through the massive, passive galaxy population found at z ∼ 1-3 (Trujillo et al. 2006;Buitrago et al. 2008;Damjanov et al. 2009).
To summarize, 850.1, 850.2, and B2 likely reside in a grouplike halo serendipitously discovered behind the core of A 1489.The two submillimeter sources represent extremes of the dust-mass-selected population at z > 4, with 850.1 being one of the most massive galaxies known at these redshifts, with a high stellar-mass-density core.The gas reservoirs in the two galaxies contribute very different fractions of their total baryonic content at this early phase of the Universe (age of 1.4 Gyr): gas fractions of 16% and 62% for 850.1 and 850.2, respectively, suggesting different evolutionary states.These differences were also reflected in the different characteristic ages of their best-fit SEDs and dust/stellar ratios, although both or more speculatively, a dusty wind emanating from the ∼1 kpc stellar core.
3. The high-quality multiband JWST and HST images available for this field also allow for a pixel-level SED analysis that provides unique information on the internal structures of the galaxies.This resolved SED analysis shows good correspondence with the integrated SED fits, but only because it showed the need to exclude the foreground/companion UV-bright sources, B1 and B2, from the integrated photometry fits for 850.1.If these sources had been included in a single SED fit, this would result in a significantly worse fit and biased derived properties.Similar issues may affect integrated photometry modeling of systems comprising a mix of lowand high-obscuration regions or components.This effect may also explain previous reports of mid-infrared power-law excesses in ∼10% of dusty galaxies that were previously attributed to buried AGN.Such resolved analyses benefit considerably from the inclusion of resolved long-wavelength constraints on the pixel-level SEDs, as provided by sensitive submillimeter interferometry.
4. 850.1 and 850.2 are only ∼100 kpc (∼ 15″) apart in the source plane and offset in velocity by only ∼ 30 km s −1 , indicating that they likely reside in a single structure, a small group at z ∼ 4.26.Without the lensing effect of the foreground cluster, they would appear as a single, blended, bright (S 850μm ∼17 mJy) source in a single-dish submillimeter survey.There are a number of similar lensed groups at z ∼ 2-4, comprising one or more dusty star-forming galaxies and other companions, reported in the literature.The frequency with which these systems are being found may reflect the correspondence between the spatial scale of group-sized halos and the size of high-magnification region of massive clusters at z  1; see also Frye et al. (2023).This adds evidence that at least a subset of the high-redshift submillimeter population resides in groupscale halos.These are expected to be the most active environments for the accretion of gas from the surrounding intergalactic medium, as well as companion galaxies, potentially providing both the fuel and the trigger needed to power their active star formation, 5.The properties of 850.1 suggest it is the central galaxy of this group.Its stellar size and central stellar-mass surface density will match those measured for massive quiescent galaxies at z ∼ 1.5 and z ∼ 0 if it follows the expected size evolution due to the accretion, through minor mergers, of satellite galaxies such as 850.2 (which itself may be an interacting system) and potentially B2.
6.The exquisite JWST imaging of the A 1489 field, combined with the natural magnification of the foreground cluster lens, provides a striking view of the stellar structure of these two strongly star-forming galaxies.The structure of 850.1 is particularly noteworthy, with features that appear to correspond to arms and a central bar in a very massive galaxy at z = 4.26.It is tempting to identify these features as a response to dynamical perturbations from the group environment, including the presence of 850.2 as a close companion.Bars with ∼1 kpc radii have been reported from high-resolution dust continuum imaging of high-redshift SMGs by Hodge et al. (2019), as well as from statistical analyses of the dust continuum shapes of the submillimeter galaxy population (Gullberg et al. 2019).These features have similar sizes to the dust continuum peak and the massive stellar core in 850.1, suggesting that this galaxy is experiencing a bar-driven inflow that is fueling a central starburst in its baryon-dominated core.
The study of high-redshift dusty star-forming galaxies with JWST/NIRCam within the PEARLS survey will be expanded beyond these two bright lensed examples to the more typical examples of the population at millijansky-flux levels with the ongoing JWST survey of the PEARLS Time Domain Field at the North Ecliptic Pole, which has been imaged with SCUBA-2 by Hyun et al. (2023).

Figure 1 .
Figure1.A three-color JWST NIRCam image of a 100″ × 100″ region (north at the top, east to the left) around the z = 0.35 cluster A 1489 (where F090W+F150W is shown as blue, F200W+F277W is green, and F356W+F444W is red).Signal-to-noise contours are shown, representing the SCUBA-2 850 μm emission (starting at 5σ, in 5σ increments, 14 5 FWHM beam), the SCUBA-2 450 μm emission (3σ and 4σ, 7 5 FWHM beam), and the NOEMA 3 mm continuum emission (10σ, ∼5″ FWHM beam).Note that 850.1 and 850.2 are the two bright sources visible at 450 μm, 850 μm, and 3 mm to the northeast of the cluster core, which is identified by several bright elliptical galaxies, as well as associated strong lensing features.The counterparts to these sources are unambiguously identified using the higherresolution SMA 880 μm observations shown in Figure2.

Figure 2 .
Figure 2. Three-color images of 9″ × 9″ regions around (a) 850.1 and (b) 850.2 constructed from JWST NIRCam imaging, with F090W+F150W as blue, F200W+F277W as green, and F356W+F444W as red.Panels (c) and (e) show ¢¢5 5″ regions using F277W, F356W, and F444W as blue, green, and red and the same filter combination is used in panels (d) and (f), which show views of 850.1 and 850.2 "de-sheared" to correct for the lensing magnification (the scale bar indicates 1 kpc).Panels (a) and (b) also show the SMA 880 μm contours (starting from 2σ in 7σ increments), these unambiguously identify the two submillimeter sources.Note that 850.1 corresponds to an extremely red, highly structured source lying close to a bright, foreground disk galaxy with several fainter sources in close proximity.In contrast, 850.2 is more isolated and appears to be a bluer, disk-like system.The various sources and features discussed in the text are identified on each panel: distinct sources have labels starting with "A," "B" for potential subcomponents, and "C" for likely emission-line features.The colors and photometric redshifts suggest that "knot" B1 near 850.1 is likely to lie at z < 4.26, while B2 is probably a close/interacting companion seen in projection to 850.1 or a less-obscured component within the galaxy.The photometric-redshift analysis also suggests that A1, A4, and A5 are likely to be cluster members at z ∼ 0.35, while A2 and A3 are probably behind the cluster, but foreground to 850.1 and 850.2.

Figure 3 .
Figure3.The NOEMA 3 mm spectra for 850.1 (left) and 850.2 (right), showing the CO(4-3), [C I](1-0), and CO(5-4) emission lines.The panels are centered on the expected frequencies of the lines at the adopted redshifts for each source from Table2and a model fit comprising a Gaussian and a constant continuum level is shown for each line.Both sources have relatively weak [C I] emission compared to their CO(4-3) and CO(5-4).A single Gaussian provided an adequate description of the emission for the lines in 850.1, but all three lines in 850.2 display a dip at the systemic velocity, suggesting a double-peaked profile, with the higher-frequency CO peak close to the redshift derived from the rest-frame UV spectroscopic features.The spectra for 850.2 have been rebinned to 40 MHz per channel, except for the [C I] line, which has been rebinned to 100 MHz per channel (and as a result is shown over a slightly wider frequency window to illustrate the off-line noise).

Figure 4 .
Figure 4. 6″ × 6″ thumbnails of 850.1and 850.2.From left to right, the panels show the HST ACS F435W, F606W, and F814W imaging and JWST NIRCam F090W, F150W, F200W, F277W, F356W, and F444W images.The striking red colors of 850.1 are clear as the bulk of the emission from the galaxy, apart from the knots B1 and B2 (we suggest that the former is likely to be foreground), is undetectable shortward of 2 μm (the F200W band).In contrast, the high surface brightness regions of 850.2 are detected down to F606W (corresponding to ∼1100 Å in the rest frame at z = 4.26), but are undetected in F435W (rest frame ∼830 Å) showing the presence of a strong Lyman break.Fitted ellipses at three isophote levels are overlaid on the F356W image of 850.1, illustrating the "bar"-like feature seen in the central regions of the galaxy.At z = 4.26, the 912-Å Lyman break falls at the extreme red end of the F435W filter transmission, while the Lyα emission line falls in F606W, Hα in F356W (strong emission from which likely explains the pointlike feature C1 in 850.2, which probably corresponds to a giant [HII] region with an intrinsic size of ∼50-100 pc), and the nebular emission lines [OIII] 4959/5007 are in F277W and [OII] 3727 lies in F200W.The various sources and features discussed in the text are identified on the F090W thumbnails.

Figure 5 .
Figure 5. (a) The spectral energy distributions (SEDs) for 850.1 (including B2) and (b) 850.2, showing the measured photometry (Tables1 & 4) and the best-fit MAGPHYS model at the redshifts determined from the NOEMA CO observations, z = 4.26.The lower two panels, (c) and (d), show an expanded view of the restframe UV-near-infrared region of the SEDs.Two SEDs are shown for each source, one derived from summing the resolved-pixel-level SED fits (Σ pix , including B2) and one from fitting to the integrated source photometry.For the latter in 850.1, the individual fits to B2 and 850.1 were summed, although the emission longward of ∼3 μm is dominated by 850.1.Wave bands used in the resolved-pixel fits are shown as solid points.For the source-integrated fits, these bands were used as well as the unresolved photometry in the wave bands marked with open symbols.Also shown for comparison are the composite SED from the MAGPHYS analysis of z 3 SMGs from the AS2UDS survey byDudzevičiūtė et al. (2020) and the SED of Arp 220, both normalized to the SMA 880 μm flux density.These illustrate that the combined SED of 850.1+B2 is similar to the AS2UDS composite between rest-frame B band (observed ∼2μm) and the near-infrared, but due to the contribution from B2 it becomes much bluer into the rest-frame UV, better matching Arp 220.850.1 alone is redder than the AS2UDS composite shortward of 2 μm.While 850.2 has an SED comparable to Arp 220 longward of ∼3 μm (rest frame ∼6000 Å) but bluer shortward of this and considerably bluer than the AS2UDS composite.In panel (c) the results of fitting to the photometry and limits for just 850.1 are shown (gray open squares), as well as a single MAGPHYS fit to the combined ("unresolved") photometry from 850.1 and B2, illustrating that this provides a poorer description of the observations, especially longward of 3 μm.In all panels, the resolved-pixel modeling, Σ pix , had only a single constraint beyond 4.4 μm from the SMA at 880 μm, yet it does at least as well out to 8 μm as the fits to the integrated source photometry (that included the IRAC 5.8 and 8.0 μm fluxes in the fit), although the absence of spatially resolved constraints at 450 μm resulted in poorer agreement in the rest-frame far-infrared region.

Figure 6 .
Figure6.6″ × 6″ thumbnails of 850.1 and 850.2 showing the JWST F200W and F444W images along with the pixel-based MAGPHYS estimates of the log-scaled stellar mass and A V .The stellar-mass surface density is per 0 18× 0 18 pixel, corresponding to ∼0.6 × 0.6 kpc in the source plane for 850.1 and ∼0.5 × 0.5 kpc for 850.2.B1 was excluded from the map for 850.1 due to the poor quality of the SED fits at z = 4.26, supporting its likely lower redshift.B2 has much lower A V and stellar-mass surface density compared to the bulk of the massive and dusty remainder of 850.1.These very distinct properties motivate the decision to fit B2 and the bulk of 850.1 separately and then sum the results.In contrast, 850.2 exhibits a smaller variation in A V and stellar-mass surface density across the source.
e and an SFR = 900 - + 300 e and SFR = 1200 - + 100 e , confirming that the baryonic masses are broadly consistent with the galaxy kinematics derived from the CO line width, given the uncertainties on the inclination corrections.If the calculation had instead used a Milky Way-like α CO = 4.3 to estimate the cold-gas mass, then the baryonic masses would rise (e).Following Alaghband-Zadeh et al. (2013), a simple photodissociation region (PDR) model can be used to assess the physical origin of these differences.The ratio [C I]/CO(4-3) traces the characteristic density of the ISM, and [C I]/L FIR indicates the strength of the radiation field, although the absolute values are sensitive to the chosen model (compare the tracks from Alaghband-Zadeh et al. 2013 and Valentino et al. 2020 in Figure 7(e)).The Kaufman et al. (1999) PDR model grids indicate well-constrained ISM densities for 850.1 and 850.2, n log 10 ( ) = 4.5 ± 0.1 cm −3 , and radiation fields of G log 10 0

Table 4
Observed Visible and Near-infrared Photometry
a b z ph =

Table 6
Intrinsic Properties of Sources