Unveiling the hidden universe with JWST: The contribution of dust-obscured galaxies to the stellar mass function at z ∼ 3 − 8

The emergence of massive, optically-faint galaxies in infrared observations has revealed that our view of the high-redshift Universe was previously incomplete. With the advent of JWST, we can for the first time probe the rest-frame optical emission of galaxies at 𝑧 > 3 with high sensitivity and spatial resolution, thus moving towards a more complete census of the galaxy population at high redshifts. To this end, we present a sample of 148 massive, dusty galaxies from the JWST/CEERS survey, colour-selected using solely JWST bands. With deep JWST/NIRCam data from 1.15 𝜇 m to 4.44 𝜇 m and ancillary HST/ACS and WFC3 data, we determine the physical properties of our sample using spectral energy distribution fitting with BAGPIPES. We demonstrate that our selection method efficiently identifies massive ( ⟨ log M ★ / M ⊙ ⟩ ∼ 10) and dusty ( ⟨ A V ⟩ ∼ 2 . 7 mag) sources, with a majority at 𝑧 > 3 and predominantly lying on the galaxy main-sequence. The main results of this work are the stellar mass functions (SMF) of red, optically-faint galaxies from redshifts between 3 < 𝑧 < 8: these galaxies make up a significant fraction of the pre-JWST total SMF at 3 < 𝑧 < 4, and dominate the high-mass end of the pre-JWST SMF at 4 < 𝑧 < 6 and 6 < 𝑧 < 8, suggesting that our census of the galaxy population needs amendment at these epochs. While larger areas need to be surveyed in the future, our results suggest already that the integrated stellar mass density at log M ★ / M ⊙ > 9 . 25 may have been underestimated by ∼ 20-25% at 𝑧 ∼ 3 − 6, and ∼ 110% at 𝑧 ∼ 6 − 8.


INTRODUCTION
For decades, observational astronomers have been on a quest to determine how the galaxy population evolves through cosmic time.The Hubble Space Telescope (HST) has pioneered the study of this ★ E-mail: rashmi.gottumukkala@gmail.com question: HST has observed high-redshift galaxies, primarily through their rest-frame ultraviolet (UV) emission.These high-redshift galaxies, usually referred to as 'normal' or Lyman-break galaxies (LBGs), have been studied extensively from z ∼ 3 to z ∼ 11, tending to have moderate star formation rates (SFR) and stellar masses, and are thought to make up the bulk of the galaxy population (e.g., Labbé et al. 2013;Schaerer et al. 2013;Bouwens et al. 2015;Finkelstein et al. 2015;Oesch et al. 2016;Faisst et al. 2020).These mostly dust un-obscured galaxies are also thought to dominate the cosmic starformation rate density (SFRD) at z > 4, while at lower redshifts the universe was dominated by obscured star-formation (e.g., Madau & Dickinson 2014;Zavala et al. 2021).While 'normal', un-obscured galaxies have been well-studied, our census of the galaxy population remains incomplete at z > 3 as rest-frame UV selections systematically miss massive, dust-obscured sources (e.g., Alcalde Pampliega et al. 2019;Wang et al. 2019).
Over the last decade, a significant population of optically undetected galaxies with relatively bright infrared (IR) or sub-millimetre (sub-mm) emission has been discovered, several in Spitzer/IRAC data and some of them with ALMA detections (e.g., Huang et al. 2011;Simpson et al. 2014;Caputi et al. 2015;Stefanon et al. 2015;Wang et al. 2016;Franco et al. 2018;Alcalde Pampliega et al. 2019;Wang et al. 2019;Yamaguchi et al. 2019;Williams et al. 2019;Dudzevičiūtė et al. 2020;Sun et al. 2021;Smail et al. 2021;Manning et al. 2022;Shu et al. 2022;Xiao et al. 2023b).They typically have very red spectral energy distributions (SEDs) and remain undetected even in deep HST -band observations -hence their name: HST-dark galaxies.Their SEDs are not well-constrained, with a few photometric detections and lack of spectroscopic redshifts, which result in very large uncertainties on their photometric redshifts, stellar masses, and SFRs (e.g., Caputi et al. 2012;Stefanon et al. 2015;Williams et al. 2019;Alcalde Pampliega et al. 2019).The physical properties of these galaxies were largely unconstrained until the arrival of the James Webb Space Telescope (JWST, Gardner et al. 2023).
JWST has revolutionised the field of optically-faint galaxies, providing for the first time reliable physical parameters (e.g., Barrufet et al. 2023;Nelson et al. 2023;Pérez-González et al. 2023;Rodighiero et al. 2023;Labbé et al. 2023b;Gómez-Guĳarro et al. 2023).With its unprecedented sensitivity and resolution in the near-IR, JWST probes the rest-frame optical emission of galaxies at  ≳ 3, allowing one to identify the Balmer break, a good redshift and mass indicator.Additionally, the SEDs of massive galaxies are typically highly dust-attenuated with characteristic red slopes in the rest-frame optical.With its extensive photometric coverage from 1 − 5 m, JWST's Near-Infrared Camera (NIRCam; Rieke et al. 2023) is the ideal instrument to identify sources based on these features.
The early JWST era has seen the puzzling emergence of two additional populations of galaxies.The first is a population of massive sources (> 10 10 M ⊙ ) at  > 7, less than 700 Myr after the Big Bang (e.g., Labbé et al. 2023b).With the currently accepted theory of hierarchical structure formation within ΛCDM cosmology, it is challenging to explain how galaxies could accumulate this much mass through mergers or accretion alone (Boylan-Kolchin 2023; Menci et al. 2022), while it might still be possible to reconcile such observations with theory (Mason et al. 2023;Dekel et al. 2023).One possibility is that these sources are actually active galactic nuclei (AGN), with one Labbé et al. (2023b) source being spectroscopically confirmed to be an AGN with broad emission lines (Kocevski et al. 2023).A deeper investigation into massive galaxies in the early universe is needed in order to determine their abundance and place constraints on mass assembly.
The second emergent population consists of massive quiescent galaxies at high redshifts, now spectroscopically confirmed up to  = 4.658 (Carnall et al. 2023).Relatively little physical insight has been provided by simulations thus far to explain the emergence of quiescent galaxies at  > 3, with simulations struggling to predict observed number densities (Valentino et al. 2023;Gould et al. 2023).While it is highly likely that sub-millimetre galaxies (SMGs) evolved into massive quiescent galaxies at  ∼ 2 (Toft et al. 2014), their num-ber densities are insufficient to explain the presence of quiescent galaxies at  ∼ 3 − 4 (Valentino et al. 2020(Valentino et al. , 2023)).Hence, an important step towards understanding the emergence of quenched galaxies is to look for previously-missed massive, dusty galaxies in the early universe and determine their stellar masses and abundances.
For the study of galaxy abundances, the stellar mass function (SMF) is an extremely useful statistical tool to quantify the evolution of the galaxy population as a function of stellar mass across cosmic history.Determining the SMF at various epochs in the history of the universe allows us to track early galaxy build-up.Several studies have so far constrained high- SMFs with ground-and space-based multiwavelength observations (e.g., Davidzon et al. 2017;Stefanon et al. 2015Stefanon et al. , 2017bStefanon et al. , 2021;;McLeod et al. 2021;Santini et al. 2021;Weaver et al. 2023b;Navarro-Carrera et al. 2023), with the shape of the total SMF being found to be accurately described by the empirically motivated Schechter (1976) function.Given that JWST is primed to find massive, dust-obscured sources that have previously been missed in the galaxy census, this raises the question of whether or not the total SMF at high- epochs requires modification.The central question we aim to address with this work is, 'How do massive, dusty galaxies selected with JWST affect the high-mass end of the galaxy stellar mass function in the early universe?' In this study, we use data from the Cosmic Evolution Early Release Science (CEERS) survey (Finkelstein et al. 2022(Finkelstein et al. , 2023)), a JWST Cycle 1 community survey in the CANDELS/EGS field.CEERS is aimed at discovering the first galaxies and observing galaxy assembly at  > 3. Given its deep photometric coverage with JWST/NIRCam from 1.15m -4.44m, CEERS is the ideal survey to look for red, IR-bright galaxies.
This paper is structured as follows: In Section 2, we discuss the photometric data used from HST and JWST and the production of the HST-JWST merged photometric catalogue.We introduce our colour selection using photometry solely from JWST.Furthermore, we describe how we create an AGN-cleaned sample of purely starforming galaxies.In Section 3, we explain the SED fitting performed using the Python tool BAGPIPES (Carnall et al. 2018).In Section 4, we discuss the physical properties of our sample, and situate our galaxies on the galaxy main-sequence (4.3).In Section 5, we discuss the methodology used to compute the SMFs (5.1) and present the SMFs of massive, dusty galaxies at 3 <  < 4, 4 <  < 6, and 6 <  < 8 (5.2).Finally, we discuss our sample in the context of other JWST studies in Section 6, and we summarise and conclude our study in Section 7.
For this work, we assume a flat ΛCDM cosmological model with H 0 = (67.8± 0.9) km s −1 Mpc −1 and Ω  = 0.308 ± 0.012 as found by the Planck Collaboration et al. (2016).All magnitudes are quoted in the AB magnitude system (Oke & Gunn 1983).Throughout this paper, we use a Kroupa (2001) initial mass function (IMF).If required for comparison, we scale mass values used in the literature from Salpeter (1955) or Chabrier (2003) to Kroupa (2001) using the scale factors quoted in Madau & Dickinson (2014).

