Evidence for a Shallow Evolution in the Volume Densities of Massive Galaxies at z = 4–8 from CEERS

We analyze the evolution of massive ( log 10 [ M å / M e ] > 10 ) galaxies at z ∼ 4 – 8 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 ) ﬁ tting 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 con ﬁ rmations from NIRSpec or NIRCam WFSS. We sample the recovered posteriors in stellar mass from SED ﬁ tting to infer the volume densities of massive galaxies across cosmic time, taking into consideration the potential for sample contamination by active galactic nuclei. We ﬁ nd that the evolving abundance of massive galaxies tracks expectations based on a constant baryon conversion ef ﬁ ciency in dark matter halos for z ∼ 4 – 8. At higher redshifts, we observe an excess abundance

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
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 and/or the efficiencies with which star formation occurs in massive dark matter halos, and are not in tension with modern cosmology.

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).
Observationally, extensive efforts have been made in attempting to understand how galaxy stellar masses evolve over time, with many studies measuring galaxy stellar mass functions (SMFs) spanning cosmic time (e.g., Muzzin et al. 2013;Duncan et al. 2014;Kelvin et al. 2014;Tomczak et al. 2014;Moffett et al. 2016;Song et al. 2016;Wright et al. 2017;Bhatawdekar et al. 2019;Leja et al. 2020;Adams et al. 2021;McLeod et al. 2021;Stefanon et al. 2021;Weaver et al. 2022a).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 formation history (SFH), 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 (∼0.5-1.8 μm) 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 observations to accurately attribute observed fluxes to their sources.
The 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 the JWST Near Infrared Camera (NIRcam; Rieke et al. 2023), observations of the rest-frame optical (including the SFH-sensitive 4000 Å break) are accessible out to z ∼ 10, allowing for the measurement of accurate stellar masses out to redshifts that were previously inaccessible.Nearly 2 yr since the launch of JWST, the data collected have already begun to revolutionize our understanding of the Universe, with new discoveries answering long-standing questions in the field of extragalactic science, while also providing unexpected results and generating new questions.
One 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;Finkelstein et al. 2022Finkelstein et al. , 2023Finkelstein et al. , 2024;;Naidu et al. 2022;Adams et al. 2024;Bouwens et al. 2023;Casey et al. 2024;Franco et al. 2023;Harikane et al. 2023;Leung et al. 2023).This population of galaxies may indicate evidence for different star formation and stellar feedback physics in the early Universe, leading to the earlier, more rapid formation of massive and luminous galaxies than previously thought possible (e.g., Dekel et al. 2023;Ferrara et al. 2023;Mason et al. 2023;Pallottini & Ferrara 2023;Shen et al. 2023;Yung et al. 2023a).While models shown tension with direct observables, simulation-inferred UV luminosities are very sensitive to the assumed IMFs, dust, and active galactic nucleus (AGN) contamination.Having reliable estimates for inferred 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 M M log 10  ( )  > 11), further indicating the possibility of more efficient star formation activity at higher redshift than expected (e.g., Labbé et al. 2023;Xiao et al. 2023), and/or a revision to Lambda cold dark matter (ΛCDM; Boylan-Kolchin 2023).
To improve our understanding of the abundance of massive (defined here as log 10 [M å /M e ] > 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 Cosmic Evolution Early Release Survey (CEERS) NIRCam field.We briefly introduce CEERS and the data reduction (Section 2.1) before discussing how the galaxies studied here were selected (Section 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-frame UV-optical spectra alongside increasing brightness at longer wavelengths, possibly hosting dust-reddened AGN (Barro et al. 2024;Greene et al. 2023;Kocevski et al. 2023;Labbe et al. 2023).We discuss how this category of sources may bias selection of massive galaxies (Section 3.5).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 (Section 5.1) and model predictions (Section 5.2).Finally, we discuss how the observed abundances of galaxies may indicate changes to the efficiency of star formation in galaxies (Section 5.3.1) or how the observed light from the galaxy is translated to mass (Section 5.3.2).Throughout this paper, we assume a cosmology of H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 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).

NIRCam and HST Catalog
We perform our analysis on the publicly released CEERS NIRCam mosaics (Data Releases 0.5 and 0.6, DR0.5 and DR0.6, respectively), which include imaging from 1 to 5 μm in F115W, F150W, F200W, F277W, F356W, F410M, and F444W and cover 10 pointings in the Extended Groth Strip (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 the 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 2022 December in parallel to NIRSpec microshutter assembly observations.Pointings 5, 7, 8, and 9 were also observed with NIRCam Wide Field Slitless Spectroscopy (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 the 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, and 1/f noise, and to improve the astrometric alignment.We direct the reader to Bagley et al. (2023) for further details on the reduction.
We make use of the photometric catalog from Finkelstein et al. (2024), which also includes the available HST imaging in this field from CEERS HST Data Release 1, including the Advanced Camera for Surveys F606W and F814W filters 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. (2023), with some updates to improve the accuracy of the 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. (2024) for full details on all cataloging procedures.

MIRI
The CEERS survey includes MIRI imaging conducted over eight pointings, with four blue pointings (P3, P6, P7, and P9) providing deep imaging with F560W and F770W and four with contiguous wavelength coverage in F1000W, F1280W, F1500W, and F1800W (P1, P2, P5, and P8).Two of these pointings (P1 and P2) also include coverage in F770W and F2100W.These observations were performed in two epochs, in 2022 June and 2022 December.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 WFSS.Details of the MIRI data reduction can be found in Yang et al. (2023).
Source Extractor v2.19.5 (Bertin & Arnouts 1996) was run on the blue pointings using F560W + F770W as detection images, and using the point-spread function (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.

Spectroscopy
CEERS also includes coverage with NIRSpec with the prism and F356W WFSS observations.The NIRSpec prism data processing is explained in detail in P. Arrabal Haro et al. (2024, in preparation).However, the main steps of the reduction follow those employed in Arrabal Haro et al. (2023), Fujimoto et al. (2023), Kocevski et al. (2023), andLarson et al. (2023), and are 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 of "snowball" events.The resulting count-rate maps are processed through stage two using calwebb_spec2 pipeline modules, where the images of the three nods are combined to extract the 1D spectra.
The WFSS spectra were extracted following a simulationbased extraction method, which was described in Pirzkal et al. (2017).To summarize, we used the CEERS master mosaic and available photometry to duly simulate each individual WFSS observation and these were used 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 2023 August.The estimated accuracy of the field dependence of the calibration is 0.2 pixels while the estimated uncertainty in the final wavelength calibration is 2.5 Å.

Selection of Massive Galaxies
Here, we detail our process of selecting a sample of massive galaxies.We first make use of eazy (Brammer et al. 2008) to determine photometric redshifts.We then use the derived redshift probability distribution functions (PDFs; 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.

Photometric Redshift Fitting
Photometric redshifts for the catalog were determined with eazy.We use the results from Finkelstein et al. (2024), briefly discussing the procedure here.eazy fits nonnegative 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 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 restframe UV colors expected from z > 9 galaxies (Larson et al. 2022).

Spectral Energy Distribution Fitting
We perform SED fitting for all sources with signal-to-noise ratio in F444W >3 with the dense basis SED fitting code, which is designed to robustly recover SFH constraints from galaxy SEDs.dense basis uses a flexible nonparametric SFH represented by a Gaussian mixture model (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 (Byler et al. 2017;Ferland 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 (2003) IMF.We impose a uniform (flat) prior on the specific star formation rate (sSFR) with limits on the sSFR (sSFR yr −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 e ä [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 z ( )  to a top-hat form by first fitting a Gaussian to the primary peak of the z ( )  , and taking the width of the top-hat 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.
While 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 higher-redshift 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 M M log 10 10  ( )>


, and (2) a parent high-redshift sample, with a best-fit photometric redshift of z > 3.5 (from dense basis), at least 70% of the integrated z ( )  from eazy to be at z > 2.5 (to reject objects with overly broad redshift distributions), and a less stringent mass cut of 2.5% of the posterior in stellar mass having M M log 10  ( )  > 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 whose posteriors may include a massive ( M M log 10  ( )  > 10) solution to properly encompass the uncertainties from the mass estimates in computing the volume densities of these galaxies.We thus manually inspect 561 sources satisfying these criteria.

Sample Cuts
We 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) of the 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 low-z solution with Δχ 2 < 4 from the high-z fit, and source ID 75728 had significant flux detected below the putative Lyα break at its fitted 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.
53 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 highredshift, 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 described above.After correcting the photometry, 50/53 sources were still considered high redshift and massive.
This 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 Section 4.1), to account for the tails of their mass posteriors, which exceeds our mass limit of M M log 10  ( )  > 10. Figure 1 shows the redshift and stellar mass distributions recovered for the entire catalog.We expect our sample to be highly complete, even at high redshifts.With the typical magnitude of our z > 6 sample at ∼25 mag in F444W, we expect >90% completeness based on completeness simulations run for the CEERS survey, as described in Finkelstein et al. (2023).
While 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 to 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.

Additional Data
The CEERS field includes coverage by both MIRI imaging (Wright et al. 2023), NIRSpec spectroscopy (Böker et al. 2023), and NIRCam WFSS data (Rieke et al. 2023).NIRSpec allows for the detection of strong [O III] emission out to z ∼ 8, 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.

Spectral Energy Distribution Fitting Including MIRI
18 high-redshift massive galaxies fall within the area covered by the CEERS MIRI imaging, extending the wavelength coverage for SED fitting to longer wavelengths; 16 were covered in the deeper blue pointings and two in the red pointings.We perform our SED fitting including MIRI constraints and rederive their physical parameters.We find that the stellar masses largely agree between SED fitting with and without MIRI observations within the errors; however, including the MIRI photometry slightly decreases the stellar masses of the sample by ΔM å ∼ 0.02 dex (Figure 2).
One source (the Larson et al. 2023 galaxy) resulted in a decrease in the measured stellar mass of 0.5 dex, similar to the 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.

Spectroscopic Redshifts
For this sample, five sources were observed with NIRSpec with the prism, and one object was also observed with the medium-resolution gratings.For 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 the observed lines, converting to redshifts using the relevant vacuum rest-frame wavelength.
82 sources were covered by F356W WFSS observations.Of these 82 sources, 22 have clearly identifiable emission lines in both the row and column dispersions.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 the Hα or [O III] λλ5007, 4959 lines.For the Hα objects, we fit a single Gaussian function with a local continuum offset to The darker shaded region denotes our redshift and mass ranges 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 mass 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 Section 3.5).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 M M log 10  ( ) > 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.
Figure 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., Arrabal Haro et al. 2023;Fujimoto et al. 2023), we do observe a slight systematic offset between the photometric redshift estimates and the confirmed spectroscopic redshifts, with all but two galaxies having confirmed spectroscopic redshifts lower than the best-fit photometric estimates.However, we find general good agreement between z phot and z spec , with only four sources with Δz > 0.5, and no catastrophic outliers.We update our SED fitting with spectroscopic redshifts when available.the emission line, and for the [O III] objects we fit two Gaussians with a fixed 2.98:1 amplitude ratio, a fixed linecenter 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.
For the majority of galaxies with spectroscopic coverage, we find remarkable agreement between the 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 (This table is available in its entirety in machine-readable form in the online article.)clear emission line detections.We again perform our SED fitting with this updated redshift information.

Extremely Red Objects
One of the early JWST discoveries is the prominence of a new class of sources: high-redshift point-like sources that are blue or flat in the short-wavelength channels in JWST but show extreme reddening toward longer wavelengths.Barro et al. (2024) analyzed these sources and found that the most likely explanation of the photometric colors are various combinations of dust-reddened AGN with bluer stellar components.Indeed, early JWST observations seem to indicate that accreting supermassive black holes are relatively common at z > 5 (e.g., Furtak et al. 2023;Harikane et al. 2023;Juodžbalis et al. 2023;Kocevski et al. 2023;Labbe et al. 2023;Larson et al. 2023;Leung et al. 2023;Maiolino et al. 2023;Matthee et al. 2023;Scholtz et al. 2023).This interpretation is supported by a number of such sources with confirmed broad-line AGN (Greene et al. 2023;Kocevski et al. 2023), though weaker than expected long-wavelength MIRI fluxes may indicate this population also includes non-AGN (Williams et al. 2023).
With 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. 2024;Kocevski et al. 2023).These degeneracies mean that estimated stellar masses may vary significantly, up to ∼3 dex, depending on the assumed galaxy components (Barro et al. 2024).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 via visual inspection for point-like sources in F444W.We show the color space of our sources in Figure 4 in comparison to the color selection presented by Barro et al. (2024), who have found that a straightforward color cut of F277W − F444W > 1.5 can effectively identify sources dominated by potential AGN emission.We find nine sources which satisfy this color cut, many of which are fit as higherredshift sources (z > 7).We find that all but one source is a point source in F444W and remove it from our ERO sample.We further add into our sample three point-like sources not identified from the color-color cut.
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 e , and the same sample, but with the red sources described here removed.
In summation, we use SED fitting to identify galaxies with M M log 10 10  ( )>  at z = 1-9 in the CEERS survey.We visually inspect sources in the high-redshift sample; 28 of these galaxies have spectroscopic redshifts.We also identify 11 extremely reddened point-like objects in F444W and explore the implications of the following results both with and without these EROs included in the sample.

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 covariant), and assign the galaxy to a redshift bin if the drawn stellar mass is greater than our intended mass limit ( M M log 10 10  ( )>


).We perform this sampling over 500 draws to recover the number of massive galaxies across the redshift bins.For each draw, we further sample from the Poisson distribution of the number of galaxies in each redshift bin.Combined, these provide 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 largescale structure, we use the approach presented in Yung et al. (2023b) and calculate the variance from 2500 random samples of a 6′ × 14 7 elongated field from five of the 2 deg 2 light cones presented in that work.In Figure 5, we report the median as the number densities and show the 68% confidence limit as the shaded region, where the cosmic variance is shown as error bars.
We show in Figure 5   These sources are bright and point-source-like in F444W, exhibiting blue restframe UV colors with redder rest-frame optical colors.Due to the brightness of these sources at λ observed > 3 μm, SED fits generally 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 AGN.Due to these degeneracies, the estimated stellar masses can vary significantly, up to ∼3 dex, depending on the 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).In 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 nonstandard galaxy IMF, discussed in Section 5.

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 semimajor axes measured at restframe 0.5 μm.Details of the GALFIT fitting will be provided in a future paper (E.J. McGrath et al. 2024, in preparation).Briefly, sources were fit using empirical PSFs (Finkelstein et al. 2024), and the ERR (error) 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.
The slope in the size-mass plane for our massive galaxies sample is consistent with the 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 toward 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.Evolution toward smaller sizes at higher redshifts would imply that the gas densities in these galaxies should also evolve to higher values at higher redshifts.Simulations predict that higher gas densities should give rise to a greater fraction of baryons converted into stars (Fukushima & Yajima 2021).Thus our observation of smaller sizes at higher redshift is  2016) halo mass function.Our results are systematically higher than previous studies though often within the uncertainties, while following a similar evolution to higher redshifts (details in Section 5) We find that the volume density of massive galaxies follows that expected given a constant baryon conversion efficiency of ò ∼ 0.14 up to z ∼ 5, before which the volume densities begin to exceed the predicted values.
Figure 6.The size-mass distribution of z > 4 massive galaxies presented here.We show the effective semimajor axes measured at rest-frame 0.5 μm (measured in F277W at z = 4 and F444W at z = 7).The points are colored by redshift and the shape indicates the 10 Myr 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 are also shown for comparison.The gray line shows the size range of a point source from z = 4 to 7. We find that our sample is primarily star forming and the slope of our sample in the size-mass plane is consistent with the 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.
consistent with an evolution in the efficiency of gas conversion, as discussed in the previous subsection.

Formation Timescales
SFHs, or the SFR of a galaxy across its lifetime, probes galaxy formation and evolution.By studying the SFHs 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 valuable insights into the early stages of galaxy formation and crucial benchmarks for galaxy formation theory.
From dense basis, we estimate the times at which each galaxy formed 25%, 50%, and 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 M M log 10 10 

( )>
 and z > 4 with EROs removed (118 sources).We separate our sample into five bins of an 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 the results found by Estrada-Carpenter et al. (2020); however the spread and errors are large, making this relationship difficult to robustly constrain.
We 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 note that the posterior SFHs often have a large spread in the 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 outshined by younger more massive stars (Conroy 2013;Giménez-Arteaga et al. 2023;Papovich et al. 2023;Song et al. 2023; though as shown in Figure 2 we do not find evidence that the stellar masses of our objects are underestimated when MIRI photometry is not available).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, 2020).Furthermore, larger samples of massive galaxies may be able to recover more robust trends with star formation timescales.While only with available photometry, the SFHs are difficult to constrain, follow-up detailed spectroscopic observations may help disentangle possible SFH models that occur in the early Universe, enabling further understanding at the timescales and the onset of star formation of massive galaxies.

