JWST CEERS probes the role of stellar mass and morphology in obscuring galaxies

In recent years, observations have uncovered a population of massive galaxies that are invisible or very faint in deep optical/near-infrared (near-IR) surveys but brighter at longer wavelengths. However, the nature of these optically dark or faint galaxies (OFGs; one of several names given to these objects) is highly uncertain. In this work, we investigate the drivers of dust attenuation in the JWST era. In particular, we study the role of stellar mass, size, and orientation in obscuring star-forming galaxies (SFGs) at $3<z<7.5$, focusing on the question of why OFGs and similar galaxies are so faint at optical/near-IR wavelengths. We find that stellar mass is the primary proxy for dust attenuation, among the properties studied. Effective radius and axis ratio do not show a clear link with dust attenuation, with the effect of orientation being close to random. However, there is a subset of highly dust attenuated ($A_V>1$, typically) SFGs, of which OFGs are a specific case. For this subset, we find that the key distinctive feature is their compact size (for massive systems with $\log (M_{*}/M_{\odot})>10$); OFGs exhibit a 30% smaller effective radius than the average SFG at the same stellar mass and redshift. On the contrary, OFGs do not exhibit a preference for low axis ratios (i.e., edge-on disks). The results in this work show that stellar mass is the primary proxy for dust attenuation and compact stellar light profiles behind the thick dust columns obscuring typical massive SFGs.


Introduction
Until JWST came online, our understanding of the early (z > 3) cosmic history of galaxies was mainly based on samples selected from their rest-frame ultraviolet (UV) light (see Madau & Dickinson 2014, for a review). However, this view is limited, as classical color selections to find distant star-forming galaxies (SFGs), such as the Lyman break galaxies (LBGs; e.g., Steidel et al. 1996;Giavalisco et al. 2004;Bouwens et al. 2012), are effective at identifying blue galaxies with low dust attenuation, but are generally biased against redder dusty SFGs (see Casey et al. 2014; Hodge & da Cunha 2020, for a review).
In recent years, Spitzer Space Telescope (Spitzer) and Atacama Large Millimeter/submillimeter Array (ALMA) observations from numerous studies have uncovered a population of massive dusty SFGs completely missed in deep optical/nearinfrared (optical/near-IR) surveys, but bright at longer farinfrared/millimeter (far-IR/mm) wavelengths (see e.g., Simpson et al. 2014;Wang et al. 2016Wang et al. , 2019Franco et al. 2018;Alcalde Pampliega et al. 2019;Yamaguchi et al. 2019;Williams et al. 2019;Romano et al. 2020;Toba et al. 2020;Umehata et al. 2020;Zhou et al. 2020;Gruppioni et al. 2020;Talia et al. 2021;Smail et al. 2021;Fudamoto et al. 2021;Manning et al. 2022; Gómez-⋆ Marie Curie Fellow ⋆⋆ NASA Hubble Fellow ⋆⋆⋆ NASA Postdoctoral Fellow Guijarro et al. 2022a;Shu et al. 2022;Enia et al. 2022;Xiao et al. 2023, so-called optically dark/faint galaxies (OFGs), or, alternatively, optically invisible galaxies, optical/near-IR-dark, HSTdark, among others). These dusty SFGs are different from the rare extreme dusty starbursts (e.g., Walter et al. 2012;Riechers et al. 2013;Marrone et al. 2018), as they are more numerous, with number densities two orders of magnitude greater, and are characterized by milder star formation rates (SFRs; e.g., Wang et al. 2019) typical of galaxies in the main sequence of SFGs (MS; e.g., Brinchmann et al. 2004;Daddi et al. 2007;Elbaz et al. 2007;Noeske et al. 2007;Whitaker et al. 2012;Speagle et al. 2014;Schreiber et al. 2015). These galaxies are believed to be a dominant contributor to the massive (M * > 10 10.3 M ⊙ ) end of the SFR density of the Universe at z ∼ 3-6 (e.g., Wang et al. 2019;Xiao et al. 2023; Barrufet et al. 2023). However, the dearth of available secure spectroscopic redshifts and determinations of their stellar mass is a hinderance to the interpretation of these sources among the SFG population. It is still unclear as to how many owe their faintness at optical/near-IR wavelengths to a dusty or a quiescent nature. In addition, being so faint or completely dark to the Hubble Space Telescope (HST) has rendered the study and understanding of the role of their stellar morphologies impossible.
The JWST recently began to be used to investigate this galaxy population, unveiling the stellar light of galaxies previously undetected by HST (HST-dark). In a general sense, it Article number, page 1 of 13 arXiv:2304.08517v2 [astro-ph.GA] 4 Sep 2023 A&A proofs: manuscript no. jwstceers_obs should be noted that there is no consensus definition of this galaxy population and the properties of its members heavily depend on the specific magnitude and color cuts (e.g., Xiao et al. 2023). Among these JWST studies, Barrufet et al. (2023), based on a red (F160W -F444W) color selection, reported their relatively high-redshift, massive, MS-like SFRs, a dusty nature (z ∼ 2-8; M * > 10 10 M ⊙ ; A V ∼ 2 mag), and their important contribution to the cosmic SFR density at z ∼ 6. Nelson et al. (2023) measured a prevalence of low axis ratios and linked their elusiveness to a dusty disk-like and edge-on nature. Pérez-González et al. (2023), based on a red (F150W -F356W) color selection, described this galaxy population as threefold, including a majority (∼ 70%) of dusty SFGs at z = 2-6, with M * = 10 9−10 M ⊙ , ∼ 20% quiescent galaxies (QGs) with M * > 10 10 M ⊙ at z = 3-4, and ∼ 10% young starbursts at z = 6-7 with M * ∼ 10 9.5 M ⊙ . Rodighiero et al. (2023) reported evolved z = 8-13 massive M * = 10 9−10 M ⊙ galaxies with high dust attenuation (A V > 5 mag). Kokorev et al. (2023) presented a spatially resolved study of a z ∼ 2.5 massive (M * ∼ 10 11.3 M ⊙ ) dark galaxy, finding very high levels of dust attenuation (A V ∼ 4 mag) all across its stellar extent.
In this work, we investigate the drivers of dust attenuation in the JWST era, focusing on understanding the nature of obscured galaxies to unveil the reason behind their elusiveness in the pre-JWST era. We study the roles of stellar mass and morphology as proxies for dust attenuation. The layout of the paper is as follows. The data used in this work are described in Sect. 2. In Sect. 3, we present the catalog and photometry, along with the derivation of stellar-based properties, including morphological measurements. Section 4 describes the sample selection, classification into LBGs and OFGs, and evaluation of the star-forming or quiescent nature of OFGs. In Sect. 5, we characterize the stellar-based properties, including redshift, stellar mass, dust attenuation, SFR, and morphological measurements. We investigate the relevance of the stellar mass and morphology in driving dust attenuation in Sect. 6. In Sect. 7, we discuss and interpret our results. We summarize our main findings and conclusions in Sect. 8.