OBSERVATIONS AND SAMPLE SELECTION
In this section, we describe the imaging data used in this work and the production of the HST and JWST merged photometric catalogue for the CEERS field.In addition, we present our sample selection criteria in order to identify massive and dusty galaxies, including the colour-selection we develop as well as the criteria used to identify and remove AGN from our final sample.

Imaging data
We use data from the Cosmic Evolution Early Release Science (CEERS) programme, one of JWST's first early-release science surveys in Cycle 1, with data collected in June and December 2022 (Finkelstein et al. 2022(Finkelstein et al. , 2023)).CEERS comprises 10 NIRCam pointings covering ∼100 arcmin2 in the Extended Groth Strip (EGS) field, a CANDELS legacy field containing a wealth of ancillary HST multiwavelength data.The NIRCam data covers a range of wavelengths from 1.15m to 4.44m in the following filters: F115W, F150W, F200W, F277W, F356W, F410M, and F444W (where W and M indicate a wide or medium band filter).Ancillary HST data from the ACS imager is available at wavelengths between 435nm to 814nm (in 3 filters: F435W, F606W, and F814W) and from the WFC3 imager at wavelengths between 1.05m to 1.60m (in 4 filters: F105W, F125W, F140W and F160W) (Koekemoer et al. 2011;Grogin et al. 2011;Stefanon et al. 2017a).
For this work, we use the version 5 images reduced with the grizli pipeline and made publicly available by G. Brammer1 , following the same steps as outlined in Valentino et al. (2023).The images include all available data over these fields taken with HST and JWST.The imaging depths as measured in circular apertures with a radius of 0.16 ′′ are listed in Table 1.They vary between 28.6 mag to 29.2 mag in the JWST wide filters and are ∼ 28.3 mag in the shortest wavelength ACS imaging.

Production of the HST-JWST photometric catalogue
We use the JWST and ancillary HST images to create photometric catalogues, taking into account the wavelength-dependent pointspread function (PSF).In the following, we briefly describe how the PSF-matched photometric catalogue used in this work was produced (see Weibel et al. in prep. for details).
We match the fluxes in all HST+JWST filters to the PSF resolution in the reddest JWST/NIRCam filter, F444W.For the NIRCam and WFC3 filters, we use the PSFs provided by G. Brammer for use with the CEERS grizli mosaics (Brammer 2018)  2 .
For the ACS filters, we derive effective PSFs from the science images by first identifying bright, but unsaturated stars without bright neighbouring sources or flagged pixels, from a preliminary SourceExtractor run (Bertin & Arnouts 1996).Then, we use the method EPSFBuilder from the python package photutils (Bradley et al. 2022) which is based on the model developed by Anderson & King (2000) to obtain the final effective PSFs.
We compute matching kernels from each ACS and NIRCam PSF to the NIRCam/F444W PSF using the software package pypher (Boucaud et al. 2016) and convolve each flux and root mean square (rms) image with the respective kernel to match the PSF resolution in F444W.
We follow a different procedure for the WFC3 filters because their PSFs are broader than the NIRCam/F444W PSF.First, we compute matching kernels from all of them and from the F444W PSF to the WFC3/F160W PSF, in the same way as described above, and produce PSF-matched flux and rms images accordingly.
Then, we run SourceExtractor in dual mode, using an inversevariance weighted stack of the unaltered F277W+F356W+F444W images as the detection image and measuring fluxes in circular apertures with a radius of 0.16 ′′ on the original images, the images that were PSF-matched to F444W as well as the images that were PSFmatched to F160W.For the final catalogue, we use the flux measurements on the original image in F444W and those on the images PSF-matched to F444W for all other filters, except the WFC3 data.
For the latter, we correct the fluxes measured on the original images to match the colour between the respective filter and F444W as measured on the images PSF-matched to F160W.We scale all fluxes to the flux measured in Kron-like apertures by SourceExtractor in F444W, obtained using the default Kron parameters 2.5 and 3.5.To account for residual flux outside the Kron aperture, we measure the fraction of the energy enclosed by a circular aperture with a radius of √   kron_radius, where ,  and kron_radius characterise the Kron-ellipse, on the theoretical F444W PSF obtained from webbpsf, and divide all fluxes by that fraction.Finally, we correct all fluxes for Milky Way foreground extinction using the extinction model from Fitzpatrick & Massa (2007) through the python package extinction, using the E(B-V) map outlined in Schlafly & Finkbeiner (2011).
To get a more realistic estimate of the rms uncertainty of our flux measurements that accounts for correlated noise, we put down circular apertures with a radius of 0.16 ′′ in 5000 random positions on the "signal-to-noise" image (i.e., the flux image divided by the rms image).We multiply the uncertainties on all fluxes, measured from the rms map respectively, by the scatter measured among those apertures.This leads to a scaling of the flux uncertainties by ∼5 -∼35% depending on the filter -the largest correction being applied to F115W and the smallest to F444W.
To identify and flag stars we used a flux ratio criterion similar to Weaver et al. (2023b).We also flag objects as artefacts that are too small to be real sources (typically left-over bad pixels).The full CEERS catalogue contains 93,922 sources.Out of these, we remove 930 sources that are either identified as stars or flagged as artefacts based on the above criteria, resulting finally in 92,992 sources. .The colour criterion F150W-F444W>2.1 mag (black dashed line) in principle identifies  > 3 sources with A V ≳ 2 mag and log(M ★ /M ⊙ ) ∼ 10, while the magnitude cut F150W>25 mag (dark red dashed line) is designed to rid the sample of low- sources while retaining the most massive galaxies in the sample (∼ 10 11 M ⊙ ).Also shown for reference is the F150W = 26 mag cut that is a proxy for identifying HST-dark sources (Pérez- González et al. 2023).We select 179 red galaxies with these criteria that theoretically restrict our sample to massive, dusty, high-redshift galaxies.

