CEERS Key Paper. IV. A Triality in the Nature of HST-dark Galaxies

Previous works have shown that searching for massive, dusty, and quiescent galaxies at z ≳ 3 using color–magnitude diagrams is a very effective technique (Mancini et al. 2009; Huang et al. 2011; Caputi et al. 2012; Labbé et al. 2013; Nayyeri et al. 2014; Stefanon et al. 2015; Wang et al. 2016; Alcalde Pampliega et al. 2019). Briefly, these methods use a long-baseline optical/NIR color to identify galaxies with very red SEDs, which would be caused either by strong Balmer breaks/D4000 spectral features or large, dust-driven UV attenuations, which are more common in evolved, massive galaxies. Mid-to-near-IR large colors might also be present in galaxies with strong (high equivalent width) emission lines (see, e.g., Smit et al. 2014, 2015; Roberts-Borsani et al. 2016; Alcalde Pampliega et al. 2019). Prior to JWST, these methods typically relied on the HST/F160W filter for the blue band and either Spitzer/IRAC [3.6] or [4.5] channels for the red (and selection) band. The resulting samples find galaxies with relatively bright IRAC magnitudes ([3.6] and or [4.5] ≤ 24.5 mag) that are often very faint, even undetected, in the H band, implying F160W > 26–27 mag (depending on the observed field). Thus, they are commonly referred to as IRAC-bright, HST-dark (or -drops) galaxies (Caputi et al. 2012; Stefanon et al. 2015; Alcalde Pampliega et al. 2019).

The new and extremely deep near- and mid-IR first data sets from JWST reach more than 2 mag deeper in F150W and F356W than the previous HST and Spitzer photometry and, consequently, it has become easier to identify and characterize both IR-bright and much fainter HST-dark galaxies (Barrufet et al. 2022; Endsley et al. 2022; Nelson et al. 2022; Labbe et al. 2022).

Here we aim to study the properties of optical/near-IR faint galaxies selected with a similar method to those in previous papers but taking advantage of the new JWST capabilities. Therefore, we use the following thresholds in color–magnitude: F150W − F356W > 1.5 mag, F356W < 27.5 mag. The red color is similar to the values used in previous works (e.g., Wang et al. 2016; Alcalde Pampliega et al. 2019), while the much deeper limiting magnitude in F356W, compared to IRAC depths, [3.6] ∼ 24.5 mag, implies that we can potentially find massive galaxies ($\mathrm{log}{M}_{\star }/{M}_{\odot }\,\gt$ 11) up to z = 8 without including too many lower-redshift galaxies with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lt$ 10, therefore keeping the overall sample focused on massive galaxies.

In order to understand better its heterogeneous nature, we divide the color-selected sample into two subsets, with magnitudes brighter and fainter than F150W = 26 mag, labeled HST-faint and HST-dark, respectively. The latter are a proxy for HST/WFC3 dropouts that were previously only detected in IRAC and/or IR/radio surveys, or completely unknown. The HST-faint sample, which consists mainly of HST/F160W-detected galaxies identified before in CANDELS and 3D-HST catalogs (Skelton et al. 2014; Stefanon et al. 2017), will be used for completeness, reference, and comparison to illustrate how the formal definition of an HST/WFC3-dropout (which depends on the depth of the survey and changes from paper to paper) affects the median redshifts and stellar population properties of an HST-dark sample (see next section). We do not consider here galaxies brighter than F150W < 23 mag, which are typically massive galaxies at z ≲ 2 and have been well-characterized by HST and ground-based telescopes, even with spectroscopy.

Figure 1 shows our selection method. The left panel depicts the F150W − F356W versus F356W color–magnitude plot for the galaxies in the CEERS sample using a gray-scale hexdiagram to indicate the regions with a higher number density. The selection of sources and integrated photometry used in this plot have been performed with a hot+cold Sextractor run (Bertin & Arnouts 1996), using Kron (1980) apertures. The details of the catalog will be presented in S. Finkelstein et al. (2023, in preparation).

Zoom In Zoom Out Reset image size

The color threshold used to select the sample is indicated with a dashed gray line, and the resulting sample of HST-faint and -dark galaxies is shown with circles and hexagons (light and dark colors, respectively). For comparison purposes, we show the sample of HST-dark galaxies from Alcalde Pampliega et al. (2019), selected in the CANDELS/GOODS region, with green circles. The overall distribution of those galaxies agrees very well with the location of our HST-dark sample at relatively bright F356W ≲ 25 mag, except for the three bright and very red galaxies in the Alcalde Pampliega et al. (2019) sample, which have no equivalent in our JWST-based selection.

To illustrate in more detail the types of galaxies and redshift ranges selected by the color–magnitude diagram and also highlight the advantages and shortcomings of this method, the middle panel of Figure 1 shows the tracks in the color–magnitude diagram of six different galaxy templates as a function of redshift. These galaxies are chosen to be (broadly) representative of the typical types in different regions of the UVJ diagram (top-right panel; Williams et al. 2009; Whitaker et al. 2010; Patel et al. 2012; Leja et al. 2019b; Zuckerman et al. 2021), i.e., they range from blue, low-mass, and low-attenuation galaxies (bottom-left region of the UVJ diagram; blue colors) to redder, more massive, dusty star-forming galaxies (SFGs; top right; purple colors) and quiescent galaxies (top left region of the UVJ diagram; red/orange colors). The average SEDs for these galaxies are shown in the bottom-right panel of Figure 1. The SEDs are normalized at ∼1.6 μm rest frame to highlight the reddening of the UV/optical emission with increasing age and/or dust attenuation (from blue to purple templates).

In Figure 1, the thin/thick line-width color–redshift tracks show the effect of scaling the templates to low and high stellar masses, $\mathrm{log}{M}_{\star }/{M}_{\odot }=11$ and 10, respectively, which are roughly traced by the F356W magnitude (with some degeneracies and scatter). The low-mass/extinction templates are shown only for $\mathrm{log}{M}_{\star }/{M}_{\odot }=10$, since high dust contents typically exist at the high-mass end.

Overall, the color–redshift tracks illustrate that red colors F150W − F356W ≳ 1.5 mag and bright F356W magnitude provide an excellent proxy to identify massive, dusty SFGs at redshifts above z ≳ 2 and massive, quiescent galaxies at slightly higher redshifts above z ≳ 3. We note also that some of the bluest quiescent galaxies (U − V ∼ 1.3 mag, orange color track) are still slightly bluer than F150W − F356W = 1.5 mag at z = 3, and were missed by color thresholds of previous HST+Spitzer-based papers. For this paper, considering composite stellar population models with some additional residual recent star formation and not adding much mass (and thus not affecting the F356W flux significantly), we use a color cut F150W − F356W = 1.5 mag.

The color–redshift tracks also show that using a limiting magnitude down to F356W ∼ 27 mag includes many more galaxies at higher redshifts and/or lower masses than previous IRAC-bright searches, which are capped at [3.6] ∼ 24.5 mag (and in relatively large sky areas), often picked the most massive galaxies ($\mathrm{log}{M}_{\star }/{M}_{\odot }\,\sim$ 11) at redshifts of z ≳ 3–4 (Wang et al. 2016; Alcalde Pampliega et al. 2019). Indeed, the $\mathrm{log}{M}_{\star }/{M}_{\odot }=10$ color tracks (thick lines) indicate that the selection would be roughly complete down to that mass for quiescent galaxies up to z ≲ 8 (red/orange) and even for the dustiest galaxies up to z ∼ 6 (purple/magenta). Similarly, at the high-mass end ($\mathrm{log}{M}_{\star }/{M}_{\odot }\,\gt$ 11; thin lines), the color–magnitude selection would identify all the quiescent or dusty galaxies up to z ∼ 6. We note that the top-right corner of the figure, where we would expect the most massive dusty galaxies at z > 5, is nearly empty. This is likely due to the relatively low number density of such sources, and very massive galaxies in general (e.g., Muzzin et al. 2013; Stefanon et al. 2015) combined with the small area surveyed by the current CEERS pointings (∼5 × smaller than for example the CANDELS/GOODS regions). We further discuss these implications in Section 6.2.