JWST data
In this work, we used JWST data from the Cosmic Evolution Early Released Science (CEERS) survey, an early release science program in the Extended Groth Strip (EGS) field (Finkelstein et al. 2022). The CEERS data include JWST/NIRCam imaging over ten pointings, the first four (CEERS1, CEERS2, CEERS3, CEERS6) observed in June 2022 and the remaining six (CEERS4, CEERS5, CEERS7-10) observed in December 2022, together covering a total area of 97 arcmin 2 . This dataset comprises observations with seven NIRCam filters F115W, F150W, F200W, F277W, F356W, F410M, and F444W, reaching an average 5σ point source depth of 29.15, 29.00, 29.17, 29.19, 29.17, 28.38, and 28.58 mag, respectively. For details about the observations and data reduction process, we refer the reader to Bagley et al. (2023). Briefly, the June pointings were processed using the JWST Calibration Pipeline version 1.7.2 and CRDS pmap 0989, while the December pointings were processed using the pipeline version 1.8.5 and CRDS pmap 1023. For both June and December pointings, the reduction steps were the same, with the raw images passing through all pipeline stages with some additional steps, including 1/ f noise subtraction and wisp removal. The mosaics were aligned to the Gaia-EDR3 astrometry (Gaia Collaboration et al. 2021), mapped into a pixel scale of 0 ′′ . 03/pixel, and background subtracted.

Ancillary data
In order to complement shorter wavelengths, we used HST imaging in the EGS field from its public release reprocessed by the CEERS team, matching the astrometry and pixel scale of the JWST data (see Koekemoer et al. 2011). This dataset comprises observations with six filters: ACS/WFC F606W, F814W, WFC3/IR F105W, F125W, F140W, and F160W, reaching an average 5σ point source depth of 28.62, 28.30, 27.11, 27.31, 26.67, and 27.37 mag, respectively. We also used the publicly available ground-based imaging from the Canada-France-Hawaii Telescope (CFHT)/MegaCam observations in the EGS field in u * , g ′ , r ′ , i ′ , and z ′ -band, reaching average 5σ point source depths of 27.1, 27.3, 27.2, 27.0, and 26.1 mag, respectively. This dataset corresponds to the D3 field of the CFHT Legacy Survey (CFHTLS) (Gwyn 2012).

Source detection and photometry
For source detection, we followed a similar approach to that used in the CANDELS catalog in the EGS field (Stefanon et al. 2017). We employed SExtractor (Bertin & Arnouts 1996) using the F444W-band as a detection image. This band was found to be optimal for our science case, because our goal is to unveil galaxies that become increasingly faint at shorter wavelengths. The F444W-band comprises the longest wavelengths of the JWST/NIRCam bands. We carried out the source detection in both the so-called cold and hot modes, that is, with tailored parameters optimized to detect from bright extended to faint compact sources. We also employed a RMS map (defined as RMS = 1/ √ WHT, where the WHT map is a weight image giving the relative weight of the output pixels) scaled to the median value of the ERR map (data array containing resampled uncertainty estimates, given as standard deviation) as a weighting image for detection. The final merged catalog comprises all cold sources and the hot sources that are not part of a cold source. Discarded hot sources are those that fall within the Kron ellipse of a cold source.
In the 13 JWST and HST bands, we measured fluxes using SExtractor in dual image mode. In this mode, one image is used for detections (F444W-band in our case), while measurements are carried out on another image. Running SExtractor with various measurement images while keeping the same detection image, one ends up with a catalog with the same sources measured through different bands. Measurements were carried out on images PSF-matched to the angular resolution of the detection F444W-band image (0 ′′ . 16 FWHM), which provides the coarsest angular resolution of the JWST/NIRCam bands and a resolution that is comparable to those of the HST WFC3 bands. PSF-matching was performed on JWST/NIRCam and HST/ACS images by convolving them with a convolution kernel created using the software PyPHER (Boucaud et al. 2016), from PSF images of the different bands constructed by stacking point sources that were selected in the mosaics using the software PSFEx (Bertin 2011). Total flux measurements and uncertainties in both the cold and hot modes in the detection F444W-band image were calculated using Kron elliptical apertures (FLUX_AUTO with PHOT_AUTOPARAMS = 2.5, 3.5). For the remaining bands, SExtractor dual mode uses the segmentation map of the detection image to measure isophotal fluxes and uncertainties (FLUX_ISO) that we corrected to total fluxes by multiplying them with the ratio FLUX_AUTO/FLUX_ISO as measured in the detection F444W-band (see Stefanon et al. 2017). HST/WFC3 fluxes and uncertainties were scaled to account for the missing flux when comparing their PSF with the F444Wband PSF. An additional aperture correction optimized for point sources was added in all bands to account for the flux missed outside of the Kron elliptical apertures. For each source, we scaled their fluxes and uncertainties by the fraction of flux outside of the area covered by its Kron ellipse as measured in the F444W-band PSF.
In the five ground-based CFHT bands, we measured fluxes using aperture photometry on images PSF-matched to the angular resolution of the u * -band (0 ′′ . 9 FWHM), which provides the coarsest angular resolution of the CFHT bands (i.e., the poorest seeing conditions at the moment of observation). We measured fluxes in 1 ′′ . 2 diameter apertures, which provides an optimal trade-off between total flux retrieval and signal-to-noise ratio (S/N; see also, e.g., Straatman et al. 2016). Total flux measurements were then calculated using aperture corrections derived tracing the u * -band PSF growth curve to account for the flux losses outside the chosen aperture.
Fluxes and uncertainties in all bands were corrected for Milky Way attenuation following a Cardelli et al. (1989) extinction function with an average E(B − V) = 0.006 in the EGS field (Schlafly & Finkbeiner 2011).

