# CEERS Key Paper. V. Galaxies at 4 < z < 9 Are Bluer than They Appear-Characterizing Galaxy Stellar Populations from Rest-frame 1 µm Imaging

Casey Papovich, Justin W. Cole, Guang Yang, Steven L. Finkelstein, Guillermo Barro, Véronique Buat, Denis Burgarella, Pablo G. Pérez-González, Paola Santini, Lise-Marie Seillé, et al.

► To cite this version:

Casey Papovich, Justin W. Cole, Guang Yang, Steven L. Finkelstein, Guillermo Barro, et al. CEERS Key Paper. V. Galaxies at 4 < z < 9 Are Bluer than They Appear-Characterizing Galaxy Stellar Populations from Rest-frame 1 µm Imaging. The Astrophysical journal letters, 2023, 949, 10.3847/2041-8213/acc948. insu-04479082

HAL Id: insu-04479082

https://insu.hal.science/insu-04479082

Submitted on 5 Mar 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CEERS Key Paper. V. Galaxies at $4 < z < 9$ Are Bluer than They Appear–Characterizing Galaxy Stellar Populations from Rest-frame ~1 µm Imaging

Casey Papovich<sup>1,2</sup>, Justin W. Cole<sup>1,2</sup>, Guang Yang<sup>3,4</sup>, Steven L. Finkelstein<sup>5</sup>, Guillermo Barro<sup>6</sup>, Véronique Buat<sup>7</sup>, Denis Burgarella<sup>7</sup>, Pablo G. Pérez-González<sup>8</sup>, Paola Santini<sup>9</sup>, Lise-Marie Seillé<sup>7</sup>, Lu Shen<sup>1,2</sup>, Pablo Arrabal Haro<sup>10</sup>, Micaela B. Bagley<sup>5</sup>, Eric F. Bell<sup>11</sup>, Laura Bisigello<sup>12,13</sup>, Antonello Calabrò<sup>9</sup>, Caitlin M. Casey<sup>5</sup>, Marco Castellano<sup>9</sup>, Katherine Chworowsky<sup>5,30</sup>, Nikko J. Cleri<sup>1,2</sup>, Luca Costantin<sup>8</sup>, M. C. Cooper<sup>14</sup>, Mark Dickinson<sup>10</sup>, Henry C. Ferguson<sup>15</sup>, Adriano Fontana<sup>9</sup>, Mauro Giavalisco<sup>16</sup>, Andrea Grazian<sup>13</sup>, Norman A. Grogin<sup>15</sup>, Nimish P. Hathi<sup>15</sup>, Benne W. Holwerda<sup>17</sup>, Taylor A. Hutchison<sup>18,31</sup>, Jeyhan S. Kartaltepe<sup>19</sup>, Lisa J. Kewley<sup>20</sup>, Allison Kirkpatrick<sup>21</sup>, Dale D. Kocevski<sup>22</sup>, Anton M. Koekemoer<sup>15</sup>, Rebecca L. Larson<sup>5,30</sup>, Arianna S. Long<sup>5,32</sup>, Ray A. Lucas<sup>15</sup>, Laura Pentericci<sup>9</sup>, Nor Pirzkal<sup>15,23</sup>, Swara Ravindranath<sup>15</sup>, Rachel S. Somerville<sup>24</sup>, Jonathan R. Trump<sup>25</sup>, Stephanie M. Urbano Stawinski<sup>14</sup>, Benjamin J. Weiner<sup>26</sup>, Stephen M. Wilkins<sup>27,28</sup>, L. Y. Aaron Yung<sup>18,31</sup>, and Jorge A. Zavala<sup>29</sup>

Received 2022 December 29; revised 2023 March 23; accepted 2023 March 24; published 2023 May 25

Abstract

We present results from the Cosmic Evolution Early Release Survey on the stellar population parameters for 28 galaxies with redshifts $4 < z < 9$ using imaging data from the James Webb Space Telescope (JWST) Mid-Infrared Instrument (MIRI) combined with data from the Hubble Space Telescope and the Spitzer Space Telescope. The JWST/MIRI 5.6 and 7.7 µm data extend the coverage of the rest-frame spectral energy distribution to nearly 1 µm for galaxies in this redshift range. By modeling the galaxies’ SEDs the MIRI data show that the galaxies have, on average, rest-frame UV (1600 Å)—I-band colors 0.4 mag bluer than derived when using photometry that lacks MIRI. Therefore, the galaxies have lower ratios of stellar mass to light. The MIRI data reduce the stellar masses by $\langle\Delta \log M_*\rangle = 0.25$ dex at $4 < z < 6$ and 0.37 dex at $6 < z < 9$. This also reduces the star formation rates (SFRs) by $\langle\Delta \log \text{SFR}\rangle = 0.14$ dex at $4 < z < 6$ and 0.27 dex at $6 < z < 9$. The MIRI data also improve constraints on the allowable stellar mass formed in early star formation. We model this using a star formation history that includes

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

<footnote n="1"><sup>1</sup> Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA; papovich@tamu.edu</footnote>

<footnote n="2"><sup>2</sup> George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242 USA</footnote>

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

<footnote n="4"><sup>4</sup> SRON Netherlands Institute for Space Research, Postbus 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="5"><sup>5</sup> Department of Astronomy, The University of Texas at Austin, Austin, TX, USA</footnote>

<footnote n="6"><sup>6</sup> Department of Physics, University of the Pacific, Stockton, CA 90340 USA</footnote>

<footnote n="7"><sup>7</sup> Aix Marseille Univ, CNRS, CNES, LAM Marseille, France</footnote>

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

<footnote n="9"><sup>9</sup> INAF—Osservatorio Astronomico di Roma, via di Frascati 33, I-00078 Monte Porzio Catone, Italy</footnote>

<footnote n="10"><sup>10</sup> NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA</footnote>

<footnote n="11"><sup>11</sup> Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI 48109-1107, USA</footnote>

<footnote n="12"><sup>12</sup> Dipartimento di Fisica e Astronomia “G.Galilei,”Università di Padova, Via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote n="13"><sup>13</sup> INAF—Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122, Padova, Italy</footnote>

<footnote n="14"><sup>14</sup> Department of Physics & Astronomy, University of California, Irvine, 4129 Reines Hall, Irvine, CA 92697, USA</footnote>

<footnote n="15"><sup>15</sup> Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA</footnote>

<footnote n="16"><sup>16</sup> University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA</footnote>

<footnote n="17"><sup>17</sup> Physics & Astronomy Department, University of Louisville, 40292 KY, Louisville, USA</footnote>

<footnote n="18"><sup>18</sup> Astrophysics Science Division, NASA Goddard Space Flight Center, 8800 Greenbelt Rd, Greenbelt, MD 20771, USA</footnote>

<footnote n="19"><sup>19</sup> Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA</footnote>

<footnote n="20"><sup>20</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="21"><sup>21</sup> Department of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA</footnote>

<footnote n="22"><sup>22</sup> Department of Physics and Astronomy, Colby College, Waterville, ME 04901, USA</footnote>

<footnote n="23"><sup>23</sup> ESA/AURA Space Telescope Science Institute, USA</footnote>

<footnote n="24"><sup>24</sup> Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>

<footnote n="25"><sup>25</sup> Department of Physics, 196 Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA</footnote>

<footnote n="26"><sup>26</sup> MMT/Steward Observatory, University of Arizona, 933 N. Cherry Ave., Tucson, AZ 85721, USA</footnote>

<footnote n="27"><sup>27</sup> Astronomy Centre, University of Sussex, Falmer, Brighton BN1 9QH, UK</footnote>

<footnote n="28"><sup>28</sup> Institute of Space Sciences and Astronomy, University of Malta, Msida MSD 2080, Malta</footnote>

<footnote n="29"><sup>29</sup> National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="30"><sup>30</sup> NSF Graduate Fellow.</footnote>

<footnote n="31"><sup>31</sup> NASA Postdoctoral Fellow.</footnote>

<footnote n="32"><sup>32</sup> NASA Hubble Fellow.</footnote>
both a “burst” at $z_f=100$ and a slowly varying (“delayed-$\tau$”) model. The MIRI data reduce the allowable stellar mass by 0.6 dex at $4 < z < 6$ and by $\approx 1$ dex at $6 < z < 9$. Applying these results globally, this reduces the cosmic stellar-mass density by an order of magnitude in the early Universe ($z \approx 9$). Therefore, observations of rest-frame $\gtrsim 1$ µm are paramount for constraining the stellar-mass buildup in galaxies at very high redshifts.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Galaxy formation (595); Cosmology (343); Galaxy colors (586); Galaxy properties (615); Galaxy masses (607); Spectral energy distribution (2129)*

*Supporting material: figure sets*

# 1. Introduction

There is growing evidence that galaxies must have started forming stars very quickly following the Big Bang. Theory predicts the first stars should form at $z \gtrsim 20$ (e.g., Barkana & Loeb 2001; Miralda-Escudé 2003; Yoshida et al. 2003; Wise et al. 2012; Visbal et al. 2020). The ionization from these sources is needed to explain observations that the hydrogen-neutral fraction of the intergalactic medium (IGM) was 50% by $z \sim 8$ (Planck Collaboration et al. 2020).<fnref n="33" /> Indeed, early observations from JWST have already identified candidates for galaxies at $z \gtrsim 15$ (Finkelstein et al. 2022b; Curtis-Lake et al. 2023; Donnan et al. 2023; Robertson et al. 2023). Spectroscopy from JWST of galaxies at $z \sim 8–9$ shows emission lines from heavy elements that appear to require metallicities of $\approx 5\%–10\% Z_\odot$ (Arellano-Córdova et al. 2022; Fujimoto et al. 2022; Heintz et al. 2023; Langeroodi et al. 2022; Matthee et al. 2022; Schaerer et al. 2022; Trump et al. 2023; Curti et al. 2023; Katz et al. 2023), implying these galaxies have experienced at least one (and probably multiple) generation(s) of previous stars. This is consistent with earlier detections of metal lines in $z > 6$ galaxies from ground-based telescopes (e.g., Stark et al. 2017; Hutchison et al. 2019). All of these results point to the fact that star formation began early and those early generations of stars enriched the Universe with heavy elements.

It is then important to consider how we may constrain the history of star formation in these early galaxies. The number of stars, and therefore the stellar mass, in galaxies appears to rise rapidly. The (comoving) stellar mass density in galaxies at $z \sim 5–6$ is already 1% of the present value (e.g., Madau & Dickinson 2014; Finkelstein 2016). Obviously, these stars must form at earlier times, and because the age of the Universe is only $\approx 1$ Gyr at these redshifts, the time to form these stars is relatively short.

Even early JWST observations find some evidence for massive galaxies at $z > 7$ (Labbé et al. 2023), with some candidate objects having masses, $\log M_*/M_\odot > 11$, as large as the stellar mass of the Milky Way today (e.g., Papovich et al. 2015). Such objects would be in tension with galaxy formation models (Boylan-Kolchin 2023) where there are not sufficient numbers of massive dark-matter halos to support these objects, even if all the baryons in the halos are in the form of stellar mass. However, the uncertainty in these measurements is in the assumed star formation histories, the contributions of emission lines to the photometric measurements from broadband data (e.g., Endsley et al. 2022; Pérez-González et al. 2023; Steinhardt et al. 2022), and the effects of young stellar populations “outshining” older stellar populations in the integrated emission of galaxies (e.g., Giménez-Arteaga et al.

2023). Clearly, there are unknown systematics in the assumptions of the data analysis, or missing physics in our theoretical understanding of stellar populations and galaxy formation, or some combination of all of these things. It is therefore crucial to understand constraints on the stellar masses (which are the integral of the star formation histories) as much as possible.

Motivated by these issues, in this paper we use new data from JWST to better constrain the stellar masses, star formation rates (SFRs), and star formation histories of galaxies during the first one and a half billion years after the Big Bang ($z > 4$). One of the problems with initial studies from JWST is that they currently rely entirely on observations from JWST’s Near-IR Camera (NIRcam), which only probes to wavelengths $\lesssim 5$ µm, or about 6000 Å rest frame for $z = 6–7$ galaxies. This complicates the ability to disentangle massive galaxies with older stellar populations from younger, dusty galaxies or galaxies with emission lines with extreme equivalent widths (see discussion in Antwi-Danso et al. (2023), and recent work by Giménez-Arteaga et al. (2023), who find evidence for older stellar populations mixed with recent bursts in spatially resolved studies using JWST/NIRCam data). To better constrain the SEDs of these galaxies requires observations at longer wavelengths. This is where JWST/MIRI is important, as it has the sensitivity to detect $z \approx 10$ galaxies at rest-frame 1 µm (see Bisigello et al. 2017). Previous work on this subject has been limited to data from the Spitzer Space Telescope, which is primarily sensitive to the emission of distant galaxies at 3.6–8.0 µm. JWST offers immense gains to Spitzer: JWST has a collecting area that is 45 times larger than that of Spitzer, and the larger aperture provides image quality (i.e., angular resolution) that is improved by a factor of order 10 (Rigby et al. 2023). These gains are especially manifest at longer wavelengths, and they make JWST 5.6 and 7.7 µm data vastly more sensitive than Spitzer.

The outline for the paper is as follows. In Section 2, we discuss the CEERS data set and the ancillary data sets used in this study. We also discuss the processes to create (and validate) the flux densities of galaxies in the CEERS JWST/MIRI data at 5.6 and 7.7 µm. In Section 2.2, we discuss the sample of $4 < z < 9$ galaxies used in this study, and we present the MIRI data for these objects. In Section 3, we discuss the analysis methods to derive constraints on the galaxy stellar populations. In Section 4, we discuss the resulting improvements that including the MIRI data provide on constraints on the galaxies’ stellar populations (specifically their stellar masses and SFRs) derived from fitting stellar population models to the observed photometry. In Section 4.4, we discuss constraints on the range of allowed stellar masses in these high-redshift galaxies by allowing for an early (“maximally old”) burst of stars at $z = 100$, and we show that adding the MIRI data improves the limit on this hypothetical population of $z = 100$ stars by a factor of 6–10 for galaxies at $4 < z < 9$. In Section 5, we discuss the implications these constraints have

<footnote n="33">³³ The Planck Collaboration et al. (2020) analysis suggests a nonzero optical depth of CMB photons scattering off free electrons at $z \approx 15$, which implies ionization of the IGM had begun by this epoch.</footnote>
for our understanding of galaxy colors, stellar populations at these high redshifts, and the evolution of the galaxy stellar mass density, in particular during the epoch of reionization, and what this could mean for future studies of galaxies at higher redshifts (from JWST). In Section 6, we present our conclusions and prospects for future studies.

Throughout this work, we use a flat cosmology with $\Omega_{m,0} = 0.315$, $H_0 = 67.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Planck Collaboration et al. 2020). All magnitudes reported here are on the Absolute Bolometric (AB) system (Oke & Gunn 1983). Throughout this work, we use the Chabrier (2003) initial mass function (IMF) for all stellar masses and SFRs. We denote magnitudes measured in the MIRI F560W and F770W bands as [5.6] and [7.7], respectively. Similarly, we denote magnitudes measured in IRAC Channel 1 (3.6 µm), Channel 2 (4.5 µm), Channel 3 (5.8 µm), and Channel 4 (8.0 µm) as [3.6], [4.5], [5.8], and [8.0], respectively.

## 2. Data and Sample

### 2.1. MIRI Catalog

