The ionizing photon production efficiency at z ∼ 6 for Lyman-alpha emitters using JEMS and MUSE Open Access

ABSTRACT

We study the ionizing photon production efficiency at the end of the Epoch of Reionization (z ∼ 5.4 − 6.6) for a sample of 30 Ly α emitters. This is a crucial quantity to infer the ionizing photon budget of the universe. These objects were selected to have reliable spectroscopic redshifts, assigned based on the profile of their Ly α emission line, detected in the MUSE deep fields. We exploit medium-band observations from the JWST Extragalactic Medium-band Survey (JEMS) to find the flux excess corresponding to the redshifted Hα emission line. We estimate the ultraviolet (UV) luminosity by fitting the full JEMS photometry, along with several HST photometric points, with Prospector. We find a median UV continuum slope of |$\beta = -2.09^{+0.23}_{-0.21}$|⁠, indicating young stellar populations with little-to-no dust attenuation. Supported by this, we derive ξ_ion,0 with no dust attenuation and find a median value of log|$\frac{\xi _{ion,0}}{\text{Hz erg}^{-1}} = 25.44^{+0.21}_{-0.15}$|⁠. If we perform dust attenuation corrections and assume a Calzetti attenuation law, our values are lowered by ∼0.1 dex. Our results suggest Ly α emitters at the Epoch of Reionization have slightly enhanced ξ_ion,0 compared to previous estimations from literature, in particular, when compared to the non-Ly α emitting population. This initial study provides a promising outlook on the characterization of ionizing photon production in the early universe. In the future, a more extensive study will be performed on the entire data set provided by the JWST Advanced Deep Extragalactic Survey (JADES). Thus, for the first time, allowing us to place constraints on the wider galaxy populations driving reionization.

1 INTRODUCTION

The Epoch of Reionization (EoR) is one of the major phase changes of the universe. It corresponds to the transition between a dark and neutral universe to an ionized one, where the intergalactic medium (IGM) became transparent to Lyman continuum (LyC) radiation. Observational evidence places the end of this epoch at redshift z ∼ 6 (Becker et al. 2001; Fan et al. 2006; Yang et al. 2020). Understanding the sources responsible for ionizing the IGM is one of the most significant, unsolved questions in modern-day astronomy. Young massive stars in galaxies are currently believed to be the main responsible culprit, producing copious amounts of hydrogen-ionizing photons (E ≥ 13.6 eV) that escape the interstellar medium (ISM) to then ionize the IGM (Hassan et al. 2018; Rosdahl et al. 2018; Trebitsch, Volonteri & Dubois 2020). Until recently, average escape fractions (f_esc) of 10 – 20 per cent were believed to be a key threshold for galaxies to be the main responsible sources of reionization (Ouchi et al. 2009; Robertson et al. 2013, 2015; Finkelstein et al. 2019; Naidu et al. 2020). High escape fractions have been seen in some LyC leaking galaxies (e.g. Borthakur et al. 2014; Bian et al. 2017; Vanzella et al. 2018; Fletcher et al. 2019; Ji et al. 2020; Izotov et al. 2021). These are galaxies from which we observe (or infer through indirect probes) LyC radiation. Understanding how LyC photons escape into the IGM and subsequently ionize it during the EoR is of utmost importance. Strong LyC leakers are not generally observed in large samples (Leitet et al. 2013; Leitherer et al. 2016; Steidel et al. 2018; Flury et al. 2022a). Based on population-wide studies, it is currently uncertain which types of galaxies dominate the budget of reionization (Finkelstein et al. 2019; Naidu et al. 2020). Promisingly, studies targeting galaxies at high redshift (up to z ∼ 9) have found that as we look further into the past of the universe, galaxies seem to be more efficient in producing ionizing photons. This can be measured as the ratio between the production rate of ionizing photons over the non-ionizing luminosity density, called the ionizing photon production efficiency (ξ_ion), which has been shown to increase with redshift (e.g. Bouwens et al. 2016; Faisst et al. 2019; Endsley et al. 2021; Stefanon et al. 2022; Tang et al. 2023). An increase of ξ_ion implies that smaller f_esc are required in galaxies to explain the reionization of the universe.

Direct observations of LyC radiation are impossible at the EoR due to the increasing absorption by neutral hydrogen in the IGM along the line of sight. At z ∼ 6 the average IGM transmission of photons with λ_rest−frame = 900 Å is virtually zero (Inoue et al. 2014). However, hydrogen recombination lines can provide indirect evidence of ionizing photons, since they are produced after photoionization has taken place. The strongest of such lines is Ly α, and due to its resonant nature, the shape of its profile is to date one of the most reliable indirect tracers of LyC leakage. In particular, when double peaks are detected, the separation between them can be related to the neutral hydrogen column density and consequently, to the escape of LyC photons (Verhamme et al. 2015, 2017; Hogarth et al. 2020; Izotov et al. 2021; Naidu et al. 2022). Moreover, LyC leakers tend to show strong Ly α emission (Nakajima et al. 2020). These kind of objects have high Ly α equivalent widths (EW(Ly α) > 20 Å; Ajiki et al. 2003; Ouchi et al. 2005, 2008), in some extreme cases reaching hundreds of angstroms (Kerutt et al. 2022; Saxena et al. 2023), and are appropriately called Ly α emitters (LAEs) (see Ouchi, Ono & Shibuya 2020, and references within). We note that the Ly α–LyC correlation has scatter, and that strong LyC leakage has also been observed in galaxies with relatively low Ly α emission (e.g. Ji et al. 2020; Flury et al. 2022b).

