UNCOVER: The rest ultraviolet to near infrared multiwavelength structures and dust distributions of sub-millimeter-detected galaxies in Abell 2744

The universe looks very different in the sub-millimeter (sub-mm) than at shorter wavelengths. Galaxies which are bright in the sub-mm are often faint or undetected at optical or near infrared (NIR) wavelengths, as their high dust content attenuates the shorter wavelength light and then re-emits this energy as thermal continuum in sub-mm to millimeter (mm) wavelengths (e.g., Casey et al. 2014). These sub-mm detected galaxies are typically distant ($z\gtrsim 1$), massive, dusty, and star-forming (e.g., Casey et al. 2014, Smith et al. 2017, Dudzevičiūtė et al. 2020; Hodge & da Cunha 2020). However, their detailed rest-frame NIR structures, which contain key information about their formation, remained uncertain with previously available observations.

The sensitivity of sub-mm and mm interferometers, particularly ALMA, has opened a window to directly study the impacts of dust well below the limits of ultraluminous infrared galaxies (ULIRGs) and classical bright sub-millimeter galaxies (SMGs), revealing dust continuum and obscured star formation even in “main sequence” star-forming galaxies. Probing the faint sub-mm population requires deep field or lensing cluster surveys, and there have been relatively limited ALMA studies in the $\sim 0.1-1\,\mathrm{mJy}$ flux density regime to date (e.g., Fujimoto et al. 2016, Aravena et al. 2016, González-López et al. 2020, Gómez-Guijarro et al. 2022, Muñoz Arancibia et al. 2023). However, these surveys have shown that the dust content in the fainter population is still significant, with $\sim 85\%$ of star formation at cosmic noon ($z\sim 1-3$) obscured by dust (e.g., Dunlop et al. 2017). Thus, a full picture of galaxy assembly at cosmic noon requires an understanding of these lower sub-mm luminosity sources.

While the detailed substructures of bright dusty objects such as SMGs have been accessible with ALMA and NOEMA (revealing merging starbursts as well as disks; e.g., Hodge et al. 2016, 2019, Rujopakarn et al. 2016, Calistro Rivera et al. 2018, Lang et al. 2019, Puglisi et al. 2019, Tadaki et al. 2020), reaching the fainter population at high spatial resolution to map the dust continuum has remained a challenge, requiring very long exposures even with ALMA’s sensitivity. Little is known about this population’s intrinsic structures, including their distribution of dust and active star formation, leaving us without a complete picture of how this faint sub-mm population builds up or whether this population represents a particular evolutionary phase in galaxy formation.

An alternative, complementary avenue to characterize the dust distribution and obscured structures in faint sub-mm sources is using high-resolution rest-optical imaging to probe the stellar distribution directly. Prior to JWST, the high spatial resolution and sensitivity Hubble Space Telescope (HST) enabled detailed structural measurements, but its limited wavelength coverage (up to $\sim 1.6\mu\mathrm{m}$) probed only as far as the rest frame optical at cosmic noon. And while Spitzer had longer wavelength coverage to probe the rest-frame near infrared (NIR) at these redshifts, it lacked the spatial resolution necessary to map stellar structures. Consequently, past structural measurements have typically been based on one HST filter at the rest-optical. However, dust complicates the interpretation of galaxies’ rest-optical structures. In particular, dust (as well as age and metallicity, which also contribute to radial color gradients) has been shown to result in dramatically different stellar size estimates depending on the dust concentration and distribution, as probed by relative differences between half-mass and half-light sizes (e.g., Wuyts et al. 2012, Liu et al. 2016, 2017, Suess et al. 2019, Miller et al. 2023, Zhang et al. 2023). Mock observations have also found dust geometry to be a major complication in structural interpretation of dusty galaxies (Popping et al. 2022, Zhang et al. 2023).