One potential concern for mass completeness would be the existence of massive galaxies at z > 5–7 with less canonical (composite) SEDs, such as the ones identified by Labbe et al. (2022). Those galaxies have mirror-flipped-L-shaped SEDs (V-shaped in f_λ ) with red colors in the rest-frame near-IR (V − J ∼ 2 mag) and blue colors in the rest-frame UV similar to those of our blue templates (black tracks and marker in Figure 1). The strong UV upturn and relatively red optical SED suggest that these are experiencing a recent burst of unobscured star formation but also have a preexisting more evolved or dust-enshrouded stellar population. However, the optical emission could be partially contaminated by the presence of high-equivalent-width emission lines, which makes the estimate of the stellar mass more uncertain (Endsley et al. 2022). The color tracks for these galaxies decline very quickly down to F150W − F356W ∼ 1 mag at z ≳ 5 instead of plateauing at around F150W − F356W ∼ 2.5 mag slightly under the quiescent and dusty tracks. As we will demonstrate in the following sections, we have several of these types of galaxies in our sample thanks to our color cut at F150W − F356W = 1.5 mag, but still, some more active systems could exist and be missed by our selection. If those galaxies are common at the highest redshifts, they will be systematically missed by the typical color–magnitude selections unless we decrease the color cut (at the expense of sample contamination by less extreme systems).

Our color–magnitude-selected sample is composed of 138 galaxies, whose median/quartiles magnitude and colors are $\langle {\rm{F}}356{\rm{W}}\rangle ={24.2}_{23.0}^{25.6}$ mag and $\langle {\rm{F}}150{\rm{W}}-{\rm{F}}356{\rm{W}}\rangle ={1.9}_{1.6}^{2.3}$ mag. Half of the sample, 50% (36%), 69 (49) galaxies, are HST-dark, F150W > 26(27) mag, the rest would qualify as HST-faint.

Among the HST-dark subsample, the fraction of sources included in HST catalogs of the EGS area (Skelton et al. 2014; Stefanon et al. 2017) is less than 28%, and all of these galaxies are brighter than F160W = 27 mag. This number goes up to 35% if we include an IRAC selection (Barro et al. 2011). Among the HST-faint subsample, 85% of galaxies were previously cataloged by HST-selected and/or IRAC-selected surveys.

Out of the 138 galaxies in our sample, 32 of them (23%) can be identified with MIPS 24 μm sources in the catalog presented in Barro et al. (2019), searching in a 15 radius region, 75% of them being in the HST-faint sample. One-sixth of them have MIPS 70 μm measurements, all with a signal-to-noise ratio (S/N) < 2. A visual inspection of the MIPS24 sources reveals that 11 of them (75% in the HST-faint sample) are very bright and isolated MIPS emitters, with an average flux of 67 μJy detected at the 9σ level and average redshift 〈z〉 = 2.8. The rest of the counterparts, 21 galaxies (65% in the HST-faint sample), are quite faint and/or possibly blended with nearby objects, presenting an average flux of 50 μJy detected at the 7σ level and average redshift 〈z〉 = 3.2.

Herschel detections (all above 4σ) by PACS at 100 μm were found for seven galaxies (5% of the sample), average flux 3.1 mJy detected on average at the 7σ level, with two of them in the HST-dark sample. For PACS 160 μm, we found six galaxies with S/N > 4 measurements, average flux 7.5 mJy detected on average at the 7σ level, and one in the HST-dark sample. For SPIRE bands, we found three possible counterparts with S/N > 4, average flux 13 mJy detected at, on average, 4σ level, and one in the HST-dark sample.

A total of 10 galaxies in our sample (7% of the total) are likely associated with SCUBA-2 sources reported in Zavala et al. (2017, we note that this survey does not cover our whole area). Nine of them are unequivocally associated with the dusty galaxies via 1.1 mm high-angular-resolution observations obtained with NOEMA (Zavala et al. 2023; L. Ciesla et al. in preparation), while the remaining four sources lie within the 15 of the 850 μm position. An additional three sources were found in the VLA 20 cm maps from Ivison et al. (2007), within a 15 radius.

Finally, we also cross-correlated our sources with the X-ray catalog presented in Nandra et al. (2015) finding two (three) association(s) within 05 (16).

The IDs, coordinates, and basic information about our sample are given in Table 1.

Table 1. Sample of Mid‐IR Bright Near‐IR Faint Galaxies Presented in this Paper

ID	R.A. (J2000)	Decl. (J2000)	F356W	F150W − F356W	z	$\mathrm{log}{M}_{\star }$	logSFR(10 Myr)	Comment
	(degrees)	(degrees)	(mag)	(mag)		(M⊙)	(M⊙ yr−1)
nircam1-349	214.95571202	+52.98342646	27.31 ± 0.06	>2.77	${6.24}_{-0.02}^{+0.03}$	9.01 ± 0.09	−0.29 ± 0.10	XELG-z6
nircam1-1085	214.95787499	+52.98030148	23.31 ± 0.03	1.91 ± 0.04	${3.52}_{-0.06}^{+0.04}$	10.46 ± 0.03	−0.76 ± 0.03	QG
nircam1-1623	214.92577403	+52.95443735	24.79 ± 0.03	2.11 ± 0.05	${3.07}_{-0.13}^{+0.06}$	9.82 ± 0.03	1.47 ± 0.06	SFG
nircam1-1744	215.01127115	+53.01358970	24.79 ± 0.03	1.78 ± 0.04	${5.17}_{-0.00}^{+0.00}$	9.96 ± 0.06	2.79 ± 0.05	SFG
nircam1-1752	215.01221476	+53.01469376	24.72 ± 0.03	1.56 ± 0.04	${3.62}_{-0.13}^{+0.03}$	9.80 ± 0.03	1.53 ± 0.03	SFG
nircam1-1867	214.99840491	+53.00461910	25.95 ± 0.03	2.03 ± 0.06	${6.52}_{-0.08}^{+0.03}$	9.63 ± 0.10	2.07 ± 0.04	XELG-z6
nircam1-2080	214.98181546	+52.99123394	22.32 ± 0.03	1.72 ± 0.04	${3.17}_{-0.24}^{+0.18}$	10.69 ± 0.01	0.52 ± 0.06	QG

Note. Basic information about the sample of galaxies in this paper: ID, coordinates, magnitude and colors used in the selection, redshift, stellar mass, SFR (averaged in the last 10 Myr), and comments (including galaxy activity type based on SED-based sSFR estimations (see text for details; we also mark quiescence results based on the UVJ diagram), and detection by MIPS, submillimeter, and/or X-ray surveys).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Photometric redshifts for each pixel were estimated with a modified version of the eazy code (Brammer et al. 2008). The modification consisted in allowing the template-fitting algorithm to use (5σ) upper limits as an input, not allowing any fit to provide brighter fluxes (achieved by penalizing the χ² calculation) than those limits, and excluding the band in the χ² calculation for templates providing lower fluxes (Mérida et al. 2023). The code was run with the template error feature disabled (which penalized IRAC bands as a default; i.e., NIRCam LW fluxes would also be affected) and not using any prior. We used v1.3 templates, which include a dusty galaxy and a high-EW emission-line galaxy spectrum.

