# Life beyond 30: Probing the −20 < *M*<sub>UV</sub> < −17 Luminosity Function at 8 < *z* < 13 with the NIRCam Parallel Field of the MIRI Deep Survey

Pablo G. Pérez-González<fnref n="1" />, Luca Costantin<fnref n="1" />, Danial Langeroodi<fnref n="2" />, Pierluigi Rinaldi<fnref n="3" />, Marianna Annunziatella<fnref n="1" />, Olivier Ilbert<fnref n="4" />, Luis Colina<fnref n="1" />, Hans Ulrik Nørgaard-Nielsen<fnref n="5" />,<fnref n="6" />, Thomas R. Greve<fnref n="5" />,<fnref n="6" />,<fnref n="7" />, Göran Östlin<fnref n="8" />, Gillian Wright<fnref n="9" />, Almudena Alonso-Herrero<fnref n="10" />, Javier Álvarez-Márquez<fnref n="1" />, Karina I. Caputi<fnref n="3" />, Andreas Eckart<fnref n="11" />, Olivier Le Fèvre<fnref n="4" />, Álvaro Labiano<fnref n="12" />, Macarena García-Marín<fnref n="13" />, Jens Hjorth<fnref n="2" />, Sarah Kendrew<fnref n="13" />, John P. Pye<fnref n="14" />, Tuomo Tikkanen<fnref n="14" />, Paul van der Werf<fnref n="15" />, Fabian Walter<fnref n="16" />, Martin Ward<fnref n="17" />, Arjan Bik<fnref n="8" />, Leindert Boogaard<fnref n="18" />, Sarah E. I. Bosman<fnref n="16" />, Alejandro Crespo Gómez<fnref n="1" />, Steven Gillman<fnref n="5" />,<fnref n="6" />, Edoardo Iani<fnref n="3" />, Iris Jermann<fnref n="5" />,<fnref n="6" />, Jens Melinder<fnref n="8" />, Romain A. Meyer<fnref n="18" />, Thibaud Moutard<fnref n="4" />, Ewine van Dishoek<fnref n="15" />, Thomas Henning<fnref n="18" />, Pierre-Olivier Lagage<fnref n="19" />, Manuel Guedel<fnref n="18" />,<fnref n="20" />,<fnref n="21" />, Florian Peissker<fnref n="11" />, Tom Ray<fnref n="22" />, Bart Vandenbussche<fnref n="23" />, Ángela García-Argumánez<fnref n="24" />,<fnref n="25" />, and Rosa María Mérida<fnref n="1" />

## Abstract

We present the ultraviolet luminosity function and an estimate of the cosmic star formation rate density at 8 < *z* < 13 derived from deep NIRCam observations taken in parallel with the MIRI Deep Survey of the Hubble Ultra Deep Field (HUDF), NIRCam covering the parallel field 2. Our deep (40 hr) NIRCam observations reach an F277W magnitude of 30.8 (5σ), more than 2 mag deeper than JWST public data sets already analyzed to find high-redshift galaxies. We select a sample of 44 *z* > 8 galaxy candidates based on their dropout nature in the F115W and/or F150W filters, a high probability for their photometric redshifts, estimated with three different codes, being at *z* > 8, good fits based on χ<fnref n="2" /> calculations, and predominant solutions compared to *z* < 8 alternatives. We find mild evolution in the luminosity function from *z* ~ 13 to *z* ~ 8, i.e., only a small increase in the average number density of ~0.2 dex, while the faint-end slope and absolute magnitude of the knee remain approximately constant, with values α = −2.2 ± 0.1, and *M* = −20.8 ± 0.2 mag. Comparing our results with the predictions of state-of-the-art galaxy evolution models, we find two main results: (1) a slower increase with time in the cosmic star formation rate density compared to a steeper rise predicted by models; (2) nearly a factor of 10 higher star formation activity concentrated in scales around 2 kpc in galaxies with stellar masses ~10⁸ *M*⊙ during the first 350 Myr of the universe, *z* ~ 12, with models matching better the luminosity density observational estimations ~150 Myr later, by *z* ~ 9.

*Unified Astronomy Thesaurus concepts: Galaxy formation (595); Galaxy evolution (594); Early universe (435); High-redshift galaxies (734); Broad band photometry (184); James Webb Space Telescope (2291)*

*Supporting material: machine-readable table*

## 1. Introduction

Some of the oldest stars discovered in the Milky Way indicate a very early onset of star formation in the universe (Tumlinson 2010). Indeed, some inferred stellar ages are older than the Hubble time or very close (Cowan et al. 2002; Creevey

<figure><img src="image_1.png" /></figure>

<footnote n="1">¹Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain</footnote>

<footnote n="2">²DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark</footnote>

<footnote n="3">³Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="4">⁴Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="5">⁵Cosmic Dawn Center (DAWN), Denmark</footnote>

<footnote n="6">⁶DTU Space, Technical University of Denmark, Elektrovej, Building 328, DK-2800, Kgs. Lyngby, Denmark</footnote>

<footnote n="7">⁷Dept. of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="8">⁸Department of Astronomy, Stockholm University, Oscar Klein Centre, AlbaNova University Centre, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="9">⁹UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="10">¹⁰Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote n="11">¹¹I. Physikalisches Institut der Universität zu Köln, Zülpicher Str. 77, D-50937 Köln, Germany</footnote>

<footnote n="12">¹²Telespazio UK for the European Space Agency, ESAC, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Spain</footnote>

<footnote n="13">¹³European Space Agency, Space Telescope Science Institute, Baltimore, MD, USA</footnote>

<footnote n="14">¹⁴School of Physics & Astronomy, Space Research Centre, Space Park Leicester, University of Leicester, 92 Corporation Road, Leicester LE4 5SP, UK</footnote>

<footnote n="15">¹⁵Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="16">¹⁶Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="17">¹⁷Centre for Extragalactic Astronomy, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="18">¹⁸Max Planck Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany</footnote>

<footnote n="19">¹⁹Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191, Gif-sur-Yvette, France</footnote>

<footnote n="20">²⁰Dept. of Astrophysics, University of Vienna, Türkenschanz Str 17, A-1180 Vienna, Austria</footnote>

<footnote n="21">²¹ETH Zürich, Institute for Particle Physics and Astrophysics, Wolfgang-Pauli-Str. 27, 8093 Zürich, Switzerland</footnote>

<footnote n="22">²²School of Cosmic Physics, Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, D02 XF86 Dublin, Ireland</footnote>

<footnote n="23">²³Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven, Belgium</footnote>

<footnote n="24">²⁴Departamento de Física de la Tierra y Astrofísica, Facultad de CC Físicas, Universidad Complutense de Madrid, E-28040, Madrid, Spain</footnote>

<footnote n="25">²⁵Instituto de Física de Partículas y del Cosmos IPARCOS, Facultad de CC Físicas, Universidad Complutense de Madrid, E-28040 Madrid, Spain

*Received 2023 February 5; revised 2023 May 12; accepted 2023 May 22; published 2023 June 27*</footnote>

<footnote>Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
et al. 2015; Schlaufman et al. 2018; Vagnozzi et al. 2022), within 100–200 Myr after (or even before!) the Big Bang. On a complementary approach, (the quite complex) stellar population modeling of stellar clusters and nearby galaxies also indicates very old formation ages for, e.g., the central regions or cores of some elliptical galaxies (Friel 1995; McDermid et al. 2015; Cappellari 2016).

From the point of view of cosmological evolution, one of the most fundamental questions is how fast the universe was able to start forming stars and galaxies. The Hubble Space Telescope (HST) was able to reach redshifts up to $z \sim 11$ (Oesch et al. 2016; Jiang et al. 2021), just 400 Myr after the Big Bang, thus providing statistical samples of high-redshift galaxy candidates from which luminosity functions at $z = 8–10$ could be built (see review Robertson 2022, and references therein; among them, Bouwens et al. 2015; Finkelstein et al. 2015; McLeod et al. 2016; Ishigaki et al. 2018; Yue et al. 2018; Bhatawdekar et al. 2019; Bowler et al. 2020; Rojas-Ruiz et al. 2020; Bagley et al. 2022).

With the launch of JWST on Christmas Day 2021, and the start of scientific operations on 2022 July, several thorough multifilter searches for the highest redshift galaxies have been carried out with the NIRCam instrument, mainly with data sets integrating for less than 2 hr per band. Tens of high-redshift galaxy candidates have been reported at $z = 10–15$, with a few at $15 < z < 18$, and even a handful at $z \sim 20$ (Finkelstein et al. 2022; Naidu et al. 2022a; Bouwens et al. 2023b; Castellano et al. 2023; Harikane et al. 2023; Robertson et al. 2023; Donnan et al. 2023; Finkelstein et al. 2023; Rodighiero et al. 2023; Whitler et al. 2023; Yan et al. 2023). However, different teams analyzing the same observations have only agreed on a few objects (Bouwens et al. 2023a). The spectroscopic confirmation of $z > 10$ candidates is starting to be carried out with JWST/NIRSpec; however, it may require exposure times of tens or even hundreds of hours (Curtis-Lake et al. 2023; Arrabal Haro et al. 2023a, 2023b; or taking advantage of lensing clusters, Roberts-Borsani et al. 2023; Williams et al. 2023), thus bringing into question the feasibility of obtaining robust samples composed of hundreds or even tens of $z > 10$ galaxies, with more promising results at $8 < z < 10$ (Fujimoto et al. 2023; Isobe et al. 2023; Sanders et al. 2023; Tang et al. 2023).

Much more important than the race to discover the highest redshift galaxy or the oldest galaxy ever formed is the exploration of the first hundred million years of galaxy evolution. Constraining the evolution of galaxies within this primordial epoch will provide constraints to simulations that still have significant uncertainties in a broad range of fundamental areas such as (super) massive black hole early growth or primordial black hole abundance, the nature and interplay between dark and baryonic matter in the Dark Ages (Lovell et al. 2014; Dayal et al. 2017), or the mechanism followed by star formation fueled by pristine or low metallicity gas (Zackrisson et al. 2011), affecting properties such as the initial mass function (IMF; Bastian et al. 2010; Marks et al. 2012), and the feedback processes in primitive galaxies (Haslbauer et al. 2022).

