# MIDIS: Strong (Hβ+[O III]) and Hα Emitters at Redshift z 7-8 Unveiled with JWST NIRCam and MIRI Imaging in the Hubble eXtreme Deep Field

P. Rinaldi, K. I. Caputi, L. Costantin, S. Gillman, E. Iani, P. G.

Pérez-González, G. Östlin, L. Colina, T. R. Greve, H. U. Noorgard-Nielsen, et al.

► To cite this version:

P. Rinaldi, K. I. Caputi, L. Costantin, S. Gillman, E. Iani, et al., MIDIS: Strong (Hβ+[O III]) and Hα Emitters at Redshift z 7-8 Unveiled with JWST NIRCam and MIRI Imaging in the Hubble eXtreme Deep Field. The Astrophysical Journal, 2023, 952, 10.3847/1538-4357/acdc27. insu-04479062

HAL Id: insu-04479062

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

Submitted on 6 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# MIDIS: Strong (Hβ+[O III]) and Hα Emitters at Redshift $z \simeq 7-8$ Unveiled with JWST NIRCam and MIRI Imaging in the Hubble eXtreme Deep Field

P. Rinaldi¹, K. I. Caputi¹,², L. Costantin³, S. Gillman²,⁴, E. Iani¹, P. G. Pérez-González³, G. Östlin⁵, L. Colina²,³, T. R. Greve²,⁴, H. U. Noorgard-Nielsen²,⁴, G. S. Wright⁶, A. Alonso-Herrero⁷, J. Álvarez-Márquez³, A. Eckart⁸, M. García-Marín⁹, J. Hjorth¹⁰, O. Ilbert¹¹, S. Kendrew⁹, A. Labiano¹², O. Le Fèvre¹¹, J. Pye¹³, T. Tikkanen¹³, F. Walter¹⁴, P. van der Werf¹⁵, M. Ward¹⁶, M. Annunziatella³,¹⁷, R. Azzollini³,¹⁸, A. Bik⁵, L. Boogaard¹⁴, S. E. I. Bosman¹⁴, A. Crespo Gómez³, I. Jermann²,⁴, D. Langeroodi¹⁰, J. Melinder⁵, R. A. Meyer¹⁴, T. Moutard¹¹, F. Peissker⁸, M. Topinka¹⁸, E. van Dishoeck¹⁵, M. Güdel¹⁹,²⁰,²¹, Th. Henning¹⁴, P.-O. Lagage²², T. Ray¹⁸, B. Vandenbussche²³, C. Waelkens²³, R. Navarro-Carrera¹, and V. Kokorev¹

¹Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV Groningen, The Netherlands; rinaldi@astro.rug.nl

²Cosmic Dawn Center (DAWN), Denmark

³Centro de Astrobiología (CAB), CSIC-INTA, Carretera de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain

⁴DTU-Space, Elektrovej, Building 328, DK-2800, Kongens Lyngby, Denmark

⁵Department of Astronomy, Stockholm University, Oscar Klein Centre, AlbaNova University Centre, SE-106 91 Stockholm, Sweden

⁶UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, UK

⁷Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain

⁸I.Physikalisches Institut der Universität zu Köln, Zülpicher Strasse 77, D-50937 Köln, Germany

⁹European Space Agency/Space Telescope Science Institute, 3700 San Martin Drive, Baltimore MD 21218, USA

¹⁰DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark

¹¹Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388, Marseille, France

¹²Telespazio UK for the European Space Agency (ESA), ESAC, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Spain

¹³School of Physics & Astronomy, Space Research Centre, Space Park Leicester, University of Leicester, 92 Corporation Road, Leicester, LE4 5SP, UK

¹⁴Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

¹⁵Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands

¹⁶Centre for Extragalactic Astronomy, Durham University, South Road, Durham, DH1 3LE, UK

¹⁷INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, I-80131 Napoli, Italy

¹⁸Dublin Institute for Advanced Studies, Astronomy & Astrophysics Section, 31 Fitzwilliam Place, Dublin 2, Ireland

¹⁹Department of Astrophysics, University of Vienna, Türkenschanzstr 17, A-1180 Vienna, Austria

²⁰Max-Planck-Institut für Astronomie (MPIA), Königstuhl 17, D-69117 Heidelberg, Germany

²¹ETH Zürich, Institute for Particle Physics and Astrophysics, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland

²²AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

²³Institute of Astronomy, KU Leuven, Celestijnenlaan 200D bus 2401, B-3001 Leuven, Belgium

Received 2023 January 30; revised 2023 May 23; accepted 2023 June 1; published 2023 July 25

## Abstract

We make use of JWST medium-band and broadband NIRCam imaging, along with ultradeep MIRI 5.6 µm imaging, in the Hubble eXtreme Deep Field to identify prominent line emitters at $z \simeq 7-8$. Out of a total of 58 galaxies at $z \simeq 7-8$, we find 18 robust candidates ($\simeq 31\%$) for (Hβ + [O III]) emitters, based on their enhanced fluxes in the F430M and F444W filters, with EW$_{0}$(Hβ + [O III]) $\simeq 87-2100$ Å. Among these emitters, 16 lie in the MIRI coverage area and 12 exhibit a clear flux excess at 5.6 µm, indicating the simultaneous presence of a prominent Hα emission line with EW$_{0}$(Hα) $\simeq 200-3000$ Å. This is the first time that Hα emission can be detected in individual galaxies at $z > 7$. The Hα line, when present, allows us to separate the contributions of Hβ and [O III] to the (Hβ + [O III]) complex and derive Hα-based star formation rates (SFRs). We find that in most cases [O III]/Hβ > 1. Instead, two galaxies have [O III]/Hβ < 1, indicating that the NIRCam flux excess is mainly driven by Hβ. Most prominent line emitters are very young starbursts or galaxies on their way to/from the starburst cloud. They make for a cosmic SFR density $\log_{10}(\rho_{SFR_{H\alpha}}/(M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3})) \simeq -2.35$, which is about a quarter of the total value ($\log_{10}(\rho_{SFR_{tot}}/(M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3})) \simeq -1.76$) at $z \simeq 7-8$. Therefore, the strong Hα emitters likely had a significant role in reionization.

*Unified Astronomy Thesaurus concepts: Galaxy formation (595); Reionization (1383); Star formation (1569); Galaxy evolution (594); Starburst galaxies (1570)*

## 1. Introduction

Quantifying the presence and properties of galaxies present at the Epoch of Reionization (EoR) is necessary to explain how this major phase transition of the universe has occurred. Over

the past decade, many studies have focused on this topic, but a few important problems complicated the selection of galaxies at this cosmic time. The increasing intergalactic medium absorption with redshift means that basically all photons blueward of the Lyα spectral line at $\lambda_{rest} = 1216$ Å cannot reach us. Indeed, it is well known that the incidence of Lyα emitters (LAEs) has a sharp drop at $z > 7$ (e.g., Fontana et al. 2010; Ono et al. 2012; Caruana et al. 2014; Pentericci et al. 2014). Therefore, other emission lines at longer wavelengths must be considered to

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** RGB mosaic of the Hubble eXtreme Ultra Deep Field. This image has been obtained by exploiting the HST and JWST images currently available in this field. In particular, the background image has been obtained by combining the HST and JWST/NIRCam filters. The zoom-in shows the region covered by MIRI/F560W. In this case, to create the RGB image, we adopted the JWST filters only where R: F560W; G: F430M, F444W, F460M, and F480M; B: F182M and F210M. The RGB images have been produced using the software Trilogy (Coe 2015).</figcaption></figure>

facilitate the search of galaxies at such high redshifts (e.g., Stark et al. 2015).

However, detecting the optical emission from atomic transitions at $z > 7$ was virtually impossible until now, given the lack of sufficiently sensitive near and mid-infrared observatories. The recent advent of the JWST is now radically changing this situation by offering, for the first time, sensitive imaging and spectroscopy at such long wavelengths. Indeed, in the first six months of operations, the JWST has enabled a number of studies of $z > 7$ galaxies, particularly on their line emission properties (e.g., Arellano-Córdova et al. 2022; Langeroodi et al. 2022; Morishita & Stiavelli 2023; Trump et al. 2023; Wang et al. 2022; Williams et al. 2023).

With imaging, the search of line emitters is facilitated by the fact that the rest-frame equivalent widths (EWs$_0$) of some of the main optical emission lines appear to increase, on average, with the redshift (e.g., De Barros et al. 2019; Matthee et al. 2023). This has allowed for the search of prominent line emitters at intermediate and high redshifts, by identifying galaxies with photometric excess in narrowband images (e.g., Khostovan et al. 2016) and even broadband images (e.g., Faisst et al. 2016; Roberts-Borsani et al. 2016; Smit et al. 2016; Caputi et al. 2017). This trend of increasing EWs$_0$ with the redshift is indicative of an evolution in the galaxy average specific star formation rates (sSFR; e.g., Faisst et al. 2016; Tang et al. 2019), as well as the conditions of their interstellar medium (ISM; e.g., Schaerer & de Barros 2009).