Selection of red, optically-dark/faint sources at z>3
Over the last decade, numerous studies of HST-dark galaxies and red galaxies have been conducted, with dropout and colour selections shown to be effective methods for selecting high-redshift sources.
Typically, these studies combine HST and Spitzer data to select massive and dusty star-forming galaxies (e.g., Huang et al. 2011;Alcalde Pampliega et al. 2019;Wang et al. 2019;Sun et al. 2021).Several unique colour cuts have been used over the last decade for efficient selections of red galaxies using HST/WFC3 bands in the optical and Spitzer/IRAC and (recently) JWST/NIRCam bands in the near-IR (e.g., Huang et al. 2011;Caputi et al. 2012;Wang et al. 2016;Alcalde Pampliega et al. 2019;Wang et al. 2019;Sun et al. 2021;Barrufet et al. 2023;Labbé et al. 2023b;Nelson et al. 2023;Rodighiero et al. 2023;Pérez-González et al. 2023;Xiao et al. 2023b;Long et al. 2023a).Here, we build on these and make a broad selection of red galaxies using solely JWST/NIRCam bands in order to fully exploit the increased sensitivity and resolution of JWST.By designing and implementing a colour selection capable of identifying the effects of the Balmer-break and reddened stellar continuum emission in a galaxy's photometry, we expect to select massive and dusty galaxies at high redshifts.For this, we use the Python tool Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (BAGPIPES, Carnall et al. 2018)  3 to investigate the evolution of colour with redshift.We generate galaxy spectra, from 3 https://bagpipes.readthedocs.io/en/latest/which we extract the photometry and compute modelled colours.We use a delayed- star-formation history, ages of 1 Gyr, an -folding time of 3 Gyr, a mass of 10 10 M ⊙ and metallicity of 0.5 Z ⊙ .We model galaxies at redshifts between  = (1., 6.) in steps of Δ = 0.1 and at discrete dust attenuation values of A V = [2.0,3.0, 4.0] mag using a Calzetti dust model (Calzetti et al. 2000) to produce the SED tracks of massive, dusty galaxies as shown by the coloured lines in Figure 1.We also model a single SED track of a 10 11 M ⊙ massive galaxy with A V = 4 mag.
As the Balmer break gets redshifted beyond 1.5m at  ≳ 3, we design a colour cut that requires galaxies to be faint in F150W in comparison to longer wavelength bands.Pérez- González et al. (2023) show with a JWST-selected sample that HST-faint sources extend to higher masses than HST-dark sources.We therefore move beyond the strict HST-dark classification by including HST-faint sources in our selection, so as not to miss the most massive and bright galaxies (HST-dark classification referenced from Pérez- González et al. 2023).We also use the F444W band to get the broadest redshift range possible (as the highest redshift sources will have their Balmer break closer to F444W).Given our choice of using the F150W and F444W bands, we identify the F150W -F444W colour at which we expect to select galaxies that are (i) high redshift ( ≳ 3), (ii) massive (log M ★ /M ⊙ ∼ 10), and (iii) dusty (A V ≳ 2 mag).In addition, from the SED-tracks shown in Figure 1, we estimate the F150W magnitude at which we rid the sample of low- sources ( ≲ 2) while retaining the most massive and dusty galaxies in our sample.
Using the SED modelling described above, we determine a selection that is optimised to identify galaxies with A V ≳ 2 mag and log M ★ /M ⊙ ∼ 10 at  ≳ 3, described in Equation 1: F150W > 25 mag. (2) Additionally, as the prominent feature of our galaxies is their redness, this suggests that they must have significant emission in the long wavelength bands.To ensure reliable detections, we require SNR > 5 in all three wide filters in the long wavelength channels: F277W, F356W, and F444W.Altogether, this colour selection is more flexible than in previous studies (i.e.Barrufet et al. 2023); we later remove the  < 3 sources after evaluating their physical properties (see Section 4.1).
We find 179 galaxies that satisfy the F150W − F444W > 2.1 mag and F150W > 25 mag criteria out of the >90,000 sources in our catalogue (see Figure 1).

Identifying and removing obscured AGN
In recent literature, there has been mounting evidence from JWST of a population of high redshift obscured AGN that displays very red colours in the NIR (Labbe et al. 2023a;Matthee et al. 2023;Barro et al. 2023;Greene et al. 2023).These so-called little red dots (LRDs) have characteristically blue rest-UV colours which possibly arise from star-forming regions, and red rest-optical colours that arise from the hot, dusty torus of the AGN (Labbe et al. 2023a;Greene et al. 2023).These sources are potential contaminants in selections of red, star-forming galaxies, and it is important to address their presence in our sample.
Based on the colour and compactness criteria outlined in Labbe et al. (2023a) and Greene et al. (2023), we identify a parent sample of 29 potential AGN candidates.In order to further identify point-like sources, we perform a two-component PSF+Sérsic fit in the F444W Postage stamps and SED fits of four selected galaxies from our sample of red galaxies.The stamps boxed in blue are from ancillary HST/ACS and WFC3 data, and the stamps boxed in red are from JWST/NIRCam imaging (each stamp is 4 × 4 arcsec 2 ).There is a variety in the morphological properties of our sample, ranging from spatially extended sources to compact ones.The lower panels display the SED fits: the maroon points represent the photometry and the downward arrows represent the flux upper limits.The orange lines are the SED fits from BAGPIPES and the photometric redshift probability density functions are inlaid in the lower right part of the graphs.The physical properties of these galaxies are quoted on the graphs.They are massive (log M ★ /M ⊙ ≳ 9.5) and dusty (A V ∼ 1.5 − 4 mag) with redshifts ranging from  ∼ 3 − 8. filter using the GalfitM4 (Häußler et al. 2013;Vika et al. 2015) software, identifying sources where the flux associated with the PSF component exceeds the flux associated with the Sérsic component (Labbe et al. 2023a).We identify 20 sources that satisfy these criteria.
We remove these 20 sources from our sample during analysis (Section 4 onwards), thus considering a purely star-forming sample of galaxies.Figure A1 in Appendix A shows the postage stamps and SED of source 6583, identified as one of the 20 AGN candidates in our sample selection.Figure A2 shows the effect of AGN on the SMF, showing that in particular the SMF at 6 <  < 8 is significantly overestimated by including AGN.