In this Letter, we address the topic of identifying and characterizing the galaxy population at $z > 8$ using new ultra-deep NIRCam observations. We observed for a factor 7× longer than previous JWST observations (e.g., in the SMACS0723, GLASS, or CEERS fields), reaching 30–31 mag 5σ limiting

magnitudes and enabling us to probe the high-redshift universe nearly a factor of 10 fainter than the previously analyzed observations in JWST cycle (1) surveys. This allowed us to constrain, for the first time, the faint-end of the ultraviolet (UV) luminosity function from $z = 8$ to $z \sim 13$.

This Letter is organized as follows. Section 2 presents our data set, the NIRCam parallel of the MIRI Deep Survey. Section 3 describes our method to select galaxies at $z > 8$. The sample is used to construct UV luminosity functions and derive the early evolution of the luminosity density in the universe, presented in Section 4. Section 5 summarizes our findings.

Throughout this Letter, we assume a flat cosmology with $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, and a Hubble constant $H_0 = 70 \text{ km s}^{-1}\text{Mpc}^{-1}$. We use AB magnitudes (Oke & Gunn 1983). All stellar mass and star formation rate (SFR) estimations assume a universal Chabrier (2003) IMF.

## 2. Description of the Data

This Letter is based on the analysis of the NIRCam data taken in parallel with the Guaranteed Time Observations of the Hubble Ultra Deep Field (HUDF) carried out in 2022 December by the MIRI European Consortium Team, Program ID 1283 (PIs Hans Ulrik Nørgaard-Nielsen R.I.P., and Göran Östlin). While MIRI was observing the HUDF with the F560W filter during 61 hr (including overheads) in our MIRI Deep Survey (MDS; G. Östlin et al. 2023, in preparation), parallel data were also being acquired by the NIRISS (20 hr) and NIRCam (40 hr) instruments.

The NIRCam observations were carried out with the F115W, F150W, F277W, and F356W filters. Nominal exposure times were 55,187 s in each filter, achieved with 2 visits per filter using 20 integrations with NGGROUP = 7, DEEP8 readout pattern, and 10 dithering positions following the MIRI-optimized CYCLING medium-size pattern with different starting points in each visit. Unfortunately, one of the 4 visits scheduled for our NIRCam observations could not be executed due to several JWST safe-mode events in 2022 December, and one more visit had a significant difference in the roll angle with respect to the first two completed observations ($23^\circ$ versus $31^\circ$). This resulted in a smaller area in common between the four filters: 8.0 arcmin$^2$ instead of the 9.7 arcmin$^2$ expected for a full NIRCam field of view (a further area cut is described in Section 3.1). The aborted observations also resulted in shorter exposure times for the data used in this paper: nominal exposure times of 55,187 s for F115W and F277W, and 27,594 s for F150W and F356W.

The typical 5σ depth of these observations (measured in a $0.\prime\prime3$ diameter circular aperture and corrected for the encircled energy fraction, ~75%) is 30.3 and 30.8 mag in the short and long wavelength channels. This is 2–3 mag deeper than other surveys such as CEERS, GLASS, or SMACS0723 (see Bouwens et al. 2023a).

The NIRCam data were calibrated with the *jwst* pipeline version 1.8.4, reference files in pmap version 1023. Apart from the standard pipeline stages (including snowball and wisp correction), we also applied a background homogenization algorithm prior to obtaining the final mosaics, including 1/f-noise removal. The whole data set was registered to the same World Coordinate System reference frame using the Hubble Legacy Field (HLF) catalog in Whitaker et al. (2019), based on Gaia Data Release 1.2 (Gaia Collaboration et al. 2016a; Gaia Collaboration et al. 2016b). All images were drizzled to the
<table><caption>**Table 1**
Table with Information for All Filters Used in This Work</caption><thead><tr><th>Filter</th><th>Area<br>(arcmin<sup>2</sup>)</th><th>5&sigma; depth<br>(mag)</th><th>Filter</th><th>Area<br>(arcmin<sup>2</sup>)</th><th>5&sigma; Depth<br>(mag)</th></tr></thead><tbody><tr><td>HST/F435W</td><td>4.2</td><td>29.2<sup>29.4</sup><sub>29.0</sub></td><td>HST/F125W</td><td>3.5</td><td>29.9<sup>30.1</sup><sub>28.9</sub></td></tr><tr><td>HST/F606W</td><td>7.2</td><td>29.9<sup>30.2</sup><sub>29.6</sub></td><td>HST/F140W</td><td>2.2</td><td>27.6<sup>27.7</sup><sub>27.5</sub></td></tr><tr><td>HST/F775W</td><td>7.0</td><td>29.5<sup>29.7</sup><sub>29.3</sub></td><td>HST/F160W</td><td>3.4</td><td>29.6<sup>30.1</sup><sub>38.1</sub></td></tr><tr><td>HST/F850LP</td><td>7.2</td><td>29.3<sup>29.5</sup><sub>29.0</sub></td><td>JWST/F115W</td><td>7.3</td><td>30.4<sup>30.6</sup><sub>30.1</sub></td></tr><tr><td>HST/F814W</td><td>7.3</td><td>30.3<sup>30.5</sup><sub>29.9</sub></td><td>JWST/F150W</td><td>7.3</td><td>30.2<sup>30.4</sup><sub>29.9</sub></td></tr><tr><td>HST/F105W</td><td>4.0</td><td>29.7<sup>30.0</sup><sub>28.5</sub></td><td>JWST/F277W</td><td>7.3</td><td>30.8<sup>31.1</sup><sub>30.4</sub></td></tr><tr><td></td><td></td><td></td><td>JWST/F356W</td><td>7.3</td><td>30.8<sup>31.0</sup><sub>30.4</sub></td></tr></tbody></table>

**Note.** We show the area in common with the MIRI European Consortium Guaranteed Time Observations, MIRI Deep Survey NIRCam parallel observations, and the 5σ depth corresponding to a point-like source measured in a 0.′′3 diameter circular aperture and corrected for the limited aperture using empirical PSFs. Median and quartiles are provided for the depth, illustrating the varying depth within the analyzed sky region.

same plate scale, namely, 0.′′03 pixel⁻¹. The FWHM of the point-spread function (PSF) is 0.′′07, 0.′′07, 0.′′12, and 0.′′14 for F115W, F150W, F277W, and F356W, respectively. All images were PSF-matched to the reddest filter.

The MDS NIRCam-parallel observations target a region known as the HUDF parallel field 2 (hereafter HUDF-P2) region, a field where the deepest HST/Advanced Camera for Surveys (ACS) F814W data exist, obtained simultaneously while observing the HUDF with WFC3. The HUDF-P2 field is located to the southeast of the HUDF, at the edge of the CANDELS footprint (Grogin et al. 2011; Koekemoer et al. 2011). A varying fraction of the NIRCam field of view is covered by HST ACS and WFC3 observations, all counting with heterogeneous depths. We used the HLF v2 images for all filters used by ACS and WFC3 (Whitaker et al. 2019). Table 1 presents depths and area coverage for all the data sets used in this paper.

### 3. Selection of $z \ge 8$ Galaxy Candidates

A robust selection of high-redshift galaxy candidates is hampered by three different systematic effects. The first one is the photometry, since we are dealing with intrinsically faint sources, which should be almost or completely undetected in some bands (e.g., those blue wards of the Lyman break or the Lyα jump), thus complicating the spectral energy distribution (SED) analysis. The comparison of results obtained by analyzing photometry in different apertures can potentially help to get rid of noise or contamination by neighboring sources in the determination of the photometric redshifts, as we will explain in Section 3.1. The second difficulty is the degeneracy in the analysis of the SED for the candidates, given that large red colors consistent with the Lyα break can also be mimicked by the Balmer or 4000 Å breaks as well as by the presence of strong (i.e., high equivalent width) emission lines or dusty starbursts (Naidu et al. 2022b; McKinney et al. 2023; Pérez-González et al. 2023; Zavala et al. 2023). The third problem, directly related to the previous one, is the limited and/or biased set of templates that different programs estimating photometric redshifts use. This translates to methodological differences in how the SED degeneracy is addressed. In order to cope with the last two difficulties, we not only used several SEDs for a given source (measured in different apertures), but we also estimated the photometric redshifts with different codes and a priori assumptions

(including the set of templates, but also how to treat low signal-to-noise data points). This is described in Section 3.2.

#### 3.1. Photometric Catalog

We selected galaxies in the F277W filter (the deepest) and measured photometry in all bands following a methodology especially developed to deal with ultra-deep data and crowded fields.

First, we masked the four bright stars present in one of the modules and the very bright spikes coming from two other stars located to the east of our pointing (outside the field of view, but the bright diffraction spikes protrude our observations). This resulted in the loss of 9% of the field of view. The final MDS-NIRCam-par survey area is 7.3 arcmin².

We then detected sources in several passes, starting from the most extended ones and progressing to smaller and fainter galaxies in a number of levels. The galaxies detected in a given level were masked completely down to the isophote corresponding to 5σ of the sky (~26.3 mag arcsec⁻² in F277W), building a first segmentation map of galaxy cores. Afterwards, the segmentation map for each galaxy was extended to the fainter outskirts by fitting isophotes constructed by dilation of the galaxy cores. The outskirts were fitted down to 1σ of the sky (around 28.0 mag arcsec⁻² in F277W) or the Kron (1980) aperture, whichever was reached first. The Kron factor was varied as we were progressing to smaller (and fainter) galaxies, starting at 2.5 for the brightest galaxies and reaching 1.2 for the faintest sources (see Finkelstein et al. 2023). These numbers were calibrated with the sky surface brightness limit mentioned above. The method allows us to select faint galaxies (or globular clusters) in the outskirts of bright extended objects and improve the background determination in the whole field. Once all galaxies in a given level were fitted, the residual image was constructed by removing the core and the outskirts fitted with the isophotal analysis, and the new frame was fed to the next level of detection. We implemented a total of 10 levels for our final catalog in this data set, starting from galaxies covering more than 10,000 pixels (10 arcsec²) in the first pass to the faintest objects with just 30 pixels (~0.03 arcsec²).

