Unveiling the nature of infrared bright, optically dark galaxies with early JWST data

ABSTRACT

Over the last few years, both Atacama Large Millimeter/submillimeter Array (ALMA) and Spitzer observations have revealed a population of likely massive galaxies at z > 3 that was too faint to be detected inHubble Space Telescope(HST) rest-frame ultraviolet imaging. However, due to the very limited photometry for individual galaxies, the true nature of these so-called HST-dark galaxies has remained elusive. Here, we present the first sample of such galaxies observed with very deep, high-resolution NIRCam imaging from the Early Release Science programme CEERS. 30 HST-dark sources are selected based on their red colours across 1.6–4.4 |$\mu$|m. Their physical properties are derived from 12-band multiwavelength photometry, including ancillary HST imaging. We find that these galaxies are generally heavily dust-obscured (A_V ∼ 2 mag), massive (log (M/M_⊙) ∼ 10), star-forming sources at z ∼ 2−8 with an observed surface density of ∼0.8 arcmin⁻². This suggests that an important fraction of massive galaxies may have been missing from our cosmic census at z > 3 all the way into the Epoch of Reionization. The HST-dark sources lie on the main sequence of galaxies and add an obscured star formation rate density of |$\mathrm{3.2^{+1.8}_{-1.3} \times 10^{-3} \,{\rm M}_{\odot }\,yr^{-1}\,Mpc^{-3}}$| at z ∼ 7, showing likely presence of dust in the Epoch of Reionization. Our analysis shows the unique power of JWST to reveal this previously missing galaxy population and to provide a more complete census of galaxies at z = 2−8 based on rest-frame optical imaging.

galaxies: high-redshift, infrared: galaxies

1 INTRODUCTION

At z > 3, the rest-frame optical light redshifts out of the view of the Hubble Space Telescope (HST). Therefore, our understanding of ‘normal’ star-forming galaxies such as Lyman break galaxies (LBGs) that were identified in HST strongly relies on rest-frame ultraviolet (UV) observations. This can result in an incomplete galaxy census at earlier times due to UV-faint galaxy populations such as quiescent or dust-obscured sources.

While a small number of extremely obscured, star-forming galaxies [e.g. submillimetre galaxies (SMGs)] have been known to exist at z > 4 for a long time, such galaxies are 100× less common than LBGs (with a sky density of only 0.01 arcmin⁻²; e.g. Riechers et al. 2013; Marrone et al. 2018). Hence, they only provide a limited contribution to the integrated star formation rate density (SFRD). Similarly, massive (⁠|$\mathrm{\gt 10^{11}\,{\rm M}_{\odot }}$|⁠) quiescent objects have only been identified to z ∼ 3−4 (e.g. Tanaka et al. 2019; Carnall et al. 2020; Santini et al. 2020; Valentino et al. 2020), but they remain relatively rare such that degree-scale survey volumes are needed for their detection. Nevertheless, due to our limited sensitivity in rest-frame optical light at z > 3, it is still possible that the early cosmic star formation and stellar mass densities are dominated by less extreme, dusty star-forming galaxies (DSFGs) or quiescent galaxies, which have been missing from rest-frame UV data sets. The key question is: how common are such red sources at z > 3?

Over the last few years, evidence has been mounting that, indeed, a population of obscured, star-forming galaxies with less extreme star formation rates (SFRs) than SMGs could be quite common in the early Universe. These sources have been undetected in current HST surveys, due to their very red colours between |$1.6$| and |$4.5\, \mu$|m (e.g. Franco et al. 2018; Wang et al. 2019), and they have thus been named HIEROs (HST to IRAC extremely red objects), H dropouts, HST-dark galaxies, or HST-faint galaxies. Very red dusty sources exist even below the sensitivity limits in IRAC, making them only detectable by ALMA. However, previous studies with ALMA relied on small areas of the sky, thus reporting only one or two galaxies in some cases (e.g. Williams et al. 2019; Fudamoto et al. 2021). Hence, HST-dark sources are typically selected in Spitzer/IRAC imaging or as serendipitous sources in ALMA submillimetre continuum data (e.g. Huang et al. 2011; Simpson et al. 2014; Caputi et al. 2015; Stefanon et al. 2015; Wang et al. 2016, 2019; Franco et al. 2018; Alcalde Pampliega et al. 2019; Williams et al. 2019; Yamaguchi et al. 2019; Umehata et al. 2020; Fudamoto et al. 2021; Sun et al. 2021; Manning et al. 2022). However, due to the limited photometric information available, their spectral energy distributions (SEDs) and redshifts remained uncertain. Nevertheless, the general consensus in the literature is that the vast majority of these galaxies are heavily dust-obscured, star-forming sources at z ∼ 3−6 that dominate the high-mass end of the galaxy population at these early cosmic times (e.g. Alcalde Pampliega et al. 2019; Wang et al. 2019).