SED FITTING WITH BAGPIPES TO DETERMINE THE PHYSICAL PROPERTIES OF GALAXIES
To calculate the physical properties of our sample, we use the Python tool BAGPIPES (Carnall et al. 2018).BAGPIPES is an SED-fitting tool capable of modelling galaxies with various star formation histories (such as delayed-, exponential, constant, bursts, etc.) and dust models (Cardelli et al. 1989;Calzetti et al. 2000;Charlot & Fall 2000, etc.), using stellar population synthesis (SPS) models (Bruzual & Charlot 2003).We choose to use a delayed- SFH, which has been shown as an effective SFH to model the bulk of the stellar population, and accurately recover stellar masses (Ciesla et al. 2017).Furthermore, this SFH has been successfully used in previous studies of HST-dark galaxies and massive galaxies (Wang et al. 2016(Wang et al. , 2019;;Alcalde Pampliega et al. 2019;Pérez-González et al. 2023;Barrufet et al. 2023).
We perform SED-fitting within a broad parameter space, allowing the code to explore the following ranges: redshifts between  = (0, 10), a delayed- SF history with  = (0.1, 9) Gyr, masses in the range log M ★ /M ⊙ = (6, 13), metallicities between Z = (0.2, 1.2) Z ⊙ , a Calzetti dust model with A V = (0.2, 4) mag, nebular emission with an ionisation parameter of log U = −2, and a velocity dispersion of 300.The models chosen have been successfully used for similar types of galaxies, being able to fit red SEDs (Wang et al. 2016(Wang et al. , 2019;;Barrufet et al. 2023).The broad parameter space in each model allows us to explore this enigmatic galaxy population and unveil their physical properties in more detail, in particular their stellar masses.
To test the suitability of our chosen A V range, we allow A V to vary from (0, 6) mag, finding that some galaxies are fit to very dusty (A V > 4 mag) solutions at low redshifts ( < 0.75).Galaxies with similar properties have been reported in Caputi et al. (2012) and more recently in Bisigello et al. (2023), where BAGPIPES is used.We compare the redshifts from BAGPIPES with redshifts derived from the Easy and Accurate Zphot from Yale (EAZY) software (Brammer et al. 2008), using the blue_sfhz template set5 .We find that the photo-'s of the A V > 4 mag sources are not in good agreement with EAZY, where EAZY typically finds higher- solutions with lower A V .This is expected, as the maximum A V that EAZY can describe is redshift dependent, reaching a maximum of A V ∼ 4 at  ∼ 3.In addition, upon visual inspection of the postage stamps, we find that several of these sources are very compact, completely dropping out of the shorter wavelength filters and thus being more likely to lie at higher redshifts than at  < 0.75.The inclusion of MIRI data could potentially rule out the low- solutions.However, this is only available over a very small portion of the field currently.We refer to Alcalde Pampliega et al. in prep.for a more detailed analysis including MIRI data.
Additionally, given that our aim is to derive accurate stellar masses in order to calculate the stellar mass function, we test whether the derived stellar masses change significantly if we use the EAZY photometric redshifts as an input to the BAGPIPES SED fitting.We find that with EAZY- as an input, ⟨log M ★ /M ⊙ ⟩ = 10.18 +0.40  −0.50 and with the BAGPIPES-, ⟨log M ★ /M ⊙ ⟩ = 10.15 +0.43  −0.50 .Both derived stellar masses follow a tight 1:1 relation with an average scatter of 0.2 dex, suggesting that the final stellar mass functions will not be strongly affected by our choice of input redshift.We finally use the BAGPIPES- in all SED fittings.
Examples of some SED fits are shown in Figure 2, which showcases the variety in galaxy morphology and physical properties.Most of our sources have very red slopes indicating high dust attenuation.We find a diversity in morphology: some sources are spatially extended, while others are extremely compact (see Figure 2).
We performed a visual inspection of SEDs and postage stamps for all sources while considering their derived physical properties.We remove 11 sources from our sample due to either clearly overestimated photometric redshifts and masses (spatially extended sources that are likely at lower redshift) or sources with deblending issues.Our final sample thus contains 148 galaxies.
To recapitulate, out of the colour-selected sample of 179 galaxies outlined in Section 2.3, we remove 20 AGN candidates (described in Section 2.4) and further remove 11 sources that have poor SED fits, resulting in a final sample of 148 galaxies.

PHYSICAL PROPERTIES OF RED, OPTICALLY-FAINT GALAXIES
JWST's outstanding sensitivity and resolution in the near-IR allow us to determine photometric redshifts and physical parameters (such as stellar masses, star formation rates, etc.) with unprecedented accuracy.This will allow us to place tighter constraints on the stellar mass build-up in the early universe.In this section, we present the photometric redshifts and physical characteristics of our galaxies as determined with BAGPIPES (see Table B1 for a list of the derived physical parameters of our full sample; AGN candidates are denoted as such but removed from the following analysis).

Photometric redshifts
We determine photometric redshifts for our sample of red galaxies using BAGPIPES (see Section 3).The redshift distribution is shown in Figure 3. ∼60% of the sample lies at  ≳ 3 and ∼ 90% at  ≳ 2, with an average redshift of  mean = 3.46.This shows that our colour selection successfully identifies high redshifts galaxies, out to  ∼ 8.The redshift is mostly in agreement with EAZY redshifts using standard templates.We note the significant number of galaxies that lie at  < 3 in our selection.We draw the reader's attention back to Figure 1, where we show with the use of arrows that masses and redshifts increase in opposing directions in the colour space of our selection.Further, at a given redshift and stellar mass, there is a scatter in stellar ages and dust which means that invariably, there is a scatter in the properties of the selected population.Therefore, in order to build the most inclusive sample and so as not to miss the most massive and dusty galaxies, it is unavoidable for low-redshift galaxies to enter our selection.
Further, we draw the reader's attention to a caveat of this selection technique, namely the two local peaks seen in the redshift distribution at  ∼ 5.5 and  ∼ 7.5 in Figure 3.The F444W detection is likely driven by the H+[NII] lines at  ∼ 4.9 − 6.6, and the [OIII]+H lines at  = 6.9 − 9.0 (see Oesch et al. 2023).The samples at these redshifts are thus qualitatively different from the bulk sample because their 'redness' comes from emission lines rather than the continuum.However, we note that our selection includes 5 detection masks in the long-wavelength filters (F277W, F356W and F444W), thus ensuring that the continuum is relatively bright over an extended wavelength range and not just in F444W.Additionally, as elaborated in the following section, all sources in this study have high A V magnitudes; therefore, even if the F444W fluxes of a few select sources are slightly boosted by emission lines, they still qualify as targets for our study.Furthermore, the red and optically faint selection criteria imply that such sources were missing from previous estimates, further justifying their inclusion in our sample.

Physical properties of red galaxies
One of JWST's most important improvements in the NIR is its increased photometric coverage at 1-5 m in comparison with its predecessor, Spitzer.This allows JWST to better probe the Balmer break and thus derive more accurate photometric redshifts than previously possible.With more accurate photometric redshifts, through SEDfitting we can additionally derive more reliable estimates of the stellar masses of galaxies and their star-formation rates.
We present the distributions of the physical properties of our sample of 148 red, optically-faint galaxies in Figure 4. We find these galaxies to be massive, with a median stellar mass of ⟨log M ★ /M ⊙ ⟩ = 10.15 +0.43  −0.50 .They also have high dust attenuations of ⟨A V ⟩ = 2.71 +0.88  −0.91 mag.Additionally, they have moderate star-formation rates, with ⟨log SFR/M ⊙ yr −1 ⟩ = 1.64 +0.43  −0.68 and ⟨sSFR/Gyr −1 ⟩ = 2.66 +4.72  −1.42 .As expected, we find our sample to be dominated by relatively massive and dusty star-forming systems.
We note that the SFRs derived in our study are based on restframe UV to optical SED fits.We are therefore not modelling the starlight that is reprocessed by dust and emitted in the FIR.While for a more complete picture of the SFR, more FIR data is needed to recover the full infrared SED (see e.g.Xiao et al. 2023a), we refer the reader to Williams et al. (2023), where they show with a selection of optically-dark galaxies that only those with the most extreme SFRs are significantly affected by the inclusion of MIR and FIR data.
To further illustrate the dusty nature of our galaxies we situate them on the widely used UVJ diagram.We classify the star-forming vs. quiescent regions on the UVJ diagram following Williams et al. (2009), and further split the star-forming region into dusty and unobscured zones following the classification in Spitler et al. (2014).Figure 5 shows the UVJ classification of our galaxies and of the full CEERS sample.Rest-frame colours for our red galaxies are determined by the best-fit SEDs from BAGPIPES, while for the full CEERS sample they are determined with EAZY due to less expensive computational time.Except for one galaxy lying in the quiescent region of the diagram, the sample lies in the star-forming region, with ∼ 75% of the sample lying particularly in the dusty region.Thus the UVJ classification further indicates the dust-obscured nature of our sample.

Red galaxies on the galaxy main-sequence
To place our galaxies within the context of galaxy evolution, we explore their position on the galaxy main-sequence.Figure 6 shows a plot of SFR vs. M ★ for our sample, comparing them to the starforming main-sequence (MS) of galaxies at = 2, 4, and 6 (from Speagle et al. 2014).As shown, our galaxies lie on the star-forming MS, indicative of the 'normal' nature of their ongoing star-formation.The three galaxies lying significantly below the SF main-sequence are candidate quiescent galaxies at  < 3.They form less than 2% of our sample.
We compare our galaxy sample at 3 <  < 5 to two studies of There is a clear tendency for less dusty sources to be at higher redshift.The majority of our sample lies on the MS at redshifts of  < 6, suggesting that they are normal star-forming galaxies with moderate SFRs.The three sources lying significantly below the MS are candidate quiescent galaxies.interest from the literature: Wang et al. (2019), that studied ALMAdetected HST-dark galaxies using HST and Spitzer, and the more recent Barrufet et al. (2023), that studied HST-dark galaxies with HST and JWST.The comparison between samples is shown in Figure 7.
The high-mass end of our sample overlaps with the Wang et al. ( 2019) sample as our colour selection is inclusive of the Wang et al. (2019) selection criteria.Additionally, the Spitzer/IRAC sensitivity is considerably lower than that of JWST/NIRCam in the same range, thus resulting in the detection of only the brightest and most massive galaxies.We also select lower-mass galaxies than Wang et al. (2019) as JWST can detect galaxies that are fainter in F444W, and our colour selection is less extreme than that used in Wang et al. (2019).
The low-mass end of our sample overlaps with the range covered by HST-dark galaxies from Barrufet et al. (2023).This study specifically looked at HST-dark galaxies (F160W > 27 mag), with JWST/NIRCam's sensitivity permitting detections of lower-mass systems.However, this magnitude cut also limits the detection of brighter, higher-mass sources.By using a less restrictive magnitude cut at 1.5 m (F150W > 25mag) our selection criteria ensure we find higher-mass galaxies than in Barrufet et al. (2023) while still including the lower-mass HST-dark galaxies in their study.In Figure 7, we show that our sample of red, optically-faint sources lie on the galaxy main-sequence, similar to HST-dark galaxies (Barrufet et al. 2023).The comparison with these select studies from the literature shows that the mass-range spanned by our sample overlaps with both pre-JWST and JWST-selected HST-dark/faint galaxies.

STELLAR MASS FUNCTIONS OF RED GALAXIES: FINDING THE MISSING SOURCES THAT DOMINATE THE HIGH-MASS END
In this section, we present the SMFs of red, optically-faint galaxies at redshifts of 3 <  < 8.We describe the method used to derive the SMFs and their uncertainties.The SMFs are then presented, discussed, and compared to studies in the literature.We note that the sample statistics quoted in the previous section were for the full sample of 148 red galaxies across the whole redshift range shown in Figure 3.For the 86 galaxies in the redshift range 3 <  < 8, the average stellar masses and dust attenuation values are ⟨log M ★ /M ⊙ ⟩ = 10.17 +0.41  −0.56 and ⟨A V ⟩ = 2.30 +1.22 −0.56 mag, setting the stage for the exploration of the SMFs of massive, dust-obscured galaxies at these epochs.

Determining SMFs
We use the step-wise method to calculate the SMFs of our sample (Bouwens et al. 2008;Santini et al. 2021).The SMFs are approximated by binning the mass distribution, calculating the number of galaxies within each mass bin and dividing this number by the differential comoving volume of the survey.The mass resolution is judiciously chosen to have reasonable statistics within individual mass bins and to have an appropriate mass resolution in order to determine the shape of the SMF.
Given that we detect sources on a stacked image of F277W+F356W+F444W and additionally select sources based on their F150W-F444W colour, we determine the area overlapped by all four filters in the CEERS survey, which is 83.3 arcmin 2 .We accordingly calculate the differential comoving volume within the considered redshift bins respectively.
The final SMFs are calculated as shown in Equation 3, where Φ ,  is the estimated number density in a redshift bin '' and mass bin ' ' per fixed mass bin Δ log M.   is the number of galaxies in the 'th mass bin,  i, comoving is the differential comoving volume determined within the 'th redshift bin and  ★ is a multiplicative factor derived from a completeness simulation used to account for missing sources in our detection catalogues (described in 5.1.1).

Completeness
We measure the source detection completeness by running a simple simulation using our custom version of the publicly available software GLACiAR2 (Carrasco et al. 2018;Leethochawalit et al. 2022).
We first select a representative 1.5 ′ ×1.5 ′ cutout approximately in the middle of the CEERS image with average depth and no contamination by bright stars.Using GLACiAR2, we inject artificial sources, spanning a range of input UV magnitudes from -24.4 to -16.2 in 35 bins at a fixed redshift of  = 6 into the cutout.We inject 500 sources per bin in batches of maximally 100 sources at a time to avoid overcrowding and run SourceExtractor with the same settings as outlined in Section 2.2.The injected galaxies follow a Gaussian distribution in the logarithm of the effective radius, centred at 0.8 kpc and with a scatter of 0.17 dex and they have Sérsic light profiles with 50% of the galaxies having a Sérsic index of 1.5, and 25% having indices of 1 and 2 respectively.We further assume a flat SED (i.e., a fixed UV-slope of  = −2), since we only wish to estimate the completeness as a function of apparent magnitude.We repeat this experiment 10 times, therefore injecting 175,000 sources in total.
To obtain the completeness of our sample, we first measure the fraction of recovered galaxies as a function of the input magnitude.Then, for each bin in apparent output magnitude, we determine the completeness as the weighted mean of the completeness values found in each input magnitude bin, weighted by the number of sources from that bin that were observed in the given output magnitude bin.Then, we additionally determine the fraction of detected sources in each apparent magnitude bin that have a measured SNR > 5 in all of F277W, F356W and F444W (cf.Section 2.3) and multiply that fraction with the detection completeness obtained in the previous step.Since all the observed galaxies considered in this paper have AB-magnitudes ≲ 27 in F444W, they are in a regime where the completeness is high and approximately constant as a function of the apparent magnitude (e.g., in F444W).From our analysis, we derive a mean completeness factor of  ★ = 0.87 by which we scale all our mass functions (see Equation 3).
To determine the mass limit above which we are 80% mass complete, we project our mass distribution onto the SNR limit of our selection (see, e.g., Pozzetti et al. 2010).Given that we select sources that are detected with a 5 certainty in F444W, F356W and F277W, the SNR limit of our selection is SNR lim = 5 √ 3. We calculate the joint SNR for all sources in our sample as SNR 2 joint = SNR 2 F277W + SNR 2 F356W + SNR 2 F444W .Assuming that stellar mass values linearly scale with source brightness, we find the hypothetical mass that each source would have if detected at SNR lim : log M hypothetical = log M ★ − log(SNR joint /SNR lim ).The 80 th percentile of the M hypothetical distribution provides the limit above which the sample is 80% mass complete, given the specific mass-to-light ratios and SEDs in our sample.We determine that the 80% mass complete limits are log M ★ /M ⊙ = 9.15 at 3 <  < 4, log M ★ /M ⊙ = 9.07 at 4 <  < 6 and log M ★ /M ⊙ = 9.21 at 6 <  < 8. Therefore, in general, we find that our sample is 80% mass complete above M ★ /M ⊙ ∼ 9.25 in all redshift bins, and therefore we plot SMFs above this conservative limit.We lose a negligible number of sources by limiting the sample in this manner (1 source each in the redshift bins 3 <  < 4 and 6 <  < 8).
To consider the completeness of our sample given the flux density limits of the telescope survey, we consider the widely used / max correction, used to test uniformity in the spatial distribution of sources (Schmidt (1968), see also Weaver et al. (2023b) for a detailed explanation) which particularly affects faint sources.This method considers the maximum redshift,  max , at which a source within a bin  low <  <  high would still be observable before falling below the detection limit.Each source is then associated with a maximum observable differential comoving volume,  max , associated with  max , and the actual differential comoving volume it is detected in, , associated with  high .If  max <  high , the source is given a weight of / max , and if  max >  high , / max = 1 (as the source would anyways have been detected in the survey, and therefore does not need to be given a higher weightage).Like the step-wise method used to calculate the SMF, the / max too is non-parametric.It assumes no functional form for the SMF, but it does assume a uniform spatial distribution of galaxies.However, Weaver et al. (2023b) show that this is problematic only at  < 1, thus not affecting our study.We apply the / max correction to our sources, finding that given the redshift bins we choose, no galaxies in our sample require this correction.This is expected, as our galaxies are red by definition and on average massive and therefore bright in F444W.The / max correction mostly affects only faint galaxies with the propensity to be detected close to the noise threshold.

Sources of uncertainty
We estimate the uncertainty of the SMFs by considering the Poisson noise  N , the uncertainty due to cosmic variance  cv , and the systematic uncertainty  sys due to SED-fitting.
Given that the calculation of the SMF is fundamentally a discrete counting process, the distribution of galaxies within a particular redshift and mass bin must follow Poissonian statistics.We calculate the uncertainty  N by using frequentist central confidence intervals 6(for details, see Maxwell 2011).
An added factor of uncertainty arises from cosmic variance, the field-to-field variation in galaxy number counts due to large-scale structure.It becomes an important source of uncertainty in narrow and deep surveys (Somerville et al. 2004), and is routinely included in uncertainty estimates of the stellar mass function (Davidzon et al. 2017;McLeod et al. 2021;Weaver et al. 2023b).To estimate the cosmic variance  cv , we use the CosmicVarianceCalculator v1.037 (Trenti & Stiavelli 2008), evaluated at the respective number density of our sample.We find relative cosmic variances for our sample to lie between 20%-30%, with the cosmic variance increasing with stellar mass.
Uncertainties on redshifts and stellar masses can give rise to a scatter,  fit , due to SED fitting.In order to estimate  fit , we generate 1000 independent realisations of the SMF by sampling from the posterior distributions of physical properties derived with BAGPIPES and calculate the variance of the number densities from these realisations.This method provides an estimate of the SMF as well as the uncertainty  fit on the SMF.2021) at  = 5 and  = 7 respectively (diamond and hexagon scatter points), derived from HST and Spitzer imaging.At 3 <  < 4, comparing our SMF to the pre- JWST Weaver et al. (2023b) SMF suggests that up to ∼ 30% of the galaxy population could have been missed at log M ★ /M ⊙ = 10.5 and up to ∼ 20% at log M ★ /M ⊙ = 11.0;similarly at 4 <  < 6, we find missed fractions of up to ∼ 25% at log M ★ /M ⊙ = 10.5 and 11.0.At 6 <  < 8, the obscured SMF exceeds the pre-JWST SMF from Weaver et al. (2023b) at log M ★ /M ⊙ = 10.375.At both 3 <  < 4 and 4 <  < 6, our SMFs dominate the dusty model SMF predicted by Long et al. (2023b) at log M ★ /M ⊙ > 9.5.
The final uncertainty  tot of the SMF is the quadrature addition of the Poisson uncertainty, cosmic variance, and the uncertainty due to SED fitting (as done in Davidzon et al. 2017), calculated via Equation 4: In the absence of detections, upper limits are calculated as the right confidence interval of the Poisson distribution.This is 1.841/(d ,comoving • Δ log M) following Gehrels (1986).

SMFs at 3 < z < 8
Figure 8 and Tables 2 and 3 present the SMFs of our sample in three redshift ranges: 3 <  < 4, 4 <  < 6, and 6 <  < 8, calculated using the method outlined in the previous section.We compare our dust obscured SMFs to the observed pre-JWST total SMFs from Weaver et al. (2023b), McLeod et al. (2021) and Stefanon et al. (2021), derived from ground-and space-based observations.We also compare our SMFs to model dust-obscured SMFs from Long et al. (2023b), which are derived from semi-empirical simulations of dusty star-forming galaxies (DSFGs).
In order to determine the previously missed fraction of the SMF, we assume that together, the selection functions of our study and pre-JWST studies produce a more complete survey than solely pre-JWST studies.Therefore, we compute upper limits on the previously missed fraction of the SMF by dividing our SMF by the sum of our SMF with the pre-JWST SMF from Weaver et al. (2023b).2021) (based on HST and groundbased observations).The 3 <  < 4 SMF deviates the most at the low-mass end but comes closest to the pre-JWST study at the highmass end, suggesting that up to ∼ 30% of the galaxy population could have been missed in the pre-JWST SMF from Weaver et al. (2023b) at log M ★ /M ⊙ ∼ 10.5 and up to ∼ 20% could have been missed at log M ★ /M ⊙ ∼ 11.0 -dusty galaxies detected with JWST therefore make up a sizeable fraction of the galaxy population at the highmass end, which suggests that galaxies at the high-mass end have been missing from our galaxy census in this epoch.Further, above log M ★ /M ⊙ ∼ 9.5, our SMF at 3 <  < 4 is significantly higher than the model SMF from Long et al. (2023b) (at  ∼ 3.0 − 3.5), with the difference being most pronounced at log M ★ /M ⊙ ∼ 10.5.Therefore, we could be seeing an emergent population of mainsequence dusty galaxies that are distinct from the widely studied DSFGs, which are typically more strongly star-forming (and which the Long et al. (2023b) simulation is based on).These results indicate that a significant population of obscured galaxies are prevalent at this redshift range.At 4 <  < 6 (central panel of Figure 8), we compare the SMF of our sample to Weaver et al. (2023b) (at  ∼ 4.5 − 5.5).For Table 2. Stellar mass function values of massive and dusty galaxies from 3 <  < 4 and 4 <  < 6, as shown graphically in the first two panels of Figure 8. Uncertainties are calculated as the quadrature addition of Poissonian noise, cosmic variance, and scatter due to SED fitting.2023b) SMF and exceeds it at log M ★ /M ⊙ = 10.375.This suggests the emergence of an extensive population of galaxies in the Epoch of Reionisation, hidden in the pre-JWST era but constituting a dominant part of the high-mass population.
We see a strong evolution in our SMFs around  ∼ 4 between masses of 9.5 ≲ log M ★ /M ⊙ ≲ 11.0.Comparing the SMFs at 4 <  < 6 and 3 <  < 4 in Table 2, we see an increase by a factor of ∼ 4 at log M ★ /M ⊙ ∼ 10.0, ∼ 9 at log M ★ /M ⊙ ∼ 10.5, and ∼ 7 at log M ★ /M ⊙ ∼ 11.0 -this shows an accelerated evolution in the knee of the SMF at this epoch, suggesting the onset of rapid dust-obscured stellar mass growth at  ∼ 4.
The SMFs of our sample between 4 <  < 6 and 6 <  < 8 show little evolution.At the high-mass end, we do not see a strong evolution across the whole redshift range, but we are heavily limited by small sample statistics, systematic uncertainties and cosmic variance, making it challenging to comment on SMF properties without a larger sample.
Globally, our analysis of red, dust-obscured galaxies shows that these sources recover a sizeable fraction of the high-mass end of the pre-JWST SMFs from Weaver et al. (2023b).Not only does this reveal the nature of the massive galaxy population, it highlights the efficiency of JWST in characterising the massive end of the galaxy SMF.