Redshift, stellar mass, and dust attenuation
We estimated photometric redshifts from spectral energy distribution (SED) fitting to the 18 bands in our catalog (u * to F444W) using the code EAZy (Brammer et al. 2008) in its updated Python version EAZy-py. We employed the set of 13 Flexible Stellar Population Synthesis (FSPS; Conroy & Gunn 2010) templates (corr_sfhz_13; see Kokorev et al. 2022, for details). These templates cover redshift-dependent star formation histories (SFHs), which disfavor SFHs that start earlier than the age of the Universe at a given epoch, span a large range in ages and dust attenuation, and include emission lines. Additionally, the set was complemented with the best-fit template to the extreme emission line galaxy 4590 at z = 8.5 with a JWST/NIRSpec spectrum (Carnall et al. 2023) rescaled to match the normalization of the FSPS templates (see e.g., Kokorev et al. 2022;Gould et al. 2023;Valentino et al. 2023, for a similar approach). EAZy-py finds the best linear combination of templates for the flux densities and uncertainties of each source. We searched for updated spectroscopic redshifts by cross-matching the source catalog (r < 0 ′′ . 5) with the ancillary spectroscopic compilation in Stefanon et al. (2017) and the MOSDEF survey in the EGS field (Kriek et al. 2015). When a high-quality measurement was available, we substituted the photometric redshift by its spectroscopic value. In addition, we used the best-quality spectroscopic redshifts to derive zero-point corrections applied iteratively to the EAZy-py run to improve photometric redshifts (N = 320, σ = 0.02, corrections < 10%).
The set of EAZy templates is best suited for redshift determination, but they are limited in terms of their physical parameters. Therefore, we estimated stellar masses and dust attenuation (A V ) using the SED fitting code FAST++ 1 , which is an updated version of the SED fitting code FAST (Kriek et al. 2009), fixing the redshifts to the values previously obtained. The stellar population models were from Bruzual & Charlot (2003) (BC03), with delayed exponentially declining SFHs (τ = 100 Myr-30 Gyr; age = 50 Myr-age of the universe at a given redshift), a Calzetti et al. (2000) dust attenuation law (A V = 0.0-5.0) -which is common in similar studies and provides a benchmark for comparisonand metallicity (Z = 0.004, 0.008, 0.02 (solar), 0.05).

Star formation rate
We estimated SFR values from a ladder of SFR estimators (see e.g., Wuyts et al. 2011;Magnelli et al. 2014;Barro et al. 2019, for a similar approach).
First, we checked whether there are counterparts in the mid-IR-to-mm bands in the "super-deblended" catalog in the EGS field (Le Bail et al., in prep). This catalog is built using a state-of-the-art de-blending methodology similar to that of similar catalogs in the GOODS-North ) and COS-MOS (Jin et al. 2018) fields. These mid-IR-to-mm bands include Spitzer/MIPS/24 µm, Herschel/PACS/100, 160 µm, Herschel/SPIRE/250, 350, 500 µm, JCMT/SCUBA2/450, 850 µm, and ASTE/AzTEC 1.1 mm. For galaxies with S/N IR > 5 (sum in quadrature of the S/N of all the bands beyond and not including Spitzer/MIPS/24 µm), we fit Draine & Li (2007) dust emission templates with the correction from Draine et al. (2014) using the code CIGALE (Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019) to derive total IR luminosity estimates (L IR , 8-1000 µm rest-frame). For galaxies with S/N IR < 5, but S/N 24µm > 5, we renormalized IR templates from Schreiber et al. (2018) to the Spitzer/MIPS/24 µm fluxes, obtaining L IR . The total SFR accounts for the contribution of the obscured star formation probed in the IR (SFR IR ) and the unobscured star formation probed in the UV (SFR UV ), following the prescription of Bell et al. (2005) (for a Chabrier (2003) IMF): where L UV = 1.5νl ν,2800 is the total UV luminosity in the range 1260-3000 Å rest frame derived from the rest frame 2800 Å luminosity density (l ν,2800 ). For galaxies with S/N IR < 5 and S/N 24µm < 5, or for those without counterparts in the mid-IR-to-mm bands in the superdeblended catalog, the total SFR is given by the SFR UV corrected from dust attenuation (SFR UV,cor ). We calculated SFR UV,cor following the same prescription above for the UV term, but with L UV,cor = 1.5νl ν,2800,cor , and with l ν,2800,cor being the rest-frame 2800 Å luminosity density corrected for dust attenuation (A 2800 ) using the Calzetti et al. (2000) law.