Accurately constraining the number densities and the physical properties of such sources is critical for our understanding of early galaxy build-up. It is conceivable that dust-obscured, HST-undetected sources dominate the SFRD at z > 4 unlike what is typically assumed (e.g. Williams et al. 2019; Gruppioni et al. 2020; Zavala et al. 2021). However, the most recent 2 mm number counts as well as IR luminosity functions measured from ALMA continuum observations out to z ∼ 4−7 indicate a less extreme scenario (see Casey et al. 2018; Barrufet et al. 2023). Nevertheless, it has become clear that our census of galaxies at z > 3 is indeed far from complete, leaving a critical gap in our understanding of early galaxy build-up.

Finding and characterizing this missing galaxy population is a key goal of extragalactic astronomy in order to inform theoretical models of galaxy evolution. In this paper, we address this important question with some of the first JWST NIRCam imaging from the Early Release Science (ERS) survey CEERS (Finkelstein et al. 2022). This allows us to identify HST-dark galaxies with the same selection that has been used in the past for HST + Spitzer samples, but now with significantly measured multiwavelength photometry in several bands to derive their redshifts and physical properties reliably for the first time. This follows previous analyses with JWST data that have already revealed massive galaxies at z > 7 (see Labbe et al. 2022; Santini et al. 2022). In a separate paper, we also analyse the morphologies of bright F444W-selected sources (Nelson et al. 2022).

This paper is organized as follows. In Section 2, we describe the observational data and sample selection. Section 3 outlines our SED modelling and the variety of physical properties that we find for HST-dark galaxies. In Section  4, we present our results and compare our galaxy sample to SMGs, before we derive their contribution to the cosmic star formation rate density (CSFRD) in Section  5 and finish with a summary and conclusions in Section  6.

Magnitudes are given in the AB system (e.g. Oke & Gunn 1983), and where necessary we adopt a Planck Collaboration XIII (2016) cosmology.

2 OBSERVATIONS AND SAMPLE SELECTION

2.1 CEERS observations and HST + JWST photometric catalogue

This analysis is based on the first four NIRCam pointings from the ERS programme CEERS (PI: Finkelstein; PID: 1345; see Finkelstein et al. 2022). These images include exposure sets of short- and long-wavelength filters: F115W +F277W, F115W+ F356W, F150W + F410M, and F200W + F444W. The typical integration time for each filter set is 0.8 h. Hence, the F115W images obtained double the exposure time compared to the remaining filters.

The data were retrieved from the MAST archive and were then WCS-matched and combined using the publicly available grizli¹ pipeline. Most notably, the images are aligned with the Gaia Data Release 3 catalogues. For details on the reduction, we refer the reader to Brammer et al. (in preparation). In the following, we work with images at 40 mas pixel scale. The 5σ depths are measured in circular apertures of 0.48 arcsec diameter. They range from 27.8 mag in the medium band filter F410M to 28.4–28.6 mag in the wide filters F200W, F277W, and F356W at |$2{\!-\!}3.5\, \mu$|m (see also table 1 in Naidu et al. 2022).

We also include the available ancillary data from HST since the CEERS survey covers the previous CANDELS/EGS field. In particular, we include a re-reduction of all the ACS and WFC3/IR imaging taken in these fields in the five filters: F606W, F814W, F125W, F140W, and F160W. The 5σ depths measured in the 0.48 arcsec circular apertures for the HST data range from 28.4 mag in F606W to 27.0 mag in F140W. While the HST near infra-red (NIR) data are shallower, they still provide additional multiwavelength constraints. Most importantly, using the F160W filter allows us to use HST-dark galaxy selections analogous to previous HST + Spitzer searches (see the next section).

We match the point spread functions (PSFs) in all but the WFC3 images to the F444W filter using the deconvolution algorithm from Lucy (1974). The multiwavelength photometric catalogues were derived with SExtractor (Bertin & Arnouts 1996), which we run in dual mode, using a PSF-matched, inverse-variance weighted stack of the F277W, F356W, and F444W images as the detection image and measuring fluxes in all bands in circular apertures of 0.48 arcsec diameter. We scale the obtained fluxes to total based on the AUTO flux measurements in the F444W image with the default Kron parameters. In order to take the PSF differences between F444W and the WFC3 filters into account, we match their PSFs to the F160W filter, which has the largest FWHM among all filters considered here. We run SExtractor again on the resulting images, using the same detection image as before. For each WFC3 filter, we then correct the fluxes measured in the original image to match the PSF-matched colour between the WFC3 filter and F444W, respectively. As a general quality cut, the catalogue contains objects with S/N > 5 in at least one of the available NIRCam-wide filters. We also removed objects near the edge and bad pixels; hence, we guarantee that all the sources of the catalogue are not spurious. Finally, to account for systematic uncertainties in the photometry we apply an error floor of 5  per cent to all flux measurements.

The final catalogue includes ∼26 000 sources with photometry in 12 bands, spanning 0.6–4.4 |$\mu$|m.

2.2 HST-dark galaxy colour selection with JWST