Integrated stellar mass density
The cosmic stellar mass density (SMD) is an efficient measure of stellar mass assembly.The total SMD is tightly coupled with the cosmic star-formation rate history, and thus could provide insights into early galaxy build-up such as previous epochs of star-formation and the stellar IMF of early stellar populations (Dickinson et al. 2003).Multiple works have observationally tracked the evolution of the SMD (Stark et al. 2009;González et al. 2010;Davidzon et al. 2017;McLeod et al. 2021;Weaver et al. 2023b), reaching up to  ∼ 8 − 10 (e.g., Weaver et al. 2023b;Stefanon et al. 2021).The observationally determined SMD, however, can be substantially affected if a significant population of high-mass galaxies have been missing in previous observations.This work in part aims to determine the fraction by which pre-JWST studies have underestimated the SMD.
We integrate the measured SMFs presented in Section 5.2 in order to get an estimate of the SMD for our galaxy sample.For each redshift bin '', we numerically integrate over the mass bins indexed by ' ' following Equation 5: where Φ ,  is the SMF value inferred via Equation 3, M  is the the central mass within each mass bin, Δ log M is the fixed mass bin size, and the limits are given by the mass range covered by our SMFs.Uncertainties on   are calculated via addition in quadrature, where the upper limits in the SMFs contribute to the upper uncertainty on   .We find that the SMD in units of [10 5 M ⊙ Mpc −3 ] is 51.6 24.8 −17.2 at 3 <  < 4, 7.5 +13.4 −3.0 at 4 <  < 6 and 4.6 +6.7 −2.4 at 6 <  < 8.The large uncertainty estimates reflect the uncertainty in the stellar mass functions where we are limited by sample size, especially at the high-mass end.Additionally, the upper limits in the highest mass bins make sizeable contributions to the upper uncertainties on the SMD.
In order to determine the missed SMD fraction in pre-JWST studies at the high-mass end, we compare our results with the Weaver et al. (2023b) study.As similarly explained in Section 5.2, we calculate an upper limit on the missed SMD fraction as the SMD of our sample divided by the sum of our SMD and the Weaver et al. (2023b) observed SMD, based on the assumption that our two studies together form a more complete survey than pre-JWST studies alone.
We integrate the observed SMFs from Weaver et al. (2023b) (shown in Figure 8) in order to estimate the observed pre-JWST total SMD.Given that the Weaver et al. (2023b) SMFs do not reach the lower mass limit of our study at 4 <  < 6 and 6 <  < 8, we expand the Weaver+23 SMFs to lower masses with their Schechter fits down to log M ★ /M ⊙ = 9.25, so as to perform a mass-consistent comparison with our sample.We find missed SMD fractions of 19 +9 −6 % at 3 <  < 4 and 15 +26 −6 % at 4 <  < 6.At 6 <  < 8, we find a missed fraction of 46 +66 −24 %, possibly doubling the SMD at this epoch.Therefore, our results indicate that the SMD could have been underestimated in pre-JWST studies, in particular significantly at  > 6.In future studies, it will be imperative to include dust-obscured galaxies at the high-mass end in order to accurately trace stellar mass build-up in the early universe.