After obtaining photometric redshifts for each pixel, we combined the probability distribution functions (zPDF) for all pixels belonging to a given galaxy using the procedure described in Dahlen et al. (2013). Briefly, different zPDFs are combined allowing for catastrophic fits as well as some degree of interdependency among the different estimations, parameterized with a quantity called α, which can take values from 1 to the number of combined zPDFs. We note that, in our case, the estimation of photometric redshifts for nearby pixels should be correlated, since the FWHM of the PSF-matched data is 016 (5 pixels), but for most galaxies, we counted with completely independent estimations given their large extension. We made some tests on the parameter used to account for this interdependence of estimations and found no significant variations in the final most probable redshifts provided for each galaxy for α > 2.

The procedure of photometric redshift estimation in two dimensions is exemplified in Figure 2 with a peculiar source, as explained next. We show one of the selected galaxies, nircam2-7122, which seems to have two stellar clumps in the JWST data. The northern part, orange in the RGB image, was included in HST catalogs (e.g., in Stefanon et al. 2017). At wavelengths redder than 2 μm, another clump starts to be visible to the south, a region that gets brighter at longer wavelengths while the northern knot gets relatively dimmer (and almost disappears by 4 μm). Our 2D analysis indicates that the northern (already known) galaxy lies at a lower redshift, as revealed by the photo-z map, which counts with 4–5 resolution elements for the northern galaxy and 6–7 times more for the new JWST source. SEDs and photo-z fits are shown in Figure 2 for the centers of the two distinct galaxies. A similar analysis was performed in the rest of the sample, just keeping the regions that present consistent photometric redshifts in the study of the stellar populations presented in the following sections.

Our approach to analyzing the photo-z maps was that in order to consider possible contamination from an interloper at a different redshift, the region must have a size of at least two resolution elements (defined in terms of the FWHM). This allows us to get rid of noise, visible in some zones of the photo-z map presented in Figure 2. In that plot, we clearly identify a different galaxy or knot with this criterion to the north (but not in other zones). We must also consider that the map is presented in terms of the most probable photo-z, but each pixel counts with a complete zPDF.

With this in mind, in order to assign redshifts to different regions of a given galaxy, we considered a probability threshold to decide when two zPDFs (constructed with several pixels as explained above) provide incompatible redshifts. This assumes that a zPDF would catch degeneracies in the photo-z determination linked to stellar population gradients throughout the galaxy (or belonging to different galaxies). We considered as separate galaxies those regions where the combined zPDF implied a probability smaller than 10% of the region being at the redshift of the rest of the galaxy.

We remark that the identification of redshift interlopers in this work implies that some regions of a given source were removed from the subsequent analysis (as in Figure 11). This means that the total masses could be underestimated if the removed zone is indeed part of the galaxy, but the SFHs (or other properties such as SFRs, sSFR, A(V), etc...) would be robust for the retained galaxy area.

The evaluation and calibration of the 2D photometric redshift determination, i.e., how to combine the information from different pixels to obtain the best redshift for a given galaxy as well as help with the segmentation of objects, would need more extended data sets in terms of area and availability of accurate measurements at higher redshifts than what current spectroscopic samples provide (and would be extended in the coming months with JWST observations).

As a first evaluation, we run the procedure through all the galaxies with spectroscopic redshifts in the four available CEERS pointings (see Section 2), just using the nine filters selected for our study of HST-dark/faint galaxies. These galaxies sum up 1.1 million pixels, for which we calculated photometric redshifts. We obtained σ_NMAD values (as defined in Brammer et al. 2008) a factor of 2 better (0.035 versus 0.018) than using the integrated SEDs.

Once we obtained a most probable redshift by combining the estimations for all pixels belonging to each galaxy, we fixed its value and analyzed the stellar populations. We fitted the pixel-by-pixel and integrated SEDs with the synthesizer code (Pérez-González et al. 2003; Pérez-González et al. 2008b) assuming a delayed exponential as the SFH, with timescale values τ between 100 Myr and 5 Gyr, ages between 1 Myr and the age of the universe at the redshift of the source, all discrete metallicities provided by the Bruzual & Charlot (2003) models, and a Calzetti et al. (2000) attenuation law with V-band extinction values, A(V) between 0 and 5 mag. We assume a universal initial mass function (IMF) that follows the Chabrier (2003) functional form. Nebular continuum and emission lines were added to the models as described in Pérez-González et al. (2008b). A Monte Carlo method was carried out in order to obtain uncertainties and account for degeneracies (see Domínguez Sánchez et al. 2016). We warn the reader that this kind of methodology would not account for systematic uncertainties such as those arising from the parameterization of the SFH, metallicity evolution, or the complexity of dust effects (in general, and as a function of position in the galaxy).

Figure 3 shows the results of the stellar population synthesis (SPS) modeling for the galaxy presented in Figure 2 after removing the northern component, which is revealed to lie at a different redshift by our 2D photo-z method. All relevant quantities have been corrected by the aperture correction mentioned in Section 4. The galaxy presents a knot with significant amounts of dust, A(V) ∼ 3 mag, separated from the stellar-mass centroid, and (integrated) SFR around 10 M_⊙ yr⁻¹. The rest of the galaxy presents a less active and low-dust-content disky structure (see also the UVJ diagram per pixel). The SFH of the galaxy shows three main episodes of star formation, one beyond z ∼ 10, another (extended) at z ∼ 5.5, and the recent burst in the center. The galaxy is not detected by MIPS, PACS, SPIRE, or (sub)millimeter wavelengths (nor in X-ray) at the 3σ level or better. This galaxy is transiting through a poststarburst phase with some residual dusty star formation near the nuclear region, presenting mass-weighted ages of up to a few hundred Myr.

Zoom In Zoom Out Reset image size

The 2D UVJ diagram on the bottom right of Figure 3 shows that the full galaxy would be located in an intermediate region between blue, low-dust content, SFGs and dusty starbursts. The pixels near the mass-weighted center of the galaxy lie toward the dust-enshrouded zone, while outer pixels beyond the effective radius tend to lie in the blue cloud and poststarburst/quiescent wedge. In order to demonstrate the reliability of this spatially resolved UVJ diagram, we show color maps directly constructed with the JWST bands that roughly probe the UVJ bands, namely, F150W, F277W, and F444W. For the latter, given that we would need some extrapolation to probe the J band, we multiply it by a factor of 1.5, so we roughly match in the color map the V − J values in the color–color plot.

Our 2D-SPS methodology was compared in García-Argumánez et al. (2022) with fits of the integrated emission of simulated galaxies in Illustris using SFHs including a single burst and two star formation events, each following a delayed exponentially decaying SFH. We demonstrated that this more classical approach could underestimate the mass-weighted ages by factors of a few (see Figure 8 in García-Argumánez et al. 2022). The paper also showed that our 2D-SPS approach could reproduce the Illustris SFHs with 20%–30% accuracies in physical properties such as mass-weighted ages and formation times. For this paper, we compared our 2D-SPS-based stellar masses with those obtained from single-burst delayed exponentially decaying SFHs fitting the integrated emission, obtaining a very small systematic offset, ${\rm{\Delta }}\mathrm{log}{M}_{\star }=0.02$ dex, and a scatter that ranges from 0.05 dex at $\mathrm{log}{M}_{\star }/{M}_{\odot }\gt 10$ to 0.3 dex for smaller masses. A detailed comparison of modeling procedures including 2D-SPS and integrated photometry with parametric and nonparametric SFHs will be presented in other papers (e.g., A. García-Argumánez et al. 2023, in preparation).