The second strongest hydrogen recombination line is Hα. Unlike Ly α, it does not scatter resonantly under normal star-forming conditions. It is commonly used as a measure of star formation rate (SFR), because its luminosity is directly related to the amount of ionizing photons that are scattered or absorbed by the medium, and thus, do not escape the ISM. If used in combination with a measure of the ultraviolet (UV) continuum, Hα can constrain ξ_ion (e.g. Bouwens et al. 2016; Chisholm et al. 2022; Stefanon et al. 2022). We note that ξ_ion has a strong degeneracy with LyC escape fractions, thus, it is common to report the values obtained under the assumption of zero LyC escape (ξ_ion,0).

In this work we aim to estimate ξ_ion for a sample of bright LAEs identified with the ground-based Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010), by measuring Hα and the UV continuum through observations obtained with the Near-Infrared Camera (NIRCam; Rieke, Kelly & Horner 2005) onboard the JWST (Gardner et al. 2006). In this first study we focus on a sample of bright LAEs, the advantage of using confirmed LAEs is the availability of reliable systemic redshifts, which can be derived from their Ly α profiles. In addition, LAEs could be the main producers of ionizing photons at the EoR, due to their connection to LyC leakage (Gazagnes et al. 2020; Nakajima et al. 2020; Izotov et al. 2022; Naidu et al. 2022; Matthee et al. 2022b). In the future, this study will be expanded upon with larger samples using the JWST Advanced Deep Extragalactic Survey (JADES; Eisenstein et al. 2017; Rieke 2019, Eisenstein (in preparation)) and a stellar-mass selected sample, to understand how ionizing radiation varies among all known galaxy types at the EoR.

The structure of this paper is the following. In Section 2 we describe the data sets and catalogues used in this work, followed by the methods used to measure Hα and the UV continuum luminosity in Section 3. In Section 4 we present our Prospector fitting method, along with the dust attenuation prescription used in this work. Using these results, we place constraints on ξ_ion in Section 5, and discuss them in Section 6. Finally, we provide concluding remarks in Section 7. Throughout this work we assume Ω₀ = 0.315 and H₀ = 67.4 km s⁻¹ Mpc⁻¹, following Planck Collaboration et al. (2020).

2 DATA

2.1 JEMS

The JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023, PID = 1963) is a JWST imaging program whose primary goal is to target Hα and the UV continuum at the EoR (z ∼ 5.4 – 6.6). Many galaxies have been identified in this redshift range through the Lyman break technique, using the original HST/ACS images (e.g. Bunker et al. 2004; Bouwens et al. 2015). Including several with spectroscopic confirmation (e.g. Bunker et al. 2003; Stanway et al. 2004). JEMS covers the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006) region in the GOODS-South field, which has been extensively studied with (among others) deep HST and MUSE observations. The survey consists of single visit observations with three filter pairs on the NIRCam: F210M-F430M, F210M-F460M, and F182M-F480M. In each pair, one module covers the entirety of the UDF while the second one points to the surrounding region, both pointings have publicly available ancillary data. Fig. 1 shows the throughputs of each filter as a function of wavelength, along with the wavelength evolution of the observed Ly α and Hα with redshift (dashed and full black lines, respectively) in the redshift range of interest. The spectral coverage of MUSE is shown as a grey shaded area. The effective wavelength, exposure time and 5σ sensitivity limit for each band is given in the caption.

The source detection and photometry leverage both the JEMS NIRCam medium band and JADES (Eisenstein et al., in preparation) NIRCam broad and medium band imaging. Mosaics of the NIRCam long wavelength filter images (F277W, F335M, F356W, F410M, F430M, F444W, F460M, and F480M) are combined into an inverse variance-weighted signal-to-noise ratio (SNR) image. Detection is performed using the photutils (Bradley et al. 2022) software package, identifying sources with contiguous regions of the SNR mosaic with signal >3σ and five or more contiguous pixels. We also use photutils to perform circular aperture photometry with filter-dependent aperture corrections based on empirical point-spread-functions measured from stars in the mosaic. The object flux within a range of aperture sizes is measured, ranging from r = 0.1″ to r = 0.5″ as well as the PSF 80 per cent encircled energy radius in each filter. Kron (1980) photometry is also performed using Kron parameters of 1.2 and 2.5. Details of the JADES source catalogue generation and photometry will be presented in Robertson et al., (in preparation).

2.2 LAEs from the MUSE ultra-deep field surveys

The MUSE is a large integral field unit on the VLT, in Chile, with a wavelength coverage of λ = 4800 – 9300 Å. The MUSE Hubble Ultra-Deep Field surveys (MUSE UDF) consist of deep observations in the HUDF area, in this work we use these surveys to find LAEs at the tail end of the EoR. The first release (Bacon et al. 2017; Inami et al. 2017) was based on two data sets: (1) a 3 × 3 arcmin² mosaic of 9 MUSE fields with an exposure of 10 h, called ‘MOSAIC’, and (2) a single 1 × 1 arcmin² field with a depth of 31 h, called ‘UDF-10’. The second release (Bacon et al. 2023) extended this survey to unprecedented depths by including an adaptive optics assisted survey in the same area, reaching 141 h of exposure time. This new survey is called the MUSE eXtremely Deep-Field (MXDF) and is the deepest spectroscopic survey in existence. The astrometry was matched to the Hubble ACS astrometry, however, there is an offset between the HST and Gaia DR2 catalogues (Dunlop et al. 2017; Franco et al. 2018) amounting to ΔRA = +0.094 ± 0.042 arcsec and ΔDEC = −0.26 ± 0.10 arcsec. The average offset for each data set is given in Bacon et al. (2023), we take this offset into account when assigning counterparts to the LAEs in JEMS.

