We analyze the evolution of massive logSUB10SUBM SUBSUBM SUBSUB gt 10 galaxies at z  14 selected from JWST Cosmic Evolution Early Release Survey CEERS. We infer the physical properties of all galaxies in the CEERS NIRCam imaging through spectral energy distribution SED fitting with dense basis to select a sample of highredshift massive galaxies. Where available we include constraints from additional CEERS observing modes including 18 sources with MIRI photometric coverage and 28 sources with spectroscopic confirmations from NIRSpec or NIRCam WFSS. We sample the recovered posteriors in stellar mass from SED fitting to infer the volume densities of massive galaxies across cosmic time taking into consideration the potential for sample contamination by active galactic nuclei. We find that the evolving abundance of massive galaxies tracks expectations based on a constant baryon conversion efficiency in dark matter halos for z  14. At higher redshifts we observe an excess abundance of massive galaxies relative to this simple model resulting in a shallower decline of observed volume densities of massive galaxies. These higher abundances can be explained by modest changes to star formation physics andor the efficiencies with which star formation occurs in massive dark matter halos and are not in tension with modern cosmology.
['Galaxy evolution', 'Galaxy formation', 'High-redshift galaxies', 'Galactic and extragalactic astronomy', '594', '595', '734', '563', 'Astrophysics - Astrophysics of Galaxies']
Evidence for a Shallow Evolution in the Volume Densities of Massive Galaxies at z  48 from CEERS
{'Evidence for a Shallow Evolution in the Volume Densities of Massive Galaxies at z = 4 to 8 from CEERS': "Katherine Chworowsky, 1, ∗ Steven L. Finkelstein, 1 Michael Boylan-Kolchin, 1 Elizabeth J. McGrath, 2 Kartheik G. Iyer, 3 Casey Papovich, 4, 5 Mark Dickinson, 6 Anthony J. Taylor, 1 L. Y. Aaron Yung, 7, 8, † Pablo Arrabal Haro, 6 Micaela B. Bagley, 1 Bren E. Backhaus, 9 Rachana Bhatawdekar, 10 Yingjie Cheng, 11 Nikko J. Cleri, 4, 5 Justin W. Cole, 4, 5 M. C. Cooper, 12 Luca Costantin, 13 Avishai Dekel, 14 Maximilien Franco, 15 Seiji Fujimoto, 16, 17 Christopher C. Hayward, 18 Benne W. Holwerda, 19 Marc Huertas-Company, 20, 21, 22 Michaela Hirschmann, 23 Taylor A. Hutchison, 7, † Anton M. Koekemoer, 8 Rebecca L. Larson, 24 Zhaozhou Li, 14 Arianna S. Long, 15, ‡ Ray A. Lucas, 25 Nor Pirzkal, 26 Giulia Rodighiero, 27, 28 Rachel S. Somerville, 29 Brittany N. Vanderhoof, 30 Alexander de la Vega, 31 Stephen M. Wilkins, 32, 33 Guang Yang, 34, 35 and Jorge A. Zavala 36  \n1 Department of Astronomy, The University of Texas at Austin, Austin, TX 78712 USA  \n2 Department of Physics and Astronomy, Colby College, Waterville, ME 04901, USA 3 Dunlap Institute for Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada 4 Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA 5 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA 6 NSF's National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA 7 Astrophysics Science Division, NASA Goddard Space Flight Center, 8800 Greenbelt Rd, Greenbelt, MD 20771, USA 8 Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA 9 Department of Physics, 196 Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269 10 European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Ca˜nada, Madrid, Spain 11 University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA 12 Department of Physics & Astronomy, University of California, Irvine, 4129 Reines Hall, Irvine, CA 92697, USA 13 Centro de Astrobiolog'ıa (CAB), CSIC-INTA, Ctra de Ajalvir km 4, Torrej'on de Ardoz, 28850, Madrid, Spain 14 Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel 15 Department of Astronomy, The University of Texas at Austin, Austin, TX, USA 16 Cosmic Dawn Center (DAWN), Jagtvej 128, DK2200 Copenhagen N, Denmark 17 Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK2100 Copenhagen Ø, Denmark 18 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA 19 Physics & Astronomy Department, University of Louisville, 40292 KY, Louisville, USA 20 Instituto de Astrof'ısica de Canarias, La Laguna, Tenerife, Spain 21 Universidad de la Laguna, La Laguna, Tenerife, Spain 22 Universit'e Paris-Cit'e, LERMA - Observatoire de Paris, PSL, Paris, France 23 Institute of Physics, Laboratory of Galaxy Evolution, Ecole Polytechnique F'ed'erale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland 24 School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA 25 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA 26 ESA/AURA Space Telescope Science Institute 27 Department of Physics and Astronomy, Universit'a degli Studi di Padova, Vicolo dell'Osservatorio 3, I-35122, Padova, Italy 28 INAF - Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122, Padova, Italy 29 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY, 10010, USA 30 Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA 31 Department of Physics and Astronomy, University of California, 900 University Ave, Riverside, CA 92521, USA 32 Astronomy Centre, University of Sussex, Falmer, Brighton BN1 9QH, UK 33 Institute of Space Sciences and Astronomy, University of Malta, Msida MSD 2080, Malta 34 Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands 35 SRON Netherlands Institute for Space Research, Postbus 800, 9700 AV Groningen, The Netherlands 36 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan", 'ABSTRACT': 'We analyze the evolution of massive (log 10 [ M ⋆ /M ⊙ ] > 10) galaxies at z ∼ 4-8 selected from the JWST Cosmic Evolution Early Release Science (CEERS) survey. We infer the physical properties of all galaxies in the CEERS NIRCam imaging through spectral energy distribution (SED) fitting with dense basis to select a sample of high redshift massive galaxies. Where available we include constraints from additional CEERS observing modes, including 18 sources with MIRI photometric coverage, and 28 sources with spectroscopic confirmations from NIRSpec or NIRCam wide-field slitless spectroscopy. We sample the recovered posteriors in stellar mass from SED fitting to infer the volume densities of massive galaxies across cosmic time, taking into consideration the potential for sample contamination by active galactic nuclei (AGN). We find that the evolving abundance of massive galaxies tracks expectations based on a constant baryon conversion efficiency in dark matter halos for z ∼ 1-4. At higher redshifts, we observe an excess abundance of massive galaxies relative to this simple model. These higher abundances can be explained by modest changes to star formation physics and/or the efficiencies with which star formation occurs in massive dark matter halos, and are not in tension with modern cosmology.  \nKeywords: Galaxy Evolution, Extragalactic Astronomy, Galaxy Formation', '1. INTRODUCTION': "Measurements of galaxy stellar masses can provide a powerful benchmark to compare to theoretical models and simulations. The ability of the Universe to create and form galaxies with large amounts of stellar mass is largely dependent on galaxy baryonic feedback processes, both internally and within a galaxy's environment. In particular, massive galaxies at high redshift provide interesting laboratories to study galaxy formation physics, as these extreme systems allow us to begin to understand how the Universe is able to rapidly build up large amounts of stellar mass in a short amount of time, and to constrain the feedback mechanisms that may be regulating these processes (e.g. Hayward et al. 2021).  \nObservationally, extensive efforts have been made in attempting to understand how galaxy stellar masses evolve over time, with many studies measuring galaxy stellar mass functions (GSMFs) spanning cosmic time (e.g. Moffett et al. 2016; Weaver et al. 2022a; Adams et al. 2021; McLeod et al. 2021; Duncan et al. 2014; Muzzin et al. 2013; Song et al. 2016; Bhatawdekar et al. 2019; Tomczak et al. 2014; Wright et al. 2017; Kelvin et al. 2014; Stefanon et al. 2021; Leja et al. 2020). The stellar mass of a galaxy is typically inferred via modeling the rest-frame ultraviolet (UV) through near-infrared (NIR) spectral energy distribution (SED), which necessarily relies on assumptions about the galaxy's star  \nformation history, metal enrichment history, dust attenuation, and stellar initial mass function (IMF). Prior to JWST , the Hubble Space Telescope ( HST ) provided the most accurate galaxy stellar mass estimates out to redshifts of z ∼ 5, however, at increasingly high redshift, the wavelength coverage of HST probes an increasingly limited wavelength range of the SED, and only UV-luminous galaxies are detected. While the Spitzer Infrared Array Camera (IRAC) provided longer wavelength space-based observations (e.g. Song et al. 2016; Bhatawdekar et al. 2019; Morishita 2021; Tacchella et al. 2022), the poor spatial resolution ( ∼ 2') and small collecting area of Spitzer IRAC with respect to HST meant that galaxy observations were often blended with neighboring sources, limiting the ability of observstions to accurately attribute observed fluxes to their sources.  \nThe launch of JWST (Gardner et al. 2023) enabled a huge leap forward in both the resolution and sensitivity of infrared imaging. With the longer wavelength coverage of JWST Near Infrared Camera (NIRcam; Rieke et al. 2023), observations of the rest-frame optical (including the star-formation history sensitive 4000 ˚ A break) are accessible out to z ∼ 10, allowing for the measurement of accurate stellar masses out to redshifts that were previously inaccessible. Nearly two years since the launch of JWST , data collected has already begun to revolutionize our understanding of the Universe, with new discoveries answering long-standing questions in the field of extra-galactic science, while also providing unexpected results and generating new questions.  \nOne of these unexpected early results is an apparent excess of UV-luminous z ≳ 8 galaxies relative to many recent theoretical models (e.g. Castellano et al. 2022;  \nNaidu et al. 2022; Finkelstein et al. 2022, 2023a,b; Casey et al. 2023; Franco et al. 2023; Leung et al. 2023; Adams et al. 2023; Bouwens et al. 2023; Harikane et al. 2023). This population of galaxies may indicate evidence for different star formation and stellar feedback physics in the early Universe, leading to earlier, more rapid formation of massive and luminous galaxies than previously thought possible (e.g. Mason et al. 2023; Shen et al. 2023; Pallottini & Ferrara 2023; Yung et al. 2023a; Dekel et al. 2023; Ferrara et al. 2023). While models shown tension with direct observables, these UV luminosities are very sensitive to assumed initial mass functions, dust, and AGN contamination. Having reliable estimates for observed stellar masses provides an alternative way of testing this putative tension with models. Excitingly, early JWST studies have claimed evidence for very massive galaxies (some exceeding log 10 ( M ⋆ /M ⊙ ) > 11), further indicating the possibility of more efficient starformation activity at higher redshift than expected (e.g. Labb'e et al. 2023; Xiao et al. 2023), and/or a revision to ΛCDM (Boylan-Kolchin 2023).  \nTo improve our understanding of the abundance of massive (defined here as log 10 [ M ⋆ /M ⊙ ] > 10) galaxies and how this population evolves over cosmic time, we present a study of massive galaxies selected over the full the 88 arcmin 2 Comic Evolution Early Release Survey (CEERS) NIRCam field. We briefly introduce CEERS and the data reduction ( § 2) before discussing how the galaxies studied here were selected ( § 3), focusing on sources at z > 4. In selecting massive galaxies, we also note the interesting prominence of a population of galaxies arising from recent JWST observations: point-like sources with flat or blue rest UV-optical spectra alongside increasing brightness at longer wavelength, possibly hosting dust-reddened active galactic nuclei (AGN Barro et al. 2023; Kocevski et al. 2023; Greene et al. 2023; Labbe et al. 2023). We discuss how this category of sources may bias selection of massive galaxies ( § 3.4). We present how the cumulative number density of massive galaxies evolves across cosmic time in Section 4 and compare our findings to both previous studies ( § 5.1) and model predictions ( § 5.2). Finally, we discuss the our observed abundances of galaxies may indicate changes to the efficiency of star-formation in galaxies ( § 5.3.1) or how the observed light from the galaxy is translated to mass ( § 5.3.2). Throughout this paper, we assume a cosmology of H 0 =70 km s -1 Mpc -1 , Ω M =0.3, Ω Λ = 0.7. When necessary, we assume a baryonic density parameter of Ω b h 2 = 0 . 02237 (Planck Collaboration et al. 2016). We use a Chabrier (2003) IMF, and all magnitudes given are in the AB system (Oke & Gunn 1983).", '2. NIRCAM & HST CATALOG': "We perform our analysis on the publicly released CEERS NIRCam mosaics (Data Releases 0.5 and 0.6), which include imaging from 1 -5 µ m in F115W, F150W, F200W, F277W, F356W, F410M and F444W and cover ten pointings in the Extended Groth Strip (EGS) (Bagley et al. 2023). The CEERS observations were completed and reduced over two epochs. Epoch 1 NIRCam pointings 1, 2, 3 and 6 were obtained in parallel to the Mid-Infrared Instrument (MIRI) imaging in June 2022. These pointings are reduced with JWST Calibration Pipeline (Bushouse et al. 2023) version 1.7.2 and Calibration Reference Data System (CRDS) pmap 0989, and are available from CEERS DR0.5. The Epoch 2 observations released with CEERS DR0.6 include NIRCam imaging over pointings 4, 5, 7, 8, 9 and 10 obtained in December 2022 in parallel to NIRSpec MSA observations. Pointings 5, 7, 8, 9 were also observed with NIRCam WFSS with MIRI imaging in parallel. These Epoch 2 images were reduced with Pipeline version 1.8.5 and pmap 1023. The two CEERS NIRCam reductions follow identical procedures, with the only difference being in Pipeline and CRDS versions. We note that the only NIRCam update between pmaps 1089 and 1023 is to the distortion reference files, and so this change will not affect the noise properties, flat fields, or flux calibration between CEERS DR0.5 and DR0.6. The CEERS NIRCam imaging reduction uses the JWST Calibration Pipeline, with custom modifications to correct for features and challenges in the data such as snowballs, wisps, 1/f noise and to improve the astrometric alignment. We direct the reader to Bagley et al. (2023) for further details on the reduction.  \nWe make use of the photometric catalog from Finkelstein et al. (2023b), which also includes the available HST imaging in this field from CEERS HST Data Release 1, including the ACS F606W, F814W and WFC3 F105W, F125W, F140W and F160W filters from HST programs including CANDELS (Grogin et al. 2011; Koekemoer et al. 2011). This catalog's methodology is similar to Finkelstein et al. (2023a), with some updates to improve the accuracy of color and total flux measurements. Flux uncertainties were calculated empirically via measurements of the flux spread in randomly placed apertures. We direct the reader to Finkelstein et al. (2023b) for full details on all cataloguing procedures.", '3. SELECTION OF MASSIVE GALAXIES': "Here, we detail our process of selecting a sample of massive galaxies. We first make use of the eazy (Brammer et al. 2008) to determine photometric redshifts. We then use the derived redshift probability distribution  \nFigure 1. The distribution in redshift and stellar mass of the full CEERS catalog (Finkelstein et al. 2023b), as measured by dense basis (with eazy priors on the dense-basis -derived redshift). The darker shaded region denotes our redshift and mass range of interest, while the dark points indicate any source which has been visually inspected, the larger of which are sources whose median stellar mass and redshift satisfy our high redshift massive criteria. We also show galaxies with MIRI coverage (red triangles) and those with spectroscopic confirmation (orange triangles). The galaxies marked with a red box are 'Extremely Red Objects' (EROs): F444W point sources with red F277W-F444W colors that are likely hosts to dust-reddened AGN (discussion in § 3.4). We find 118 galaxies with median stellar mass and redshift values in our region of interest, including objects up to z ∼ 7.5. Notably, we find no sources at log 10 ( M ⋆ /M ⊙ ) > 11 or at z > 8, and the highest redshift sources we do find are all classified as EROs, thus reliable mass measurements are not possible.  \n<!-- image -->  \nfunctions (PDFs, or P ( z )) as a prior to dense basis (Iyer & Gawiser 2017; Iyer et al. 2019), which we use to perform Bayesian SED fitting to recover physical properties (i.e. stellar mass) and their corresponding posterior distributions.", '3.1. Photometric Redshift Fitting': "Photometric redshifts for the catalog were determined with eazy . We use these results from Finkelstein et al. (2023b), briefly discussing the procedure here. eazy fits non-negative linear combinations of stellar population synthesis templates through a user-defined grid of redshifts, and at each redshift finds the best fitting synthetic template spectrum by minimizing χ 2 . eazy derives a P ( z ) based on the quality of fit to the observed photometry. The templates used include the 'tweak fsps QSF 12 v3' set of 12 FSPS (Conroy & Gunn 2010) templates recommended by the eazy documentation as well as an additional six templates designed to encompass bluer rest-frame UV colors expected from z > 9 galaxies (Larson et al. 2022).", '3.2. SED fitting': "We perform spectral energy distribution (SED) fitting for all sources with signal-to-noise in F444W > 3 with  \nthe dense basis SED fitting code, which is designed to robustly recover star formation history (SFH) constraints from galaxy SEDs. dense basis uses a flexible non-parametric SFH represented by a Gaussian Mixture Model (GMM; Iyer et al. 2019), and stellar templates generated from FSPS (Conroy et al. 2009; Conroy & Gunn 2010) including implementation of nebular emission lines using CLOUDY (Ferland et al. 2017; Byler et al. 2017). For this work, we define three 'shape' parameters that describe the SFH: t 25 , t 50 , and t 75 (requiring the recovered SFH of the galaxy to form ' x ' fraction of its total mass by time t x ). All sources were fit with dense basis assuming a Calzetti et al. (2000) dust law and a Chabrier initial mass function (IMF Chabrier 2003). We impose a uniform (flat) prior on the specific star formation rate (sSFR) with limits on the sSFR (sSFRyr -1 ∈ [ -14 , -7]), an exponential prior on the dust attenuation over a wide range of values ( A V ∈ [0 , 4]), and a uniform (in log-space) prior on the metallicity ( Z/Z ⊙ ∈ [0 . 01 , 2 . 0]). For redshift, dense basis is only able to take a prior in the form of a top-hat function. Therefore, in order to input the photometric redshift information as recovered by eazy , we modify the eazy P ( z ) to a top-hat form by first fitting a Gaussian to the primary peak of the P ( z ), and taking the width  \nof the tophat function as the 2 σ width of the Gaussian, if σ > 1 . 5, or the 5 σ width of the Gaussian, if σ ≤ 1 . 5 (in the case of a very narrowly peaked Gaussian). For each source, we thus have a redshift prior based on the eazy recovered redshift PDF.  \nWhile we perform SED fitting over the entire catalog and all available redshifts, as the emphasis of this work is on the evolution at z > 4, we further vet that higherredshift sample to ensure accurate results. Therefore, we split our massive sample into two categories: 1. a low redshift sample of 1796 sources with best fit redshift z < 3 . 5 and recovered median stellar masses of log 10 ( M ⋆ /M ⊙ ) > 10, and 2. a parent high redshift sample, with a photometric redshift of z > 3 . 5 (from dense basis ), at least 70% of the integrated eazy P ( z ) greater than 2.5, and a less stringent mass cut of 2.5% of the posterior in stellar mass having log 10 ( M ⋆ /M ⊙ ) > 10. The parent high redshift sample is intentionally chosen to include sources which have median stellar masses below our intended mass limit, to ensure we perform visual inspection on all sources where posteriors may include a massive (log 10 ( M ⋆ /M ⊙ ) > 10) solution to assess the reliability of the photometry. We thus manually inspect 561 sources satisfying these criteria.  \nWe remove a total of 297 sources from our sample due to unreliable photometry (e.g. too close to a detector edge, bad pixels, spurious sources from Source Extractor (Bertin & Arnouts 2010), and too close to a neighboring bright source), the majority (168/297) sources removed have z phot > 8, therefore we are confident in their spurious nature. Two particularly high redshift sources were determined to be incorrectly fit at high redshift (IDs 52092 and 75728 at z phot ∼ 15 and z phot ∼ 12, respectively). Source ID 52092 has a lowz solution with ∆ χ 2 < 4 from the highz fit, and source ID 75728 had significant flux detected below the putative Lymanα break at its fit redshift. Both sources were thus refit with dense basis at the low redshift solution from eazy . After refitting, one source was removed from the sample for not satisfying our mass limit.  \n53 additional sources were 'oversplit', where multiple sources in the photometric catalog appear, based on visual inspection, to be part of a single large galaxy (this is because the Source Extractor run was optimized to find high redshift, faint sources). For these galaxies, we further inspect the segmentation map to determine the neighboring components that are likely to be from the same galaxy. The photometry of each component was combined if the individual best-fitting redshift was within ∆ z ≤ 0 . 5. The errors on flux are added in quadrature and the updated photometry was then refitted both with eazy and dense basis as de-  \nbed above. After correcting the photometry, 50/53 sources were still considered high redshift and massive.  \nThis resulted in 261 sources with robust photometry which satisfied our selection. While the mass limit in this initial selection implies that many of these galaxies in fact have median stellar masses estimated to be below our intended mass limit, we draw stellar mass posteriors from this larger sample when measuring the volume densities of massive galaxies (details in § 4.1), to account for the tail of their mass posterior which exceeds our mass limit of log 10 ( M ⋆ /M ⊙ ) > 10. Figure 1 shows the redshift and stellar mass distribution recovered here.  \nWhile we expect there to exist similar contaminants and spurious sources in the lower redshift sample, we elect not to inspect the complete sample as this would be unfeasible for such a high number of sources, we do however inspect a random subsample of the low redshift sample and found a very low spurious rate, therefore the expected contamination is not expected contribute significantly to the derived statistics from the sample due to the large number of real sources. Therefore the following discussion primarily focuses on the visually inspected high redshift subsample of 261 massive galaxies presented here.", '3.3. Additional Data': 'The CEERS field includes coverage by both MIRI imaging (Wright et al. 2023), NIRSpec spectroscopy (Boker et al. 2023), and NIRCam Wide field slitless spectroscopy (WFSS; Rieke et al. 2023). NIRSpec allows for detection of strong [O iii ] emission out to z < 9 . 6, while the CEERS WFSS observations with the F356W filter are sensitive to [O iii ] redshifts at 5 < z < 7, allowing spectroscopic confirmation for sources with strong emission lines in our sample.', '3.3.1. MIRI': 'Eighteen high redshift massive galaxies fall within the area covered by the CEERS MIRI imaging, extending the wavelength coverage for SED fitting to longer wavelengths. The MIRI imaging is conducted over eight pointings, four blue pointings (P3, P6, P7, P9) providing deep imaging with F560W and F770W, and four with contiguous wavelength coverage in F1000W, F1280W, F1500W, and F1800W (P1, P2, P5, P8). Two of these pointings (P1, P2) also include coverage in F770W and F2100W. These observations were performed in two epochs, in June 2022 and December 2022. Pointings 1, 2, 3, and 6 were observed in Epoch 1 as prime observations with parallel NIRCam imaging, and pointings 5, 7, 8, and 9 were observed in parallel to NIRCam Wide Field Slitless Spectroscopy. Details of the MIRI data reduction can be found in Yang et al. (2023).  \nFigure 2. Left: Estimated stellar masses for galaxies with and without MIRI photometry. With one exception (the Larson et al. (2023) galaxy, which also has spectroscopic confirmation), including MIRI photometric constraints does not significantly alter the recovered stellar masses (average ∆ M ⋆ < ∼ 0 . 1 dex), though it does reduce the median error by ∼ 0 . 05 dex. Right: Comparison of photometric to measured spectroscopic redshifts for sources covered by NIRSpec (five sources) and WFSS (22 sources). Similar to previous analyses of JWST spectroscopy (e.g. Fujimoto et al. 2023; Arrabal Haro et al. 2023), we do observe a slight systematic offset between photometric redshift estimates and confirmed spectroscopic redshifts, with all but two galaxies having confirmed spectroscopic redshifts lower than the best fit photometric estimate. However, we find general good agreement between z phot and z spec , with only four sources with ∆ z > 0 . 5, and no catastrophic ouliers. We update our SED fitting with spectroscopic redshifts when available.  \n<!-- image -->  \n<!-- image -->  \nz  \nSource Extractor v2.19.5 (Bertin & Arnouts 1996a) was run on the blue pointings using F560W+F770W as detection images, and using the PSF-matched F560W and F770W images and their associated RMS maps for flux measurements (Papovich et al. 2023). In the red MIRI pointings, photometry extraction was determined using TPhot v2.1 (Merlin et al. 2016) following Yang et al. (2021). TPhot was used for the red pointings instead of Source Extractor due to the large wavelength coverage of these pointings and to mitigate quality loss by PSF matching across filters. Of the galaxies presented here, 18 were covered in MIRI (16 in the deeper blue pointings, and 2 in the red pointings). We perform our SED fitting including MIRI constraints and re-derive their physical parameters. We find that the stellar masses largely agree between SED fitting with and without MIRI observations within errors, however, including MIRI photometry slightly decreases the stellar masses of the sample by ∆ M ⋆ ∼ 0 . 02 dex (Figure 2), with 1/18 source (the Larson et al. 2023 galaxy) having a decrease in measured stellar mass of 0.5 dex, similar to results found by Papovich et al. (2023). The inclusion of these longer wavelength observations does provide marginally tighter constraints on the estimated stellar masses, reducing the median error bar by ∼ 0 . 05 dex.', '3.3.2. Spectroscopy': "For this sample, five sources were observed with NIRSpec with the prism, and one object was also observed  \nwith the medium-resolution gratings. The data processing is explained in detail in Arrabal Haro et al. (in prep.). However, the main steps of the reduction follow those employed in Arrabal Haro et al. (2023); Fujimoto et al. (2023); Kocevski et al. (2023); Larson et al. (2023), and is briefly summarized here. These NIRSpec data were reduced using JWST Calibration Pipeline v1.8.5 and CRDS mapping 1045. Custom parameters were utilized to improve correction on 'snowball' events. The resulting count-rate maps are processed through stage two using calwebb\\_spec2 pipeline modules, the images of the three nods are combined to extract the one-dimensional spectra.  \nFor all six sources, strong emission lines were observed, including [O iii ] and Hα , providing clear spectroscopic redshift confirmations. One of these sources was previously presented in Larson et al. (2023), thus we use their measurements (the source with medium-resolution grating observations). For the remaining five we measure spectroscopic redshifts by fitting Gaussian functions to observed lines, converting to redshifts using the relevant vacuum rest-frame wavelength.  \nEighty-two sources were covered by F356W WFSS observations. Of these 82 sources, 22 have clearly identifiable emission lines in both the row and column dispersions. The WFSS spectra were extracted following the Simulation Based Extraction method which was described in Pirzkal et al. (2017). To summarize, we used the CEERS master mosaic and available photometry to  \n)  \nf  \n(  \nx  \nu  \nl  \nf  \nFigure 3. Examples of spectroscopic observations of sources within our sample of log 10 ( M ⋆ /M ⊙ ) > 10 galaxies. The top five panels show spectra for the five sources with NIRSpec PRISM observations (one additional source, presented in Larson et al. 2023, has a NIRSpec grating redshift). Each galaxy shows clearly detected [O iii ] doublet+H β and H α emission. The bottom eight panels show line detections in the GRISM, with the top row showing H α and the bottom row shows the [O iii ] doublet+H β . The red (blue) shows the row (column)-dispersion. Here we smooth the GRISM data using a Gaussian filter with a width of one spectral element for clarity. 28 sources have clear emission lines providing robust spectroscopic redshift confirmation of many of the sources in our sample.  \n<!-- image -->  \nrest wavelength ( m)  \nduly simulate each individual WFSS observations and these served to estimate and correct for spectral contamination in the 2D WFSS wavelength calibrated spectra, before these were combined and used to produce a pair of independent (GRISMR and GRISMC) 1D spectra using optimal extraction. We used the latest, most up-to-date version of the NIRCAM WFSS calibration products which were finalized in August 2023. The estimated accuracy of the field dependence of the calibration is 0.2 pixel while the estimated uncertainty in the final wavelength calibration is 2.5 ˚ A. We inspected each of the objects by eye and identified objects with obvious emission lines present in both the row and column dispersion data. We conservatively require detections in both observations to ensure no spurious lines were included; several sources contain partial observations in only one dispersion direction, thus future work may increase the yield of spectroscopic observations. We used a combination of the photometric redshifts and the presence or lack of a clear line doublet to classify the strong emission lines as H α or [O iii ] λλ 5007,4959 lines. For the H α objects, we fit a single Gaussian function with a local continuum offset to the emission line, and for the [O iii ] objects we fit two Gaussians with a fixed 2.98:1 amplitude ratio, a fixed line-center ratio of 5006.843:4958.911, a common Gaussian width, and a local continuum offset. Using a Monte Carlo technique, we fit the appropriate model to each object 1000 times, and varied the flux data of the object spectrum using the spectral error data between each run. We use the median line centers and their standard deviations from these Monte Carlo runs to determine spectroscopic redshifts with uncertainties for each object. We report these values in Table 1.  \nFor the majority of galaxies with spectroscopic coverage, we find remarkable agreement between photometric and spectroscopic redshifts, with four sources with ∆ z > 0 . 5 and only one source with ∆ z > 1 between the photometric and spectroscopic redshifts (Figure 2, right panel). Notably, we find no catastrophic outliers, with all galaxies remaining in our sample with their spectroscopic redshifts, providing confidence in the validity of the full sample of massive galaxies. We show observed spectra for the five NIRSpec prism sources and eight example NIRCam WFSS detections in Figure 3, all showing clear emission line detections. We again perform our SED fitting with this updated redshift information.", '3.4. Extremely Red Objects': "One of the early JWST discoveries is the prominence of a new class of sources: high redshift point-like  \nFigure 4. Color selection of 'extremely red objects' (EROs) (Barro et al. 2023) colored by redshift. These sources are bright and point-source like in F444W, exhibiting blue restUV colors with redder rest-optical colors. Due to the brightness of these sources at λ observed > 3 µ m, SED fits generally recover large stellar masses. However, these galaxies have photometric colors consistent (and degenerate) with different combinations of dust-obscured massive star-forming galaxies, less massive non-dust-obscured galaxies, and dust reddened AGNs. Due to these degeneracies, the estimated stellar masses can vary significantly, up to ∼ 3 dex, dependent on assumed galaxy components. Therefore, we remove these sources from our fiducial analysis; however, we still show cumulative number densities found inclusive of these objects for comparison (Figure 5).  \n<!-- image -->  \nsources that are blue or flat in the short wavelength channels in JWST but show extreme reddening towards longer wavelengths. Barro et al. (2023) analyzed these sources and found that the most likely explanation of the photometric colors are various combinations of dustreddened AGN with bluer stellar components. Indeed, early JWST observations seem to indicate that accreting super-massive black holes are relatively common at z > 5 (e.g. Kocevski et al. 2023; Larson et al. 2023; Leung et al. 2023; Harikane et al. 2023; Labbe et al. 2023; Matthee et al. 2023; Furtak et al. 2023; Juodˇzbalis et al. 2023). This interpretation is supported by a number of such sources with confirmed broad-line AGN (Kocevski et al. 2023; Greene et al. 2023), though weaker than expected long-wavelength MIRI fluxes may indicate this population also includes non-AGNs (Williams et al. 2023).  \nWith only photometric colors available, it is extremely difficult to accurately determine the light contributed by the AGN component of these galaxies, making photometric stellar mass estimates for these sources extremely uncertain (Barro et al. 2023; Kocevski et al. 2023). These degeneracies mean that estimated stellar masses may vary significantly, up to ∼ 3 dex depen-  \nTable 1. Properties of Massive Galaxy Sample  \nNote -Table reporting redshifts and stellar masses of z > 4 massive galaxy candidates, sorted by mass. Sources with auxiliary data are marked accordingly, the full table can be found in the electronic version of this paper.  \ndent on assumed galaxy components (Barro et al. 2023). These uncertainties in stellar mass would bias our results in determining the cumulative number densities of massive galaxies observed; therefore we remove these sources from our fiducial sample with a straightforward color cut of F277W-F444W > 1 . 5, which is found to effectively identify sources dominated by potential AGN emission (Barro et al. 2023). We show the color-space of our sources in Figure 4, as well as the proposed colorcut to remove these objects. We find nine sources which satisfy this color cut, many of which are fit as higher redshift sources ( z > 7). We further visually inspect these and find that they are indeed point-like in F444W. With this in mind, we perform the following analysis on two samples, one inclusive of all sources determined to have stellar masses of > 10 10 M ⊙ , and the same sample, but with the red sources described here removed.  \nIn summation, we use SED fitting to identify galaxies with log 10 ( M ⋆ /M ⊙ ) > 10 at z = 1 -9 in the CEERS survey. We visually inspect sources in the high redshift sample which are all required to have median z phot > 3 . 5 and P( z > 2 . 5) > 70%; 28 of these galaxies have spectroscopic redshifts. We also identify 9 extremely reddened objects via their F277W-F444W colors and explore the implications of the following results both with and without EROs included in the sample.", '4.1. Cumulative Number Densities of Massive Galaxies': "To determine the cumulative number densities of massive galaxies across cosmic time, we separate our sample into redshift bins of ∆ z = 1. At z < 4, we take the median fitted physical parameters for each source recovered from dense basis . Above z > 4 we draw from the recovered posterior distribution for each individual (visually inspected) source. For each source, we conduct a draw from the posterior SED to obtain a stellar mass and redshift (stellar mass and redshift sampled as a pair as they are co-variant), and assign the galaxy to a redshift bin if the drawn stellar mass is greater than our intended mass limit (log 10 ( M ⋆ /M ⊙ ) > 10). We perform this sampling over 500 draws to recover the number of massive galaxies across redshift bins. For each draw, we further sample from the Poisson distribution of the number of galaxies in each redshift bin. Combined, these provides uncertainties in the number densities inclusive of the spread in redshift, stellar mass, and Poisson uncertainties. To account for cosmic variance, or the variance that arises from the effects of large-scale structure, we use the approach presented in Yung et al. (2023c) and calculate the variance from 2500 random samples of a 6' × 14.7' elongated field from 5 of the 2-deg 2 light-  \nresented in that work. In Figure 5, we report the median as the number densities and show the 68% confidence limit as the shaded region, the CV is shown as error bars.  \nWe show in Figure 5 the cumulative number densities of log 10 ( M ⋆ /M ⊙ ) > 10 galaxies as recovered from this analysis in both our entire galaxy sample, and the same sample with red point-sources removed ( § 3.4) in the CEERS field, as a function of redshift. For comparison, we also plot curves of cumulative number densities of log 10 ( M ⋆ /M ⊙ ) > 10 galaxies given various values of constant baryon conversion efficiency, or the rate of conversion from baryonic to stellar mass using halo mass functions from the Bolshoi-Planck and MultiDarkPlanck ΛCDM cosmological simulations (Rodr'ıguezPuebla et al. 2016). Assuming a cosmic baryon fraction f baryon = 0.16, which is consistent with the value established by observations of the cosmic microwave background (Spergel et al. 2003; Jarosik et al. 2011), we convert the underlying halo mass function to a stellar mass function. We then integrate the converted stellar mass function of log 10 ( M ⋆ /M ⊙ ) > 10 sources assuming varying values of efficiency, ϵ .  \nIn CEERS, we see that above a redshift of z ∼ 4, while the number density of massive galaxies continues to decrease with increasing redshift, there is an increased abundance over that expected with a constant baryon conversion efficiency as we progress to higher redshifts. This increased abundance could be explained by a difference in star-formation physics in the early Universe, resulting in an increased global baryon conversion efficiency and a non-standard galaxy initial mass function, discussed in § 5.", '4.2. Sizes of Massive Galaxies': "In Figure 6, we present the size-mass distribution of our massive sample from 4 < z < 8. Sizes were measured with galfit and are effective semi-major axes measured at rest-frame 0.5 µ m. Details of the galfit fitting will be provided in a future paper (McGrath et al., in prep.). Briefly, sources were fit using empirical PSFs (Finkelstein et al. 2023b), and the ERR array was used as the input sigma image, which includes noise elements from the detector, sky, and source. Neighboring galaxies were either fit simultaneously or masked during fitting, depending on brightness and distance from the primary source. Only galaxies whose best-fit models yielded magnitudes consistent with our photometry and which did not reach a constraint limit are shown (83/118 sources). We also show the best-fit size-mass relations from van der Wel et al. (2014) extrapolated to this redshift range for comparison.  \nFigure 5. The redshift evolution of massive (log 10 ( M ⋆ /M ⊙ ) > 10) galaxies in CEERS. The black solid (dotted) line shows the volume density for log 10 ( M ⋆ /M ⊙ ) > 10 galaxies exclusive (inclusive) of EROs, with the darker (lighter) shaded region showing the 68% confidence interval (CI) inclusive of sampling the redshift and stellar mass posteriors as well as the Poisson error, the errorbars show the error from cosmic variance. We show a variety of preJWST results (green-blue points) and comparison to the massive galaxies presented in Labb'e et al. (2023), all converted to a Chabrier IMF. The orange-red curves show the expected cumulative number densities above our mass threshold for a given baryon conversion efficiency ( ϵ ) calculated using the Rodr'ıguez-Puebla et al. (2016) halo mass function. Our results are systematically higher than previous studies though often within uncertainty, while following a similar evolution to higher redshifts (details in § 5) We find that the volume density of massive galaxies follow that expected given a constant baryon conversion efficiency of ϵ ∼ 0 . 14 up to z ∼ 5, before which the volume densities begin to exceed predicted values.  \n<!-- image -->  \nThe slope in the size-mass plane for our massive galaxies sample is consistent with results for star-forming galaxies at lower redshift, but when we extrapolate the evolving intercept of the size-mass relation from lower redshift to z > 4, we find that galaxies at 4 < z < 8 are smaller on average by a factor of 1 . 5 -2 than predicted from previous work (van der Wel et al. 2014). This is consistent with the results of Ormerod et al. (2023) and Ward et al. (2023), who also found evolution towards smaller sizes at higher redshifts. Furthermore, we see evidence for continued redshift evolution, with z = 6 -7 galaxies smaller at fixed mass than z = 4 -5 galaxies by a factor of 1.67, implying that the processes responsible for driving this evolution are already in place at very early times in the universe.", '4.3. Formation Timescales': 'Star-formation histories, or the SFR of a galaxy across its lifetime, probes galaxy formation and evolution. By studying the star-formation histories of massive galaxies in particular, we can attempt to understand how these galaxies manage to form surprisingly large amounts of stellar mass over relatively short timescales, and provide  \nvaluable insights into the early stages of galaxy formation and crucial benchmarks for galaxy formation theory.  \nFrom dense basis , we estimate the times at which each galaxy formed 25%, 50%, 75% of its total stellar mass at the time of observation and the associated spread in each timescale. Figure 7 shows the distribution of τ 75 (the difference between the time of observation compared to the time when the galaxies formed 75% of their stellar mass) for our sample with median log 10 ( M ⋆ /M ⊙ ) > 10 with EROs removed (120 sources). We separate our sample into 5 bins of equal number of galaxies across stellar mass and show the median τ 75 and 68% spread in each bin. We notice that there may be a slight dependence on τ 75 at which 75% of the observed stellar mass formed with respect to current observed stellar mass, similar to results found by EstradaCarpenter et al. (2020); however the spread and errors are large, making this relationship difficult to robustly constrain.  \nWe show some example fitted SFHs in Figure 8. With dense basis , we recover a diversity of SFH shapes, including periods of bursts, and rising/falling SFHs. We  \nFigure 6. The size-mass distribution of z > 4 massive galaxies presented here. We show the effective semi-major axes measured at rest-frame 0.5 µ m. The points are colored by redshift and the shape indicates the 10Myr averaged sSFR values from dense-basis . The best-fit size-mass relations for quiescent and star-forming galaxies from van der Wel et al. (2014) extrapolated to z = 4 -8 is also shown for comparison. The grey line shows the size range of a point source from z = 4 -7. We find that our sample is primarily star-forming and the slope of our sample in the size-mass plane is consistent with results for star-forming galaxies at lower redshifts, however, they are smaller by a factor of ∼ 2 than predicted when the relation is extrapolated out to higher redshifts.  \n<!-- image -->  \nnote that the posterior SFHs often have a large spread in recovered shapes, with some sources exhibiting unconstrained SFHs, and others which are more constrained. Without spectral information, it is difficult to disentangle the emission/absorption lines from the continuum for each galaxy; therefore the exact stellar populations also become difficult to fully constrain, particularly for older stellar populations which are often out-shined by younger more massive stars (Papovich et al. 2023; Conroy 2013). While recent bursts of star-formation activity can be indicated by photometry, spectroscopic observations are necessary to robustly constrain the full mass assembly histories of these galaxies (Estrada-Carpenter et al. 2020, 2021). Furthermore, larger samples of massive galaxies may be able to recover more robust trends with star-formation timescales.', '5.1. Comparison to Literature Results': "Here we compare our measured cumulative number densities to those inferred from previously published stellar mass functions, converting all results to a Chabrier IMF where needed (Figure 5). At lower redshifts, ( z < 4), we compare to results from Tomczak et al. (2014), who surveyed an area of 316 arcmin 2 distributed over three independent fields, using observations from the FourStar Galaxy Evolution Survey (ZFOURGE). Our results show slightly higher number  \nFigure 7. The estimated formation time of 75% of the stellar mass (in terms of lookback time, τ 75 ) for each galaxy in our massive, high redshift sample. The points are colored by redshift, we bin galaxies into bins of stellar mass, each bin with an equal number of sources; the open black squares show the median of each bin, with the error bars denoting the 18% and 84% spread. There is a slight trend with τ 75 with respect to the observed stellar mass, with galaxies of larger mass forming the bulk of their stars at an earlier time, but this is not significant with constraints on these quantities from the current data. Due to the difficulty constraining star-formation histories with only photometry, we note there exist large errors associated with τ 75 . SFHs are often difficult to constrain based on photometry alone, with follow-up spectroscopy often necessary to properly disentangle the underlying stellar populations. Furthermore, at very high redshifts, our sample is small ( ∼ a few galaxies at z > 7); larger surveys will be able to provide the sample sizes needed to properly infer trends with redshift.  \n<!-- image -->  \ndensities that agree within the quoted uncertainties. We also show results from the K s selected catalog of the 1.62 deg 2 COSMOS/UltraVISTA field from 0 . 2 < z < 4 . 0 from Muzzin et al. (2013). Our results show higher number densities; however, we note that beyond z ≥ 2 . 0, the analysis from Muzzin et al. (2013) is limited to galaxies with log 10 ( M ⋆ /M ⊙ ) > 10 . 51, reducing its measured number densities. We also show recent work from Leja et al. (2020), who compute an updated SMF from SED fitting with non-parametric SFHs of ∼ 10 5 galaxies in the 3D-HST and COSMOS-2015 surveys to infer updated stellar masses, and constructs a continuity model that directly fits for the redshift evolution of the SMF.  \nAt higher redshifts, our results are overall higher than the stellar mass functions from Great Observatories Origins Deep Survey Re-ionization Era Wide-area Treasury from Spitzer program (GREATS: Stefanon et al. 2021). We are within the error of the stellar mass function presented by Song et al. (2016), calculated from a restframe UV selected sample of ∼ 4500 galaxies, over the ∼ 280 arcmin 2 CANDELS/GOODS HUDF fields, and  \nFigure 8. This figure set contains example plots of the recovered SFHs from dense-basis of log ( M ⋆ /M ⊙ ) > 10 and z > 4 galaxies. The blue (orange, red) panels show the time at which the galaxy formed 25% (50%, 75%) of the observed stellar mass. The black line shows the best-fitting SFH, and the thinner blue lines show other draws for the SFH; the opacity scales with χ 2 of each fit, where darker representing a lower χ 2 (better fit). The bottom two panel shows the added constraints including MIRI observations provides. The SFHs of highz massive galaxies are varied, and in certain cases, difficult to tightly constrain when determined from photometry alone.  \n<!-- image -->  \nDuncan et al. (2014) who measure the SMF in the CANDELS GOODS-South field from 4 < z < 7. Weaver et al. (2022a) studied the galaxy SMF from redshifts 0 . 2 < z ≤ 7 . 5, leveraging the 1.27 deg 2 coverage of the COSMOS2020 catalog (Weaver et al. 2022b). We find overall higher number densities across all redshifts, while following the same general trend as we push to higher redshifts.  \nWhile resulting cumulative number densities in this work are generally higher than previous works, our results make use of deeper data covering both the restframe UV and optical at these redshifts. The improvement in angular resolution compared to Spitzer/IRAC largely corrects the source blending seen with IRAC , resulting in higher completeness. We note that in particular, the number density inclusive of EROs we find here (if extrapolated) agrees within the uncertainties with previous studies of massive galaxies in the CEERS field by Labb'e et al. (2023).", '5.2. Comparisons to Models': 'We show in Figure 9 comparisons between the number densities of massive galaxies in CEERS and a variety of models, including hydrodynamical simulations Millennium-TNG (MTNG, Kannan et al. 2023; Kannan et al. 2022; Springel et al. 2018) and FLARES (Lovell et al. 2021; Wilkins et al. 2023), the empirical model Universe-Machine (Behroozi et al. 2019), and semi-analytic models DRAGONS (Wang et al. 2016) and the Santa Cruz SAM (SC-SAM, Somerville et al. 2015; Yung et al. 2019a,b), and the Feedback-Free Starbursts (FFB, Dekel et al. 2023) model.  \nWhile the models agree well at low redshift, they tend to diverge at higher redshifts due to a combination of effects from physical assumptions and the volume probed. Below a redshift of z ∼ 5 (i.e. z < 5), our results agree well with SC-SAM, MTNG, and DRAGONS, however we find higher abundances of massive galaxies than that predicted by FFB, which is made by construction to converge with the Universe-Machine at low redshifts. At z ≳ 5, we find a higher median cumulative number density than almost all of the models that we compared with except for FFB, which argues for high star-formation efficiency due to FFB in massive galaxies at cosmic dawn, where the gas density in star-forming clouds is above a critical value (Dekel et al. 2023). The shaded area corresponds to a range of maximum efficiency from ϵ max = 0 . 2 (bottom) to 1 . 0 (top). Our results at z > 6 most closely follow the FFB model with ϵ max ∼ 0 . 3.', '5.3. Evolution of Massive Galaxies': 'JWST is revealing a high abundance of bright galaxies in the early Universe (e.g. Finkelstein et al. 2023a;  \nFigure 9. Comparison of the number density of log 10 ( M ⋆ /M ⊙ ) > 10 . galaxies found in CEERS vs. models. Here we show the sampled number densities with EROs removed (details in § 3.4), the grey shaded region in our 95% CI. We note that at lower redshifts ( z < 4), models seem to agree well; however, they diverge at higher redshifts. Overall, we find good agreement with DRAGONS, MTNG, and SC-SAM at z < 4 . 5. At z > 5 we find higher median number density than predicted by all models shown aside from FFB, where our results are best fit by the FFB predictions with ϵ max = 0 . 3 at z > 6.  \n<!-- image -->  \nCastellano et al. 2022; Casey et al. 2023; Franco et al. 2023; Leung et al. 2023; Naidu et al. 2022; Adams et al. 2023; Bouwens et al. 2023; Harikane et al. 2023). These observations may be explained by changing the physical processes that govern the relationship between observed galaxy light and the host dark matter halo mass. In this analysis we find a similar trend of a high abundance of massive galaxies in the early Universe. Here, we discuss two possible physical processes that can explain this excess: 1. a higher baryon conversion efficiency at high redshifts, and 2. a changing mass-to-light ratio, where the assumptions made in converting the observed light to stellar mass grounded in the local Universe no longer hold at high redshifts. Other physical processes may also be invoked to explain excesses in UV luminous and massive galaxies, such as modifying cosmology, or varying dust attenuation and/or levels of stochasticity of star formation, though we note that the latter two are encompassed in the performed SED modelling and should not effect the recovered stellar masses significantly; we refer the reader to works by Mason et al. (2023); Finkelstein et al. (2023b); Pallottini & Ferrara (2023); Shen et al. (2023); Sun et al. (2023); Ferrara et al. (2023); Dressler et al. (2023); Boylan-Kolchin (2023); Parashari  \n&Laha (2023); Maio & Viel (2023) for further discussion on these particular topics.', '5.3.1. A Redshift Dependence in Global Baryon Conversion Efficiency': "The baryon conversion efficiency ( ϵ ) describes how efficiently the baryonic gas accreted onto dark matter halos is converted into stars. It is defined as:  \nM ⋆ = ϵf baryon M halo (1)  \nThis efficiency, as defined here, is averaged over the star-formation history of the galaxy, and is determined to strongly correlate with both halo mass as well as the stellar mass of the central galaxy (Moster et al. 2010; Behroozi et al. 2010). In the local Universe, low values of ϵ ( ≲ 0 . 15) are inferred (Bregman 2007; Zhang et al. 2022). Energetic feedback from a combination of effects including stellar winds, supernovae, AGN, mergers/shock heating and morphological quenching are thought to quench and prevent large amounts of efficient star-formation from occurring in the local Universe (e.g. Kauffmann & Charlot 1998; Grimes et al. 2009; Kondapally et al. 2023; Martig et al. 2009; Somerville & Dav'e 2015; Naab & Ostriker 2017, etc.). Dekel et al. (2023) have argued that at extremely high redshifts ( z ≳ 10), ϵ may approach unity because in the early Universe, the combination of high ISM density and low metallicities implies that molecular cloud freefall times are shorter than the time it takes for stars to develop winds and supernovae. When they explode in higher density gas, supernovae explosions become less efficient at heating and driving winds, perhaps causing weaker stellar feedback in the early universe (Walch et al. 2015).  \nTo explore how the baryon conversion efficiency may have had to evolve in order to match the observed number densities for each redshift bin, we find the cumulative number density for each value of ϵ that most closely matches the observed cumulative number density in CEERS. From this, we show how the baryon conversion efficiency would evolve with redshift to explain the surplus of massive galaxies in the early Universe. We find that, to match the observed number of log 10 ( M ⋆ /M ⊙ ) > 10 galaxies in this study, the efficiency remains relatively constant ( ϵ ∼ 0 . 14) from z = 1 to z ∼ 4, at which point it steadily increases, reaching ϵ ∼ 0 . 3 at z ∼ 7, a factor of 3 × more efficient than the average efficiency expected in the local Universe (Bregman (2007); Zhang et al. (2022); Figure 10; panel a).  \nIf the light observed from EROs is indeed primarily contributed by their stellar components, and thus the estimated masses are robust, we show that there  \nFigure 10. a) The evolution in baryon conversion efficiency given a constant mass-to-light ratio needed to match our observations of massive galaxies ( § 5.3.1). Or, b) the factor change in mass-to-light required to match the number densities of log 10 ( M ⋆ /M ⊙ ) > 10 expected assuming a constant efficiency of ϵ = 0 . 14. We find that the observed cumulative number densities in CEERS follow that expected from a constant baryon conversion efficiency of ∼ 0 . 14 up to a redshift of z ∼ 4, before which we find an excess of log 10 ( M ⋆ /M ⊙ ) > 10 galaxies, potentially indicating that the global baryon conversion efficiency is higher or that the mass-to-light ratio is lower in the early Universe than at lower redshift.  \n<!-- image -->  \nwould need to be a dramatic increase in ϵ , particularly at z > 6, where they represent a significant portion of massive galaxies in our sample. Including these sources the baryon conversion efficiency would have to have reached a value of 0 . 4 by ∼ 7. At even higher redshifts Dekel et al. (2023) predicts an even higher efficiency value (close to unity at z ∼ 10). While we do not probe to such high redshifts in this work due to our mass limit combined with the cosmic volume probed (the area coverage of CEERS is not expected to find rare, extremely massive high redshift sources Santini et al. 2021).  \nOur inference of ϵ as defined in Eq. 1 requires knowledge of the underlying dark matter halo mass function. We use mass functions obtained directly from large cosmological simulations (Bolshoi-Planck and MultiDarkPlanck) run with cosmological parameters that are consistent with the Planck cosmology (see Sec. 4.1). If the true halo mass function deviates from these simulation results, the inferred star formation efficiencies will change accordingly. However, this would only affect our results at a quantitative rather than qualitative level unless the background cosmology itself deviates from the flat Planck cosmology. For a recent detailed investigation of halo mass functions at high redshift, see Yung et al. (2023a).", '5.3.2. Changing Mass-to-Light Ratio': 'Another explanation for the implied excess of massive galaxies may be due to an incorrect assumption about the mass-to-light ratio at high redshifts. The determination of the stellar mass function (SMF) is a cornerstone in astrophysics, for the stellar mass distribution determines the evolution, surface brightness, chemical enrichment, and baryonic content of galaxies. SED fitting codes estimate the mass of the galaxy from its observed brightness, dependent on an assumed initial mass function (IMF), or the assumption of the abundance of stars at any given mass populating the galaxy. Various functional forms of the IMF are widely used (e.g. Salpeter 1955; Scalo & Miller 1979; Kroupa & Boily 2002; Chabrier 2003). The determination of the initial mass function is not direct; but instead is dependent on the transformation of the observed light from objects into mass reliant on theories of stellar evolution. Current IMFs are calibrated on a Galactic scale, where individual stars are resolvable, while the IMF characteristics of early star formation at large redshift remain undetermined. However, changes in the characteristic stellar mass for an IMF may be expected in the early Universe (Steinhardt et al. 2023; Sneppen et al. 2022). The IMF depends upon several properties of the starforming clouds, with the high temperature and low gas metallicities in the early Universe motivating the formation of more massive stars. This can lead to the decrease in mass-to-light ratio by a factor of several (e.g. Bromm et al. 1999; Raiter et al. 2010), which could also explain the high observed abundances of UV bright galaxies at z > 10 from recent JWST observations (Finkelstein et al. 2023a; Harikane et al. 2023; Yung et al. 2023a).  \nIn Figure 10, we show how the mass measured from our SED fitting would need to be adjusted in order for the cumulative number densities of massive galaxies observed to match those expected from a constant baryon conversion efficiency ( ϵ = 0 . 14). Our reported stellar masses are determined using a Chabrier (2003) IMF, thus the change in the mass-to-light ratio shown here is the deviation from such. To determine this, we repeat the procedure as detailed in § 4.1, drawing from the posteriors of the SED fits, and bin them according to mass and redshift. For each redshift bin, we then divide all masses by some value δ (starting at 1.0 with increments of 0.0001) until the number of galaxies with log 10 ( M ⋆ /M ⊙ ) > 10 is within 10% of that which is expected given an assumed constant baryon conversion efficiency of ϵ = 0 . 14 (Figure 10; panel b). This analysis shows that the mass-to-light ratio should be ∼ 2 × lower by z = 6 . 5 and ∼ 3 × lower by z = 7. These ratios are consistent with those expected of top-heavy IMFs, re-  \ng et al. (2023a,b) semi-analytic model predictions based on the GUREFT simulation suite shows that a 3 × decrease in the mass-to-light ratio can be invoked to match the excess of UV luminous galaxies and match the observed UVLF at z ∼ 11.', '6. CONCLUSION': "In this work, we study the evolution of massive (log 10 ( M ⋆ /M ⊙ ) > 10) galaxies in the CEERS survey. The launch of JWST has brought with it an unparalleled ability to measure stellar masses at redshifts beyond z > 5, which was historically limited due to the abilities of previous telescope facilities to constrain restframe optical emission.  \nTo select massive galaxies, we perform SED fitting with eazy to obtain photometric redshift probabilities and with dense basis to fit flexible non-parametric SFHs and recover posteriors on stellar masses. We perform visual inspection of all high redshift ( z > 3 . 5) sources for which the 97 th percentile posterior in stellar mass is log 10 ( M ⋆ /M ⊙ ) > 10 (561 sources), allowing us to further vet photometry and their resulting fits. We remove a total of 300 sources from our final high redshift sample and note that this significant number of contaminant sources speaks to the necessity of careful selection for these rare galaxies. We remove likely AGN contaminants from our sample, as the stellar mass estimates cannot be accurately determined without further information on light contribution from the central AGN. We then sample the posterior fit for each galaxy and take the median ± 1 σ spread and arrive at a cumulative number density inclusive of errors from SED fitting, Poisson, and cosmic variance. From this analysis, we find that the volume densities of massive galaxies in CEERS follow that expected from a constant baryon conversion ratio of ϵ ∼ 0 . 14 up to z ∼ 5, before which the volume densities exceed a constant ϵ and appear to indicate higher ϵ at higher redshifts.  \nWe find that the number of z > 4 massive galaxies in CEERS exceeds most theoretical predictions, particularly at redshifts above z ∼ 5. While we do find a higher than expected cumulative number density at z > 4, our findings are entirely consistent with the ΛCDM cosmological model if we posit an increased efficiency of conversion of baryons into stars in massive dark matter halos with increasing redshift. We show that a redshift-dependent global baryon conversion efficiency or a change to the mass-to-light ratios of galaxies' stellar populations, e.g., due to the variation in the IMF can both reproduce the number of observed massive galaxies at high redshifts. To disentangle these scenarios, further observations to measure galaxy clustering may  \nprovide insights into star-formation efficiencies (Mu˜noz et al. 2023) while deep spectroscopic observations to detect very massive stars via strong hydrogen and helium recombination lines alongside a lack of heavy metal features to provide constraints on assumed IMFs (Trinca et al. 2023; Yung et al. 2023a; Venditti et al. 2023).  \n- K.C. and S.L.F. acknowledge support from the Na1\n- tional Science Foundation through grant AST-2009905 2\n- and from NASA through STScI award JWST-ERS-1345. 3\n- This material is based upon work supported by the 4  \n5  \nNational Science Foundation Graduate Research Fel-  \n- lowship under Grant No. DGE 2137420. MBK ac6\n- knowledges support from NSF CAREER award AST7\n- 1752913, NSF grants AST-1910346 and AST-2108962, 8\n- NASA grant 80NSSC22K0827, and HST-AR-15809, 9\n- HST-GO-15658, HST-GO-15901, HST-GO-15902, HST10\n- AR-16159, HST-GO-16226, HST-GO-16686, HST-AR11\n- 17028, and HST-AR-17043 from the Space Telescope 12\n- Science Institute, which is operated by AURA, Inc., un13\n- der NASA contract NAS5-26555. 14", 'REFERENCES': 'Adams, N. J., Bowler, R. A. A., Jarvis, M. J., Haußler, B., & Lagos, C. D. P. 2021, MNRAS, 506, 4933,  \ndoi: 10.1093/mnras/stab1956  \nAdams, N. J., Conselice, C. J., Austin, D., et al. 2023, arXiv e-prints, arXiv:2304.13721, doi: 10.48550/arXiv.2304.13721  \nArrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, ApJL, 951, L22, doi: 10.3847/2041-8213/acdd54  \nAstropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33,  \ndoi: 10.1051/0004-6361/201322068  \nAstropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f  \nBagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12, doi: 10.3847/2041-8213/acbb08  \nBarro, G., Perez-Gonzalez, P. G., Kocevski, D. D., et al. 2023, arXiv e-prints, arXiv:2305.14418,  \ndoi: 10.48550/arXiv.2305.14418  \nBehroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182  \nBehroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379, doi: 10.1088/0004-637X/717/1/379  \nBertin, E., & Arnouts, S. 1996a, A&AS, 117, 393,  \ndoi: 10.1051/aas:1996164  \n- -. 1996b, A&AS, 117, 393, doi: 10.1051/aas:1996164\n- -. 2010, SExtractor: Source Extractor, Astrophysics Source Code Library, record ascl:1010.064. http://ascl.net/1010.064  \nBhatawdekar, R., Conselice, C. J., Margalef-Bentabol, B., & Duncan, K. 2019, MNRAS, 486, 3805, doi: 10.1093/mnras/stz866  \nBoker, T., Beck, T. L., Birkmann, S. M., et al. 2023, PASP, 135, 038001, doi: 10.1088/1538-3873/acb846  \nBouwens, R. J., Stefanon, M., Brammer, G., et al. 2023,  \nMNRAS, 523, 1036, doi: 10.1093/mnras/stad1145  \nBoylan-Kolchin, M. 2023, Nature Astronomy, 7, 731,  \ndoi: 10.1038/s41550-023-01937-7  \n- Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786  \nBregman, J. N. 2007, ARA&A, 45, 221,  \ndoi: 10.1146/annurev.astro.45.051806.110619  \nBromm, V., Coppi, P. S., & Larson, R. B. 1999, ApJ, 527, L5, doi: 10.1086/312385  \nBushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, 1.10.2, Zenodo, Zenodo, doi: 10.5281/zenodo.7829329  \nByler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66 Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692  \nFacilities: HST(STIS), JWST (STSci)  \nSoftware: astropy (Astropy Collaboration et al. 2013, 2018), EAZY (Brammer et al. 2008), Source Extractor (Bertin & Arnouts 1996b), Dense-Basis (Iyer et al. 2019; Iyer & Gawiser 2017), Colormaps from Crameri (2023)  \n- Venditti, A., Graziani, L., Schneider, R., et al. 2023, Monthly Notices of the Royal Astronomical Society, 522, 3809, doi: 10.1093/mnras/stad1201\n- Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238, doi: 10.1093/mnras/stv1975\n- Wang, L., Spurzem, R., Aarseth, S., et al. 2016, Monthly Notices of the Royal Astronomical Society, 458, 1450, doi: 10.1093/mnras/stw274\n- Ward, E. M., de la Vega, A., Mobasher, B., et al. 2023, arXiv e-prints, arXiv:2311.02162,\n- doi: 10.48550/arXiv.2311.02162\n- Weaver, J. R., Davidzon, I., Toft, S., et al. 2022a, arXiv e-prints, arXiv:2212.02512, doi: 10.48550/arXiv.2212.02512\n- Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022b, ApJS, 258, 11, doi: 10.3847/1538-4365/ac3078\n- Wilkins, S. M., Vijayan, A. P., Lovell, C. C., et al. 2023, MNRAS, 519, 3118, doi: 10.1093/mnras/stac3280\n- Williams, C. C., Alberts, S., Ji, Z., et al. 2023, arXiv e-prints, arXiv:2311.07483,\n- doi: 10.48550/arXiv.2311.07483\n- Wright, A. H., Robotham, A. S. G., Driver, S. P., et al. 2017, MNRAS, 470, 283, doi: 10.1093/mnras/stx1149\n- Wright, G. S., Rieke, G. H., Glasse, A., et al. 2023, PASP, 135, 048003, doi: 10.1088/1538-3873/acbe66\n- Xiao, M., Oesch, P., Elbaz, D., et al. 2023, Massive Optically Dark Galaxies Unveiled by JWST Challenge Galaxy Formation Models.\n- https://arxiv.org/abs/2309.02492\n- Yang, G., Papovich, C., Bagley, M. B., et al. 2021, ApJ, 908, 144, doi: 10.3847/1538-4357/abd6c1\n- Yang, G., Papovich, C., Bagley, M., et al. 2023, CEERS MIRI Imaging: Data Reduction and Quality Assessment. https://arxiv.org/abs/2307.14509\n- Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Dav\'e, R. 2019a, MNRAS, 483, 2983, doi: 10.1093/mnras/sty3241\n- Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Wilkins, S. M., & Gardner, J. P. 2023a, arXiv e-prints, arXiv:2304.04348, doi: 10.48550/arXiv.2304.04348\n- Yung, L. Y. A., Somerville, R. S., Nguyen, T., et al. 2023b, arXiv e-prints, arXiv:2309.14408,\n- doi: 10.48550/arXiv.2309.14408\n- Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019b, MNRAS, 490, 2855, doi: 10.1093/mnras/stz2755\n- Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., et al. 2023c, MNRAS, 519, 1578, doi: 10.1093/mnras/stac3595\n- Zhang, Z., Wang, H., Luo, W., et al. 2022, A&A, 663, A85, doi: 10.1051/0004-6361/202142866'}