Previous HST-dark galaxy studies were based on analyses of HST + Spitzer/IRAC data (e.g. Wang et al. 2016; Alcalde Pampliega et al. 2019; Wang et al. 2019). These selections were limited by two factors: (1) the very broad PSF of the IRAC data with FWHM of ∼1.7 arcsec and (2) the limited depth of available IRAC imaging. Hence, typical HST-dark samples were limited to very bright sources with [4.5] ≲ 24.5 mag. Now, with the unparalleled high-resolution, deep JWST data, we can push such analyses to much fainter sources. Nevertheless, in this first paper, we focus on sources that are identified using very similar selection methods as previous HST + Spitzer samples in order to reveal their nature.

We use a similar colour cut to the ‘traditional’ HST-dark galaxies selection colour adopted in the literature (e.g. Caputi et al. 2012; Wang et al. 2016). This colour cut is specifically designed to identify dusty galaxies at z ∼ 3–6. Indeed, previous analyses have tentatively confirmed that such red HST-dark sources lie at these redshifts on average. However, the available photometry for these sources is very poor, with significant detections only in the IRAC 3.6 and 4.5 |$\mu$|m bands, and mostly upper limits at HST wavelengths.

In this first paper, we substitute the IRAC data for NIRCam data. To find an appropriate colour cut, we compute theoretical fluxes using Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (BAGPIPES; Carnall et al. 2018). We perform similar modelling that we use on the SED fitting (see Section 3). We set to a delayed SFH (e-folding time τ = 3 Gyr), a stellar mass of |$\mathrm{{\it M}_{*}=10\ {\rm M}_{\odot }}$| and a metallicity of Z = 0.5 Z_⊙. We model the fluxes F160W and F444W with a dust attenuation set to A_V = 2.0, 3.0,  4.0 mag, respectively, tracking the redshift evolution with steeps Δz = 0.1 between 2 < z < 6. We find that applying a colour cut of F160W–F444W > 2.0 selects galaxies at z > 3 at A_V = 2.0 (see Fig. 1).

Figure 1.

Colour–magnitude diagram (F160W–F444W versus F444W) used to select the HST-dark galaxies. The dashed line shows the location of H = 27 mag, which corresponds to the average 3σ limit of the shallowest parts of the H band in the HST data. The black line marks the adopted colour cut H-F444W > 2.0 and the red points correspond to the galaxies included in our sample. Of the 39 HST-dark galaxies selected, 19 have lower limit fluxes in the H band (grey arrows). The continuous lines show the theoretical evolution in galaxy colours at redshifts of 2 < z < 6 for three different levels of reddening A_v = 2.0, 3.0, and 4.0 mag (yellow, orange, and purple lines, respectively). The grey dots represent the density of sources from the total CEERS catalogue with detections in the F444W band (∼26 000 sources).

We thus apply a colour cut of F160W–F444W > 2.0 to our multiwavelength catalogue to select HST-dark galaxies. After also limiting the sample to F160W > 27 mag, we identify a total of 39 galaxies in the ∼40 arcmin² down to an F444W = 26.4 mag (red diamonds in Fig. 1). We remove six galaxies that are at the edge of the image and do not have a clear detection in the longer wavelength of the wide NIRCam filters (F444W, F356W, F277W). And three more galaxies for which the photometry is extremely noisy, likely due to remaining detector artefacts. Hence, our final sample consists of 30 HST-dark galaxies.

3 SED FITTING AND VARIETY OF HST-DARK GALAXIES

We derive the physical properties of the 30 HST-dark galaxies by using BAGPIPES (Carnall et al. 2018). This code generates a complex model galaxy spectra to fit these to arbitrary combinations of spectroscopic and photometric data using the MultiNest nested sampling algorithm (Feroz, Hobson & Bridges 2009).

The assumptions by Wang et al. (2016, 2019) have proven effective to characterize HST-dark galaxies and this work follows the same SED modelling. Specifically, we use a delayed star formation history (SFH) with an e-folding time from τ = 0.1–9 Gyr and a uniform prior. Hence, we include relatively short bursts as well as, effectively, a constant SFH.

The stellar models are based on the 2016 updated version of the Bruzual & Charlot (2003) library with a Kroupa initial mass function. We allow for a broad range of metallicities from 0.2 to 2.5 Z_⊙. Nebular continuum and line emission are added self-consistently based on the photoionization code cloudy(Ferland et al. 1998) with the ionization parameter set as logU = −2. Finally, we adopt Calzetti et al. (2000) dust attenuation allowing for very heavily attenuated SEDs by setting a range |$\mathrm{{\it A}_{v} = 0{\!-\!}6\, mag}$|⁠, again using a uniform prior.

The above settings adequately reproduce the range of SEDs of our sample. We find good fits for all our sources. Fig. 2 shows postage stamps and SED fits of four of these sources at z > 3.3. They represent the variety of galaxy SEDs that are identified as HST-dark sources, both in terms of physical properties and in terms of morphology. We also show the cut-outs of previous IRAC images from the S-CANDELS programme (Ashby et al. 2015). The outstanding increase in depth and spatial resolution at |$3{\!-\!}5\, \mu$|m can be appreciated. With the exception of a handful of our sample galaxies, none of these would have been easily selected based on previously available data. Hence, the majority of our sample was completely missing from our previous cosmic census of HST + Spitzer. In the following, we use our sample of 30 sources to analyse the physical properties of HST-dark galaxies. The main parameters are summarized in Table  A1.