Once the procedure was completed for the selection band, the photometry was measured in all other filters by fixing centers, shapes, and sizes. The final catalog is composed of 40,526 sources, of which 33,905 were detected at the 5σ level or higher in F277W. The photometric catalog included color measurements in the Kron (1980) aperture as well as fixed-
diameter circular apertures with sizes 0"2, 0"3, 0"4, and 0"5. The SEDs corresponding to these fixed apertures were scaled to the Kron aperture to obtain total integrated magnitudes by applying a constant offset to all bands in a galaxy-by-galaxy basis.

Photometric errors were estimated by measuring the background noise locally around each source with a procedure similar to that presented in Pérez-González et al. (2008), devised to take into account correlated noise. We remark that estimating realistic uncertainties not affected by drizzling effects (i.e., not underestimated; see, e.g., Labbé et al. 2003; Quadri et al. 2007; Pérez-González et al. 2008; or Finkelstein et al. 2023) and consistently between the short-wavelength (SW) and long-wavelength (LW) channels (which are drizzled in a different way given their distinct nominal pixel sizes, see below) is of utmost importance for our selection and analysis of dropout sources, especially given that the observations for the bluest channels just provide upper limits for the flux.

We built artificial apertures composed of randomly selected noncontiguous pixels, adding up the number of pixels of the photometric aperture whose uncertainty we were trying to estimate. We also imposed that each chosen pixel should be more than three pixels away from any other pixel entering the analysis, this distance being a good approximation of the area that contributes to the flux of each pixel in the final mosaics for any band (after drizzling). With this additional procedure added to a random selection strategy, we ensured that all the pixels included in the noise calculation are independent, thus minimizing the effects of correlated noise. We remark that, in the selection of high-redshift galaxy candidates, we analyzed five different SEDs for each source using the fluxes and noise calculations in the mentioned apertures.

We also calculated photometric errors using randomly distributed circular apertures (as in Harikane et al. 2023, for example). With this method, we find values of the rms of the sky that are 60% and 3% smaller for the long and short wavelength channels, respectively, compared to our method based on randomly selected noncontiguous pixels. We interpret this difference as the effect of correlated noise, which biases the noise calculation in randomly distributed circular apertures because they include contiguous pixels whose signal comes from the same original pixel (i.e., with nominal size). Our random-selection method avoids contiguous pixels, so the effect of correlated noise should be minimized. We note that the F277W and F356W images are drizzled to half the nominal pixel size (65 mas pixel⁻¹), while the nominal pixel size for the F115W and F150W filters (32 mas pixel⁻¹) is very similar to the value used in our final mosaics (30 mas pixel⁻¹). Therefore, we should expect a larger effect of correlated noise in the rms calculations with circular apertures for the long wavelength bands, in agreement with the results of our test.

## 3.2. Photometric Redshift Estimation

The SEDs for each source were fitted with three different codes and techniques, each providing a redshift probability distribution function (zPDF).

Our fiducial photometric redshift estimation came from running the eazy code (Brammer et al. 2008) using the default FSPS templates (Conroy & Gunn 2010), plus a dusty galaxy template (Muzzin et al. 2013). We added the templates presented by Larson et al. (2022) to optimize the analysis of

high-redshift galaxies, which include emission-line galaxies with high equivalent widths (known to contaminate JWST high-redshift samples, see Naidu et al. 2022b; Zavala et al. 2023) and a variety of UV slopes. We used a flat prior in F277W magnitude and no template error, allowing the redshift to take values between z = 0, and z = 20. The data points with low (<3) signal-to-noise ratio (S/N) were treated in two different ways. We either used directly measured fluxes (even if negative) or we also changed those values to 5σ upper limits and used them in a modified version of EAZY. The modification consisted of not allowing any fit to provide brighter fluxes (achieved by penalizing the $\chi^2$ calculation) than those limits and excluding the band in the $\chi^2$ calculation for templates providing lower fluxes (Mérida et al. 2023; Pérez-González et al. 2023).

We also used the prospector SED-fitting algorithm (Johnson et al. 2021) to further verify the estimated photometric redshifts of our candidates. We adopted a setup similar to that used in Langeroodi et al. (2022) for fitting the stellar populations of spectroscopically confirmed $z \sim 8$ galaxies. In brief, the free parameters include the total formed stellar mass, the stellar metallicity, the nebular metallicity, the nebular ionization parameter, the dust-attenuation (modeled with three free parameters, adopting the dust model of Kriek & Conroy 2013), and the optical depth of intergalactic medium (see Langeroodi et al. 2022, for details). Unlike Langeroodi et al. (2022), where the redshift was fixed to the spectroscopically measured value, here we treated redshift as a free parameter with a flat prior in the range $z = 0-20$. The star formation history (SFH) was modeled nonparametrically in five temporal bins, with the last bin spanning 0–10 Myr in lookback time and the rest evenly spaced in log(lookback time) up to $z = 35$. We sampled the parameter space using the dynesty sampler (Speagle 2020; Koposov et al. 2022), a dynamic nested sampler based on Higson et al. (2019).

Finally, we used the template-fitting code LePhare (Arnouts et al. 2002; Ilbert et al. 2006) to help with the selection of our candidates. We adopted two configurations of this program. For the first one, we followed Rinaldi et al. (2022, 2023). Briefly, we made use of galaxy templates with the following set of SFHs: a standard exponentially declining (SFR(t) ∝ exp(-(t-t₀)/τ)) and an instantaneous burst by adopting a simple stellar population model (SFR(t) ∝ δ(t)). In particular, for the standard exponentially declining models, we used the following e-folding timescales (τ) in Gyr: 0.01, 0.1, 0.3, 1, 3, 5, 10, 15. We employed the stellar population synthesis models of Bruzual & Charlot (2003) based on a Chabrier (2003) IMF, considering two different values of the metallicity: solar metallicity and one-fifth of solar metallicity. To account for the effects of internal dust extinction, we convolved the model templates with a modified version of the reddening law of Calzetti et al. (2000), where we adopted the extrapolation provided by Leitherer et al. (2002) at shorter wavelengths. For that purpose, we allowed the color excess $E(B-V)$ to range from 0 to 1.5 mag (with steps of 0.1), reaching $A_V \simeq 6$ mag. A second configuration was explored, this time following the one described in Ilbert et al. (2015), Kauffmann et al. (2022). Since our aim is to select rare sources at $z > 8$, we explore a large range of dust-attenuation in order to identify and reject possible contaminants. We consider two dust-attenuation curves as a free parameter (Calzetti et al. 2000; Arnouts et al. 2013) with $E(B-V)$ varying from 0 to 2 mag
(reaching $A_V = 8$ mag). In this execution, we also adopted the recipes of Saito et al. (2020) to include the emission lines, allowing for an additional free rescaling by a factor two of all line fluxes. We note that LePhare's fits are performed based on the fluxes and not the magnitudes, which do not require the need of introducing upper limits, considering that the flux uncertainties at low signal-to-noise are still meaningful (following Laigle et al. 2016). We rejected sources with $\chi^2 > 100$, which corresponds to an unreliable fit and sources that are best fitted by a brown dwarf stellar template.

Given that we had two executions of *eazy* and two of LePhare, in order to avoid a bias toward any of these codes, we considered the highest redshift among the two executions for each code to select galaxies as explained in the next subsection.

Summarizing, we searched for high-redshift galaxy candidates using five different SEDs and three photometric redshift codes and techniques.

### 3.3. Selection of Candidates

We constructed a sample of $z \gtrsim 8$ galaxy candidates by searching for F115W and F150W dropouts in our NIRCam data and applying a similar methodology to that presented in Finkelstein et al. (2023). Distinctively, we explored the effect of the photometric aperture and the peculiarities of different photometric redshift estimation codes. Indeed, for each galaxy, we used five SEDs derived from different aperture sizes and three redshift estimation codes, as outlined in the previous section, working with weighted averages for all relevant quantities used in the selection. We applied twice the weight of any other SED to that corresponding to 0.′′3 diameter aperture (typically used in this type of studies, see Naidu et al. 2022a; Adams et al. 2023; Bouwens et al. 2023a; Donnan et al. 2023; Harikane et al. 2023), which maximized the S/N of the measurements while keeping PSF-related aperture correction below 25%.

First, we selected all galaxies whose median F277W and F356W flux, averaged across all five SEDs, had a S/N > 5. This cut resulted in a sample of 16,133 sources. Then, we restricted the catalog to all sources with a median S/N < 3 in F115W and/or F150W (in both the PSF-matched and original images), resulting in a sample of 3883 1.15 µm dropouts and 2547 1.50 µm dropouts.

We then selected the galaxies whose most probable (as found by integrating the zPDF) and peak (that providing the smallest $\chi^2$ value) redshifts were above $z > 8$ for any of the 3 codes (1942 sources). We further cut the sample to keep only those galaxies with a >50% cumulative probability of lying at $z > 8$ (1831 sources). The final cut was in $\chi^2$ values, ensuring that the difference between the minimum value of the zPDF, at $z < 8$ and $z > 8$, $\Delta \log \chi^2$ was larger than 0.4 dex (as in Finkelstein et al. 2023; but not as large as in Harikane et al. 2023), leaving 160 sources. The additional cut used in Finkelstein et al. (2023) based on the probability in $\Delta z$ intervals being higher at $z > 8$ than in any other unity-redshift interval did not add anything to our method. All these cuts were based on our fiducial *eazy* run.

We then inspected each of the 160 sources visually in the NIRCam, WFC3, and ACS bands, as well as the SED fits from the different codes. We vetted galaxies that were not affected by artifacts such as (surviving) spikes and/or contamination by nearby objects. We only kept in our final sample the sources for

which at least two out of the three photometric redshift codes and three of the five photometric apertures provided a peak redshift above $z > 8$.