To assemble our sample, we use the catalogues presented in Bacon et al. (2023), in which the data sets from the first release were reprocessed with the same tools as the MXDF. The redshifts assignment and confidence are of particular interest for this work. Briefly, at z > 2.9 the redshift is estimated through the Ly α profile. The redshift is based on the peak of the Ly α line, known to be offset from the systemic redshift due to its resonant nature (Shapley et al. 2003; Song et al. 2014). In order to estimate the systemic redshift, statistical corrections from Verhamme et al. (2018) are applied in the catalogue. In summary, two corrections are taken into account: one based on the separation of the peaks when double peaks are detected, and one based on the full-width-half-maximum (FWHM) of the line. For the redshifts range of interest, the confidence of the detected line being Ly α (ZCONF) is a number between zero and three, where

• ZCONF = 0. No redshift solution can be found

• ZCONF = 1. Low confidence arising for example from low S/N of the lines or the existence of other redshift solutions

• ZCONF = 2. Good confidence. The Ly α line is detected with a S/N > 5, and has a width and asymmetry that are compatible with Ly α line profiles

• ZCONF = 3. High confidence. If Ly α is the only detected line, then it has to have a S/N > 7 with the expected line shape (i.e. a red asymmetrical profile and/or a blue bump, or double peaked profile)

A group of experts iterates until they converge on a classification for each object. Often the difference between ‘good’ and ‘high’ confidence are subtle and objects in either group can be considered as having a certain redshift assignment.

2.3 Sample

We focus on LAEs at z ∼ 5.4 − 6.6 that have ZCONF = 2 or ZCONF = 3, and for which Hα falls in one of the following JEMS filters: F430M, F460M, or F480M. This leaves us with a sample of 30 LAEs with spectroscopic redshift, shown as horizontal line markers on Fig. 1. After correcting for the astrometry offset we find counterparts in JEMS by searching within a 0.5″ radius of the coordinates, this radius is favourable due to the known possible offset between the Ly α and the stellar continuum emission (Hoag et al. 2019; Claeyssens et al. 2022). We then visually examine each source to confirm the counterpart assignment. Table 1 shows the coordinates of JEMS survey counterparts, along with their systemic redshift, Ly α flux, and JEMS photometry. We note that we see a flux excess that corresponds to the redshifted Hα location in these objects. All the Ly α profiles and medium band cutouts are shown in the Appendix  A. For reference, seven of these galaxies were identified previously as Lyman Break Galaxies (LBGs) in Bunker et al. (2004), and four of them are part of the Hα emitting sample studied in Helton et al. (2023) they have been highlighted in Table 1.

3 Hα MEASUREMENT

There has been extensive work in the literature to retrieve Hα fluxes by combining multiband photometry to estimate both the line flux and the local continuum (e.g. Bunker et al. 1995; Shim et al. 2011; Labbé et al. 2013; Stark et al. 2013). It is a common approach to estimate Hα emission from photometry, and then use models to interpret the results (e.g. Sobral et al. 2012; Vilella-Rojo et al. 2015; Faisst et al. 2016). The JEMS survey has proven to be useful to measure emission lines in high redshift galaxies (7 < z < 9; Laporte et al. 2023), in this work we use JEMS photometry to measure Hα luminosities. Specifically, depending on the redshift of the sources, we interpret colour excess in the F430M, F460M, or F480M bands as the presence of Hα emission. The luminosity can then be estimated by using the remaining two filters to measure the local continuum and subtract it. We first study the possible contamination of [N ii] in the photometry, then describe the steps that were followed to measure Hα, and correct it for dust attenuation.

3.1 [N ii] contribution to Hα photometry

Extracting Hα emission from broad-band photometry can be troublesome, mainly due to the nearby [N ii] emission lines at rest-frame 6548 and 6583 Å, although also potentially affected by other metal lines in the vicinity of the continuum band measurement. To study the effect this contamination might have in our medium band observations, we create appropriate photoionization models and simulate photometry for F182M, F210M, F430M, F460M, and F480M, for the redshifts of the LAEs in the sample. It is important to note that at the redshifts of interest for this work (i.e. z ∼ 5.4 − 6.6) we expect metal-poor galaxies with high ionization parameters (log〈U〉) to dominate (Harikane et al. 2020; Katz et al. 2022; Sugahara et al. 2022). Although recently some tentative evidence has been shown for a galaxy with solar metallicity at z > 7 (Killi et al. 2023), this is not the norm. Promisingly, Cameron et al. (2023) produce emission line diagnostic diagrams for 26 LBGs at z ∼ 5.5 − 9.5 using spectroscopy obtained with the JWST/NIRSpec micro-shutter assembly (Ferruit et al. 2022; Jakobsen et al. 2022), and find that these galaxies populate a parameter region comparable to extreme local ones, having low metallicity and high ionization parameters (Curti et al. 2023; Nakajima et al. 2023).

To investigate the relative intensity of [N ii] compared to Hα we produce, as in Simmonds, Schaerer & Verhamme (2021), photoionization models following a sample of 5607 star-forming galaxies (Izotov et al. 2014) from the Sloan Digital Sky Survey (SDSS) Data Release 14 (Abolfathi et al. 2018), compiled by Y. Izotov and collaborators. The properties of these galaxies are discussed in previous literature (Guseva et al. 2020; Ramambason et al. 2020) but of main importance to this work are their high [O iii]/Hβ ratios, and their low metallicities (mean metallicity, 12+log(O/H) = 7.97|$^{+0.14}_{-0.17}$|⁠).