Figure 2.

Examples of four HST-dark galaxies at z ∼ 3.3, 3.6, 4.6, 8 that represent the variety of physical properties, isolation/clustering and morphologies that we find in our sample (see also Nelson et al. 2022). Top panels: Postage stamps in the optical/NIR for four examples of HST-dark galaxies, where the sources show no significant detections in the HST filters. The first row shows HST images with the bands F606W, F814W, F850LP, F125W, F140W, F160W, and the Spitzer bands IRAC1, IRAC2 at 3.6  and  4.5 |$\mu$|m, respectively. The second row shows new NIRCam JWST images: F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The extraordinary improvement in both depth and resolution of JWST at |$3-5\, \mu$|m compared to the IRAC images is obvious. Most of these sources would not have been easily selected in previous data and they would have potentially been missing from our cosmic census. Bottom panels: Posterior SEDs (yellow; median with 16–84 percentile uncertainties) with the photometry (red dots). Orange points indicate the expected model fluxes from the posterior SEDs. The probability distribution function of the photometric redshift is shown in the lower right part of the panel. In general, HST-dark sources are found to be dust-obscured, massive, star-forming galaxies at z ∼ 2−8.

4 HST-DARK GALAXIES: THE HIGH-Z EXTENSION OF SMGs?

4.1 HST-dark galaxy properties

Thanks to the very sensitive NIRCam data, we can now – for the first time – measure accurate photometric redshifts and SEDs for HST-dark galaxies. Hence, we can now derive significantly more reliable physical parameters such as SFRs and stellar masses compared to previous IRAC-selected HST-dark sources. In this section, we thus present the physical properties derived from SED fitting for our final sample of 30 sources (see Figs  2 and 3).

Figure 3.

Histograms of the main physical properties analysed: photometric redshift, SFR (first row), stellar masses, and dust attenuation (second row). The histograms show the generally dusty high-z nature of the sample with most of the galaxies at z > 3.3 and a median attenuation of A_v ∼ 2 mag. The SFRs are moderate, mostly at |$\mathrm{SFR \lt \, 50 \,{\rm M}_{\odot }\,yr^{-1} }$|⁠, whereas the stellar masses are large with |$\mathrm{log({\it M}_{*}/{\rm M}_{\odot }) \sim 9.5{\!-\!}10}$|⁠. Our analysis thus shows that previous stellar mass function estimates at z > 3 were likely underestimated at high masses.

We find that an HST-dark galaxy selection includes a variety of sources with the vast majority of the sample being high redshift dust-obscured, massive, star-forming systems. The photometric redshifts are accurate for the galaxies at z > 3.3 with small uncertainties (Δz ∼ 0.1−0.2) with only one solution in the redshift posterior distribution function (see Table A1 and Fig. 2). The subset of galaxies that is placed at z < 3.3 presents more considerable uncertainties and typically shows two possible redshift solutions.

As expected, we confirm the dusty nature of these colour-selected HST-dark galaxies. We have allowed the dust attenuation to vary broadly in the modelling, but the values only reach up to A_v = 4 mag for almost all HST-dark galaxies with a median value for the sample of A_v ∼ 2.1 mag. An alternative procedure to classify among dusty and quiescent galaxies is the commonly accepted UVJ diagrams (see e.g. Spitler et al. 2014). Fig. 4 shows the rest-frame UV colours for the total CEERS sample, corroborating the dusty nature of HST-dark galaxies.

Figure 4.

UV against VJ rest-frame colours diagram. HST-dark galaxies are presented as big circles with the H-F444W colour on the third axis, whereas the grey dots are the total CEERS sample. The black-dashed lines show the usual classification of galaxies in quiescent, star-forming and dusty from Spitler et al. (2014). Our sample of HST-dark galaxies lies mostly in the dusty colour space with part of the sample in the star-forming galaxy region. There are no quiescent galaxies in our selection due to the colour cut; however, given the photometric uncertainties, we cannot exclude that some of the galaxies are quiescent.

Interestingly, the sample reaches redshifts z ∼ 8 for several galaxies. In total, we identify eight galaxies at z > 6, i.e. in the epoch of reionization, whereas previous HST-dark galaxies selected with Spitzer were typically limited to z < 6. The high attenuation for the galaxies above z > 6, with A_V = 1.7 ± 0.4 suggests the presence of dust in the reionization epoch (see also Fudamoto et al. 2021). We refer to a future paper for more detailed analyses of these sources.