DISCUSSION
In this section, we discuss the results of our work in the context of similar studies conducted with JWST's first year of observations on dusty galaxies.We additionally discuss the abundance of massive galaxies that is suggested by our dust-obscured SMFs, compare our SMD estimates with the SMDs estimated from integrating the Weaver et al. (2023b) Schechter functions, discuss the move towards redder selection functions, and place this in the context of past work and future studies on galaxy censuses.

Comparison of sample to recent literature in CEERS
JWST's pilot year has seen the output of a great amount of science, with several papers and teams already providing novel insights into obscured galaxies at  > 3 (e.g., Barrufet et al. 2023;Nelson et al. 2023;Pérez-González et al. 2023;Rodighiero et al. 2023;Labbé et al. 2023b;Akins et al. 2023).Additionally, it was shown that very dusty galaxies can sometimes contaminate extremely high redshift selections (e.g., Naidu et al. 2022;Zavala et al. 2023;Arrabal Haro et al. 2023).Here, we discuss our sample in comparison with some select studies in the CEERS field: Barrufet et al. (2023), Pérez-González et al. (2023), Labbé et al. (2023b), andNaidu et al. (2022).Barrufet et al. (2023) studied HST-dark galaxies in the CEERS field, identifying massive, obscured galaxies at  > 3 and into the Epoch of Reionisation.Of the 30 HST-dark sources in their study, we identify 12 in our sample, likely due to the different colour selection.Our SMF results support the findings of Barrufet et al. (2023) that suggest that a significant fraction of massive, obscured sources were previously missing from our galaxy census at  > 3.
Pérez- González et al. (2023) studied HST-dark and -faint galaxies in the first four NIRCam pointings of the CEERS field, using a selection based on F150W-F356W colours.Out of their sample of 138 HST-dark galaxies, we identify 65 sources in our sample.
Comparing their total sample to our study, we find similar redshift ranges (⟨⟩ = 3.68 +1.60  −1.00 in their study, ⟨⟩ = 3.46 +2.04 −1.35 in ours) and stellar masses (⟨log M ★ /M ⊙ ⟩ = 10.20 +0.46  −0.73 in their study, ⟨log M ★ /M ⊙ ⟩ = 10.15 +0.43  −0.50 in ours).We note however that our redshift distribution has a longer high-end tail, where we find more sources at  ≳ 6 than the Pérez- González et al. (2023) study.This is most likely because we use the longer wavelength F444W filter in our colour selection, where we are possibly picking up the [OIII] line at  ∼ 7.
Using a selection based on blue rest-UV and red restoptical colours, Labbé et al. (2023b) found six massive galaxies (M ★ /M ⊙ > 10 10 ) at 7.4 <  < 9.1.We identify two of their sources in our sample (IDs 48444 and 67066).We most likely do not select the remaining four sources in Labbé et al. (2023b) due to their blue rest-UV colour selection.Additionally, one of the Labbé et al. (2023b) sources originally identified as a massive galaxy at  = 8.13 has now been spectroscopically determined to be a likely AGN candidate at  = 5.64 (Kocevski et al. 2023); we do not find this source in our sample.Naidu et al. (2022) proposed a luminous candidate  ≈ 17 or  ≈ 5 galaxy, dubbed "Schrodinger's Galaxy", now confirmed to be an obscured source at  = 4.912 ± 0.001 (Arrabal Haro et al. 2023).We find this galaxy in our sample (ID 81918) at  = 4.79 +0.05 −0.08 with a dust attenuation of A V = 1.74 +0.11  −0.17 mag.Such studies show that there is increasing evidence for a population of massive, obscured galaxies at high redshifts, close to and into the Epoch of Reionisation (see also Fudamoto et al. 2021).