Morphology
In order to estimate morphological structural parameters, we used surface brightness profile fitting to the JWST/NIRCam F444W images by employing the code statmorph (Rodriguez-Gomez et al. 2019). The F444W-band is the reddest NIRCam filter available, probing 0.52 < λ rf < 1.1 µm in the redshift range of interest (3 < z < 7.5, see Sect. 4) at ∼ 0 ′′ . 16 angular resolution (PSF FWHM). Therefore, the F444W-band probes optical-tonear-IR rest-frame morphologies, mitigating the morphological differences that can occur at shorter wavelengths, especially in the rest-frame UV (see e.g., Suess et al. 2022;Miller et al. 2023). statmorph performs 2D Sérsic profile fitting and nonparametric morphological diagnostics of a given background-subtracted image and associated segmentation map with the sources of interest. The code includes a flag for the fit quality, where f lag = 0 indicates a good fit and f lag ≥ 2 indicates a problem with the measurements (e.g., artefacts, foreground stars, or incompletely masked secondary sources). Following the documentation, we restricted our analysis to f lag ≤ 1 and f lag_sersic = 0 when considering morphological measurements.

General parent SFGs sample
In order to include new insight from JWST/NIRCam into the rest-frame UV, we introduced a lower limit in the redshift for our study. We focused our work on the redshift range 3 < z < 7.5. The lower redshift limit was imposed to cover the rest-frame UV (λ rf < 4000 Å) at all redshifts with at least two NIRCam bands (F115W and F150W), while the higher redshift limit was imposed because of the scarcity of galaxies at higher redshifts. We also introduced a lower limit in the stellar mass of log(M * /M ⊙ ) ≥ 9.0 to focus our work on intermediate-to-massive galaxies.
We classified this redshift-and stellar mass-limited sample into SFGs and QGs through a UV J classification. The UV J plane is a classical color diagram widely used to make this discrimination (e.g., Wuyts et al. 2007;Williams et al. 2009;Patel et al. 2012). We employed the definition by Williams et al. (2009). The rest-frame colors were obtained from FAST++ and the resulting best-fit SEDs. Although the UV J diagram has proven to be a useful tool for SFGs and QGs discrimination, it is not free of some plausible contamination. In order to mitigate contaminants, we checked that the galaxies classified as QGs by the UV J criteria are undetected (< 2σ) in mid-IR-to-mm bands (from Herschel/PACS/100 µm to ASTE/AzTEC 1.1 mm) in the super-deblended catalog in the EGS field (Le Bail et al., in prep); otherwise we classified them as SFGs.
As a result, we defined a redshift-(3 < z < 7.5) and stellar mass-limited (log(M * /M ⊙ ) ≥ 9.0) sample of UV J-selected SFGs (parent SFGs sample), which serves as a benchmark for comparisons. In total, the parent SFG sample comprises 1490 galaxies.

Lyman break galaxies
The Lyman limit (λ rf = 912 Å) is a signature that allows the identification of high-redshift SFGs. SFGs exhibit a prominent break beyond their Lyman limit owing to absorption by intervening hydrogen. These galaxies are called LBGs. An LBG will appear very faint or invisible in all bands bluewards of its Lyman limit and prominent again redwards of its Lyman limit. The more redshifted the Lyman limit, the more distant the galaxy. Therefore, the Lyman break technique is an efficient way of selecting high-redshift SFGs (e.g., Steidel et al. 1996;Giavalisco et al. 2004;Bouwens et al. 2012). The technique was used up to z ∼ 10 in the pre-JWST era (e.g., Bouwens et al. 2015).
In this work, we label LBGs in the parent SFGs sample galaxies; these are selected by applying the LBG color criteria as defined by Bouwens et al. (2020): faint in all optical and near-IR bands up to and including the Hband in the deepest cosmological fields (typical 5σ point source depth H > 27 mag) but bright at longer near-IR bands ([3.6] or [4.5]-band). Their redshifts, stellar masses, dust attenuation, and star-forming or quiescent nature heavily depend on the specific magnitude and color cuts. A typical selection consists in targeting sources that are Hband dropouts (with studies defining them below a given magnitude limit or as absent in catalogs) but bright at 3.6/4.5 µm, or galaxies with very red colors (e.g., H − 3.6/4.5 > 2-4 mag) (see e.g., Wang et al. 2016Wang et al. , 2019Alcalde Pampliega et al. 2019;Pérez-González et al. 2023). This selection typically yields z > 3 massive and dust obscured galaxies (see also Bisigello et al. 2023, for low-mass dusty SFG selection criteria). Xiao et al. (2023), with the purpose of bridging more extreme optically dark galaxies with more common lower-mass, less-dust attenuated galaxies, such as LBGs, formally define OFGs based on the following criteria: (1) H > 26.5 mag; (2) [4.5] < 25 mag. The H > 26.5 cut selects both the extreme optically dark and intrinsically fainter galaxies with less dust attenuation. This cut also weeds out typical massive (log(M * /M ⊙ ) > 10) QGs. The [4.5] < 25 mag cut selects not only massive galaxies, but also galaxies with intermediate stellar masses (see Xiao et al. 2023, for details).
In this work, we applied the magnitude cuts as defined by Xiao et al. (2023) using the new insight from the JWST/NIRCam bands: (1) F150W > 26.5 mag; (2) F444W < 25 mag. This parent OFG sample comprises 35 galaxies. There are 25/35 (71%) in the parent SFG sample and 3/35 (8.6%) QGs based on the UV J classification. The remaining 7/35 (20%) galaxies belong to a lower-redshift tail (1.5 < z < 3). This tail is left out of the parent SFG sample when we impose the lower redshift limit (z > 3) to cover the rest-frame UV (λ rf < 4000 Å) at all redshifts with at least two NIRCam bands (F115W and F150W). Nevertheless, the OFG criteria appear to be an effective filter for relatively high-redshift (3 < z < 7.5), intermediate-to-massive (M * > 10 9 M ⊙ ), dusty (A V ≳ 1) SFGs (25/35, 71%), with little contamination from QGs (3/35, 8.6%). Hereafter, we continued the analysis focusing on the OFGs sample that is part of the parent SFGs sample (see Sect. 5, for details on their redshift, stellar mass, and dust attenuation distributions). Figure 1 shows the (F150W − F444W) color versus F444W magnitude distribution, including the parent SFGs, LBGs, and OFGs samples. The distribution of LBGs follows the bulk of the parent SFGs population, although we see that the LBG sample starts to miss galaxies as they become redder. The OFG sample is complementary to the LBG sample by selecting bright F444W < 25 mag galaxies missed by the LBG criteria. We note that LBGs and OFGs are not mutually exclusive galaxy types; at the edges of the OFGs selection region, seven galaxies are both LBG and OFG. An example OFG is shown in Fig. 2.
In light of the small percentage of QG contaminants in the OFG sample, we find that most of the z > 3 massive galaxies missed by the LBG criteria in the pre-JWST era were missed mainly because of high dust attenuation rather than because of old stellar populations. In order to illustrate the effect of dust attenuation in the selection criteria for LBGs and OFGs, we performed the following experiment: First, we took the intrinsic (before applying the dust attenuation law) best-fit SEDs of the parent SFG sample from FAST++. We used the SEDs of the parent SFG sample to build a representative set of SEDs, as opposed to taking all the SED templates of FAST++. We then progressively applied increasing A V values following the Calzetti et al. (2000) dust attenuation law. In every A V step, we obtained the magnitudes in the relevant filters for the selection criteria for LBGs and OFGs by convolving their transmission curves with the best-fit SEDs. In each A V step, we calculated the fraction of galaxies that meet the selection criteria for LBGs and OFGs, that is, f LBGs = N LBGs /N SFGs , where N SFG is the number of galaxies in the parent SFG sample; and f OFGs = N OFGs /N SFGs,bright , where N SFG,bright is the number of galaxies in the parent SFG sample with F444W < 25 mag. Figure 3 presents these fractions as a function of increasing A V . This illustrates how the Lyman break technique naturally misses galaxies as a function of increasing A V (∼ 50% for A V = 2 mag and ∼ 0% for A V > 3 mag). On the contrary, the higher the A V , the higher the fraction of OFGs (∼ 50% for A V = 1 mag and ∼ 100% for A V > 5 mag).