With the launch of JWST (Gardner et al., 2023), it is now possible to obtain high resolution imaging over an expanded wavelength range from $\lesssim 1\mu\mathrm{m}$ out to $\sim 4.4\mu\mathrm{m}$ (corresponding to the rest-frame $\sim 0.3-1.4\mu\mathrm{m}$ at $z\sim 1-3$), providing information across the full stellar spectral energy distribution (SED). Critically, this coverage includes the rest-frame NIR probing the largely unattenuated stellar continuum. In addition to enabling robust measurements of sizes and stellar densities, comparing the rest-frame NIR with the rest-frame UV/optical (“color gradients”) allows us to constrain radial variations in mass-to-light ratios. These color gradients can be driven by non-uniform dust distributions or variations in relative stellar ages and/or metallicities; therefore, color gradients also indicate where active growth occurs in galaxies now relative to their past (e.g., Suess et al. 2021, 2022, 2023). Early results from JWST point to dust distribution as a major contributor to observed color profiles (e.g., Miller et al. 2022; also Suess et al. 2022) and without the simultaneous rest-NIR coverage, intrinsic size estimates may be biased towards larger values (e.g., Zhang et al. 2023, among others).

By combining JWST/NIRCam imaging with existing shorter wavelength HST imaging, we can now perform detailed multiwavelength structural studies of the distribution of stars, dust, and star formation in extremely dusty distant galaxies, revealing new insights about the nature and formation of this population. A number of studies have already leveraged early JWST observations to constrain the structures of distant dusty galaxies, including both sub-mm detected galaxies (including Cheng et al. 2022, 2023, Gillman et al. 2023, Kokorev et al. 2023, Rujopakarn et al. 2023, Smail et al. 2023, Wu et al. 2023, Boogaard et al. 2024) and other dusty galaxies (e.g., Gómez-Guijarro et al. 2023, Le Bail et al. 2024, Liu et al. 2023, Nelson et al. 2023, Magnelli et al. 2023). Building upon these works, we study the structures of a sample of galaxies within a cluster field, where the strong lensing boost allows us to examine the relation between total dust content and detailed multiwavelength structures for a variety of galaxies pushing to faint intrinsic fluxes and compact physical sizes.

Here we present a multiwavelength ($0.5-4.4\mu\mathrm{m}$) structural analysis of 11 $z\sim 0.9-3.5$ sub-mm detected ($S_{\mathrm{1.2mm,int}}\sim 0.2-1.5\,\mathrm{mJy}$) galaxies in the Abell 2744 cluster field as measured with HST and JWST. This sample includes all galaxies detected in $1.2\,\mathrm{mm}$ ALMA continuum observations of the Abell 2744 primary cluster core from the ALMA Lensing Cluster Survey (ALCS; Fujimoto et al. 2023a) and ALMA-Hubble Frontier field (ALMA-HFF; Muñoz Arancibia et al. 2023) programs. We visually examine the rest-optical and NIR structures, then quantify the multiwavelength structures by fitting single component models and deriving residual-corrected flux profiles in every filter for the full sample. We also calculate the ratio of sizes between long ($>4.4\mu\mathrm{m}$) and short ($\sim 1.5-2.7\mu\mathrm{m}$) wavelengths, which is sensitive to color gradient strengths probing radial changes in dust, age, or metallicity. Dust is expected to be an important driver of these color gradients, given this sample is sub-mm detected and a subset exhibits signatures consistent with dust lanes or patchy dust distributions. Finally, we examine the structural measurements together with the intrinsic, lensing-corrected 1.2mm flux densities, to probe the impact of both dust quantity and geometry on these galaxies’ multiwavelength structures.