We then use the version 17.03 of the photoionization code cloudy (Ferland et al. 2017), combined with stellar population models from the Binary Population and Spectral Synthesis version 2.2.1 (BPASS; Eldridge et al. 2017). We run a grid of simple cloudy photoionization models to convergence, using BPASSv2.2.1 models with binary interactions. The initial mass function (IMF) only slightly changes the results, therefore, we choose to use a Salpeter IMF (α = −2.35, with M_*,max = 100 M_⊙). At redshifts close to the EoR, we expect galaxies to have young and metal-poor stellar populations, so we decide to use 10 Myr old stellar population models with a constant star formation history (SFH). We use two different metallicities, one that follows the bulk of the the SDSS sample (Z = 0.006; 12+log(O/H) = 8.1), and one to explore lower metallicity galaxies (Z = 0.001; 12+log(O/H) = 7.3), shown as upwards and downwards triangles in Fig. 2, respectively. Most importantly, we vary the ionization parameter (log〈U〉), in a range between −3.5 and −0.5. This is a dimensionless ratio of ionizing over non-ionizing hydrogen photon densities. A spherical cloud with the same gas-phase metallicity as the stellar SEDs is assumed in each case. The left panel of Fig. 2 shows one of the classical low-ionization emission line diagrams (BPT diagram; Baldwin, Phillips & Terlevich 1981; Kewley et al. 2006). The star-forming SDSS sample is plotted as grey dots. These galaxies populate virtually the same region of the BPT diagram that is expected of galaxies at z > 6, that is, log〈U〉 > −3.0 and log [O iii]/Hβ >0.4 (see figs 5 and 6 of Sugahara et al. 2022), which has recently been confirmed by observations of z ∼ 6 LBGs in Cameron et al. (2023), making them appropriate low-redshift analogues to galaxies at the EoR.

The relation between the [N ii]/Hα intensities ratio and log〈U〉 is shown for our cloudy models in the middle panel of Fig. 2, and are described by the following cubic equations

and

where the subscripts represent the intensity ratios for the cloudy models corresponding to metallicities Z = 0.001 (dashed curve, 12+log(O/H) = 7.3) and Z = 0.006 (filled curve, 12+log(O/H) = 8.1), respectively, and ‘x’ represents log〈U〉.

We find that the Hα emission dominates over [N ii] in our photoionization models, in agreement with previous studies (Maiolino & Mannucci 2019; Onodera et al. 2020). In addition, [N ii] appears to be virtually non-existent in a large sample of galaxies at z = 5.1 − 5.5, for which grism spectra covering the Hα emission line is available (Sun in preparation). At its maximum, our models indicate [N ii] can account for ∼45 per cent of the band flux, but this is only in cases where log〈U〉 = –3.5. At z ∼ 6 higher ionization parameters are predicted (e.g. Sugahara et al. 2022). Moreover, log〈U〉 is directly related to the [O iii]/[O ii] ratio, which has been observed to be elevated at z > 6 (Cameron et al. 2023), and therefore, we assume negligible [N ii] contamination in our next steps.

3.2 Hα luminosity estimation

We define four redshift bins in our sample, depending on the expected wavelength of Hα, and for each bin use a slightly different prescription to estimate Hα luminosities. We find that the local continuum is more reliably measured when we use the full JEMS photometry to fit a line in logarithmic space, and use this fit to estimate the local continuum around Hα. We note that for our objects there is no evidence of any strong emission lines (other than Hα) contaminating the JEMS photometry, and thus, for the continuum estimation we only exclude the band that contains Hα in each bin, as follows

• z ≤ 5.75: f(Hα) falls in F430M

• 5.75 < z ≤ 6: f(Hα) falls in F460M

• 6 < z ≤ 6.15: f(Hα) falls in both F460M and F480M

• z > 6.15: f(Hα) falls in F480M

The JEMS photometry flux densities in units of nJy can be found in Table 1, while the derived Hα luminosities are presented in Table 2. We have enforced a floor error of 5 per cent in the Hα measurements. Our measurements (when available) are consistent with those reported in Helton et al. (2023). It must be noted, however, that this simple prescription in some cases underestimates the flux of Hα emission. In particular, there are two LAEs (JADES-GS+53.17829-27.77728 and JADES-GS 53.18744-27.77804) in which the F430M photometry yields negative fluxes. In these cases we ignore this filter when estimating Hα, resulting in lower limits for ξ_ion,0 (circles with black edges in Fig. 5). In the case of JADES-GS 53.16674-27.80424, the Hα emission line falls in the wavelength range with low transmission of the F460M filter. The Hα line flux derived from the NIRCam grism observations through the First Reionization Epoch Spectroscopic Complete Survey (FRESCO, PI: Oesch; Oesch et al. 2023) is approximately six times higher than our photometric measurement (private communication from F. Sun). For consistency, we estimate ξ_ion,0 using the underestimated Hα for this source (shown as a circle with a black edge and a plus sign in Fig. 5), but note that this value is a lower limit.

3.3 The Ly α/Hα ratio

The synergy between MUSE and NIRCam is exceptional to study the Ly α/Hα ratio, since there are no slit loss effects and MUSE is currently the best instrument to observe Ly α at the redshift range of this work. Fig. 3 compares the estimated Hα to the measured Ly α luminosities. If Case B recombination is assumed, Ly α/Hα =8.7 (filled grey line, Osterbrock & Ferland 2006). In LAEs this ratio has been observed to fluctuate between ∼0.1 − 10, the resonant nature of the Ly α emission line complicates the interpretation of this fluctuation, but it has been proposed that dust attenuation plays a big role by reducing the Ly α emission (see Scarlata et al. 2009; Cowie, Barger & Hu 2011; Finkelstein et al. 2011; Hayes et al. 2013, 2014). Therefore, our range of Hα fluxes are expected and in agreement with previous studies.