Stellar-based properties and morphologies
In this section, we characterize the stellar-based properties, including redshift, stellar mass, dust attenuation, SFR, and morphological measurements, as derived in Sect. 3. We compare the samples defined as in the previous section with regard to these stellar properties. In Fig. 4, we present the redshift, stellar mass, and dust attenuation distributions of the parent SFG, LBG, and OFG samples to illustrate the parameter space covered by these galaxies. All samples exhibit a similar redshift distribution. In terms of stellar mass and dust attenuation distributions, the main difference is in the OFG sample. These galaxies are more massive and dust attenuated than LBGs and the parent sample of SFGs. This reflects the OFG criteria, which are designed to select intermediate to massive galaxies with moderate to high levels of dust attenuation. In Fig. 5, we place our different samples in the stellar mass versus SFR plane, along with the MS of SFGs as defined by Schreiber et al. (2015), assuming that this MS definition can be extrapolated to z > 4. The SFR values are scaled to a common redshift, corresponding to the median value of the parent SFGs sample (z med = 3.95). The scaled SFR values do not represent the true SFR values, but the overall distribution of the different samples in the SFR-M * plane, maintaining the ∆MS = SFR/SFR MS each galaxy would have if plotted against the MS associated with its redshift. The parent SFG sample has a median ∆MS = 0.91 +0.74 −0.40 (where the uncertainties are the 16% and 84% percentiles of the distributions), tracing the MS of SFGs. Similarly, LBGs have ∆MS = 0.95 +0.67 −0.37 and therefore trace the bulk of the MS of SFGs. In the case of OFGs, ∆MS = 0.38 +0.66 −0.19 locates them slightly below the MS of SFGs. Previous studies revealed that the SFR of the OFGs (or comparable selections) is consistent with the typical MS SFGs as constrained from stacked mid-IR-to-mm photometry or deep mm individual detections from ALMA (e.g., Wang et al. 2019;Gómez-Guijarro et al. 2022b;Xiao et al. 2023). SFR could be underestimated in OFGs in the absence of counterparts in the mid-IR-to-mm that typically offer the best SFR constraints. In this work, none of the OFGs have counterparts in the mid-IR-to-mm bands in the super-deblended catalog in the EGS field, and so SFRs would be underestimated, likely because of an underestimated dust attenuation, which already reaches saturation in the SED fitting of such highly dustattenuated systems. The extreme cases are the OFGs that are still undetected in their rest-frame UV by JWST (9/25), for which dust attenuation constraints are heavily dependent on extrapolations of the best-fit SED. These galaxies do not show differences in terms of redshift when compared to the OFG sample (∼ 5% difference in the medians) and exhibit 1.7 times higher stellar masses and 1.2 times higher dust attenuation. Pérez-González et al. (2023) reported a similar SFR behaviour compared to the MS for HST-dark galaxies, which is mitigated when performing 2D SED-fitting.
In terms of morphological parameters, in Fig. 6, we present the (circularized) effective radius (R e ), axis ratio (q = b/a), and Sérsic index (n) distributions of the parent SFG, LBG, and OFG samples (flagged as good fits). We recall that these were measured in the JWST/NIRCam F444W-band and that we restricted our analysis to f lag ≤ 1 and f lag_sersic = 0 when considering morphological measurements. Among the OFGs sample, 6/25 galaxies were flagged as poor fits, in all cases being found to be consistent with point sources. Therefore, while it was not possible in these cases to measure their morphologies, a limit on their effective radii is PSF FHWM /2 < 0 ′′ . 08. At first order, the distributions are rather similar in the different samples. The main difference is in the lower scatter around small effective radii and Sérsic index n ∼ 1 in the case of the OFGs sample compared to LBGs and the parent SFGs sample. In addition, the distribution of axis ratios appears skewed toward larger values in the case of OFGs. While not displayed here, the upper limits on the sizes for the 6/25 OFGs classified as point sources are also consistent with these sources being clustered around smaller effective radii. Their compact sizes could also explain the skewness toward large axis ratios, as the latter are more difficult to constrain when the angular resolution starts to be comparable with the intrinsic size of the galaxy (marginally resolved sources).