Throughout this work we adopt a flat $\Lambda$CDM cosmology with $\mathrm{H_{0}=70\,km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{M}=0.3$, and $\Omega_{\Lambda}=0.7$.

To study the rest-frame UV to near infrared (NIR) structures of sub-mm detected galaxies at cosmic noon ($z\sim 1-3$), we leverage the combination of low spatial resolution ALMA continuum observations (beam FWHM $\sim 0.9\arcsec\times 0.7\arcsec$; Fujimoto et al. 2023a) with deep, high spatial resolution JWST/HST coverage available for the Abell 2744 cluster field. We select galaxies from ALCS (Kokorev et al. 2022, Sun et al. 2022, Fujimoto et al. 2023a; including ALMA-HFF observations; González-López et al. 2017, Muñoz Arancibia et al. 2023), as presented in Fujimoto et al. (2023a). We include both sources in the ALCS 1.2mm continuum blind ($\mathrm{SNR_{natural\,map}}\geq 5$; $N=6$) and IRAC-prior ($4.0\leq\mathrm{SNR_{natural\,map}}<5$ and $\mathrm{SNR_{IRAC,CH2}}\geq 5$; $N=5$) detection catalogs, for a total of 11 galaxies.

For this analysis, we use all public JWST/NIRCam (F115W, F150W, F200W, F277W, F356W, F410M, F444W) and HST/ACS (F435W, F606W, F814W) and WFC3 (F105W, F125W, F140W, F160W) imaging covering Abell 2744. Mosaics for all filters are produced with a common WCS grid using grizli (1.8.16.dev12; Brammer 2019, Kokorev et al. 2022) and astrodrizzle (Gonzaga et al., 2012), with JWST/NIRCam short wavelength (SW; F115W, F150W, F200W) filters drizzled onto a 20mas scale, while JWST/NIRCam long wavelength (LW; F277W, F356W, F410M, F444W) and HST filters are drizzled onto a 40 mas scale.¹¹1https://s3.amazonaws.com/grizli-v2/JwstMosaics/v7/index.html See Bezanson et al. (2024) for full details of the imaging reduction and mosaicing. Background subtraction and modeling and subtraction of bright cluster galaxies is then performed, as detailed in Weaver et al. (2024).

The ALCS sources are matched to JWST-detected sources from the UNCOVER photometric catalog (DR2; Weaver et al. 2024) with a search radius of $0.5\arcsec$. NIR counterparts for all ALCS sources are identified. However, in two cases (UNCOVER IDs 9018, 24143; denoted throughout with an asterisk), visual inspection reveals the galaxies are shredded into multiple sources in the UNCOVER catalog as a result of the detection deblending thresholds, which were optimized to detect faint, high-redshift sources. As the photometry and segmentation maps are used for simultaneous modeling or masking of neighboring objects in the structural modeling (see Sec. 3), for these two objects we instead perform an alternative detection (using SEP; Barbary 2018) with less aggressive deblending parameters (KERNEL=3.5 pixel FWHM Gaussian, MINAREA=3 pixels, $\texttt{THRESH=0.8}\,\sigma$, DEBLEND_NTHRESH=8, DEBLEND_CONT=0.0001, CLEAN=Y, CLEAN_PARAM=1.0), and use the resulting segmentation maps and photometric parameters (derived using the same photometric pipeline as Weaver et al. 2024 except for the parameter changes) for structural modeling of these two galaxies. We calculate the magnifications ($\mu$) for each object from the updated lensing model by Furtak et al. (2023, v1.1), adopting the best-available redshifts. We also estimate stellar masses from the HST and JWST photometry, as well as the 1.2mm ALMA fluxes (from the most recent DUALZ catalog; Fujimoto et al. 2023b) using the same procedure as the UNCOVER DR2 stellar population synthesis (SPS) catalog (Wang et al., 2024), with the following modifications. First, we fix the redshifts to $z_{\rm spec}$ or $z_{\rm gris}$ when available. Second, we use the alternative, less deblended photometry measured for the aforementioned two objects. For the single object with only a $z_{\rm phot}$ (16840), we use the redshift and stellar mass directly from the DR2 SPS catalog (Wang et al., 2024). For clarity, we emphasize that all physical properties (i.e., stellar mass and physical size) and the intrinsic 1.2mm flux densities are corrected for gravitational lensing using the magnification $\mu$.

We use deep Chandra X-ray maps of Abell 2744 (2.1Ms, presented in Bogdán et al. 2024, considering both blind and IRAC-prior detections) to determine if any of these objects host AGN. We find that the majority of our sample is X-ray detected (7/11: 24143*, 26131, 43148, 24852, 21972, 14034, 24823), likely hosting AGN. This suggests that the high AGN fraction amongst somewhat brighter sub-mm selected galaxies (e.g., Serjeant et al., 2010; Shim et al., 2022) continues even towards lower sub-mm fluxes. This X-ray detection fraction is even higher than found by previous studies for brighter sub-mm-detected galaxies (e.g., Laird et al., 2010; Johnson et al., 2013; Wang et al., 2013).²²2We defer any comparison of AGN and X-ray-detection fractions among faint ($\lesssim 1\,\mathrm{mJy}$) sub-mm-detected galaxies compared to brighter sub-mm-detected populations to future analysis. Objects which are X-ray detected are noted in figures throughout this work. As discussed later (see Sec. 3, Appendix A), we determine that the host galaxy dominates the light even out to $4.4\mu\mathrm{m}$ in 5 of these objects (24143*, 26131, 43148, 24852, 21972); we therefore infer that these objects host highly obscured AGN, such that we can still analyze the wavelength-dependent host galaxy morphology. We infer that the other 2 objects (14034, 24823) host dust-reddened AGN, which start to appreciably contribute to the total flux at $\gtrsim 2\mu\mathrm{m}$.

Color composite JWST images for our sample are shown in Fig. 1, and the sample IDs, positions, best-available redshifts, magnifications, and stellar masses are listed in Table 1. For structural modeling, we use empirical PSFs constructed for each band from unsaturated stars within the mosaic, renormalized so the energy enclosed within $4\arcsec$ aligns with typical calibration levels (see Weaver et al. 2024 for full details).

To measure the multiwavelength structural parameters of our sample, we perform single Sérsic component fits in the image plane using GALFIT (Peng et al., 2002, 2010) on $8\arcsec\times 8\arcsec$ cutouts (from the original resolution, non-PSF matched images) in all filters where the integrated galaxy flux has $\mathrm{SNR}\geq 10$. First, we fit the parameters in F444W, using initial values based on the detection catalog and adopting parameter limits on the Sérsic index ($n_{S}=[0.2,8]$), major axis effective radius ($R_{E,S}=[0.3,300]\,\mathrm{pixels}$), and axis ratio ($b/a=[0.05,1]$), while requiring the total magnitude to be within $\pm 3\,\mathrm{mag}$ and the center position to within $\pm 10$ pixels of the initial values. We simultaneously fit any neighboring objects with $\mathrm{min}(\Delta\mathrm{mag_{X}})<1.5$ (considered over all filters X) and with any part of their segmentation maps falling within $1.5\arcsec$ of the primary object, and mask all other detected objects are using the detection segmentation maps. We then fit in all other filters in which the galaxy has $\mathrm{SNR}\geq 10$, fixing the position and PA of all objects to the F444W values. For all filters, we flag the fit if GALFIT returns a flag $>0$ or if one of the parameters reaches the enforced limit (with the exception that we do not flag fits where $n_{S}=8$ or $b/a=1$). Additionally, as some of our objects are nearly point sources in the reddest filters, we do not flag fits where only $R_{E,S}$ has numerical stability issues and $R_{E,S}<1\,\mathrm{pixel}$.

We modify this procedure if $(b/a)_{\mathrm{F444W}}>0.75$, as the position angle (PA) may be poorly constrained. In this case, we fix the positions to the best-fit $(x_{c},y_{c})_{\mathrm{F444W}}$ and perform an intermediate fit on the shortest-wavelength JWST filter with $\mathrm{SNR}\geq 10$ to determine the galaxy’s PA. After this SW filter fit, we refit F444W with the positions and PA of all objects fixed, and then fit all remaining filters as in the normal procedure. We show example short and long wavelength GALFIT model results for each object (F150W if fit and unflagged, otherwise F200W/F277W; and F444W) in the first three panels of the left and right columns in Fig. 2.

As single Sérsic component fits do not capture galaxy substructures (which are clearly present in many of the filters for galaxies in our sample), we also compute residual-corrected flux profiles and corrected effective radii $R_{E}$ for all filters following Szomoru et al. (2010, 2012). We compute residual flux profiles from the GALFIT model residual images (masking all neighboring objects) using elliptical annular apertures with fixed PA and axis ratio based on each band’s best-fit values. We then determine uncertainties using derived per-band empty aperture scaling relations (to capture correlated noise; as in Skelton et al. 2014). We show these residual-corrected and GALFIT Sérsic fit surface brightness profiles and respective effective radii ($R_{E}$, $R_{E,S}$) for example short and long wavelength bands (F150W or F200W and F444W) of each object in the fourth panel of each column in Fig. 2. We note that the parametric Sérsic $R_{E,S}$ are in fairly good agreement with the corrected $R_{E}$ (median $R_{E}$/$R_{E,S}$ = 1.01) suggesting that to at least first order the single-component fits capture the average light distributions fairly well, even in the presence of the observed substructures.

Because structural fitting is performed in the image plane, we must account for the lensing magnification on sizes. As our sample has only modest magnifications of $\mu\sim 2-4$, we correct for lensing with a factor of $1/\sqrt{\mu}$. Throughout this paper, all physical sizes have been corrected for both cosmology and this lensing factor. Angular sizes, as measured directly in the image plane, are uncorrected for lensing.

As an alternative to the parametric Sérsic index, we also derive the empirical concentration $C$ as defined in Conselice (2003) for all filters. We calculate the concentration from the residual-corrected flux profiles, in order to mitigate the impact of PSF-broadening (compared to measuring $C$ from the observed images).

We also investigate the impact of dust-reddened or dust-obscured AGN on the morphologies of the X-ray detected objects in our sample (7/11). We perform a second, alternative set of morphological fits, using the same methodology as above but adding a point source component. The point source-to-extended flux ratios from these fits are presented in Appendix A, and also discussed in Sec. 4. For the analysis that follows, we exclude objects with appreciable ($\gtrsim 10\%$) estimated point-source AGN flux contribution, and use morphological parameters derived from the single-component Sérsic fits in the analysis.

Our multiwavelength $0.5-4.4\mu\mathrm{m}$ imaging of these 11 ALMA-detected galaxies reveals a diversity of structures. Visual inspection (i.e., Figs. 1 & 2) reveals galaxies that appear to be large disks, both edge-on (9018*, 16790, 26131, 43148; see also Kokorev et al. 2023 for other analysis of 9018*/ID319, which Kokorev et al. report may be in the early phases of a merger with the compact companion to the West given its undisturbed disk morphology) and more face-on (24143*; see also Wu et al. 2023). The sample also includes other smaller galaxies that appear disk-like and edge-on (16840, 24852), a small galaxy that appears to be a moderately inclined disk with a misaligned dust lane but may be irregular (21972), and small objects which are irregular at short wavelengths but appear compact and round at long wavelengths (43086; also 14034, 24823; see Sec. 3, Appendix A, and below).

Quantitatively, the sample also exhibits a wide range structural parameter values. In Fig. 3 we show $R_{E,\mathrm{F444W,kpc}}$ versus inferred stellar masses for our sample. Our sample ranges from low ($\sim 10^{9.5}{\rm{M_{\odot}}}$) to high ($\sim 10^{11}{\rm{M_{\odot}}}$) stellar mass. We find some of the objects visually that are identified as large disks lie above the mass-size distribution of “normal” star-forming galaxies at $z\sim 1-3$ from Suess et al. (2022), while the remainder of edge-on and moderately-inclined disks fall within the scatter of the Suess et al. sample. Unresolved/round objects include one galaxy (with the lowest stellar mass) just barely below the population scatter.

The two most extremely small objects (14034, 24823) host dust-reddened X-ray detected AGN. Based on alternative Sérsic+point source morphological fits, these objects have appreciable ($\gtrsim 10\%$) AGN flux contributions at long wavelengths ($\gtrsim 2\mu\mathrm{m}$; see Sec. 3 and Appendix A), and are therefore excluded from the later structural analysis. The other 5 X-ray-detected objects appear to have fully obscured AGN even at $4.4\mu\mathrm{m}$ (i.e., minimal point source flux contribution), with the HST and JWST-derived morphology are dominated by the host galaxy (see Sec. 3 and Appendix A). These 5 objects are retained in the sample, and the morphological parameters derived from the single-Sérsic fits are adopted for the remainder of this analysis.

In Fig. 4, we consider measured structural parameters (effective radius, axis ratio, and Sérsic index) as a function of wavelength. We find sizes $R_{E}$ spanning $1-2$ orders of magnitude at a given wavelength, axis ratios $b/a$ spanning from $\sim 0.1-1$, and Sérsic indices ranging from $n_{S}\sim 0.2-1$ in HST/ACS F435W to $n_{S}\sim 0.5-8$ in JWST/NIRCam F444W. The three galaxies with the smallest $R_{E}$ in F444W ($\lesssim 0.1\arcsec$) have sufficiently concentrated long-wavelength light distributions (as evidenced by their relatively high $n_{S}$) to exhibit visible diffraction spikes in the log-scale color images (Fig. 1). We do note that the presented structural parameters are sensitive to a number of additional systematic uncertainties. Previous work has shown systematic uncertainties can be important for GALFIT-derived morphological parameters (e.g., van der Wel et al. 2012), likely leading to underestimates of the total uncertainties presented in this work. We attempt to account for some of these systematic uncertainties (due to non-smooth galaxy light profiles) by using residual corrections to derive $R_{E}$, though this does not provide corrections for other, directly-measured structural parameters (i.e., $b/a$, $n_{S}$).

We also observe common trends in the variation of the structural parameters with wavelength, as seen in Fig. 4 (showing both the measured structural parameters, top, and the trends normalized to the values in F444W, bottom). With increasing wavelength, we find $R_{E}$ decreases and $n_{S}$ generally increases, with each galaxy appearing more compact and concentrated at longer wavelengths. (We find a similar trend of increasing $C$ with wavelength, in agreement with the $n_{S}$ trend.) We also find the axis ratios do not drastically change with wavelength (i.e., varying only by $\sim 0.1-0.2$), with the exception of three objects that are extremely compact at long wavelengths (43086; and the two hereafter excluded dust-reddened AGN, 14034, 24823).

Changes in morphology with wavelength, or equivalently color gradients, reflect non-uniform dust distributions as well as relative age or metallicity changes. Color gradients thus provide a powerful probe of how galaxies assemble their stellar mass profiles (e.g., Suess et al. 2021, 2023). In these sub-mm detected galaxies, we infer that dust is a primary driver of the observed color gradients and morphological changes with wavelength, as we will discuss later. Thus to estimate the radial distribution of dust, we calculate a size ratio $\log_{10}(R_{E,\mathrm{LW}}/R_{E,\mathrm{SW}})$ between a long ($>4\mu\mathrm{m}$) and and short wavelength ($\sim 1.5-2.7\mu\mathrm{m}$) band, probing the rest-frame NIR and optical, respectively, for this sample. We use F444W as the LW band, and preferentially use F150W (if the galaxy has sufficient SNR and an unflagged fit in that band), or otherwise F200W or F277W, as the SW band. These size ratios are shown as a function of F444W $R_{E}$ in kpc (corrected for lensing using $1/\sqrt{\mu}$ as in Fig. 3), $n_{S}$, and $C$ in Fig. 5. Our sample has size ratios that range from relatively mild values ($\sim-0.2$), similar to those observed in $z\sim 1-3$ HST-selected star-forming galaxies by Suess et al. (2022), down to steeper ratios ($\sim-0.5$) for the more compact objects in F444W.

While we observe a correlation between $\log_{10}(R_{E,\mathrm{LW}}/R_{E,\mathrm{SW}})$ and $R_{E,\mathrm{F444W}}$, where the smallest objects have the most extreme gradients, we find more clear anti-correlations between $n_{S,\mathrm{F444W}}$ and $C_{\mathrm{F444W}}$ and the size ratio. The Suess et al. (2022) sample also show slight tails of stronger gradients towards lower $R_{E,\mathrm{F444W}}$ and higher $n_{S,\mathrm{F444W}}$, though the observed trend for the sub-mm detected galaxies extends far beyond the region spanned by “typical” SFGs.

We also consider the size ratios versus stellar mass, $\log_{10}(M_{*}/M_{\odot})$ (fourth panel). We find that most (7/9) of our sample (excluding the two X-ray detected objects hosting dust-reddened AGN that contribute to the long-wavelength flux) roughly overlaps with the anti-correlation trend found by Suess et al. (2022), but the remaining two galaxies have lower size ratios at a fixed stellar mass. This suggests that stellar mass alone does not set galaxies’ size ratios and that measures of long wavelength central concentration (i.e., $C_{\mathrm{F444W}}$ or $n_{S,\mathrm{F444W}}$) correlate more strongly with size ratio (in contrast to Gómez-Guijarro et al. 2023, who find dust concentration correlates most strongly with mass).

Finally, we examine the size ratios as a function of the integrated observed F150W–F356W color, which approximates the rest-frame V–J color for our sample redshift ($z\sim 2$; Fig. 5, fifth panel). Again, this is consistent with the anti-correlation of $\log_{10}(R_{E,\mathrm{LW}}/R_{E,\mathrm{SW}})$ with rest-frame V–J found by Suess et al. (2022). Our sample extends that trend to redder objects. In Fig. 5 we show that most (7/9) of our sample continues the $\log_{10}(R_{E,\mathrm{LW}}/R_{E,\mathrm{SW}})$ – V–J trend shown by the Suess et al. contours; however, our galaxy sample also includes two objects that have lower size ratios than expected given their color (as also seen with stellar mass).

Rest-frame UVJ colors have been proposed to help distinguish between color gradients that are driven by dust versus age (e.g., Miller et al., 2022, 2023), though Leja et al. 2019 caution that interpretation of spatially-resolved UVJ colors with respect to dust gradients requires further testing. This analysis relies on the fact that more dust tends to produce redder U–V colors and redder V–J colors, where as older ages tend to produce redder U–V colors but bluer V–J colors. To test our rest-frame color gradients, we use EAzY-Py (Brammer et al. 2008; using tweak_fsps_QSF_12_v3 templates) to interpolate, or for 4 objects, slightly extrapolate, the residual-corrected deconvolved profiles (as in Fig. 2) to construct rest-frame UVJ color profiles. We find that all galaxies in our sample have both redder U–V and V–J colors in the centers than the outskirts, consistent with increasing dust attenuation towards their centers. Combined with the sub-mm detection, the visible signatures consistent with dust lanes or patchy dust distributions in some objects, and the emergence of a stellar bar at long wavelengths for at least one object (43148, which cannot be explained by stellar age gradients), these signatures are consistent with high central dust concentrations being an important driver of the observed color gradients. Accurately detangling the contribution of dust attenuation and age to the observed color profiles will require pixel-to-pixel resolved SED fitting. Future high-resolution sub-mm observations (e.g., with a forthcoming ALMA Cycle 10 program; PI: V. Kokorev) and spatially-resolved spectroscopy will also provide valuable constraints on the dust attenuation and age profiles.

This sample demonstrates that the population of sub-mm detected galaxies is not homogeneous in structure. Among our 11 objects, we observe disk-like structures over a range of inclinations and sizes, as well as smaller, more irregular objects that also host extremely compact central components at the longest wavelengths. In the rest-frame NIR (see Fig. 2), these structures are consistent with stellar disks (with a range of central concentrations; 8/11), small spheroidal or face-on stellar disks (1/11), and NIR continuum dominated by heavily dust-reddened AGN that have appreciable AGN flux contributions at long wavelengths (2/11; excluded from the full analysis).

Even with this morphological diversity, we find evidence that within individual galaxies the dust and NIR continuum light distributions are similar. Specifically, we find the size ratios — and thus color gradients — correlate with size and concentration in F444W. As we infer that the color gradients are primarily driven by dust, based on their sub-mm detection, visible dust lane and patchy distributions, and inferred rest-frame UVJ color profiles (Sec. 4), this correlation suggests that the objects with smaller and more concentrated NIR continuum light (assumed to trace the stars, for all but the 2 excluded dust-reddened X-ray-detected AGN) also have smaller and more compact dust distributions, while objects with less concentrated NIR light also have less concentrated dust. We see further evidence for self-similarity between the NIR continuum and dust distributions from the trend of axis ratios with wavelength: for most objects there is little change in $b/a$ with wavelength, suggesting similar dust geometry to the underlying stellar distribution. For the two excluded dust-reddened X-ray detected AGN, the axis ratio changes reflect the emergence of these point-source AGN as they become less obscured at longer wavelengths, also suggesting distribution self-similarity in order to obscure the central AGN.

Finally, the distribution of 1.2mm flux densities (using updated fluxes from the recently-released DUALZ catalog by Fujimoto et al. 2023b; Fig. 6) versus size ratios suggests that both dust geometry and total dust quantity (as potentially evidenced by sub-mm flux) are important in setting these galaxies’ multiwavelength structures. For fixed dust-star geometry, increasing the total dust will tend to result in stronger color gradients and more extreme short-to-long wavelength size ratios. Increasing the total dust for a fixed distribution (i.e., increasing the dust column densities) will preferentially attenuate the inner regions compared to the outskirts (assuming dust column density $\sim A_{\rm V}$; Salim & Narayanan 2020), increasing the effective radii at short wavelengths while minimally impacting the long-wavelength sizes (assuming the dust is not optically thick at those wavelengths). This would tend to produce an anti-correlation between the sub-mm flux densities and size ratios.

Variations in dust-star (or dust-star-AGN) geometry — including patchiness in the dust distribution — could change the relation between total dust quantity and color gradient, introducing scatter between the sub-mm flux densities and size ratios. Furthermore, if a galaxy is optically thick even in the rest-frame NIR, dust attenuation would also increase the long-wavelength effective radii relative to the stellar half-mass radii. This could result in larger LW-to-SW size ratios (as both the SW and LW effective radii are larger than the stellar half-mass radius), compared to objects where only the SW is attenuated while the LW is relatively unattenuated.

Our sample exhibits both some degree of anti-correlation between $S_{\rm 1.2mm,int}$ and $\log_{10}(R_{E,\mathrm{LW}}/R_{E,\mathrm{SW}})$ and notable scatter (including a strong outlier, 9018*, which may be optically thick even in F444W, potentially explaining the unexpectedly high size ratio given the total sub-mm flux density as discussed above). Taken together, we infer that both total dust quantity and dust geometry are important in explaining size ratios. However, we note that non-dust drivers of color gradients (e.g., age gradients) could also introduce scatter. Future analysis of high-resolution maps of the dust continuum, and spatially-resolved stellar population synthesis modeling, will reveal the relative impact of dust quantity and geometry on multiwavelength morphologies and color gradients, and disentangle the impact of non-dust effects.

Furthermore, high (low) sub-mm fluxes can also reflect high (low) incident radiation from star formation or AGN activity. The uncertain AGN contributions to the 1.2mm flux density for the 5 X-ray detected but highly obscured AGN host galaxies retained in this analysis (24143*, 26131, 43148, 24852, 21972) thus complicates the above inferences regarding dust content and dust geometry. Nonetheless, the non-X-ray-detected objects (9018*, 16790, 16840, 43086) span much of the parameter space that might be expected for both total dust and dust geometry contributions to setting the observed size ratios — suggesting that the conclusion that both dust quantity and geometry are important holds even with the high AGN fraction in this sample. High-resolution dust continuum maps will additionally help address this issue, by spatially detangling the central AGN contribution from the extended emission driven by star formation radiation.

This analysis highlights the power of using high-resolution, deep NIR imaging from JWST to unravel the multi-wavelength, multi-component structures of extremely dusty galaxies, which were previously shrouded from view. More detailed observations and modeling to directly recover the stellar mass and dust distributions are necessary to better understand the structures of sub-mm detected galaxies. This type of spatially-resolved spectral energy distribution (SED) modeling would benefit dramatically from higher-resolution ALMA continuum observations, which are forthcoming in an approved ALMA Cycle 10 program (2023.1.00626.S; PI V. Kokorev). Resolved NIR spectroscopy (either from a slit or IFU) would also be critical to directly probe the reddening and thus disentangle the dust, age, and metallicity color degeneracies. Future work with these observations will also determine how well dust and stellar geometries are captured through multiwavelength structural measurements similar to those presented in this work. Finally, JWST and ALMA will continue to map larger fields, including the recently completed ALMA Cycle 9 program DUALZ covering the full extended Abell 2744 cluster (Fujimoto et al., 2023b), yielding significantly larger samples that fully span the demographics of dust-obscured/sub-mm-detected galaxies across cosmic time. Larger samples would enable full explorations of sub populations and constrain their respective evolutionary pathways. Combined with parallel efforts to characterize structures of less extreme galaxy populations, these studies will continue bridging the gap between optical studies of relatively low dust galaxies and sub-mm studies of the most dusty galaxies.