We have also analysed both the stellar masses and the SFRs from the SED fitting. The SFRs are moderate, when compared to SMGs, with |$\mathrm{\sim 80~{{\ \rm per\ cent}}}$| of our sample presenting |$\mathrm{SFR \,\lt\, 50 {\rm M}_{\odot }\,yr^{-1}}$|⁠. The stellar masses of our sample are |$\mathrm{log({\it M}_{*}/{\rm M}_{\odot })\,\gt \, 8.8 }$| and for the eight galaxies at z > 6 we find |$\mathrm{log({\it M}_{*}/{\rm M}_{\odot })= 9.9 \pm 0.3}$|⁠. Hence, we can confirm that HST-dark galaxies are massive with relatively low star formation.

4.2 HST-dark galaxies compared to SMGs

Given their dusty nature, it is interesting to compare the HST-dark galaxies with a classical SMG sample in order to put them in a broader context of massive high-z dusty galaxies. We find that HST-dark galaxies have higher redshifts compared to SMGs (see the histogram in Fig. 5). Interestingly and probably not unexpectedly, HST-dark galaxies have lower SFRs than SMGs. Similarly, we also find HST-dark galaxies have more than an order of magnitude lower stellar masses than the massive SMGs that weigh |$\mathrm{{\it M}_{*} \sim 10^{12} {\rm M}_{\odot }}$|⁠.

Figure 5.

HST-dark galaxies compared with SMGs. Top panel: The redshift histogram shows that HST-dark galaxies (red line) lie mostly at z > 3.2 and extend to higher redshift compared to SMGs (orange bars). Bottom panel: Photometric redshifts against SFRs. On average, HST-dark galaxies lie at higher redshift and have significantly lower SFRs [orange squares; ALESS sample from da Cunha et al. (2015)]. The significantly smaller error bars in the photometric redshifts are due to the improvement coming from the deep, high-resolution JWST/NIRCam data.

In Fig. 6, we further show the location of HST-dark and SMGs on the SFR–M_* diagram with respect to the main sequence (MS) of galaxies at z ∼ 4. As can be seen, HST-dark sources at these redshifts follow the MS. They populate its lower mass end compared to both SMGs and IRAC-selected HST-dark galaxies. Specifically, we have compared our results with two previous studies of dusty galaxies. The ALESS sample from da Cunha et al. (2015) shows a larger number of starburst galaxies that lie significantly above the MS. In between our sample and SMGs, we find more massive HST-dark galaxies from Wang et al. (2019). This difference could be due to the selection criteria and lower mass HST-dark galaxies could have been missing from Spitzer-selected samples.

Figure 6.

SFR against stellar mass for our sample corresponding to 3 < z < 5 (red diamonds) and the total sample (empty diamonds). The black line shows the MS at z = 4 from Schreiber et al. (2015). HST-dark galaxies lie on the MS as other optically dark massive galaxies selected with Spitzer data from Wang et al. (2019) (yellow dots). SMGs from da Cunha et al. (2015) are above the MS (orange squares).

In a nutshell, HST-dark galaxies appear to be the higher redshift, lower SFR extension of SMGs. They have likely just not been detected in previous submillimetre single-dish studies due to their faintness. The synergy between JWST and ALMA data will be needed to confirm this rather interesting outcome.

5 COSMIC SFR: IS HST MISSING A SIGNIFICANT CONTRIBUTION AT HIGH Z?

The cosmic SFH is key to explaining the mass assembly and the metal enrichment of the galaxy population both vital to understanding galaxy evolution. Although the CSFRD is well constrained up to z ∼ 8 for UV-selected galaxies (i.e. LBGs), the contribution from obscured, IR rest-frame selected galaxies is still highly uncertain and only poorly constrained at z > 3 (e.g. Casey, Narayanan & Cooray 2014; Casey et al. 2019; Zavala et al. 2021). A central question is: when does the dust-obscured star formation start to dominate the CSFRD? And, more specifically, what is the contribution to the CSFRD from HST-dark galaxies? In this section, we approach this question and put our results in the context of the most relevant high- z galaxy populations.

We calculate the SFRD for our sample of HST-dark galaxies at z > 3.5 in four redshift equidistant bins centred at z = 3.3, 4.5, 5.7, and 6.9, respectively. We find a contribution of |$\mathrm{SFRD = 1.4^{+1.1}_{-0.8} \times 10^{-3}\, {\rm M}_{\odot }\,yr^{-1}\,Mpc^{-3}}$| at z = 4.5 that remains constant within the uncertainty up to z ∼ 7. Specifically, we measure |$\mathrm{SFRD \sim 2.0^{+2.0}_{-1.1} \times 10^{-3}\, {\rm M}_{\odot }\,yr^{-1}\,Mpc^{-3}}$| at z ∼ 6, and |$\mathrm{SFRD \sim 3.2^{+1.8}_{-1.3} \times 10^{-3}\, {\rm M}_{\odot }\,yr^{-1}\,Mpc^{-3}}$| at z ∼ 7.