Abundance of red galaxies at the high-mass end of SMFs
The SMFs of JWST-detected dust-obscured galaxies in our study point toward an abundance of galaxies at the massive end of the pre-JWST SMF, possibly leading to an excess of the galaxy population with respect to the pre-JWST determined SMF.This abundance is even more pronounced with respect to the Schechter fits from Weaver et al. (2023b) (solid lines in Figure 8).Comparing our SMD estimates to those found by integrating the Schechter fits from Weaver et al. (2023b) down to log M ★ /M ⊙ = 9.25, we find missed SMD fractions of 19 +9 −6 % at 3 <  < 4 and 18 +32 −7 % at 4 <  < 6.At 6 <  < 8, we find a missed SMD fraction of 52 +76 −27 %, effectively doubling the SMD at this epoch.This excess with respect to the Schechter fit at the massive end of the SMFs at  ∼ 3 − 5 was shown in Weaver et al. (2023b) with a sample of 2m-selected sources from the COSMOS2020 dataset, with hints that this population could be star-forming, dusty galaxies.
It is evident from past work that selecting sources deeper into the NIR results in stronger constraints on the high-mass end of SMFs.At its time, the Weaver et al. (2023b) study represented some of the reddest SMFs in comparison with earlier studies (e.g.Davidzon et al. 2017;Stefanon et al. 2017bStefanon et al. , 2021)).The effect of this is evident from Figure 8 2021) SMF at 6 <  < 8 at the high mass end.Now, with JWST/NIRCam allowing us to move even deeper into the NIR regime, our study represents the natural next step in the move towards more complete selections: with sources selected based on their 1.5m-4.44mcolour, our study enables a more complete characterisation of massive and dusty galaxies than was possible with previous SMF studies.
Our results reinforce the conclusion that dust-obscured galaxies contribute significantly to the high-mass end of the SMF.In future galaxy censuses with JWST, it will be critical to explore how the SMF measurements at the high-mass end compare with the Schechter formalism of our description of galaxy evolution.

SUMMARY AND CONCLUSION
In this work, we used data from the JWST/CEERS survey (Finkelstein et al. 2022(Finkelstein et al. , 2023) ) in the CANDELS/EGS field to identify red, optically-faint galaxies at high redshifts in order to determine the obscured stellar mass function at various epochs in the first two billion years of the history of the universe.Some key results are summarised in the following: • Using a colour criterion designed to select red, optically-faint galaxies, we show that we efficiently select massive and dusty galaxies (⟨log M ★ /M ⊙ ⟩ = 10.15 +0.43  −0.50 and ⟨A V ⟩ = 2.71 +0.88 −0.91 mag) with a majority lying at  > 3 (see Figures 3 and 4).
• Our sample contains predominantly star-forming galaxies, largely lying on the star-forming main-sequence.They therefore represent a normal population of galaxies without extreme starburst properties (see Figures 5 and 6).Our sample overlaps with the Wang et al. (2019) sample at the high-mass end and the Barrufet et al. (2023) sample at the low-mass end, showing that our sample of red galaxies has similar star-forming properties to that of HST-dark galaxies (see Figure 7).
• Our analysis of the obscured galaxy SMF (see Figure 8) shows that in the pre-JWST era, we have missed a significant fraction of galaxies, particularly at the high-mass end of the SMF at redshifts of  > 3. The SMFs of red, optically-faint galaxies suggest a missed fraction of ≳ 20% of the galaxy population in the 3 <  < 4 and 4 <  < 6 epochs (at log  ★ / ⊙ ≥ 10.5).At 6 <  < 8, our SMF overtakes the pre-JWST SMF from Weaver et al. (2023b) around log M ★ /M ⊙ ∼ 10.375.
• Our results at 6 <  < 8 highlight the importance of accounting for massive, dust-obscured galaxies in the final stages of the Epoch of Reionisation.
• Our SMFs show a strong evolution at  ∼ 4 at masses of 9.5 ≲ log M ★ /M ⊙ ≲ 11.0, suggesting the onset of rapid dustobscured stellar mass assembly in this epoch.
• The derived stellar mass density of our sources at log M ★ /M ⊙ ≥ 9.25 suggests that the missed SMD fraction could be a factor of ∼15-20% at  ∼ 3 − 6.We find a missed fraction of ∼45% at  ∼ 6 − 8, possibly doubling the SMD at this epoch.
These findings point towards an emergent population of massive, obscured galaxies from  ∼ 3 up to and into the Epoch of Reionisation, supporting the findings of early JWST studies (e.g., Barrufet et al. 2023;Labbé et al. 2023b;Akins et al. 2023).The strong evolution of the SMF at  ∼ 4 suggests that this is a period of rapid stellar mass growth in obscured galaxies.Interestingly,  ∼ 4 is also roughly when the obscured SFRD is thought to overtake the un-obscured SFRD, dominating the cosmic star-formation history at later epochs (e.g., Zavala et al. 2021;Bouwens et al. 2020Bouwens et al. , 2021)).
Our results indicate that obscured stellar mass assembly occurred as early as  ∼ 8, suggesting that the build-up of dusty galaxies could begin close to 600 Myrs after the Big Bang.To further explore the beginning of obscured stellar mass assembly and push the observable redshift boundary farther back, studying the SMF by collating all public JWST surveys is critical.Including surveys such as COSMOS-Web (Casey et al. 2023), PRIMER (Dunlop et al. 2021), UNCOVER (Bezanson et al. 2022) and PANORAMIC (Williams et al. 2021) will satisfy the need of the hour: larger sample sizes.These surveys, and others to come with JWST, will surely result in us establishing a complete census of the massive, dust-obscured galaxy population in the early universe.
support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51515.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.MS acknowledges support from the CIDEGENT/2021/059 grant, from project PID2019-109592GB-I00/AEI/10.13039/501100011033 from the Spanish Ministerio de Ciencia e Innovación -Agencia Estatal de Investigación.MS also acknowledges the financial support from the MCIN with funding from the European Union NextGenerationEU and Generalitat Valenciana in the call Programa de Planes Complementarios de I+D+i (PRTR 2022) Project (VAL-JPAS), reference ASFAE/2022/025.The work of CCW is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

APPENDIX A: AGN IDENTIFICATION AND EFFECT ON THE SMF
Given that our study focuses on star-forming galaxies, it is of importance to remove AGN from our sample.We identify and remove AGN candidates, the so-called little red dots (LRDs) as described in Section 2.4. Figure A1 shows the postage stamps and SED of one such AGN candidate, galaxy 6583.This is a very compact source, as is characteristic of LRDs, with a red slope beyond 2m and a blue slope below this.As shown, BAGPIPES does not fit the short-wavelength end of the slope well, possibly because BAGPIPES cannot perform multi-component SED-fitting and additionally does not contain AGN templates.This results in inaccurate photometric redshifts and derived physical properties of LRDs.Further, Figure A2 shows the effect of LRDs on the stellar mass function of our sample.While the SMF of the full sample overlaps neatly with the AGN-cleaned sample at 3 <  < 4, the difference between SMFs is more pronounced at 4 <  < 6 and differs the most at 6 <  < 8.This highlights the importance of addressing the presence of LRDs in our sample, so as not to overestimate the SMFs and stellar mass density at high-redshifts.The AGN-removal has the least effect on the SMF at 3 <  < 4, a larger effect at 4 <  < 6 and the most pronounced effect at 6 <  < 8.This highlights the importance of identifying and removing AGN candidates at high redshifts in order to derive accurate mass functions and SMD estimates.

