Delving deep: a population of extremely dusty dwarfs observed by JWST

TL;DR

This study uses JWST NIRCam data to identify extremely dusty, low-mass galaxies at redshifts below 2, revealing unexpectedly high dust extinction levels in these faint objects.

Contribution

It introduces a new selection method for dusty galaxies using JWST data, uncovering low-mass, highly obscured galaxies at z<2 with unexpectedly high dust extinction.

Findings

Most dusty galaxies are at z<2 with median stellar mass of 10^7.3 M_sun.

No strong dusty galaxy candidates found at z>6.5.

Galaxies exhibit median dust extinction of Av=4.9 mag.

Abstract

We take advantage of the NIRCam photometric observations available as part of the Cosmic Evolution Early Release Science survey (CEERS) to identify and analyse very red sources in an effort to discover very dusty star forming galaxies. We select red galaxies as objects with a S/N>3 at 4.4 $\mu$m and a S/N<2 in all JWST and HST filters at $\lambda\leq2\mu$m, which corresponds to [F200W]-[F444W]>1.2 considering CEERS depths. This selection is ideal to identify very dusty (Av>1 mag) galaxies with stellar masses between $10^6$ to $10^{10}\, \rm M_{\odot}$ at z<5, more massive dusty galaxies at z=5-18 and galaxies at z>18 due to the Lyman absorption, independently of their dust extinction. Our sample of F200W-dropouts contains no strong candidates at z>6.5, instead it consists almost completely (~81%) of z<2 low-mass galaxies, with a median stellar mass of $10^{7.3} \rm M_{\odot}$. These galaxies show an exceptional dust extinction with median value of Av=4.9 mag, completely unexpected given their low stellar mass. The remaining galaxies, which are at z<6.5, show similar large dust extinction (Av>1), but they are generally more massive $>10^{7.5}\rm M_{\odot}$.

Tables

Table 1. Table 1: Summary of the filters and 5 σ 𝜎 \sigma depths of the catalogue considered in this work.

JWST	depth(5$\sigma$)	HST	depth(5$\sigma$)
F115W	29.2	F606W	28.6
F150W	29.0	F814W	28.3
F200W	29.2	F105W	27.1
F277W	29.2	F125W	27.3
F356W	29.2	F140W	26.7
F410M	28.4	F160W	27.4
F444W	28.6

Full text

*Methods. *We select red galaxies as objects with a $S/N>3$ at 4.4 $\mu$m and a $S/N<2$ in all JWST and HST filters at $\lambda\leq 2\mu$m, which corresponds to $\rm[F200W]-[F444W]>1.2$ considering CEERS depths. This selection is ideal to identify very dusty ($A_{V}>1$ mag) galaxies with stellar masses between $10^{6}$ to $10^{10}\,\rm M_{\odot}$ at $z<5$, more massive dusty galaxies at $z=5-18$ and galaxies at $z>18$ due to the Lyman absorption, independently of their dust extinction.

*Results. *Our sample of F200W-dropouts contains no strong candidates at $z>6.5$, instead it consists almost completely ($\sim 81\%$) of $z<2$ low-mass galaxies, with a median stellar mass of $10^{7.3}\rm M_{\odot}$. These galaxies show an exceptional dust extinction with median value of $A_{V}=$4.9 mag, completely unexpected given their low stellar mass. The remaining galaxies, which are at z¡6.5, show similar large dust extinction ($A_{V}>1$), but they are generally more massive $>10^{7.5}\rm M_{\odot}$.

1 Introduction

Dust is an ubiquitous component of galaxies. It has a significant impact on the observed spectra as well as on the derived physical properties and, in turn, on our understanding of galaxy evolution. Indeed, around half of the cumulative energy emitted by stars at all epochs is absorbed by the surrounding dust and re-emitted at infrared (IR) wavelengths (Hauser & Dwek, 2001). Therefore, the dust content and its evolution with redshift remains one of the main topic in extragalactic astrophysics.

For example, dust attenuation is one of the major uncertainties when deriving the cosmic star-formation-rate density (SFRD), particularly, but not only, at $z>3$ (Madau & Dickinson, 2014). On the one hand, studies based on rest-frame ultra-violet (UV) data (e.g. Bouwens et al., 2015; Oesch et al., 2018), which are largely affected by dust but have been possible up to $z=10$, have reported a sharp decline of the SFRD at $z>2$. On the other hand, far-IR or radio observations (Rowan-Robinson et al., 2016; Novak et al., 2017; Gruppioni et al., 2020; Enia et al., 2022), which have been generally limited to $z<7$ due to observational challenges, have shown that the SFRD may show a less drastic decline, or even a plateau, at $z>2$. At the same time, the comparison between IR and rest-frame UV and optical observations have highlighted the presence of sources that were completely missed by even the deepest HST observations (e.g. Simpson et al., 2014; Franco et al., 2018; Wang et al., 2019; Talia et al., 2021; Gruppioni et al., 2020). The presence of such sources show how the absence of deep IR observations may have biased our past view of the stellar mass assembly of the Universe. Moreover, dusty galaxies are interlopers when searching for $z>10$ galaxies, as their red continuum could be misinterpreted as a Lyman-break at ultra-high redshift (Wilkins et al., 2014; Holwerda et al., 2015, 2020; Zavala et al., 2022).