Comparison to Literature Results
Here we compare our measured cumulative number densities to those inferred from previously published SMFs, 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.Our results show slightly higher number densities 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 M M log 10.51 10 

( )>
 , reducing its measured number densities.We also show recent work from Leja et al. (2020), who compute an updated SMF from SED fitting with nonparametric SFHs of ∼10 5 galaxies in the 3D-HST and COSMOS-2015 surveys to infer updated stellar masses, and construct a continuity model that directly fits for the redshift evolution of the SMF.
At higher redshifts, our results are overall higher than the SMFs from Great Observatories Origins Deep Survey Reionization Era Wide-area Treasury from Spitzer program (Stefanon et al. 2021).We are within the error of the SMF presented by Song et al. (2016), calculated from a rest-frame UV-selected sample of ∼4500 galaxies, over the ∼280 arcmin 2 CANDELS/GOODS HUDF fields, and Duncan 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.
While the resulting cumulative number densities in this work are generally higher than in previous works, our results make use of deeper data covering both the rest-frame 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 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 the constraints on these quantities from the current data.Due to the difficulty constraining SFHs 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 (of order a few galaxies at z > 7); larger surveys will be able to provide the sample sizes needed to properly infer trends with redshift.
previous studies of massive galaxies in the CEERS field by Labbé et al. (2023).
While 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, The bottom two panels show the added constraints that including MIRI observations provide.The SFHs of high-z massive galaxies are varied, and in certain cases, difficult to tightly constrain when determined from photometry alone.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 0.2 max =  (bottom) to 1.0 (top).Our results at z > 6 most closely follow the FFB model with 0.3 max  .