At $z > 7$ both the $H\beta\lambda4861\AA$ and [O III] $\lambda\lambda4959$, 5007 emission lines are shifted into the JWST's NIRCam (Rieke et al. 2005) wavelength range, making that these lines together can produce a flux excess in the NIRCam filters at $\simeq 4-5 \mu$m. In turn, the ($H\alpha\lambda6563 + [NII]\lambda\lambda6548, 6583 + [SII]\lambda\lambda6716, 6730$) complex appears in the MIRI (Rieke et al. 2015; Wright et al. 2015) wavelength domain at observed $> 5 \mu$m.

In this paper, we make use of publicly available NIRCam images in the Hubble eXtreme Deep Field (XDF) to search for ($H\beta + [OIII]$) emitters at $z \simeq 7-8$. In most of this field, we also benefit from ultradeep MIRI 5.6 $\mu$m imaging, which we analyze to search for the presence of $H\alpha$ emission in the same galaxies. This is the first time that the $H\alpha$ line can be detected and quantified in individual galaxies at $z > 7$. This paper is organized as follows: in Section 2 we describe the data sets, photometric measurements and spectral energy distribution (SED) fitting that allows us to select galaxies at $z \simeq 7-8$. In Section 3 we explain our methodology to identify strong ($H\beta + [OIII]$) and $H\alpha$ emitters among these galaxies. We present all our results in Section 4 and our conclusions in Section 5. Throughout this paper, we consider a cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_{\Lambda} = 0.7$. All magnitudes are total and refer to the AB system (Oke & Gunn 1983). A Chabrier (2003) initial mass function (IMF) is assumed.

## 2. Data Sets, Photometry, and SED Fitting

### 2.1. Data Sets

The Hubble XDF (Illingworth et al. 2013; see Figure 1) is a small field of the sky with the deepest Hubble Space Telescope (HST) observations ever taken since this telescope started operations more than 30 years ago. This field has been the main window to study the early universe before the JWST advent, with numerous works scientifically exploiting its unique possibilities. Now in the JWST era, the HST data in the XDF and surroundings are being enhanced with deep imaging and spectroscopy obtained with the JWST/NIRCam and MIRI, extending the wavelength coverage of high-spatial-resolution observations to the mid-infrared.
### 2.1.1. JWST/NIRCam

In this work, we made use of the recent JWST/NIRCam images collected by Williams et al. (2023) in a General Observers Cycle-1 program across the Hubble eXtreme Ultra Deep Field (HUDF; PID: 1963; PI: Christina C. Williams). Observations have been taken in five JWST/NIRCam medium bands: F182M, F210M, F430M, F460M, and F480M. In particular, 7.8 hr of the total integration time have been dedicated to F182M, F210M, and F480M. Instead, only 3.8 hr of observations have been collected for F430M and F460M. In order to complement these data sets, we also made use of the imaging data taken as part of The First Reionization Epoch Spectroscopic Complete Survey (FRESCO; Oesch et al. 2021, 2023, PID: 1895; PI: Pascal Oesch). On the one hand, this GO program allowed us to add more depth to F182M and F210M; on the other hand, it gave us the opportunity to include F444W in our analysis.

All JWST/NIRCam images have been reduced by adopting a modified version of the official JWST pipeline<sup>24</sup> (based on *jwst*1.8.2 and Calibration Reference Data System pipeline mapping (CRDS; pmap) 1018). More detailed information about the reference files is available on the official STScI/CRDS website.<sup>25</sup>

Compared to the official JWST pipeline, our version includes different procedures, following some of the ideas presented in Bagley et al. (2023), to deal with the unresolved problems that still affect the official software. In our data reduction, we minimized the impact of the so-called “snowballs,” the 1/f noise, the “wisps,”<sup>26</sup> and the residual cosmic rays. After reducing all the JWST/NIRCam images from Williams’s and FRESCO programs, we drizzled all the NIRCam calibrated files to 0.03 pixel<sup>-1</sup>, as the final pixel scale we adopted in this work. All the final images have been aligned to the Hubble Legacy Fields (HLF) catalog.<sup>27</sup>

As a sanity check, we compared the photometry for the brightest sources (<24 mag) in all the NIRCam filters. To do that, we produced two versions of our final images, with and without the extra steps we employed in our modified version of the official pipeline. Then, we extracted the sources by using the software SOURCE EXTRACTOR (SEXTRACTOR; Bertin & Arnouts 1996) and compared their photometry. This test demonstrated that our extra steps do not introduce any kind of systematic effect in the photometry.

### 2.1.2. JWST/MIRI

We complemented the JWST/NIRCam observations with the MIRI 5.6 µm imaging from the JWST Guaranteed Time Observations (GTO) program: MIRI Deep Imaging Survey (MIDIS; PID: 1283, PI: Göran Östlin). The MIRI observations were carried out in 2022 December and targeted with the broadband filter F560W the HUDF for a total amount of 50 hr (≈41 hr on-source), covering an area of about 4.7 arcmin<sup>2</sup>. By reaching a median depth of 29.15 mag (5σ, r = 0.15), this set of observations represents the deepest imaging available at 5.6 µm to date. A complete description of the data collection and reduction, as well as the source statistics on these 5.6 µm images, will be presented by Östlin et al. (G. Östlin et al. 2023,

in preparation). Here we only summarize the basic information of this data processing.

As in the case of the NIRCam imaging, we adopted a modified version of the official JWST pipeline to reduce the MIRI data. In fact, the final products that can be obtained by running the JWST pipeline are still affected by strong patterns (e.g., vertical striping and background gradients) that impact the scientific quality of the images (e.g., Iani et al. 2022). To overcome these problems, we added to the pipeline some extra steps at the end of stages 2 and 3 that allowed us to significantly mitigate the intensity of the striping, the background inhomogeneities as well as the noise of the output image. A comparison between the F560W magnitude of the brightest galaxies (<24 mag) measured in MIRI images obtained with and without the extra steps ensured that our modified version of the pipeline did not introduce any systematic offset.

Finally, we drizzled the final MIRI image to the same pixel scale adopted for the JWST/NIRCam images and registered its astrometry to the HLF catalog.

### 2.1.3. Ancillary HST Data

We obtained all of HST images over the HUDF from the Hubble Legacy Field GOODS-S (HLF-GOODS-S).<sup>28</sup> The HLF-GOODS-S provides 13 HST bands covering a wide range of wavelengths (0.2–1.6 µm), from the UV (WFC3/UVIS F225W, F275W, and F336W filters), optical (ACS/WFC435W, F606W, F775W, F814W, and F850LP filters), to near infrared (WFC3/IR F098M, F105W, F125W, F140W, and F160W filters). See Whitaker et al. (2019) for more detailed information on these observations.

## 2.2. Photometric Analysis

We used the software SEXTRACTOR to detect the sources and measure their photometry in all the 20 filters available from the HST and JWST, covering a wide range of wavelengths (0.2–5.6 µm). We used SEXTRACTOR in dual-image mode adopting a superdetection image that we created by combining photometric information from different bands. In order to maximize the number of the detected sources, we opted to use a hot-mode extraction, as presented in Galametz et al. (2013), which is well suited to find very faint sources.

We combined aperture photometry, adopting circular apertures (i.e., MAG_APER) of 0.5 diameter, and Kron apertures (i.e., MAG_AUTO, Kron 1980) following the same prescription we adopted in Rinaldi et al. (2022, see Section 3.2). We chose a circular-aperture flux over a Kron flux when the sources were fainter than a given magnitude. In this case, as we were dealing with very deep images, we decided to consider mag<sub>lim</sub> = 27 as our faint limit for the Kron aperture. This final decision has been taken after several tests we performed with the HST photometry, comparing our fluxes with the HLF photometric catalog from Whitaker et al. (2019). We corrected the aperture fluxes to the total. For the HST, these corrections are well known.<sup>29,30,31</sup> For the JWST, instead, we estimated the aperture corrections using the software WEBBPSF.<sup>32</sup>

<footnote n="24"><sup>24</sup> The pipeline is available at the following link.</footnote>

<footnote n="25"><sup>25</sup> http://jwst-crds.stsci.edu</footnote>

<footnote n="26"><sup>26</sup> More information about these artifacts at the following link.</footnote>

<footnote n="27"><sup>27</sup> The HLF catalog is available at the following link.</footnote>

<footnote n="28"><sup>28</sup> The HST images (0.03 pixel<sup>-1</sup>) have been downloaded from the following link.</footnote>

<footnote n="29"><sup>29</sup> Aperture corrections for HST/ACS.</footnote>