The dust attenuation present in a star-forming galaxy has a positive correlation with the stellar mass (e.g. Pannella et al., 2015; McLure et al., 2018), with more massive galaxies being in general more dust attenuated and galaxies below 10${}^{8.5}\,\rm M_{\odot}$ having a little dust attenuation. This is not surprising given that a large stellar mass generally correspond to intense star-formation activity (Brinchmann et al., 2004; Noeske et al., 2007) and, therefore, to large dust content, as the dust is produced by stellar activity such as pulsating moderate mass stars and supernovae (e.g. Scalo & Slavsky, 1980). However, the search for dusty low-mass galaxies requires very deep IR observations, which are not possible with HST or Spitzer, as they are extremely faint in the rest-frame optical given the combination of low stellar mass and high dust extinction. In addition, the mentioned relation between stellar mass and dust extinction seems not to evolve with redshift up to $z\sim 3$ (e.g., Whitaker et al., 2017; Shapley et al., 2022), with a possible decrease at higher redshifts (e.g. Fudamoto et al., 2020).

The launch of JWST has opened, as expected, a new observational window in the Universe. On the one hand it has allowed for identifying a large number of $z>10$ candidates just in the first few months of observations (e.g. Adams et al., 2022; Carnall et al., 2022; Castellano et al., 2022; Donnan et al., 2022; Finkelstein et al., 2022b; Naidu et al., 2022; Pérez-González et al., 2023), some of which have shown a puzzling, high content of dust (Rodighiero et al., 2023). On the other hand, it has enabled statistical investigation of the nature of optically-faint and optically-dark galaxies (Barrufet et al., 2022; Pérez-González et al., 2022), which was impossible to do with previous near-IR observatory such as Spitzer. In particular, such analysis can be pushed now to intermediate and small stellar masses, without being limited to the most massive and IR-luminous galaxies.