Evolution of Massive Galaxies
JWST is revealing a high abundance of bright galaxies in the early Universe (e.g., Castellano et al. 2022;Naidu et al. 2022;Adams et al. 2024;Bouwens et al. 2023;Casey et al. 2024;Finkelstein et al. 2023;Franco et al. 2023;Harikane et al. 2023;Leung 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, varying dust attenuation, and/or levels of stochasticity of star formation, though we note that the latter two are encompassed in the performed SED modeling and should not effect the recovered stellar masses significantly; we refer the reader to the works by Boylan-Kolchin (2023),

A Redshift Dependence in the 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: where M å is the stellar mass, f baryon is the fraction of baryons available, and M halo is the mass of the underlying dark matter halo.This efficiency, as defined here, is averaged over the SFH of the galaxy, and is determined to strongly correlate with both halo mass as well as the stellar mass of the central galaxy (Behroozi et al. 2010;Moster 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;Martig et al. 2009;Somerville & Davé 2015;Naab & Ostriker 2017;Kondapally et al. 2023, 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 higherdensity gas, supernova explosions become less efficient at heating and driving winds, perhaps causing weaker stellar feedback in the early Universe (Walch et al. 2015).
To 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 M M log 10 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)).
If the light observed from EROs is indeed primarily contributed by their stellar components, and thus the estimated masses are robust, we show that there would need to be a dramatic increase in ò, particularly at z > 7, where they represent a significant portion (5/13) of massive galaxies in our sample.Including these sources the baryon conversion efficiency would have to have reached a value of 0.4 by z ∼ 7.At even higher redshifts Dekel et al. (2023) predict an even higher efficiency value (close to unity at z ∼ 10).We do not probe to such high redshifts in this work due to our mass limit  galaxies found in CEERS vs. models.Here we show the sampled number densities with EROs removed (details in Section 3.5), the gray shaded region in our 95% CI.We note that at lower redshifts (z < 4), the 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 a higher median number density than predicted by all models shown aside from FFB, where our results are best fit by the FFB predictions with 0.3 max =  at z > 6.
combined with the cosmic volume probed (the area coverage of CEERS is not expected to find rare, extremely massive highredshift sources; Santini et al. 2021).Our inference of ò as defined in Equation (1) requires knowledge of the underlying dark matter halo mass function.We use mass functions obtained directly from large cosmological simulations (Bolshoi-Planck and MultiDark-Planck) run with cosmological parameters that are consistent with the Planck cosmology (see Section 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).