<footnote n="30"><sup>30</sup> Aperture corrections for HST/WFC3-IR.</footnote>

<footnote n="31"><sup>31</sup> Aperture corrections for HST/WFC3-UVIS.</footnote>

<footnote n="32"><sup>32</sup> The software WEBBPSF is available at the following link.</footnote>
Moreover, we adopted a minimum error of 0.05 mag for all the HST photometry because SExTRACTOR typically underestimates photometric errors (e.g., Sonnett et al. 2013). We decided to adopt this minimum error value for the JWST images as well to account for possible uncertainties in the NIRCam and MIRI flux calibrations.

Finally, all our fluxes have been corrected for Galactic extinction. Those values have been estimated adopting a python package called `DUSTMAPS`.<fnref n="33" /> As a sanity check, we compared the correction factors for the HST filters with Schlafly & Finkbeiner (2011), finding an excellent agreement with the values we can recover following their prescription, as expected.

### 2.3. SED Fitting

We performed the SED fitting and derived the properties of our sources by making use of the code LEPHARE (Arnouts & Ilbert 2011). We constructed the libraries for LEPHARE by adopting the same configuration we used in Rinaldi et al. (2022, see Section 4). Briefly, we considered the stellar population synthesis (SPS) models proposed by Bruzual & Charlot (2003, hereafter BC03), based on the Chabrier IMF (Chabrier 2003). We made use of two different star formation histories (SFHs): a standard exponentially declining SFH (known as "τ-model") and an instantaneous burst adopting a simple stellar population (SSP) model. In particular, we adopted two distinct metallicity values, a solar metallicity ($Z_⊙ = 0.02$) and a fifth of solar metallicity ($Z = 0.2Z_⊙ = 0.004$). Moreover, to take the strong contribution from the nebular emission lines that can occur at very young ages into account, we also considered STARBURST99 templates (Leitherer et al. 1999, hereafter SB99) for young galaxies (age $\le 10^7$ yr) with constant star formation histories. We considered the Calzetti et al. (2000) reddening law in combination with Leitherer et al. (2002) to better constrain wavelengths below 912 Å. In particular, we adopted the following color excess values: $0 \le E(B-V) \le 1.5$, with a step of 0.1. We also decided to run LEPHARE between $z=0$ and $z=20$, by considering the following steps: $\Delta z=0.04$ between $z=0$ and $z=6$ and $\Delta z=0.1$ between $z=6$ and $z=20$ (291 steps in total). We summarize the parameters we adopted to perform the SED fitting in Table 1.

We estimated upper limits for each source that SExTRACTOR was not able to detect. To do that, around each source, we placed random circular apertures ($0.5''$ diameter) to estimate the background rms ($1\sigma$). For LEPHARE, we opted to use the $3\sigma$ upper limit for the flux in those filters where we did not have a detection. Finally, for all those sources for which we did not have any photometric information (e.g., the MIRI/F560W and NIRCam coverage areas are different), we simply ignored those filters during the SED fitting (i.e., we used $-99$ as input flux in LEPHARE).

### 3. Selection of Strong ($H\beta + [OIII]$) and $H\alpha$ Emitters at $z \simeq 7-8$

LEPHARE returns the best-fit SED and derived parameters for each source. We performed two different runs with LEPHARE, one adopting BC03 models only and the other one adopting SB99 models only. Therefore, we created the final catalog choosing for each source the best $\chi^2$ between

<table><caption>**Table 1**
Parameters Used to Perform the SED Fitting with LEPHARE by Adopting BC03 and SB99 Models</caption><thead><tr><th colspan="3">Parameter</th></tr></thead><tbody><tr><td>Templates</td><td>Bruzual & Charlot (2003)</td><td>Leitherer et al. (1999)</td></tr><tr><td>e-folding time (τ)</td><td>0.01-15 (8 steps) + SSP</td><td>Constant SFH</td></tr><tr><td>Metallicity (Z)</td><td>0.004; 0.02 (=Z<sub>☉</sub>)</td><td>0.008; 0.001</td></tr><tr><td>Age (Gyr)</td><td>0.001-13.5 (49 steps)</td><td>0.001-0.1 (6 steps)</td></tr><tr><td colspan="3">Common values</td></tr><tr><td>Extinction laws</td><td>Calzetti et al. (2000) + Leitherer et al. (2002)</td><td></td></tr><tr><td>E(B-V)</td><td>0-1.5 (16 steps)</td><td></td></tr><tr><td>IMF</td><td>Chabrier (2003)</td><td></td></tr><tr><td>Redshift</td><td>0-20 (291 steps)</td><td></td></tr><tr><td>Emission lines</td><td>Yes</td><td></td></tr><tr><td>Cosmology (H<sub>0</sub>, Ω<sub>0</sub>, Λ<sub>0</sub>)</td><td>70, 0.3, 0.7</td><td></td></tr></tbody></table>

**Note.** For the run with SB99 models, we used the same configuration as for the BC03 models for the extinction law, $E(B-V)$, IMF, redshift interval, and cosmology. Moreover, for the run with SB99, we opted for only six steps in age because nebular emission lines only matter for very young ages.

the BC03 and SB99 solutions. Finally, we cleaned our catalog of possible stars. To do so, we first cross-matched our catalog with Gaia Data Release 3 (DR3; Babusiaux et al. 2023). Then, we looked at the stellarity parameter (i.e., CLASS_STAR) we have from SExTRACTOR. In particular, we applied the same criterion adopted in Caputi et al. (2011, Section 3.1). We removed all those sources that have CLASS_STAR > 0.8 and occupy the stellar locus in the (F435W - F125W) versus (F125W - F444W) color-color diagram. In total, less than $\simeq 1\%$ sources have been discarded from our full catalog because they have been classified as stars (eight of them have been identified in GAIA DR3).

As our goal is to look for potential ($H\beta + [OIII]$) and ($H\alpha + [NII]+[SII]$) emitters in the XDF at $z \simeq 7-8$, we only focused on those sources for which the best photometric redshift falls in that redshift range.

For each candidate, we created postage stamps to make a careful visual inspection in order to exclude all those galaxies that either fall on stellar spikes or are heavily contaminated by the light of the nearby sources. After this visual inspection, we were left with 58 robust galaxy candidates at $z \simeq 7-8$.

Among these sources, we searched for ($H\beta + [OIII]$) and $H\alpha$ emitters. We first analyzed if they show a flux excess in the following three bands: NIRCam/F430M, NIRCam/F444W, and MIRI/F560W. The first two filters have been used to look at the flux enhancement produced by ($H\beta + [OIII]$). In turn, MIRI/F560W has been used to look at the flux excess produced by $H\alpha$.

To convert the flux excess into an $EW_0$ we followed the canonical approach described by Mármol-Queraltó et al. (2016). Following that procedure, we know that

$$EW_0 = \frac{W_{\text{rec}}}{1+z} (10^{(-0.4\Delta\text{mag})} - 1), \quad (1)$$

where $W_{\text{rec}}$ is the rectangular width of the filter containing the emission line in question, in our case ($H\beta + [OIII]$) or $H\alpha$, and $\Delta\text{mag}$ is the difference between the observed magnitude in that filter and the synthetic magnitude<fnref n="34" /> from the SED fitting (i.e.,

<footnote n="33">³³ The DUSTMAPS python package is available at the following link.</footnote>

<footnote n="34">³⁴ For each galaxy, LEPHARE returns the synthetic magnitude in each filter (i.e., $mag_{syn}$) for the best-fit model along with the stellar parameters.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Postage stamps ($5'' \times 5''$) of one of our ($H\beta$ + [O III]) candidates (ID: 9434, $z_{\text{best}} = 7.68_{-0.01}^{+0.03}$). The last postage stamp refers to the stacked image we adopted as the detection map with *SEXTRACTOR*. Here we show all the bands we used in our analysis, from 0.2 to 5.6 µm. The green circle has been placed to only guide the eye on the source. In particular, from these postage stamps, there is a clear excess at 4.3 µm. This source shows an excess in MIRI/F560W as well.</figcaption></figure>

$\Delta\text{mag} = \text{m}_{\text{obs}} - \text{m}_{\text{syn}}$) that we adopted as a proxy for the continuum emission.

Therefore, to estimate the flux excess, we assumed that the continuum flux was well described by the synthetic NIRCam/F460M obtained from the best-fit template for each galaxy. In particular, we selected all those sources for which $|\text{mag}_{\text{obs}}(\text{F460M}) - \text{mag}_{\text{syn}}(\text{F460M})| \le 2 \times \text{mag}_{\text{err}}(\text{F460M})$, where $\text{mag}_{\text{obs}}$ and $\text{mag}_{\text{syn}}$ are the observed and best-fit synthetic magnitudes, respectively. This condition ensures that the continuum at 4.6 µm can be considered flat within the error bars. We also double-checked if this condition was satisfied in NIRCam/F480M.

