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

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 λ_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 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₀) 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₀ 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β λ4861 Å and [O iii] λ λ4959, 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 ≃4–5 μm. In turn, the (Hα λ6563 + [N ii] λ λ6548, 6583 + [S ii] λ λ6716, 6730) complex appears in the MIRI (Rieke et al. 2015; Wright et al. 2015) wavelength domain at observed >5 μm.

In this paper, we make use of publicly available NIRCam images in the Hubble eXtreme Deep Field (XDF) to search for (Hβ + [O iii]) emitters at z ≃ 7–8. In most of this field, we also benefit from ultradeep MIRI 5.6 μm imaging, which we analyze to search for the presence of Hα emission in the same galaxies. This is the first time that the Hα 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 ≃ 7–8. In Section 3 we explain our methodology to identify strong (Hβ + [O iii]) and Hα 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₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 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.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.

Zoom In Zoom Out Reset image size

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 ²⁴ (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. ²⁵

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," ²⁶ 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 003 pixel⁻¹, as the final pixel scale we adopted in this work. All the final images have been aligned to the Hubble Legacy Fields (HLF) catalog. ²⁷

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². By reaching a median depth of 29.15 mag (5σ, r = 015), 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). ²⁸ 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/WFC F435W, 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

2.3. SED Fitting

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 }_{\nu }^{2}$ between 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 ≃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β + [O iii]) and (Hα+[N ii]+[S ii]) emitters in the XDF at z ≃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 ≃ 7–8.

Among these sources, we searched for (Hβ + [O iii]) and Hα 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β + [O iii]). In turn, MIRI/F560W has been used to look at the flux excess produced by Hα.

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

$$\begin{eqnarray}&&{\mathrm{EW}}_{0}=\displaystyle \frac{{W}_{\mathrm{rec}}}{1+z}({10}^{(-0.4{\rm{\Delta }}\mathrm{mag})}-1),\end{eqnarray} \tag{ 1 }$$

where W_rec is the rectangular width of the filter containing the emission line in question, in our case (Hβ + [O iii]) or Hα, and Δmag is the difference between the observed magnitude in that filter and the synthetic magnitude ³⁴ from the SED fitting (i.e., Δmag = m_obs - m_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 ∣mag_obs(F460M) − mag_syn(F460M)∣ ≤ 2 × mag_err(F460M), where mag_obs and mag_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: Δmag = (mag_X − mag_cont), where mag_X represents the magnitude in one of the filters we chose to select (Hβ+[O iii]) or Hα, 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: Δmag < −0.2. Note that a Δmag = −0.2 in NIRCam/F430M corresponds to a EW₀ ≃58 Å at z = 7, while in NIRCam/F444W it would imply an EW₀ ≃270 Å. For MIRI/F560W, the same Δmag would correspond to an EW₀ ≃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 Δ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β + [O iii]) emitters. These emitters constitute ≃31% of our total galaxy sample at z ≃ 7–8 (see Figure 2 where we show the multiwavelength images of an example source). The derived EW₀ values cover a wide range that goes from a minimum of ${87.5}_{-27}^{+30}\,{\rm{\mathring{\rm{A}} }}$ to a maximum value of ${2140.4}_{-154}^{+970}\,{\rm{\mathring{\rm{A}} }}$, with a median $\left\langle {\mathrm{EW}}_{0}\right\rangle \simeq {943}_{-194}^{+737}\,{\rm{\mathring{\rm{A}} }}$ (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 ≃ 8 galaxy candidates. Out of the 18 (Hβ + [O iii]) emitters, 83% have a best-fit SED with subsolar (0.2 Z_⊙) metallicity and the remaining ≃17% with solar (Z_⊙) metallicity.

Zoom In Zoom Out Reset image size

Among the 18 (Hβ + [O iii]) emitters at z ≃ 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α+[N ii]+[S ii]) line complex at z ≃ 7–8. To obtain the net value of the Hα EW₀, we applied the correction recipes provided by Anders & Fritze-v. (2003), as follows: f(Hα)=0.63f(Hα + [N ii] + [S ii]) for a solar metallicity, and f(Hα)=0.81f(Hα + [N ii] + [S ii]) for a 0.2 Z_⊙ 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β + [O iii]) and Hα 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.