The left panel of Figure 4 shows the photometric redshift versus stellar-mass diagram for galaxies in the HST WFC3/F160W-selected CANDELS catalog (gray-scale density map) overlapping with the CEERS region, and the galaxy sample selected by our JWST color–magnitude criterion. We also distinguish between JWST galaxies that are brighter and fainter than F150W = 26 mag, namely, HST-faint and -dark galaxies (light and dark colors, hexagons and circles), respectively. At first glance, the main difference in the distribution of galaxies selected with HST and the new JWST-selected sample is that the latter reveals a population of red, massive ($\mathrm{log}{M}_{\star }/{M}_{\odot }\,=$ 9–11) galaxies at z ≳ 2, many of which were previously undetected (HST-dark) even in the deepest HST/WFC3 surveys such as CANDELS.

Zoom In Zoom Out Reset image size

In more detail, the HST-faint sample selects predominantly massive galaxies ($\mathrm{log}{M}_{\star }/{M}_{\odot }\gtrsim 10$) over a redshift range around z ∼ 3. This is consistent with the distribution of points for HST-faint galaxies in the color–magnitude diagram presented in Figure 1. As illustrated by the template color tracks, similarly massive galaxies at redshifts z ≲ 2 are not in the HST-faint sample because they have bluer colors (F150W − F356W < 1.5 mag) and/or brighter F150W mag, while massive galaxies at higher redshifts tend to be redder and fainter, and thus fall into the selection region of the HST-dark sample.

Consequently, the HST-dark sample exhibits a higher median redshift z ∼ 4.4 and a broader distribution that extends up to z ∼ 7. The redshift overlap between HST-faint and -dark samples at z ∼ 3.5 (central panel of Figure 4) is a direct consequence of the magnitude limit used to define the two sets (F150W = 26 mag). If we choose a fainter magnitude, many of the HST-dark galaxies at higher-z will move to the HST-faint sample, increasing its median redshift and lowering the median redshift of the HST-dark sample. We note that while F150W = 26 mag is an adequate limit to identify WFC3 dropouts in the shallower CANDELS mosaics in EGS, the deeper WFC3 imaging in GOODS/UDF fields (see, e.g., Guo et al. 2013; Barro et al. 2019) means that HST-dark galaxies in those fields will be fainter in F160W (by 0.5 to 1 mag) and will have, on average, higher redshifts (as presented exactly in those fields in, e.g., Alcalde Pampliega et al. 2019).

The redshift histogram for HST-dark galaxies reveals a peak around z ∼ 6.5. The sources in this peak are peculiar, as we will also show in the coming sections. Here we emphasize that these sources typically present a blue F356W − F410M color, on average −0.3 ± 0.2 mag, compared with an average of +0.2 ± 0.2 mag for the rest of the sample. Some of them are also brighter in the F410M compared to both F356W and F444W.

In terms of the stellar-mass distribution (right panel of Figure 4), the HST-dark sample has a lower average stellar mass than the HST-faint sample, $\mathrm{log}{M}_{\star }/{M}_{\odot }\sim 9.8$ versus $\mathrm{log}{M}_{\star }/{M}_{\odot }\sim \,10.4$, although both distributions have galaxies with masses up to $\mathrm{log}{M}_{\star }/{M}_{\odot }\sim 11.0$. The HST-dark sample also extends to lower stellar masses than the HST-faint sample. This is primarily due to the faint limiting magnitude of the selection in F356W, which, to first order, correlates with the stellar mass of the galaxies. As illustrated by the template color tracks in Figure 1, some of the galaxies in the HST-dark sample with magnitudes around F356W ∼ 27 mag will have $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 10$ at redshifts ranging from z = 3 to 7.

Another consequence of having a fainter limiting magnitude is that the average stellar mass of our HST-dark sample is also lower than those of pre-JWST HST-dark studies, which selected galaxies using Spitzer/IRAC and thus were restricted to much brighter limiting magnitudes of [3.6] ∼24.0–24.5 mag. Therefore, these samples had virtually no galaxies with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 10$. For example, Wang et al. (2016) or Alcalde Pampliega et al. (2019) find average stellar masses of the order of $\mathrm{log}{M}_{\star }/{M}_{\odot }=10.5$ and 10.8, respectively, for their HST-dark, IRAC-bright samples. A different effect that can contribute to the larger stellar masses of these IRAC-based studies is that they cover a larger area than the current CEERS survey by up to a factor of 4–5 and, therefore, they are more likely to find a higher fraction of rare galaxies with very large stellar masses. Furthermore, the number of those galaxies detected over the relatively small area surveyed so far by CEERS is more susceptible to cosmic variance.

In this section, we study the star formation properties and rest-frame colors of our sample in the context of the star formation main sequence and the UVJ diagram. Based on this analysis, we will divide the sample into star formation activity subsets that will be explored in the upcoming sections to analyze their stacked properties (e.g., star formation or mass-assembly histories). In particular, we will classify the sample into three classes: (1) quiescent sources; (2) SFGs, many of them detected at far-infrared and submillimeter wavelengths (note that the coverage of the CEERS footprint by submillimeter surveys is very inhomogeneous, so this is not a complete sample); and (3) peculiar galaxies at z ∼ 6. We refer the reader to Section 6.2.3 for more details.

6.2.1. Recent Star Formation Properties and the Star-forming Main Sequence

Figure 5 shows the SFR versus stellar-mass plot for our sample, known as the star formation main-sequence (SFMS) diagram. In order to further understand the properties of our sample, which derive directly from the selection method, we divide the sample into several redshift bins in the following plots. The selected redshift bins are: one dominated by HST-faint galaxies at 2 < z < 3, one dominated by HST-dark galaxies at 4 < z < 7, and an intermediate range with a similar number of galaxies from the two subsamples, 3 < z < 4. SFRs have been derived by calculating the average value from the SFHs of each galaxy in different time intervals, ranging from the last 5 to the last 100 Myr (namely, 5, 10, 25, 50, 75, and 100 Myr). In this figure, we use the SFR defined as that averaged over the last 10 Myr. Similar trends are observed in plots using other definitions of the SFR.

Zoom In Zoom Out Reset image size

In Figure 5, we see our sample forming a band in the SFR–stellar-mass plane around the SFMS; most galaxies would qualify as SFGs and a small fraction presents very small specific SFRs (a 0.1 Gyr⁻¹ limit is marked in the plot), expediting the classification as quiescent/dormant galaxy. Although our SFR and stellar-mass estimations do follow a sequence, many of the galaxies lie below the "classical" SFMS (this is especially clear when using SFRs averaged on the last 75–100 Myr), typically obtained by estimating stellar masses from integrated-aperture SED fitting and SFRs from observed rest-frame UV luminosities corrected for attenuation with recipes based on UV slopes, or with SFR tracers linked to dust emission. At first sight, this kind of behavior (which has also been observed when comparing the observed SFMS with galaxy formation simulations such as Illustris; see Sparre et al. 2015; Donnari et al. 2019) might be interpreted as the SED-fitting analysis underestimating the SFRs due to high dust contents.