Once we selected all those sources that survive the condition described above, we estimated the flux excess in the following way: $\Delta\text{mag} = (\text{mag}_X - \text{mag}_{\text{cont}})$, where mag$_X$ represents the magnitude in one of the filters we chose to select ($H\beta$+[O III]) or H$\alpha$, and mag$_{cont}$ refers to F460M$_{syn}$. We highlight that this selection is purely based on the photometric excess we considered above. None of our derivations is based on emission lines modeled by LEPHARE. For a conservative approach, we only considered those galaxies for which the flux excess with respect to the stellar continuum satisfies the following condition: $\Delta\text{mag} < -0.2$. Note that a $\Delta\text{mag} = -0.2$ in NIRCam/F430M corresponds to a $EW_0 \approx 58$ Å at $z = 7$, while in NIRCam/F444W it would imply an $EW_0 \approx 270$ Å. For MIRI/F560W, the same $\Delta\text{mag}$ would correspond to an $EW_0 \approx 239$ Å at the same redshift.

We inspected again the postage stamps of the 58 possible candidates, after estimating the flux excess in each band (NIRCam/F430M, NIRCam/F444W, and MIRI/F560W), to make a cross-match between the values we got for $\Delta\text{mag}$ and the visual inspection of the sources themselves. We also examined the best-fit SED for each galaxy. This safely allowed us to conclude that 18 sources can be securely classified as ($H\beta$ + [O III]) emitters. These emitters constitute $\approx 31\%$ of our total galaxy sample at $z \simeq 7-8$ (see Figure 2 where we show the multiwavelength images of an example source). The derived $EW_0$ values cover a wide range that goes from a minimum of

$87.5_{-27}^{+30}$ Å to a maximum value of $2140.4_{-154}^{+970}$ Å, with a median $\langle EW_0 \rangle \simeq 943_{-194}^{+737}$ Å (lower and upper errors refer to the 16th and 84th percentiles). This value is higher, but still marginally consistent with the error bars, than that derived by Labbé et al. (2013) from Spitzer Space Telescope observations of bright $z \simeq 8$ galaxy candidates. Out of the 18 ($H\beta$ + [O III]) emitters, 83% have a best-fit SED with subsolar (0.2 $Z_\odot$) metallicity and the remaining $\simeq 17\%$ with solar ($Z_\odot$) metallicity.

Among the 18 ($H\beta$ + [O III]) emitters at $z \simeq 7-8$, a total of 16 lie on the ultradeep MIRI 5.6 µm coverage field. Out of them, 12 show a significant 5.6 µm flux excess with respect to the continuum (as defined above), which we interpret as the presence of the ($H\alpha$+[N II]+[S II]) line complex at $z \simeq 7-8$. To obtain the net value of the $H\alpha$ $EW_0$, we applied the correction recipes provided by Anders & Fritze-v. (2003), as follows: $f(H\alpha) = 0.63f(H\alpha + [N\text{II}] + [\text{S II}])$ for a solar metallicity, and $f(H\alpha) = 0.81f(H\alpha + [N\text{II}] + [\text{S II}])$ for a 0.2 $Z_\odot$ metallicity. Note that with this procedure we are assuming that the stellar and gas metallicities are similar in these galaxies.

We also compared the derived stellar properties from the SED fitting between the ($H\beta$ + [O III]) and H$\alpha$ emitters and nonemitters. Performing the two-sample Kolmogorov–Smirnov test, we do not find any significant difference between the two samples in terms of ages, $E(B-V)$, metallicity, and stellar mass. Regarding the SFR$_{best}$ distributions, we see a difference between the two populations (SFR$_{best}$ for the emitters tend to be higher than SFR$_{best}$ for the nonemitters) that might reflect the fact that we are looking at strong emitters (i.e., SFR is higher). We show these distributions in Figure 3.

## 4. Results

Once we estimated the stellar properties of our candidates by performing the SED fitting with LEPHARE, we analyzed the properties of these sources by comparing our results with the recent literature at high redshifts. Before doing that, we first ensured that the stellar masses we inferred with LEPHARE were not affected by the presence of the flux excess we estimated in
<figure><img src="image_5.png" /></figure>

<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 3.** Comparison of the best-fit properties for emitters ((Hβ + [O III]) and Hα) and nonemitters at z = 7–8: stellar mass and age (upper row); star formation rate and color excess (middle row); and metallicity (bottom row). No significant differences have been noticed between the two populations for most of the stellar parameters, as determined by performing a two-sample Kolmogorov–Smirnov test. Differences in SFRbest between emitters and nonemitters might be explained by the fact that we are only looking at strong emitters that show a higher SFRbest.
<figure><img src="image_10.png" /><figcaption>**Figure 4.** Best-fit SEDs for two examples of line-emitter candidates at $z \simeq 7-8$. On the left panel, we show a source at $z_{phot} \simeq 7.19$ (ID = 9432). On the right panel, we show another source at $z_{phot} \simeq 7.68$ (ID = 9434; shown in Figure 2). Both panels show how well-constrained the best-fit SEDs and the derived photometric redshifts are, which is evident by simple inspection of the best-fit templates and their probability density functions (PDF(z)). In each case, we notice the clear presence of an excess in F430M, F444W, and F560W, which we adopted as the criterion to select our sample of ($H\beta$ + [O III]) and $H\alpha$ emitters.</figcaption></figure>