Drivers of dust attenuation
In this section, we study the relevance of stellar mass and morphology in driving dust attenuation in a general sense and, specifically, in the case of highly dust attenuated galaxies, such as OFGs. We note that drivers in this case refer to parameters that show correlations, and these can therefore be interpreted as predictors of the dust attenuation. For a discussion on the relevant physical processes involved, we refer to Sect. 7. Figure 7 shows dust attenuation as a function of stellar mass. Dust attenuation correlates with stellar mass in SFGs, as reported before in numerous studies (e.g., Garn & Best 2010;Zahid et al. 2013;Heinis et al. 2014;Pannella et al. 2015;Álvarez-Márquez et al. 2016). The parent sample of SFGs establishes the general trend. LBGs exhibit a shallower trend, while OFGs display a steeper trend populating a much more dust attenuated regime.

Stellar mass
The dust attenuation and stellar mass correlation arises from the mass-metallicity relation (e.g., Tremonti et al. 2004;Erb et al. 2006;Sánchez et al. 2013;Genzel et al. 2015;Ma et al. 2015). Massive SFGs are capable of more efficiently producing and retaining their metals than low-mass SFGs, which tend to expel them in galactic winds. Therefore, massive SFGs tend to be more dusty than low-mass SFGs. Stellar mass acts as a primary proxy for dust attenuation in SFGs. However, it remains unclear as to whether or not additional parameters can drive dust attenuation. In the following section, we study the role of morphology.

Morphology
We studied whether dust attenuation shows signs of correlation with basic morphological structural parameters. In Fig. 7, we show dust attenuation as a function of effective radius and axis ratio. The behavior of A V for the parent sample of SFGs and LBGs is generally flat across the range of effective radii and axis ratios. The main difference arises in the case of OFG sizes. The effective radii of OFGs are clustered around small sizes (as shown in Fig. 6) that are associated to high dust-attenuation values. In addition, we see an increase in A V when moving toward smaller effective radii. On the contrary, we do not see an increase in A V toward lower axis ratios (i.e., edge-on galaxies).
Generally, there is no clear link -as in the one involving stellar mass-between the studied basic structural parameters and dust attenuation. However, galaxy size appears to be related in the most dust attenuated galaxy type, the OFG sample. The question that rises is whether or not OFGs are indeed smaller when compared to the typical SFG sizes at similar stellar mass and redshift. Figure 8 shows the ratio between median (circularized) effective radius of the different galaxy types and the median (circularized) effective radius of the parent sample of SFGs. For better statistics, we focus on a single redshift bin (3 < z < 5) of massive galaxies (log(M * /M ⊙ ) > 10). In addition, we also move be- yond the selection criteria for LBGs and OFGs by simply defining galaxies with higher dust attenuation (A V > 1 sample), of which OFGs would be specific cases, and galaxies with lower dust attenuation (A V < 1 sample), which would comprise the majority of LBGs. The A V > 1 and OFG samples are indeed smaller than the parent sample of SFGs at the same redshift and stellar mass bin. In the case of axis ratios, these appear generally consistent with the values of the parent sample of SFGs, albeit slightly skewed toward larger values, which is consistent with the picture in which axis ratios are more difficult to constrain when the angular resolution starts to be comparable with the intrinsic size of the galaxy. This affects more galaxy types that are intrinsically smaller, such as the A V > 1 and OFG samples and, generally, galaxies at higher redshifts following the decreasing sizes of SFGs with increasing redshift (e.g., van der Wel et al. 2014). We note that the upper limits in the sizes for the 6/25 OFGs classified as point sources are not included here. If included, they would lower the size ratio even further, in line with A V > 1 and OFG samples exhibiting smaller sizes than the parent sample of SFGs.  2015) displayed as a solid blue line. Its 1σ scatter associated with 0.5 < ∆MS < 2 (∼ 0.3 dex) is represented as a shaded blue area, with a more extended typical scatter of 0.33 < ∆MS < 3 (∼ 0.5 dex) in lighter blue. We note that the values are scaled to a common redshift (z med = 3.95) as explained in the main text.
While we studied a given galaxy type compared with the parent sample of SFGs in the same redshift and stellar mass bin, there are still differences in the redshift and stellar mass distributions. In Fig. 8, measurements (filled squares) are placed at the median stellar mass of the galaxy type. These galaxy types do not have the exact same median stellar mass (or redshift) between them and they do not have the exact same median stellar mass (or redshift) as the parent SFG sample. This could have an impact on the comparison. However, we know from literature studies that SFG sizes decrease with increasing redshift and increase with increasing stellar mass (e.g., Shen et al. 2003;van der Wel et al. 2014;Barro et al. 2017;Suess et al. 2019;Mowla et al. 2019). While these studies were carried out in the pre-JWST era, and therefore they are typically limited to z < 3 for stellar sizes in the rest-frame optical/near-IR, if we assume that the same trends can be extrapolated to higher redshifts, we can correct for the differential effect that redshift and stellar mass have on the sizes (see also Gómez-Guijarro et al. 2022a, for a similar analysis). Following this approach, we corrected the effective radii of each galaxy to a common redshift and stellar mass as given by the median values of the 3 < z < 5 massive (log(M * med /M ⊙ ) > 10) parent sample of SFGs (open squares in Fig. 8). After this exercise, A V > 1 and OFG samples are a factor 1.3 smaller than the parent sample of SFGs and a factor 1.5 smaller than the A V < 1 sample at the same redshift and stellar mass bin.