If we assume Case-B recombination, the ratio between the Ly α and Hα luminosities can give us constraints on the Ly α escape fractions, using the simple relation given by atomic physics: f_esc(Ly α)|$= \frac{1}{8.7} \frac{\text{L(Ly$\alpha $)}}{\text{L(H$\alpha $)}}$|⁠. We provide these estimations, excluding the cases in which Hα is a lower limit, in Table 2.

4 UV LUMINOSITIES AND DUST ATTENUATION

4.1 SED fitting with Prospector

To obtain the UV continuum slope, β, we use the galaxy spectral energy distribution (SED) fitting code Prospector (Johnson et al. 2019, 2021). This tool infers stellar population parameters from the UV to IR wavelengths using photometry and/or spectroscopy as inputs. For this work we use the JEMS photometry, in addition to observations in the HST ACS bands: F435W (λ_eff = 0.432 µm), F606W (λ_eff = 0.578 µm), F775W (λ_eff = 0.762 µm), F815W (λ_eff = 0.803 µm), F850LP (λ_eff = 0.912 µm), and the HST WFC3 IR bands: F105W (λ_eff = 1.055 µm), F125W (λ_eff = 1.249 µm), F140W (λ_eff = 1.382 µm), and F160W (λ_eff = 1.537 µm). The same circular aperture of radius 0.15″ is used to extract both the JEMS and HST photometry, all photometry has been aperture corrected. We fix the redshifts to the ones provided in Bacon et al. (2023), then allow the dust attenuation and properties of the stellar populations (stellar and gas-phase metallicity, mass, ionization parameter, and SFH) to vary following the prescription in Tacchella et al. (2022a). Specifically, for the dust attenuation we use a two component dust model that accounts for the increased effect on young stars, embedded in dust, and the nebular emission (<10 Myr, as described in Conroy, Gunn & White 2009), with a variable dust index. We allow the stellar metallicity to vary between log(Z/Z_⊙) = −3.0 and 0.19, adopting a top-hat prior in logarithmic space. The continuum and line emission properties of the SEDs are provided by the Flexible Stellar Population Synthesis (FSPS) code (Byler et al. 2017), which are based on interpolated models from cloudy (version 13.03; Ferland et al. 2013), limiting the ionization parameter to a maximum value of –1. Regarding SFHs, we adopt a non-parametric SFH (continuity SFH; Leja et al. 2019). In brief, this option describes the SFH as a combination of six different SFR bins with the bursty-continuity prior (Tacchella et al. 2022b). The latter consists of five free parameters controlling the ratios of the amplitudes of the adjacent bins. Finally, we include both nebular and dust emission. The nebular emission is particularly important because we incorporate the JEMS photometry that contain the Hα emission line in our fits.

Once the best fit spectrum has been produced for each galaxy, we fit a straight line in logarithmic space between rest-frame λ = 1250 − 2600 Å (Calzetti, Kinney & Storchi-Bergmann 1994), in the form f_λ ∝ λ^β. The results are given in Table 2. In addition to providing constraints on the dust attenuation, this fitted line also allows us to predict the monochromatic UV luminosity density at rest-frame λ = 1500 Å. A representative example is shown on Fig. 4.

4.2 Dust attenuation estimation

The observed Hα luminosities must now be corrected for dust attenuation. Since we do not have access to Hβ, and therefore, cannot measure the Balmer decrement, we can use the rest-frame UV continuum slope (β; Calzetti et al. 1994) obtained with Prospector to estimate the effect of dust attenuation in the galaxies (e.g. Reddy et al. 2018). We note that the values obtained following this method are consistent with the dust attenuation obtained with Prospector.

The stellar and nebular continuum colour excess, E(B – V), can be estimated using

given in Reddy et al. (2018). The nebular excess is then simply given by E(B − V)_neb = 2.27 × E(B − V)_stellar (Calzetti et al. 2000). The optical depth of Hα is given as a function of E(B − V)_neb in Reddy et al. (2020) as τ ∼ 3 × E(B − V)_neb assuming a Calzetti attenuation law. The attenuation law at high redshifts is not well understood (Gallerani et al. 2010; Ma et al. 2019). Moreover, the amount of dust attenuation that affects emission lines versus the stellar continuum is highly uncertain. At z ∼ 2, nebular lines have been observed to be attenuated by the same amount as the stellar continuum (Calzetti et al. 2000; Erb et al. 2006; Reddy et al. 2010, 2015). In addition, the lack of understanding of the geometry and amount of dust attenuation at early galaxy phases complicate the situation even further (Bowler et al. 2018, 2022). It has been shown that a steeper attenuation curve like SMC is more appropriate than a Calzetti law for young high-redshift galaxies (Shivaei et al. 2020). The choice to use Calzetti over SMC in this work is due the increased effect it has in lowering ξ_ion,0 (e.g. log ξ_ion (SMC) ∼ log ξ_ion (Calzetti)+0.3; Shivaei et al. 2018). And to thus highlight the increased ionizing photon production we find in our sample.