We use the data release DR0.5 images produced by the CEERS team for the MIRI 3 and MIRI 6 fields (see Finkelstein et al. 2022b).<fnref n="34" /> These data were acquired in 2022 June 21 and 22. The properties of the data and the data reduction are discussed elsewhere (G. Yang et al. 2023, in preparation), but we provide a summary here. The data were processed using the JWST Calibration Pipeline (v1.7.2) using the default parameters for stage 1 and 2. We then removed the backgrounds with a custom routine that combines images taken in the same bandpass but from different fields and/or dither positions, rejecting pixels in each image that contain galaxies, in order to create a “super-background” image. We then removed this background from each image and applied an astrometric correction to each image prior to processing them with stage 3 of the pipeline. This produced the final science images (extension i2d), rms images (extension rms, which account for Poisson, readout, and correlated pixel noise; see G. Yang et al. (2023, in preparation)), and weight maps (wht) for each field with a pixel scale of 0.′′09, registered astrometrically to the existing HST/CANDELS v1.9 WFC3 and ACS images (see Koekemoer et al. 2011; Bagley et al. 2023). Our tests find that the MIRI images achieve limiting 5σ depths in the MIRI images of 26.5 and 27.1 AB mag measured in 0.′′45 diameter apertures.

For the purpose of this study, we are interested in sources detected in the MIRI data, so we create a catalog of sources derived from these images. Prior to object detection, we convolved the 5.6 µm image to match the image quality of the 7.7 µm image. For this step, we constructed an “effective” point-source function (ePSF) for each image by identifying unblended stars using the photutils (v1.5.0) detection task and modeling them with the photutils psf task. This produced model ePSFs with measured full width at half maxima (FWHM) of 0.′′24 and 0.′′28, for the F560W and F770W images, respectively. This is consistent with the expected image quality, but takes into account the exact dithering and reduction steps for the CEERS MIRI data. We then used PyPHER (Boucaud et al. 2016b) to construct a convolution kernel to match the image quality of the model

<table><caption>**Table 1**
CEERS MIRI F560W and F770W SExtractor Parameter Settings</caption><thead><tr><th>SEExtractor Parameter<br>(1)</th><th>Value<br>(2)</th></tr></thead><tbody><tr><td>DETECT_MINAREA</td><td>10 pixels</td></tr><tr><td>DETECT_THRESH</td><td>1.3</td></tr><tr><td>ANALYSIS_THRESH</td><td>1.3</td></tr><tr><td>FILTER_NAME</td><td>Gauss_2.5_5x5<sup>a</sup></td></tr><tr><td>WEIGHT_TYPE</td><td>MAP_WEIGHT,MAP_RMS</td></tr><tr><td>DEBLEND_NTHRESH</td><td>32</td></tr><tr><td>DEBLEND_MINCONT</td><td>0.005</td></tr><tr><td>MAG_ZEROPOINT</td><td>25.701<sup>b</sup></td></tr><tr><td>PIXEL_SCALE</td><td>0.′′09</td></tr><tr><td>BACK_TYPE</td><td>AUTO</td></tr><tr><td>BACK_FILTERSIZE</td><td>5 pixels</td></tr><tr><td>BACK_SIZE</td><td>32 pixels</td></tr><tr><td>BACKPHOTO_THICK</td><td>8</td></tr><tr><td>BACKPHOTO_TYPE</td><td>LOCAL</td></tr><tr><td>SEEING_FWHM</td><td>0.′′3</td></tr></tbody></table>

**Notes.** SE was run using the weighted sum of the PSF-matched F560W and F770W images for detection, and using the images separately for photometry. All other SE parameters are set to the program defaults (for SExtractor v.2.19.5).

a This is a Gaussian kernel with $\sigma = 2.5$ pixels and size $5 \times 5$ pixel² used to filter the image for source detection.

b The AB magnitude zero point for the images, converting from the JWST default of MJy sr⁻¹ to µJy pixel⁻¹ at the 0.′′09 pixel⁻¹ scale.

ePSFs. We applied these kernels to each F560W image, creating a “PSF-matched” image. Our tests on point sources in the PSF-matched F560W and F770W images show that we measure the same fraction of light to better than 2% in fixed circular apertures of radii larger than 0.′′35.

We then created a detection image constructed from the sum of the MIRI F560W and F770W science images (using the extension sci) weighted by the appropriate weight image (using the extension wht). We also created a detection weight map as the sum of the weights for these images. We then created F560W and F770W catalogs using Source Extractor (SE, version 2.19.5; Bertin & Arnouts 1996) in “dual-image” mode using the detection image and its weight map for object detection, where we measured photometry in the PSF-matched F560W and F770W image. We used the parameters in Table 1. We then measured fluxes and magnitudes using circular apertures of 0.′′9 diameter, and we scaled these to a total aperture (MAG_AUTO) measured for each source in the detection image. Uncertainties for each object are measured from the rms image in the same apertures, and scaled to a total magnitude in the same way.

We compared the MIRI flux densities for sources with F560W and F770W detections against those for bright objects from existing IRAC 5.8 and 8.0 µm catalogs (Stefanon et al. 2017). For bright objects ([5.8] or [8.0] ≤ 22 AB mag) in the IRAC data, we measure small offsets of $\Delta m = [5.8] - [5.6] = 0.16$ mag between the IRAC 5.8 and MIRI 5.6 data, and $\Delta m = [8.0] - [7.7] = 0.07$ mag between the IRAC 8.0 and MIRI 7.7 data (i.e., the MIRI flux densities are slightly brighter). Most of these offsets can be explained by differences in the shape of the MIRI and IRAC passbands and because of differences in the angular resolution of the instruments (MIRI has a PSF FWHM smaller by a factor of more than seven). These tests are

<footnote n="34">³⁴ https://ceers.github.io/releases.html</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustration of galaxy spectra (in relative units of erg s⁻¹ cm⁻² Å⁻¹) compared to the broadband data for the observations in this work. The bottom panel shows the relative transmission functions for the HST/ACS and WFC3 filters (ACS F606W and F814, and WFC3 F125W and F160W), Spitzer/IRAC 3.6 µm (Ch1) and 4.5 µm (Ch2), and JWST/MIRI F560W and F770W. The top panel shows model spectra of star-forming galaxies at z = 5.1 and z = 8.7, which match two galaxies with spectroscopic redshifts in this sample. Key emission lines and features are labeled. The MIRI data probe the shape of the galaxy spectral energy distributions to the 8000 Å rest frame, even for galaxies with z = 9.</figcaption></figure>

discussed more fully in G. Yang et al. (2023, in preparation), but this gives us confidence that the MIRI data are calibrated to better than ≈0.1–0.2 mag.

## 2.2. Galaxy Sample

For this study, we use galaxies identified in the CEERS/ MIRI first-epoch fields with photometric redshifts 4.3 < z < 10 from Finkelstein et al. (2022a). The lower redshift bound is selected to ensure the HST photometric data probe the redshifted Lyman break. The upper redshift limit includes the highest-redshift galaxies detectable by HST/WFC3 data (Bouwens et al. 2019; Finkelstein et al. 2022a). This is illustrated in Figure 1, which shows that, for galaxies around z ~ 5 and z ~ 9, the HST/ACS and WFC3 data constrain this break. This improves the quality of the sample, as compared to, for example, using galaxies at z ≈ 3 where the Lyman break shifts blueward of the HST/ACS F606W band.

As mentioned above, the parent sample for our study is the catalog from Finkelstein et al. (2022a), which uses the existing HST/ACS, WFC3, and Spitzer/IRAC 3.6 and 4.5 µm data to select photometric samples of galaxies at these high redshifts. The data include both the imaging from the original CANDELS survey (HST/ACS F606W and F814W, and WFC3 F125W and F160W) and additional imaging from WFC3 F140W (see Footnote 2). Finkelstein et al. (2022a) then use these data to measure photometric redshifts and redshift probability distribution functions, P(z) for each object.

We matched objects from the catalog from Finkelstein et al. (2022a) with objects in our MIRI catalogs that are detected with S/N > 3 in either the MIRI 5.6 or 7.7 µm data (Section 2.1) using a matching radius of 0.5. This sample includes 29 objects, though one object was later identified as a foreground star and removed. Table 2 lists the observed

properties of the 28 galaxies in the sample, their ID numbers from Finkelstein et al. (2022a), their HST/F160W flux densities, and their MIRI 5.6 and 7.7 µm flux densities. The table includes the photometric redshifts derived by Finkelstein et al. (2022a), including the 16th and 84th percentile range from the P(z) used for object selection. The table also includes the amount of the integrated P(z) contained between Δz = ± 0.5 of the stated redshift, for example,