Changing Mass-to-light Ratio
Another explanation for the implied excess of massive galaxies may be due to an incorrect assumption about the massto-light ratio at high redshifts.Determination of the SMF is a cornerstone in astrophysics, for the stellar mass distribution determines the evolution, surface brightness, chemical enrichment, and baryonic content galaxies.SED fitting codes estimate the mass of a galaxy from its observed brightness, dependent on an assumed 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 IMF 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 (Sneppen et al. 2022;Steinhardt et al. 2023).The IMF depends upon several properties of the star-forming clouds, with the high temperatures and low gas metallicities in the early Universe motivating the formation of more massive stars.This can lead to a decrease in the 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. 2023;Harikane et al. 2023;Yung et al. 2023a).
In Figure 10, we show how the masses 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 Section 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 M M log 10 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; e.g., the results from Yung et al. (2023aYung et al. ( , 2024) ) are SAM predictions based on the GUREFT simulation suite, which show that a 3× decrease in the mass-tolight ratio can be invoked to match the excess of UV-luminous galaxies and match the observed UV luminosity function at z ∼ 11.

Conclusion
In this work, we study the evolution of massive ( M M log 10 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 rest-frame optical emission.
To select massive galaxies, we perform SED fitting with eazy to obtain photometric redshift probabilities and with dense basis to fit flexible nonparametric SFHs and recover posteriors on stellar masses.We perform visual inspection of all high-redshift (z > 3.5) sources for which the 97th percentile posterior in stellar mass is M M log 10 10  ( )>  (561 sources), allowing us to further vet the photometry and their resulting fits.We remove a total of 300 sources from our final highredshift 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 the 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  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 M M log 10 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.
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.
We 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 variations 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 provide insights into star formation efficiencies (Muñoz 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. 2024;Venditti et al. 2023;Yung et al. 2023a).

Figure 1 .
Figure1.The distributions in redshift and stellar mass of the full CEERS catalog(Finkelstein et al. 2024), as measured by dense basis (with eazy priors on the dense basis-derived redshift).The darker shaded region denotes our redshift and mass ranges 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 mass 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 Section 3.5).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 M M log 10