To probe the possible existence of highly embedded, optically thick star-forming knots that cannot be characterized with rest-frame UV-to-near-IR SEDs and to test how the 2D approach handles this classical problem, we discuss the effects of dust on the SFMS plot in Figure 6. For this purpose, we calculate the total energy absorbed by dust according to our pixel-by-pixel SPS analysis, adding up the contributions of all pixels for each galaxy to calculate the stellar dust-absorbed luminosity. Assuming an energy balance, this luminosity should be converted into dust emission in the mid- to far-IR. We then translate the dust-absorbed stellar luminosity to an IR-based SFR, applying the calibration in Kennicutt (1998). We remark that this exercise should provide results that are completely comparable to the regular technique used to construct the SFMS plot at the high-mass end. For low-mass galaxies, with lower dust contents, SFRs are typically based on rest-frame UV luminosities, maybe added to absorbed emission with hybrid SFR tracers (see Buat et al. 2005; Calzetti et al. 2007) or dust-corrected using UV slopes (with more or less complicated recipes, see Barro et al. 2019). Therefore, for our galaxies, which span a relatively wide range of masses, the SFRs calculated with an energy-budget technique might be underestimated at the low stellar-mass end if significant amounts of nonobscured star formation exist.

Zoom In Zoom Out Reset image size

Figure 6 demonstrates that the 2D-based, global dust attenuations for the galaxies in our sample do reproduce the dust content and emission observed in galaxies and are used for the construction of the SFMS. Even though the 5–100 Myr time-averaged SFRs fall below the MS, the fair comparison with observations behaves much better. On average, and down to $\mathrm{log}{M}_{\star }/{M}_{\odot }\sim 9.7$, our energy-budget-based SFRs reproduce the SFMS, with the median sSFR lying within 0.1 dex of the literature curves. However, the scatter of our SFMS is 2× larger than the typical 0.2–0.3 dex reported in the literature. Whether this behavior is intrinsic or related to observational uncertainties is beyond the scope of this paper.

In Figure 6 (and the right panel of Figure 5), we mark the galaxies in the redshift peak at z ∼ 6 reported in the previous section (green triangles). All of these sources present very large SFRs for their stellar mass, well above the SFMS, consequently qualifying as starburst galaxies. The SFRs derived with the energy-budget calculation are also large, indicating the presence of significant amounts of dust. We remark here that the starburst galaxies seem to follow a separate sequence, very similar to the bimodality reported in Rinaldi et al. (2022). We cannot claim any robust interpretation for this effect due to the small number of galaxies we have in our sample, but our modeling points to systematic effects linked to age (rather, for example, than for extinction) since most of the galaxies located far away from the SFMS present very young bursts, many of them belonging to the XELG-z6 type. This would be consistent with the presence of strong emission lines, most probably [O iii]+Hβ, in these starburst systems, similar to the sources reported in Endsley et al. (2022), Matthee et al. (2022), and Laporte et al. (2022), but they also present relatively high dust contents, with average attenuations around A(V) = 2 mag. Indeed, Matthee et al. (2022) report not negligible amounts of dust in some of the (NIRCam-grism-selected) spectroscopically confirmed [O iii] emitters.

We conclude that our 2D-SPS method based on JWST+HST data does provide robust dust estimations at spatial resolutions around 200 pc and that SFRs from SED fitting should be compared with SFRs from classical tracers with caution.

6.2.2. Rest-frame UVJ Colors

Figure 7 shows the overall distribution of all the CEERS galaxies (gray-scale density map) and the color-selected sample, divided into HST-faint and -dark in the rest-frame UVJ diagrams at different redshifts. The distribution of rest-frame colors is roughly consistent with what we see in the SFMS plot: The bulk of the sample is massive, dusty SFGs with a small fraction of quiescent galaxies starting at z ≳ 3. The majority of the galaxies identified as quiescent by sSFR (red squares) are located in the quiescent region of the UVJ diagram, showing that the two methods are quite consistent.

Zoom In Zoom Out Reset image size

At 2 < z < 3, nearly all the color-selected galaxies are HST-faint, dusty and star-forming, with high average extinctions 〈A_V 〉 = 2.2 mag and relatively large stellar masses of $\langle \mathrm{log}{M}_{\star }/{M}_{\odot }\rangle =10.4$. This is consistent with the color tracks in Figure 1 that indicate that only dusty galaxies are redder than the selection threshold at z < 3. That is, quiescent galaxies at these redshifts do not make the color cut. Similarly, since these galaxies are also relatively massive, they tend to be on the brighter end of both the F356W and F150W mag, and consequently, the majority are HST-faint. Interestingly, many of these galaxies are located in the very high dust-obscuration region of the UVJ diagram with rest-frame colors exceeding V − J > 2.0 mag. Previous works have pointed out that this dusty corner of the UVJ diagram seems to exhibit a higher fraction of objects with late-type morphologies (or low Sérsic 1968 indices) and closer to edge-on inclinations that can increase the reddening along the line of sight due to dust within the disks, which lead to values of A_V ≳ 2 mag (e.g., Patel et al. 2012; Wang et al. 2018).

At 3 < z < 4, there is roughly an even split between HST-dark and -faint galaxies. The average stellar mass is roughly the same, but the average extinction is much lower, 〈A_V 〉 = 1.4 mag, than at z ∼ 2.5. Similarly, the average V − J ∼ 1.5 mag is smaller and fewer galaxies are found in the dusty corner V − J > 2 mag. Notable exceptions are the submillimeter/radio galaxies with V − J > 1.5 mag. The lower average attenuation and V − J color are likely driven by the inclusion of a larger number of quiescent galaxies (∼17%), which enter the color selection at z > 3. We note also that at z ∼ 3.5 some of the younger, perhaps recently quenched, quiescent galaxies with ages ∼1 Gyr have bluer U − V ∼ 1.3 mag colors, and some of them will not make the color threshold of the selection, F150W − F356W > 1.5 mag, even at z = 3 (orange template in Figure 1). To emphasize this point, in Figure 7 we show 6 out of the 17 galaxies in the recent Carnall et al. (2023) paper, which finds massive quiescent galaxies at z > 3 in the CEERS field. Three out of these six galaxies (green circles) are in our UVJ quiescent region but are not selected by our color criterion, as they are too blue in F150W − F356W. The other 11 galaxies with redder colors and higher redshifts are included in our sample and 9 of them are also classified as quiescent by our criteria.

At 4 < z < 7, nearly all selected galaxies are HST-dark. The average extinction and color continue their downward trend relative to the lower-redshift bins with 〈A_V 〉 = 1.2 mag, 〈V − J〉 = 1.2 mag, but there is also a noticeable decrease in the average stellar mass to $\langle \mathrm{log}{M}_{\star }/{M}_{\odot }\rangle =10$. Since stellar mass and dust are strongly correlated (e.g., Fang et al. 2018), a decrease in the average mass will also lower the attenuation. The lower average mass is caused by a decrease in the number of very massive (and dusty) galaxies (e.g., $\mathrm{log}{M}_{\star }/{M}_{\odot }\gtrsim 10.8\mbox{--}11$) in our sample at redshift z ≳ 4 (see Figure 4). This is most likely caused by a combination of the steep decline in the number density of massive galaxies with redshift and our small surveyed area. For example, Alcalde Pampliega et al. (2019; see also Muzzin et al. 2013; Stefanon et al. 2015) find densities of $\mathrm{log}(N)\sim -4.75$ or −5.00 [Mpc⁻³] at 3 < z < 6, even after accounting for HST-dark galaxies. In the small area covered by CEERS, that density translates to less than ∼1 galaxy (we find none), and even at $\mathrm{log}{M}_{\star }/{M}_{\odot }\gt 10$ the number of galaxies is expected to be of the order of ∼20 based on pre-JWST studies (e.g., extrapolating Muzzin et al. 2013, $\mathrm{log}(N)\sim -3.75$ [Mpc⁻³] at z ∼ 3.5). Consequently, more massive, dusty SFGs like those identified in submillimeter surveys (Wang et al. 2019) will likely be found in larger area surveys. In fact, two (of the five) submillimeter NOEMA galaxies at this redshift are among the most massive and higher-extinction galaxies. In particular, galaxy (a) was also recently discussed in Zavala et al. (2023) as a source that might be misidentified with very high-z (z ≳ 10) candidates.