$$P(z = z_c) = \int_{z_c - 0.5}^{z_c + 0.5} P(z') \, dz', \quad (1)$$

in bins with central redshifts of z_c = 4, 5, 6, 7, 8, and 9 (see Finkelstein et al. 2022a). These integrated probabilities indicate a likelihood that a given galaxy is within the redshift bin to which it is assigned. For our analysis, we rederive the photometric redshifts below (from SED modeling that includes the new 5.6 and 7.7 µm MIRI data) but include this here because we use the P(z) in Table 2 as a prior likelihood on the SED fitting (discussed in Section 3 below).

Three of the galaxies in our sample have spectroscopic redshifts. This includes a previously known galaxy (ID 6811 in our catalog) with z = 8.683 from Zitrin et al. (2015), as well as two new redshifts obtained by the CEERS and WERLS collaborations from observations with Keck/DEIMOS and Keck/LRIS. The latter two sources are ID 37653 with z = 4.899 measured by (S. Stawinski et al. 2023a, in preparation) and ID 19180 with z = 5.077 (S. Stawinski et al. 2023b, in preparation). In all of these cases, the spectroscopic redshifts are consistent with the photometric redshifts in Table 2, and we fix the redshift to the value of the spectroscopic redshift in our analysis of the spectral energy distributions below.
<table><caption>Table 2
Observed Properties of the Galaxy Sample</caption><thead><tr><th>ID</th><th>R.A. (J2000)<br>(deg)</th><th>Decl. (J2000)<br>(deg)</th><th>F160<br>nJy</th><th>E160<br>(nJy)</th><th>F560<br>(nJy)</th><th>E560<br>(nJy)</th><th>F770<br>(nJy)</th><th>E770<br>(nJy)</th><th>z<sub>phot</sub></th><th>z<sub>16</sub></th><th>z<sub>84</sub></th><th>z<sub>spec</sub></th><th>P(z = 4)</th><th>P(z = 5)</th><th>P(z = 6)</th><th>P(z = 7)</th><th>P(z = 8)</th><th>P(z = 9)</th></tr></thead><tbody><tr><td>5090</td><td>215.04973</td><td>52.89656</td><td>102.0</td><td>12.5</td><td>158.8</td><td>17.9</td><td>148.0</td><td>14.4</td><td>4.38</td><td>4.18</td><td>4.56</td><td>-1.000</td><td>0.68</td><td>0.27</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>11329</td><td>215.04086</td><td>52.90623</td><td>49.4</td><td>6.9</td><td>63.5</td><td>20.4</td><td>44.5</td><td>15.7</td><td>4.47</td><td>4.15</td><td>4.65</td><td>-1.000</td><td>0.55</td><td>0.39</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>34813</td><td>214.97597</td><td>52.920297</td><td>82.2</td><td>16.2</td><td>67.0</td><td>16.7</td><td>43.9</td><td>19.4</td><td>4.52</td><td>4.17</td><td>4.59</td><td>-1.000</td><td>0.52</td><td>0.36</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>7600</td><td>215.04415</td><td>52.898731</td><td>564.7</td><td>29.4</td><td>3283.1</td><td>32.4</td><td>4101.7</td><td>24.8</td><td>4.57</td><td>4.13</td><td>4.65</td><td>-1.000</td><td>0.62</td><td>0.38</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>13389</td><td>215.03793</td><td>52.909329</td><td>453.9</td><td>29.7</td><td>578.3</td><td>26.9</td><td>661.0</td><td>22.0</td><td>4.57</td><td>4.43</td><td>4.61</td><td>-1.000</td><td>0.38</td><td>0.61</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>37703</td><td>214.99094</td><td>52.924279</td><td>116.7</td><td>11.8</td><td>278.9</td><td>30.4</td><td>61.0</td><td>17.5</td><td>4.60</td><td>4.47</td><td>4.70</td><td>-1.000</td><td>0.22</td><td>0.78</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>15445</td><td>215.02689</td><td>52.907215</td><td>124.5</td><td>14.4</td><td>187.8</td><td>20.2</td><td>116.6</td><td>14.9</td><td>4.60</td><td>4.21</td><td>4.84</td><td>-1.000</td><td>0.24</td><td>0.62</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>41564</td><td>214.98708</td><td>52.912734</td><td>266.4</td><td>13.4</td><td>481.6</td><td>20.8</td><td>276.8</td><td>15.7</td><td>4.63</td><td>4.56</td><td>4.71</td><td>-1.000</td><td>0.03</td><td>0.97</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>45145</td><td>215.00878</td><td>52.919869</td><td>101.2</td><td>15.2</td><td>108.6</td><td>29.6</td><td>69.0</td><td>19.3</td><td>4.68</td><td>4.37</td><td>5.20</td><td>-1.000</td><td>0.15</td><td>0.74</td><td>0.03</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>14913</td><td>215.02081</td><td>52.901523</td><td>78.3</td><td>14.1</td><td>72.8</td><td>17.8</td><td>9.4</td><td>8.9</td><td>4.68</td><td>4.47</td><td>4.94</td><td>-1.000</td><td>0.19</td><td>0.81</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>42638</td><td>214.97762</td><td>52.90349</td><td>547.0</td><td>31.8</td><td>711.4</td><td>28.5</td><td>596.3</td><td>23.4</td><td>4.71</td><td>4.57</td><td>4.77</td><td>-1.000</td><td>0.02</td><td>0.90</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>41375</td><td>215.00283</td><td>52.924312</td><td>109.2</td><td>15.4</td><td>47.2</td><td>18.1</td><td>119.1</td><td>14.8</td><td>4.74</td><td>4.60</td><td>5.13</td><td>-1.000</td><td>0.07</td><td>0.92</td><td>0.01</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>35896</td><td>214.98522</td><td>52.924266</td><td>96.3</td><td>10.2</td><td>114.5</td><td>14.8</td><td>110.0</td><td>15.3</td><td>4.76</td><td>4.35</td><td>5.33</td><td>-1.000</td><td>0.06</td><td>0.73</td><td>0.07</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>13179</td><td>215.04260</td><td>52.911968</td><td>736.4</td><td>23.7</td><td>1355.2</td><td>26.8</td><td>1519.1</td><td>21.5</td><td>4.78</td><td>4.72</td><td>4.91</td><td>-1.000</td><td>0.00</td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>19180</td><td>215.03170</td><td>52.919632</td><td>773.3</td><td>35.6</td><td>749.1</td><td>26.0</td><td>787.7</td><td>19.0</td><td>4.93</td><td>4.86</td><td>4.95</td><td>5.077<sup>a</sup></td><td>0.00</td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>37653</td><td>214.99190</td><td>52.925053</td><td>103.1</td><td>12.5</td><td>154.0</td><td>22.2</td><td>108.4</td><td>18.5</td><td>4.95</td><td>4.72</td><td>5.59</td><td>4.899<sup>b</sup></td><td>0.06</td><td>0.71</td><td>0.22</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>18449</td><td>215.02314</td><td>52.912683</td><td>69.0</td><td>10.0</td><td>185.0</td><td>24.2</td><td>165.3</td><td>17.6</td><td>5.05</td><td>4.30</td><td>5.57</td><td>-1.000</td><td>0.14</td><td>0.59</td><td>0.20</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>41545</td><td>215.00312</td><td>52.924103</td><td>384.9</td><td>18.3</td><td>749.2</td><td>21.5</td><td>656.6</td><td>17.5</td><td>5.19</td><td>5.02</td><td>5.44</td><td>-1.000</td><td>0.00</td><td>0.90</td><td>0.10</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>7818</td><td>215.02758</td><td>52.887744</td><td>289.9</td><td>21.4</td><td>421.4</td><td>24.5</td><td>399.5</td><td>20.9</td><td>5.27</td><td>4.74</td><td>5.46</td><td>-1.000</td><td>0.02</td><td>0.86</td><td>0.12</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>12773</td><td>215.03057</td><td>52.902606</td><td>56.4</td><td>12.2</td><td>17.0</td><td>7.5</td><td>71.0</td><td>14.1</td><td>5.35</td><td>4.87</td><td>5.43</td><td>-1.000</td><td>0.00</td><td>0.90</td><td>0.09</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>24007</td><td>214.95185</td><td>52.928275</td><td>47.1</td><td>4.9</td><td>107.8</td><td>18.8</td><td>56.7</td><td>12.6</td><td>5.64</td><td>5.16</td><td>5.79</td><td>-1.000</td><td>0.00</td><td>0.46</td><td>0.51</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>49365</td><td>215.00994</td><td>52.910669</td><td>143.9</td><td>14.1</td><td>408.5</td><td>26.8</td><td>541.0</td><td>20.6</td><td>5.64</td><td>5.03</td><td>5.82</td><td>-1.000</td><td>0.01</td><td>0.38</td><td>0.50</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>18441</td><td>215.03209</td><td>52.918972</td><td>97.2</td><td>11.5</td><td>199.6</td><td>28.2</td><td>169.9</td><td>22.0</td><td>6.76</td><td>5.88</td><td>7.24</td><td>-1.000</td><td>0.00</td><td>0.01</td><td>0.36</td><td>0.48</td><td>0.07</td><td>0.00</td></tr><tr><td>39096</td><td>214.98901</td><td>52.919652</td><td>64.1</td><td>8.7</td><td>159.5</td><td>17.0</td><td>57.7</td><td>11.2</td><td>6.82</td><td>6.67</td><td>7.08</td><td>-1.000</td><td>0.00</td><td>0.00</td><td>0.02</td><td>0.97</td><td>0.01</td><td>0.00</td></tr><tr><td>12514</td><td>215.03716</td><td>52.906712</td><td>71.6</td><td>13.2</td><td>142.3</td><td>17.8</td><td>50.4</td><td>13.6</td><td>7.57</td><td>6.90</td><td>8.41</td><td>-1.000</td><td>0.00</td><td>0.00</td><td>0.02</td><td>0.39</td><td>0.44</td><td>0.13</td></tr><tr><td>7364</td><td>215.03561</td><td>52.892208</td><td>292.6</td><td>17.5</td><td>575.7</td><td>28.4</td><td>611.4</td><td>23.1</td><td>8.52</td><td>7.47</td><td>8.68</td><td>-1.000</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.17</td><td>0.52</td><td>0.31</td></tr><tr><td>6811</td><td>215.03538</td><td>52.890666</td><td>314.5</td><td>13.3</td><td>426.8</td><td>21.6</td><td>404.0</td><td>16.9</td><td>8.93</td><td>8.60</td><td>9.05</td><td>8.683<sup>c</sup></td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.08</td><td>0.92</td></tr><tr><td>26890</td><td>214.96754</td><td>52.932966</td><td>133.3</td><td>8.9</td><td>101.9</td><td>21.4</td><td>60.8</td><td>16.8</td><td>9.16</td><td>8.46</td><td>9.24</td><td>-1.000</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.06</td><td>0.83</td></tr></tbody></table>

Notes. References for spectroscopic redshifts.

a S. Stawinski et al. 2023, in preparation

b WERLS.

c Zitrin et al. (2015).
<table><caption>**Table 3**
Parameter Settings for BAGPIPES</caption><thead><tr><th>Model</th><th>Parameter</th><th>Prior</th><th>Limits</th></tr></thead><tbody><tr><td rowspan="3">Star Formation History (1): Delayed-τ, Ψ ∝ (t/τ)exp(-t/τ)</td><td>e-folding timescale, τ / Gyr</td><td>Uniform</td><td>(0.01, 10)</td></tr><tr><td>age, t / Gyr</td><td>Uniform</td><td>(0.01,15)</td></tr><tr><td>stellar mass, log(M*/M<sub>☉</sub>)</td><td>Uniform</td><td>(5, 12)</td></tr><tr><td rowspan="2">Star Formation History (2): Burst at z<sub>f</sub> = 100 and delayed-τ model from (1)</td><td>burst age, t<sub>burst</sub> / Gyr</td><td>Fixed</td><td>t<sub>burst</sub> = Age(z) - Age(z<sub>f</sub> = 100)</td></tr><tr><td>burst stellar mass, log(M*/M<sub>☉</sub>)</td><td>Uniform</td><td>(0, 13)</td></tr><tr><td rowspan="5">Additional parameters for all models</td><td>dust attenuation law</td><td>...</td><td>Calzetti (2001)</td></tr><tr><td>dust attenuation, A(V) / mag</td><td>Uniform</td><td>(0, 3)</td></tr><tr><td>metallicity, Z/Z<sub>☉</sub></td><td>Uniform</td><td>(0,1)</td></tr><tr><td>ionization parameter, log U</td><td>Fixed</td><td>-2</td></tr><tr><td>redshift<sup>a</sup>, z</td><td>EAZY P(z)</td><td>(3, 15)</td></tr></tbody></table>

**Notes.**

<sup>a</sup> For galaxies with photometric redshifts the redshift prior is the posterior from the photometric redshift. For galaxies with spectroscopic redshifts, the redshift is fixed at the spectroscopic redshift.

Figure 2 shows the HST/ACS, HST/WFC3, Spitzer/IRAC, and JWST/MIRI imaging for all the objects in our sample, with the objects ordered by increasing redshift (the full figure set of all 28 objects is available online). In all cases, the galaxies show prominent “Lyman breaks” at the location of the redshifted Lyman limit and/or Lyα. In some cases, the flux density appears to be much brighter in a given passband compared to the adjacent band (for example, galaxy ID = 12514 at z = 7.6 shows evidence of enhanced emission at MIRI 5.6 µm, indicative of strong redshifted Hα). Figure 1 illustrates how the bandpasses are sensitive to different features in the SED of galaxies (using z = 5.1 and z = 8.7 as examples because these are similar to two of the objects with spectroscopic redshifts in our sample). We will return to these cases below, when we explore constraints on the galaxy stellar populations by modeling their SEDs (Section 3).

The 4 < z < 9 galaxies in our sample have MIRI [5.6]–[7.7] colors largely consistent with expectations: most objects have relatively flat ([5.6] – [7.7] ≈ 0 mag) or blue colors ([5.6]–[7.7] ≲ 0 mag). This implies that, in most cases, the MIRI data sample the continuum of galaxies. In some cases, the MIRI colors suggest very blue colors, [5.6]–[7.7] ≲ − 0.5 to − 1 mag. This could indicate the presence of an emission line in the F560W band that boosts the brightness in this band. We will explore evidence for this interpretation below.

### 3. Spectral Energy Distribution Modelling

We model the spectral energy distributions (SEDs) of the galaxies in our sample using stellar population synthesis models. Our goal is to test how the inferred properties—in particular, the stellar masses and the SFRs—of the stellar populations in the high-redshift galaxies change when we include the JWST/MIRI 5.6 and 7.7 µm data. Previous work before the launch of JWST showed that MIRI data are able to recover these quantities accurately (Bisigello et al. 2017), and here we test how they improve the constraints on the stellar population parameters. This is critical for galaxies at higher redshifts, z ≳ 4, where the rest-frame optical features shift to longer wavelengths (rest-frame 4000 Å corresponds to > 2 µm), which probe light from longer-lived stars.

Perhaps more problematic are the effects of nebular emission lines, which can litter the optical portion of the SED (see

Figure 1). There are observations that z > 2 galaxies have a higher incidence of “extreme” emission lines with rest-frame EWs up to ≈ 1000 Å (e.g., van der Wel et al. 2011; Tang et al. 2019; Tran et al. 2020; Boyett et al. 2022; Matthee et al. 2022; Pérez-González et al. 2023; Sun et al. 2022), consistent with inferences made from the > 3 µm colors of z > 6 galaxies (Smit et al. 2015; Roberts-Borsani et al. 2016; Castellano et al. 2017; Hutchison et al. 2019; Endsley et al. 2021). Because the EW scales with redshift as (1 + z), this implies these lines have a stronger impact for high-redshift galaxies for bandpasses of fixed wavelength width (e.g., Papovich et al. 2001; Burgarella et al. 2023).

We model each galaxy by fitting HST/ACS and WFC3, Spitzer/IRAC, and JWST/MIRI data with stellar population models using BAGPIPES (Carnall et al. 2018). BAGPIPES is a Bayesian SED-fitting code that models the multiband photometry (flux densities) with stellar population synthesis models formed over a wide range of user-defined parameters. The code has flexibility with regard to the type of stellar population synthesis models, star formation history, dust attenuation, and nebular emission. It has the ability to incorporate prior knowledge on parameters. The code then computes a probability density for model parameters (i.e., posteriors) given the data by calculating a likelihood weighted by priors on the parameters, and it samples the posteriors for the parameters using the MultiNest nested sampling algorithm (see Feroz et al. 2009; Carnall et al. 2018).

Table 3 lists the range of parameters considered for the SED fitting in this study. For all models, we use stellar population synthesis models from Bruzual & Charlot (2003) formed with a Chabrier IMF. The table defines the parameters and their range of parameter values we explored. In most cases, we adopt uniform priors on these parameters, as listed in Table 3, with two exceptions. The first is related to the nebular ionization parameter, which controls the strength of the nebular emission features. Current evidence from spectroscopy (e.g., Le Fèvre et al. 2015; Oesch et al. 2015; Stark et al. 2015; Sanders et al. 2016; Laporte et al. 2017; Stark et al. 2017; Sanders et al. 2020; Backhaus et al. 2022; Papovich et al. 2022), including recent JWST spectroscopy (Brinchmann 2022; Schaerer et al. 2022; Trump et al. 2023), shows that emission lines are common in star-forming galaxies at z ≳ 1 (and the strength appears to
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Montage of images of a subset of the galaxies used in this study, ordered by increasing redshift. Each row shows images for one galaxy (labeled by galaxy ID and redshift). The images are (left to right), ACS F606W, F814W, WFC3 F125, F160W, IRAC Ch1 (3.6 µm) and Ch2 (4.5 µm), and MIRI F560W and F770W. The images above include the three galaxies with spectroscopic redshifts (ID 37653, 19180, and 6811). The images are 6" × 6", centered on the galaxy in each bandpass, as labeled along the top row.</figcaption></figure>

(The complete figure set (29 images) is available.)
increase with increasing redshift). Therefore, we fix the ionization parameter to a high, physically plausible value of $log U = -2$, as this is representative of the values used in previous studies when fitting the SEDs of galaxies (see, for example, the discussion in Whitler et al. 2023). We plan to explore how variations in the ionization parameter impact the constraints on the stellar populations using future data that can include spectroscopy, e.g., from JWST/NIRSpec.

The other exception is for galaxies with photometric redshifts, where we use the photometric redshift posterior, $P(z)$, derived from EAZY as the prior on the redshift (see Chworowsky et al. 2023). For galaxies with spectroscopic redshifts, we force the fit to the spectroscopic redshift value listed in Table 2.

For the galaxy star formation histories (SFHs), we test two possibilities. First, we primarily employ “delayed-$\tau$” models, where SFR $\propto (t/\tau) \times \exp(-t/\tau)$ for age, $t$, and star formation $e$-folding timescale, $\tau$. These models allow for SFHs that rise with time (when $t/\tau \ll 1$) (as is expected for high-redshift galaxies; Finlator et al. 2011; Papovich et al. 2011) as well as for exponentially declining models (when $t/\tau \gg 1$), and these have the flexibility to broadly reproduce the evolution of galaxies over long time periods (e.g., Larson & Tinsley 1978; Tinsley 1980; Carnall et al. 2019; García-Argumánez et al. 2023).

Second, following Papovich et al. (2001), we also consider more extreme SFHs that include both the delayed-$\tau$ model (above) and an early burst of stars that formed at $z_f = 100$. When a burst forms, the stellar population immediately begins aging. A burst at $z_f = \infty$ has the oldest possible age at any subsequent time, and the smallest amount of light at a given mass, and therefore it would have a maximal M/L at any observed epoch. As such, an early burst adds the maximum amount of stellar mass, with the minimum impact on the observed galaxy SED. In contrast, the slowly evolving delayed-$\tau$ model provides a fit to the light from the more recently formed stellar populations that dominate the rest-frame UV and optical light. Therefore, it is conceivably possible to hide significant amounts of stellar mass formed at earlier times that are “lost in the glare” of the luminous, more recently formed stars. This is also referred to as the outshining effect (e.g., Papovich et al. 2001; Dickinson et al. 2003; Papovich et al. 2006; Finkelstein et al. 2010; Maraston et al. 2010; Pforr et al. 2012; Conroy 2013).

The choice of $z_f = 100$ is motivated by the fact that current models expect stars to be forming by $z = 20–30$ Barkana & Loeb (2001), with galaxy candidates identified at $z \approx 15$ (see Section 1). A redshift of $z = 100$ is essentially “immediately” in the history of the Universe, as it corresponds to an age of only $\approx 17$ Myr after the Big Bang for the assumed cosmology. The time from $z = 100$ to $z = 20$ spans less than 200 Myr, during which little stellar evolution occurs for the longer-lived stars that dominate the stellar mass. Therefore, $z_f = 100$ seems to be a reasonable upper bound to ensure we capture the earliest time when stars could plausibly form. This formation redshift allows for a “maximally old” stellar population, and we constrain any star formation that may have occurred at the earliest times, which could have the highest possible M/L (for canonical stellar populations).<fnref n="35" />

By studying the SEDs of galaxies to longer wavelengths, we can constrain the amount of light in this population. For example, using K-band data, Papovich et al. (2001) and Dickinson et al. (2003) found that galaxies at $z \sim 3$ could hide as much as 75%-90% of their stellar mass in early bursts formed at $z_f = 100$. Indeed, at the risk of foreshadowing, we find that including the MIRI 5.6 and 7.7 µm data reduces the amount of possible stellar mass formed in such maximally old bursts by up to an order of magnitude (see Section 4.4).

# 4. Results

## 4.1. Analysis of Galaxy SEDs

Figure 3 shows the BAGPIPES SED fits and one-dimensional (1D) posterior likelihoods for select parameters of the SED fit. The figure shows six galaxies as an example. The online version of the paper includes a figure set with these plots for the full sample. For each galaxy, the plots compare the SED fits with and without the MIRI 5.6 and 7.7 µm data. Tables 4 and 5 provide the medians (50th percentiles), 16th and 84th percentiles derived from these posterior likelihoods for the stellar masses, SFRs, and photometric redshifts for all galaxies in the sample, both with and without using the MIRI data, respectively. Throughout, we take the “stellar mass” to be the “mass formed,” which is equivalent to the integral of the star formation history.

To test the robustness of the stellar masses and SFRs derived from the BAGPIPES fits, we have refit all the galaxies in our sample using several independent SED-fitting codes (CIGALE, Boquien et al. 2019; FAST, Kriek et al. 2009; and the codes of Santini et al. 2022 and Pérez-González et al. 2008). Comparing the stellar masses and SFRs, we find they agree in the mean (with bias, $\mu \approx 0$ dex) and an inter-method scatter of $\sigma = 0.23$ dex in stellar mass and $\sigma = 0.27$ dex in SFR. This scatter is typical in comparisons of SED-fitting results (e.g., Mobasher et al. 2015). We therefore interpret the scatter as representative of the systematic uncertainties on the stellar masses and SFRs here.

In some cases, adding the MIRI data has a small effect on the median values of the stellar mass and SFR, but it does tighten the allowable range of these parameters. Panels (a) and (b) of Figure 3 show galaxy ID 7364 (at $z = 8.1–8.2$) and galaxy ID 6811 (with $z_{sp} = 8.683$; Zitrin et al. 2015). These have MIRI data that support the interpretation inferred using only the HST and Spitzer data. However, in both of these cases, adding the MIRI data tightens the allowed range of models, and thus improves the constraints on the stellar population parameters. In both cases shown here, the favored range of stellar mass and SFR is improved significantly when MIRI data are included (with improvements in the inter-68th-percentile range by more than a factor of two). Below the SED fit for each galaxy, Figure 3 shows the posterior probability densities for the SFR, mass-weighted age (Age$_{MW}$), stellar mass, and dust attenuation. Adding the MIRI data typically produces narrower posteriors for SFR and stellar mass. This is a result of improved constraints on the dust attenuation ($A(V)$), and this forces the models to a narrower range of SFR and stellar mass. This is the case for 25% of the sample here (7 of the 28 galaxies), based on our visual inspection of the SEDs and 1D posteriors in Figure 3.

In other cases, adding the MIRI data changes the interpretation of the galaxy stellar populations dramatically. Panels (c)

<footnote n="35">³⁵ Using a lower formation redshift, $z_f < 100$, would lower the upper limit on the stellar masses that could form in the bursts as these would be younger, with lower M/L.</footnote>
<figure><img src="image_6.png" /></figure>

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

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

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

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

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

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

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

**Figure 3.** Examples of SED fits for galaxies in the sample. The green data points show the measured flux densities and uncertainties on the HST/WFC3, Spitzer/IRAC, and JWST/MIRI bands. The red-shaded regions show the model fit to all the data points (the shaded region shows the inner 68% range of models; the solid red line shows the median). The black-shaded region shows the fit to the data points excluding the MIRI bands. The small red and gray points show the median-model photometry. The inset text gives the median and 68% tile uncertainties on the stellar masses and SFRs inferred from the fits. Below each SED plot, the panels show the posteriors (probability density) for the SFR, mass-weighted (MW) age, stellar mass, and dust attenuation for each galaxy (using MIRI and excluding MIRI data). The dashed lines denote the 5%, 50%, and 95% intervals. The examples include galaxies where adding the MIRI data yields similar stellar masses and SFRs, but with tighter constraints (panels (a) and (b)). Other examples show galaxies where adding the MIRI data greatly reduces both the stellar masses and SFRs. Typically, this results from contamination in the IRAC data or because of strong emission lines impacting the IRAC data (or both; see panels (c) and (d)), or because the stellar continua appear very blue (see panels (e) and (f)).

(The complete figure set (28 images) is available.)
<figure><img src="image_14.png" /><figcaption>Figure 3 (Continued.)</figcaption></figure>

<table><caption>Table 4
Stellar Masses, SFRs, and Redshifts Derived Including the MIRI [5.6] and [7.7] Data</caption><thead><tr><th rowspan="2">ID<br>(1)</th><th colspan="3">Redshift</th><th colspan="3">log M<sub>*</sub>/M<sub>☉</sub> (dex)</th><th colspan="3">logSFR/M<sub>☉</sub> yr<sup>-1</sup> (dex)</th><th rowspan="2">"Burst" Mass<br>log M<sub>*</sub>/M<sub>☉</sub> (dex)<br>(11)</th></tr><tr><th>z<sub>50</sub><br>(2)</th><th>z<sub>16</sub><br>(3)</th><th>z<sub>84</sub><br>(4)</th><th>logM<sub>50</sub><br>(5)</th><th>logM<sub>16</sub><br>(6)</th><th>logM<sub>84</sub><br>(7)</th><th>log SFR<sub>50</sub><br>(8)</th><th>log SFR<sub>16</sub><br>(9)</th><th>log SFR<sub>84</sub><br>(10)</th></tr></thead><tbody><tr><td>5090</td><td>4.33</td><td>4.19</td><td>4.45</td><td>9.01</td><td>8.84</td><td>9.16</td><td>0.59</td><td>0.50</td><td>0.69</td><td>9.55</td></tr><tr><td>6811</td><td>...</td><td>...</td><td>...</td><td>9.76</td><td>9.57</td><td>9.88</td><td>1.60</td><td>1.49</td><td>1.72</td><td>10.34</td></tr><tr><td>7364</td><td>8.07</td><td>7.35</td><td>8.55</td><td>9.91</td><td>9.72</td><td>10.04</td><td>1.70</td><td>1.57</td><td>1.82</td><td>10.61</td></tr><tr><td>7600</td><td>4.41</td><td>4.22</td><td>4.56</td><td>10.54</td><td>10.36</td><td>10.65</td><td>2.09</td><td>1.97</td><td>2.20</td><td>10.93</td></tr><tr><td>7818</td><td>5.19</td><td>4.93</td><td>5.39</td><td>9.51</td><td>9.30</td><td>9.64</td><td>1.19</td><td>1.10</td><td>1.29</td><td>10.23</td></tr><tr><td>11329</td><td>4.45</td><td>4.31</td><td>4.59</td><td>8.48</td><td>8.14</td><td>8.73</td><td>0.27</td><td>0.17</td><td>0.37</td><td>9.17</td></tr><tr><td>12514</td><td>7.69</td><td>7.03</td><td>8.21</td><td>8.78</td><td>8.47</td><td>9.02</td><td>0.70</td><td>0.53</td><td>0.85</td><td>9.45</td></tr><tr><td>12773</td><td>5.08</td><td>4.91</td><td>5.24</td><td>8.01</td><td>7.78</td><td>8.24</td><td>0.08</td><td>-0.18</td><td>0.30</td><td>8.87</td></tr><tr><td>13179</td><td>4.79</td><td>4.73</td><td>4.86</td><td>10.07</td><td>9.86</td><td>10.20</td><td>1.71</td><td>1.60</td><td>1.83</td><td>10.57</td></tr><tr><td>13389</td><td>4.53</td><td>4.45</td><td>4.60</td><td>9.61</td><td>9.42</td><td>9.75</td><td>1.25</td><td>1.15</td><td>1.37</td><td>10.16</td></tr><tr><td>14178</td><td>4.89</td><td>4.76</td><td>5.08</td><td>8.54</td><td>8.29</td><td>8.79</td><td>0.54</td><td>0.35</td><td>0.67</td><td>9.16</td></tr><tr><td>15445</td><td>4.59</td><td>4.46</td><td>4.71</td><td>8.92</td><td>8.67</td><td>9.11</td><td>0.63</td><td>0.54</td><td>0.73</td><td>9.36</td></tr><tr><td>18441</td><td>6.54</td><td>5.90</td><td>7.11</td><td>9.32</td><td>9.12</td><td>9.46</td><td>1.03</td><td>0.90</td><td>1.17</td><td>10.00</td></tr><tr><td>18449</td><td>4.97</td><td>4.43</td><td>5.32</td><td>9.17</td><td>8.93</td><td>9.32</td><td>0.81</td><td>0.67</td><td>0.95</td><td>9.83</td></tr><tr><td>19180</td><td>...</td><td>...</td><td>...</td><td>9.72</td><td>9.48</td><td>9.90</td><td>1.54</td><td>1.38</td><td>1.68</td><td>10.14</td></tr><tr><td>24007</td><td>5.36</td><td>5.17</td><td>5.52</td><td>8.76</td><td>8.48</td><td>8.91</td><td>0.41</td><td>0.32</td><td>0.53</td><td>9.59</td></tr><tr><td>26890</td><td>8.80</td><td>8.61</td><td>8.97</td><td>8.79</td><td>8.46</td><td>9.04</td><td>0.84</td><td>0.50</td><td>0.97</td><td>9.45</td></tr><tr><td>34813</td><td>4.38</td><td>4.27</td><td>4.50</td><td>8.26</td><td>7.97</td><td>8.62</td><td>0.33</td><td>0.02</td><td>0.46</td><td>8.97</td></tr><tr><td>35896</td><td>5.13</td><td>4.82</td><td>5.40</td><td>8.92</td><td>8.67</td><td>9.08</td><td>0.63</td><td>0.51</td><td>0.76</td><td>9.44</td></tr><tr><td>37653</td><td>...</td><td>...</td><td>...</td><td>8.95</td><td>8.58</td><td>9.13</td><td>0.66</td><td>0.53</td><td>0.83</td><td>9.57</td></tr><tr><td>37703</td><td>4.63</td><td>4.53</td><td>4.71</td><td>8.41</td><td>8.20</td><td>8.86</td><td>0.47</td><td>0.23</td><td>0.65</td><td>9.21</td></tr><tr><td>39096</td><td>6.78</td><td>6.68</td><td>6.87</td><td>8.81</td><td>8.48</td><td>8.98</td><td>0.63</td><td>0.52</td><td>0.75</td><td>9.41</td></tr><tr><td>41375</td><td>4.79</td><td>4.66</td><td>4.96</td><td>8.54</td><td>8.27</td><td>8.77</td><td>0.51</td><td>0.34</td><td>0.61</td><td>9.15</td></tr><tr><td>41545</td><td>5.13</td><td>5.00</td><td>5.25</td><td>9.80</td><td>9.62</td><td>9.89</td><td>1.38</td><td>1.29</td><td>1.49</td><td>10.32</td></tr><tr><td>41564</td><td>4.66</td><td>4.59</td><td>4.74</td><td>9.36</td><td>9.08</td><td>9.52</td><td>1.05</td><td>0.93</td><td>1.17</td><td>9.78</td></tr><tr><td>42638</td><td>4.70</td><td>4.64</td><td>4.76</td><td>9.65</td><td>9.39</td><td>9.80</td><td>1.34</td><td>1.20</td><td>1.45</td><td>10.12</td></tr><tr><td>45145</td><td>4.72</td><td>4.55</td><td>4.96</td><td>8.59</td><td>8.28</td><td>8.83</td><td>0.51</td><td>0.34</td><td>0.61</td><td>9.52</td></tr><tr><td>49365</td><td>5.31</td><td>5.14</td><td>5.52</td><td>9.63</td><td>9.43</td><td>9.75</td><td>1.29</td><td>1.18</td><td>1.42</td><td>10.23</td></tr></tbody></table>

Note. (1) Galaxy ID; (2)–(4) redshift median (50th percentile), 16th percentile, and 84th percentile; galaxies with no redshift use the spectroscopic redshift in Table 2; (5)–(7) stellar mass median (50th percentile), 16th percentile, and 84th percentile in the delayed-$\tau$ model; (8)–(10) SFR median (50th percentile), 16th percentile, and 84th percentile (all SFRs are averaged over the past 100 Myr) in the delayed-$\tau$ model; (11) maximum stellar mass allowed in the “burst” formed at $z_f = 100$, where these satisfy the BIC criteria (Equation (2) in Section 4.4) and correspond approximately to a $3\sigma$ upper limit.
<table><caption>**Table 5**
Stellar Masses, SFRs, and Redshifts Derived Excluding the MIRI Data</caption><thead><tr><th rowspan="2">ID<br>(1)</th><th colspan="3">Redshift</th><th colspan="3">log M<sub>k</sub>/M<sub>☉</sub> (dex)</th><th colspan="3">logSFR/M<sub>☉</sub> yr<sup>-1</sup> (dex)</th><th rowspan="2">"Burst" Mass<br>log M<sub>k</sub>/M<sub>☉</sub> (dex)<br>(11)</th></tr><tr><th>z<sub>50</sub><br>(2)</th><th>z<sub>16</sub><br>(3)</th><th>z<sub>84</sub><br>(4)</th><th>log M<sub>50</sub><br>(5)</th><th>log M<sub>16</sub><br>(6)</th><th>log M<sub>84</sub><br>(7)</th><th>log SFR<sub>50</sub><br>(8)</th><th>log SFR<sub>16</sub><br>(9)</th><th>log SFR<sub>84</sub><br>(10)</th></tr></thead><tbody><tr><td>5090</td><td>4.32</td><td>4.17</td><td>4.46</td><td>9.08</td><td>8.75</td><td>9.26</td><td>0.61</td><td>0.50</td><td>0.74</td><td>10.00</td></tr><tr><td>6811</td><td>-1.00</td><td>-1.00</td><td>-1.00</td><td>9.85</td><td>9.47</td><td>10.12</td><td>1.70</td><td>1.43</td><td>1.98</td><td>11.17</td></tr><tr><td>7364</td><td>8.23</td><td>7.42</td><td>8.61</td><td>9.77</td><td>9.48</td><td>9.97</td><td>1.60</td><td>1.38</td><td>1.80</td><td>10.85</td></tr><tr><td>7600</td><td>4.42</td><td>4.26</td><td>4.57</td><td>10.53</td><td>10.32</td><td>10.68</td><td>2.10</td><td>1.97</td><td>2.24</td><td>11.17</td></tr><tr><td>7818</td><td>4.83</td><td>4.67</td><td>5.17</td><td>9.83</td><td>9.57</td><td>10.02</td><td>1.45</td><td>1.28</td><td>1.61</td><td>10.51</td></tr><tr><td>11329</td><td>4.40</td><td>4.24</td><td>4.56</td><td>8.82</td><td>8.50</td><td>9.05</td><td>0.40</td><td>0.24</td><td>0.55</td><td>9.84</td></tr><tr><td>12514</td><td>7.38</td><td>6.77</td><td>8.13</td><td>9.30</td><td>8.92</td><td>9.64</td><td>1.10</td><td>0.77</td><td>1.44</td><td>10.41</td></tr><tr><td>12773</td><td>4.97</td><td>4.83</td><td>5.15</td><td>8.96</td><td>8.55</td><td>9.22</td><td>0.58</td><td>0.41</td><td>0.78</td><td>10.23</td></tr><tr><td>13179</td><td>4.76</td><td>4.71</td><td>4.81</td><td>10.25</td><td>10.03</td><td>10.38</td><td>1.81</td><td>1.70</td><td>1.96</td><td>10.89</td></tr><tr><td>13389</td><td>4.52</td><td>4.44</td><td>4.59</td><td>9.82</td><td>9.65</td><td>9.97</td><td>1.37</td><td>1.24</td><td>1.51</td><td>10.41</td></tr><tr><td>14178</td><td>4.78</td><td>4.65</td><td>4.92</td><td>9.27</td><td>8.99</td><td>9.45</td><td>0.90</td><td>0.72</td><td>1.07</td><td>10.08</td></tr><tr><td>15445</td><td>4.58</td><td>4.45</td><td>4.70</td><td>9.03</td><td>8.70</td><td>9.25</td><td>0.68</td><td>0.55</td><td>0.85</td><td>10.07</td></tr><tr><td>18441</td><td>6.61</td><td>6.07</td><td>6.87</td><td>10.18</td><td>9.90</td><td>10.42</td><td>1.94</td><td>1.66</td><td>2.16</td><td>11.11</td></tr><tr><td>18449</td><td>5.00</td><td>4.43</td><td>5.41</td><td>9.24</td><td>8.92</td><td>9.47</td><td>0.90</td><td>0.64</td><td>1.10</td><td>10.13</td></tr><tr><td>19180</td><td>-1.00</td><td>-1.00</td><td>-1.00</td><td>10.52</td><td>10.31</td><td>10.70</td><td>2.09</td><td>1.90</td><td>2.25</td><td>11.52</td></tr><tr><td>24007</td><td>5.28</td><td>5.06</td><td>5.47</td><td>8.93</td><td>8.59</td><td>9.27</td><td>0.58</td><td>0.38</td><td>0.88</td><td>10.40</td></tr><tr><td>26890</td><td>8.92</td><td>8.67</td><td>9.17</td><td>9.39</td><td>9.10</td><td>9.67</td><td>1.28</td><td>1.05</td><td>1.51</td><td>10.53</td></tr><tr><td>34813</td><td>4.34</td><td>4.23</td><td>4.47</td><td>8.78</td><td>8.40</td><td>9.05</td><td>0.48</td><td>0.39</td><td>0.59</td><td>9.64</td></tr><tr><td>35896</td><td>4.88</td><td>4.66</td><td>5.20</td><td>9.23</td><td>8.94</td><td>9.47</td><td>0.93</td><td>0.72</td><td>1.12</td><td>10.08</td></tr><tr><td>37653</td><td>-1.00</td><td>-1.00</td><td>-1.00</td><td>9.55</td><td>9.26</td><td>9.75</td><td>1.18</td><td>1.00</td><td>1.36</td><td>10.39</td></tr><tr><td>37703</td><td>4.58</td><td>4.48</td><td>4.68</td><td>9.16</td><td>8.90</td><td>9.35</td><td>0.75</td><td>0.62</td><td>0.88</td><td>9.94</td></tr><tr><td>39096</td><td>6.73</td><td>6.59</td><td>6.84</td><td>9.04</td><td>8.68</td><td>9.46</td><td>0.80</td><td>0.57</td><td>1.21</td><td>10.03</td></tr><tr><td>41375</td><td>4.73</td><td>4.60</td><td>4.87</td><td>9.15</td><td>8.79</td><td>9.34</td><td>0.73</td><td>0.58</td><td>0.91</td><td>10.19</td></tr><tr><td>41545</td><td>5.12</td><td>4.98</td><td>5.26</td><td>9.85</td><td>9.61</td><td>10.03</td><td>1.45</td><td>1.31</td><td>1.60</td><td>10.49</td></tr><tr><td>41564</td><td>4.64</td><td>4.57</td><td>4.71</td><td>9.56</td><td>9.34</td><td>9.72</td><td>1.15</td><td>1.04</td><td>1.29</td><td>10.18</td></tr><tr><td>42638</td><td>4.68</td><td>4.62</td><td>4.75</td><td>9.85</td><td>9.62</td><td>10.01</td><td>1.47</td><td>1.32</td><td>1.64</td><td>10.55</td></tr><tr><td>45145</td><td>4.66</td><td>4.52</td><td>4.86</td><td>9.08</td><td>8.71</td><td>9.30</td><td>0.69</td><td>0.53</td><td>0.88</td><td>9.98</td></tr><tr><td>49365</td><td>5.43</td><td>5.23</td><td>5.63</td><td>9.45</td><td>9.14</td><td>9.65</td><td>1.09</td><td>0.89</td><td>1.29</td><td>10.31</td></tr></tbody></table>

Note. (1) Galaxy ID, (2)–(4) redshift median (50th percentile), 16th percentile, and 84th percentile; galaxies with no redshift use the spectroscopic redshift in Table 2; (5)–(7) stellar mass median (50th percentile), 16th percentile, and 84th percentile in the delayed-$\tau$ model; (8)–(10) SFR median (50th percentile), 16th percentile, and 84th percentiles (all SFRs are averaged over the past 100 Myr) in the delayed-$\tau$ model; (11) maximum stellar mass allowed in the “burst” formed at $z_f = 100$, where these satisfy the BIC criteria (Equation (2) in Section 4.4) and correspond approximately to a $3\sigma$ upper limit.

and (d) of Figure 3 show galaxies with ID 19180 (with $z_{sp} = 5.077$) and ID 18441 (at $z = 6.5–6.6$). In both of these cases, without MIRI data, the HST to IRAC 3.6 and 4.5 µm data implied very red rest-UV-to-optical colors, leading to high dust attenuation values ($A(V) \gtrsim 1$ mag) with high SFRs ($\gtrsim 90–100$ M$_{\odot}$ yr$^{-1}$). The stellar masses in these cases are also elevated primarily because the higher dust attenuation increases the mass-to-light ratio (M/L) of the models, and therefore increases the stellar mass. Including the MIRI data changes the favored stellar population models to ones with much bluer rest-UV-to-optical colors. As a result, the dust attenuation, SFR, and stellar mass are decreased—in some cases, by an order of magnitude. The decrease in stellar mass is more than that of the SFR, implying the specific SFR increases slightly as well. This impacts 36% of the sample here (10 of the 28 galaxies).

Yet in other cases, adding the MIRI data forces the constraints on the stellar populations to be bluer than expected based on the HST and Spitzer data. Panels (e) and (f) of Figure 3 show two galaxies that demonstrate these situations. In both the cases of galaxy ID 12514 (at $z = 7.4–7.7$) and ID 26890 (at $z = 8.8–8.9$), the MIRI data favor very blue UV/optical colors. In the case of galaxy 12514, there are indications that the IRAC and MIRI 5.6 µm data are boosted by the presence of strong emission lines. Having the 7.7 µm

photometry favors a lower stellar continuum. As a result, the SFR and stellar mass are reduced (the presence of redshifted $H\alpha$ in the MIRI F560W band is apparent even in the galaxy image in Figure 2, which shows the 5.6 µm flux density to be noticeably brighter than the 7.7 µm flux density). For this reason, the SED fit favors a slightly higher photometric redshift to accommodate the $H\alpha$ emission in the F560W bandpass; see Section 4.2. This is the case for 39% of the sample here (11 of the 28 galaxies). In the latter category, galaxy ID 26890 is noteworthy in itself because it is the highest-redshift galaxy in the sample, and because the HST/WFC3-to-JWST/MIRI colors are $H_{160} – [5.6] \approx –0.5$ mag and indicative of stellar populations with very low M/L ratios. The IRAC 4.5 µm emission shows indications of enhancement, possibly owing to redshifted $H\alpha+[O III]$. The MIRI data rein in the allowable range of stellar population parameters, favoring models with lower SFRs than the constraints that lack MIRI data.

Therefore, the MIRI data favor bluer rest-UV/optical colors compared to the IRAC data. In part, this may have a physical explanation. At the redshifts of the galaxies under consideration, the IRAC data contain strong emission lines, including $H\alpha+[O III]$ at $5 < z < 8$, $[O II]+[Ne III]$ at $7 < z < 12$, and $H\alpha +[N II]$ at $4 < z < 7$ (Labbé et al. 2013; Smit et al. 2015; De Barros et al. 2019; Roberts-Borsani et al. 2020; Endsley et al.
<figure><img src="image_15.png" /><figcaption>**Figure 4.** Comparison of the redshifts for galaxies in our sample with and without the MIRI data. The data points show the relative difference between the photometric redshifts (derived from our *BAGPIPES* fits) for the galaxies excluding MIRI ($z_{\text{no MIRI}}$) and when including MIRI ($z_{\text{MIRI}}$) as a function of the prior photometric redshift from Finkelstein et al. (2022a). Annotated points indicate objects with $|z_{\text{no MIRI}} - z_{\text{MIRI}}| > 0.2$. These objects indicate shifts in their photometric redshift, where this is at least in part because of nebular emission in one or more bands.</figcaption></figure>

2021). At certain redshifts, both the 3.6 and 4.5 µm bands can be affected (e.g., Smit et al. 2014; Roberts-Borsani et al. 2020). Without the SED fits, the models that fit the data may include both models with strong emission lines and/or models with redder colors (indicative of strong dust attenuation) or both. The models in Figure 3 show that the presence of strong emission lines in the bands augments the flux densities. This could account for part of the difference between the SED fits with and without MIRI: the MIRI data provide evidence in these cases that the stellar populations have bluer colors, and the elevated IRAC flux densities then would require strong emission line EWs to reproduce them.

However, crowding between sources in the IRAC data may be another reason for the increased IRAC flux densities. The lower angular resolution of the IRAC data can cause blending from bright, nearby galaxies, and this can lead to additional uncertainties in the flux density. As the images in Figure 2 illustrate, some galaxies (e.g., ID 19180 and 6811) have bright objects within 3". The light from these objects makes deblending more difficult and could potentially bias the flux-density measurements (see, e.g., Laidler et al. 2007; Skelton et al. 2014; Merlin et al. 2016). For this reason, the elevated IRAC flux densities may include systematic measurement uncertainties from blended objects. In the Appendix, we test for the effect of source blending by excluding objects that have any bright neighbor within 3" (we show that source blending does not bias our interpretation for the stellar masses and SFRs, nor their evolution, that we infer for the galaxy population). However, this does emphasize the benefit of having data with the enhanced angular resolution available from JWST imaging.

## 4.2. The Impact of MIRI on Redshifts and Implications for Emission Lines

Figure 4 compares the redshifts for the galaxies in our sample derived from our *BAGPIPES* fits for the galaxies excluding the MIRI data ($z_{\text{no MIRI}}$) and when including the MIRI data ($z_{\text{MIRI}}$) as a function of the prior photometric redshift

from Finkelstein et al. (2022a). For most galaxies, there is good agreement. For completeness, we include all 28 galaxies in this comparison, using the photometric redshifts even for galaxies with spectroscopic redshifts.

In several instances, we see the median redshift from the posterior shifts appreciably when including the MIRI data (with $\Delta z > 0.2$; these are annotated in Figure 4). In all cases, the shift in the redshift probability density appears to be related to the effects of one or more emission lines in the IRAC or MIRI passbands. Figure 5 illustrates the shift in $P(z)$ from the *BAGPIPES* fits for these galaxies.

Two of these galaxies lie at $4 < z < 6$ (galaxy IDs 35896 and 7818). Adding the MIRI data favors having strong nebular emission in both of the IRAC bands. This occurs at $z \gtrsim 5$ when $H\alpha+[O III]$ enters the IRAC Ch1 bandpass (at 3.6 µm) and $H\alpha$ enters the IRAC Channel 2 bandpass (at 4.5 µm; see Figure 1). Inspection of the galaxy SEDs (see Figure Set. 3) shows that the IRAC-to-MIRI colors ([3.6]–[5.6] and [4.5]–[5.6]) are blue for both of these galaxies. The *BAGPIPES* fits that include the MIRI data favor models in which the galaxy has strong nebular emission in the IRAC bands to account for the observed color. This increases (decreases) the redshift probability density at higher (lower) redshifts for these galaxies.

The remaining galaxy, ID 12514, lies at $z \approx 7.7$. Figure 5 shows the redshift probability densities for this galaxy. The MIRI F560W image shows evidence of “boosted” flux compared to the MIRI F770W image (see the images in Figure 2), where the MIRI color is [5.6]–[7.7] = $-1.1 \pm 0.3$ mag (see Table 2). The most likely explanation for this boosted emission appears to be from redshifted $H\alpha+[N II]$ at $z = 7.6–8.0$ in the MIRI F560W bandpass (see Figure 3). The strength of the flux in the F560W bandpass decreases the probability density for redshifts with $z \lesssim 7$, as these would place $H\alpha$ at wavelengths shorter than those covered by the bandpass (see Figure 1).

Therefore, nebular emission lines appear to be responsible for the larger shifts in the photometric redshifts. However, these cases are generally exceptions. For most galaxies in the sample, the redshift posteriors are consistent; Figure 4 shows the differences are negligible in median redshifts with and without including the MIRI data. We expect that many (even most) of the galaxies in the sample also exhibit strong emission lines, given the preponderance of evidence from other galaxies at these redshifts (see Section 4.1) and by inspection of the SEDs in Figure 3. In these cases, adding the MIRI data supports the redshifts, either because the emission lines make less of an impact on the redshift likelihoods or because the MIRI data reinforce them.

## 4.3. The Impact of MIRI on the Stellar Masses, Dust Attenuation, and SFRs

Figure 6 compares the stellar masses and SFRs derived for the galaxies in our sample using the simple delayed-$\tau$ models with *BAGPIPES* for the case where we include the MIRI F560W and F770W data versus the case where we exclude it. In general, including the MIRI data reduces the stellar masses and SFRs for the galaxies.

We study the offsets in stellar mass and SFR for our galaxies in two redshift bins, $4 < z < 6$ and $6 < z < 10$, where the median offsets in stellar mass and SFR are indicated in the Figure 6 as large rectangles, and are listed in Table 6. For this comparison, we
<figure><img src="image_16.png" /><figcaption>**Figure 5.** Comparison of redshift posterior probability densities, $P(z)$, derived from the BAGPIPES fits to three of the galaxies in the sample (as labeled). These three galaxies each have a change of more than 0.2 between the respective median redshifts when including the MIRI F560W and F770W data in the fit vs. when they are excluded. For the case of galaxy IDs 35896 and 7818, the IRAC–MIRI colors are blue, indicating redshifted emission lines (likely $H\alpha$+[O III] and $H\alpha$) inhabit both IRAC bands while the MIRI data probe the galaxy continua. In the case of galaxy ID 12514, the MIRI F560W band clearly shows boosted emission, likely from redshifted $H\alpha$ at $z = 7.6–8$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 6.** Comparison of stellar masses and SFRs derived from the SED modeling for galaxies including the MIRI F560W and F770W data and those without the MIRI data. The top set of panels show the comparison for the stellar mass. The bottom set of panels show the comparison for the SFR. In each row, the left panel shows the direct comparison, where the dashed line shows the one-to-one relation. The right panel shows the logarithmic difference, defined as $\Delta \log M_*^{\text{with MIRI}} - \log M_*^{\text{without MIRI}}$ (and similarly for the SFR). The symbols are color-coded by redshift (using the right panel). In the right panel, the large rectangular boxes and error bars show the median value (and scatter) in two bins of redshift ($4 < z < 6$ and $6 < z < 10$). The values for these are given in Table 6. Adding the MIRI data generally lowers the SFRs and stellar masses of these galaxies, though the scatter is significant. The median offset is larger for galaxies at higher redshifts.</figcaption></figure>
<table><caption>**Table 6**
Offsets in Stellar Mass and SFRs for 4 < z < 10 Galaxies when Including the MIRI F560W and F770W Data</caption><thead><tr><th rowspan="3">Sample</th><th colspan="2">Stellar mass offsets</th><th colspan="2">SFR offsets</th></tr><tr><th colspan="2">&lang;&Delta; log M<sub>*</sub>&rang; =</th><th colspan="2">&lang;&Delta;logSFR&rang; =</th></tr><tr><th>Median</th><th>Scatter</th><th>log SFR<sub>no MIRI</sub> / SFR<sub>MIRI</sub></th><th>Median</th><th>Scatter</th></tr></thead><tbody><tr><td>4 &lt; z &lt; 6</td><td>0.25 dex</td><td>0.28 dex</td><td>0.15 dex</td><td>0.12 dex</td></tr><tr><td>6 &lt; z &lt; 10</td><td>0.38 dex</td><td>0.44 dex</td><td>0.29 dex</td><td>0.27 dex</td></tr></tbody></table>

**Notes.** The quantities with the subscript “noMIRI” denote values derived without the MIRI data. Quantities with the subscript “MIRI” denote the values derived with the 5.6 and 7.7 µm MIRI data. The scatter is the inter-68th-percentile interval derived from the median absolute deviation.

have defined $\Delta \log M_* = \Delta \log M_*^{\text{with MIRI}} - \log M_*^{\text{without MIRI}}$ and $\Delta \log \text{SFR} = \log \text{SFR}^{\text{with MIRI}} - \log \text{SFR}^{\text{without MIRI}}$. We have added their uncertainties in quadrature, but this assumes the stellar masses and SFRs derived with and without MIRI are independent. Strictly speaking, this is not true, as both values use the same HST and Spitzer/IRAC data. Therefore, it is likely that the error bars are overestimated, assuming the values are correlated. Nevertheless, this does not have a major impact on our work, because we are interested in the medians and scatter of the sample.

Formally, the offsets in stellar mass are $\Delta \log(M_*) = 0.25$ dex for $4 < z < 6$ and 0.38 dex for $6 < z < 10$. The inter-68th-percentile intervals are 0.28 and 0.44 dex, respectively (estimated using the median absolute deviation, $\sigma_{\text{NMAD}}$). For the SFRs, the offsets are $\Delta \log(\text{SFR}) = 0.15$ dex for $4 < z < 6$ and 0.29 dex for $6 < z < 10$, with an inter-68th-percentile interval of 0.12 dex and 0.27 dex, respectively. Because the impact of MIRI is larger on the stellar mass than the SFR, the specific SFRs will be increased by approximately $\Delta \log \text{SFR} - \Delta \log M_* \approx 0.1$ dex (see Table 6).

The reasons for the offsets are similar to those discussed for the individual galaxies in Section 4.1. Including the MIRI data generally favors stellar populations with bluer rest-UV/optical colors, with lower M/L ratios. This forces the constraints to lower stellar-mass models. For the SFRs, because the colors are bluer, there is less dust attenuation favored in the models, which lowers the SFRs compared to models with higher dust attenuation (at fixed observed galaxy luminosity). The MIRI data also remove some of the degeneracy between models with redder stellar populations versus those with strong emission lines impact select bands. These have the combined effect of favoring lower stellar masses and SFRs when including the MIRI data.

There is significant scatter in the offsets of stellar mass and SFR for individual objects. In Figure 6, the error bars on the large rectangles show the inter-68th-percentile range (i.e., the difference between the 16th and 84th percentiles) for both the $4 < z < 6$ and $6 < z < 10$. For some galaxies, the offsets are insignificant, with $\Delta \log M_* \approx 0$ and $\Delta \log \text{SFR} \approx 0$. In other words, in some cases, the MIRI data support the range of stellar population parameters favored by the fits to the HST/ACS + WFC3 and Spitzer/IRAC data, but the MIRI data improve the accuracy of the measurements, typically by a factor of order two.

In most cases, the MIRI data reduce the stellar mass and SFRs—in some cases, substantially so. Figure 6 shows that

this is generally the case, where the MIRI data decrease the average stellar mass and SFR for galaxies in our sample, typically by a factor of order two. We will explore the implications this has on the evolution of the galaxy stellar-mass density in Section 5.

### 4.4. The Impact of MIRI on the Inferred Star Formation History (and the Mass Formed in Bursts)

Using the results from Section 3, we also explored the possibility that an early burst of star formation could contribute stellar mass to the galaxies in our sample, as well as how the JWST/MIRI data can improve the constraints on this population. For this study, we used the SED fits to the galaxies in our sample, using a star formation history that includes both the delayed-$\tau$ model and the burst of stars formed at $z_f = 100$. The details of the other model parameters are listed in Table 3.

To quantify the amount of allowed stellar mass formed in the burst for each galaxy, we use the Bayesian Information Criterion (BIC; Bailer-Jones 2017). The BIC provides a criterion for model selection in the case that one introduces a model with an additional parameter. In our case, we are selecting between two models: one with a star formation history that has a delayed-$\tau$ model only, and one with both the delayed-$\tau$ model and an early burst at $z_f = 100$. When comparing the two models, the BIC applies a penalty term for the additional parameter to determine if the additional parameter improves the fit or is overfitting the data.

We use the BIC defined as (e.g., Buat et al. 2019)

$$ \text{BIC} = \chi_0^2 \times k \ln N, \quad (2) $$

where $\chi_0^2$ is the minimum chi-squared from the SED-fitting, k is the number of independently fitted parameters (we use $k=7$), and N is the number of data points ($N=8$ or 10, depending on if the MIRI [5.6] and [7.7] data are included or not).

Our goal is to quantify the upper limit on the stellar mass formed in an early burst that could exist in our galaxies. To do this, we select models with bursts ($z_f = 100$) that are not excluded by the BIC criterion. That is, for a given galaxy, we identify all models that satisfy $\chi^2 \le BIC$, where the BIC is defined in Equation (2) and $\chi^2$ is the fit for a given model (this is similar to the approach adopted by Buat et al. 2019). For each galaxy, we select the model with the highest stellar mass in the $z_f = 100$ burst from the subset of models that satisfy the BIC. We then take this as an upper limit on the stellar mass permitted in the burst. Comparing this upper limit on the stellar mass formed in bursts to the range of stellar masses from the fits, we find that the limiting stellar mass from models that satisfy the BIC criteria corresponds roughly to a 99.7% upper limit on the mass (i.e., the values we report correspond approximately to a $3\sigma$ upper limit on the stellar mass). Tables 4 and 5 list the upper limit on the stellar mass in the burst component for all galaxies in the sample for the cases where we include and exclude the MIRI data, respectively.

Figure 7 shows example SED fits for galaxies in our sample, both with and without the bursts, for both the cases where we include or exclude the MIRI 5.6 and 7.7 µm data. The complete figure set (28 images) is available online. The examples shown in the figure include a galaxy where the IRAC data are bright relative to the MIRI data (ID 18441). In this case, without the
<figure><img src="image_18.png" /><figcaption>**Figure 7.** Example fits to SEDs for the galaxies in our sample, comparing the results from fits that include an early burst of star formation at $z_f = 100$. Data points and upper limits have the same definitions as in Figure 3. In each panel, the shaded regions show the stellar population fit to the SED using the total model (cyan-shaded = delayed-$\tau$ plus the burst) and delayed-$\tau$ model only (gray-shaded; these are identical to those in Figure 3, and there is almost no difference between the cyan- and gray-shaded models). The tan-shaded region shows the light permitted in the burst component. The labels indicate the amount of stellar mass in each component. Each row shows the results for one galaxy, where the left panel shows the results that exclude the MIRI F560W and F770W data, and the right panel shows the results including the MIRI F560W and F770W data.</figcaption></figure>

(The complete figure set (28 images) is available.)
<figure><img src="image_19.png" /><figcaption>**Figure 8.** Comparison of stellar masses derived for galaxies with and without early bursts (at $z_f = 100$). The top row shows the results that lack MIRI data. The top left panel shows the stellar masses derived from the delay-$\tau$ models only (labeled “no burst”) compared to the models that include the burst (labeled “allowed in burst”). The top right panel shows the difference between the stellar masses as a function of redshift. The bottom row shows the same results for the galaxies including the MIRI data. In the left panels, the dashed lines show the one-to-one relation and the solid lines show the median offsets. In the right panels, the large rectangles show the medians in the two bins of redshift ($4 < z < 6$ and $6 < z < 9$) these are given in Table 7. Adding the MIRI data reduces the amount of stellar mass allowed in the burst components.</figcaption></figure>

MIRI data, the allowed stellar mass in the burst can reach $\log M_*/M_\odot = 11.0$, but adding MIRI reduces this by nearly an order of magnitude. For galaxy ID 7364, the MIRI data favor red IRAC-to-MIRI colors. Nevertheless, because the MIRI data constrain the models at longer wavelengths, they also lower the amount of stellar mass allowed in the burst: without the MIRI data, the burst can include $\log M_*/M_\odot = 10.9$; when MIRI data are included, the mass in the burst declines by 0.3 dex (a factor of two). For galaxy ID 6811, the MIRI data show that, because they constrain the SED at longer wavelengths, the amount of stellar mass allowed in the burst is reduced by a factor of 0.6 dex (nearly a factor of five).

The addition of the MIRI F560W and F770W data reduces the amount of stellar mass that can form in early bursts.

Figure 8 shows the change in stellar mass for the case that the star formation histories include only delayed-$\tau$ models (labeled “no burst” in the figure) compared to when an early burst of star formation at $z_f = 100$ is included (labeled “allowed in burst” in the figure). Figure 8 shows the results for both the case that the MIRI data are excluded (top row) and when the MIRI data are included (bottom row). For the fits that lack MIRI data, the amount of stellar mass allowed in the burst is nearly an order of magnitude higher than constrained in the delayed-$\tau$ models: in this case, the median differences in the log of the stellar mass of the models with early bursts and those with only delayed-$\tau$ models is 0.9 dex at $4 < z < 6$ and 1.1 dex at $6 < z < 9$. For the fits that include the MIRI data, the amount of stellar mass allowed in the bursts is significantly reduced: the
<figure><img src="image_20.png" /><figcaption>**Figure 9.** Comparison of the median and relative SED for each galaxy, both for the case that MIRI data are used (top panel) and when the MIRI data are excluded (middle panel). The individual lines are the median SED model fit to each galaxy in the sample, shifted to the rest-frame. The large data points show the median rest-frame magnitude at 1600 Å, 2800 Å, and *U*, *B*, *V*, and *I*. The error bars show the scatter in the sample. All models have been normalized to the 1600 Å magnitude (which accounts for the lack of scatter at that wavelength). The bottom panel shows the difference in color ($\Delta m$) between the models with and without the MIRI data. The error bars show the range of the 16th–84th percentiles of the sample. The change in the color at the reddest wavelengths probed (about the rest-frame *I*-band) corresponds to a $\Delta m \approx 0.4$ mag.</figcaption></figure>

median differences in this case are 0.6 dex at $4 < z < 6$ and 0.7 dex at $6 < z < 9$. These values are listed in Table 7.

# 5. Discussion

## 5.1. On the Colors, Stellar Masses, and Nebular Emission in Early Galaxies

One of main findings in this paper is that the MIRI data favor bluer colors in galaxies at $4 < z < 9$. Figure 9 shows this by comparing the relative SED for each galaxy, both in the case of including the MIRI data and that without the MIRI data. Including the MIRI data reduces the derived rest-frame *I*-band light by approximately $\Delta m_{1600\AA} - I \approx 0.4$ mag. Upon inspecting Figure 9, this appears to result from many galaxies favoring bluer SEDs when the MIRI data are included. In other words, without the MIRI data, the SED is unconstrained at longer wavelengths, and this allows for a greater range of SED shape (where the median favors a solution that is redder on average). Adding the MIRI data shifts the likelihood to bluer populations for many galaxies. This has a major impact on the implied M/

L, as the blue rest-frame colors imply younger ages, lower dust attenuation, or both. The fact that adding the MIRI data makes the galaxies bluer largely explains the differences in the derived stellar masses and SFRs observed in Figure 6, where adding the MIRI data lowers the stellar masses and SFRs compared to when the MIRI data are excluded.

Therefore, our interpretation of the MIRI data is that galaxies at high redshifts ($z > 4$) are bluer than inferred from previous studies. This adds to other studies that find that galaxies at high redshifts must have (very) blue colors. Studies from the pre-JWST era argued that galaxies at $z > 4$ show indications of declining (i.e., steepening) UV spectral slopes with increasing redshift (e.g., Bouwens et al. 2012; Finkelstein et al. 2012; Wilkins et al. 2016; Bhatawdekar & Conselice 2021). These conclusions have been reinforced by early JWST imaging that shows very blue colors among UV-selected galaxies (e.g., Nanayakkara et al. 2023; Topping et al. 2022). We find here that including the MIRI 5.6 and 7.7 µm data provides strong evidence that the stellar populations are very blue in their rest-frame UV-to-*I*-band colors, seemingly more so than inferred from these previous studies (as in some cases the galaxies in our sample are identical to those in these other studies). Therefore, the MIRI data favor lower stellar masses and lower SFRs in distant galaxies.

Adding the MIRI data also changes the interpretation of the strength of nebular emission in the galaxies in the sample. Many of the galaxies show red colors between the HST/WFC3 and Spitzer/IRAC data (Figure 3). In the absence of MIRI data, the SED fitting interprets these colors as a combination of higher dust obscuration, older ages, and strong nebular emission. This has been reported previously for individual galaxies and for stacked samples at $z > 6$ (e.g., Castellano et al. 2017; Stark et al. 2017; De Barros et al. 2019; Hutchison et al. 2019; Endsley et al. 2021; Stefanon et al. 2022). The fact that, without the MIRI data, the models allow for higher dust obscuration and older stellar populations increases the M/L and the stellar masses. The higher dust obscuration also leads to higher SFRs. However, including the MIRI data changes this interpretation for the galaxy population. Nebular emission lines appear to be the primary explanation for the red HST/WFC3 to Spitzer/IRAC colors (while there are some galaxies where the MIRI data does not change [e.g., galaxy IDs 6811 and 7364 in Figure 3], this does not change the main result in general). This means (1) the nebular emission lines for high-redshift galaxies must be strong and (2) the stellar populations are blue. Early JWST results show emission lines remain strong in high-redshift galaxies and impact the reddest (4–5 µm) bandpasses in NIRCam (e.g., Endsley et al. 2022; Giménez-Arteaga et al. 2023; Santini et al. 2023; Topping et al. 2022; Whitler et al. 2023). The MIRI data appear necessary to extend the rest-frame wavelength coverage to >7000 Å, past the strongest of the nebular emission features in the rest-frame optical.

## 5.2. Implications for Early Star Formation and Stellar Masses in Galaxies

The difference between the delayed-$\tau$ star formation history and the model that includes bursts, while simplistic, arguably spans the gamut of available star formation histories of galaxies. Simulations show that galaxies experience many discrete bursts, but when averaged over long time baselines, the evolution is mostly smooth (e.g., Diemer et al. 2017; Iyer et al. 2019; Leja et al. 2019). Therefore, the smoothly evolving
<table><caption>**Table 7**
Ratio of the Stellar Mass Allowed in Models That Include an Early Burst of Star Formation (at $z_f = 100$) to Those That Include Only a Delayed-$\tau$ Model</caption><thead><tr><th rowspan="3">Sample</th><th colspan="4">log M<sub>*</sub>(with burst) – log M<sub>*</sub>(no burst)</th></tr><tr><th colspan="2">With MIRI data</th><th colspan="2">No MIRI data</th></tr><tr><th>Median</th><th>Scatter</th><th>Median</th><th>Scatter</th></tr></thead><tbody><tr><td>4 &lt; z &lt; 6</td><td>0.59 dex</td><td>0.16 dex</td><td>0.87 dex</td><td>0.23 dex</td></tr><tr><td>6 &lt; z &lt; 10</td><td>0.69 dex</td><td>0.04 dex</td><td>1.11 dex</td><td>0.13 dex</td></tr></tbody></table>

**Notes.** The values labeled “with burst” denote the upper limit on the stellar mass for models that include bursts. The values labeled “no burst” denote the stellar masses for models that assume only a delayed-$\tau$ star formation history. The values “with MIRI” include the 5.6 and 7.7 µm flux densities from MIRI.

delayed-$\tau$ model represents the slowly evolving evolution of galaxy star formation histories. This can reproduce the bluest colors (and lowest M/L) of the stellar population. Of course, galaxies can experience a host of stochastic bursts through changes in gas accretion or events that can cause sudden changes in the star formation, possibly as a result of mergers (e.g., Kartaltepe et al. 2023), or from the onset of strong feedback from an AGN (e.g., Wagner et al. 2016). “Bursts” of star formation from these events will add stellar mass, but if these occur at $z < 100$, then they will be younger and will have M/L lower than a model with a burst at $z_f = 100$. As such, the models with a burst at $z_f = 100$ have the maximum M/L and the oldest possible ages for a stellar population at the observed redshift of each galaxy. For these reasons, the models with this burst represent a maximum stellar mass possibly formed in these galaxies.

An advantage of focusing on galaxies at high redshifts is that the amount of *time* for discrete episodes of star formation (i.e., many individual bursts) is small, given the age of the Universe is 1 Gyr at $z = 5.7$ and 500 Myr at $z = 9.6$ for the adopted cosmology. This is shorter than the lifetimes of stars of spectral type A and later. The short age of the Universe limits the maximum M/L for stellar populations in these galaxies. Combined with longer wavelengths probed by the MIRI 5.6 and 7.7 µm data, this allows us to place tighter constraints on the mass from earlier bursts in high-redshift galaxies than has been possible previously.

Using the MIRI data, we find that the maximum stellar mass from bursts for the galaxies in our sample is lower by 0.6 dex at $4 < z < 6$ and 0.7 dex at $6 < z < 9$ (see Table 7). Furthermore, our results show that including MIRI reduces the amount of stellar mass allowed in these models—by an order of magnitude, in some cases. This is interesting in itself, as nearly all galaxies show no direct evidence for such early star formation. Comparing the median stellar masses of galaxies when fit by the delayed-$\tau$ models only and those with the delayed-$\tau$ models and the early burst at $z_f = 100$ shows they differ by a median($\Delta M_*$) ≈ 0.1 dex, and the 16th-to-84th percentile range also is similar (see the plots in Figure 7). We find no convincing cases in our sample where the galaxies require a burst at $z_f = 100$ to better fit their SEDs. This would imply that galaxies do not experience early bursts of star formation (or at least that such bursts do not form sufficient mass that we require them). We can compare our results for $4 < z < 9$ galaxies to those at lower redshifts. For example, Papovich et al. (2001) found, with K-band data for $z \sim 2–3.5$

galaxies, that early bursts could contain between three and three times the stellar mass of the younger populations that dominated the SED in those galaxies. The K band is roughly the same rest-frame wavelength for galaxies at $z \sim 2$ as MIRI is for galaxies at $z \sim 8$. With MIRI, we find that the stellar mass in an early burst can be ~4–5× the stellar mass of the younger population in $4 < z < 9$ galaxies. The early burst mass is modestly lower than the results from Papovich et al. (2001). The reason for this is a combination of effects. At $4 < z < 9$, the galaxies are blue, so they can hide more stellar mass in an early burst. However, because the Universe is younger at $4 < z < 9$, any early burst of stars is also younger, with a lower M/L and therefore a lower stellar mass. These effects offset, leading to only a modest difference in the amount of stellar mass that could be contained in bursts. Nevertheless, this highlights the importance of probing the SEDs of galaxies at rest-frame wavelengths longer than ~5000–7000 Å to constrain these populations.

Part of our findings could be impacted by biases and limitations. First, there is selection bias: all the galaxies studied here were selected in HST/WFC3 data and therefore required some rest-frame UV emission above the HST/WFC3 detection limit. It will be important to test if there is evidence for early episodes of star formation in galaxies with future JWST observations. This will be potentially very important for red galaxies, including JWST/NIRCam-selected galaxies that lack HST counterparts (Glazebrook et al. 2023; Pérez-González et al. 2023). Second, very recent work shows evidence for older stellar populations in the spatially resolved colors of $6 < z < 9$ galaxies (Giménez-Arteaga et al. 2023). As these issues become better constrained, then the ages of the stellar populations in distant galaxies could begin to inform us about when galaxies first form stars (see Whitler et al. 2023).

### 5.3. Implications for Galaxy Growth

The question of how much mass is contained in galaxies is related to the integral of the galaxies’ star formation histories. This is important because it contains the integrated record of how rapidly galaxies acquire their baryons and how efficiently they convert these into stars. This has already been discussed as an impossibly early galaxy problem, where galaxies may have acquired too much stellar mass: in typical JWST surveys, galaxies should have less stellar mass than a few times $10^{11} M_\odot$ (Behroozi & Silk 2018; Boylan-Kolchin 2023).

The effect of adding the MIRI data already is that the stellar masses and SFRs derived for galaxies tend to be overestimated when such data are excluded. The median offsets for the galaxies in our sample are 0.15 dex at $4 < z < 6$ and rise to ≈0.3 dex at $6 < z < 9$ (Figure 6 and Table 6). Assuming that these (median) offsets apply to previous estimates of galaxy stellar masses, then it implies that measurements of the cosmic SFR density (SFRD, which is the average SFR in all galaxies per comoving volume element) are similarly biased to higher values at these redshifts.

Figure 10 shows the impact of these lower SFRs on the cosmic stellar-mass density, $\rho_*$: First, the figure shows $\rho_*$ derived from the integral of the SFRD for two empirical models calibrated against measurements from the literature (Madau & Dickinson 2014; Finkelstein 2016). Finkelstein (2016) shows that these models are consistent with a compilation of measurements of $\rho_*$ from the literature at $4 < z < 10$ prior to
<figure><img src="image_21.png" /><figcaption>**Figure 10.** Evolution of the cosmic stellar mass density, $\rho_*$, in galaxies from $4 < z < 10$. The lines show pre-JWST constraints from Finkelstein (2016) and Madau & Dickinson (2014). The data points show recent measurements of $\rho_*$ at $z > 6$ from the literature (Kikuchihara et al. 2020 [Ki20], Bhatawdekar & Conselice 2021 [Bh21], Stefanon et al. 2021 [St21], and Santini et al. 2022 [Sa22]), which largely follow the pre-JWST constraints. The shaded regions show maximally allowable stellar mass density assuming galaxies experience a burst at $z = 100$ followed by “normal” star formation. Constraints lacking JWST/MIRI coverage to rest-frame 1 µm allow for a stellar mass density that is up to 0.8 dex higher at $z = 4$ and 1.4 dex higher at $z = 10$. Including MIRI 5.6 and 7.7 µm data lowers the maximum allowed by up to a factor of 5.</figcaption></figure>

the launch of JWST, (Duncan et al. 2014; Oesch et al. 2014; Grazian et al. 2015; Song et al. 2016). The thick line in Figure 10 labeled “MIRI corrected” shows the empirical model of Finkelstein corrected by the offsets of the SFRs derived in Table 6. To derive these corrections, we have interpolated the results from Table 6 assuming median redshifts of $z = 5$ and 8 for the derived offsets in the two redshift bins. Although we use the offsets derived from the SFRs, using those for the stellar masses changes the results by ≈0.1 dex. We note, however, that because the MIRI data imply offsets in the SFRs of galaxies, they similarly lower the values of the cosmic SFRD at $z > 4$.

Second, the MIRI data improve the constraints on the amount of stellar mass possible in early bursts of star formation (Section 4.4). This is illustrated by the shaded regions in Figure 10. To derive the area in the shaded swaths, we applied the ratio between the mass permitted in early bursts at $z_f = 100$ to the fiducial value (listed in Table 7) and interpolated them over redshift as above. Therefore, the effects of adding the MIRI data both lower SFR (and stellar masses) and limit the total stellar mass allowed in early bursts. The combination of these effects reduces the upper bound on the total cosmic stellar mass density allowed by the data by 0.4 dex at $z = 4$ and by 1.0 dex at $z = 9$. As illustrated in Figure 10, this implies that the JWST/MIRI data have constrained the stellar mass in galaxies at $z = 9$ to be less than 0.1% that of the present-day value ($\rho_*(z=0)$).

## 5.4. Implications for the Nature of Early Galaxies

The fact that galaxies have bluer colors than inferred from previous constraints has other consequences (i.e., they are “bluer than they appear”). Our findings dovetail with other recent results from JWST NIRCam imaging. The results in this paper are likely only the first foray into the properties of distant galaxies using the longer-wavelength data available from MIRI. Future studies will be able to combine both NIRCam and MIRI imaging, providing JWST-quality data from 0.8 to 10 µm, with NIRSpec spectroscopy, which will improve the constraints presented here. Already there are indications, using JWST/NIRCam data only, that the number density of luminous galaxies at $z > 10$ may be much higher than predictions (e.g., Finkelstein et al. 2022b; Bouwens et al. 2023; Donnan et al. 2023; Harikane et al. 2023; Naidu et al. 2022; Robertson et al. 2023). One interpretation of these discoveries is that the UV luminosity per unit of stellar mass (the UV “efficiency”) may be higher than models predict. This could be a result of changes in the stellar populations such that they produce more ionizing photons per unit mass. This is expected both for the case that the galaxies contain high-mass, very metal-poor stars (Olivier et al. 2022) or if their stellar populations have an IMF weighted toward high-mass stars (Raiter et al. 2010). Importantly, either of these effects would further decrease the M/L ratio of the stellar population, and thus lead to even lower stellar masses than we have measured with the MIRI data. Therefore, it will be important to constrain the nature of the stars in these galaxies. In the meantime, the MIRI data have better constrained the available light in stars at these early epochs, revealing that they contain at least three times less “light” at rest-frame 0.7–1 µm than previously known.

## 6. Conclusions

In this paper, we have presented results from CEERS on the stellar population parameters for 28 galaxies with redshifts $4 < z < 9$ using new imaging data from JWST/MIRI at 5.6 and 7.7 µm. Our galaxy sample was detected in deep data from HST/WFC3 and ACS and has observations from Spitzer/IRAC at 3.6 and 4.5 µm. The MIRI 5.6 and 7.7 µm data extend the coverage of the rest-frame spectral energy distribution to nearly 1 µm for galaxies in this redshift range. We use these data to study the improvements in the stellar masses and SFRs of the galaxies at these redshifts when the MIRI data are included. Our main results are the following.

1. Galaxies at $4 < z < 9$ have bluer rest-frame UV-$I$-band colors ($m_{1600\AA} - I$) when using the MIRI data compared to when they are excluded. Using the MIRI data, we model the SEDs using stellar population synthesis models (with BAGPIPES). When we compare the average galaxy SED (Figure 9), we find that models that include the MIRI data are (on average) $\Delta(m_{1600\AA} - I) \approx 0.4$ mag bluer in their rest-frame colors compared to models that exclude the MIRI data.

2. Galaxies generally have lower stellar masses and SFRs when the MIRI data are included. For the majority of the galaxies (Figure 6), adding the MIRI data reduces the derived stellar masses by 0.25 dex at $4 < z < 6$ and by 0.38 dex at $6 < z < 9$. Similarly, including the MIRI data reduces the SFRs by 0.15 dex at $4 < z < 6$. Because the impact is larger on the stellar masses than on the SFRs,
the specific SFRs will be increased by approximately
$\Delta\log\text{SFR} - \Delta\log M_* \approx 0.1$ dex when the MIRI data are
included. There are multiple reasons the stellar masses
and SFRs are lowered when we include the MIRI data.
The first reason is that the galaxies are blue, and the fits
favor models with lower dust attenuation and models
with lower M/L in general. The second reason is that, in
many cases, the IRAC 3.6 and 4.5 µm data probe the rest-
frame optical, and these show indications of containing
light from strong emission lines (e.g., redshifted $H\alpha$ +
$[O\,III]$, $H\alpha+[N\,II]$, $[O\,II]$, etc.). These boost the flux in
these bands. In the absence of MIRI data, the models
cannot determine if the red rest-frame UV-optical colors
are a result of dust attenuation, older stellar populations,
or strong emission lines (or all of them). The parameter
constraints then give more weight (probability density) to
models with higher stellar masses and SFRs. When the
MIRI data are included, they probe more of the stellar
continuum at $>7000$ Å rest frame. The model fits that
include the MIRI data then show the dominant effect in
the rest-frame optical to be strong nebular emission lines.
This problem will persist for models that use data from
NIRCam, as it also is limited to wavelengths less than
$5$ µm, but this can be tested with forthcoming data sets.

3. The amount of stellar mass that could have formed in early bursts is lower. We estimated the amount of stellar mass formed by using a star formation history that includes an early burst (at $z_f=100$) in addition to a smoothly evolving component. A stellar population formed at this early time would fade and redden with time, and it would have the highest M/L at any subsequent time. Therefore, it represents an upper limit on the amount of mass that could exist in these galaxies. The MIRI data improve the constraint on this stellar population by probing longer wavelengths where the impact from this stellar population is more pronounced. Figure 8 shows that, without the MIRI data, the amount of stellar mass in this population can be as much as 0.9 dex higher at $4 < z < 6$ and as much as 1.1 dex higher at $6 < z < 9$. Including the MIRI data, these drop to 0.6 dex at $4 < z < 6$ and 0.7 dex at $6 < z < 9$. Therefore, adding the MIRI reduces the amount of mass in early bursts by a factor of order 2 compared to when no MIRI data are used.

4. Our analysis of the MIRI 5.6 and 7.7 µm therefore provides evidence that there is less star formation in distant galaxies than found in previous studies, because the SFRs and stellar masses are lowered. The MIRI data also reduce the limits on the amount of stellar mass possibly formed at early times. The combination of these results has implications for the evolution of the cosmic stellar-mass density, $\rho_*$*. We showed (Figure 10) that applying our results lowers the amount of stellar mass density in galaxies at $z=9$ to be less than 0.1% of the present-day, $z=0$, value. This is an order of magnitude lower than implied by previous studies (i.e., pre-JWST).

We wish to thank everyone that brought JWST to fruition.
We also thank our other colleagues in the CEERS collaboration
for their hard work and valuable contributions on this project.
The authors offer their sincerest thanks to the anonymous
referee whose comments and suggestions improved the quality

and rarity of this work. C.P. thanks Marsha and Ralph Schilling
for generous support of this research. Portions of this research
were conducted with the advanced computing resources
provided by Texas A&M High Performance Research
Computing (HPRC, http://hprc.tamu.edu). This work bene-
fited from support from the George P. and Cynthia Woods
Mitchell Institute for Fundamental Physics and Astronomy at
Texas A&M University. This work acknowledges support from
the NASA/ESA/CSA James Webb Space Telescope through
the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy,
Incorporated, under NASA contract NAS5-03127. Support
for program No. JWST-ERS01345 was provided through a
grant from the STScI under NASA contract NAS5-03127.

*Software:* AstroPy (Astropy Collaboration et al. 2013), BAGPIPES(Carnall et al. 2018), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), photutils (Bradley et al. 2020), PyPHER (Boucaud et al. 2016a), SE(Bertin & Arnouts 1996), SciPy (Virtanen et al. 2020), Seaborn (Waskom 2021).

Appendix
Impact of Crowded Sources in IRAC Data

One source of potential bias relates to the photometry of our
sources in the Spitzer/IRAC data. As illustrated in the images
(Figure 2), some objects have bright neighbors. In the case of
the IRAC images, the light from the wings of these objects can
blend with that for our sources. There is a large body of
literature on the subject of performing crowded source
photometry (Laidler et al. 2007; Labbé et al. 2013; Merlin
et al. 2015, 2016). We have used the catalog from Finkelstein
et al. (2022a) who used the HST/F160W image as a prior for
the locations of sources and neighbors. Source photometry is
then carried out using TPHOT (Merlin et al. 2015), which
estimates the source flux from objects simultaneously when
measuring photometry. While this method is theoretically
robust, residuals from poorly modeling ePSFs and changes in
galaxy morphology with wavelength (the “morphological K-
correction”; see Papovich et al. 2005) can lead to systematic
uncertainties in source photometry.

To test if our results are impacted by blended sources in the
IRAC bands, we did the following. We first searched around
each of the galaxies in our sample and identified galaxies that
have neighbors in the MIRI 5.6 µm catalog within a radius of
$r \le 3''$ and a magnitude of $[5.6] \le 26.7$ mag (near the flux
limit). We selected neighbors in the MIRI 5.6 µm image as the
central wavelength is closest to that of IRAC for our data set
(see Figure 1). The IRAC ePSF has an FWHM of $\approx 2''$, so any
source within $3''$ in the MIRI data could therefore have IRAC
light blended with our source of interest.

From our sample, we identified 11 galaxies that have a
neighbor within 3" in the MIRI 5.6 µm image. To estimate their
effect on our study, we removed these objects from the sample
and recomputed the offsets in stellar mass and SFR for the
results that include the MIRI 5.6 and 7.7 µm data versus the
results that exclude the MIRI data. These results are shown in
Figure 11. Contrasting this figure with the one for the full
sample (Figure 6) shows there is little change in the median
offsets in stellar mass and SFR. The galaxy sample used in this
appendix is obviously smaller, but the median values do not
change appreciably. For the sample that excludes blended
objects, the offsets in stellar mass are $\Delta \log M_* = 0.21$ dex for
<figure><img src="image_22.png" /><figcaption>**Figure 11.** Testing the impact of sources with “crowded” IRAC photometry. The plots in this figure are similar to those in Figure 6, and they compare the stellar masses and SFRs derived from the SED modeling for galaxies including the MIRI F560W and F770W data and without the MIRI data. In both panels, the results show the difference between the mass (SFR) derived without MIRI data and the mass (SFR) derived including the MIRI data. In this figure, we have excluded objects that have a neighbor with $r \le 3''$ and MIRI [5.6] $\le 26.7$ mag. This eliminates 11 objects and allows us to test if crowding in the IRAC data (which has lower angular resolution) impacts object photometry in the IRAC bands. We do not observe any significant offset compared to the results in Figure 6: the median offsets in stellar mass and SFR change by $\approx 0.1$ dex. Therefore, we conclude that blended IRAC photometry does not significantly impact the results here.</figcaption></figure>

$4 < z < 6$ and $0.53$ for $6 < z < 9$ (though the latter now includes only three galaxies). The offsets in SFR are $\Delta\log\text{SFR} = 0.13$ dex for $4 < z < 6$ and $0.42$ dex for $6 < z < 9$. These are within $\approx 0.1$ dex of the values reported for the full sample (in Figure 6).

Similarly, we also investigated how the IRAC data for sources with close neighbors impact our finding that the stellar populations of the galaxies in our sample are generally “bluer” when the MIRI data are included in the analysis (see Section 5.1 and Figure 9). We repeated our analysis of the rest-frame colors in Section 5.1 with our sample of galaxies that excludes those objects with a neighboring MIRI 5.6 µm source with $[5.6] < 26.7$ mag and within $r \le 3''$. We find that, in this case, the relative rest-frame colors change only slightly. The rest-frame far-UV-$I$ color become bluer by $0.015$ mag (to have a total rest-frame (blue) color of $\Delta(m_{1600\AA} - I) \approx 0.42$ mag) when the objects with crowded IRAC photometry are excluded.

Therefore, we conclude that our results are not dominated by photometry from sources crowded in the IRAC data. Obviously, future studies using JWST/NIRCam will be valuable to testing the IRAC photometry (see, e.g., Bagley et al. 2023).

## ORCID iDs

Casey Papovich https://orcid.org/0000-0001-7503-8482

Justin W. Cole https://orcid.org/0000-0002-6348-1900

Guang Yang https://orcid.org/0000-0001-8835-7722

Steven L. Finkelstein https://orcid.org/0000-0001-8519-1130

Guillermo Barro https://orcid.org/0000-0001-6813-875X

Véronique Buat https://orcid.org/0000-0003-3441-903X

Denis Burgarella https://orcid.org/0000-0002-4193-2539

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

Paola Santini https://orcid.org/0000-0002-9334-8705

Lise-Marie Seillé https://orcid.org/0000-0001-7755-4755

Lu Shen https://orcid.org/0000-0001-9495-7759

Pablo Arrabal Haro https://orcid.org/0000-0002-7959-8783

Micaela B. Bagley https://orcid.org/0000-0002-9921-9218

Eric F. Bell https://orcid.org/0000-0002-5564-9873

Laura Bisigello https://orcid.org/0000-0003-0492-4924

Antonello Calabrò https://orcid.org/0000-0003-2536-1614

Caitlin M. Casey https://orcid.org/0000-0002-0930-6466

Marco Castellano https://orcid.org/0000-0001-9875-8263

Katherine Chworowsky https://orcid.org/0000-0003-4922-0613

Nikko J. Cleri https://orcid.org/0000-0001-7151-009X

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

M. C. Cooper https://orcid.org/0000-0003-1371-6019

Mark Dickinson https://orcid.org/0000-0001-5414-5131

Henry C. Ferguson https://orcid.org/0000-0001-7113-2738

Adriano Fontana https://orcid.org/0000-0003-3820-2823

Mauro Giavalisco https://orcid.org/0000-0002-7831-8751

Andrea Grazian https://orcid.org/0000-0002-5688-0663

Norman A. Grogin https://orcid.org/0000-0001-9440-8872

Nimish P. Hathi https://orcid.org/0000-0001-6145-5090

Benne W. Holwerda https://orcid.org/0000-0002-4884-6756

Taylor A. Hutchison https://orcid.org/0000-0001-6251-4988

Jeyhan S. Kartaltepe https://orcid.org/0000-0001-9187-3605

Lisa J. Kewley https://orcid.org/0000-0001-8152-3943

Allison Kirkpatrick https://orcid.org/0000-0002-5537-8110

Dale D. Kocevski https://orcid.org/0000-0002-8360-3880

Anton M. Koekemoer https://orcid.org/0000-0002-6610-2048

Rebecca L. Larson https://orcid.org/0000-0003-2366-8858

Arianna S. Long https://orcid.org/0000-0002-7530-8857

Ray A. Lucas https://orcid.org/0000-0003-1581-7825

Laura Pentericci https://orcid.org/0000-0001-8940-6768

Nor Pirzkal https://orcid.org/0000-0003-3382-5941

Swara Ravindranath https://orcid.org/0000-0002-5269-6527
Rachel S. Somerville https://orcid.org/0000-0002-6748-6821

Jonathan R. Trump https://orcid.org/0000-0002-1410-0470

Stephanie M. Urbano Stawinski https://orcid.org/0000-0001-8169-7249

Benjamin J. Weiner https://orcid.org/0000-0001-6065-7483

Stephen M. Wilkins https://orcid.org/0000-0003-3903-6935

L. Y. Aaron Yung https://orcid.org/0000-0003-3466-035X

Jorge A. Zavala https://orcid.org/0000-0002-7051-1100

## References

Antwi-Danso, J., Papovich, C., Leja, J., et al. 2023, ApJ, 943, 166

Arellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, ApJL, 940, L23

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Backhaus, B. E., Trump, J. R., Cleri, N. J., et al. 2022, ApJ, 926, 161

Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12

Bailer-Jones, C. A. L. 2017, Practical Bayesian Inference (Cambridge: Cambridge Univ. Press)

Barkana, R., & Loeb, A. 2001, PhR, 349, 125

Behroozi, P., & Silk, J. 2018, MNRAS, 477, 5382

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bhatawdekar, R., & Conselice, C. J. 2021, ApJ, 909, 144

Bisigello, L., Caputi, K. I., Colina, L., et al. 2017, ApJS, 231, 3

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Boucaud, A., Bocchio, M., Abergel, A., et al. 2016a, PyPHER: Python-based PSF Homogenization kERnels, Astrophysics Source Code Library, ascl:1609.022

Boucaud, A., Bocchio, M., Abergel, A., et al. 2016b, A&A, 596, A63

Bouwens, R., Illingworth, G., Oesch, P., et al. 2023, MNRAS, Advance Access

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJ, 754, 83

Bouwens, R. J., Stefanon, M., Oesch, P. A., et al. 2019, ApJ, 880, 25

Boyett, K., Mascia, S., Pentericci, L., et al. 2022, ApJL, 940, L52

Boylan-Kolchin, M. 2023, NatAs, Advanced Online Publication

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: 1.0.0, Zenodo, doi:10.5281/zenodo.4044744

Brinchmann, J. 2022, arXiv:2208.07467

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

Buat, V., Ciesla, L., Boquien, M., Malek, K., & Burgarella, D. 2019, A&A, 632, A79

Burgarella, D., Theulé, P., Buat, V., et al. 2023, A&A, 671, A123

Calzetti, D. 2001, PASP, 113, 1449

Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, ApJ, 873, 44

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Castellano, M., Pentericci, L., Fontana, A., et al. 2017, ApJ, 839, 73

Chabrier, G. 2003, PASP, 115, 763

Chworowsky, K., Finkelstein, S. L., Spilker, J., et al. 2023, arXiv:2305.06309

Conroy, C. 2013, ARA&A, 51, 393

Curti, M., D'Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, NatAs, Advanced Online Publication

De Barros, S., Oesch, P. A., Labbé, I., et al. 2019, MNRAS, 489, 2355

Dickinson, M., Papovich, C., Ferguson, H. C., & Budavári, T. 2003, ApJ, 587, 25

Diemer, B., Sparre, M., Abramson, L. E., & Torrey, P. 2017, ApJ, 839, 26

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

Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960

Endsley, R., Stark, D. P., Chevillard, J., & Charlot, S. 2021, MNRAS, 500, 5229

Endsley, R., Stark, D. P., Whitler, L., et al. 2022, arXiv:2208.14999

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Finkelstein, S. L. 2016, PASA, 33, e037

Finkelstein, S. L., Bagley, M., Song, M., et al. 2022a, ApJ, 928, 52

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

Finkelstein, S. L., Papovich, C., Giavalisco, M., et al. 2010, ApJ, 719, 1250

Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, ApJ, 756, 164

Finlator, K., Oppenheimer, B. D., & Davé, R. 2011, MNRAS, 410, 1703

Fujimoto, S., Ouchi, M., Nakajima, K., et al. 2022, arXiv:2212.06863

Giménez-Arteaga, C., Oesch, P. A., Brammer, G. B., et al. 2023, ApJ, 948, 126

Glazebrook, K., Nanayakkara, T., Jacobs, C., et al. 2023, ApJL, 947, L25

Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96

García-Arégumáñez, Á., Pérez-González, P. G., Gil de Paz, A., et al. 2023, ApJ, 944, 3

Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5

Heintz, K. E., Giménez-Arteaga, C., Fujimoto, S., et al. 2023, ApJL, 944, L30

Hunter, J. D. 2007, CSE, 9, 90

Hutchison, T. A., Papovich, C., Finkelstein, S. L., et al. 2019, ApJ, 879, 70

Iyer, K. G., Gawiser, E., Faber, S. M., et al. 2019, ApJ, 879, 116

Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJL, 946, L15

Katz, H., Saxena, A., Cameron, A. J., et al. 2023, MNRAS, 518, 592

Kikuchihara, S., Ouchi, M., Ono, Y., et al. 2020, ApJ, 893, 60

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

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Labbé, I., Oesch, P. A., Bouwens, R. J., et al. 2013, ApJL, 777, L19

Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Natur, 616, 266

Laidler, V. G., Papovich, C., Grogin, N. A., et al. 2007, PASP, 119, 1325

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

Laporte, N., Nakajima, K., Ellis, R. S., et al. 2017, ApJ, 851, 40

Larson, R. B., & Tinsley, B. M. 1978, ApJ, 219, 46

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Maraston, C., Pforr, J., Renzini, A., et al. 2010, MNRAS, 407, 830

Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2022, arXiv:2211.08255

Merlin, E., Bourne, N., Castellano, M., et al. 2016, A&A, 595, A97

Merlin, E., Fontana, A., Ferguson, H. C., et al. 2015, A&A, 582, A15

Miralda-Escudé, J. 2003, ApJ, 597, 66

Mobasher, B., Dahlen, T., Ferguson, H. C., et al. 2015, ApJ, 808, 101

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

Nanayakkara, T., Glazebrook, K., Jacobs, C., et al. 2023, ApJL, 947, L26

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2014, ApJ, 786, 108

Oesch, P. A., van Dokkum, P. G., Illingworth, G. D., et al. 2015, ApJL, 804, L30

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

Olivier, G. M., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 938, 16

Papovich, C., Dickinson, M., & Ferguson, H. C. 2001, ApJ, 559, 620

Papovich, C., Dickinson, M., Giavalisco, M., Conselice, C. J., & Ferguson, H. C. 2005, ApJ, 631, 101

Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, J. M., & Giavalisco, M. 2011, MNRAS, 412, 1123

Papovich, C., Labbé, I., Quadri, R., et al. 2015, ApJ, 803, 26

Papovich, C., Moustakas, L. A., Dickinson, M., et al. 2006, ApJ, 640, 92

Papovich, C., Simons, R. C., Estrada-Carpenter, V., et al. 2022, ApJ, 937, 22

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

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

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Pforr, J., Maraston, C., & Tonini, C. 2012, MNRAS, 422, 3285

Raiter, A., Schaerer, D., & Fosbury, R. A. E. 2010, A&A, 523, A64

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

Roberts-Borsani, G. W., Bouwens, R. J., Oesch, P. A., et al. 2016, ApJ, 823, 143

Roberts-Borsani, G. W., Ellis, R. S., & Laporte, N. 2020, MNRAS, 497, 3440

Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, NatAs, Advanced Online Publication

Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 816, 23

Sanders, R. L., Shapley, A. E., Reddy, N. A., et al. 2020, MNRAS, 491, 1427

Santini, P., Castellano, M., Fontana, A., et al. 2022, ApJ, 940, 135

Santini, P., Fontana, A., Castellano, M., et al. 2023, ApJL, 942, L27

Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Smit, R., Bouwens, R. J., Franx, M., et al. 2015, ApJ, 801, 122

Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Stark, D. P., Ellis, R. S., Charlot, S., et al. 2017, MNRAS, 464, 469

Stark, D. P., Richard, J., Charlot, S., et al. 2015, MNRAS, 450, 1846

Stefanon, M., Bouwens, R. J., Illingworth, G. D., et al. 2022, ApJ, 935, 94

Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29

Stefanon, M., Yan, H., Mobasher, B., et al. 2017, ApJS, 229, 32

Steinhardt, C. L., Kokorev, V., Rusakov, V., Garcia, E., & Sneppen, A. 2022, arXiv:2208.07879

Sun, F., Egami, E., Pirzkal, N., et al. 2022, arXiv:2209.03374

Tang, M., Stark, D. P., Chevillard, J., & Charlot, S. 2019, MNRAS, 489, 2572
Tinsley, B. M. 1980, *FCPh*, **5**, 287

Topping, M. W., Stark, D. P., Endsley, R., et al. 2022, *ApJ*, **941**, 153

Tran, K.-V. H., Forrest, B., Alcorn, L. Y., et al. 2020, *ApJ*, **898**, 45

Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2023, *ApJ*, **945**, 35
van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, *CSE*, **13**, 22

van der Wel, A., Straughn, A. N., Rix, H. W., et al. 2011, *ApJ*, **742**, 111

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, *NatMe*, **17**, 261

Visbal, E., Bryan, G. L., & Haiman, Z. 2020, *ApJ*, **897**, 95

Wagner, A. Y., Bicknell, G. V., Umemura, M., Sutherland, R. S., & Silk, J. 2016, *AN*, **337**, 167

Waskom, M. L. 2021, *JOSS*, **6**, 3021

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

Wilkins, S. M., Bouwens, R. J., Oesch, P. A., et al. 2016, *MNRAS*, **455**, 659

Wise, J. H., Turk, M. J., Norman, M. L., & Abel, T. 2012, *ApJ*, **745**, 50

Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, *ApJ*, **592**, 645

Zitrin, A., Labbé, I., Belli, S., et al. 2015, *ApJL*, **810**, L12