Given that the galaxies in this work show evidence of blue UV slopes (β ∼ −2.2) we believe the dust attenuation effects to be small. We therefore present our results uncorrected by dust in Table 2. However, for illustrative purposes, in Fig. 5 we show both uncorrected (circles) and dust-corrected (crosses) values, assuming a Calzetti attenuation law, and adopting local relations for the nebular and stellar continuum attenuation. This prescription lowers our resulting ξ_ion,0 by ∼0.1 dex in the few cases where β suggests non-negligible dust attenuation. We note that if an SMC law is assumed instead, the effect is even smaller.

5 CONSTRAINTS ON ξ_ion

The dust-corrected Hα luminosity is related to the amount of ionizing photons (N(H⁰)) that are being emitted. If we assume Case B recombination (f_esc = 0), a temperature of 10⁴ K and an electron density of 100 cm⁻³, this relation is described in Osterbrock & Ferland (2006) as

where N(H⁰) is in units of Hz and L(Hα) in units of erg s⁻¹. Since we expect our sources to have a non-zero LyC escape fraction, the Case B recombination assumption yields a conservative estimation of the amount of ionizing photons. If these galaxies are LyC leakers, which we will study in a future work, the ionizing photon production efficiency would be boosted with respect to the values presented here. We use this equation to estimate the production rate of ionizing photons, ξ_ion,0, using

where the zero subscript indicates the assumption of zero escape of ionizing photons into the IGM, L_UV is the observed monochromatic luminosity density in units of erg s⁻¹ Hz⁻¹ at rest-frame λ = 1500 Å.

Given the uncertainty in the effects of dust on the Hα emission line and the UV continuum, in addition to the blue continuum slopes found in our sample (median values of |$\beta = -2.09^{+0.23}_{-0.21}$|⁠, corresponding to |$L_{UV} = 1.58^{+1.09}_{-0.60} \times 10^{43}$| erg s⁻¹), in this work we present the ξ_ion,0 values uncorrected for dust in Table 2. We show both corrected (crosses) and uncorrected (circles) ξ_ion,0 estimations in Fig. 5, and note that if a Calzetti law is applied, ξ_ion,0 is reduced by ∼0.1 dex for the few cases where dust is non-negligible.

We find that the median value of ξ_ion,0 (assuming negligible dust correction for Hα and the UV continuum) for our sample is log|$\frac{\xi _{ion,0}}{\text{Hz erg}^{-1}} = 25.44^{+0.21}_{-0.15}$|⁠. The individual estimations can be found in Table 2. This value is in broad agreement with those found for other LAEs at similar redshifts (e.g. Ning et al. 2023; Prieto-Lyon et al. 2023). A compilation of results from the literature is given in Stefanon et al. (2022), in particular, they provide a linear fit (in logarithmic space) of the ξ_ion observations up to z = 9 from Stark et al. (2015, 2017), Bouwens et al. (2016), Mármol-Queraltó et al. (2016), Nakajima et al. (2016), Matthee et al. (2017), Harikane et al. (2018), Shivaei et al. (2018), De Barros et al. (2019), Faisst et al. (2019), Lam et al. (2019), Tang et al. (2019), Emami et al. (2020), Naidu et al. (2020), Nanayakkara et al. (2020), Endsley et al. (2021), Atek et al. (2022). This extensive compilation illustrates a clear trend of increasing ξ_ion,0 with redshift (see fig. 4 of Stefanon et al. 2022). As shown in Fig. 5, our estimations scatter around this relation. It must be noted that this relation is not specific to LAEs, and is mostly tracing the ionizing photon production of LBGs. Thus, the increased ionizing photon production in a few of our sources is likely linked to the population being studied. Matthee et al. (2022a) investigate this further by comparing ξ_ion for different galaxy populations, finding non-LAEs tend to have lower ξ_ion than LAEs. The ionizing photon production efficiency has not been extensively studied in the redshift range of this work, for comparison, we show the median value obtained for 22 massive galaxies at z ∼ 7 from (black diamond; Endsley et al. 2021), and the 7 LAEs from (stars; Ning et al. 2023).

Fig. 6 shows a comparison of ξ_{ion, 0} and Ly α escape fractions, excluding the cases where the measured Hα flux is a lower limit (i.e. JADES-GS+53.16674-27.80424, JADES-GS+53.17829-27.77728, and JADES-GS 53.18744-27.77804). We observe a large variation of escape fractions for a given ξ_ion,0, with a tentative anticorrelation between f_esc(Ly α) and ξ_ion,0. The f_esc(Ly α) values are in agreement with other studies at similar redshifts (Jung et al. 2023; Ning et al. 2023; Tang et al. 2023, Saxena in preparation).

6 DISCUSSION

Studying ξ_ion,0 at the EoR is interesting to understand the escape fractions required for galaxies to be the main sources responsible for the reionization of the universe. The redshift range investigated here (z ∼ 5.4 − 6.6) has not been extensively studied in the past. This study has been made possible due to the deep imaging provided by NIRCam onboard JWST, in combination with the ground-based multi-unit spectrograph MUSE. The latter provided reliable systemic redshifts, that were used to find Hα through flux excesses in the JEMS bands F430M, F460M, and F480M.

One of the caveats in this work arises from the Hα measurements. As discussed in Section 3.1, we expect these galaxies to be metal poor, and thus, that contamination in the JEMS filters from other emission lines (especially [N ii]) is negligible. This might not be the case for all the objects in our sample, for example, Killi et al. (2023) suggest evidence of a solar metallicity galaxy at z ∼ 7. Thus, it is possible that some of the Hα luminosities here provided are upper limits. The Hα luminosities directly relate to the ionizing photon production through equation (4), where the conversion between Hα luminosity and ionizing photons produced (N(H⁰)) assumes Case B recombination. In other words, an escape fraction of zero since all the ionizing photons produced are being used to power the nebular emission. This assumption is adopted for simplicity, but it has to be noted that we expect non-zero escape fractions for these objects.