F430M, F444W, and F560W. To do that, we rerun LEPHARE following the methodology explained by Caputi et al. (2017). This time, for each source, we turned off those bands (NIRCam/F430M, NIRCam/F444W, and MIRI/F560W) in which we found a flux excess (i.e., -99 following LEPHARE's prescription). Moreover, for this run, we fixed the redshifts adopting the photometric ones we estimated from the original run. Doing this test allows us to ensure that our stellar mass estimates are not affected by any emission line that falls in one of those filters. We found a good agreement within $2\sigma$. Finally, we also inspected that the stellar continuum was well described by inspecting the best-fit SEDs we obtained from LEPHARE. In Figure 4 we show two examples (ID: 9432, 9434) of the best-fit SEDs for the candidates we have in our sample.

## 4.1. Emission Line EW Versus Stellar Mass and Age in Galaxies at $z \simeq 7-8$

Having calculated the ($H\beta$ + [O III]) $EW_0$ for the prominent line emitters, we can compare their best-fit SED properties with those of the other $z \simeq 7-8$ galaxies in our sample. In Figure 5, we show the derived ($H\beta$ + [O III]) $EW_0$ versus the best-fit age. From this plot, we can see that all except three of the ($H\beta$ + [O III]) emitters are characterized by young best-fit ages ($\le 10^8$ yr), which indicates that these objects may be in their first major star formation episode. The remaining three objects are older ($>10^8$ yr), with two having almost the age of the universe at their redshifts. This fact suggests that these galaxies could be having a rejuvenation episode, as is known to happen at lower redshifts (Rosani et al. 2020), as it is unlikely that they could have sustained their high instantaneous SFR values for all of their lifetimes.

The gray triangles in Figure 5 refer to the EW upper limits that we estimated for all those galaxies at $z \simeq 7-8$ that do not have a significant flux excess in the NIRCam/F430M band. In contrast to the ($H\beta$ + [O III]) emitters, the nonemitters span different possible ages at those redshifts, without any bias toward young/old ages.

We also compared our results with the recent literature. In particular, Endsley et al. (2021) studied a sample of 20 rest-frame ultraviolet (UV) bright ($H\beta$ + [O III]) emitters at

<figure><img src="image_11.png" /><figcaption>**Figure 5.** Age vs. ($H\beta$ + [O III]) $EW_0$ for our line emitters at $z = 7-8$. The gray triangles refer to the $EW_0$ upper limits that we estimated for all the “nonemitter” galaxies in our sample. We also show the data points from Endsley et al. (2021) to make a comparison with the recent literature at high redshift, albeit in a much higher luminosity regime. The curves refer to the evolution of the $H\beta$ $EW_0$ as a function of age expected from SB99 models, corresponding to the two metallicities (solar and subsolar) that we have considered in our work and for two different SFHs. A clear anticorrelation between $EW_0$ and age is evident in this plot, which is in line with previous findings in the literature at lower redshifts (Reddy et al. 2018).</figcaption></figure>

$z \simeq 6.8-7$ that have been selected over a wide sky area (2.7 deg$^2$ in total). Endsley et al. (2021) found this rare population of very strong ($H\beta$ + [O III]) emitters with an $EW_0$ >1200 Å. The fact that we find similarly high ($H\beta$ + [O III]) $EW_{S0}$ among faint galaxies in a much smaller area of the sky indicates that prominent ($H\beta$ + [O III]) emitters were much more common at the EoR than what can be inferred from the brightest galaxies.

Finally, the solid and dashed lines in Figure 5 show the expected variation of the $H\beta$ (only) $EW_0$ versus age for SB99 model galaxies. These theoretical tracks are based on a Chabrier IMF with a stellar mass cutoff of 100$M_\odot$ and were obtained both for a solar and a subsolar metallicity (0.2$Z_\odot$),
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Stellar mass vs. (Hβ + [O III]) EW₀. The gray triangles refer to the upper limits we estimated for all those galaxies we classified as “nonemitters” during our selection. Also in this case, we report data points from Endsley et al. (2021) to make a comparison with the recent literature at high redshift. The EW₀ broadly anticorrelates with stellar mass, similarly to what has been reported by Reddy et al. (2018) and Endsley et al. (2021) at lower redshifts.</figcaption></figure>

each for a single burst and constant SFH. As expected, our data points are located nicely above these curves, following the trend of the models with constant star formation histories, albeit with higher EW₀, due to the [O III] contribution.

Over the past decades, the recombination line equivalent widths have been used as proxies for stellar population age in star-forming galaxies. The ratios between the fluxes of the recombination line, which are sensitive to the instantaneous star formation rates (SFRs), and the fluxes of the continuum, which are sensitive to the previous average SFR, are indeed what we define as recombination line equivalent widths (Stasińska & Leitherer 1996). In particular, Reddy et al. (2018) found a very strong anticorrelation between (Hβ + [O III]) EW₀ and young ages at z ≈ 1.8–3.8, which does not evolve as a function of redshift at that range of cosmic time. By looking at Figure 5, we can see that this anticorrelation is evident also at z ≈ 7–8 where strong (Hβ + [O III]) emitters prefer young ages, which is in line with what has been found at lower redshifts. We double-checked this result by estimating the Spearman’s rank correlation coefficient, finding that those two quantities anticorrelate (i.e., Spearman’s coefficient ≈ −0.5) with a p-value ≈ 0.03. Therefore, we can conclude that there is evidence of a moderate anticorrelation between age and EW₀(Hβ + [O III]).

We repeated the same exercise looking, this time, at the derived (Hβ + [O III]) EW₀ versus stellar mass for our (Hβ + [O III]) emitters (Figure 6). Also, in this case, the stellar masses come directly from the best-fit SED obtained with LEPHARE. As we can see from Figure 6, our (Hβ + [O III]) emitters have a stellar mass that ranges from a minimum value of log₁₀(M*/M☉) ≈ 7.5 to a maximum value of log₁₀(M*/M☉) ≈ 9. In previous works, it has been shown that the normalization of the (Hβ + [O III]) EW₀ versus stellar mass relation should increase with redshift (e.g., Reddy et al. 2018). Here we find a broad anticorrelation between the two quantities. The gray triangles in Figure 6 refer to the upper limits that we estimated for the (Hβ + [O III]) EW₀ for the z ≈ 7–8 galaxies that are not classified as emitters from a NIRCam flux excess.

<figure><img src="image_13.png" /><figcaption>**Figure 7.** (Hβ + [O III]) EW₀ vs. Hα EW₀. Here we color coded our data points for the [O III]λ5007/Hβ flux ratio. Most galaxies have [O III]λ5007/Hβ > 1, and they mostly lie on or below the identity line. Among the four galaxies that lie above the identity line, two have [O III]λ5007/Hβ < 1, i.e., these line ratios are dominated by Hβ. Instead, the other two galaxies above the identity line have [O III]λ5007/Hβ > 1 and correspond to cases with nonnegligible dust extinction.</figcaption></figure>

Finally, for all those galaxies that show an “Hα excess,” we compare their (Hβ + [O III]) EW₀ versus Hα EW₀, where the “Hα excess” has been corrected to only take the real Hα flux into account following Anders & Fritze-v. (2003). We show this comparison in Figure 7. In particular, we also plot the recent results from Prieto-Lyon et al. (2023) where they inferred those quantities studying a sample of galaxies at z ≈ 3–7. We see that our sample is in good agreement with the expected correlation that has been found in Prieto-Lyon et al. (2023) as well. As a matter of fact, as we can derive the Hα line flux from our data, we can also infer the Hβ line flux independently and separate the contributions of Hβ and [O III] for each galaxy, considering the following:

$$f(\mathrm{H}\beta) = f(\mathrm{H}\alpha) \times 10^{-0.4 \times 1.27 E(B-V)/2.86}, \quad (2)$$

where f(Hα) and f(Hβ) refer to the observed fluxes and E(B-V) is the color excess obtained from the best-fit SED model. The denominator 2.86 corresponds to assuming case-B recombination (e.g., Osterbrock & Ferland 2006), while the factor −1.27 = k(Hα) − k(Hβ) is obtained from the Calzetti et al. (2000) reddening law.

Once we know the Hβ flux for each source, we can independently work out the [O III]λλ 4959, 5007 fluxes for all those emitters that show an Hα excess.

The data points in Figure 7 are color coded according to each galaxy’s [O III]λ5007/Hβ ratio. From that figure, we see that most line emitters have [O III]/Hβ > 1, indicating the predominance of [O III], which is consistent with recent literature findings at similar redshifts. Instead, two galaxies have [O III]/Hβ < 1, i.e., the Hβ line flux is larger than the [O III] line flux for them. These two galaxies are well above the identity line in Figure 7, as expected. We separate the Hβ and [O III]λ5007 line fluxes simply assuming the case-B recombination Hα/Hβ = 2.86 ratio and the corresponding color excess mentioned above. The [O III]/Hβ < 1 values could indicate very low metallicities, but this would need to be confirmed with a spectroscopy follow up of these sources.
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Stellar mass vs. SFR. Here we show the SFR-M* plane populated by the SFR directly inferred from the “Hα excess.” To make a comparison with the recent literature at high redshifts, we plot data points from Rinaldi et al. (2022) that give us the opportunity to populate this plane with very low-mass galaxies at z ≃ 2.8–6.5. We also show data points from Endsley et al. (2021), who studied a sample of 20 bright (Hβ + [O III]) emitters at z ≃ 6.8–7, and indicate the starburst zone, as defined by Caputi et al. (2017, 2021). We also plot the expected MS of galaxies at z ≃ 7–8 from Speagle et al. (2014). Our data points are color coded by their [O III]/Hβ ratio. We see no correlation between this ratio and the position of sources on the SFR-M* plane.</figcaption></figure>

## 4.2. Hα-derived SFR and the Location of Galaxies on the SFR-M* Plane

For all the 12 Hα emitters at z ≃ 7–8 as determined from the MIRI 5.6 µm imaging, we estimated their SFRs from their inferred Hα luminosities.

After we obtained the net observed Hα flux for each source, we converted those fluxes into the intrinsic ones by simply applying the Calzetti reddening law. We then estimate the luminosity for the Hα emission line and apply the following formula from Kennicutt (1998) to obtain the corresponding SFR(Hα):

$$ \text{SFR}(M_{\odot} \text{ yr}^{-1}) = 7.9 \times 10^{-42} L_{H\alpha} (\text{erg s}^{-1}). \quad (3) $$

As the aforementioned formula has been originally calibrated for a Salpeter IMF over (0.1–100) M* (Salpeter 1955), we applied a conversion factor (Madau & Dickinson 2014; i.e., 1.55) to rescale it to a Chabrier IMF (Chabrier 2003).

We then placed our sources on the SFR-M* plane, as we show in Figure 8. To make a comparison with the recent literature, we also populated this plane with star-forming galaxies at z ≃ 3.0–6.5 from Rinaldi et al. (2022) and (Hβ + [O III]) emitters at z ≃ 6.8–7 (Endsley et al. 2021). We also indicate the starburst (SB) zone as determined in Caputi et al. (2017, 2021), which empirically defined as starburst galaxies all those sources with sSFR > 10⁻⁷.⁶⁰ yr⁻¹.

We see that five (≈42%) of the galaxies that show an “Hα excess” lie in the starburst zone, while only two are located on the star formation main sequence (MS; Brinchmann et al. 2004; Noeske et al. 2007; Peng et al. 2010; Speagle et al. 2014; Rinaldi et al. 2022). The remaining five galaxies appear close, but slightly below the starburst envelope, in what has been defined in Caputi et al. (2017) as the star formation valley (SFV), i.e., in between the starburst cloud and the MS, suggesting that they are on the way to/from a starbursting phase. The fact that the vast majority of emitters are in or close

<figure><img src="image_15.png" /><figcaption>**Figure 9.** (a) Comparison between SFRUV and SFRHα. The error bars reflect the usual scatter that has been observed with the Kennicutt’s relations we used to derive those two quantities. (b) The ratio of SFRHα and SFRUV as a function of stellar mass. Both SFRs have been corrected by adopting the same reddening curve (Calzetti et al. 2000). The horizontal line indicates a one-to-one ratio. The pale blue shade refers to Atek et al. (2022) results at lower redshifts.</figcaption></figure>

to the starburst zone is consistent with the findings of Endsley et al. (2021) for brighter galaxies, as it can be seen in Figure 8.

We also color coded our Hα emitters according to their [O III]λ5007/Hβ ratios. We find no correlation between these ratios and the position of galaxies on the SFR-M* plane.

For the Hα sample in Figure 9(a) we show the comparison between the two different SFR indicators that we considered in this paper (UV and Hα luminosities). From that plot, we clearly see differences between those two indicators (SFRHα and SFRUV). This finding is not surprising as it has been already pointed out in the literature (e.g., Flores Velázquez et al. 2021; Atek et al. 2022; Patel et al. 2023).

Differences between these two SFR tracers (Figure 9(a)) may be partly explained by uncertainties in the dust-extinction correction, which mostly affect the UV continuum fluxes, and by our assumption that the dust extinction of the continuum and emission lines is the same and only depends on wavelength. However, part of the scatter observed in the SFRHα and SFRUV plane may be real and due to the following:
1. In very young galaxies (below 100 Myr), SFR<sub>UV</sub> underestimates the real value because the UV luminosity associated with star formation is still growing. Indeed, when comparing SFR<sub>H&alpha;</sub> and SFR<sub>UV</sub>, one has to also take age effects into account. UV traces typically 1500–2000 &Aring; (i.e., nonionizing photons), while H&alpha; traces directly <912 &Aring; photons. For example, UV-bright regions without H&alpha; emission trace the presence of star-forming clumps dominated by B-type stars and where most massive O-type have already evolved;

2. Different ionizing photon production efficiencies (e.g., Nanayakkara et al. 2020; Endsley et al. 2023; P. Rinaldi et al. 2023, in preparation).

Finally, by exploiting the FIRE simulations (Hopkins et al. 2014), Sparre et al. (2017) showed that the Hα measurement of the SFR over a short timescale can fluctuate significantly, up to a factor of ten, compared to the UV indicator.

Following Atek et al.'s (2022) procedure at lower redshifts, in Figure 9(b) we inspected the ratio between SFR<sub>Hα</sub> and SFR<sub>UV</sub> as a function of the stellar mass. We find similar results as Atek et al. (2022, see their Figure 8) where the ratio of SFR<sub>Hα</sub>/SFR<sub>UV</sub> seems to be generally higher for the low-mass galaxies. Similarly, Faisst et al. (2019) found that more than 50% of their sample has SFR<sub>Hα</sub> in excess compared to SFR<sub>UV</sub>, particularly in low-mass galaxies. However, there are still uncertainties in determining the ratio of SFR<sub>Hα</sub>/SFR<sub>UV</sub> and how it changes with different galaxy parameters. As we know from the literature, the SFR indicators use conversion factors from Hα and UV luminosities, which assume that the SFR is constant. Nonetheless, this assumption may not be that accurate for different SFHs, especially when we consider cases of bursty star formation.

### 4.3. The Role of the Hα Emitters in the Cosmic Star Formation History at $z \simeq 7-8$

With the SFR values derived in the previous section, we computed the contribution of the prominent Hα emitters to the cosmic star formation rate density (SFRD) at $z \simeq 7-8$. To do that, we sum up the individual SFRs ($SFR_{H\alpha,\text{total}} \simeq 51.39 M_{\odot} \text{yr}^{-1}$) and then divide the total by the comoving volume<fnref n="35" /> encompassed by the area ($A \simeq 4.7 \text{ arcmin}^2$) and redshift bin (i.e., $z \simeq 7-8$) analyzed in this work ($V_{\text{sky}} \simeq 11580.26 \text{ Mpc}^3$). We obtain that, at these redshifts, the Hα emitters make for $\log_{10}(\rho_{SFR_{H\alpha}}/(M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3})) \simeq -2.35 \pm 0.3$.