Fig. 7 shows that our results have similar SFRD to previously selected HST-dark galaxy samples at z ∼ 4−6. There are no current observations of HST-dark galaxies at z > 6. Our results, however, suggest that the SFRD of HST-dark galaxies remains constant, within the uncertainty, or becomes even slightly larger at z ∼ 7. This could suggest a remarkable presence of dust in the Epoch of Reionization. As a caveat, we notice that our selection criteria tend to select HST-dark galaxies at z ∼ 7 due to the bright F444W flux that selects [O iii] lines in the SED. Also, our redshift bins currently only contain a few sources each and future studies with larger samples are needed to corroborate this result.

Figure 7.

SFRD against redshift for the HST-dark galaxies (red diamonds) compared to a compilation of previous studies. The orange area represents the dust-obscured SFH from semiempirical models presented in Zavala et al. (2021), whereas the blue area shows the unobscured star formation (Bouwens et al. 2022). The contribution of our HST-dark galaxies remains almost constant from redshift 3.5 < z < 6. The HST-dark contribution is similar to SMGs at z ∼ 5−6 (orange squares; da Cunha et al. 2015). Our SFRD is in agreement with previous HST-dark studies at z ∼ 4−6 (orange dots; Wang et al. 2019), but we extend the result to z ∼ 7.

Our analysis shows an HST-dark contribution to the SFRD similar to SMGs at 4 < z < 5, but larger at z > 5. An evident explanation is that the number of SMGs drops at z > 5 (da Cunha et al. 2015). Even if their SFRs are an order of magnitude higher than the average of HST-dark galaxies, their SFRD is lower at high z.

Our study shows that HST-dark galaxies contribute an order of magnitude less than 2 mm selected DSFGs at z ∼ 4 from the MORA survey presented in Casey et al. (2021). The reason probably is the lower SFRs and this suggests that fainter dusty galaxies could have been missed by shallow ALMA surveys at z ∼ 4–7.

Our study and previous studies mentioned above indicate that the dusty population at high z is larger than what was expected. In particular, more extensive samples are needed from deeper surveys and multiwavelength observations to understand how the disparate objects seen today relate and contribute to the overall cosmic SFRD. A very important result from our analysis is that HST-dark galaxies significantly contribute to the CSFRD, in particular at z > 5.

6 SUMMARY AND CONCLUSIONS

In this work, we have exploited some of the first JWST/NIRCam data to identify red high- z galaxies that were missing from rest-frame UV selections based on HST. These HST-dark galaxies were selected based on H-F444W colours. We studied 30 galaxies in the four ERS CEERS NIRCam pointings and used SED fitting to derive their properties and put them into context by comparison to other dusty high-redshift galaxy populations.

The main results are summarized as follows:

• The average properties of HST-dark galaxies show their high-z dusty nature with |$\mathrm{80~{{\ \rm per\ cent}}}$| of the sample at 3 ≲ z ≲ 8. Surprisingly, we identify eight galaxies at z > 6 with a median dust attenuation of A_v ∼ 1.9 ± 0.4 mag. This indicates that massive, heavily obscured galaxies have been missing from our cosmic census in the Epoch of Reionization.

• We find that HST-dark galaxies have moderate SFRs compared to SMGs with |$\mathrm{ SFR \,\lt \,50\, {\rm M}_{\odot }\,yr^{-1} }$|⁠, but they extend to higher redshifts. HST-dark galaxies lie on the MS of galaxies with lower masses than SMGs, but still with a median of |$\mathrm{log({\it M}_{*}/{\rm M}_{\odot })= 9.9 \pm 0.3}$| for the eight galaxies at z > 6.

• We find that the contribution of HST-dark galaxies to the CSFRD is |$\mathrm{SFRD \sim 3.2^{+1.8}_{-1.3} \times 10^{-3} \,{\rm M}_{\odot }\,yr^{-1}\,Mpc^{-3}}$| at z ∼ 7 at z ∼ 7 remaining almost constant over redshift 3.3 < z < 7. Our SFRD results are in broad agreement with previous HST-dark galaxy studies, but the JWST-selected sample indicates that HST-dark galaxies existed up to higher redshifts than previously assumed, i.e. z ∼ 8.

Thanks to the outstanding improvement of NIRCam data over the previous Spitzer/IRAC imaging, this work provides more reliable photometric redshifts, SFRs, and stellar masses compared to previous HST-dark galaxy studies. Hence, we can finally conclude that HST-dark galaxies have a high-z nature and that they are numerous at z > 3, with a surface density of ∼0.8 arcmin⁻². The high number of dusty galaxies with z > 6 suggests that the presence of dust is relevant in the Epoch of Reionization.

The SFRD contribution of our sample of HST-dark galaxies is in agreement with previous H-dropout studies at z < 6. We also find a significant SFRD at z ∼ 7 that needs to be corroborated with larger samples of HST-dark galaxies. But already now, with these first NIRCam images, JWST has demonstrated its power to identify more complete samples of rest-frame optically selected galaxies out to very high redshifts. We have shown that dusty galaxies at z > 3 are likely to be far more numerous than expected and that such HST-dark galaxies provide an important, but hitherto underestimated contribution to the high-mass end of the galaxy mass function and the cosmic SFRD up to z ∼ 7−8.