Aside from the average trends, the fraction of quiescent galaxies is a bit smaller compared to the previous redshift bin, ∼10%, but we note an increase in the overall scatter such that there are more galaxies with more extreme V − J colors (e.g., V − J > 2 or V − J < 0.25). The majority of these galaxies are all at z > 6, in the redshift peak reported in previous sections. They also have peculiar SEDs, in terms of F356W − F410M color, and SEDs that are red at λ > 2 μm but get bluer at shorter wavelengths, similar to the SEDs of the galaxies presented in Labbe et al. (2022). At those redshifts, the V − J color is not probed by the observed SED, and therefore, it is extrapolated from the best-fit SED template. Most of these galaxies are typically a combination of a dusty stellar population (hence the red V − J) plus a low mass, un-extinguished burst that causes the blue UV rest-frame color and a high, recent SFR (see SFH section). The true stellar properties of these galaxies are unclear and still debated, most remarkably whether the red color is caused predominantly by strong emission lines. Our 2D SED fits reveal strong (i.e., high equivalent width, EW) [O iii]λ4959,5007 and Hβ emission boosting the F356W flux (and sometimes the F140M emission), linked to a very young starburst on top of an older population. For this reason, we will identify these sources as extreme emission-like galaxies at z ∼ 6 (XELG-z6, hereafter). They are a relatively numerous population; in fact, they dominate our sample at z > 6, and we might still be missing some, since their colors decrease with redshift, ranging from F150W − F356W = 3 mag down to our limit 1.5 mag (see postage stamps in Figure 8). In fact, the majority of galaxies from Labbe et al. (2022) at z > 7 are bluer and fainter at F356W, so we only have three in common with their sample. Given the selection of our sample, we are biased toward galaxies with emission lines entering the F356W passband, but if similar high-EW ELGs exist at lower redshifts, the lines could lie on other filters and affect the interpretation of JWST data. This is the case for one of the objects discussed in Zavala et al. (2023), for which z = 4.6 and z = 16.7 (Donnan et al. 2023) redshift solutions are obtained, depending strongly on the usage of templates with high-EW emission lines. In our 2D study, that exact galaxy, with ID nircam2-2159 in our catalog, is assigned with a redshift z = 4.59, with all its 74 pixels preferring the low-redshift solution over the z ∼ 16 value. Apart from that, and interestingly, all XELGs-z6 share identical morphologies/appearances (see Section 6.3).

Zoom In Zoom Out Reset image size

6.2.3. Selecting Subsets of Galaxies Based on Star Formation Activity

6.3.1. Visual Appearances

6.3.2. Morphological Parametric Measurements

Figure 9 shows the average global SFHs derived for the galaxies in our sample divided by star formation activity, namely, QGs, SFGs, and XELGs-z6. Global SFHs were computed by summing up the measurements from all pixels in a galaxy (each one fitted with a delayed exponentially decaying parameterization). This SFH was scaled up to reproduce the total stellar mass of the object by multiplying by the average (using all photometric bands) flux aperture correction obtained from the comparison of the sum of pixels arriving at our 2D stellar population analysis and the integrated-light elliptical aperture. After that, we transformed all SFHs to the common time frame of the age of the universe. Then, we averaged the functions in 100 Myr time intervals using all galaxies lying at a redshift below that corresponding to the given universe age. In order to make the average SFH curves more meaningful (in terms of combining coetaneous galaxies with similar natures), we divided the SFG and QG samples into two regimes, one at lower redshifts, z < 4 and z < 3.5, respectively, and one at higher, with roughly the same number of galaxies in each bin.

Zoom In Zoom Out Reset image size

The average SFH of our sample (black line) is roughly constant, at a 4 M_⊙ yr⁻¹ level, from z ∼ 30 down to z ∼ 3, when star formation bursts occur peaking 3–10 times above the steadier previous state. This translates to forming around 10¹⁰ M_⊙ of stars above z > 3, in 2 Gyr, and an extra 30% of that mass in more recent starbursts. Our sample is thus dominated by (dusty) a population of SFGs that smoothly (when averaging their properties) assemble their mass until there is an intense star-forming event with high dust content.

If we restrict the analysis to SFGs, a similar SFH is observed for the z < 4 sample (which dominates the sample), with a slightly larger average level. For z > 4 SFGs, we see a steady state at a similar level, 5 M_⊙ yr⁻¹, followed by a starburst at z = 4–6, with peak SFRs six times larger than the previous plateau. We note that these starbursts are averaged out when considering the full sample, which implies that they must be short and stochastic from galaxy to galaxy.

The combined interpretation for the SFHs we derive for SFGs point to a relatively constant SFR for around 1 Gyr broken by a first major starburst event beyond z = 4–6 (maybe only in some galaxies) and another one below z < 3. The SMGs in our sample are caught in the first burst, implying that the first stage with less extreme SFRs and lasting for 1 Gyr was capable of producing enough amounts of metals to feed the dusty starburst revealed by the submillimeter/millimeter observations.

The left panel of Figure 9 also shows the SFH for the QGs in our sample, again divided by redshift (z < 3.5 and z > 3.5 to have similar numbers of galaxies in each subsample). For the z < 3.5 QGs, after some constant activity at a 10–15 M_⊙ yr⁻¹ level, larger than the value observed for other galaxy types, we see a very prominent star formation peak, reaching five times larger SFRs values, around 30 M_⊙ yr⁻¹, occurring at z = 5–7, at a similar epoch when high-z SFGs present their first peak. Then, the SFH decays very rapidly (a factor of 4 in 500 Myr) and is established at the 5–10 M_⊙ yr⁻¹ level for 1 more Gyr.

For high-z, QGs the SFH peak occurs even before, z ∼ 7–10, it is not as marked, decays slightly more rapidly (300–400 Myr), and reaches almost null activity afterwards. We note that the high-z and low-z subsamples of QGs have different median masses (see Table 2), which explain the average levels of the SFH in Figure 9 (also true for other subsamples).

Table 2. Statistical Properties of the Stellar Populations in the Different Types of Galaxies Presented in the Paper

Sample	z	$\mathrm{log}{M}_{\star }$	$\mathrm{log}\,\mathrm{SFR}(5\,\mathrm{Myr})$	$\mathrm{log}\,\mathrm{SFR}(25\,\mathrm{Myr})$	$\mathrm{log}\,\mathrm{SFR}(100\,\mathrm{Myr})$	$\mathrm{log}\,\mathrm{SFR}({\rm{E}}-\mathrm{budget})$	A(V)	Agem−w	z5%	$\mathrm{log}{{\rm{\Sigma }}}_{{M}_{\star }}$
		(M⊙)	(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)	(M⊙ yr−1)	(mag)	(Gyr)		(M⊙ kpc−2)