In Figure 10 we show the redshift evolution of the SFRD as proposed by Lilly et al. (1996) and Madau et al. (1996), the so-called “Lilly–Madau diagram.” In this plot, we show our own estimation of the SFRD, along with a compilation of recent results from the literature based on different SFR tracers. In particular, we also show the SFRD values that have recently been obtained by Bouwens et al. (2023) using JWST data, tracing the SFR directly from the UV continuum emission at $z \simeq 9$ to $z \simeq 15$ as well as Pérez-González et al. (2023) results at $z \simeq 8-13$ from ultradeep NIRCam images in HUDF-P2 (PID proposal: 1283, PI: Göran Östlin). Our inferred SFRD appears to be in good agreement with what has been found in the literature at similar redshifts. We also find a very good agreement with the predictions from theoretical models (e.g.,

IllustrisTNG; Springel et al. 2018). In particular, in Figure 10 we also show the total SFRD, which has been estimated from both Hα emitters and nonemitters at $z \simeq 7-8$ ($\log_{10}(\rho_{SFR_{\text{tot}}}/(M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3})) \simeq -1.76 \pm 0.3$).<fnref n="36" /> For the nonemitters, the SFR has been obtained from the rest-frame UV continuum luminosity at 2000 Å and adopting the conversion formula from Kennicutt (1998).

### 4.4. The Evolution of the Rest-frame EW(Hα) As a Function of the Redshift

Finally, our derived values of the Hα EW<sub>0</sub> allow us to extend the study of the redshift evolution of this parameter to $z \simeq 7-8$. In Figure 11 we present our results along with the most recent determinations from the literature (for sources at $z \simeq 0.5-6$, Erb et al. 2006; Shim et al. 2011; Fumagalli et al. 2012; Stark et al. 2013; Sobral et al. 2014; Mármol-Queraltó et al. 2016; Smit et al. 2016; Reddy et al. 2018; Lam et al. 2019; Atek et al. 2022; Boyett et al. 2022; Sun et al. 2022; Ning et al. 2023) and a stacking analysis measurement by Stefanon et al. (2022) at $z \simeq 8$. These previous works made use of different methods and techniques to determine the Hα EW, such as medium/high-resolution spectroscopy, low-resolution grism spectroscopy, and narrowband and broadband photometry combined with SED modeling, as we did in this paper.

Our sample of strong line emitters at $z \simeq 7-8$ allows us to populate a virtually unexplored part of parameter space. At those redshifts ($z \simeq 8$), only Stefanon et al. (2022) previously obtained an estimate of the average Hα EW<sub>0</sub>, by median stacking 102 Lyman-break galaxies (LBG) in the 3.6, 4.5, 5.8, and 8.0 µm bands from the Spitzer Infrared Array Camera (IRAC).

We also incorporate in the analysis the empirical prescriptions from Fumagalli et al. (2012) and Faisst et al. (2016), who predict that the EW<sub>0</sub>(Hα) should evolve differently below and above $z \simeq 2$. In particular, according to the recent literature, at $z < 2$, the EW<sub>0</sub>(Hα) should evolve as $\propto(1+z)^{1.8}$, while at $z > 2$ it should evolve as $\propto(1+z)^{1.3}$.

By inspection of Figure 11, we can see that JWST observations at $z \ge 6$ (i.e., Sun et al. 2022; Ning et al. 2023; this present work) suggest that the break proposed at $z \simeq 2$ in the past literature does not really hold up to such high redshifts (thin, black, and dashed line). For that reason, we fit the evolution of EW<sub>0(Hα)</sub> as a function of redshift again by considering the recent JWST observations at $z \ge 6$ as well. In this case, we find that EW<sub>0</sub>(Hα) $\propto (1+z)^{2.1}$ (bold, dark red, and dashed line in Figure 11). However, larger galaxy samples are needed to confirm this finding.

Our data points are in good agreement with the stacking estimate obtained by Stefanon et al. (2022). Some of these values are well above the empirical median extrapolation at those redshifts, while others are consistent with it. The prominent line emitters we analyze here constitute almost a quarter of all the MIRI-detected galaxies at $z \simeq 7-8$. The remaining MIRI sources at those redshifts should lie below the extrapolation of the empirical determination. This very large variation in the Hα EW<sub>0</sub> at $z \simeq 7-8$ suggests that, even at these very high redshifts, galaxies may be at different stages of their evolution, as we discuss in the next section.

<footnote n="35">³⁵ We estimated the comoving volume for the entire sky at $z \simeq 7-8$ by using the Cosmology calculator at the following link.</footnote>

<footnote n="36">³⁶ For the nonemitters, $A = 5.25 \text{ arcmin}^2$ (corresponding to the NIRCam coverage) and $V_{\text{sky}} = 12999.90 \text{ Mpc}^3$.</footnote>
<figure><img src="image_16.png" /></figure>