Figure 3 .
Figure 3. Examples of spectroscopic observations of sources within our sample of M M log 10 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 at λ rest = 5007 Å and 6563 Å, respectively.The bottom eight panels show line detections made with 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.
as recovered from this analysis in both our entire galaxy sample, and the same sample with red point sources removed (Section 3.5) in the CEERS field, as a function of redshift.For comparison, we also plot curves of cumulative number densities of M constant baryon conversion efficiency, or the rate of conversion from baryonic to stellar mass using halo mass functions from the Bolshoi-Planck and MultiDark-Planck ΛCDM cosmological simulations(Rodríguez-Puebla 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;

Figure 4 .
Figure 4. Color-color space of EROs(Barro et al. 2024) colored by redshift, where sources with red squares are selected as ERO based on visual inspection.These sources are bright and point-source-like in F444W, exhibiting blue restframe UV colors with redder rest-frame optical colors.Due to the brightness of these sources at λ observed > 3 μm, SED fits generally 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 AGN.Due to these degeneracies, the estimated stellar masses can vary significantly, up to ∼3 dex, depending on the 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 (Figure5).
Jarosik et al. 2011), we convert the underlying halo mass function to an SMF.We then integrate the converted SMF of M assuming varying values of the efficiency, ò.

Figure 5 .
Figure 5.The redshift evolution of massive ( M M log 10  ( ) 

Figure 7 .
Figure 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, where we bin galaxies into bins of stellar mass, with each bin having 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 the constraints on these quantities from the current data.Due to the difficulty constraining SFHs 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 (of order a few galaxies at z > 7); larger surveys will be able to provide the sample sizes needed to properly infer trends with redshift.

Figure 8 .
Figure8.This figure set contains example plots of the recovered SFHs from dense basis of log (M å /M e ) > 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 starting time is defined as the start of the Universe.The black line shows the best-fitting SFH, and the thinner blue lines show other draws for the SFH; the opacity scales with the χ 2 of each fit, where darker represent a lower χ 2 (better fit).The bottom two panels show the added constraints that including MIRI observations provide.The SFHs of high-z massive galaxies are varied, and in certain cases, difficult to tightly constrain when determined from photometry alone.

Figure 9 .
Figure 9.Comparison of the number density of M M log 10. 10  ( )>

Figure 10 .
Figure 10.(a) The evolution in baryon conversion efficiency given a constant mass-to-light ratio (M/L) needed to match our observations of massive galaxies (Section 5.3.1).(b) The factor change in the mass-to-light ratio required to match the number densities of M M log 10 10  ( )>

Table 1
Properties of the Massive Galaxy Sample Table reporting the 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.