Zoom In Zoom Out Reset image size

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 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σ. 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.

Zoom In Zoom Out Reset image size

4.1. Emission Line EW Versus Stellar Mass and Age in Galaxies at z ≃ 7–8

Having calculated the (Hβ + [O iii]) EW₀ for the prominent line emitters, we can compare their best-fit SED properties with those of the other z ≃ 7–8 galaxies in our sample. In Figure 5, we show the derived (Hβ + [O iii]) EW₀ versus the best-fit age. From this plot, we can see that all except three of the (Hβ + [O iii]) emitters are characterized by young best-fit ages (≤10⁸ yr), which indicates that these objects may be in their first major star formation episode. The remaining three objects are older (>10⁸ 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.

Zoom In Zoom Out Reset image size

The gray triangles in Figure 5 refer to the EW upper limits that we estimated for all those galaxies at z ≃ 7–8 that do not have a significant flux excess in the NIRCam/F430M band. In contrast to the (Hβ + [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β + [O iii]) emitters at z ≃ 6.8–7 that have been selected over a wide sky area (2.7 deg² in total). Endsley et al. (2021) found this rare population of very strong (Hβ + [O iii]) emitters with an EW₀ >1200 Å. The fact that we find similarly high (Hβ + [O iii]) EWs₀ among faint galaxies in a much smaller area of the sky indicates that prominent (Hβ + [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β (only) EW₀ versus age for SB99 model galaxies. These theoretical tracks are based on a Chabrier IMF with a stellar mass cutoff of 100M_⊙ and were obtained both for a solar and a subsolar metallicity (0.2Z_⊙), 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.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&f({\rm{H}}\beta )=f({\rm{H}}\alpha )\times {10}^{-0.4\times 1.27E(B-V)}/2.86,\end{eqnarray} \tag{ 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.

Zoom In Zoom Out Reset image size

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.

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α):

$$\begin{eqnarray}&&\mathrm{SFR}\ ({M}_{\odot }\,{\mathrm{yr}}^{-1})=7.9\times {10}^{-42}\,{L}_{{\rm{H}}\alpha }(\mathrm{erg}\,{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 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^−7.60 yr⁻¹.

Zoom In Zoom Out Reset image size

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 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 (SFR_Hα and SFR_UV). 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 SFR_Hα and SFR_UV plane may be real and due to the following:

• 1.  
  In very young galaxies (below 100 Myr), SFR_UV underestimates the real value because the UV luminosity associated with star formation is still growing. Indeed, when comparing SFR_Hα and SFR_UV, one has to also take age effects into account. UV traces typically 1500–2000 Å (i.e., nonionizing photons), while Hα traces directly <912 Å photons. For example, UV-bright regions without Hα 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_Hα and SFR_UV 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_Hα /SFR_UV 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_Hα in excess compared to SFR_UV, particularly in low-mass galaxies. However, there are still uncertainties in determining the ratio of SFR_Hα /SFR_UV 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 ≃ 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 ≃ 7–8. To do that, we sum up the individual SFRs (SFR_Hα,total ≃51.39 M_⊙ yr⁻¹) and then divide the total by the comoving volume ³⁵ encompassed by the area ($A\simeq 4.7\,{\mathrm{arcmin}}^{2}$) and redshift bin (i.e., z ≃ 7–8) analyzed in this work (V_sky ≃ 11580.26 Mpc³). We obtain that, at these redshifts, the Hα emitters make for ${\mathrm{log}}_{10}({\rho }_{{\mathrm{SFR}}_{{\rm{H}}\alpha }}/({{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}))\,\simeq -2.35\pm 0.3$.

Zoom In Zoom Out Reset image size

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 ≃ 9 to z ≃ 15 as well as Pérez-González et al. (2023) results at z ≃ 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 ≃ 7–8 (${\mathrm{log}}_{10}({\rho }_{{\mathrm{SFR}}_{\mathrm{tot}}}/({M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}))\simeq -1.76\pm 0.3$). ³⁶ 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).

Zoom In Zoom Out Reset image size

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

Finally, our derived values of the Hα EW₀ allow us to extend the study of the redshift evolution of this parameter to z ≃ 7–8. In Figure 11 we present our results along with the most recent determinations from the literature (for sources at z ≃ 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 ≃ 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.

Zoom In Zoom Out Reset image size

Our sample of strong line emitters at z ≃ 7–8 allows us to populate a virtually unexplored part of parameter space. At those redshifts (z ≃ 8), only Stefanon et al. (2022) previously obtained an estimate of the average Hα EW₀, 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₀(Hα) should evolve differently below and above z ≃ 2. In particular, according to the recent literature, at z < 2, the EW₀(Hα) should evolve as ∝(1 + z)^1.8, while at z > 2 it should evolve as ∝(1 + z)^1.3.

By inspection of Figure 11, we can see that JWST observations at z ≥ 6 (i.e., Sun et al. 2022; Ning et al. 2023; this present work) suggest that the break proposed at z ≃ 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_0(Hα) as a function of redshift again by considering the recent JWST observations at z ≥ 6 as well. In this case, we find that EW₀ ${({\rm{H}}\alpha )\propto (1+{\rm{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 ≃ 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₀ at z ≃ 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.

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β+[O iii]) and Hα emitters at z ≃ 7–8. This is the first time the Hα 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β+[O iii]) emitters at z ≃ 7–8, as determined from their F430M and F444W flux excess. These 18 galaxies constitute ≃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α emission line. The (Hβ+[O iii]) EWs₀ that we derive range from $\simeq {87.5}_{-27}^{+30}\,{\rm{\mathring{\rm{A}} }}$ to ${2140.4}_{-154}^{+970}\,{\rm{\mathring{\rm{A}} }}$, with a median value of ${943}_{-194}^{+737}\,{\rm{\mathring{\rm{A}} }}$. For most of these galaxies, we find [O iii]/Hβ > 1, but a few have [O iii]/Hβ < 1. The two line fluxes can be separated by making use of the independent Hα emission line measurement. This is telling us that some of the prominent (Hβ+[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α emitters show an EW₀ that ranges from a few hundred to a few thousand Angstroms. Some of these values are substantially above the expected median Hα EW at these redshifts, as extrapolated from lower redshift determinations. We also report that, by considering the recent JWST findings at z ≥ 6 (including this present work), EW₀(Hα) as a function of the redshift should evolve as follows: EW₀(Hα) ∝ (1 + z)^2.1. However, larger samples of galaxies are needed to confirm this result. We note, however, that the prominent Hα emitters only constitute about a quarter of all the MIRI-detected galaxies at z ≃ 7–8. For the remaining galaxies, the Hα EW₀ should lie below the expected median trend. As the Hα 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 z ≃ 7–8, which are relatively evolved galaxies, with best-fit ages >10⁷–10⁸ yr and stellar masses >10⁸ M_⊙.

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⁷ 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⁸ M_⊙), 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 ≃ 7–8. We found ${\mathrm{log}}_{10}({\rho }_{{\mathrm{SFR}}_{{\rm{H}}\alpha }}/({M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}))$ ≃ −2.35 ± 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 ≃ 7–8. We also considered the SFR_UV for all the other galaxies at z ≃ 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 ≃ 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.

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 expertize 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 Opt-IR 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 (SNSA). A.A.H. acknowledges support from PID2021-124665NB-I00 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-I00.

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

Facilities: HST - Hubble Space Telescope satellite, JWST - James Webb Space Telescope.

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).