ACKNOWLEDGEMENTS

We acknowledge support from the Swiss National Science Foundation through project grant 200020_207349 (LB, PAO, and AW). We also acknowledge the anonymous referee for the constructive feedback that resulted in an improved version of this manuscript. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant no. 140. YF acknowledges support from NAOJ ALMA Scientific Research grant number 2020-16B. VG acknowledges support by the ANID BASAL projects ACE210002 and FB210003. KEH acknowledges support from the Carlsberg Foundation Reintegration Fellowship Grant CF21-0103. RS acknowledges an STFC Ernest Rutherford Fellowship (ST/S004831/1). MS acknowledges support from the CIDEGENT/2021/059 grant, and from project PID2019-109592GB-I00/AEI/10.13039/501100011033 from the Spanish Ministerio de Ciencia e Innovacion–Agencia Estatal de Investigacion. RJB and MS acknowledge support from NWO grant TOP1.16.057.

Cloud-based data processing and file storage for this work are provided by the AWS Cloud Credits for Research programme.

This work is based on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with programme #1345.

Facilities: JWST, HST.

Software: matplotlib (Hunter 2007), numpy (Oliphant 2015), scipy (Virtanen et al. 2020), jupyter (Kluyver et al. 2016), astropy (Astropy Collaboration 2013, 2018), grizli (Brammer 2018), SExtractor(Bertin & Arnouts 1996), and bagpipes (Carnall et al. 2018).

DATA AVAILABILITY