Random forest analysis
The results above indicate that stellar mass is a primary proxy for dust attenuation. In a general sense, effective radius and axis ratio do not play a clear role in dust attenuation. However, when considering massive galaxies at similar redshifts, highly dust-attenuated galaxies exhibit more compact (smaller effective radii) stellar light profiles. We check these results by performing a random forest analysis (see e.g., Ellison et al. 2020;Baker et al. 2022;Bluck et al. 2022, for a similar analysis).
A random forest is a machine learning technique that identifies the most important parameter in driving a given quantity, especially when several parameters are inter-correlated. A target quantity is selected and removed from the data, leaving a dataset containing the features that the target quantity imprinted on it. These two datasets (target and features) are then split into a training and a test sample. The algorithm is employed in the training sample, which sorts the data in different nodes into the different decision trees with the goal being to minimize the Gini Impurity (Pedregosa et al. 2011). The final model is then applied to the test sample in order to check its performance.
We used a random forest regressor in order to identify the most predictive parameter of A V between M * , SFR, R e , and q. These parameters are chosen to reflect different aspects of galaxies, going from the total amount of stars (M * ), the rate of creating new stars (SFR), the extent of the stellar light R e , and the way the stellar light is oriented (q). We also added a random uniform variable (Rand) as a control variable. Using a randomized cross-search validation method we fine-tuned the hyperparameters. We checked the performance of the fit by calculating the mean squared errors (MSEs) between model and test sample. Figure 9 shows the relative importance of the studied parameters in predicting A V from the random forest regressor. Stellar mass is the most predictive parameter, followed by SFR. The morphological parameters of effective radius and axis ratio do not appear as strong predictors of dust attenuation, with the axis ratio being close to random. These results are in line with the results shown in Sect. 6.