A second important point to keep in mind is the dust corrections applied to the Hα emission line and the UV stellar continuum at rest-frame λ = 1500 Å. Attenuation laws are highly uncertain at high redshift, as well as how much their relative contribution is to nebular lines versus to continuum attenuation. We show both our corrected and uncorrected results based on this uncertainty, we note that the blue slopes found with Prospector indicate young objects with little-to-no dust attenuation (β ∼ −2.2). In studies such as Bouwens et al. (2016), galaxies with UV slopes of β < −2.23 are left uncorrected for dust (A_UV = 0 according to a Calzetti et al. 2000, law). This is a reasonable approach since a steep β slope indicates a young stellar population with low dust attenuation.

We find slightly enhanced ionizing photon production efficiency values for some of the LAEs in our sample, in comparison to the average ξ_ion,0 at the redshift range studied. We again note that after correction for dust attenuation, assuming a Calzetti law and local relations between nebular and continuum attenuation, our values are lowered by ∼0.1 dex in the few cases where β suggests non-negligible dust attenuation. However, even after dust correction some of these objects remain above average for their redshift. This is not unheard of in LAEs, for example, Maseda et al. (2020) find enhanced ionizing photon production in a sample of 30 LAEs at z ∼ 4 − 5, with high Ly α EW (median value of 250 Å). Moreover, enhanced ξ_ion,0 values have been observed in LyC leakers at z ∼ 0.3, thought to be local analogues to high-redshift (z ∼ 7) galaxies (Schaerer et al. 2016). All these studies point towards LAEs being key to the reionization of the universe. Future observations of indirect tracers of LyC escape for this galaxy population at high redshift will allow us to understand how ξ_ion relates to LyC (and Ly α) escape. Regarding the latter, we find a tentative anticorrelation between f_esc(Ly α) and ξ_ion,0, which has two interesting implications. First, the escape of Ly α photons has been shown to have a similar behaviour with Ly α peak separations (e.g. Gazagnes et al. 2020), which are related to the porosity of the ISM. A more porous ISM favours low density channels through which Ly α (and LyC) photons can escape. Secondly, the fact that a higher ionizing photon production is not linked to higher Ly α escape fractions, points toward a time delay between the production and (potentially) escape of ionizing photons, in agreement with simulations (Barrow et al. 2020; Katz et al. 2020). This time delay is related to the creation of low-density channels in galaxies (for example, through SNe feedback).

Somewhat interestingly, we do not find any significant correlation between ξ_ion,0 and the galactic properties derived with Prospector (table provided in Appendix  B), however, the median values for our sample indicate it is dominated by metal-poor (log₁₀(Z/Z_⊙) |$= -2.2^{+0.4}_{-0.4}$|⁠), low-mass (log₁₀(M_*/M_⊙) |$= 8.3^{+0.7}_{-0.6}$|⁠) galaxies. We find that the lookback time at which half of their mass was assembled is given by |$t_{50} = 0.2^{+0.1}_{-0.1}$| Gyr, and that they have relatively high-ionization parameters (log〈U〉 |$= -2.4^{+0.4}_{-0.3}$|⁠). Interestingly, the LAE with the highest ξ_ion,0, JADES-GS 53.13859-27.79024, also has the highest stellar mass (log₁₀(M_*/M_⊙) ∼10), shown as the circle with a red edge in Fig. 5. This LAE was flagged as an AGN candidate due to its high stellar mass in Helton et al. (2023). Moreover, slitless spectroscopy reveals broad Hα emission, indicative of AGN activity. This object will be discussed in detail in a future work.

The redshift distribution of the sources in this work is shown on the top panel of Fig. 5. Given the small footprint of the MUSE coverage, groups of sources close in redshift space suggest clustering. In particular, we can distinguish two such groups, at z ∼ 5.4 and z ∼ 5.9, respectively. The former was identified and presented in Helton et al. (2023). It is not surprising that the LAEs we observe at the tail end of the EoR seem to be clustered, since there is growing evidence for these kind of objects to lie in overdense environments (e.g. Jung et al. 2020; Tang et al. 2023). Moreover, a recent study by Witten et al. (2023), which combines the power of simulations with observational data, finds that LAEs deep into the EoR necessarily reside in clustered environments. The clustering acts as a mechanism that aids in the creation of ionized bubbles, which favour Ly α escape as the IGM becomes more opaque. It is possible that the environment of these LAEs also has an effect on ξ_ion,0, but further studies are needed in order to draw solid conclusions.

Finally, we call attention to the biased nature of our sample. The IGM transmission makes Ly α detection at high redshifts a difficult task. Here we study the brightest LAEs at the end of the EoR, in other words, they are bright enough that we can reliably identify the Ly α emission line, which is later used to infer systemic redshifts. LAEs have been shown to be LyC leakers Gazagnes et al. (2020); Nakajima et al. (2020); Izotov et al. (2022), making them interesting objects to analyse. Our sample is not necessarily representative of the entire galaxy population at the EoR; this is a first exploratory work to showcase the science that can be performed using JEMS. However, as previously mentioned, there is increasing evidence that LAEs might be key drivers of the reionization of the universe. A follow-up study in the future will focus on a more complete sample using the entire JADES data base. Of particular interest, a recent study by Ning et al. (2023) focuses on the same class of objects at z ∼ 6, shown as black arrows in Fig. 5. Their sample was constructed by selecting spectroscopically confirmed LAEs observed with the fiber-fed, multi-object spectrograph Michigan/Magellan Fiber System (M2FS; Mateo et al. 2012), in a program that aims to build a sample of high-z galaxies (Jiang et al. 2017). Therefore, even though it was built in an analogous way to the sample presented here, it differs in the instrumentation used. Promisingly, our findings are broadly consistent with theirs.