Figure 1 .
Figure 1.Colour-magnitude diagram of F150W-F444W vs. F444W showing our selection method.The grey scatter points show all the sources in the CEERS catalogue while the orange scatter points are our selected galaxies.The coloured lines are SED tracks for various dust attenuation values (A V is indicated in boxes), with coloured numbers indicating the redshift -solid lines correspond to 10 10 M ⊙ and the dashed line corresponds to 10 11 M ⊙ .The grey arrows show upper limits for the sources with F150W mag lower than 2 (median errors are shown by the cross in the upper right of the figure).The colour criterion F150W-F444W>2.1 mag (black dashed line) in principle identifies  > 3 sources with A V ≳ 2 mag and log(M ★ /M ⊙ ) ∼ 10, while the magnitude cut F150W>25 mag (dark red dashed line) is designed to rid the sample of low- sources while retaining the most massive galaxies in the sample (∼ 10 11 M ⊙ ).Also shown for reference is the F150W = 26 mag cut that is a proxy for identifying HST-dark sources (Pérez-González et al. 2023).We select 179 red galaxies with these criteria that theoretically restrict our sample to massive, dusty, high-redshift galaxies.
Figure2.Postage stamps and SED fits of four selected galaxies from our sample of red galaxies.The stamps boxed in blue are from ancillary HST/ACS and WFC3 data, and the stamps boxed in red are from JWST/NIRCam imaging (each stamp is 4 × 4 arcsec 2 ).There is a variety in the morphological properties of our sample, ranging from spatially extended sources to compact ones.The lower panels display the SED fits: the maroon points represent the photometry and the downward arrows represent the flux upper limits.The orange lines are the SED fits from BAGPIPES and the photometric redshift probability density functions are inlaid in the lower right part of the graphs.The physical properties of these galaxies are quoted on the graphs.They are massive (log M ★ /M ⊙ ≳ 9.5) and dusty (A V ∼ 1.5 − 4 mag) with redshifts ranging from  ∼ 3 − 8.

Figure 6 .
Figure 6.Star formation rate vs. stellar mass for our sample of red, opticallyfaint galaxies, coloured by photo- (circular scatter points).Uncertainties are given by the 16 th and 84 th percentile of the posterior distribution from SED-fitting with BAGPIPES.The galaxy main-sequence lines shown at  = 2, 4, and 6 are from Speagle et al. (2014) (solid coloured lines with scatter).The majority of our sample lies on the MS at redshifts of  < 6, suggesting that they are normal star-forming galaxies with moderate SFRs.The three sources lying significantly below the MS are candidate quiescent galaxies.

Figure 7 .
Figure7.Star formation rate vs. stellar mass for the subset of our sample of red, optically-faint galaxies at 3 <  < 5 (circular scatter points).We show the galaxy main-sequence at  = 4. fromSpeagle et al. (2014) (solid line with scatter).We compare our sample to HST-dark galaxies fromWang et al. (2019) with  median = 4 (empty diamonds), and a sample subset fromBarrufet et al. (2023) at 3 <  < 5 (red squares).Our sample overlaps with theBarrufet et al. (2023) sample at the lower-mass end and with theWang et al. (2019) sample at the high-mass end, showing that our study covers the mass-range spanned by HST-dark/faint galaxies in both the pre-JWST and JWST era.

Figure 8 .
Figure 8. SMFs of our sample of massive and dusty galaxies in three redshift ranges: 3 <  < 4, 4 <  < 6, and 6 <  < 8. Uncertainties shown are derived from Poisson statistics, cosmic variance, and scatter due to SED fitting.Upper bounds (downward arrows) are derived from the right confidence interval of the Poisson distribution.The fixed size of the mass bins are shown in the lower left of each panel.We compare our results to the observed SMFs of the pre-JWST total galaxy population from Weaver et al. (2023b) (blue scatter points and shaded area) derived from COSMOS2020 observations (Weaver et al. 2023a) and to the model SMFs from Long et al. (2023b) (solid grey line) derived from semi-empirical simulations for dusty star-forming galaxies.For reference, the Schechter fits from Weaver et al. (2023b) are also shown (solid blue line).The SMF at 3 <  < 4 is additionally compared to McLeod et al. (2021) (square scatter points), derived from ground-based observations, and the 4 <  < 6 and 6 <  < 8 SMFs are compared to Stefanon et al. (2017b) and Stefanon et al. (2021) at  = 5 and  = 7 respectively (diamond and hexagon scatter points), derived from HST and Spitzer imaging.At 3 <  < 4, comparing our SMF to the pre-JWST Weaver  et al. (2023b)  SMF suggests that up to ∼ 30% of the galaxy population could have been missed at log M ★ /M ⊙ = 10.5 and up to ∼ 20% at log M ★ /M ⊙ = 11.0;similarly at 4 <  < 6, we find missed fractions of up to ∼ 25% at log M ★ /M ⊙ = 10.5 and 11.0.At 6 <  < 8, the obscured SMF exceeds the pre-JWST SMF fromWeaver et al. (2023b) at log M ★ /M ⊙ = 10.375.At both 3 <  < 4 and 4 <  < 6, our SMFs dominate the dusty model SMF predicted byLong et al. (2023b) at log M ★ /M ⊙ > 9.5.
The left panel of Figure8shows our SMF at 3 <  < 4 in comparison with Weaver et al. (2023b) (at  ∼ 3.0−3.5)and McLeod et al. (2021) (at  ∼ 3.25).At all masses shown in this redshift range, the SMF of our sample lies below the pre-JWST SMF from Weaver et al. (2023b) (based on COSMOS2020 observations; see Weaver et al. 2023a) and McLeod et al. ( reference, theStefanon et al. (2017b) mass function for LBGs at  = 5 is shown.Of particular interest in this epoch is the comparison of our sample to the pre-JWST SMF at the high-mass end, where at log M ★ /M ⊙ ∼ 10.5 and 11.0, we find that up to ∼ 25% of the SMF could have been missed in theWeaver et al. (2023b) mass function.In addition, like at 3 <  < 4, we find an SMF at 4 <  < 6 which is significantly higher than the dust-obscured model SMF predicted by theLong et al. (2023b) simulation at  ∼ 4.0 − 6.0.At 6 <  < 8 (right panel of Figure8), we compare our SMF toWeaver et al. (2023b) (at  ∼ 6.5 − 7.5).For reference, theStefanon et al. (2021) mass function at  = 7 is shown.Above log M ★ /M ⊙ ∼ 10.0, our sample overtakes theWeaver et al. ( , where the Weaver et al. (2023b) SMF overtakes the LBG-based Stefanon et al. (2017b) SMF at 4 <  < 6 and the Stefanon et al. (

Figure A1 .
Figure A1.Postage stamps and SED of galaxy 6583, identified as a potential AGN and removed from our final sample.The postage stamps show the compactness of the source.The SED shows a characteristic red slope above 2m, and a blue slope below 2m.The blue part of the slope is poorly fit with BAGPIPES, thus resulting in an inaccurate photometric redshift and subsequently inaccurate derived physical properties.

Figure A2 .
Figure A2.SMFs shown as in Figure8.In addition, open circles represent the SMFs of the full sample of 168 sources before AGN-removal.The AGN-removal has the least effect on the SMF at 3 <  < 4, a larger effect at 4 <  < 6 and the most pronounced effect at 6 <  < 8.This highlights the importance of identifying and removing AGN candidates at high redshifts in order to derive accurate mass functions and SMD estimates.

Table 3 .
Stellar mass function values of massive and dusty galaxies from 6 <  < 8, as shown graphically in the third panel of Figure8.Uncertainties are calculated as the quadrature addition of Poissonian noise, cosmic variance, and scatter due to SED-fitting.

Table B1 -
continued S.No.ID RA Dec  phot log M ★ /M ⊙ log SFR/M ⊙ yr −1 A V /mag