We estimated the fraction of low-redshift interlopers and completeness of our selection method by simulating galaxies with the SEDs provided by the JAGUAR mock catalog (Williams et al. 2018). We added noise to the catalog based on the depths of our survey and applied the selection technique (based on aggregated probabilities, goodness of fits at low and high redshift, and photometric redshift estimation with 3 codes) to the resulting list of sources. This analysis provided a completeness of 80% down to F277W magnitude 31, and a contamination of 3%, all interlopers and missed sources presenting magnitudes fainter than 30.5 mag. No correction from interlopers is applied in the luminosity function calculations presented in the following section.

The final sample of $z > 8$ galaxy candidates consists of 44 objects, listed in Table 2. We show examples in Figures 1 and 2, four dropout sources in the F115W and F150W filters, with the first figure centered on brighter galaxies (F277W < 30 mag) and the second one showing fainter examples (F277W > 30 mag).

We compared our selected $z > 8$ candidate galaxies with those reported in Austin et al. (2023) for the common area with the NGDEEP survey (Bagley et al. 2023). Only 2 sources in Austin et al. (2023) lie within our field of view; another 3 (2) are located in the edges (gap between detectors) of our imaging data and only covered by 2 (0) bands (i.e., they could not arrive to our selection in any case). Of the 2 sources in common, our MDS011049 ($z \sim 9.5$ in this Letter) is NGD-z11a ($z \sim 11.1$ for Austin et al. 2023). We assign it with a lower redshift mostly because it is well detected in F150W (although near our image edge). The other common source is NGD-z11b, MDS007481 in our parent catalog; it did not enter our final selection because we detected two redshift solutions, one at $z \sim 6.5 \pm 0.2$, and another one at $z = 11.1 \pm 0.5$; the goodness of fit difference between them lied below our cut, and the low-$z$ solution was favored by the photometry measured in some of our aperture.

Figure 3 presents the distribution of F277W (Kron-aperture) magnitudes and photometric redshifts of the final sample, compared to general catalogs and other similar compilations of $z > 8$ candidates and confirmed galaxies in the first surveys carried out by JWST in 2022. The figure shows that our MDS-NIRCam-par data reach 2 mag fainter than the surveys of SMACS0723, GLASS, and CEERS, with a histogram peaking at F277W ~ 30 mag, compared to 28 mag for the first data sets. The 2022 shallow surveys add up ~8 times the area of our observations, but we increase by almost a factor of 10 the number of detected sources per unit magnitude bin down to 31 mag.