The calibrated (stage 2) data used here are available at Mikulski Archive for Space Telescopes (https://mast.stsci.edu). Further data products are available from the authors upon request.

Footnotes

1

https://github.com/gbrammer/grizli/

APPENDIX

See Table A1.

Table A1.

HST-dark galaxy sample and derived properties.

ID	RA	Dec.	z	Log SFR	Log Mass	AV
	(deg)	(deg)		(M⊙ yr−1)	(M⊙)	(mag)
9795	214.9984037	53.0046235	8.06 |$^{+0.13}_{-0.10}$|	2.30 |$^{+0.10}_{-0.09}$|	10.24 |$^{+0.11}_{-0.09}$|	2.69 |$^{+0.18}_{-0.17}$|
18038	214.8147216	52.8036165	8.05 |$^{+0.17}_{-0.11}$|	1.99 |$^{+0.09}_{-0.08}$|	9.94 |$^{+0.10}_{-0.09}$|	1.77 |$^{+0.13}_{-0.15}$|
2921	215.0084889	52.9779736	8.01 |$^{+0.11}_{-0.08}$|	1.61 |$^{+0.08}_{-0.07}$|	9.54 |$^{+0.08}_{-0.08}$|	1.63 |$^{+0.14}_{-0.12}$|
7194	214.9272387	52.9338962	7.94 |$^{+0.08}_{-0.06}$|	1.90 |$^{+0.06}_{-0.06}$|	9.88 |$^{+0.10}_{-0.08}$|	1.64 |$^{+0.09}_{-0.10}$|
4836	214.9233733	52.9255929	7.28 |$^{+0.05}_{-0.07}$|	1.86 |$^{+0.06}_{-0.07}$|	9.79 |$^{+0.08}_{-0.07}$|	2.00 |$^{+0.09}_{-0.10}$|
32479	214.7944062	52.8373983	7.07 |$^{+0.07}_{-0.04}$|	1.68 |$^{+0.08}_{-0.06}$|	9.71 |$^{+0.15}_{-0.13}$|	1.46 |$^{+0.11}_{-0.07}$|
17951	214.8867950	52.8553796	6.83 |$^{+0.07}_{-1.33}$|	2.22 |$^{+0.15}_{-0.13}$|	10.46 |$^{+0.14}_{-0.23}$|	2.37 |$^{+0.23}_{-0.18}$|
15222	214.8493871	52.8118241	6.55 |$^{+0.20}_{-0.33}$|	1.60 |$^{+0.10}_{-0.13}$|	10.03 |$^{+0.09}_{-0.09}$|	1.49 |$^{+0.13}_{-0.09}$|
29279	214.8383981	52.8851900	5.99 |$^{+0.35}_{-0.17}$|	1.54 |$^{+0.08}_{-0.07}$|	9.47 |$^{+0.09}_{-0.08}$|	1.74 |$^{+0.06}_{-0.07}$|
23985	214.8644947	52.8709797	5.36 |$^{+0.07}_{-0.05}$|	2.06 |$^{+0.09}_{-0.07}$|	10.01 |$^{+0.10}_{-0.09}$|	2.13 |$^{+0.14}_{-0.09}$|
23008	214.8588182	52.8620152	5.18 |$^{+0.02}_{-0.01}$|	1.44 |$^{+0.20}_{-0.16}$|	9.36 |$^{+0.20}_{-0.15}$|	3.19 |$^{+0.57}_{-0.42}$|
11011	214.9319470	52.9589227	5.17 |$^{+0.05}_{-0.05}$|	0.91 |$^{+0.05}_{-0.06}$|	8.86 |$^{+0.07}_{-0.07}$|	1.14 |$^{+0.10}_{-0.08}$|
25718	214.8400370	52.8606498	4.72 |$^{+0.93}_{-0.13}$|	1.33 |$^{+0.11}_{-0.19}$|	9.60 |$^{+0.12}_{-0.19}$|	1.66 |$^{+0.34}_{-0.25}$|
29538	214.7665813	52.8315228	4.59 |$^{+0.16}_{-0.19}$|	1.10 |$^{+0.11}_{-0.10}$|	9.58 |$^{+0.06}_{-0.09}$|	1.74 |$^{+0.18}_{-0.15}$|
11555	214.9155599	52.9484970	4.47 |$^{+0.13}_{-0.14}$|	0.90 |$^{+0.08}_{-0.10}$|	9.59 |$^{+0.05}_{-0.05}$|	0.77 |$^{+0.12}_{-0.13}$|
22858	214.8475516	52.8533710	4.36 |$^{+1.68}_{-0.52}$|	1.60 |$^{+0.42}_{-0.36}$|	10.03 |$^{+0.34}_{-0.35}$|	3.95 |$^{+0.50}_{-0.47}$|
24502	214.8081702	52.8322122	4.22 |$^{+0.12}_{-0.13}$|	1.54 |$^{+0.10}_{-0.10}$|	10.10 |$^{+0.05}_{-0.07}$|	2.07 |$^{+0.12}_{-0.10}$|
30673	214.8482898	52.8847859	3.84 |$^{+0.24}_{-2.09}$|	1.06 |$^{+0.28}_{-0.97}$|	9.50 |$^{+0.20}_{-0.45}$|	2.01 |$^{+0.61}_{-0.28}$|
9594	214.9508404	52.9668643	3.74 |$^{+0.46}_{-1.99}$|	0.91 |$^{+0.21}_{-0.95}$|	9.51 |$^{+0.17}_{-0.70}$|	2.28 |$^{+0.79}_{-0.28}$|
9784	214.9091114	52.9372131	3.65 |$^{+1.95}_{-1.90}$|	1.39 |$^{+0.64}_{-0.96}$|	9.95 |$^{+0.46}_{-0.63}$|	3.26 |$^{+1.69}_{-0.41}$|
27303	214.7680265	52.8164017	3.56 |$^{+0.13}_{-0.15}$|	2.14 |$^{+0.14}_{-0.14}$|	10.68 |$^{+0.09}_{-0.12}$|	3.34 |$^{+0.18}_{-0.20}$|
16342	214.8871236	52.8453781	3.56 |$^{+0.69}_{-2.02}$|	0.72 |$^{+0.27}_{-0.97}$|	9.30 |$^{+0.29}_{-0.68}$|	2.06 |$^{+0.73}_{-0.39}$|
20994	214.8556793	52.8487133	3.50 |$^{+0.20}_{-1.74}$|	0.33 |$^{+0.17}_{-0.65}$|	8.96 |$^{+0.09}_{-0.38}$|	1.48 |$^{+0.43}_{-0.19}$|
24835	214.8894314	52.8922071	3.30 |$^{+0.10}_{-0.15}$|	0.43 |$^{+0.10}_{-0.15}$|	8.92 |$^{+0.10}_{-0.11}$|	0.87 |$^{+0.17}_{-0.14}$|
2439	214.9440430	52.9297453	2.04 |$^{+1.79}_{-0.95}$|	0.38 |$^{+0.73}_{-0.83}$|	9.20 |$^{+0.32}_{-0.53}$|	2.68 |$^{+1.10}_{-0.44}$|
22807	214.8476086	52.8534052	1.86 |$^{+2.01}_{-0.23}$|	0.32 |$^{+0.95}_{-0.20}$|	9.15 |$^{+0.62}_{-0.23}$|	4.53 |$^{+0.61}_{-1.18}$|
21771	214.8560277	52.8546738	1.76 |$^{+1.85}_{-0.11}$|	0.81 |$^{+1.05}_{-0.22}$|	9.67 |$^{+0.35}_{-0.17}$|	3.63 |$^{+0.45}_{-0.61}$|
4710	215.0045567	52.9835262	1.73 |$^{+1.86}_{-0.20}$|	1.02 |$^{+0.98}_{-0.25}$|	9.99 |$^{+0.62}_{-0.19}$|	4.40 |$^{+0.49}_{-1.16}$|
28170	214.7598238	52.8334124	1.68 |$^{+1.90}_{-0.12}$|	-0.09 |$^{+0.88}_{-0.22}$|	8.80 |$^{+0.50}_{-0.22}$|	3.33 |$^{+0.45}_{-0.69}$|
10881	214.9257706	52.9544456	1.32 |$^{+0.72}_{-0.22}$|	0.14 |$^{+0.54}_{-0.31}$|	9.14 |$^{+0.40}_{-0.23}$|	4.04 |$^{+0.77}_{-0.80}$|