ALL	${3.7}_{3.1}^{4.6}$	${10.2}_{9.7}^{10.6}$	${1.4}_{0.3}^{2.1}$	${0.7}_{-0.1}^{1.4}$	${0.5}_{-0.1}^{1.2}$	${1.5}_{0.8}^{2.1}$	${0.7}_{0.2}^{1.5}$	${0.77}_{0.39}^{1.04}$	${12.0}_{9.0}^{15.2}$	${8.4}_{8.0}^{8.8}$
HST-faint	${3.3}_{2.5}^{3.7}$	${10.4}_{10.2}^{10.7}$	${1.3}_{0.6}^{2.1}$	${0.8}_{0.4}^{1.4}$	${0.7}_{0.3}^{1.3}$	${1.6}_{1.1}^{2.2}$	${0.7}_{0.2}^{1.4}$	${0.77}_{0.46}^{1.11}$	${10.5}_{7.8}^{13.5}$	${8.4}_{7.9}^{8.8}$
HST-dark	${4.4}_{3.7}^{5.7}$	${9.8}_{9.4}^{10.2}$	${1.4}_{0.2}^{2.1}$	${0.3}_{-0.8}^{1.2}$	${0.2}_{-0.7}^{0.9}$	${1.4}_{0.5}^{1.9}$	${0.9}_{0.3}^{1.7}$	${0.65}_{0.25}^{0.88}$	${13.5}_{10.1}^{16.0}$	${8.5}_{8.2}^{8.9}$
SFG	${3.6}_{2.9}^{4.3}$	${10.2}_{9.8}^{10.6}$	${1.5}_{0.9}^{2.1}$	${0.9}_{0.4}^{1.5}$	${0.7}_{0.3}^{1.3}$	${1.7}_{1.2}^{2.2}$	${0.9}_{0.3}^{1.6}$	${0.71}_{0.33}^{1.04}$	${12.1}_{9.0}^{15.1}$	${8.4}_{8.0}^{8.8}$
SFG@z < 4	${3.3}_{2.6}^{3.7}$	${10.4}_{9.9}^{10.7}$	${1.4}_{0.7}^{1.9}$	${1.0}_{0.6}^{1.5}$	${0.9}_{0.5}^{1.4}$	${1.7}_{1.1}^{2.1}$	${0.9}_{0.3}^{1.6}$	${0.74}_{0.35}^{1.11}$	${11.7}_{8.8}^{14.8}$	${8.3}_{8.0}^{8.8}$
SFG@z > 4	${5.0}_{4.4}^{5.3}$	${9.9}_{9.5}^{10.3}$	${1.8}_{1.1}^{2.5}$	${0.6}_{-0.0}^{1.2}$	${0.4}_{-0.0}^{0.8}$	${1.8}_{1.3}^{2.3}$	${1.0}_{0.2}^{2.1}$	${0.48}_{0.01}^{0.77}$	${13.5}_{9.5}^{15.1}$	${8.8}_{8.4}^{9.1}$
QG	${3.6}_{3.3}^{4.0}$	${10.3}_{9.5}^{10.5}$	$-{0.6}_{-1.0}^{0.2}$	$-{0.6}_{-1.0}^{0.1}$	$-{0.5}_{-0.8}^{0.1}$	${0.4}_{-0.4}^{0.8}$	${0.2}_{0.0}^{0.5}$	${0.82}_{0.77}^{1.04}$	${11.7}_{9.2}^{13.6}$	${8.4}_{8.0}^{8.9}$
QG@z < 3.5	${3.2}_{2.6}^{3.3}$	${10.4}_{9.9}^{10.6}$	${0.0}_{-1.4}^{0.2}$	$-{0.1}_{-1.3}^{0.1}$	$-{0.0}_{-1.2}^{0.1}$	${0.7}_{-0.1}^{0.9}$	${0.2}_{0.1}^{0.6}$	${0.82}_{0.71}^{1.23}$	${9.7}_{8.6}^{13.9}$	${8.2}_{7.9}^{8.7}$
QG@z > 3.5	${3.9}_{3.7}^{4.4}$	${9.8}_{9.4}^{10.4}$	$-{0.8}_{-1.0}^{-0.0}$	$-{0.8}_{-1.0}^{-0.5}$	$-{0.7}_{-0.8}^{-0.4}$	${0.3}_{-0.4}^{0.7}$	${0.1}_{0.0}^{0.5}$	${0.82}_{0.77}^{0.82}$	${12.0}_{11.0}^{13.5}$	${8.6}_{8.4}^{9.0}$
XELG-z6	${6.5}_{6.3}^{6.6}$	${9.6}_{9.0}^{9.8}$	${2.2}_{1.9}^{2.8}$	${0.5}_{-1.0}^{1.3}$	$-{0.4}_{-1.0}^{0.8}$	${1.8}_{1.5}^{2.4}$	${1.4}_{0.8}^{2.4}$	${0.003}_{0.002}^{0.400}$	${14.4}_{8.6}^{17.8}$	${9.0}_{8.5}^{9.4}$

Note. Median and quartiles for relevant properties of the different galaxy subsamples presented in the paper, namely, all galaxies, HST-faint, HST-dark, SFGs (all, and divided into two redshift regimes), quiescent galaxies QGs (all, and divided into two redshift regimes), and extreme emission-line galaxies at z ∼ 6 XELG-z6. These properties are stellar mass; SED-based SFRs averaged in the last 5, 25, and 100 Myr; SFR obtained from the energy-budget calculations presented in Section 6.2.1; V-band attenuation; mass-weighted age; redshift at which the galaxies formed 5% of their current mass; and stellar-mass surface density (on a pixel-by-pixel basis).

Download table as:  ASCIITypeset image

Finally, we also show in this panel the SFH for XELG-z6, which starts at a very low level and increased very rapidly by a factor of 8 at z ∼ 6–7, similar to the behavior observed for z < 3.5 QGs and at a very similar epoch.

The right panel of Figure 9 shows the evolution of galaxies in an sSFR versus age of the universe plot. We remark that the tracks followed by galaxies roughly align with the line representing a mass-doubling time equal to the age of the universe (shown in black). This means that, on average, the SFH of galaxies transit through a steady star formation phase (the constant SFR region in the left panel) for several hundreds of Myr (even up to 1 Gyr). Indeed, the CANDELS sample (shown in yellow) nicely concentrates around that line. At z ≳ 10, however, our SFHs imply very quick mass-doubling times, which points to a very active phase in the young universe, consistent with the relatively high number of very high-z sources being detected in JWST surveys, several times above expectations (see Figure 5 in Finkelstein et al. 2022 and pink stars in our plot). We must warn the reader that the behavior of the SFHs at very high z can be affected by the very nonlinear relationship between time and redshift. By z ∼ 7, we start distinguishing the individual galaxies, which move to a starburst phase, a quiescent state, or remain in the mass-doubling line.

The right panel of Figure 9 also shows some interesting galaxy samples to help understand the nature of our sources. Apart from the general CANDELS sample shown in yellow, which follows the mass-doubling equal to the universe age line and the SFH tracks plotted in the figure, we also depict the galaxies detected with HST versus IRAC colors and presented in papers such as Wang et al. (2016) or Alcalde Pampliega et al. (2019), with the SFR and mass (quite uncertain due to data limitations prior to JWST) estimation they were able to assign. Most of these galaxies present high specific SFRs from redshifts z = 2–6, similar to some of our dusty starbursts. Another less numerous subpopulation presents specific SFRs compatible with quiescent systems, extending even at higher redshifts than our quiescent galaxy candidates.