Our deeper data allow us to find $z > 8$ galaxy candidates at fainter magnitudes, probing lower luminosities than previous works (see Section 4). Our sample, plotted in red in Figure 3, is compared with other photometrically selected candidates and the four first spectroscopic confirmations at $z > 10$ provided by JWST (Curtis-Lake et al. 2023; Arrabal Haro et al. 2023a, 2023b) as well as faint lensed galaxies above $z \sim 7$ confirmed with NIRSpec spectroscopy (Langeroodi et al. 2022; Carnall et al. 2023; Roberts-Borsani et al. 2023; Williams et al. 2023). We have seven candidates in the magnitude regime probed in the 2022 shallow surveys due to our limited area, with the bulk of our sample concentrating around 30 mag (median and quartiles F277W = 30.2$^{30.7}_{29.8}$ mag). Our sample doubles the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Examples of F115W and F150W dropout galaxies in our sample. For each source, in the top left, we show the SED (for the 0.1"3 diameter aperture; filled dots) and template fits to low- and high-redshift solutions. Arrows indicate 1σ upper limits. The top right panel shows the photometric redshift probability distribution function for the different codes used in this work and the derived redshift and uncertainties. The bottom of the figure for each source presents 1.1"5 × 1.1"5 postage stamps in ACS, WFC3, and NIRCam bands, with the source marked with a 0.1"15 radius red circle.</figcaption></figure>

number of $z \sim 9$ galaxy candidates and covers an unexplored magnitude regime up to $z \sim 13$.

We note that our data are well suited to get $z \sim 9$ galaxy candidates with well-covered SEDs to identify the break and obtain photometric redshift with relatively small (formal) uncertainties. However, at $z \gtrsim 13$, we cannot constrain

the position of the break since we do not have the F200W filter. This translates also to a higher uncertainty in the photometric redshift estimates at $z \gtrsim 10$, when the breaks start to enter the F150W filter, and we only have the F277W –F356W color to constrain the redshift further. This is the explanation for the larger error bars at those redshifts
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Examples of F115W and F150W dropout galaxies in our sample; in this case, for sources fainter than F277W = 30 mag (distinctively probed by our data set). Same information as in Figure 1 given.</figcaption></figure>

(especially for the $z \sim 19$ candidate, shown in Figure 2) and the gap at $z \sim 10$, even when the galaxies are not significantly fainter. This was tested and confirmed with JAGUAR simulations similar to those mentioned in Section 3. We compared results using our NIRCam filter set with those achieved when adding F200W fluxes (with depths similar to

the F115W filter). The inclusion of the F200W filter decreased the typical photometric redshift errors by a factor of >3 and also provided larger completeness levels, reaching values beyond 90% (compared to 80% without that filter); the interloper fraction decreased from 3% to 1% when adding F200W to the simulations.
<figure><img src="image_4.png" /></figure>

**Figure 3.** Magnitude in the NIRCam F277W filter vs. photometric redshift for $z > 8$ galaxies. General galaxy samples extracted from the catalogs released by G. Brammer for the SMACS0723, GLASS, and CEERS data sets gathered in 2022 are plotted with blue small dots to illustrate the magnitude limits of the surveys and the number of high-redshift galaxies in unsupervised and unvetted photometric catalogs. Our parent galaxy sample is plotted with red small dots, and the final sample of $z > 8$ galaxy candidates with red hexagons. We depict $1\sigma$ uncertainties plotted in both axes and the data point referring to the zPDF peak redshift. We also plot the photometric high-redshift galaxy candidates of Bouwens et al. (2023a), Naidu et al. (2022a), Harikane et al. (2023), Donnan et al. (2023), and Finkelstein et al. (2023) with large black circles (we do not remove repeated sources with different measurements), as well as spectroscopically confirmed galaxies from JADES (Curtis-Lake et al. 2023; Robertson et al. 2023), lensing clusters (Langeroodi et al. 2022; Roberts-Borsani et al. 2023; Williams et al. 2023; Carnall et al. 2023, not corrected for magnification), and CEERS sources (Arrabal Haro et al. 2023a, 2023b), all in green. The top panel shows the histogram of redshifts for these samples of high-redshift galaxy candidates. The right panel shows the magnitude distribution of the different data sets mentioned before, with the three 2022 surveys plotted with transparency. Our deeper data (peaking around 2 mag fainter than previous surveys) probe the 29–31 mag regime, providing candidates at the previously unexplored faint-end of the luminosity functions.

# 4. Results

## 4.1. UV Luminosity Functions

We divided the sample in 3 redshift bins, namely, $8 < z < 10$, $10 < z < 11.5$, and $11.5 < z < 13$, which correspond to universe ages ~540, ~420, and ~350 Myr. We calculated absolute UV magnitudes ($M_{UV}$, averaged in a $0.01~\mu$m window around $0.15~\mu$m) from SED fitting, obtaining uncertainties by repeating the stellar population modeling after varying the photometry according to their errors and considering the redshift probability distribution functions.

The luminosity functions were constructed using a $V_{max}$ formalism and the stepwise maximum likelihood method (Efstathiou et al. 1988). The completeness was estimated by inserting artificial sources extracted from the real data in the detection image, covering the magnitude range between 26 and 32 mag, and then repeating the multilayer source detection presented in Section 3.1. We find that our catalog is 80% complete at F277W = 30.4 mag, 50% complete at F277W = 30.6 mag, and 10% complete at F277W = 31.0 mag.

Using the implementation to determine cosmic variance effects presented by Trapp & Furlanetto (2020), consistent with the predictions from Bhowmick et al. (2020), Ucci et al. (2021), we estimate a ~25% (18%) uncertainty in our luminosity function estimation at the bright (faint) end for $z \sim 9$, the values increasing to 34% (28%) at the highest redshifts. These uncertainties, similar to those presented in Trenti & Stiavelli (2008), are included in the following figures and discussion.

The luminosity functions are presented in Figure 4 for the three redshift bins mentioned above, comparing with previous estimates in the literature. Our results are given in Table 3. The typical systematic offsets and scatter of previous calculations are as high as 0.4 dex, with only a minor degradation as we move to higher redshifts, and some larger differences (0.6 dex level) at the bright end, $M_{UV} < -21$ mag, at $z \sim 9$ and $z \sim 12$.

We fit all these estimates (our data and other data points) with Schechter (1976) functions without any prior using a Markov Chain Monte Carlo method. The results are given in Table 4.

We note that several of the luminosity function data points are based on the same data sets, although not on the same sources, and the methodologies to select candidates are different. This might affect the uncertainties in the fits.

To test this possibility, we repeated the fits at $z \sim 9$ only including our data points and those in Donnan et al. (2023), Harikane et al. (2023), which cover all the public JWST fields (including Stephan's quintet, not present in Bouwens et al. 2023a) as well as the COSMOS field for the bright end. We found completely consistent results and slightly larger uncertainties: $(\alpha, M^*, \phi^*) = (-2.30_{-0.24}^{+0.16}, -21.00_{-0.45}^{+0.34}, 1.8_{-1.7}^{+1.1})$ (normalization in units of $10^{-5}$ Mpc$^{-3}$ mag$^{-1}$) for the whole data set, compared to $(-2.36_{-0.25}^{+0.17}, -20.96_{-0.64}^{+0.40}, 2.0_{-2.3}^{+1.2})$ for the limited set.

For $z \sim 11$, we obtained best-fitting values $(-2.14_{-0.38}^{+0.24}, -20.74_{-0.54}^{+0.57}, 3.9_{-3.1}^{+4.0})$ for all data points (same units as above), and $(-2.39_{-0.22}^{+0.16}, -20.24_{-0.29}^{+0.45}, 5.3_{-2.5}^{+4.1})$ only using Donnan et al. (2023) and our data points.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** From left to right, UV luminosity functions at $z \sim 9$, $z \sim 11$ and $z \sim 12$. We plot our results with red hexagons (uncertainties smaller than data points in some cases, see Table 3), and compare with literature estimations from Bouwens et al. (2023a), Naidu et al. (2022a), Donnan et al. (2023), Finkelstein et al. (2023), Harikane et al. (2023), and Castellano et al. (2023). We fit all points (ours and those in the literature) to a Schechter (1976) function, which is plotted with red continuous lines in each panel, with the fits for other redshifts shown with dashed lines.</figcaption></figure>

Finally, for $z \sim 12$, we obtained ($-2.19^{+0.26}_{-0.39}$, $-20.81^{+0.77}_{-0.67}$, $1.6^{+1.6}_{-0.9}$) using all estimations (same units as above), and ($-2.05^{+0.40}_{-0.46}$, $-20.25^{+0.95}_{-0.82}$, $2.2^{+2.1}_{-1.3}$) when only using Harikane et al. (2023) and our data points.

Overall, we find little evolution between $z \sim 13$ and $z \sim 8$; the fits run roughly in parallel. The faint-end slope is consistent across all redshifts down to $M_{UV} \sim -17$ mag, with a value $\alpha = -2.2 \pm 0.1$ and no indications of a steepening (in contrast with previous findings, e.g., Bouwens et al. 2023a). We also obtain very similar knee absolute magnitudes, consistent within uncertainties, also considering the degeneracy with the other two Schechter parameters. The main difference between the three redshift bins is a small average number density evolution, increasing by 0.2–0.3 dex from $z \sim 13$ to $z \sim 8$ (between the last and first bins considered in this Letter), i.e., in $\sim 200$ Myr, although this difference is of the same order as the current uncertainties for individual luminosity function estimations.

## 4.2. The Cosmic UV Luminosity and SFR Density

We integrated the luminosity functions presented in Figure 4 down to the same absolute magnitude for all redshift bins, $M_{UV} \sim -17$ mag, obtaining the (directly observed) UV luminosity densities at $8 < z < 13$. In Figure 5, we compare our results with estimates from the literature. This figure includes pre-JWST era measurements (Oesch et al. 2018), as well as more recent estimates obtained with JWST data (Bouwens et al. 2023a; Donnan et al. 2023; Harikane et al. 2023).

Figure 5 also shows theoretical predictions from a range of state-of-the-art models in the literature covering a variety of approaches. The simulation compilation includes the UniverseMachine semiempirical model (Behroozi et al. 2020), the Santa Cruz semianalytical models (Yung et al. 2019, 2020), and the fiducial cold dark matter model in Maio & Viel (2023). We also compare with the hydrodynamical simulations Illustris/TNG (Vogelsberger et al. 2014; Nelson et al. 2015; Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018; Nelson et al. 2019; Pillepich et al. 2019), EAGLE (Crain et al. 2015; Schaye et al. 2015), and SIMBA (Davé et al. 2019). We also compare with predictions from resimulations of multiple regions, such as THESAN (two realizations are shown,

THESAN-1 and THESAN-SDAO-2, both based on IllustrisTNG galaxy formation model; Kannan et al. 2022) and the First Light And Reionisation Epoch Simulations (FLARES, which adopt the AGNdT9 variant of EAGLE; Lovell et al. 2021; Vijayan et al. 2021). All these simulations probe a large box size, spanning from 50 cMpc in TNG50 to 250 cMpc in UniverseMachine and 500 cMpc for FLARES. These models have been validated by matching to observations at lower redshifts than the epochs probed in this Letter.

For this comparison with simulations, we directly took UV luminosity functions from the mentioned papers or added the SFRs (typically averaged over 100 Myr periods) of all or certain subsamples (see below) of simulated galaxies in given simulation snapshot, which can be converted to a luminosity density assuming a constant star formation event following a given IMF, Chabrier (2003) in our case.

We note that most of these simulations provide properties integrated in galaxy regions that are larger than the actual measurements from our NIRCam observations. Indeed, the median and quartile values of the Kron apertures for our (slightly resolved) sample are $0.′′46_{0.36}^{0.58}$ (in diameter), i.e., observations typically refer to regions of $0.′′5$, which translates to $\sim 2$ kpc physical size (2.4 kpc for $z=8$, 1.7 kpc for $z=13$). This means that some of the vertical offsets seen between our and literature's measurements and the model predictions can be interpreted in terms of aperture corrections.

We exemplify this issue for the TNG100 simulation, shown with blue large-dot lines for whole-galaxy SFRs. The galaxies in the TNG100 simulation at these redshifts present typical half-mass radii around 3 kpc, 6 kpc diameter, 3 times the typical photometric aperture sizes we use in this work. In Figure 5, we also show the predictions for $0.′′5$ apertures (small dots). While the former curve (large dots, whole-galaxy) runs closer to the observational data points, the latter (small dots, limited aperture), which is more directly comparable to the data, is consistent with observations at $z \sim 8$ but predicts 0.7 dex smaller luminosity densities at $z \sim 13$. This means that the star formation is occurring in more compact regions than what these simulations predict. The importance of galaxy size and morphology in JWST surveys is indeed paramount to understanding galaxy evolution in the early universe (Costantin et al. 2023).
<table><caption>**Table 2**
Sample of $z > 8$ Galaxy Candidates Presented in This Letter</caption><thead><tr><th>ID</th><th>R.A. (J2000)<br>(degrees)</th><th>Decl. (J2000)<br>(degrees)</th><th>F277W<br>(mag)</th><th>Redshift</th><th>No. Codes z > 8</th></tr></thead><tbody><tr><td>MDS025593</td><td>53.23575778</td><td>-27.84791165</td><td>30.25 ± 0.13</td><td>19.6<sup>+4.2</sup><sub>-8.1</sub></td><td>2</td></tr><tr><td>MDS007258</td><td>53.23712648</td><td>-27.85429400</td><td>29.89 ± 0.07</td><td>9.0<sup>+1.6</sup><sub>-1.6</sub></td><td>2</td></tr><tr><td>MDS013931</td><td>53.23809134</td><td>-27.86844857</td><td>30.18 ± 0.14</td><td>11.5<sup>+7.0</sup><sub>-1.8</sub></td><td>2</td></tr><tr><td>MDS007507</td><td>53.24247528</td><td>-27.85509211</td><td>30.24 ± 0.08</td><td>9.0<sup>+1.6</sup><sub>-1.2</sub></td><td>2</td></tr><tr><td>MDS021337</td><td>53.24305327</td><td>-27.85711165</td><td>30.21 ± 0.12</td><td>8.9<sup>+0.7</sup><sub>-3.0</sub></td><td>3</td></tr><tr><td>MDS006697</td><td>53.24608867</td><td>-27.86033795</td><td>29.87 ± 0.14</td><td>10.6<sup>+1.5</sup><sub>-3.2</sub></td><td>3</td></tr><tr><td>MDS010233</td><td>53.24634046</td><td>-27.84432908</td><td>27.79 ± 0.05</td><td>9.2<sup>+2.0</sup><sub>-2.6</sub></td><td>3</td></tr><tr><td>MDS009900</td><td>53.24674102</td><td>-27.84563609</td><td>29.03 ± 0.04</td><td>9.3<sup>+0.6</sup><sub>-4.8</sub></td><td>3</td></tr><tr><td>MDS028153</td><td>53.24926367</td><td>-27.84894921</td><td>30.60 ± 0.12</td><td>9.3<sup>+2.8</sup><sub>-2.2</sub></td><td>2</td></tr><tr><td>MDS008711</td><td>53.24965269</td><td>-27.85225902</td><td>29.60 ± 0.06</td><td>10.3<sup>+3.6</sup><sub>-4.5</sub></td><td>3</td></tr><tr><td>MDS015081</td><td>53.25034986</td><td>-27.87040839</td><td>29.95 ± 0.08</td><td>8.7<sup>+1.8</sup><sub>-2.1</sub></td><td>3</td></tr><tr><td>MDS022582</td><td>53.25129541</td><td>-27.85826615</td><td>30.57 ± 0.15</td><td>9.0<sup>+1.8</sup><sub>-5.2</sub></td><td>2</td></tr><tr><td>MDS018975</td><td>53.25346260</td><td>-27.86461612</td><td>29.34 ± 0.09</td><td>12.2<sup>+3.1</sup><sub>-2.6</sub></td><td>2</td></tr><tr><td>MDS005128</td><td>53.25431438</td><td>-27.87161787</td><td>29.88 ± 0.07</td><td>9.1<sup>+1.9</sup><sub>-1.6</sub></td><td>2</td></tr><tr><td>MDS023789</td><td>53.25460386</td><td>-27.85749411</td><td>29.39 ± 0.07</td><td>9.3<sup>+1.0</sup><sub>-0.6</sub></td><td>2</td></tr><tr><td>MDS028519</td><td>53.25468485</td><td>-27.85043609</td><td>30.62 ± 0.12</td><td>8.3<sup>+0.8</sup><sub>-2.4</sub></td><td>3</td></tr><tr><td>MDS026779</td><td>53.25653543</td><td>-27.85376201</td><td>29.59 ± 0.08</td><td>11.9<sup>+1.8</sup><sub>-3.8</sub></td><td>2</td></tr><tr><td>MDS014361</td><td>53.25797759</td><td>-27.87485511</td><td>30.14 ± 0.21</td><td>9.0<sup>+3.6</sup><sub>-0.5</sub></td><td>2</td></tr><tr><td>MDS029633</td><td>53.25929715</td><td>-27.85077767</td><td>29.52 ± 0.08</td><td>8.5<sup>+2.1</sup><sub>-4.1</sub></td><td>2</td></tr><tr><td>MDS030229</td><td>53.25954821</td><td>-27.85017797</td><td>30.62 ± 0.11</td><td>9.3<sup>+1.0</sup><sub>-3.2</sub></td><td>2</td></tr><tr><td>MDS027948</td><td>53.25957616</td><td>-27.85306371</td><td>29.62 ± 0.12</td><td>9.1<sup>+1.1</sup><sub>-1.3</sub></td><td>2</td></tr><tr><td>MDS020295</td><td>53.26008969</td><td>-27.86493963</td><td>30.87 ± 0.18</td><td>8.8<sup>+2.4</sup><sub>-3.6</sub></td><td>2</td></tr><tr><td>MDS018332</td><td>53.26100195</td><td>-27.86847883</td><td>30.33 ± 0.17</td><td>10.8<sup>+2.4</sup><sub>-4.7</sub></td><td>2</td></tr><tr><td>MDS017690</td><td>53.26517306</td><td>-27.87116585</td><td>30.01 ± 0.11</td><td>8.3<sup>+1.3</sup><sub>-5.2</sub></td><td>2</td></tr><tr><td>MDS011049</td><td>53.26722949</td><td>-27.84890925</td><td>27.66 ± 0.04</td><td>9.4<sup>+0.1</sup><sub>-0.2</sub></td><td>2</td></tr><tr><td>MDS006210</td><td>53.28683244</td><td>-27.87802853</td><td>29.63 ± 0.09</td><td>9.1<sup>+2.3</sup><sub>-4.2</sub></td><td>3</td></tr><tr><td>MDS021311</td><td>53.28746863</td><td>-27.87372630</td><td>31.12 ± 0.22</td><td>8.4<sup>+1.7</sup><sub>-1.6</sub></td><td>3</td></tr><tr><td>MDS022349</td><td>53.29153766</td><td>-27.87366009</td><td>30.36 ± 0.15</td><td>8.8<sup>+1.3</sup><sub>-2.3</sub></td><td>3</td></tr><tr><td>MDS008116</td><td>53.29301264</td><td>-27.87124117</td><td>28.91 ± 0.11</td><td>12.6<sup>+2.4</sup><sub>-2.6</sub></td><td>2</td></tr><tr><td>MDS004915</td><td>53.29351197</td><td>-27.88750272</td><td>29.39 ± 0.10</td><td>10.7<sup>+3.0</sup><sub>-3.5</sub></td><td>3</td></tr><tr><td>MDS020574</td><td>53.29521802</td><td>-27.87761048</td><td>30.81 ± 0.18</td><td>10.7<sup>+2.7</sup><sub>-3.5</sub></td><td>3</td></tr><tr><td>MDS005199</td><td>53.29802771</td><td>-27.88750404</td><td>29.52 ± 0.08</td><td>8.7<sup>+1.5</sup><sub>-1.9</sub></td><td>3</td></tr><tr><td>MDS005247</td><td>53.29901733</td><td>-27.88760795</td><td>29.75 ± 0.08</td><td>9.1<sup>+0.7</sup><sub>-1.4</sub></td><td>2</td></tr><tr><td>MDS006765</td><td>53.30019902</td><td>-27.88012852</td><td>29.29 ± 0.10</td><td>9.0<sup>+1.7</sup><sub>-1.8</sub></td><td>3</td></tr><tr><td>MDS018069</td><td>53.30742315</td><td>-27.88628994</td><td>30.78 ± 0.14</td><td>8.2<sup>+2.1</sup><sub>-3.5</sub></td><td>2</td></tr><tr><td>MDS007889</td><td>53.30748690</td><td>-27.87763847</td><td>29.84 ± 0.08</td><td>8.8<sup>+2.7</sup><sub>-1.6</sub></td><td>2</td></tr><tr><td>MDS008635</td><td>53.30804862</td><td>-27.87443982</td><td>29.98 ± 0.13</td><td>8.9<sup>+1.9</sup><sub>-4.7</sub></td><td>2</td></tr><tr><td>MDS008261</td><td>53.30820570</td><td>-27.87625163</td><td>28.73 ± 0.06</td><td>9.2<sup>+2.5</sup><sub>-4.0</sub></td><td>3</td></tr><tr><td>MDS006079</td><td>53.30921148</td><td>-27.88711824</td><td>28.91 ± 0.05</td><td>10.7<sup>+3.1</sup><sub>-3.3</sub></td><td>3</td></tr><tr><td>MDS025580</td><td>53.31027610</td><td>-27.87573262</td><td>29.17 ± 0.10</td><td>10.1<sup>+2.6</sup><sub>-3.8</sub></td><td>3</td></tr><tr><td>MDS005234</td><td>53.31030222</td><td>-27.89189564</td><td>29.27 ± 0.08</td><td>12.4<sup>+8.6</sup><sub>-4.0</sub></td><td>2</td></tr><tr><td>MDS025608</td><td>53.31030283</td><td>-27.87569216</td><td>30.10 ± 0.23</td><td>10.7<sup>+3.3</sup><sub>-3.5</sub></td><td>2</td></tr><tr><td>MDS006529</td><td>53.31597314</td><td>-27.88737382</td><td>30.01 ± 0.09</td><td>8.9<sup>+0.9</sup><sub>-6.0</sub></td><td>3</td></tr><tr><td>MDS025422</td><td>53.31811452</td><td>-27.87887496</td><td>30.54 ± 0.19</td><td>11.3<sup>+4.0</sup><sub>-3.8</sub></td><td>2</td></tr></tbody></table>

Note. We provide coordinates, integrated magnitudes in the F277W filter, and photometric redshifts estimated with EAZY. The comment gives information about how many of the 3 photometric redshift codes used in this Letter agree in the $z > 8$ determination.

(This table is available in machine-readable form.)

Another point to consider is the mass range covered by observations and simulations. The typical stellar masses (assuming standard IMF and SFHs recipes) for our sample are $\log(M_*/M_\odot) = 7.8_{-0.2}^{+6.6}$ (median and quartiles, P. G. Pérez-González et al. 2023, in preparation). For some models, such as EAGLE (in green), UniverseMachine (in orange), and Illustris/TNG (comparing TNG300 versus other resolutions), we are able to distinguish between masses. These 3 sets of

models show differences of 0.2–0.6 dex when considering a cut in $M_*/M_\odot$ instead of $M_*/M_\odot$, and they predict 0.4–0.8 dex smaller densities when cutting at $M_*/M_\odot > 10^8 M_\odot$, these curves being more comparable to our sample.

Disregarding the difficulty in comparing with models due to aperture and mass effects (not always traceable, since some simulations do not count with the adequate resolution and/or information), most simulations predict an increase of
<figure><img src="image_6.png" /></figure>

**Figure 5.** Evolution of the UV luminosity density, transformed into SFR density on the right vertical axis (assuming a Chabrier 2003 IMF). Our results are plotted with a red shaded region. Estimates from the literature are plotted with gray symbols (Oesch et al. 2018; Bouwens et al. 2023a, 2023b; Donnan et al. 2023; Harikane et al. 2023). The densities at $z > 16$ are shown in light gray since they were mostly based on a $z \sim 5$ redshift interloper (Arrabal Haro et al. 2023a; Harikane et al. 2023). We also plot predictions from several galaxy formation simulations, listed in the following lines. The Illustris-1, TNG50, TNG100, TNG300 results are plotted with blue continuous, dashed, dotted, and dashed-dotted lines, respectively, all referring to SFRs measured in $0''/5$ diameter apertures for $M_* > 10^7 M_\odot$ galaxies, except for TNG300, cut at $M_* > 10^8 M_\odot$. The thick dotted line refers to TNG100 results using whole-galaxy measurements (Vogelsberger et al. 2014; Nelson et al. 2015; Naiman et al. 2018; Marinacci et al. 2018; Nelson et al. 2018; Pillepich et al. 2018, 2019; Springel et al. 2018; Nelson et al. 2019). The EAGLE simulations (Schaye et al. 2015; Crain et al. 2015) are plotted in green for galaxy samples cut in stellar masses $M_* > 10^{6.7,8} M_\odot$ in continuous, dotted, dashed lines. The SIMBA predictions are shown with a brown continuous line for all galaxies, dashed line for F277W < 31 mag sources (Davé et al. 2019). The THESAN-1 and THESAN-SDAO-2 simulations (Kannan et al. 2022) are shown in magenta with continuous and dotted lines, respectively (see text for differences). The FLARES predictions are shown in cyan (Lovell et al. 2021; Vijayan et al. 2021). The UniverseMachine semiempirical model (Behroozi et al. 2020) is shown in orange with the same mass cuts and line styles mentioned for EAGLE. Finally, the models in Maio & Viel (2023) are shown in yellow, and the predictions from the Santa Cruz semianalytic models (Yung et al. 2019, 2020) are depicted in orchid.

<table><caption>**Table 3**
Luminosity Function Data Points Obtained in This Work</caption><thead><tr><th>M<sub>UV</sub><br>AB mag</th><th>log φ<br>Mpc<sup>-2</sup> mag<sup>-1</sup></th><th>M<sub>UV</sub><br>AB mag</th><th>log φ<br>Mpc<sup>-2</sup> mag<sup>-1</sup></th><th>M<sub>UV</sub><br>AB mag</th><th>log φ<br>Mpc<sup>-2</sup> mag<sup>-1</sup></th></tr></thead><tbody><tr><td>8 &lt; z &lt; 10</td><td></td><td>10.0 &lt; z &lt; 11.5</td><td></td><td>11.5 &lt; z &lt; 13.0</td><td></td></tr><tr><td>-19.5</td><td>-4.17<sup>+0.26</sup><sub>-0.61</sub></td><td>-18.5</td><td>-3.57<sup>+0.21</sup><sub>-0.36</sub></td><td>-18.3</td><td>-3.47<sup>+0.26</sup><sub>-0.43</sub></td></tr><tr><td>-18.5</td><td>-3.53<sup>+0.15</sup><sub>-0.23</sub></td><td>-17.5</td><td>-2.82<sup>+0.21</sup><sub>-0.40</sub></td><td>-17.3</td><td>-2.87<sup>+0.28</sup><sub>-0.53</sub></td></tr><tr><td>-17.5</td><td>-2.77<sup>+0.11</sup><sub>-0.17</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>-16.5</td><td>-2.30<sup>+0.16</sup><sub>-0.27</sub></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 4**
Results for the Schechter (1976) Parameterization Fits to the Luminosity Functions Presented in Figure 4</caption><thead><tr><th>Parameter</th><th>8 &lt; z &lt; 10</th><th>10.0 &lt; z &lt; 11.5</th><th>11.5 &lt; z &lt; 13.0</th></tr></thead><tbody><tr><td>α</td><td>-2.30<sup>+0.16</sup><sub>-0.24</sub></td><td>-2.14<sup>+0.24</sup><sub>-0.38</sub></td><td>-2.19<sup>+0.26</sup><sub>-0.39</sub></td></tr><tr><td>M* [AB mag]</td><td>-21.00<sup>+0.34</sup><sub>-0.45</sub></td><td>-20.74<sup>+0.57</sup><sub>-0.54</sub></td><td>-20.81<sup>+0.77</sup><sub>-0.67</sub></td></tr><tr><td>φ* [10<sup>-5</sup> Mpc<sup>-3</sup> mag<sup>-1</sup>]</td><td>1.8<sup>+1.1</sup><sub>-1.7</sub></td><td>3.9<sup>+4.0</sup><sub>-3.1</sub></td><td>1.6<sup>+1.6</sup><sub>-0.9</sub></td></tr><tr><td>log ρ<sub>UV</sub> [10<sup>25</sup> erg s<sup>-1</sup> Hz<sup>-1</sup> Mpc<sup>-3</sup>]</td><td>1.51<sup>+1.56</sup><sub>-0.39</sub></td><td>1.34<sup>+0.45</sup><sub>-0.29</sub></td><td>0.50<sup>+0.21</sup><sub>-0.18</sub></td></tr></tbody></table>

**Note.** The last row shows the integrated luminosity for absolute magnitudes $M_{UV} < -17$ mag.
luminosity density of almost 2 orders of magnitude in the 300 Myr from $z \sim 13$ to $z \sim 8$. The observations differ from this increase rate. We note that our JWST-based results at these redshifts present a striking agreement considering the different approaches used to identify high-redshift galaxies. All these estimates together, jointly with the first results of the UV luminosity density beyond $z=13$ shown in Figure 5, suggest a shallower increase or a more constant behavior of the earliest phases of galaxy formation from $z=18$ until $z \sim 9$, followed by a more pronounced increase of star formation activity at $z \lesssim 9$. Quantitatively, the luminosity density evolution follows a $(1+z)^{-4.5 \pm 1.0}$ law according to the data at $8 < z < 14$, $(1+z)^{-2.8 \pm 0.9}$ at $8 < z < 18$, while most models predict steeper evolution, with an average slope around $-10$, ranging from $-20$ for SIMBA or $-15$ for EAGLE to $-8$ for Illustris/TNG.

Overall, most simulations present tension with these first measurements provided by JWST, as they typically under-predict the abundance of high-redshift luminous galaxies and/or their SFHs, stellar masses, and compactness. Among the simulations, the results of THESAN seem in agreement with observed log($\rho_{SFR}$) values in the redshift range $10 < z < 13$ (although with no information about the spatial extension of the stellar cores we are probing with JWST), but with faster evolution (compare the steeper magenta line compared to red region). In this case, however, there is a significant difference depending on the reionization histories. In THESAN-1 (continuous magenta line), the duration of reionization is short, where low mass halos have a significant contribution. In contrast, the duration is longer in THESAN-SDAO-2 (dotted magenta line, closer to the data), since it assumes nonstandard dark matter models, with high mass halos being the primary drivers of reionization. Furthermore, at $z > 8$, the total SFR density in THESAN-1 is dominated by the lowest mass halos ($10^8 < M_{halo}/M_{\odot} < 10^9$), which are usually unresolved in the medium-resolution simulations, thus reducing the SFR density by almost an order of magnitude.

From the observational perspective, the lack of spectroscopic confirmation (still quite demanding and/or unfeasible for a statistically significant sample) still hampers our physical interpretation of the early stage in galaxy formation (see Harikane et al. 2023). The spectroscopy is needed not only to confirm candidates but also to train our photometric redshift algorithms, necessary to probe the faint-end of the luminosity function and to improve statistics and control cosmic variance effects. On the theoretical side, various possibilities are on the table, like a variable (temperature-dependent) stellar IMF (Steinhardt et al. 2020a, 2022b), early dark energy models (Smith et al. 2022), dark matter properties (Dayal et al. 2017), or star formation efficiencies ~15%–30% higher than expected (Inayoshi et al. 2022).

## 5. Summary and Conclusions

We present a sample of 44 $z > 8$ galaxy candidates identified in the NIRCam parallel observations of the MDS (centered in the HUDF and parallel pointing 2, respectively). These data reach around 2 mag deeper than those from previous surveys carried out during Early Release Observations, Early Release Science, and cycle (1) programs (e.g., the SMACS0723 ERO survey, CEERS, and GLASS), reaching $5\sigma$ depths around 30.8 mag. The median and quartiles of the magnitude distribution of our sample is $\langle F_{277W} \rangle = 30.2_{29.8}^{30.7}$ mag.

Aiming to reduce the effects of photometric uncertainties and a priori assumptions in the determination of redshifts, the selection of our $z > 8$ galaxy candidates is based on the analysis of SEDs measured in a variety of apertures and with different photometric redshift codes. We also apply restringent cuts in the goodness of fit and the relevance of high- compared to low-redshift solutions, following standard procedures in the literature.

The sample probes the absolute magnitude regime $-19.5 < M_{UV} < -16.5$ mag of the luminosity function at $8 < z < 13$, corresponding to universe ages between 350 and 540 Myr. We constrain the faint-end, obtaining a constant slope of $-2.2 \pm 0.1$ in the full redshift range. Jointly with previous results based on shallower larger-area surveys probing the bright-end of the luminosity function, we find that the regime we explore accounts for nearly 50% of the total luminosity density (integrated through all magnitudes brighter than $-17$ mag).

Our estimates of the integrated UV luminosity density, which is a good proxy for the cosmic SFR density (via assumptions in relevant star formation properties such as the IMF, the metallicity, binary star fraction, or the SFH or burstiness behavior) are compared with a wide range of predictions from state-of-the-art galaxy evolution models. The two main results are as follows: (1) we find a shallower slope in the early stages in galaxy formation at $8 < z < 13$, where many models predict nearly a factor of 10 less star formation activity; (2) although some models predict similar values of the cosmic SFR density to those measured in this paper and in recent publications based on JWST data, when taking into account the sizes (as measured by the aperture diameters used in the photometric extractions) and the stellar masses corresponding to the apparent magnitudes (assuming typical stellar mass calculation recipes and a standard universal IMF), the simulations typically identify the location of star formation in larger, less massive galaxies. For the models where the sizes and masses can be constrained to compare more fairly with observations, systematic differences with our results about the cosmic UV luminosity (or SFR) density are found at the $\times 4-10$ level (depending on the specific model), indicating that the observations reveal a much more active universe in the production of photons (through star formation or nuclear activity) in the first 500 Myr on $\sim 2$ kpc scales, especially at $z \gtrsim 11$ (i.e., first 350 Myr of the universe).

## Acknowledgments

The authors dedicate this Letter to the memory of the members of the European Consortium MIRI Team Hans Ulrik Nørgaard-Nielsen and Olivier Le Fèvre, R.I.P. We thank the anonymous referee for their constructive comments that certainly helped to the robustness of this work. This work has made use of the Rainbow Cosmological Surveys Database, which is operated by the Centro de Astrobiología (CAB), CSIC-INTA, partnered with the University of California Observatories at Santa Cruz (UCO/Lick, UCSC). P.G.P.-G. and L.C. acknowledge support from Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through grant PGC2018-093499-B-I00. L.C. acknowledges financial support from Comunidad de Madrid under Atracción de Talento grant 2018-T2/TIC-11612. L.C. and J.Á.-M. acknowledge support by grant PIB2021-127718NB-100. G.S.W. acknowledges funding from the UK Science and
Technology Facilities Council, and the UK Space Agency. D.L. and J.H. were supported by a VILLUM FONDEN Investigator grant to J.H. (project number 16599). A.A.H. acknowledges support from PID2021-124665NB-100 funded by the Spanish Ministry of Science and Innovation and the State Agency of Research MCIN/AEI/10.13039/ 501100011033 and by ERDF A way of making Europe. T.G. and S.G. acknowledge the support of the Cosmic Dawn Center of Excellence (DAWN) funded by the Danish National Research Foundation under the grant 140. T.R.G. acknowledges support from the Carlsberg Foundation (grant No. CF20- 0534). K.I.C. acknowledges funding from the Dutch Research Council (NWO) through the award of the Vici grant VI. C.212.036. J.P.P. and T.V.T. acknowledge funding from the UK Science and Technology Facilities Council, and the UK Space Agency. O.I. acknowledges the funding of the French Agence Nationale de la Recherche for the project iMAGE (grant ANR-22-CE31-0007) and the support of the Centre National d'Etudes Spatiales (CNES). B.V. thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA). A.C.G. acknowledges support by grant PIB2021-127718NB-100. R.A.M. acknowledges support from the ERC Advanced grant 740246 (Cosmic_Gas).

Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via doi:10.17909/je9x-d314.

## ORCID iDs

Pablo G. Pérez-González https://orcid.org/0000-0003-4528-5639

Luca Costantin https://orcid.org/0000-0001-6820-0015

Danial Langeroodi https://orcid.org/0000-0001-5710-8395

Pierluigi Rinaldi https://orcid.org/0000-0002-5104-8245

Marianna Annunziatella https://orcid.org/0000-0002-8053-8040

Olivier Ilbert https://orcid.org/0000-0002-7303-4397

Luis Colina https://orcid.org/0000-0002-9090-4227

Hans Ulrik Nørgaard-Nielsen https://orcid.org/0000-0003-1693-2868

Thomas R. Greve https://orcid.org/0000-0002-2554-1837

Göran Östlin https://orcid.org/0000-0002-3005-1349

Gillian Wright https://orcid.org/0000-0001-7416-7936

Almudena Alonso-Herrero https://orcid.org/0000-0001-6794-2519

Javier Álvarez-Márquez https://orcid.org/0000-0002-7093-1877

Karina I. Caputi https://orcid.org/0000-0001-8183-1460

Andreas Eckart https://orcid.org/0000-0001-6049-3132

Olivier Le Fèvre https://orcid.org/0000-0001-5891-2596

Álvaro Labiano https://orcid.org/0000-0002-0690-8824

Macarena García-Marín https://orcid.org/0000-0003-4801-0489

Jens Hjorth https://orcid.org/0000-0002-4571-2306

Sarah Kendrew https://orcid.org/0000-0002-7612-0469

John P. Pye https://orcid.org/0000-0002-0932-4330

Paul van der Werf https://orcid.org/0000-0001-5434-5942

Fabian Walter https://orcid.org/0000-0003-4793-7880

Martin Ward https://orcid.org/0000-0003-1810-0889

Arjan Bik https://orcid.org/0000-0001-8068-0891

Leindert Boogaard https://orcid.org/0000-0002-3952-8588

Sarah E. I. Bosman https://orcid.org/0000-0001-8582-7012

Alejandro Crespo Gómez https://orcid.org/0000-0003-2119-277X

Steven Gillman https://orcid.org/0000-0001-9885-4589

Edoardo Iani https://orcid.org/0000-0001-8386-3546

Iris Jermann https://orcid.org/0000-0002-2624-1641

Jens Melinder https://orcid.org/0000-0003-0470-8754

Romain A. Meyer https://orcid.org/0000-0001-5492-4522

Thibaud Moutard https://orcid.org/0000-0002-3305-9901

Ewine van Dishoek https://orcid.org/0000-0001-7591-1907

Thomas Henning https://orcid.org/0000-0002-1493-300X

Manuel Guedel https://orcid.org/0000-0001-9818-0588

Florian Peissker https://orcid.org/0000-0002-9850-2708

Tom Ray https://orcid.org/0000-0002-2110-1068

Bart Vandenbussche https://orcid.org/0000-0002-1368-3109

Ángela García-Argumánez https://orcid.org/0000-0002-8365-5525

Rosa María Mérida https://orcid.org/0000-0001-8115-5845

## References

Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755

Arnouts, S., Le Floc'h, E., Chevallard, J., et al. 2013, A&A, 558, A67

Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023a, arXiv:2303.15431

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023b, arXiv:2304.05378

Austin, D., Adams, N. J., Conselice, C. J., et al. 2023, arXiv:2302.04270

Bagley, M. B., Finkelstein, S. L., Rojas-Ruiz, S., et al. 2022, arXiv:2205.12980

Bagley, M. B., Pirzkal, N., Finkelstein, S. L., et al. 2023, arXiv:2302.05466

Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339

Behroozi, P., Conroy, C., Wechsler, R. H., et al. 2020, MNRAS, 499, 5702

Bhatawdekar, R., Conselice, C. J., Margalef-Bentabol, B., & Duncan, K. 2019, MNRAS, 486, 3805

Bhowmick, A. K., Somerville, R. S., Di Matteo, T., et al. 2020, MNRAS, 496, 754

Bouwens, R., Illingworth, G., Oesch, P., et al. 2023a, MNRAS, 523, 1009

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023b, MNRAS, 523, 1036

Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Cappellari, M. 2016, ARA&A, 54, 597

Carnall, A. C., Begley, R., McLeod, D. J., et al. 2023, MNRAS, 518, L45

Castellano, M., Fontana, A., Treu, T., et al. 2023, ApJL, 948, L14

Chabrier, G. 2003, PASP, 115, 763

Conroy, C., & Gunn, J. E. 2010, FSPS: Flexible Stellar Population Synthesis, Astrophysics Source Code Library, ascl:1010.043

Costantin, L., Pérez-González, P. G., Vega-Ferrero, J., et al. 2023, ApJ, 946, 71

Cowan, J. J., Sneden, C., Burles, S., et al. 2002, ApJ, 572, 861

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Creevey, O. L., Thévenin, F., Berio, P., et al. 2015, A&A, 575, A26

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, NatAs, 7, 622

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Dayal, P., Choudhury, T. R., Pacucci, F., & Bromm, V. 2017, MNRAS, 472, 4414

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Efstathiou, G., Ellis, R. S., & Peterson, B. A. 1988, MNRAS, 232, 431

Finkelstein, S. L., Bagley, M. B., Arrabal Haro, P., et al. 2022, ApJL, 940, L55

Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJL, 946, L13

Finkelstein, S. L., Ryan, R. E. J., Papovich, C., et al. 2015, ApJ, 810, 71

Friel, E. D. 1995, ARA&A, 33, 381

Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023, ApJL, 949, L25

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016b, A&A, 595, A2

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016a, A&A, 595, A1

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Harikane, Y., Nakajima, K., Ouchi, M., et al. 2023, arXiv:2304.06658
Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5

Haslbauer, M., Kroupa, P., Zonoozi, A. H., & Haghi, H. 2022, ApJL, 939, L31

Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2019, S&C, 29, 891

Ilbert, O., Arnouts, S., Floc'h, E. L., et al. 2015, A&A, 579, A2

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Inayoshi, K., Harikane, Y., Inoue, A. K., Li, W., & Ho, L. C. 2022, ApJL, 938, L10

Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73

Isobe, Y., Ouchi, M., Nakajima, K., et al. 2023, arXiv:2301.06811

Jiang, L., Kashikawa, N., Wang, S., et al. 2021, NatAs, 5, 256

Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22

Kannan, R., Garaldi, E., Smith, A., et al. 2022, MNRAS, 511, 4005

Kauffmann, O. B., Ilbert, O., Weaver, J. R., et al. 2022, A&A, 667, A65

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Koposov, S., Speagle, J., Aarbary, K., et al. 2022, joshspeagle/dynesty: v1.2.3, v1.2.3, Zenodo, doi:10.5281/zenodo.6609296

Kriek, M., & Conroy, C. 2013, ApJL, 775, L16

Kron, R. G. 1980, ApJS, 43, 305

Labbé, I., Franx, M., Rudnick, G., et al. 2003, AJ, 125, 1107

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Langeroodi, D., Hjorth, J., Chen, W., et al. 2022, arXiv:2212.02491

Larson, R. L., Hutchison, T. A., Bagley, M., et al. 2022, arXiv:2211.10035

Leitherer, C., Li, I. H., Calzetti, D., & Heckman, T. M. 2002, ApJS, 140, 303

Lovell, C. C., Vijayan, A. P., Thomas, P. A., et al. 2021, MNRAS, 500, 2127

Lovell, M. R., Frenk, C. S., Eke, V. R., et al. 2014, MNRAS, 439, 300

Maio, U., & Viel, M. 2023, A&A, 672, A71

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Marks, M., Kroupa, P., Dabringhausen, J., & Pawlowski, M. S. 2012, MNRAS, 422, 2246

McDermid, R. M., Alatalo, K., Blitz, L., et al. 2015, MNRAS, 448, 3484

McKinney, J., Finnerty, L., Casey, C., et al. 2023, ApJL, 946, L39

McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812

Mérida, R. M., Pérez-González, P. G., Sánchez-Blázquez, P., et al. 2023, ApJ, in press, arXiv:2303.16234

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJS, 206, 8

Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022b, arXiv:2208.02794

Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022a, ApJL, 940, L14

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Nelson, D., Pillepich, A., Genel, S., et al. 2015, A&C, 13, 12

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105

Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2023, ApJL, 946, L16

Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Quadri, R., Marchesini, D., van Dokkum, P., et al. 2007, AJ, 134, 1103

Rinaldi, P., Caputi, K. I., Costantin, L., et al. 2023, arXiv:2301.10717

Rinaldi, P., Caputi, K. I., van Mierlo, S. E., et al. 2022, ApJ, 930, 128

Roberts-Borsani, G., Treu, T., Chen, W., et al. 2023, Natur, 618, 480

Robertson, B. E. 2022, ARA&A, 60, 121

Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, NatAs, 7, 611

Rodighiero, G., Bisigello, L., Iani, E., et al. 2023, MNRAS, 518, L19

Rojas-Ruiz, S., Finkelstein, S. L., Bagley, M. B., et al. 2020, ApJ, 891, 146

Saito, S., de la Torre, S., Ilbert, O., et al. 2020, MNRAS, 494, 199

Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023, arXiv:2301.06696

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schechter, P. 1976, ApJ, 203, 297

Schlaufman, K. C., Thompson, I. B., & Casey, A. R. 2018, ApJ, 867, 98

Smith, T. L., Lucca, M., Poulin, V., et al. 2022, PhRvD, 106, 043526

Speagle, J. S. 2020, MNRAS, 493, 3132

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Steinhardt, C. L., Sneppen, A., Hensley, H., et al. 2022a, ApJ, 934, 22

Steinhardt, C. L., Sneppen, A., Mostafa, B., et al. 2022b, ApJ, 931, 58

Tang, M., Stark, D. P., Chen, Z., et al. 2023, arXiv:2301.07072

Trapp, A. C., & Furlanetto, S. R. 2020, MNRAS, 499, 2401

Trenti, M., & Stiavelли, M. 2008, ApJ, 676, 767

Tumlinson, J. 2010, ApJ, 708, 1398

Ucci, G., Dayal, P., Hutter, A., et al. 2021, MNRAS, 506, 202

Vagnozzi, S., Pacucci, F., & Loeb, A. 2022, JHEAp, 36, 27

Vijayan, A. P., Lovell, C. C., Wilkins, S. M., et al. 2021, MNRAS, 501, 3289

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Whitaker, K. E., Ashas, M., Illingworth, G., et al. 2019, ApJS, 244, 16

Whitler, L., Endsley, R., Stark, D. P., et al. 2023, MNRAS, 519, 157

Williams, C. C., Curtis-Lake, E., Hainline, K. N., et al. 2018, ApJS, 236, 33

Williams, H., Kelly, P. L., Chen, W., et al. 2023, Sci, 380, 416

Yan, H., Ma, Z., Ling, C., Cheng, C., & Huang, J.-S. 2023, ApJL, 942, L9

Yue, B., Castellano, M., Ferrara, A., et al. 2018, ApJ, 868, 115

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., et al. 2020, MNRAS, 496, 4574

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019, MNRAS, 483, 2983

Zackrisson, E., Rydberg, C.-E., Schaerer, D., Östlin, G., & Tuli, M. 2011, ApJ, 740, 13

Zavala, J. A., Buat, V., Casey, C. M., et al. 2023, ApJL, 943, L9