Discussion
The results in this work show that stellar mass is a primary predictor of dust attenuation. The correlation between dust attenuation and stellar mass was studied in the pre-JWST era, especially at z < 4 using statistical galaxy samples (e.g., Garn & Best 2010;Zahid et al. 2013;Heinis et al. 2014;Pannella et al. 2015;Álvarez-Márquez et al. 2016), with some studies peering into the z > 4 Universe using smaller samples (e.g., Fudamoto et al. 2020;Bogdanoska & Burgarella 2020). In the JWST era, some works started to investigate the relation between dust attenuation and stellar mass at z > 3 (e.g., Shapley et al. 2023), with models predicting similar trends (e.g., Yung et al. 2019). Stellar mass is a tracer of the integrated history of star formation in galaxies and therefore the history of metal and dust production. The stellar mass defines the depth of the galaxy potential well and therefore massive galaxies are able to keep their metals and dust more efficiently than low-mass galaxies, which are prone to lose them through galactic winds. The mass-metallicity relation is strong observational proof of these physical processes (e.g., Tremonti et al. 2004;Erb et al. 2006;Sánchez et al. 2013;Genzel et al. 2015;Ma et al. 2015).
Once the stellar mass has been accounted for, SFR appears as another strong predictor of dust attenuation. The two are correlated in SFGs, as shown by the existence of the MS of SFGs (e.g., Brinchmann et al. 2004;Daddi et al. 2007;Elbaz et al. 2007;Noeske et al. 2007;Whitaker et al. 2012;Speagle et al. 2014;Schreiber et al. 2015). At fixed stellar mass, a number of studies in the literature explain the dust attenuation scatter through variations in the SFR (e.g., Garn & Best 2010;Wuyts et al. 2011;Zahid et al. 2013). However, one interesting aspect of the random forest analysis in Sect. 6.3 is that stellar mass appears as a stronger predictor of dust attenuation than SFR. While the mass-metallicity relation can be considered a first-order relation in galaxy evolution, the fundamental metallicity relation C. Gómez-Guijarro et al.: JWST CEERS probes the role of stellar mass and morphology in obscuring galaxies (FMR; Mannucci et al. 2010) comes as a second-order indicator of the metal content in galaxies. The FMR indicates that, at fixed stellar mass, galaxies with higher SFR have lower metallicities. Therefore, the FMR would favor stellar mass as a stronger predictor of dust attenuation than SFR.
Orientation has been studied in the literature as a proxy for dust attenuation. Disk galaxies can exhibit higher dust attenuation levels for higher inclinations, owing to the thicker projected dust columns (e.g., Wild et al. 2011;Patel et al. 2012;Zuckerman et al. 2021). Zuckerman et al. (2021) explains the dust attenuation scatter at fixed stellar mass as due to inclination effects in thin dust disks, but also shows that when the dust disk is thick, inclination becomes irrelevant. The results presented in this work show that at z > 3, inclination has a negligible effect on dust attenuation, and is almost indistinguishable from a random variable. Highly dust attenuated systems in this work do not exhibit a preference for more inclined orientations (i.e., edge-on disks). In particular, OFGs are not biased toward edge-on disks, contrary to the interpretation of Nelson et al. (2023), who links the obscuration of similar galaxies to an edge-on nature. A very recent study by Lorenz et al. (2023) also reaches the same conclusion at lower redshifts, showing that none of the dust properties studied (Balmer decrement, A V , rest-frame UV slope) vary with axis ratio in a sample of 308 spectroscopically characterized SFGs at 1.3 < z < 2.6. The authors propose an interpretation in which already at z > 1 the dust attenuation in SFGs would be dominated by dust in dense star-forming regions rather than by the dust in the diffuse ISM, which would mean a negligible effect of inclination on dust attenuation. Some examples of dust swept in dense clouds and depleted diffuse components have also been observed in the local Universe (Holwerda & Keel 2013;Keel et al. 2014).
While stellar mass and SFR come as strong predictors of dust attenuation for the parent sample of SFGs, our results show that there is a subset of galaxies that, even at fixed stellar mass, SFR, and redshift, are fainter than the average SFG at rest-frame UV wavelengths with higher levels of dust attenuation (A V > 1). OFGs are a specific case of these galaxies, and exhibit a factor 2.4 higher A V than the parent sample of SFGs (for 3 < z < 5 and log(M * /M ⊙ ) > 10). We find that the key distinctive feature is their compact size (smaller effective radii) compared to the average SFG at the same stellar mass and epoch (as given by the parent sample of SFGs). The physical interpretation of this is that the galaxy size is the parameter that modulates how the density profile of the galaxy is distributed. More compact stellar profiles M * SFR R e q Rand Relative importance 0.00 0.25 0.50 Fig. 9. Relative importance of the different parameters explored in predicting dust attenuation from a random forest regressor analysis. The bar height displays the median value of the distribution of 10 000 realizations of training and test samples, while the error bars represent the 16% and 84% percentiles of the distributions.
could therefore contribute to more efficiently keeping metals and dust in a galaxy compared to a more extended stellar profile. For galaxies with similar levels of dust mass (determined at first order by their redshift, stellar mass, and SFR), the dust column density would be thicker when embedded in more compact stellar profiles than in more extended stellar profiles.
More compact stellar profiles in OFGs are in line with the view of ALMA-detected OFGs; they exhibit enhanced SFR surface densities associated with compact dust morphologies with high dust covering fractions. Therefore, these galaxies are capable of forming stars in concentrated areas -which would result in the compact stellar light profiles we see in this work-and of obscuring the stellar light behind thick dust columns (e.g., Smail et al. 2021;Xiao et al. 2023;Kokorev et al. 2023). The compact stellar sizes of OFGs are comparable to those of QGs. Massive (log(M * /M ⊙ ) > 10) OFGs in the redshift bin 3 < z < 4 have a median (circularized) effective radius of R e = 0.75 ± 0.25 kpc (where the uncertainty is given by the median absolute deviation). As a reference, QGs at z = 3 of similar stellar mass have a typical effective radius of R e = 0.73 kpc (van der Wel et al. 2014). In the classical view of galaxy structural evolution, SFGs evolve mostly as extended star-forming disks, while QGs are typically more compact than SFGs at a fixed stellar mass and redshift (e.g., Shen et al. 2003;van der Wel et al. 2014;Barro et al. 2017;Suess et al. 2019;Mowla et al. 2019). The buildup of a central stellar core appears to be linked to the quenching of star formation (e.g., Kauffmann et al. 2003;Whitaker et al. 2017;Barro et al. 2017;Gómez-Guijarro et al. 2019;Suess et al. 2021). SFGs with compact star forming regions (e.g., Toft et al. 2014;Gómez-Guijarro et al. 2018, 2022bFranco et al. 2020;Puglisi et al. 2021;Magnelli et al. 2023) and SFGs with a developed compact stellar core (e.g., Barro et al. 2013;Nelson et al. 2014;Williams et al. 2014;van Dokkum et al. 2015;Gómez-Guijarro et al. 2019) have been proposed as the link between the more extended SFGs and the more compact QGs. The compact sizes (smaller effective radius) in OFGs when compared to the general SFG population shown in this work and the higher SFR surface densities in OFGs shown in literature studies are in line with a progressive buildup of the compact stellar profiles of QGs. Therefore, OFGs (at a higher redshift bin) could be progenitors of massive compact QGs (at a lower redshift bin).

Summary and conclusions
In this work, we investigate the drivers of dust attenuation in the JWST era, with a special focus on understanding why the galaxies known as optically faint/dark galaxies (OFGs) are so faint at optical/near-IR wavelengths. To this end, we use data from JWST to unveil their nature and characterize their basic stellar and morphological properties. Our findings can be summarized as follows: -Stellar mass is a primary proxy of dust attenuation, among the properties studied. Effective radius and axis ratio do not show a clear link with dust attenuation, with the effect of orientation being found to be close to random. -There is a subset of highly dust attenuated (A V > 1) SFGs, of which OFGs are a specific case. Even at the same stellar mass, SFR, and redshift, this subset of galaxies are still fainter than the average SFG at rest-frame UV wavelengths with high levels of dust attenuation. We find that their key distinctive feature is their compact size, with 30% smaller effective radius than the average SFG at the same stellar mass and redshift (for massive systems with log(M * /M ⊙ ) > 10). These galaxies do not show a preference for low axis ratios (i.e., edge-on disks). -The OFG selection criteria (F150W > 26.5 mag; F444W < 25 mag) are an effective filter for relatively high-redshift (3 < z < 7.5), intermediate to massive (log(M * /M ⊙ ) ≥ 9.0), dusty (A V ≳ 1) SFGs, with little contamination from QGs (∼ 8.6%). Compared to the general SFG population, OFGs exhibit smaller effective radii, with Sérsic indices of n ∼ 1 and axis ratios of around q ∼ 0.5.
The results in this work show that stellar mass is a primary proxy for dust attenuation and compact stellar light profiles behind the thick dust columns obscuring typical massive SFGs.