In this work we take advantage of some of the first photometric observations obtained by JWST to study some of the reddest galaxies observable. Such red colours could be caused by a high redshift or large dust content. In particular, we consider the first observations of the Cosmic Evolution Early Release Science (CEERS111\hrefhttps://ceers.github.io/https://ceers.github.io/), which cover an area of $\sim$36 arcmin2 that is larger than the area of the early release observations on SMACS0723 ($\sim 9\,\rm arcmin^{2}$) or the GLASS-JWST early release science (18 arcmin2; Treu et al., 2022), thus allowing the selection of rare sources and increasing the number statistics of previous surveys. The structure of the paper is the following. In Section 2 we describe the used data and the adopted selection criterion. In Section 3 we describe the method used to derive the physical property associated to each object and we analyse them in Section 4. We then summarise our work and report our conclusive remarks in Section 5. Throughout the paper, we assume a Chabrier initial mass function (Chabrier, 2003) and a $\Lambda$CDM cosmology with $H_{0}=70\,{\rm km}\,{\rm s}^{-1}{\rm Mpc}^{-1}$, $\Omega_{\rm m}=0.3$ and $\Omega_{\Lambda}=0.7$. All magnitudes are in the AB system (Oke & Gunn, 1983).

2 Data and sample selection

2.1 CEERS survey and catalogues

CEERS is a JWST early release science (ERS) program (Proposal ID $\#$1345, Finkelstein et al. in prep.). CEERS covers approximately 100 arcmin2 in the Extended Groth Strip field (EGS) with both imaging and spectroscopic observations. In particular, it consists of 10 pointings with the Near Infrared Camera (NIRCam) combined with 6 pointings with the Near InfraRed Spectrograph instrument (NIRSpec) and 8 pointings with the Mid InfraRed Instrument (MIRI).

In the same field there are also archival HST observations from the All-wavelength Extended Groth Strip International Survey (Davis et al., 2007), the Cosmic Assembly Deep Extragalactic Legacy Survey (CANDELS; Grogin et al., 2011; Koekemoer et al., 2011) and 3D-HST (Momcheva et al., 2016). These catalogues have been updated (v1.9222\hrefhttps://ceers.github.io/releases.htmlhttps://ceers.github.io/releases.html) aligning their astrometry with that of Gaia DR3.

For this work we consider the CEERS NIRCam catalogue presented in Finkelstein et al. (2023). In particular, the catalogue is based on the NIRCam images of the first CEERS public data release (Data Release 0.5; Bagley et al., 2022), taken as part of CEERS on the 21 June 2022. These images, for which no complete overlapping MIRI observations are available, consists of four pointings (1, 2, 3 and 6) in three short wavelength filters (i.e. F115W, F150W, and F200W) and four long-wavelength bands (i.e. F277W, F356W, F410M, and F444W). Source detections were based on the inverse-variance-weighted sum of the F277W and F356W images, after matching them by their point-spread-function (PSF) and removing the background. PSFs were derived stacking observed stars present in the field.

Photometric measurements were performed after matching all images to the F444W PSF, which has the largest full-width at half-maximum (FWHM) among the considered JWST filters. In particular, we made use of the PyPHER routine (Boucaud et al., 2016) to derive the kernel necessary to match each PSFs to the F444W one and then we convolved each image with its respective kernel. This process was applied also to HST/ACS F606W and F814W imaging, as their PSF FWHMs are smaller than the F444W one. We instead did not convolve the HST/WFC3 band as their PSF FWHMs are larger than F444W one, but we applied an additional aperture correction to take this into account.

More details on the catalogue creation and image reductions are presented in Finkelstein et al. (2023) and Bagley et al. (2022), respectively, while we list in Table 1 the filters included in the catalogue and the respective observational depths. The latter were derived, as explained in details in Finkelstein et al. (2022a, 2023), by extracting from each image the fluxes inside non-overlapping circular apertures of different radius, from 0.1” to 3”, all placed avoiding real sources. These fluxes were then used to derived how the noise varies as a function of the number of pixels in each aperture. The final 5$\sigma$ depths were estimated using the mentioned noise function in an aperture of diameter 0.2”, then corrected to total flux using the PSF.

2.2 Sample selection

We select galaxies that have a $\rm S/N>3$ in the F444W filter and a $\rm S/N<2$ in every filter equal or below 2 $\mu$m (i.e. HST/F606W, HST/F814W, HST/F105W, HST/F125W, HST/F140W, HST/F160W, JWST/F115W, JWST/F150W, and JWST/F200W)333Using a magnitude limit equivalent to a S/N$=2$ in all filters at $\lambda\leq 2\mu$m, instead of a S/N cut, results in the selection of a different galaxy sample, because of local noise fluctuation. However, we verify with pointing 3 that the general properties of the samples selected using a S/N cut and a magnitude cut are consistent.. Considering the CEERS observational depths this corresponds roughly to a colour cut [F200W]-[F444W]¿1.2. We checked that results are consistent if we consider a $\rm S/N>5$ in the F444W filter, even if the resulting sample is reduced in number. We highlight here that this method does not select equally red galaxies at all magnitudes, as objects that are brighter in the F444W filter must be intrinsically redder in the [F200W]-[F444W] colour to meet our S/N limits. The selection method has been however chosen to include $z>10$ galaxies, without limiting the sample only to blue, relatively dust-free, objects.

We also limit our analysis to galaxies with a magnitude [F444W]$<29$ mag, which roughly corresponds to the 3$\sigma$ limit of the four pointings (Bagley et al., 2022). For this selection we consider the photometry derived with the Kron aperture, with a Kron factor and a minimum radius set to 1.1 and 1.6, respectively. These values were chosen via simulations to maximise the signal-to-noise value. We also correct all fluxes in order to retrieve the total flux, using larger apertures (see Finkelstein et al., 2023, for further details). We then visually checked each selected source to remove artefacts, for example, due to the image edges or to the presence of saturated stars checking also each co-added F277W, F356W, and F444W band images. Finally, we co-added the F115W, F150W and F200W band images and exclude every source with a S/N$\geq$2 in it. Since the original catalogue has been selected in the F277W+F356W co-added image, we checked for additional sources detected in the F444W filter and not present in catalogue by Finkelstein et al. (2023), using the same extraction methodology and selection criteria.

The final sample of F200W-dropouts based on the catalogue by Finkelstein et al. (2023) consists of 135 objects, of which 15, 13, 48, and 59 are present in pointing 1, 2, 3 and 6, respectively. Among these sources, 128 have $S/N>3$ in at least two filters, indicating that they are not spurious sources. In addition, after removing artefacts (e.g., edges, cosmic rays, bad pixels) and sources with contaminated photometry, we found nine sources (three in pointing 2 and 3, two in pointing 1 and one in pointing 6) not present in the catalogue based on the F277W+F356W co-added image as they are faint in one or both of these bands. However, four of these sources are detected only in the F444W filter, so we do not attempt to derive any physical property for them. We call this sub-sample F200W-dropouts-extra hereafter. The total sample of F200W-dropout, including the F200W-dropouts-extra sub-sample, is therefore made of 144 objects. In the next sections we will analyse all 144 objects when looking at their colours, but we will consider only the 133 of them detected in at least two filters when performing the SED fitting.

Differences in numbers in the four pointings may be due to the presence of over-densities, as visible in Figure 1. However, these possible over-densities do not correspond to known galaxy groups from X-ray observations Erfanianfar et al. (2013) or density analysis (Coogan et al., 2023; Jin et al., 2022), therefore additional spectroscopic observations are necessary to confirm them.

Some differences are visible in the number of sources detected in the two NIRCam detectors, i.e. 32 in one detector and 112 in the other. First, as verified in Bagley et al. (2022), the photometric residuals from detector-to-detector are at the level of $\sim 2-5\%$ in all filters and sources do not overlap with any feature visible in the flat fields of the long filters. Second, if we apply a higher magnitude cut, i.e. [F444W]¡27.8 corresponding to a 10$\sigma$ level, the difference is still present, as four galaxies are present in one detector and ten in the other one. We can therefore exclude that the observed difference is due to different noise level and calibration.

2.3 What do we expect from this selection?

Before moving on analysing the observed sample, to investigate our selection criterion we derived the expected JWST fluxes for all the spectral energy distribution (SED) templates used later to derive physical properties, but considering a coarser grid in the different physical parameters. We then use these fluxes to understand which model is compatible with the F200W-dropout selection described in the previous paragraph. In particular, we included Bruzual & Charlot (2003) stellar population models with stellar metallicity from 0.005 solar up to solar (i.e., 0.005 Z*⊙, 0.02 Z⊙, 0.2 Z⊙, 0.5 Z⊙, Z⊙). We allow the redshift to vary up to $z=20$ ($z=$0.1, 0.3, 0.5 and from 1 to 20 with steps of 1) and the stellar mass to range from 10${}^{6}\,\rm M_{\odot}$ up to 10${}^{12.5}\,\rm M_{\odot}$, with steps of $\rm log(M_{*}/M_{\odot})=0.5$. Nebular continuum and emission lines are included based on a pre-computed model grid, already included in the SED fitting code (see Sec. 3.1), based on the photoionisation code CLOUDY (Ferland et al., 2017) assuming a ionisation parameter (i.e. dimensionless ratio of densities of ionising photons to hydrogen) of $\log_{10}U\in{-4,-3,-2}$, a Hydrogen density of 100 atoms cm-3* (see Carnall et al., 2018, for more details). We considered the same reddening law (i.e., Calzetti et al., 2000) for both the stellar continuum and the nebular emission lines with dust extinction up to $A(V)=6$, in steps of 0.5 mag. Finally, we considered a delayed exponential (i.e. SFR$\propto t\,e^{-(t/\tau)}$) star-formation history with ages ranging from 1 Myr to the age of the Universe (0.001, 0.01, 0.1, 1, 3, 5, 7, and 15 Gyr) and $\tau=$0.01 to 10 Gyr (i.e., 0.01, 0.1, 1, 3, 5, and 10 Gyr). To identify SED templates not detected at $\lambda\leq 2\,\mu$m, we consider 5$\sigma$ limits of 29.2, 29.0 and 29.2 mag for the F115W, F150W and F200W filters (Table 1).

The physical properties of the templates corresponding to the F200W-dropout selection are reported in Figure 2. First, the considered selection is expected to identify every galaxy candidate at $z\geq 18.2$ (considering the wavelength range where the F200W throughput is at least 50$\%$ of the maximum value) massive enough to be observed at 4.4$\mu$m, as there is no flux at $\lambda\leq 2\mu$m due to the inter-galactic medium absorption (IGM). Second, at $13.4\leq z<18.2$ the flux in the F200W is only partially absorbed from intergalactic hydrogen, so we expect to observe galaxies with stellar masses above 109 M*⊙* and various dust-extinction. At $z<13.4$ only some combinations of stellar mass and dust extinction can make a galaxy visible at 4.4$\mu$m, but not at $\lambda\leq 2\mu$m. For example, at $z<2$, this happens only to galaxies with $\rm M_{*}<10^{8.5}\,M_{\odot}$ and $A_{V}>3$, while galaxies with with $\rm M_{*}>10^{11}\,M_{\odot}$ are possible at $z>3.7$, as otherwise they would be detected at $\lambda\leq 2\mu$m. Third, we may select passive galaxies with $\rm log(sSFR/(M_{\odot}/yr))\leq-11$, only up to $z=6$. Overall, the chosen selection criteria is biased toward very high-$z$ galaxies and toward dusty sources.

3 Physical property derivation

3.1 SED fitting

We use the Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (BAGPIPES; Carnall et al., 2018) to derive both the redshift and galaxy physical properties (e.g. stellar mass, SFR) using the models anticipated in Section 2.3, but with a finer grid of values for each physical parameter. We include all filters listed in Table 1 and for fluxes with low S/N values (i.e., $S/N<3$), we do not impose any hard upper limit in the SED fitting, but we include the measured photometry with the associated uncertainties. This SED fitting procedure was performed using the photometry derived with the Kron aperture, but in Appendix A we repeat the analysis using a 0.2”-radius aperture. This change in the input photometry sometimes results in different physical properties for the individual galaxies, but the general results for the overall F200W-dropout sample remain unchanged. In Figure 3 we report, as example, the SED fitting procedure applied to one galaxy in the F200W-dropout sample. In the entire sample, 132/133 objects results in a good SED fit with a $\chi^{2}<10$. For comparison, the same SED fitting procedure described for the F200W-dropout was applied also to the entire sample in the CEERS.

To verify the robustness of our results we repeated the BAGPIPES fitting changing some key parameters. The details of these additional results are reported in Appendix B. First, we assumed the reddening law observed in the Small Magellanic Cloud, instead of the reddening law by Calzetti et al. (2000), but this change has a minor impact on the results. Second, in our SED fitting we assumed a Draine & Li (2007) dust emission models, but Leja et al. (2017) pointed out that this model sometimes over-predicts the 3.3 $\mu$m polycyclic aromatic hydrocarbon (PAH) emission line strength. This spectral feature is inside the F444W filter indicatively between $z=0.1$ and $z=0.5$ and we therefore decided to repeat the BAGPIPES fit limiting the mass associated to PAH to only 10% of the total dust mass. Even in this case the impact on our results is minor. Finally, we have changed the star-formation history to an exponentially declining one (i.e. SFR$\propto e^{-(t/\tau)}$), which result on differences in the physical properties of each galaxy, while leaving the general results of the paper unchanged.

3.2 Are we delving too shallow?

In this section we explore the possibility that some of the objects in our F200W-dropout sample are not extra-galactic sources but, instead, galactic brown dwarfs. To perform this test we fit our candidates F200W-dropouts with L and T dwarf models from Burrows et al. (2006). In particular, these templates span metallicities between [Fe/H]=-0.5 and 0.5, effective temperatures between 700 K and 2300 K and gravities between $10^{4.5}$ and $10^{5.5}\,\rm cm\,s^{-2}$. Of the full F200W-dropout sample, only two objects have a $\chi^{2}$ derived with the brown dwarf templates that is better than the $\chi^{2}$ derived with the galaxy models. However, the normalisation of the best brown dwarf templates for these objects indicates that they should be all within 10 pc from the Sun.

In the solar neighbourhood, McKee et al. (2015) found a surface density of brown dwarfs of $\sim 1.2\,\rm M_{\odot}/pc^{2}$, which corresponds to a number density of 16-40 $pc^{-2}$ if we associate to them a stellar mass between 0.03 and 0.075 $\,M_{\odot}$. This implies that we would expect between 3 and 7 brown dwarfs in the area covered by four NIRCam pointings, which is consistent with the found numbers. However, additional near-IR spectroscopic data are necessary to verify the presence of the absorption features characteristic of brown-dwarf and confirm their nature. These two possible stellar objects are further discussed in Appendix C.

4 Results

In this Section we analyse the observed colours and the physical properties of our F200W-dropouts sample, comparing them with other literature samples and with the full CEERS sample.

4.1 Colours

As a first test, in Figure 4 we report the [F444W] magnitude and the [F356W]-[F444W] colour of our F200W-dropout sample, divided into the four NIRCam pointings. We also show the full CEERS galaxy sample for comparison. We see no major differences between the F200W-dropouts in the four pointings, excluding any strong bias in the photometry. At the same time, F200W-dropouts have generally redder [F356W]-[F444W] colours with respect to the full CEERS sample, i.e. a median value of 0.55 (0.56, including objects in the F200W-dropout-extra sub-sample) against -0.15, respectively. Moreover, only 4/144 objects in the F200W-dropout sample have blue $\rm[F356W]-[F444W]<0.0$ colours. Analysing the SED models presented in Section 2.3, a blue [F356W]-[F444W] colour in our F200W-dropout sample is generally associated to $z>16$ templates with almost no dust extinction. However, once the photometric errors in both filters are taken into account all four objects are consistent with red colours within 1$\sigma$. All the nine objects in the F200W-dropouts-extra sub-sample have $\rm[F356W]-[F444W]>1$ showing that they are the red tail of our F200W-dropout sample.

Looking at the [F277W]-[F444W] colours (Figure 5), our F200W-dropout sample shows again generally redder colours than the full CEERS sample, with median values of 0.83 (0.85 including the F200W-dropout-extra sub-sample) against -0.19 mag, respectively. However, these values are affected by large uncertainties, given than 46$\%$ of our F200W-dropout sample has a S/N$<$3 in the F277W filter. Again, the nine objects in the F200W-dropouts-extra sub-sample are the red tail of our F200W-dropout sample, having $\rm[F277W]-[F444W]>0.6$. We do not show any additional colour, as the F410M observations are shallower than the ones in other filters and indeed only 28$\%$ of the F200W-dropout sample has a S/N$>$3 in this filter. At the same time, the remaining filters at $\lambda<2\,\mu$m are by construction associated to S/N$<$2.

Some galaxies in the general CEERS sample are very red in the [F277W]-[F444W] and [F356W]-[F444W] colours. These objects are not included in the F200W-dropout sample because their are bright at $\lambda\leq 2$, i.e. $S/N>2$, or have $S/N<3$ in the F444W filter.

In both Figure 4 and 5 we can also see that our F200W-dropouts occupy the same color-magnitude space of the bulk of the F200W-dropouts presented in Rodighiero et al. (2023). Only some of the sources in our F200W-dropout-extra sample has [F356W]-[F444W] and [F277W]-[F444W] colours similar or redder than one of their two candidates at $z>10$ (i.e. KABERLABA) and only a few of our F200W-dropout sources have colours and magnitudes similar to the other $z>10$ candidate (i.e. PENNAR). However, the only galaxies with [F356W]-[F444W] colour similar to KABERLABA is detected only in the F444W filter, making any further comparison impossible. Among the remaining four objects with [F356W]-[F444W] colours redder than KABERLABA, only the brightest one (N2-202826) in the F444W filter is detected in more than one filter (i.e. F444W and F410M), with a marginal detection in the F356W filter.

We also compare our F200W-dropout sample with the F200W-dropout samples by Yan et al. (2023a, b). Their sample selection is slightly different from the one adopted in our work, as they applied a colour cut [F200W]-[F356W]$\geq$0.8 mag, a S/N$<$2 only in the F150W band and a S/N$>$5 in the F277W filter. This selection results in galaxies with both [F277W]-[F444W] and [F356W]-[F444W] colours bluer than our sample.

4.2 Physical properties: a population of dusty dwarfs at $z<2$

We now move on analysing the physical properties derived for our F200W-dropout sample, taking into account, however, that these are results based on few observational points, given the large number of non-detections, and additional data are necessary to give more robust constrains. In Figure 6 we report the redshift, stellar mass and dust extinction for our F200W-dropout sample. In agreement with the test performed with the SED templates, our F200W-dropout sample is composed at 81$\%$ by galaxies at $z<2$, with a median stellar mass of 10${}^{7.3}\,\rm M_{\odot}$ and a median dust extinction of 4.9 mag. The entire sample of F200W-dropouts has A${}_{V}>$0.9 showing an extremely large dust content.

The pointing N3 and N6 have a larger number of F200W-dropouts and they show an excess in the redshift distribution at $\rm z\sim 1.5$. This secondary redshift peak is possibly due to over-densities present in these pointings, as previously shown in Figure 1, but additional spectroscopic data will be necessary to confirm it.

In our F200W-dropout sample we see a lack of galaxies between $z=4$ and $z=5.5$ (see left panel in Fig. 6), but this is neither motivated by the selection applied, as SED templates at these redshifts and with the colours of our F200W-dropout galaxies are indeed possible, neither by the presence of strong nebular emission lines, as the Hα nebular emission line is inside one of the red filters (i.e. F277W, F356W, and F444W) from $z\sim 2.5$ to $z\sim 6.5$ and should in case boost the flux of these filters. This is for example the case for two of the sources in the F200W-dropout-extra sub-sample, as the excess flux observed in these galaxies in the F444W filter could be explained by the Hα nebular line at $z=5.5-6.5$.

We note that the galaxy properties of our F200W-dropout sample is extremely different respect to the optical-dark populations presented in previous works using the same CEERS data (i.e. Pérez-González et al., 2023; Barrufet et al., 2022). Those optical-dark galaxies are selected based on their [F150W]-[F444W] or [F150W]-[F356W] colours and correspond to a population of dust-obscured galaxies at $z>2$ with intermediate to high stellar masses ($M>10^{9}M_{\odot}$). While all galaxies in the F200W-dropout sample are faint at 1.5$\mu$m with [F150W]$>$28.5 mag, only $\leavevmode\nobreak\ 40\%$ have [F150W]-[F444W]$>2$ mag or [F150W]-[F356W]$>1.5$ mag, as used by Barrufet et al. (2022) and Pérez-González et al. (2023), respectively. These colours criteria are however merely indicative for our sample, given that galaxies in our F200W-dropout sample are completely undetected (i.e. S/N$<2$) in the F150W filter. At the same time, galaxies in the F200W-dropout sample are also quite faint in the F444W filter with magnitudes from 29.0 to 27.5, while the sample by Barrufet et al. (2022) has [F444W]$<$26.4 mag and the one by Pérez-González et al. (2023) has [F444W]4$<$27.5 mag.

In Figure 7 we show the distribution of our F200W-dropout sample in the SFR vs. stellar mass plane. The majority of our galaxies are consistent with the main-sequence of star-forming galaxies (MS,e.g., Brinchmann et al., 2004; Noeske et al., 2007; Whitaker et al., 2012; Steinhardt et al., 2014; Santini et al., 2017; Bisigello et al., 2018), as extrapolated from the parametrisation by Speagle et al. (2014).

If we concentrate on galaxies at $z<2$, we obtain that their median SFR is 0.01 M${}_{\odot}/yr$. This value, if converted using the relation by Kennicutt (1998), results in a total IR luminosity (8-1000 $\mu$m) of L${}_{IR}\sim 9\times 10^{7}\,\rm L_{\odot}$, generally below Herschel capability beyond the local Universe. This shows why such galaxy population has been missed up to now: they are in general too red for HST and to faint for Spitzer and Herschel.

Looking into more details in the SFR-M*∗* plane (Fig. 7), the fraction of starburst, roughly defined as galaxies above 0.6 dex (Rodighiero et al., 2011) from the MS, in our F200W-dropout sample is of 26$\%$ at $z\leq 1$, 27$\%$ at $1<z\leq 2$ with only few other cases around $z\sim 3.5$. These results suggest that a burst of star-formation is probably not the cause of the large dust extinction. It is however necessary to take these results with caution given that the estimated SFR is based on the SED analysis of rest-frame UV-to-optical or optical-to-near-IR data, depending on the redshift. In the absence of far-IR observations these SFR estimates may be underestimated (e.g. Puglisi et al., 2017; Elbaz et al., 2018; Xiao et al., 2023), particularly for the most dusty objects. At the same time, it is necessary to take into account that intrinsic correlations are present between the SFR and the stellar mass, given that both are derived from the SED fitting. Moreover, the assumed delayed SFH and the minimum allowed age of 1 Myr limit the sSFR to values below $10^{-7.9}yr^{-1}$, as observed also in other works (e.g. Ciesla et al., 2017). Finally, the galaxy physical properties are based on a limited number of photometric points with large S/N values.

Of the full F200W-dropout, only five objects (4% of the sample) have a SFR at least four times below the MS, showing therefore an exceptional low SFR. However, all these objects are well described by the SED templates of star-forming galaxies when considering the aperture photometry, indicating that a second not-passive solution is also plausible.

As previously shown in Figure 6 (right panel), our F200W-dropout sample show in general a very high dust extinction. This is even more striking if we compare it with galaxies of similar stellar mass (Fig. 8). Indeed when considering all redshifts, while our F200W-dropouts have median A${}_{V}=4.7$, all the galaxies in CEERS have a median value of A${}_{V}=0.48$. The striking difference is shown in Figure 8, where we also report the median AV for the entire CEERS sample in different stellar mass bins. This median values are in agreement with the average relation between dust extinction and stellar mass (e.g., McLure et al., 2018), showing that the large dust extinction of our F200W-dropout is not a direct artefact of the adopted SED fitting procedure, but is dictated by their red colours. Overall, the dusty dwarf population we identify represent 0.5$\%$ of the population of galaxies with similar stellar mass (i.e. $10^{5.8}-10^{7.2}M_{\odot}$) at $z<1$ and 0.8$\%$ of galaxies with $10^{6.8-8.2}M_{\odot}$ at $z=1-2$, assuming that the completeness of the entire CEERS sample and of the F200W-dropout at those stellar masses are the same. This show that these exceptional galaxies clearly are rare and not representative of the general dwarf population at $z<2$.

In the same Figure it is possible to see that there are galaxies with $A_{V}>1$ that are not present in our F200W-dropout sample. On the one hand, the absence of galaxies with stellar masses above $\sim 10^{8}\,\rm M_{\odot}$ at $z<2$ in our F200W-sample is purely due to a selection effect, as mentioned in Section 2.3. As an additional check, we search for dusty galaxies with $>10^{8}\,\rm M_{\odot}$ in the entire NIRCam pointing 2 and we found 167 objects with A${}_{V}>$1 and 29 with A${}_{V}>$3, the majority of which have $S/N>2$ in the F200W filter.

On the other hand, there are different reason beyond the non selection of some galaxies with stellar masses below $10^{8}\,\rm M_{\odot}$ and $A_{V}>1$. Around 44$\%$ of these objects are simply too faint in the F444W filter, as they have [F444W]$>29$ mag or a $S/N<3$. The remaining objects have a S/N$>$2 in at least one of the filter at $\lambda>2\mu$m. In general, all these not-selected low-mass dusty sources show in general bluer colour than our F200W-dropout sample, with median [F277W]-[F444W]=-0.21 and median [F356W]-[F444W]=-0.16. These results on lower stellar masses at a fixed dust extinction.

4.3 On the absence of $z>10$ galaxies

As visible in Figure 6, no galaxies are present at $z>6.5$ ($z>9$ considering aperture photometry, see Appendix A) in our F200W-dropout sample. On the one hand, some blue $z>10$ galaxies, as the candidates presented in Finkelstein et al. (2023), are discarded from our sample as they have $S/N>3$ in the F200W filter, given that the flux at these wavelength is not totally absorbed by the IGM. On the other hand, observations considered in this work cover an area four times larger than the observations in SMACS0723 analysed in Rodighiero et al. (2023) and would therefore expect to observe at least few dusty galaxies at $z>10$. However, it is necessary to take into account that one of their two $z>10$ candidates, i.e. KABERLABA, is extremely red and faint, e.g. [F356W]=30.5 mag, and it is therefore below the detection limit of CEERS in all bands except F444W. On the contrary, a galaxy similar to the other $z>10$ candidate, i.e. PENNAR, would be bright enough to be detected in at least three filters, even without the magnification of the SMACS0723 cluster.

Even if no galaxy are associated to $z>10$, 93/133 objects have a probability larger than 10$\%$ of being at this high-z, as visible from the average redshift probability distribution in Figure 9. By integrating the probability distribution at $z>10$ of each of these objects and summing up all of them, we expect our sample to statistically include 17 galaxies at $z>10$, which however would require a spectroscopic follow-up to be identified and confirmed. Such secondary solutions correspond to a median stellar mass of 10${}^{9.5}\,\rm M_{\odot}$, median $\rm log(SFR)=1.4$ and median dust extinction of $\rm A_{V}=2$ mag, which are in line with the physical properties of the $z>10$ candidate PENNAR from Rodighiero et al. (2023).

5 Summary and conclusions

In this work we have taken advantage of the first four NIRCam pointings available as part of the CEERS survey. Available observations include six NIRCam broad-band filters, from 1.15 $\mu$m to 4.4$\mu$m, the medium F410M NIRCam filter and five HST filters, covering from 0.6$\mu$m to 1.6$\mu$m. Such photometric observations were considered to explore a population of extremely red galaxies, defined as objects with $S/N>3$ at 4$\mu$m and $S/N<2$ in all filters at $\lambda\leq 2\mu$m and in the co-added image of the JWST blue filters (F115W, F150W and F200W). We verify, using available SED models, that such colour selection is particularly suited for identifying dusty (A${}_{V}>$1) galaxies at $z<18.2$, but in principle also $z\geq 18.2$ galaxies, independently by their dust extinction.

This search results in a sample of 144 F200-dropouts, with eight of them detected only in the F444w filter. These objects have in general both [F277W]-[F444W] and [F356W]-[F444W] colours redder than the general galaxy population and we found no galaxies with secure blue colours, excluding the presence of very high-z galaxies. Indeed, when performing an SED-fitting analysis we found no candidate galaxy at z$>$9, showing the rarity of dusty high-z sources, which have been instead identified in other fields (Rodighiero et al., 2023). We however speculated, using the redshift probability distribution associated to each galaxy, that our sample could statistically include 17 galaxies with $z>10$ misinterpreted as low-z one.

At the same time, $\sim 81\%$ of the F200W-dropout sample corresponds to extremely dusty dwarf galaxies (median A${}_{V}=$4.9 mag and median stellar mass of 10${}^{7.3}\,\rm M_{\odot}$) at $z<2$. These are mainly star-forming galaxies, both inside and above the MS, with median SFR$=0.01\,\rm M_{\odot}/yr$. Such extreme dust content is striking, given that the dust content of galaxies, which correlates with the stellar mass, is expected to be negligible for galaxies with stellar masses below 10${}^{8.5}\,\rm M_{\odot}$. However, these population of extremely dusty dwarfs have a minor contribution on the overall galaxy population, being less than 1$\%$ of the galaxy sample at similar stellar masses and redshifts, but they could contribute to the diffuse extragalactic background light (Driver et al., 2016; Koushan et al., 2021; Kramer et al., 2022).

We can speculate about different possible scenario causing this excess in dust extinction. First, these galaxies could be more compact than galaxies with similar size, causing a large dust extinction even if the overall dust content remain similar. However, galaxies with such stellar masses are expected to be too small to be resolved even with JWST, extrapolating the stellar-mass size relation (van der Wel et al., 2014), and we can therefore not test this scenario. Second, these galaxies could been experiencing a burst of star-formation, having therefore an increase in dust production. Our preliminary analysis shows that only a fraction of these galaxies are considerably above the main-sequence of star-formation, but such estimates are heavily limited by the low-number of photometric data available and the absence of far-IR observations. Third, they could be the remnants of more massive galaxies, but again this scenario can not be tested with the available data. Finally, these galaxies could host an active galactic nuclei, but we can not again verify this possibility with the few photometric point available. Therefore, we need to wait for future spectroscopic observations or data at longer wavelengths to investigate further these galaxy population and confirm their exceptional dust content.

Acknowledgements.

LB and GR acknowledge the support from grant PRIN MIUR $2017-20173ML3WW$˙$001$. GG acknowledge M. Giulietti for the helpful discussions. LC acknowledges financial support from Comunidad de Madrid under Atracción de Talento grant 2018-T2/TIC-11612 and Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through grant PGC2018-093499-B-I00

Appendix A Comparison with aperture photometry

As reported before, the sample selection and the derivation of the physical properties are based on the Kron aperture photometry, however, given the faintness of our sources, the photometric extraction could strongly influence our results. In this Appendix we analyse the reliability of our results repeating the analysis by deriving the flux on a 0.2”-radius aperture, corrected to total flux.

In Figure 10 we compare the photometry derived with the two apertures. The ratios between fixed and Kron apertures is, on average, within 20$\%$ for filters at $\lambda>2\mu m$. It is larger for the filters at shorter wavelengths, where however all objects have a S/N$<$2. No clear difference is present between the different pointings.

We then repeat the SED fitting described in Section 3.1 using the 0.2”-radius aperture and the comparison between redshift, stellar mass, SFR, A(V) and age is reported in Figure 11. Looking at the redshift, we have that 57$\%$ have a normalised redshift difference (i.e., $\delta z=|z_{kron}-z_{0.2"}|/(1+z_{kron})$) above 0.15 and 30$\%$ have a redshift difference above $|z_{kron}-z_{0.2"}|=1$. These are mainly galaxies for which a double or large peak is present in the probability redshift distribution, so small differences on the photometry results change the probability of the different peaks. For galaxies for which the redshift is consistent within the two photometric sets, the stellar mass is consistent within 0.29 dex, the SFR within 0.49 dex and the dust extinction within 0.10 mag.

Even considering the differences in the physical properties of each single object derived changing from Kron-aperture to 0.2”-radius photometry, we have that the median properties of the F200W-sample remain similar and the conclusion of our work are unaffected. Indeed, the redshift distribution of the F200W-dropouts galaxies is similar to the one presented in the Section 4, with 81$\%$ of them being at $z<2$ with a median stellar mass of 10${}^{6.9}\,\rm M_{\odot}$ and a median dust extinction of 5.0 mag (Fig. 12). At the same time, the majority of the F200W-dropout galaxies are star-forming galaxies, with still 26$\%$ of SB at $z<2$, compared to 27$\%$ derived with the Kron photometry.

We can therefore conclude that, when changing the aperture to calculate the fluxes from a Kron ellipse to a fixed 0.2”-radius circle aperture, there are difference, even large ones, between the physical properties of the objects in our F200W-dropout sample. However, the median properties of the sample are not largely affected by this fluxes changes, un-affecting the main findings of our work.

Appendix B Testing the impact of the dust extinction law, PAH strength and star-formation history

In this Appendix we investigate the impact of the chosen dust extinction law, i.e. Calzetti et al. (2000), the PAH strength and star-formation history (SFH) on the results of our paper. For this reason we repeat the BAGPIPES run described in Section 3.1, but first assuming a Small Magellanic Cloud (SMC) extinction curve (Gordon et al. 2003), second fixing the dust mass associated to PAH to 10% and third assuming an exponentially declining (i.e. SFR$\propto e^{-(t/\tau)}$) SFH.

We show the comparison of the derived redshift, stellar mass, SFR and AV in figure 14 for the test on the dust extinction law, in figure 15 for the test on the PAH strength and in figures 16 and 17 for the test on the SFH. For 86/133 ($\sim 65\%$) of the objects the change in the dust extinction law has a minor impact on the redshift estimation, as $\delta z=|z_{Calzetti}-z_{SMC}|/(1+z_{Calzetti})<0.15$. For these objects the other physical properties are also quite stable, with absolute differences in the stellar mass, SFR and AV of 0.25 dex, 0.36 dex and 0.12 mag, which are generally equal or smaller than the differences induced by changing from Kron to 0.2”-radius aperture photometry (see Appendix A). The median value of the optical dust extinction remains large, varying from 4.9 mag using a Calzetti et al. (2000) reddening law to 4.5 mag considering a SMC reddening law. At the same time, almost all galaxies, 130/133, are still very dusty sources as A${}_{V}\geq$1.

Imposing a low mass content associated to PAH as a similar low impact on the results. Indeed, 85/133 objects have a redshift that change less than $\delta z=|z_{free}-z_{10\%}|/(1+z_{free})<0.15$. The median optical dust extinction is 4.2 mag, with only two objects with A${}_{V}\leq$1.

Finally, we change the SFH considering separately an exponentially declining one and a burst one. The first has a larger impact on the redshift, with only 65/133 objects with a redshift change less than $\delta z=|z_{del}-z_{dec}|/(1+z_{del})<0.15$. However, the sample is still composed at 78$\%$ by galaxies with $z<2$ with median stellar mass of $10^{7.0}\,\rm M_{\odot}$ and median dust extinction $\rm A_{V}=4.7$. At the same time, 131/133 galaxies are still very dusty sources with A${}_{V}\geq$1. The available data is not enough to discriminate between the two SFHs, as the difference in the $\chi_{red}^{2}$ is below one for all objects.

When we consider a single burst of star-formation, the fit is noticeably worst (i.e. $\Delta\chi^{2}_{red}>1$) than the case with a delayed SFH for 26% of the sample and generally results on a larger $\chi^{2}$ for 86% of the sample. Therefore, for at least part of the sample, a burst SFH is generally disfavoured by the data. In Figure 17 we compare the physical parameters for the objects for which the two SFH produce comparably good fits, i.e., $\Delta\chi^{2}_{red}<1$. As for the previous case, there is a large scatter on the derived physical properties of each single objects, however the sample still remains dominated by $z<2$ galaxies, i.e. 84$\%$ of the sample even when considering a burst SFH. These galaxies at $z<2$ have a median dust extinction A${}_{V}=4.8$ and median stellar mass $log(\rm M/M_{\odot})=7.0$, with only two of these galaxies having a A${}_{V}<$1.

Overall, we can conclude that the physical properties of each galaxy are quite uncertain, but the results presented in the paper are statistically robust.

Appendix C Possible brown dwarfs

In this appendix we report the photometry, cutouts, best models and redshift probability distributions of the two galaxies for which the brown dwarf templates correspond to a lower $\chi_{2}$ with respect to the galaxy ones. Both these galaxies are in pointing 2 (N2-2160 and N2-4539). The best brown dwarf fits of all two objects correspond to a model with gravity of $\rm 10^{5.5}cm\,s^{-2}$ and a metallicity of [Fe/H]=-0.5. The effective temperature are instead 1000 K and 1300 K for the N2-2160 and N2-4539, respectively. At the same time, the normalisation of such models would place them at 2.7 and 4.4 pc from the Sun, respectively.

Looking instead at the fit with galaxy models, the N2-2160 has a main peak at $z\leq 1$, considering both the flux derived with the Kron aperture and the one with a fixes 0.2”-radius aperture. The possible redshift of N2-4539 is more uncertain, as the two photometric sets give two different best solutions, i.e. $z=$0.7 and $z=$5.9.

Both objects are unresolved, but this is consistent both the brown-dwarf solution and the galaxy one, considering the stellar-mass-size relation by van der Wel et al. (2014) and their stellar masses (i.e.$10^{6}-10^{7}\rm M_{\odot}$).

Appendix D SED fits and cutouts

In this section we report the SED fits and cutouts images for all F200W-dropout images. The full list of images are available online.

{longfigure}

cc

&

Same as Figure 3, but for the F200W-dropout galaxies not previously shown.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.