**Figure 10.** Cosmic star formation rate density as a function of the redshift. The large red circle at $log_{10}(\rho_{SFR_{H\alpha}}/(M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3})) \simeq -2.35$ indicates our estimate at $z \simeq 7-8$, which only takes the prominent $H\alpha$ emitters into account, i.e., it should be considered a lower limit to the real SFRD value at these redshifts. The red diamond at $log_{10}(\rho_{SFR_{tot}}/(M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3})) \simeq -1.76$, instead, refers to the total SFRD that we estimated accounting for both the $H\alpha$ emitters and nonemitters at $z \simeq 7-8$. For the nonemitters, the SFR directly comes from the UV continuum emission. Other symbols refer to the recent SFRD determinations from the literature, based on different SFR tracers (Schiminovich et al. 2005; Kajisawa et al. 2010; Burgarella et al. 2013; Gruppioni et al. 2013; Sobral et al. 2014; Bouwens et al. 2015; Rowan-Robinson et al. 2016; Caputi et al. 2017; Ishigaki et al. 2018; Bouwens et al. 2020; Gruppioni et al. 2020; Loiacono et al. 2021; Bouwens et al. 2023; Rinaldi et al. 2022; Pérez-González et al. 2023). The different curves correspond to theoretical predictions. Dashed line: Madau & Dickinson (2014). Solid line: Pillepich et al. (2018). All the SFRD values in this figure correspond to a Chabrier (2003) IMF.

**5. Summary and Conclusions**

In this paper, we have taken advantage of the publicly available medium-band and broadband NIRCam imaging in the XDF, combined with the deepest MIRI 5.6 µm imaging existing in the same field, to search for prominent ($H\beta$+[O III]) and $H\alpha$ emitters at $z \simeq 7-8$. This is the first time the $H\alpha$ emission line can be detected and its flux measured in individual galaxies at such high redshifts. This has been possible thanks to the unprecedented sensitivity of JWST observations, particularly those conducted with MIRI, for which the sensitivity gain is of more than an order of magnitude with respect to previous instruments operating at similar wavelengths (Iani et al. 2022).

We found 18 galaxies which are robust candidates to be prominent ($H\beta$+[O III]) emitters at $z \simeq 7-8$, as determined from their F430M and F444W flux excess. These 18 galaxies constitute $\simeq 31\%$ of all the galaxies that we find in the XDF in the same redshift range. Among them, 16 lie on the MIRI coverage area and 12 out of 16 have a clear flux excess in the MIRI/F560W filter, indicating the simultaneous presence of a prominent $H\alpha$ emission line. The ($H\beta$+[O III]) $EW_{S_0}$ that we derive range from $\simeq 87.5_{-27}^{+30} \AA$ to $2140.4_{-154}^{+970} \AA$, with a median value of $943_{-194}^{+737} \AA$. For most of these galaxies, we find [O III]/

$H\beta > 1$, but a few have [O III]/$H\beta < 1$. The two line fluxes can be separated by making use of the independent $H\alpha$ emission line measurement. This is telling us that some of the prominent ($H\beta$+[O III]) emitters likely have hard radiation fields typical of low-metallicity galaxies, but not all of them. Some are strong line emitters simply because they are intensively forming stars.

The identified $H\alpha$ emitters show an $EW_0$ that ranges from a few hundred to a few thousand Angstroms. Some of these values are substantially above the expected median $H\alpha$ EW at these redshifts, as extrapolated from lower redshift determinations. We also report that, by considering the recent JWST findings at $z \ge 6$ (including this present work), $EW_0(H\alpha)$ as a function of the redshift should evolve as follows: $EW_0(H\alpha) \propto (1+z)^{2.1}$. However, larger samples of galaxies are needed to confirm this result. We note, however, that the prominent $H\alpha$ emitters only constitute about a quarter of all the MIRI-detected galaxies at $z \simeq 7-8$. For the remaining galaxies, the $H\alpha$ $EW_0$ should lie below the expected median trend. As the $H\alpha$ EW is a good proxy for the sSFR, the lower EW values could indicate that these other galaxies (the nonemitters) have either relatively low SFRs, or a more important underlying stellar population producing a higher continuum. This is likely the case for the nonemitters at
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Evolution of the EW$_{0}$(Hα) as a function of redshift. After re-fitting all the data points, including our own, we find that the EW$_{0}$(Hα) evolution can be described by a single law: $EW_{0} (H\alpha) \propto (1+z)^{2.1}$ (bold, dark red, and dashed line). Our data points are color coded for sSFR. We also report the recent literature regarding the evolution of the EW$_{0}$(Hα) as a function of the redshift (thin, black, and dashed line). The gray shade represents a median estimate of the error bars of the data points from the literature. The orange shade represents the redshift window where the JWST is starting to detect these kinds of sources more systematically (e.g., Boyett et al. 2022; Sun et al. 2022; Ning et al. 2023). Note that the data point at $z \simeq 8$ from Stefanon et al. (2022) has been obtained by median stacking a sample of 102 Lyman-break galaxies in the Spitzer/IRAC bands from 3.6 to 8 µm.</figcaption></figure>

$z \simeq 7-8$, which are relatively evolved galaxies, with best-fit ages $>10^{7}-10^{8}$ yr and stellar masses $>10^{8} M_{\odot}$.

In turn, most of the prominent (Hβ+[O III]) and Hα emitters are characterized by higher sSFRs, with basically all of them being starburst galaxies or on the way to/from the starburst cloud. The majority of the prominent (Hβ+[O III]) emitters are very young galaxies (best-fit ages <10$^{7}$ yr), so they might be in their first major star formation episode. A few others are almost as old as the universe at their redshifts and have already built significant stellar mass ($>10^{8} M_{\odot}$), suggesting that they may be experiencing a rejuvenation effect.

Therefore, the overall conclusion of this work is that the galaxies present at the EoR are likely at different stages of their evolution. Furthermore, strong line emission is present in a minor, but significant fraction of sources.

Considering the Hα fluxes inferred for the prominent Hα emitters, we estimated their contribution to the cosmic SFRD at $z \simeq 7-8$. We found $\log_{10}(\rho_{SFR_{H\alpha}}/(M_{\odot} \text{ yr}^{-1} \text{ Mpc}^{-3})) \simeq -2.35 \pm 0.3$, in excellent agreement with independent measurements from the literature based on rest-frame UV luminosities, and with theoretical predictions and empirical extrapolations from lower redshifts. We note, however, that this estimated SFRD must be considered a lower limit, as it only takes into account the most prominent Hα emitters at $z \simeq 7-8$. We also considered the SFR$_{UV}$ for all the other galaxies at $z \simeq 7-8$ to obtain a total SFRD value at that redshift interval. We concluded that the strong Hα emitters produced about a quarter of the total SFRD at $z \simeq 7-8$, which suggests that they likely have had a significant role in the process of reionization. In a future paper, we will conduct a more detailed investigation of these sources, in order to better understand their nature.

## Acknowledgments

In memoriam to the MIRI European Consortium members Hans-Ulrik Noorgard-Nielsen and Olivier Le Fèvre.

The authors would like to acknowledge an anonymous referee for a careful reading and useful comments on this manuscript. This work is based on observations made with the NASA/ESA/CSA James WEBB Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for the JWST. These observations are associated with programs GO #1963, GO #1895, and GTO #1283. The authors acknowledge the team led by coPIs C. Williams, M. Maseda and S. Tacchella, and PI P. Oesch, for developing their respective observing programs with a zero-exclusive-access period. Also based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The specific observations analyzed can be accessed via the Hubble Legacy Fields at MAST, DOI: 10.17909/T91019. The work presented here is the effort of the entire MIRI team and the enthusiasm within the MIRI partnership is a significant factor in its success. MIRI draws on the scientific and technical expertise of the following organizations: Ames Research Center, USA; Airbus Defence and Space, UK; CEA-Irfu, Saclay, France; Centre Spatial de Liège, Belgium; Consejo Superior de Investigaciones Científicas, Spain; Carl Zeiss Optronics, Germany; Chalmers University of Technology, Sweden; Danish Space Research Institute, Denmark; Dublin Institute for Advanced Studies, Ireland; European Space Agency, Netherlands; ETCA, Belgium; ETH Zurich, Switzerland; Goddard Space Flight Center, USA; Institute d'Astrophysique Spatiale, France; Instituto Nacional de Técnica Aeroespacial, Spain; Institute for Astronomy, Edinburgh, UK; Jet Propulsion Laboratory, USA; Laboratoire d'Astrophysique de Marseille (LAM), France; Leiden University, Netherlands; Lockheed Advanced Technology Center (USA); NOVA OptIR group at Dwingeloo, Netherlands; Northrop Grumman, USA; Max-Planck Institut für Astronomie (MPIA), Heidelberg, Germany; Laboratoire d'Etudes Spatiales et d'Instrumentation en Astrophysique (LESIA), France; Paul Scherrer Institut, Switzerland; Raytheon Vision Systems, USA; RUAG Aerospace, Switzerland; Rutherford Appleton Laboratory (RAL Space), UK; Space Telescope Science Institute, USA; Toegepast- Natuurwetenschappelijk Onderzoek (TNO-TPD), Netherlands; UK Astronomy Technology Centre, UK; University College London, UK; University of Amsterdam, Netherlands; University of Arizona, USA; University of Cardiff, UK; University of Cologne, Germany; University of Ghent; University of Groningen, Netherlands; University of Leicester, UK; University of Leuven, Belgium; University of Stockholm, Sweden; Utah State University, USA.