7 SUMMARY AND CONCLUSIONS

We study a sample of 30 bright LAEs at the end of the EoR (z ∼ 5.4 − 6.6), with spectroscopic redshifts provided by the MUSE team (Bacon et al. 2023). The Hα recombination line is shifted into wavelengths observed by JEMS in this redshift range. We aim to retrieve the Hα intrinsic luminosity using a combination of the F430M, F460M, and F480M bands, while using Prospector to fit JEMS and selected HST photometry to constrain the UV continuum slope. We finally use our measurements to estimate the ionizing photon production efficiency ξ_ion,0. The blue UV continuum slopes of our sample indicate little-to-no dust attenuation. Additionally, due to the uncertainties in dust corrections of nebular emission lines and stellar continuum, we present uncorrected values in Table 2. In summary

• we use photoionization models to quantify the importance of [N ii] contamination in our photometric bands, and find that this effect would only be significant for cases where log〈U〉 ≤−3.0, which is not expected from galaxies at z ∼ 6 (Sugahara et al. 2022). Moreover, Prospector predicts the median log〈U〉 of the sample to be log〈U〉 |$= -2.3^{+0.2}_{-0.3}$|⁠.

• We divide our sample into four redshift bins and estimate the Hα flux using combinations of JEMS and find the median luminosity of the sample is log|$\frac{L(\rm{H\alpha })}{\text{erg s}^{-1}}$||$= 41.41^{+0.44}_{-0.15}$|

• We use Prospector to fit the SEDs of our sample and find blue UV slopes obtained by fitting a straight line in the rest-frame λ = 1250 − 2600 Å spectral region. The median of the sample is given by |$\beta = -2.09^{+0.23}_{-0.21}$|⁠, and is consistent with young stellar populations with little-to-no dust. We use the fitted lines to estimate the monochromatic luminosity density at rest-frame 1500 Å

• We use our measurements to estimate the ionizing photon production efficiency assuming no dust attenuation and find log|$\frac{\xi _{ion,0}}{\text{Hz erg}^{-1}} = 25.44^{+0.21}_{-0.15}$|⁠. If a Calzetti attenuation law is assumed instead, with local relations between nebular and continuum relative attenuation, this value is reduced by ∼0.1 dex for the few cases where dust attenuation is expected to be non-negligible. Our ξ_ion,0 estimations are in broad agreement with other LAEs at similar redshifts (Ning et al. 2023), and are slightly enhanced compared to their expected value based on all the measurements from literature (from z ∼ 0 − 9). The highest ξ_ion,0 is found in JADES-GS53.13859-27.79024, which has the highest Ly α emission and shows indications of possible AGN activity, this galaxy will be analysed in detail in a future work.

In this work we focus on LAEs, making us biased towards higher ξ_ion,0 values. This is an interesting galaxy population to study because LAEs are likely to be the main drivers of the reionization of the universe. In the future we will produce a similar study on a larger sample, through a stellar mass selection in the JADES data.

ACKNOWLEDGEMENTS

These observations are associated with JWST Cycle 1 GO program number 1963. Support for program JWST-GO-1963 was provided in part by NASA through a grant from the Space Telescope Science Institute, which is operated by the Associations of Universities for Research in Astronomy, Incorporated, under NASA contract number NAS 5–26555. The authors acknowledge the FRESCO team led by PI Pascal Oesch for developing their observing program with a zero-exclusive-access period. The authors acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315. The research of CCW is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. WB, JW, and LS acknowledge support from the ERC Advanced Grant 695671, ‘QUENCH’, and the Fondation MERAC. AJB and AS acknowledge funding from the ‘FirstGalaxies’ Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement Number 789056). BER acknowledges support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015. CW thanks the Science and Technology Facilities Council (STFC) for a PhD studentship, funded by UKRI grant number 2602262. ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1). EE, BDJ, and FS acknowledge support from the JWST/NIRCam contract to the University of Arizona NAS5-02015. DJE is supported as a Simons Investigator and by JWST/NIRCam contract to the University of Arizona, NAS5-02015. KB acknowledges support from the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. RM and CS acknowledge support by the Science and Technology Facilities Council (STFC) and by the ERC through Advanced Grant number 695671 ‘QUENCH’, and by the UKRI Frontier Research grant RISEandFALL. RM also acknowledges funding from a research professorship from the Royal Society. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton–Kavli Junior Fellowship. RS acknowledges support from a STFC Ernest Rutherford Fellowship (ST/S004831/1).

DATA AVAILABILITY STATEMENT

The data underlying this article will be shared on reasonable request to the corresponding author.

References

APPENDIX A: PROSPECTOR FITS

As in Fig. 4, we show the Prospector best fit SED and selected photometry for the remainder of the sample, including JEMS 30 × 30 pixel² (0.9 × 0.9 arcsec²) cutouts. For clarity, a Gaussian interpolation has been applied to the cutouts. The information relevant to the UV continuum slope fitting is also provided.

APPENDIX B: PROSPECTOR DERIVED PROPERTIES

Here we present the galactic properties derived with Prospector. All errors represent the 25 per cent and 75 per cent percentiles of each parameter.