Another interesting sample of galaxies we include in the right panel of Figure 9 are ALMA sources (see caption for references). Especially interesting are those at 4 < z < 6, which present similar specific SFR values to some of our z > 4 SFGs (some confirmed SMGs) and XELG-z6. The SMGs in our sample and the ALMA sources concentrate around the epoch where our QGs experienced their last intense starburst event before quenching.

Finally, we also include on the right panel of Figure 9 the high-redshift massive galaxies reported by Labbe et al. (2022) and the z > 9 galaxy candidates presented by several groups based on the first analyses of the first JWST data sets (Finkelstein et al. 2022; Naidu et al. 2022). In both cases, some of our SFHs would be able to reproduce their properties, but the bulk of our galaxies would be linked to progenitors presenting smaller specific SFRs than these high-z sources, implying also that we are just starting to probe the tip of the iceberg of z > 9 populations.

Figure 10 shows, on the left panel, the evolution of the stellar masses of our galaxies as a function of the universe's age, jointly with the same comparison samples depicted in the previous figure. Virtually all of our sources experience very quick mass assembly at z ≳ 10, followed by a less pronounced mass increment, very similar for all galaxies except for very few of them quickly increasing their mass by more than 1 dex below z = 7, the XELG-z6.

Zoom In Zoom Out Reset image size

The HST-faint and HST-dark galaxies selected in this paper correspond to a high mass end of the distribution of masses and redshifts of the CANDELS sample. Many of our galaxies have stellar masses similar to those of HST-dark systems selected with IRAC versus HST data (pink crosses), but those samples seem to be biased toward massive systems, $\mathrm{log}{M}_{\star }/{M}_{\odot }\gt 10$. This might be linked to limitations in the area of our CEERS survey, not comparable (yet) to the data sets used by those previous works. It might also be linked to confusion problems in the IRAC observations since many of our galaxies are indeed surrounded by neighbors and only the JWST data have now been able to separate them in the mid-IR. And, obviously, IRAC samples are biased toward brighter objects, compared to the new JWST data.

Compared to ALMA sources, our sample presents similar stellar masses at z > 4. Probably the selections are converging in this epoch. We also remark that some of the ALMA sources at 7 < z < 10 could well be progenitors of our lower-redshift galaxies because they lie within the locus occupied by our SFHs.

It is worth noting that the sample of XELG-z6 presents an average track, implying that they are quickly assembling a significant fraction of their stellar mass at z ∼ 7.

The last two samples we compare with are the Labbe et al. (2022) massive high-z galaxies and z > 9 JWST galaxy candidates (Naidu et al. 2022; Finkelstein et al. 2022). The former significantly deviates from our SFHs; they would seem to correspond more, should their stellar mass and redshift be accurate, to progenitors of very massive HST-dark galaxies, which we have not yet detected with JWST, maybe due to survey area limitations. In fact, for the galaxies in that paper that we have in our sample, we obtain similar redshifts but up to 10× smaller stellar masses. The masses and redshifts for the z > 9 samples compare well with the tracks followed by many of our SFHs, although they typically lie around 1 dex below the average SFH for our samples, implying that we are not seeing comparable populations. Furthermore, z > 9 detections would be more compatible with the $\mathrm{log}{M}_{\star }/{M}_{\odot }\lt 10$, 3 < z < 5 galaxies in the CANDELS or our sample (including quiescent objects).

The right panels of Figure 10 present the histograms of formation redshifts of our sample. On the bottom panel, we show when our red galaxies reached masses 10^8–10 M_⊙. Typically, 10⁸ M_⊙ of stars were assembled by z ∼ 15, with a large scatter. Our galaxies formed their first 10⁹ M_⊙ by z ∼ 12 (on average) and reached masses as high as 10¹⁰ M_⊙ at z = 7. Complementary, we also show the histograms of formation redshifts for 5%, 10%, and 25% of the current stellar mass of the galaxies in our sample. Overall, our SFHs, especially those for quiescent galaxies, imply a very active universe beyond z > 10 and up to z ∼ 25, with some quite massive galaxies ($\mathrm{log}{M}_{\star }/{M}_{\odot }\gt 10$) at even z = 15–20 (or, quite probably, several progenitors to merge later). We remark that all this star formation activity can occur in situ or ex situ in distinct progenitors at high redshifts.

In this section, we focus on the analysis of how the red galaxies in our sample have assembled, investigating the stellar-mass density 2D distribution, the ages of stellar populations present in different parts of the galaxy, and the location of the dusty starbursts that many of our sources present.

Figure 11 shows the distribution of the stellar population properties for all pixels in our analysis as a function of galactocentric distance (in absolute units and relative to the effective radius) and for different galaxy subsamples. Median values and quartiles are marked and provided in Table 2. We note three interesting results from this figure.

Zoom In Zoom Out Reset image size

First, focusing on the top-right panel, the z > 4 SFGs (including high-z SMGs) subsample (depicted in blue) is identified with the highest dust attenuations, typically A(V) > 2 mag. For z < 4 SFGs, there is a tail of systems harboring regions (pixels) with similarly large dust contents, many of which belong to the MIPS24 emitters in our sample. We remark that the submillimeter and millimeter surveys in the CEERS field are heterogeneous in area coverage and depth, so our z > 4 SFG sample could still include additional SMGs. Distinctively, quiescent galaxies are composed of regions with A(V) < 1 mag dust attenuations, the median being around 0.3 mag. The regions with the highest attenuations are preferentially located in the inner 2 kpc, or 2 × r_e,M , with a significant gradient toward larger radii, similar to what has been reported at lower redshifts (see Wang et al. 2017; Miller et al. 2022), but with larger nuclear attenuations in the case of SMGs (compared to main-sequence galaxies at z ∼ 2; see, e.g., Tacchella et al. 2018). XELG-z6 present a quite flat attenuation distribution, which extends from unobscured to A(V) = 3 mag, i.e., the starburst we see as an emission line is also dusty, but these are very compact sources, as we mentioned earlier.

The second result we highlight is that all the galaxies in our sample present evolved stellar populations in terms of mass-weighted age. The frequency plots for QGs and SFGs peak at similar ages, around 1 Gyr, although the distribution for quiescent galaxies is skewed to older values, the contrary can be said for SFGs. The distribution of mass-weighted ages for regions in z > 4 SFGs is different from that for z < 4 SFGs—the former tends to have more regions with ages between 1 and 10 Myr. Most of these very young regions are in the inner 2.5 kpc, within 1.5 × r_e,M . Even though the histograms of mass-weighted ages for QGs and SFGs (especially those at the lowest redshifts) peak at similar values, with the latter starting forming earlier, there is a significant number of pixels (dominating in number but not in mass) that had already formed 5% of their current mass at 5 < z < 10. This contrasts with the properties of XELG-z6—most of their pixels are assembling their mass quickly at z = 6–7.

The last result we extract from Figure 11 refers to the stellar-mass surface density. QGs and z < 4 SFGs present very similar distributions at the dense end, at values higher than 10^8.5 M_⊙ kpc⁻². Relatively, the regions with the highest densities, above 10^9.5 M_⊙ kpc⁻², are more common in quiescent galaxies and at distances within 2 kpc or one effective radius. We see a tail at the high-density end for z > 4 SFGs and XELG-z6, with larger values than those observed by QGs, but this is probably an artifact introduced by the different redshifts and spatial resolutions of the subsample. The two higher-redshift subsamples would be more comparable, and XELG-z6 present more pixels with higher densities, compared to z > 4 SFGs.