K.I.C. and E.I. acknowledge funding from the Netherlands Research School for Astronomy (NOVA). K.I.C., R.N.C., and V.K. acknowledge funding from the Dutch Research Council (NWO) through the award of the Vici Grant VI.C.212.036. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. L.C. acknowledges financial support from Comunidad de Madrid under Atracción de Talento grant 2018-T2/TIC-11612. S.G. acknowledges the
support of the Cosmic Dawn Center of Excellence funded by
the Danish National Research Foundation under grant 140. G.
Ö., A.B., and J.M. acknowledge support from the Swedish
National Space Administration (NSSA). A.A.H. acknowledges
support from PID2021-124665NB-100 funded by the Spanish
Ministry of Science and Innovation and the State Agency of
Research MCIN/AEI/10.13039/501100011033. J.H. and D.L.
were supported by a VILLUM FONDEN Investigator grant to
J.H. (project number 16599).

J.A.M. and A.C.G acknowledge support by grant PIB2021-127718NB-100 by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe.”

PGP-G acknowledges support from Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through grant PGC2018-093499-B-100.

J.P.P. and T.V.T. acknowledge funding from the UK Science and Technology Facilities Council, and the UK Space Agency.
*Facilities: HST, JWST.*

*Software:* ASTROPY (Astropy Collaboration et al. 2018), LEPHARE (Arnouts & Ilbert 2011), NUMPy (Harris et al. 2020), PANDAS (pandas development team 2020) PHOTUTILS (Bradley et al. 2021), SCIPy (Virtanen et al. 2020) SOURCE EXTRACTOR (Bertin & Arnouts 1996), TOPCAT (Taylor 2005).

ORCID iDs

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

R. Azzollini https://orcid.org/0000-0002-0438-0886

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

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

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

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

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

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

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

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

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

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

E. van Dishoeck https://orcid.org/0000-0001-7591-1907

M. Güdel https://orcid.org/0000-0001-9818-0588

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

R. Navarro-Carrera https://orcid.org/0000-0001-6066-4624

V. Kokorev https://orcid.org/0000-0002-5588-9156

References

Anders, P., & Fritze-v., A. U. 2003, A&A, 401, 1063

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

Arnouts, S., & Ilbert, O. 2011, LEPHARE: Photometric Analysis for Redshift Estimate, Astrophysics Source Code Library, ascl:1108.009

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Atek, H., Furtak, L., Oesch, P., et al. 2022, MNRAS, 511, 4464

Babusiaux, C., Fabricius, C., Khanna, S., et al. 2023, A&A, 674, A32

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

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

Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112

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

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

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

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2021, astropy/photutils: v1.0.2, Zenodo, doi:10.5281/zenodo.4453725

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

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

Burgarella, D., Buat, V., Gruppioni, C., et al. 2013, A&A, 554, A70

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

Caputi, K. I., Caminha, G. B., Fujimoto, S., et al. 2021, ApJ, 908, 146

Caputi, K. I., Cirasuolo, M., Dunlop, J. S., et al. 2011, MNRAS, 413, 162

Caputi, K. I., Deshmukh, S., Ashby, M. L. N., et al. 2017, ApJ, 849, 45

Caruana, J., Bunker, A. J., Wilkins, S. M., et al. 2014, MNRAS, 443, 2831

Chabrier, G. 2003, PASP, 115, 763

Coe, D. 2015, Trilogy: FITS image conversion software, Astrophysics Source Code Library, ascl:1508.009

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

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

Endsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, in press

Erb, D. K., Steidel, C. C., Shapley, A. E., et al. 2006, ApJ, 647, 128

Faisst, A. L., Capak, P., Hsieh, B. C., et al. 2016, ApJ, 821, 122

Faisst, A. L., Capak, P. L., Emami, N., Tacchella, S., & Larson, K. L. 2019, ApJ, 884, 133

Flores Velázquez, J. A., Gurvich, A. B., Faucher-Giguère, C.-A., et al. 2021, MNRAS, 501, 4812

Fontana, A., Vanzella, E., Pentericci, L., et al. 2010, ApJL, 725, L205

Fumagalli, M., Patel, S. G., Franx, M., et al. 2012, ApJL, 757, L22

Galametz, A., Grazian, A., Fontana, A., et al. 2013, ApJS, 206, 10

Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 436, 2875

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hopkins, P. F., Kereš, D., Oňorbe, J., et al. 2014, MNRAS, 445, 581

Iani, E., Caputi, K. I., Rinaldi, P., & Kokorev, V. I. 2022, ApJL, 940, L24

Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6

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

Kajisawa, M., Ichikawa, T., Yamada, T., et al. 2010, ApJ, 723, 129

Kennicutt, R. C. 1998, ARA&A, 36, 189

Khostovan, A. A., Sobral, D., Mobasher, B., et al. 2016, MNRAS, 463, 2363

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

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

Lam, D., Bouwens, R. J., Labbé, I., et al. 2019, A&A, 627, A164

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

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

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Lilly, S. J., Le Fevre, O., Hammer, F., & Crampton, D. 1996, ApJL, 460, L1

Loiacono, F., Decarli, R., Gruppioni, C., et al. 2021, A&A, 646, A76

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

Madau, P., Ferguson, H. C., Dickinson, M. E., et al. 1996, MNRAS, 283, 1388

Mármol-Queraltó, E., McLure, R. J., Cullen, F., et al. 2016, MNRAS, 460, 3587

Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2023, ApJ, 950, 67

Morishita, T., & Stiavelli, M. 2023, ApJL, 946, L35

Nanayakkara, T., Brinchmann, J., Glazebrook, K., et al. 2020, ApJ, 889, 180

Ning, Y., Cai, Z., Jiang, L., et al. 2023, ApJL, 944, L1
Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43

Oesch, P., Bouwens, R., Brammer, G., et al. 2021, JWST Proposal, Cycle 1, #1895

Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, arXiv:2304.02026

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

Ono, Y., Ouchi, M., Mobasher, B., et al. 2012, ApJ, 744, 83

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei (2nd ed.; Sausalito, CA: University Science)

Pandas Development Team, T. 2020, pandas-dev/pandas: Pandas, vlatest, Zenodo, doi:10.5281/zenodo.3509134

Patel, S. G., Kelson, D. D., Abramson, L. E., Sattari, Z., & Lorenz, B. 2023, ApJ, 945, 93

Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Pentericci, L., Vanzella, E., Fontana, A., et al. 2014, ApJ, 793, 113

Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJL, 951, L1

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Prieto-Lyon, G., Strait, V., Mason, C. A., et al. 2023, A&A, 672, A186

Reddy, N. A., Shapley, A. E., Sanders, R. L., et al. 2018, ApJ, 869, 92

Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584

Rieke, M. J., Kelly, D., & Horner, S. 2005, Proc. SPIE, 5904, 1

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

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

Rosani, G., Caminha, G. B., Caputi, K. I., & Deshmukh, S. 2020, A&A, 633, A159

Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100

Salpeter, E. E. 1955, ApJ, 121, 161

Schaerer, D., & de Barros, S. 2009, A&A, 502, 423

Schiminovich, D., Ilbert, O., Arnouts, S., et al. 2005, ApJL, 619, L47

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Shim, H., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 738, 69

Smit, R., Bouwens, R. J., Labbé, I., et al. 2016, ApJ, 833, 254

Sobral, D., Best, P. N., Smail, I., et al. 2014, MNRAS, 437, 3516

Sonnett, S., Meech, K., Jedicke, R., et al. 2013, PASP, 125, 456

Sparre, M., Hayward, C. C., Feldmann, R., et al. 2017, MNRAS, 466, 88

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

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

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

Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129

Stasińska, G., & Leitherer, C. 1996, ApJS, 107, 661

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

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

Tang, M., Stark, D. P., Chevallard, J., & Charlot, S. 2019, MNRAS, 489, 2572

Taylor, M. B. 2005, in ASP Conf. Ser. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert (San Francisco, CA: ASP), 29

Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2023, ApJ, 945, 35

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

Wang, X., Cheng, C., Ge, J., et al. 2022, arXiv:2212.04476

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

Williams, C. C., Tacchella, S., Maseda, M. V., et al. 2023, arXiv:2301.09780

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

Wright, G. S., Wright, D., Goodson, G. B., et al. 2015, PASP, 127, 595