Extreme emission line galaxies detected in JADES JWST/NIRSpec – I. Inferred galaxy properties

ABSTRACT

Extreme emission line galaxies (EELGs) exhibit large equivalent widths (EW) in their rest-optical emission lines ([O iii]|$\lambda 5007$| or H |$\alpha$| rest-frame EW|$\gt 750$| Å) which can be tied to a recent upturn in star formation rate (SFR), due to the sensitivity of the nebular line emission and the rest-optical continuum to young (⁠|$\lt 10$| Myr) and evolved stellar populations, respectively. By studying a sample of 85 star-forming galaxies (SFGs), spanning the redshift and magnitude interval |$3 \lt z\lt 9.5$| and |$-16\gt $|M|$\rm _{UV}\gt -21$|⁠, in the JWST Advanced Deep Extragalactic Survey (JADES) with NIRSpec/prism spectroscopy, we determine that SFGs initiate an EELG phase when entering a significant burst of star formation, with the highest EWs observed in EELGs with the youngest luminosity-weighted ages (⁠|$\lt 5$| Myr) and the highest burst intensity (those with the greatest excess between their current and long-term average SFR). We spectroscopically confirm that a greater proportion of SFGs are in an EELG phase at high redshift in our UV-selected sample (⁠|$61\pm 4~{{\ \rm per\ cent}}$| in our |$z\gt 5.7$| high-redshift bin, compared to |$23^{+4}_{-1}\%$| in our lowest redshift bin |$3\lt z\lt 4.1$|⁠) due to the combined evolution of metallicity, ionization parameter, and star formation histories with redshift. We report that the EELGs within our sample exhibit a higher average ionization efficiency (⁠|$\log _{10}(\xi _\mathrm{ion}^\mathrm{HII}/{\rm erg^{-1}Hz}) =25.5\pm 0.2$|⁠) than the non-EELGs. High-redshift EELGs therefore comprise a population of efficient ionizing photon producers. Additionally, we report that 53 per cent (9/17) of EELGs at |$z\gt 5.7$| have observed Ly |$\alpha$| emission, potentially lying within large ionized regions. The high detection rate of Ly |$\alpha$| emitters in our EELG selection suggests that the physical conditions associated with entering an EELG phase also promote the escape of Ly |$\alpha$| photons.

1 INTRODUCTION

Extreme emission line galaxies (EELGs) exhibit significant line emission relative to their stellar continuum which manifest as large emission line equivalent widths (EWs). These large emission line equivalent widths are most commonly seen in the H |$\alpha$| and [O iii]|$\lambda 5007$| lines in the rest-frame optical. The nebular emission lines are often driven by the ionizing photons produced in massive and short-lived O and B stars (or active galactic nuclei), whereas the surrounding rest-optical continuum includes the contribution from longer-lived and less massive stars (e.g. Eldridge & Stanway 2022). Hence, the EW tells us about the ratio of the very young to the older stellar population, reflecting any change in star formation over time. Selecting EELGs within a galaxy population may identify galaxies going through an upturn or burst in star formation (e.g. Endsley et al. 2024), exhibiting large specific star formation rates (sSFR).

Over the last decade EELGs have been observed locally (e.g. the ‘green pea’ population, Cardamone et al. 2009; Amorín, Pérez-Montero & Vílchez 2010; Izotov, Guseva & Thuan 2011; Brunker et al. 2020; Kumari et al. 2024), at moderate redshifts (⁠|$z\sim 2$|⁠, Atek et al. 2011; van der Wel et al. 2011; Amorín et al. 2015; Maseda et al. 2018; Tang et al. 2019; Du et al. 2020; Onodera et al. 2020; Boyett et al. 2022a) and during the epoch of reionization (EOR; Smit et al. 2015; De Barros et al. 2019; Endsley et al. 2021; Simmonds et al. 2023). These large EW systems are identified either directly from spectroscopy or from a photometric excess of flux between two adjacent imaging filters due to the contribution of line emission to one of them. We note that there is no set definition for what EW threshold identifies an EELG, with adopted thresholds varying between studies ([O iii]|$\lambda 5007$| EW thresholds ranging from |$\sim 100\, {\text {to}}\, 1000$| Å, e.g. der Wel et al. 2011; Amorín et al. 2015; Tang et al. 2019; Du et al. 2020; Tang et al. 2021a). Studies over a wide range of redshifts have suggested that EELGs are more abundant at higher redshift (Boyett et al. 2022a), perhaps due to changes in the characteristic star formation histories (SFHs) and other properties such as metallicity (e.g. Matthee et al. 2023). However, prior to the JWST, studies of EELGs at high redshifts (⁠|$z\gt 4$|⁠) had been limited to photometric excess measurements, largely from Spitzer/IRAC, as the rest-optical H |$\alpha$| and [O iii]|$\lambda 5007$| lines moved out of the wavelength coverage of the Hubble Space Telescope (HST) and ground-based spectroscopy (due to atmospheric opacity), and hence these measurements are highly dependent on assumptions of the underlying continuum shape and line ratios.

The recent launch of JWST enables sensitive near-infrared spectroscopy out to 5.2 |$\mu$|m with NIRSpec (Jakobsen et al. 2022) permitting the direct measurement of the EWs of rest-optical emission lines out to high redshift (H |$\alpha$||$z\lt 7$|⁠, [O iii]|$\lambda 5007$||$z\lt 9.5$|⁠). This allows the unambiguous identification of EELGs and characterization of their physical properties. As spectroscopic samples grow, we will be able to constrain the fraction of the galaxy population in an EELG phase, and any evolution with redshift.

In the first cycle of JWST observations, studies with a variety of instrument modes have been used to identify EELGs at high redshifts. Spectroscopic studies in slitless mode using NIRISS (Willott et al. 2022) and NIRCam (Rieke et al. 2023a) have identified EELGs (Boyett et al. 2022b; Kashino et al. 2023; Sun et al. 2023b). The use of NIRCam medium-band imaging to measure the flux excess between adjacent filters has also enabled EELG identification (Tacchella et al. 2023; Withers et al. 2023; Endsley et al. 2024). These programmes are beginning to confirm that the galaxy properties of individual high-redshift EELGs match what had been observed at cosmic noon (z|$\sim 2$|⁠) where studies have shown that galaxies with large EWs (rest-frame [O iii]|$\lambda 5007$| EW > 750 Å, for instance) are observed to exhibit lower masses, more compact morphologies, younger stellar populations, higher ionizing photon production efficiencies, higher ionization parameters, and lower gas-phase oxygen abundances (metallicity) than typical star-forming galaxies (SFGs) found at the same epoch (e.g. Tang et al. 2019).

These properties make SFGs in an EELG phase effective ionizing photon producers. However, at moderate redshifts, their abundance has been measured to be too low to dominate the ionizing output of the star-forming galaxy population (Boyett et al. 2022a), with the majority of the ionizing output coming from SFGs in more typical modes of star formation. The bursty SFHs expected at higher redshifts and lower masses (Ceverino, Klessen & Glover 2018; Faucher-Giguère 2018; Ma et al. 2018; Tacchella, Forbes & Caplar 2020), which may make EELGs more common than at lower redshifts, suggest the potential for high escape fractions of ionizing photons¹ (⁠|$f_\mathrm{esc}$|⁠; Katz et al. 2023). To understand the nature and significance of EELGs at high redshift, it is critical to study a large spectroscopic sample.

Our sample originates from the JWST Advanced Deep Extragalactic Survey (JADES; Bunker, NIRSPEC Instrument Science Team & JAESs Collaboration 2020; Rieke 2020; Eisenstein et al. 2023), which is obtaining spectra of thousands of galaxies between cosmic noon to within the EoR. In this paper, we use the JADES first data release (Bunker et al. 2023a; Rieke et al. 2023b), which targets |$\sim 250$| galaxies in the Hubble Ultra Deep Field (HUDF) region (Beckwith et al. 2006) and the surrounding GOODS-South field (Giavalisco et al. 2004). We use the JADES NIRSpec and NIRCam data to look at SFGs over a redshift range |$3\lt z\lt 9.5$| to investigate the evolution of the abundance of these EELGs with redshift and study their physical properties in detail to determine what physical conditions initiate an EELG phase. We focus on the [O iii]|$\lambda 5007$| line as this typically has the highest EW in our spectroscopy, but we also consider cases where H |$\alpha$| (accessible at |$z\lesssim 7$|⁠) may also have extreme EWs.

This paper is laid out as follows. We discuss the JADES observations and the selection of the parent spectroscopic sample in Section 2. We measure the EWs of emissions lines and identify a subsample of galaxies in an EELG phase in Section 3. In Sections 4, we examine any trends between the EWs of [OIII]|$\lambda 5007$| and galaxy properties. Finally, in Section 5, we discuss what initiates an EELG phase and how the fraction of SFGs in an EELG phase evolves as a function of redshift. Where applicable, we use a standard |$\Lambda$|CDM cosmology with parameters H|$_{0}=70$| km s⁻¹ Mpc⁻¹, |$\Omega _{m}=$|0.3, and |$\Omega _{\wedge }=$|0.7. All magnitudes are in the AB system (Oke & Gunn 1983). When quoted, all EWs are in the rest-frame.

2 OBSERVATIONS

2.1 Spectroscopic and imaging data

The spectroscopic and photometric imaging data used in this study come from the DEEP tier of the JWST NIRSpec and NIRCam observations obtained as part of the JADES survey (programme ID 1210, 1180 P.I.: D. Eisenstein). The spectroscopic data used here have already been presented in Bunker et al. (2023a) and we refer the reader to that paper for a detailed description of the target selection, observation strategy and data reduction (further reduction details will be provided in Carniani in preparation). The NIRCam imaging used in this paper is described in Rieke et al. (2023b). We also refer the reader to Eisenstein et al. (2023) for an overview of the JADES programme.

Briefly, the spectroscopic observations were obtained using NIRSpec in the multiobject spectroscopy (MOS) mode (Ferruit et al. 2022), with a three shutters nodding sequence. Three microshutter assembly (MSA) configurations were constructed, each with three nod positions. Targets were assigned to the configurations following the prioritization discussed in Bunker et al. (2023a). The highest priority targets were included in all three MSA configurations and received the maximum exposure time, while the remaining targets assigned to only one or two configurations received only 1/3 or 2/3 of the total observing time. This approach maximizes the total number of galaxies observed while retaining the maximum exposure time for the highest priority sources. The prioritization of targets that were selected and allocated may introduce a selection bias to our sample, which we discuss further in Section 2.2.

The configurations were observed with both the low-dispersion prism² (total integrated exposure time of 100 ks) and the |$R\sim 1000$| medium resolution³ G140M/F070LP, G235M/F170LP, and G395M/F290LP gratings (each with a total integrated exposure time of 25 ks). The MSA was configured to avoid overlap in the prism spectra. To avoid spectra overlapping in the medium-resolution grating (where the dispersed light covers a greater pixel range on the detector) for the highest priority targets, a number of lower priority targets were removed from the configuration. The spectroscopic data were reduced using the pipeline developed by the ESA NIRSpec Science Operations Team and the NIRSpec GTO Team [Carniani in preparation, see Bunker et al. (2023a) for specific details].

For the imaging in this study, we utilise the JADES (Rieke et al. 2023b) F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W NIRCam observations. Each of these have integrated exposure times between 24.7 and 60.5 ks (see Eisenstein et al. 2023; Rieke et al. 2023b). The JADES F444W NIRCam imaging is supplemented with extended F444W imaging from the FRESCO survey (Oesch et al. 2023) and covers a larger footprint than the other filters. The JADES imaging over the HUDF is complimented by the JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023) F182M, F210M, F430M, F460M, and F480M observations. These have integrated exposure times between 13.9 and 27.8 ks. We note that the footprint in the different filters varies, and a subset of the galaxies targeted with NIRSpec have coverage only in F444W, while the majority are covered by all imaging filters.

2.2 Sample selection and selection bias

Across the three MSA configurations, a total of 252 unique galaxies were observed in the prism, with a subset of 197 also observed in medium-resolution gratings [we refer the reader to Bunker et al. (2023a) for a detailed break-down of the spectroscopic survey].

Of these, a total of 152 galaxies have identifiable redshifts based on at least one emission line detected at |$\gt 5\sigma$| in the prism spectra or have an identifiable Lyman spectral break. For the grating, 103 out of 197 galaxies observed had an identifiable redshift.

In this study, we will focus on the emission line sample at |$z\gt 3$|⁠, where the H|$\beta$| and [O iii]|$\lambda 4959,5007$| emission lines are deblended at the prism resolution. We define a sample of 86 galaxies which are determined to have multiple detected lines at |$\gt 5\sigma$|⁠, including [O iii], in the prism spectra. These galaxies will become our sample for this study (hereafter ‘full sample’).

In Fig. 1, we present the redshift distribution for the JADES targeted galaxies (presenting NIRSpec spectroscopic redshifts when determined otherwise photometric redshifts, see Eisenstein et al. 2023; Bunker et al. 2023a) and overlay the subset of galaxies which have multiple |$\gt 5\sigma$| emission line detections at |$z\gt 3$|⁠, which form the full sample in this study. We note that [O iii]|$\lambda 5007$| passes out of the NIRSpec wavelength coverage at |$z\gt 9.5$|⁠, providing an upper redshift limit on our study.

This sample was preselected for observation through a target prioritization strategy (as detailed in Bunker et al. 2023a). The selection of the highest redshift candidates, which were initially identified as dropout candidates from HST imaging at |$z\gt 5.7$| (these include I|$_{775}$|⁠, Y|$_{105}$|⁠, and J|$_{125}$|-band dropouts), used a rest-frame UV magnitude threshold longwards of the Ly |$\alpha$| break (rest-frame UV detection band AB magnitude |$\lesssim 29$|⁠) with prioritization based on their photometric redshift. Over a subset of the area covered by the MSA, JWST NIRCam imaging was also available prior to the MSA design and we used this to refine the selection with better S/N at longer wavelengths.⁴ We note that, as expected, all the high-redshift (⁠|$z\gt 5.7$|⁠) EELGs we identify come from targets selected to be in the photometric redshift range |$5.7\lt z\lt 8.5$| [priority classes 4, 6.1, and 6.2 in Bunker et al. (2023a) depending on their HST or JWST/NIRCam magnitude] or |$z\gt 8.5$| (priority class 1 in Bunker et al. 2023a), and of the galaxies targeted in each of these classes roughly half proved to be an EELG (see Section 3.1). As we are essentially selecting objects on the rest-UV continuum, in a spectral region without strong line emission (filters used for the rest-frame UV typically do not encompass C iii]|$\lambda 1909$|⁠), we are not biasing our magnitude limited sample towards strong line emitters, although we consider the selection effects due to SFH in Section 4.3.1.

At lower redshifts (⁠|$z\lt 5.7$|⁠), priority was given to obtaining a stellar mass-limited sample, by selecting on the longest wavelength broad-band filter with good imaging data. In the case of HST selection this was the F160W HST/WFC3 (AB |$\lt 29$| mag), but for a subset of the NIRSpec footprint NIRCam imaging was available prior to designing the MSA configuration, enabling selection on the rest-frame optical F444W (AB |$\lt 27.5$| mag). The majority of our identified low-redshift EELG sample comes from targets selected to be in classes 7.5, 7.6, and 7.7 in Bunker et al. (2023a) (corresponding to redshift slices |$4.5\lt z\lt 5.7$|⁠, |$3.5\lt z\lt 4.5$|⁠, and |$2.5\lt z\lt 3.5$|⁠, respectively), except for one object from class 6.2 (ID: 8113 at |$z=4.90$|⁠) where the spectroscopic redshift was slightly below the priority class expectation, one object in class 7.2 (ID: 4270, |$z=4.02$|⁠, upweighted for scientific interest), which was selected as a potential quiescent galaxy from the photometry SED but we now can see was driven by line contamination rather than a strong Balmer/4000 Å break, and one object below the stated magnitude selection which was included as a filler (class 9, ID: 10001916 at |$z=4.28$|⁠). From the class 7.5, 7.6, and 7.7 targeted galaxies, the fraction of EELGs we detect are 8/23, 6/31, and 3/45 (see Section 3.1).

We note that the selection on the rest-frame optical means that over specific redshift ranges, the broad-band flux may be boosted by the contribution from strong line emission (H |$\alpha$|⁠, [O iii] + H|$\beta$|⁠) which could bias the lower redshift sample towards EELGs at faint magnitudes. We will discuss in Section 4.3.1 how we can apply an apparent rest-UV magnitude cut to this subsample to address this selection effect.

We also note that the rest-optical selection (at |$z\lt 5.7$|⁠) will include galaxies with less active star formation, which may account for the lower fraction of galaxies with detected emission lines at |$z\lt 6$|⁠, seen in Fig. 1.

3 MEASUREMENT OF STRONG LINE EQUIVALENT WIDTHS

To determine the EW for emission lines detected in the prism and R1000 grating spectroscopy, we require measurements of the line flux and the continuum flux density at the wavelength of the line.

We use the line fluxes reported in Bunker et al. (2023a) for the JADES DEEP spectroscopy, which relied on Gaussian modelling of the emission lines (with the width tied for blended lines) after continuum subtraction. All lines detected at a |$\gt 5\sigma$| significance are reported in the publicly available catalogue,⁵ and where necessary (e.g. the line diagnostic diagrams in Section 4), we additionally make use of line fluxes below this significance threshold, measured the same way.

The continuum flux density at the wavelength of the lines can be determined either directly from the spectroscopy or from aperture photometry of the NIRCam imaging. We choose to measure the continuum flux density from NIRCam photometry as this provides a better constraint on the continuum for faint sources, where the continuum is often below the sensitivity of the spectroscopy. To determine the continuum flux density from the NIRCam imaging, we first wish to avoid requiring a path loss correction (to account for slit losses associated with the NIRSpec MSA) and hence the continuum flux densities are derived from apodized photometry, with the aperture set to match the placement of the MSA shutter over the target. A point spread function (PSF) correction is then applied to account for the difference in PSF between the NIRCam and NIRSpec observations. Bunker et al. (2023b) show the consistency between the spectral flux through the microshutter and the apodized multiple filter NIRCam photometry over the same wavelength ranges. The continuum flux density around the wavelength of the emission line is derived assuming a flat in |$f_\nu$| continuum (matching the typical observed |$\beta$| slopes in our full sample, see Section 4.3.1) after removing any contribution from spectroscopically detected line emission to the broad-band filter, accounting for the wavelength-dependent filter transmission. The observed-frame EW for each line is then taken as the ratio of the individual line flux to the continuum flux density (⁠|$f_\lambda$|⁠) at the location of the line. This is then corrected to the rest-frame EW value by dividing by|$(1+z)$|⁠.

The use of broad-band photometry to determine the continuum flux density allows EWs to be determined for faint galaxies where constraints cannot be placed directly from spectroscopy. For a subset of continuum bright objects, where the continuum is well detected in the spectra, we check whether both methods recover consistent EW values, and present this work in Appendix  A. From this, we determine consistent EWs measurements for galaxies with a signal to noise per pixel in continuum greater than 3. It is clear that the uncertainty in the spectrum-derived EW increases as the continuum becomes fainter, and this supports our choice to use the broad-band photometry to provide tighter constraints on the EWs.

When computing the continuum flux density from the broad-band photometry, we use the filter that contains the selected emission line. However, in a small number of cases, we do not have filter coverage at the wavelength of investigation. This accounts for 10/86 of our [O iii]|$\lambda 5007$| EW measurements. In these cases, we determine an estimate of the continuum flux density at the wavelength of the line using the nearest available filter. We find that using an adjacent filter provides a consistent EW estimate, albeit with a larger scatter driven by the variation of UV spectral slopes (⁠|$\beta$|⁠) in our sample. We discuss this further in Appendix  A.

Finally, we compare the EWs derived using line fluxes determined from the two spectral resolution modes (prism and R1000 grating). We determine good agreement for our EWs for lines detected in both the prism and R1000 spectra, with the R1000 providing higher resolving power for blended emission lines (e.g. the [O iii]|$\lambda 4959, 5007$| doublet becomes blended below |$z\lesssim 5.2$| in the prism). However, we note a small systematic offset between the two modes with the R1000 fluxes roughly 8 per cent larger than the prism fluxes. This offset, due to the flux calibration of the spectroscopy, has been observed previously by Bunker et al. (2023b) who find that the prism flux calibration better replicates the broad-band flux densities. Therefore in this study, we will adopt the prism EWs as our primary method. We note that because of this choice, we limit our sample to |$z\gt 3$| where the H |$\beta$| and [O iii] lines are unblended in the prism.

For one target (ID: 7624), two galaxies fell within the NIRSpec MSA shutter (at |$z=2.7\ \&\ z=4.8$|⁠). Despite observing strong emission lines at the expected wavelengths, the significant contribution of continuum flux from both systems means the emission line EW cannot be determined from either the spectrum or the photometry. We therefore remove this source from our full sample (leaving 85 galaxies) and our stacking analysis.

3.1 Identification of extreme emission line galaxies

There is no set EW threshold (or even specific emission line) for defining a star-forming galaxy as being in an EELG phase, with adopted thresholds varying from study to study. Commonly taken thresholds can be based on the EW of individual lines (e.g. [O iii]|$\lambda 5007$| or H |$\alpha$|⁠) or complexes (e.g. [O iii] + H |$\beta$|⁠) depending on the spectroscopic resolution of observations or if relying on broad-band flux-excess. Threshold values for [O iii]|$\lambda 5007$| EW in literature have ranged from |$\sim 100\,{\!-\!}\,1000$| Å  (e.g. der Wel et al. 2011; Amorín et al. 2015; Tang et al. 2019, 2021a; Du et al. 2020 ).

In this paper, we will adopt [O iii]|$\lambda 5007$| EW |$\gt 750$| Å as our threshold, based on work by Boyett et al. (2022a) and Tang et al. (2019). This choice is motivated by |$z\sim 2$| studies which show significantly different interstellar medium (ISM) characteristics (high ionization parameter and ionization efficiency as well as potentially high escape fractions of ionizing radiation) above this threshold compared to the general SFG population at this epoch (Tang et al. 2019).

We identify 35 galaxies that meet this threshold for [O iii]|$\lambda 5007$| in our sample, with one additional galaxy consistent within |$1\sigma$| of the [O iii]|$\lambda 5007$| threshold (ID: 10001892) which also has a H |$\alpha$| EW above 750 Å  and hence is an EELG. This makes up 42 per cent (36/85) of our sample of emission line galaxies at |$z\gt 3$|⁠. This subsample of 36 galaxies will hereafter be referred to as the EELG sample. All other galaxies in our full sample below this threshold will therefore be referred to as the non-EELG sample (49 galaxies). The emission line EWs for these EELGs are presented in Table 1. We note that one EELG, ID: 8083 has an underlying broad component of H |$\alpha$| (Maiolino et al. 2023) and likely hosts AGNs. One non-EELG in our full sample, ID: 10013704 similarly shows broad H |$\alpha$| (Maiolino et al. 2023). In Table 1, we flag these, and in Figs 12 and 14 we plot these two with different symbols. However, we retain these in our sample, since the broad-line contribution to the line fluxes is small, and the narrow lines and galaxy continuum may well be dominated by star formation rather than an AGN (Maiolino et al. 2023).

3.2 Distribution of equivalent widths

In Fig. 2, we present the EW distribution of the [O iii]|$\lambda 5007$| and H |$\alpha$| EWs. We split our full sample into two redshift bins, those above and below |$z=5.7$|⁠, dividing the EELG sample into two at the end of the EoR. We expect that the galaxies within our sample at high redshift (⁠|$z\gt 5.7$|⁠), which were selected based on the rest-UV, should avoid any selection bias towards high EW galaxies, whereas lower redshift galaxies within our sample are selected on the rest-optical and hence may preferentially include galaxies with large EWs in their rest-optical emission lines (see Section 2.2). However, although there are even numbers of EELGs from the high and lower redshift samples, the lower redshift sample contains a larger total number of galaxies (i.e. the EELG fraction is lower than in the higher redshift sample), suggesting significant evolution in the EELG fraction with redshift (discussed further in Section 5.3).

The EW distribution of H |$\alpha$| and [O iii]|$\lambda 5007$| show similar profiles (see Fig. 2). The high-redshift sample exhibits a flatter EW distribution which extends to high EW (EW|$\gt 2000$| Å), whereas at low-redshift the distribution peaks around low EW and while there is a high-EW tail this represents only a small fraction of the total population.

4 GALAXY PROPERTIES OF EELGS

In this section, we are now interested in characterizing the physical properties of our EELG sample, both individually and averaged over two redshift bins split at |$z=5.7$|⁠. We will compare these characteristics against our non-EELG sample to examine what attributes may be responsible for a galaxy entering an EELG phase. We will first discuss the combination of several spectra (stacking) in Section 4.1, then examine the galaxy properties in Sections 4.2–4.4.

4.1 Stacked spectra of EELGs

In this section, we stack the prism spectra for four different subsamples of our data set, split by redshift (⁠|$z=5.7$|⁠) and whether a galaxy meets the threshold to be considered an EELG ([O iii]|$\lambda 5007$| or H |$\alpha$| EW |$\gt 750$| Å). There are 17 EELGs in our high-redshift sample and 19 in our low-redshift sample. We additionally have 11 high-redshift non-EELGs and 38 low-redshift non-EELGs (excluding ID:7624 as described in Section 3) that we can use as a control sample, although we note that these still require at least two emission lines detected at |$\gt 5\, \sigma$|⁠.

We note that in our low-redshift non-EELG sample there is one galaxy (ID 4493, |$z=3.59$|⁠) with an extremely red UV spectral slope of |$\beta =+0.65\pm 0.03$|⁠, the only galaxy with a positive |$\beta$|⁠, which is an outlier at |$\gg 3\sigma$| from the full sample. We exclude this highly unusual galaxy from our stack, leaving 37 galaxies in our low-redshift non-EELG stack.

We de-redshift each individual prism spectrum to the rest-frame (with the redshift determined from the detected [OIII]|$\lambda 5007$| emission line), interpolating onto a common wavelength grid of 10  Å pixels, and normalize to the integrated 1400–1600 Å  flux (⁠|$f_\lambda$|⁠), where possible determined from the NIRSpec prism spectrum (if the S/N is sufficient). For redshifts below |$z=4$| the NIRSpec spectrum was not sensitive to this portion of the rest-frame UV so we instead used the rest-UV flux inferred from the broad-band photometry. At higher redshifts (⁠|$z\gt 4$|⁠), if the spectroscopic continuum is only weakly detected in the rest-UV (integrated S/N|$\lt 3$| in the range 1400–1600 Å) we again adopt the M|$\rm _{UV}$| inferred from the broad-band imaging, as described in the Appendix  B, this only applied to two of our galaxies (ID 10009693 and 10009320). We then average these UV normalized spectra; we note that because of the different redshifts of each galaxy, we do not have uniform coverage of all wavelengths (as shown in Fig. 4) and this is taken into account in the averaging. The top panel of Fig. 3 presents the four subsets over the full wavelength range, and we show individual zooms around the Ly |$\alpha$| and [O iii] + H |$\beta$| wavelengths in the lower panels. In Appendix  C, we consider a different stacking where we weight by the rest-frame UV luminosity and find consistent results.

To measure the emission line fluxes from the stacks, we note that the combination of galaxies at different redshifts, and hence at different resolving power, means a Gaussian profile no longer provides an ideal fit. We instead measure the integrated line flux rather than attempting to fit a Gaussian profile. We report the galaxy properties of the stacks in Table 2 and we will discuss the properties in the following subsections, along with a discussion of the distribution of individual galaxy measurements.

4.2 ISM properties from emission line diagnostics

In this section, we look at the physical properties of the individual galaxies (both for EELGs and non-EELGs) in our full sample and our stacked subsamples. In Section 4.2.1 we examine the detected hydrogen emission lines and we consider how the dust properties, [O iii]|$\lambda 5007$|/H |$\alpha$| flux ratio, H |$\alpha$| and H|$\beta$| EWs, and the detection of Ly |$\alpha$| emission vary across our sample. In Section 4.2.2, we examine how the ionization parameter, metallicity and ionizing photon production efficiency change with [O iii]|$\lambda 5007$| EW.

4.2.1 Hydrogen lines and the Balmer decrement

In many cases in our full sample (61/85), we detect both H |$\alpha$| and H |$\beta$| at |$\rm S/N\,\gt\, 5$|⁠, and we will first use these lines to estimate the reddening of the nebular emission using the Balmer decrement, assuming an intrinsic H |$\alpha$|/H |$\beta =2.86$| (Osterbrock & Ferland 2006 with electron temperature and density, |$T_e=10^4\,{\rm K}$| and |$n_e=10^2\,{\rm cm}^{-3}$|⁠, respectively). At the highest redshifts (⁠|$z\gt 7.1$|⁠), we lose H |$\alpha$| from our wavelength coverage and instead we consider the H|$\beta /$|H |$\delta$| decrement⁶ (with an intrinsic value of 3.88, Osterbrock & Ferland 2006). This applies to six galaxies in our full sample of which four are EELGs. For two of these galaxies (both EELG IDs: 1003905 and 10058975), we determine a H|$\beta /$|H|$\delta$| decrement consistent with no dust, but for the other four no other Balmer lines are detected leaving the dust attenuation in these galaxies unconstrained (hence we will exclude these four targets from any dust corrected analysis).

We use the Balmer decrement to infer the attenuation at other wavelengths, and we adopt the Calzetti et al. (2000) attenuation law (with |$R_v=4.05$|⁠) assuming that the reddening of the stellar continuum is related to that of the nebular emission lines by |$E(B-V)_{\mathrm{stellar}}=0.44\, E(B-V)_{\mathrm{nebular}}$|⁠. For the galaxies within our full sample where we can constrain the dust reddening, we measure the mean dust attenuation at 1600 Å  to be A|$_{1600}= 0.48\pm 0.07$|⁠. The EELG subsample exhibits a marginally lower level of attenuation of A|$_{1600}=0.38\pm 0.08$|⁠. In Section 4.3.1, we compare our measured Balmer decrements with the rest-frame UV spectral slopes, another potential indicator of dust attenuation, and we discuss the dust properties of our sample in that section.

We note that if the reddening of the nebular lines and stellar continuum were identical at the same wavelength, the EW measurements would not be affected by the dust reddening, unlike the observed line ratio diagnostics. However, in a scenario in which very young stars (preferentially powering nebular line emission) are still enshrouded in their birth clouds (e.g. Charlot & Fall 2000), the EW would be dependent on the differential reddening where nebular lines are attenuated to a greater extent than the surrounding stellar continuum (e.g. in the Calzetti et al. 2000 law). In our sample we do not correct the EW for this differential extinction, but we consider the potential effect of this in Section 4.3.1.

We now consider the relative strength of rest-optical emission line fluxes within our sample. In our sample of EELGs, it is predominantly the [O iii]|$\lambda 5007$| line which has the highest equivalent width in these extreme systems. H |$\alpha$| is usually the second strongest line in our NIRSpec spectroscopy. For the subset of 61 galaxies where both H |$\alpha$| and H |$\beta$| are detected above |$5\, \sigma$| and a de-reddening correction can be applied, we plot the ratio of [O iii]|$\lambda 5007$|/H |$\alpha$| line fluxes in Fig. 5. As can be seen, in our sample [O iii]|$\lambda 5007$| is nearly always stronger than H |$\alpha$|⁠, wi th an average flux ratio for the full sample of |$f_{\rm {[O\,{\small III}]}\lambda 5007}/f_{\rm H\,\alpha }=2.01\pm 0.07$| after correction for dust reddening. Splitting our full sample into EELGs (coloured blue in Fig. 5) and non-EELGs (coloured orange), we see only a marginally higher ratio in the EELG subsample with a |$f_{\rm {[O\,{\small III}]}\lambda 5007}/f_{\rm H\,\alpha }$| flux ratio of |$2.18\pm 0.15$|⁠, compared to |$1.93\pm 0.08$| for our non-EELG sample.

The high [O iii]/H |$\alpha$| flux ratio we observe at |$z\gt 3$| is in contrast to work at lower redshifts, where H|$\alpha$| is often more luminous than [O iii]|$\lambda 5007$|⁠. For example, Domínguez et al. (2013) at |$z=0.75{\!-\!}1.5$| find a [O iii]/H |$\alpha$| flux ratios spanning |$1.3-0.5$|⁠, with the ratio dropping at higher H |$\alpha$| luminosities (grey data points in Fig. 5). The significant increase we observe in the [O iii]|$\lambda 5007$|/H|$\alpha$| flux ratio at higher redshifts is likely attributable to decreasing metallicity and a higher ionization parameter (e.g. Kewley & Dopita 2002). Indeed, our observed flux ratio is consistent with the rising trend with redshift of [O iii]/H|$\beta$| (which we correct to H|$\alpha$| using the Case-B H |$\alpha$|/H |$\beta$| = 2.86) at |$z\gt 3$| predicted in Cullen et al. (2016), who suggest that increasing ISM pressure and ionization parameter is driving this.

We can additionally consider how the EW of these strong Balmer rest-optical emission lines varies with [O iii]|$\lambda 5007$| EW. The hydrogen and [O iii] nebular emission lines are powered by the same sources, typically the hot ionizing UV emission from hot massive stars. However, the EW of these lines will also be moderated by the stellar continuum at their respective rest-frame wavelengths. In the top two panels of Fig. 6 we present the H |$\alpha$| and H |$\beta$| EWs against the [O iii]|$\lambda 5007$| EWs. The EW of both Balmer lines show a clear positive trend with the [O iii]|$\lambda 5007$| EW, with best-fitting relations:

and

We additionally find that the EELG stacks (defined as having [O iii]|$\lambda 5007$|EW|$\gt 750$|Å) at high and low redshift exhibit consistent H |$\alpha$| EWs, which are both considerably larger than those exhibited by the non-EELG stacks (see Fig. 6). Here, we see that the high sSFRs that drive extreme EWs in the hydrogen nebular emission lines (e.g. Mármol-Queraltó et al. 2016) also contribute to the high [O iii]|$\lambda 5007$| EWs in our sample. The high H |$\alpha$| and H |$\beta$| EWs exhibited in our EELG sample likewise reflect that the stellar populations in these galaxies will likely be young (the EW of these hydrogen lines starts decreasing within 10 Myr for an instantaneous burst, e.g. Leitherer et al. 1999).

Finally, we turn to the Ly |$\alpha$| rest-UV emission line, which our NIRSpec prism spectroscopy has sensitive coverage of down to |$z\gt 4$|⁠. Although Ly |$\alpha$| is known to be resonantly scattered and highly suppressed in most galaxies, it has been suggested that extreme galaxies such as EELGs (with potentially high ionization of the ISM, and starburst-driven outflows) might be able to clear escape channels for Ly |$\alpha$| photons (Östlin et al. 2009; Herenz et al. 2017). We now look for any evidence of Ly |$\alpha$| emitters (LAEs) in our stacks and in individual EELGs, and we report strong evidence that Ly |$\alpha$| does indeed escape from EELGs at high redshift.

In the middle panel of Fig. 3, we zoom in on the rest-frame UV in the stacked data to study the emergence or otherwise of Ly |$\alpha$|⁠. We first consider the two EELG stacks. The high-redshift EELG stack has the strongest Ly |$\alpha$| emission (Ly|$\alpha$| EW = |$39\pm 3$| Å), with the low-redshift EELG stack also showing Ly |$\alpha$| in emission, but at a slightly weaker level (Ly |$\alpha$| EW = |$21\pm 5$| Å). However, we note that the lower sensitivity of NIRSpec at shorter wavelengths (⁠|$\lt 1\,\mu \rm m$|⁠) means that in our lower redshift stack we are less sensitive to weak Ly |$\alpha$| emission close to the spectral break (see Jones et al. 2024).

Turning to the non-EELG stacks, there is no Ly |$\alpha$| emission and both redshift bins have the appearance of a damping wing (Miralda-Escudé 1998). It is interesting that we observe a similar softening of the Lyman-break in both non-EELG redshift stacks. In our high-redshift stack, where we sample galaxies during the EoR, we may anticipate Lyman-damping due to a high hydrogen neutral fraction in the IGM, however, in the low-redshift stack we would not expect damping from the IGM because of the near-zero neutral fraction during this epoch (⁠|$z\lt 5.7$|⁠). Softening of the spectral break may instead be due to either damped Ly |$\alpha$| from H i in the surrounding environment [ISM and circumgalactic medium (CGM)] of the galaxies (e.g. Heintz et al. 2024) or be an artefact, resulting from the softening of the Lyman-break due to convolution with the line spread function of NIRSpec (see Jones et al. 2024).

Individually, we find that out of 17 EELGs at |$z\gt 5.7$|⁠, 53 per cent (9/17) have detected Ly |$\alpha$| emission, as determined by Jones et al. (2024), Witstok et al. (2024), and Saxena et al. (2024). This would suggest that selections based on high [O iii]|$\lambda 5007$| EW galaxies are also preferentially selecting LAEs. This LAE fraction is higher than among the high-redshift non-EELGs, where only 18 per cent (2/11) have detected Ly  |$\alpha$| emission.⁷ These two non-EELG Ly |$\alpha$| emitters show moderate [O iii]|$\lambda 5007$| EWs (⁠|$\sim 300$| Å), in both cases the relative strength of the nebular emission lines (being driven by recent |$\lt 10$| Myr star formation) has potentially decreased due to strong rest-optical continuum emission from an older stellar population.

Below |$z\sim 4$|⁠, where we have four EELGs in our sample, Ly |$\alpha$| does not lie within the wavelength window where we have sensitive spectral coverage with NIRSpec. Between |$4.0\lt z\lt 5.7$| we have a sample of 15 EELGS with 4 showing Ly |$\alpha$| in the NIRSpec prism spectroscopy. One of these is the potential AGN identified by Maiolino et al. (2023) to have a broad component to the H |$\alpha$| emission line (Galaxy ID: 8083). We additionally look at the MUSE Data release II IFU spectroscopy of the HUDF region (Bacon et al. 2023) and find one additional EELG galaxy (Galaxy ID 7938, |$z=4.815$|⁠) where MUSE detects Ly |$\alpha$| emission which is not visible in our low dispersion prism spectrum. This leaves us with an LAE fraction of 5/15 (33 per cent) of the EELG galaxies in the redshift range |$4.0\lt z\lt 5.7$|⁠, nominally lower than the fraction at higher redshifts (⁠|$z\gt 5.7$|⁠) of 53 per cent (9/17), although we note that the combination of sensitivity and spectral resolution of the prism means that NIRSpec’s ability to detect Ly |$\alpha$| decreases towards the short wavelength end (affecting the lower redshifts). Hence, the LAE fraction among our low-redshift EELG sample may be a lower limit and the difference in the fractions of LAEs in the low- and high-redshift EELG samples may not reflect genuine evolution. However, what we can say is that at high redshift, roughly half of EELGs are Ly |$\alpha$| emitters, whereas for the non-EELGs at these high redshifts only 18 per cent of our sample have Ly |$\alpha$| in the NIRSpec spectroscopy. This is a remarkable result which we will discuss further in Section 5.2.

Finally, we note 2/19 (11 per cent) non-EELGs at |$4\lt z\lt 5.7$| show Ly |$\alpha$| emission, the lowest LAE fraction in any of our subsamples. Like the high-z non-EELG LAEs these two.⁸ also show rest-optical continuum emission from an evolved stellar population.

4.2.2 Ionization parameter, excitation state, and metallicity

To investigate how the ionization state of the ISM changes with [O iii]|$\lambda 5007$| EW within our full sample we employ the O32⁹ and Ne3O2¹⁰ diagnostics (e.g. see Tang et al. 2019), which are sensitive to the ionization state of the ISM which is shaped by the ionizing flux and hardness of the ionizing spectrum. O32 provides a flux ratio between two different ionization states of the same elements, while Ne3O2 provides a flux ratio between two α-elements in different ionization states. The Ne3O2 diagnostic utilizes two lines close in wavelength and is therefore insensitive to the dust attenuation and does not require any correction, whereas the large baseline of O32 (3727–5007  Å) means we must take into account reddening in our measurements (as described in Section 4.2.1). Although the Ne3O2 ratio is sensitive to the Ne/O abundance ratio, these are both α-elements and this abundance ratio is expected to vary only minimally (e.g. Kobayashi, Karakas & Lugaro 2020). Together they provide a complimentary estimate of the ionization parameter.

In the bottom-left panel of Fig. 6 we plot the O32 diagnostic against the [O iii]|$\lambda 5007$| EW, with galaxies corrected for dust reddening. We note that when H |$\alpha$| and H |$\beta$| are detected and we can measure the Balmer decrement, the average reddening suggests minimal dust attenuation between 3727 and 5007Å. We additionally plot the Ne3O2 ratio against the measured [O iii]|$\lambda 5007$| EW of the full sample in the lower right panel. Both O32 and Ne3O2 are related to the ionization parameter, and in both these Figure panels we additionally provide the ionization parameter derived following the Witstok et al. (2021) relations.

We find a clear positive correlation in both O32 and Ne3O2 with [O iii]|$\lambda 5007$| EW, obtaining best-fitting relations:

and

This correlation suggests that galaxies with greater [O iii]|$\lambda 5007$| EW have more highly ionized ISMs. This result is not surprising considering young stellar populations are responsible for the nebular emission. The youngest stellar populations retain the most massive stars with the hardest ionizing spectra, and have higher ionizing photon fluxes, making the ISM highly ionized. EELGs show the greatest nebular emission relative to the continuum and therefore are likely to be the youngest systems (see Section 4.4)

When we turn to our stacked subsamples in Fig. 6, the high-redshift EELG stack has the highest O32 and Ne3O2 ratios, with the low-redshift EELG stack showing marginally lower values. Both EELG stacks show significantly higher values for all quantities in Fig. 6 compared to the non-EELG stacks.

To further examine the state of the ISM, we employ the dust-corrected O32 and R23¹¹ diagnostic diagram (Fig. 7), where galaxies of different ionization and excitation states occupy distinct regions. For our full sample, we determine the O32 and R23 emission line diagnostics for galaxies, which we correct for reddening (see Section 3). Where individual lines are undetected at |$3\sigma$|⁠, they are represented as |$3\sigma$| limits. These are plotted in Fig. 7. We find that our full sample (⁠|$3\lt z\lt 9.5$|⁠) lies above the local (⁠|$z\lt 0.1$|⁠) Sloan Digital Sky Survey (SDSS) galaxy locus,¹² to higher O32 and R23 values. This is the region numerous studies have found local ‘high-redshift analogue’ galaxies to occupy (e.g. Green peas and Blueberries, Yang et al. 2017a, b; Cameron et al. 2023). We find a trend that galaxies with larger [O iii]|$\lambda 5007$| EWs tend to sit further from the local SDSS galaxy distribution, and have much higher O32.

We observe no clear turnover in the R23 diagnostic for the highest EWs, suggesting that while these galaxies are all metal-poor (subsolar), they have not entered the extremely metal-poor branch of the R23 diagnostic (Curti et al. 2020) where the oxygen lines lose strength relative to the hydrogen recombination lines due to the extremely low range of oxygen abundances. It is therefore unsurprising that the R23 diagnostic does not reveal a clear trend with EW as we are probing the turnover region of the R23 diagnostic and our resulting measurements of metallicity will be poorly constrained by R23.

More robust oxygen abundance measurements for our full sample have been presented in Curti et al. (2024). These measurements rely on multiple strong-line calibrations utilizing the rest-optical emission lines, and metallicity estimates from spectral model fitting. For 10 objects in our full sample (including 1 EELG, ID 2651) metallicity constraints could not be provided due to too few key lines being detected (typically only [O iii] and H |$\alpha$| were detected in these cases). The Curti et al. (2024) metallicities for our EELGs lie in the range |$7.3 \lt 12 + \log _{10}({\rm O/H}) \,\lt\, 8$| (where solar metallicity is 8.69 in these units, Asplund et al. 2009), suggesting that while these systems are all subsolar (all |$\lt 0.2Z_\odot$|⁠) none show extremely low metallicity (below |$\lt 0.04Z_\odot$|⁠). The metallicities for our sample are shown in the lower panel of Fig. 7 against the [O iii]|$\lambda 5007$| EW.

Within our sample we find that galaxies with higher [O iii]|$\lambda 5007$| EWs have typically lower metallicities, albeit with some degree of scatter. This trend is additionally seen in the stacked subsamples in the lower panel of Fig. 7, with the metallicity measured to be lower in the stacks with the highest [O iii]|$\lambda 5007$| EW.

This trend is expected since [O iii]|$\lambda 5007$| is a collisionally excited line and the luminosity increases with the electron temperature, related to the metallicity (due to the limitation of the metal line cooling mechanism), allowing the EW to increase at subsolar oxygen abundance (up to a point when there are too few oxygen atoms – i.e. we expect monotonic behaviour at very low metallicity where the EW will decrease). We note that the highest EW EELGs in our sample have metallicities of |$\sim 0.1Z_\odot$|⁠, with maximum EW found in galaxies with measured metallicities above and below this value never reaching as high. This maximum-EW turnover with metallicity is consistent with photoionization models, such as in Tang et al. (2021b), that predict only marginal increase in the [O iii] + H |$\beta$| EW between 0.3 and 0.1|$Z_\odot$|⁠, with a turnover around |$0.1Z_\odot$|⁠, after which the EW begins to decrease with further decrease in metallicity. This turnover may explain the few extremely metal-poor galaxies with low EWs in Fig. 7.

4.3 UV properties from spectral modelling

Going beyond the rest-optical emission lines, the rest-frame UV properties of a galaxy are sensitive to recent star formation and the presence of hot O/B type stars whose emission peaks at shorter wavelengths. The UV continuum is sensitive to star formation on 100 Myr time-scales, whereas nebular lines are sensitive to shorter time-scales of 10Myr (see Fig. 17). To examine the non-ionizing UV continuum we model the rest-frame UV spectrum for each target in our full sample, using our NIRSpec prism spectra (which have been corrected for slit losses).

4.3.1 |$\beta$| slopes and M|$\rm _{UV}$|

We determine the absolute magnitude in the rest-frame UV (M|$\rm _{UV}$| at 1500 Å) for our sample, which we measure over a rest-frame 50 Å-wide boxcar filter centred on 1500 Å in our prism spectra. In Appendix  A, we also consider the photometry from the NIRCam filter best approximating the rest-UV (uncontaminated by the Lyman-break) using the Kron magnitude to get an indication of the total flux, and we find in most cases good agreement in M |$\rm _{UV}$|⁠. We also determine the UV spectral slope (⁠|$\beta$|⁠, where |$f_{\lambda }\propto \lambda ^{\beta }$|⁠), from our prism spectra over a rest-frame range 1340 to 2400 Å  (we mask the regions contaminated by the C iv|$\lambda 1550$|⁠, He ii|$\lambda 1640$|⁠, [O iii]|$\lambda 1660$|⁠, and C iii]|$\lambda 1990$| emission lines). We note that, as described in Bunker et al. (2023a), we have employed wavelength dependent slit loss corrections, accounting for the differing PSF and placement of the galaxy in the MSA shutter, ensuring we are not biasing UV slope measurements to bluer values with preferentially more flux lost at longer wavelengths. The uncertainties on the UV slope are measured via a Monte Carlo–based approach, whereby the power-law fit is performed 500 times with the uncertainties on each pixel randomly sampled from a Gaussian distribution. The resulting standard deviation of the UV slope measurement from 500 iterations is used as the uncertainty, and we note that this error reflects the statistical uncertainty but does not include systematic uncertainty. We note that for 24 galaxies in our full sample at |$z\lt 4$|⁠, the rest-frame 1500 Å  moves to the extreme blue end of our spectral coverage where the sensitivity reduces. In these cases, we move our |$\beta$| fitting spectral range to slightly longer rest-frame wavelengths [2000–3500 Å] and use this fit to estimate M|$\rm _{UV}$| at 1500 Å. For these objects, we additionally measure M|$\rm _{UV}$| from the JADES aperture photometry¹³ of the legacy HST/ACS F775W imaging (see Rieke et al. 2023b), which covers the rest-frame 1500 Å at |$z\sim 3.5-4$|⁠. There were two galaxies at |$z\gt 4$| (ID: 10009693 and 10009320) with a S/N|$\lt 3$| in the integrated spectrum over 1400–1600Å (rest-frame), and as in Section 4.1 we use instead the NIRCam photometry in the filter best approximating 1500 Å (rest) to determine M|$\rm _{UV}$|⁠.

The M|$\rm _{UV}$| and |$\beta$| slope measurements for our full sample are shown against the [O iii]|$\lambda 5007$| EW in Figs 8 and 9. Our full sample covers a broad M|$\rm _{UV}$| range, between |$-16 \gt $|M|$\rm _{UV} \gt -20$|⁠. As noted in Section 2.2, at |$z\gt 5.7$| the HST-Deep spectroscopy is essentially a rest-UV-based selection (a proxy for star formation rate, hereafter SFR) with a magnitude cut of |$AB\sim 29$| mag in the filter covering wavelengths just redwards the Ly |$\alpha$| break, but uncontaminated by the break. We plot in Fig. 8 the corresponding M|$\rm _{UV}$| selection at |$z\gt 5.7$| as a dashed line, and as can be seen our spectroscopic measurements of the UV luminosity broadly agree with the range expected from initial pre-selection of spectroscopic targets. Our selection is sensitive to all galaxies brighter than M|$\rm _{UV}=-18$| over our entire redshift range for the EELG selection (out to |$z=9.5$|⁠), and we are sensitive to galaxies with fainter UV luminosities at lower redshifts (as we impose a cut in apparent magnitude).

At lower redshifts, the spectroscopic target selection was based on redder rest-frame wavelengths (a better proxy for stellar mass rather than SFR). This means that the full sample does not have an M|$\rm _{UV}$| cut, although in Fig. 8 we show where such a cut would have fallen adopting the same rest-UV apparent magnitude limit as in the |$z\gt 5.7$| selection (dotted line). About 12 per cent of the |$z\lt 5.7$| sample (7/57) are UV faint (⁠|$m_{\rm {UV}}\gt 29$|⁠), and can be excluded to make a fairer comparison with the higher redshifts. We note that applying an apparent magnitude UV cut will not recover galaxies with extremely blue UV spectral slopes |$\beta \ll -2$|⁠, which would have been selected in a UV selection but may have been lost in our longer-wavelength selection.

We now consider the beta slopes and present these in Fig. 9. The full sample covers a broad range of |$\beta$| slopes (⁠|$-0.5 \gt \beta \gt -3$|⁠). The average observed |$\beta$| slope of the 36 EELGs is |$\beta =-2.2$| with a standard deviation of |$\sigma =0.3$|⁠, so the average is |$\beta _{\rm EELG}=-2.2\pm 0.1$| (where we quote the standard error on the mean). This UV spectral slope is steeper than that for the 49 non-EELGs, which have |$\beta _{\rm nonEELG}=-1.9\pm 0.1$| (again, quoting the standard error on the mean from the standard deviation of |$\sigma =0.5$|⁠). Hence, the EELGs have marginally more blue UV spectral slopes at the |$\approx 3\sigma$| level.

The UV spectral slope can depend on the level of dust attenuation and the SFH. Nebular continuum emission (e.g. Topping et al. 2022; Cameron et al. 2024) may also affect the measured UV slope, particularly with the very high EW lines in EELGs, making UV slopes redder than expected from very young ionizing stars alone. We correct our |$\beta$| slope measurements for dust reddening using the measured Balmer decrement (see Section 4.2.1). In Appendix  D, we present a simple relation for the dust correction to the observed |$\beta$| slope as a function of dust attenuation (⁠|$A_{1600}$|⁠). We apply this relation to our sample to correct the UV slope |$\beta$| for dust and we compare the trends with the [O iii]|$\lambda 5007$| EW in the central panel of Fig. 9. As noted in Section 4.2.1, we are unable to constrain the dust attenuation for four galaxies and these objects are excluded from the dust corrected panels in the figure.

We find the EELG sample occupies a smaller and slightly bluer range of de-reddened |$\beta$| slopes than the non-EELGs. One non-EELG (ID:10013704) with a reddening corrected |$\beta \lt -3.5$| has been identified as a potential AGN (Maiolino et al. 2023) which may drive the observed blue slope |$\beta =-2.14\pm 0.01$| and large measured Balmer decrement.

Interestingly, we find many objects with very blue slopes (⁠|$\beta \lt -2.5$|⁠) that are not classified as EELGs. These objects are also expected to host a young stellar population and therefore produce copious quantities of ionizing photons; however, they show low [O iii]|$\lambda 5007$| EWs (⁠|$\lt 500$|Å) and also low EWs of the Balmer lines |$H\alpha$| and H|$\beta$|⁠. In the right panel of Fig. 9 we additionally correct our [O iii]|$\lambda 5007$| EWs for dust reddening to understand whether this population of blue |$\beta$| slopes at low EW may be driven by differential dust attenuation of the emission line relative to the continuum (e.g. Calzetti et al. 2000). We find that while more galaxies exhibit dust corrected [O iii]|$\lambda 5007$| EWs |$\gt 750$| Å, there still remains a subsample of galaxies with blue slopes and low EWs. Galaxies with young stellar populations but low EWs have been previously identified at high redshift (⁠|$z\sim 5{\!-\!}7$|⁠; e.g. Endsley et al. 2023), where the lack of high EWs has been attributed to either extremely low oxygen abundance or a high escape fraction of ionizing radiation. Extremely low metallicities and high |$f_\mathrm{esc}$| can both reduce the luminosity of the [O iii]|$\lambda 5007$| line. However, we do not see strong H |$\alpha$| emission in these galaxies with low [O iii]|$\lambda 5007$| EWs and blue |$\beta$| slopes, which disfavours the extremely low metallicity scenario. Hence, these particular galaxies may be ‘Lyman Continuum leakers’ with non-negligible escape fractions of ionizing photons.

To further examine the effect of dust and SFH on the observed |$\beta$| slope, we compare the |$\beta$| slope with the Balmer decrement for both the EELG and non-EELG samples in Fig. 10. We determine best-fitting relations

for the EELG sample and

for the non-EELG sample. For the EELG sample (shown in blue in Fig. 10), we find an anticorrelation between the observed |$\beta$| slope and the Balmer decrement, which closely tracks that expected for an intrinsic (i.e. de-reddened) UV spectral slope of |$\beta =-2.3$| (Wilkins et al. 2011, shown in black in Fig. 10). However, the non-EELG population displays a higher scatter and typically redder observed |$\beta$| slopes that do no exhibit such a tight anticorrelation with the Balmer decrement. We find this apparent in the measure Spearman rank coefficients for the two subsamples, with the EELG sample showing a strong anti-correlation with a Spearman rank coefficient of −0.74 (p-value of 2e−4) while the non-EELG sample correlation is weaker with a Spearman rank coefficient of only −0.23 (p-value of 0.16). This is consistent with the picture where galaxies with modest emission line EWs have a variety of SFHs leading to a range of intrinsic beta slopes, whereas to have extreme EWs a galaxy must be dominated by recent intense star formation leading to a blue intrinsic |$\beta$| slope.

This is seen clearly in the stacks in Fig. 10, where the high-redshift EELG stack (shown in red) exhibits the steepest UV spectral slope of |$\beta =-2.36\pm 0.03$| while the shallowest slope is found in the low-redshift non-EELG stack (in green) with a |$\beta =-1.87\pm 0.02$|⁠. This trend reflects the older stellar ages (and/or potentially more dust attenuation) in the lower redshift non-EELG sample.

In the rest-optical, the non-EELG low-redshift stack is the only one to show evidence for a Balmer break (see Fig. 3), with an excess of continuum flux above 3645 Å, confirming the change in galaxy properties between the samples (see also Section 4.4). This is not unexpected as the low-redshift sample comes from a mass-limited survey which will include less actively star-forming galaxies. For this low redshift sample, the Universe is older than the high-redshift sample, so there has been more opportunity for older stellar populations to form and contribute flux at rest-optical wavelengths, reducing the EW of the emission lines.

4.3.2 Ionization efficiency

We determine the ionization efficiency (⁠|$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$|⁠, the hydrogen ionizing photon production rate per unit UV luminosity at 1500 Å, corrected for dust attenuation) of our full sample using the dust corrected H |$\alpha$| flux and M |$\rm _{UV}$| for each galaxy. We follow the method of Tang et al. (2019) to determine |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$|⁠. Here, we define

where we determine the UV luminosity from the measured M|$\rm _{UV}$|⁠, and determine the hydrogen ionizing photon production rate N(H|$^0$|⁠)[s|$^{-1}$|] from the measured H|$\alpha$| luminosity L(H|$\alpha$|⁠)[erg s|$^{-1}$|] using the linear scaling factor (1.36|$\times 10^{-12}$| erg) from Osterbrock & Ferland (2006). We split our sample by redshift (⁠|$z=5.7$|⁠) and plot ionization efficiency as a function of [O iii]|$\lambda 5007$| EW, along with the best-fitting relations, on Fig. 11. We over plot the |$z\sim 2$| best fit from Tang et al. (2019) and we find a comparable relation with larger |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$| being found at higher [O iii]|$\lambda 5007$| EWs. This correlation should not be surprising as the calculation of |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$| and EW are closely related quantities: |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$| comes from the ratio of H |$\alpha$| flux to the rest-UV continuum and EW is related to ratio of the [O iii]|$\lambda 5007$| flux to the rest-optical continuum. We note that the best-fitting relation for our lower redshift sample lies at higher |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$| than our high-redshift sample. If we consider |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$| simply as the total number of photons with energies higher than 13.6 eV, then at higher metallicities, these photons will be absorbed by many other metals such as C, N, Ne, Ar, etc., in addition to O, leading to a weaker increase in [O iii]|$\lambda 5007$| EW with increasing |$\xi _\mathrm{ion}^\mathrm{HII}$|⁠. Likewise, Curti et al. (2017) find that the [O iii]/H |$\beta$| ratio is strongly dependent on metallicity. Therefore, evolution of the oxygen abundance with redshift would mean that the expected [O iii]|$\lambda 5007$| EW would vary for the same |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$| between our two subsamples. We determine a best-fitting relation

for the high-redshift sample and

for the low-redshift sample.

We find that the average |$\xi _\mathrm{ion}^\mathrm{H\,{\small II}}$| in the EELG sample (high- and low-redshift combined) is |$\log _{10}(\xi _\mathrm{ion}^\mathrm{H\,{\small II}}/{\rm erg^{-1}Hz})$| of |$25.5\pm 0.2$|⁠, which is higher (although within |$1\sigma$|⁠) than the average of |$25.2\pm 0.3$| in the non-EELG sample, consistent with the harder ionizing spectra (bluer UV spectral slopes) of the EELGs (in line with results from Simmonds et al. 2024 using the JADES photometric data set). We find this pattern is repeated in Fig. 11 in the stacks of EELGs (both high and low redshift) which are either consistent or have higher ionization efficiencies than the non-EELG stacks.

This result reaffirms that EELGs, including those at high redshift, are efficient ionizing photon producers, an important characteristic as a potential source of the UV-ionizing background during the EoR.

4.4 Inferred galaxy properties from spectral fitting

In addition to studying the emission line diagnostics and rest-UV spectral properties, we can also model the full spectrum of each galaxy using spectral evolution synthesis codes to infer the stellar mass and the SFH. In this section we are interested in what characteristics of the modelled SFHs can be associated with a galaxy that has entered an EELG phase.

4.4.1 BEAGLE SED fitting

We model the full prism spectrum (⁠|$0.6\lesssim \lambda _{\rm obs}/\mu{\rm m}\lesssim 5.2$|⁠) for each galaxy using beagle (version 0.27.10, Chevallard & Charlot 2016) to constrain the galaxies’ physical properties and SFHs. Briefly, for each galaxy the modelled redshift is set to the spectroscopic value, and the maximum formation redshift is set to |$z=30$| (corresponding to an age of the Universe of 100 Myr for the adopted cosmology). We utilise the photoionization models from Gutkin, Charlot & Bruzual (2016), which are computed with Cloudy (Ferland et al. 2017) and are based on a more recent version of the Bruzual & Charlot (2003) stellar population synthesis models (for more details, see Chevallard in preparation). We adopt the Charlot & Fall (2000) two-component dust model, fixing to 0.4 the fraction of dust attenuation arising in stellar birth clouds. The specific allowed ranges for the variable beagle parameters are laid out in Table 3.

We employ a two-component SFH model. This choice reflects that EELGs may be experiencing an upturn in SFR or a short-lived starburst. Our model SFH is comprised of a delayed exponential component¹⁴ and a constant SFR burst of variable duration, where the two components are truncated, i.e. they are not simultaneous. The choice of a recent burst component is related to the nebular emission properties (continuum and line), which are powered by massive O-type stars which typically move off the main sequence before they reach an age of |$\sim 10$| Myr, and the disruption time-scale of stellar birth clouds (Murray, Quataert & Thompson 2010; Murray 2011). For the single stellar models used here¹⁵, more than 99  per cent of the ionizing flux is emitted by stars younger than 10 Myr.

The full beagle data products used in this paper will be presented in Chevallard (in preparation). We note that when observing the rest-optical caution should be taken as bright emission from young stars may hide an older and fainter/redder stellar population.

We determine galaxy properties for 84/85 galaxies in our full sample, with only one galaxy¹⁶ unable to be fit by beagle due to low signal-to-noise ratio in the continuum. We present the full range of SFH property constraints in Table E1 in Appendix  E. In Fig. 12, we present the beagle best-fitting total stellar mass (including stars and stellar remnants), and the mass-weighted- and luminosity-weighted (v-band) ages of the stellar population as a function of the [O iii]|$\lambda 5007$| EW.

We highlight the region identifying EELGs (shaded region in Fig. 12) and find that the EELG sample on average exhibits lower masses (typically |$\lesssim 10^8\,{\rm M}_\odot$|⁠) and lower ages (⁠|$\lesssim 5{\!-\!}10$| Myr, weighted either by mass or luminosity) than the non-EELG sample. While both EELG and non-EELG subsamples cover a broad range of masses and luminosities, the dynamic range of ages in the EELG subsample is considerably smaller than in the non-EELG sample.

To consider this further we present in Fig. 13 the sSFR for our sample as a function of EW, utilizing the SFR average over the most recent 3 Myr. Galaxies with high EWs typically have high sSFR, with EELGs having sSFR typically above 10–20 Gyr|$^{-1}$| (meaning mass-doubling times of |$\lesssim 100$| Myr). This relation is fundamental as the high sSFR is needed in EELGs for the nebular lines to dominate over the rest-optical stellar continuum (to produce a large EW), which in turn drives the mass-weighted-age to low values as seen in Fig. 12.

On the other hand, the larger dynamic range among non-EELGs reflects the greater variety of SFHs in this subsample. Likewise, we find that this pattern is repeated in the luminosity-weighted-ages, with EELGs showing ages below |$\sim 5$| Myr.¹⁷ Here, the young, massive stars necessary to produce the ionizing radiation that power the nebular emission lines outshine the light from older stellar populations in the rest-frame V band, heavily skewing the luminosity-weighted ages. We find greater discrimination between EELGs and non-EELGs in the luminosity-weighted ages than in the mass-weighted figure (due to the high luminosity of the highest mass and shortest lived O- and B-type stars), and we note that all bar one non-EELG is consistent within |$1\sigma$| with an age |$\gt 5$| Myr and all EELGs are consistent within |$1\sigma$| with ages younger than this.

4.4.2 SFH and SFR duty cycle

We can use the SFH of each galaxy, returned by beagle, to assess their duty cycle of star formation (i.e. what fraction of time a galaxy spends in an active star-forming phase), which is a measure of the ‘burstiness’ of star formation. We are interested in whether EELGs show an upturn or burst in their SFR in the last few million years, the period corresponding to the short lifetime of the massive O-type stars which provide the source of the ionizing radiation which power the nebular emission. Some of this information is encapsulated in the sSFR (Fig. 13), but we caution that the stellar masses returned in the fitting may miss very old stellar populations which do not significantly contribute to the observed light. Hence, to analyse the recent variability (burstiness) of the star formation, that the observations are more sensitive to than the overall stellar mass, in Fig. 14 we show the ratio of the intensity of the very recent SFR¹⁸ to the average SFR over a longer time-scale preceding this (extending to 100 Myr¹⁹), and compare this with the equivalent width of [O iii]|$\lambda 5007$|⁠. We adopt a comparison period of 100 Myr because this is the time-scale over which the rest-frame UV continuum can be used as an SFR indicator (Section 4.3), and hence if star formation is intrinsically bursty on shorter time-scales than this it will be revealed through a comparison of SFR|$_{100}$| to the SFR inferred from the emission lines (which respond to star formation within |$\lt 10$| Myr).

The majority of EELGs have recent SFRs which are more than twice the longer-term average (over the past |$\sim$|100 Myr), with many showing a factor of 10 increase in SFR, implying that a recent sharp upturn in the SFH is required to produce an EELG. In these panels we draw on a horizontal line at the SFR|$_{\rm short-term}/$|SFR|$_{\rm long-term} = 2$|⁠, this roughly splits the EELG and non-EELG samples albeit with some scatter. We remind the reader that the 750 Å [O iii]|$\lambda 5007$| EW is arbitrary and we do not necessarily expect two fully distinct samples to emerge (i.e. the population is continuous rather than bimodal).

As can be seen from Fig. 14, galaxies with high short-term/long-term SFR ratios are typically EELGs, but we now briefly discuss those with high ratios and lower EWs. In the lower panels of Fig. 14, we repeat the 3 Myr/3–100 Myr-average SFR ratio and colour code by the total mass (left panel) and the luminosity-weighted age (right panel) of the galaxies. Here, we notice that the non-EELGs with SFR ratios above our horizontal line are low mass and have young luminosity-weighted-ages, yet exhibit low EWs. These highly bursty but low-EW outliers could be due to unaccounted escape of ionizing photons (which is not included in the BEAGLE fits, where |$f_{\rm esc}$| is set to 0), or alternatively we are catching them very shortly after an EELG phase (⁠|$\sim 10$| Myr) when the emission lines are dying away.

We note considerable scatter around these relations, some of which results from the signal-to-noise ratio of the observations but also potentially the parametric form of the SFH we impose for the beagle fitting which may not be fully appropriate for a more complicated actual SFH. The scatter in Fig. 14 indicates that it may not be the SFH alone which governs the EW of the rest optical lines, but other factors such as dust attenuation and metallicity may play a role.

The high SFR|$_{\rm {short-term}}$|/SFR|$_{\rm {long-term}}$| ratio we see in our EELG sample (with a median value of 4.2 compared to 0.8 for non-EELGs for the 3 Myr short-term time-scale) suggests that we are catching the EELG galaxies in an ‘on phase’ of active star formation which is short compared to the ‘off phase’ of less-active star formation (where the galaxy will not be an EELG in this part of the duty cycle). The high excess of SFR at recent times compared to the long-term average suggests we are witnessing a short duration burst in EELGs. This is supported by the high sSFR we observe in our EELG sample in Fig. 13.

To examine the SFH as a function of EW in our sample, we stack the SFHs over the last 100 Myr of each galaxy into four [O iii]|$\lambda 5007$| EW bins. To produce each stack we take the beagle SFH with a 1 Myr time-step and normalize each individual SFH so the mass formed in the last 100 Myr is equal, then for each EW bin we stack the subsets performing a median average. We re-normalize the stacks to have equal mass and present these stacks in Fig. 15. We find that the stacks of the higher the EW bins have a greater proportion of star formation occurring in the most recent 10Myr, as expected for the time period over which the nebular emission occurs. For the bins of lower EW, the stars are on average older with star formation in the last 10 Myr lower compared to activity at earlier times. Interestingly, we note that the moderate EW stacks (⁠|$250\,\lt \,{\rm EW}\,\lt\, 500$| and |$500\,\lt \,{\rm EW}\,\lt 750$|⁠) exhibit a peak of star formation just before a lookback time of 10Myr. Such a bump in an SFH at this lookback time would generate a blue UV spectral slope |$\beta$| yet yield a low EW of nebular emission lines. This SFH may reflect a scenario where star formation has recently been quenched in these galaxies. Alternatively perhaps, these galaxies may be actively star forming at this current epoch but may have high escape of ionizing photons and hence have lower EW.

Finally, using our beagle SFHs we can measure a duty cycle for each galaxy. Here, we define the duty cycle as the fraction of the time (over the last 100 Myr) where the SFR of the galaxy exceeds our approximate threshold for an EELG (when the measured SFR exceeds the long-term average SFR by a factor of 3). The fraction of the time spent in this ‘active phase’ is shown in Fig. 16. What is apparent is that EELGs tend to have short duty cycles, spending an average of 5–10 per cent of the last 100 Myr in an active phase. The non-EELGs show a marginally more extended distribution of duty cycles, although we cautionthat the parametric form of the SFH we adopt (a recent burst, preceded by a delayed exponential) may be less appropriate for systems which have undergone a recent down-turn in star formation (due perhaps to quenching). Hence, the duty cycle inferred for the non-EELGs is less reliable. It is likely that a subset of non-EELGs in our sample might actually be observed to be EELG if we were to catch them at the right time in their SFH, but the more extended distribution duty cycles might imply there may also be a population there that has a much smoother variation of SFR (although again we caution that the non-EELG duty cycle determination may be subject to uncertainty due to the parametric form of the SFH adopted). However, overall our results indicate that star formation at high redshift has a significant contribution from bursty systems.

4.5 EELGs as bursty star formers: a repeating burst model

The BEAGLE fits we performed involved an ongoing starburst of variable duration following a delayed exponential SFH. At high-redshift bursty star formation plays an increasing role (Looser et al. 2023; Sun et al. 2023a; Dressler et al. 2024) to assess whether EELGs could be such bursty galaxies caught in an on-phase we introduce a toy model for the SFH, whereby there are repeated identical ‘top hat’ bursts of star formation that are regularly spaced in time. We vary the duration and frequency of these bursts, and consider the equivalent width of the [O iii]|$\lambda 5007$| emission as a function of time, assessing for each burst duration and frequency the fraction of the time that our EELG criterion of |${\rm EW}\gt 750$| Å would be met.

In Fig. 17, we show an example of the evolution of the EW([O iii]) for a burst of 10 Myr duration which repeats every 40 Myr (i.e. a star formation duty cycle of 25 per cent). We adopt a 10 per cent solar metallicity and |$\log _{10}(U)=-2$| fiducial model (aligning closely to the properties of our EELG sample, see Section 4.2.2, the lower panel of Fig. 7 and the lower left panel of Fig. 6). We find that galaxies following this SFH continue to meet the EELG threshold even after several successive bursts. However, the maximum EW achieved in each successive burst decreases due to the older stellar population starting to build up (as reflected in the sSFR in Fig. 17) and significantly contributing to the continuum around 5007 Å, diluting the EW of the [O iii].

We consider different burst durations and frequencies up to an age of 500 Myr (which is about the age of a galaxy at the median redshift of our EELG sample, |$z\approx 5.7$|⁠, which formed at |$z\approx 10$|⁠), and in the top panel of Fig. 18 we show the fraction of the repeating bursts up to 500 Myr which reach an EELG phase. The bottom panel of Fig. 18 shows the total fraction of time over the 500 Myr of repeating bursts that a galaxy would be observed as an EELG.

For a star-forming duty cycle (DC = T|$_{\rm duration}$|/T|$_{\rm repeat}$|⁠) of less than 0.1, every burst of star formation results in a short lived EELG phase. For higher DCs than this, the accumulation of low-mass stars (that contribute to the stellar continuum at 5007 Å  but not to the line emission) prohibits the galaxy from reaching high EW in subsequent bursts, meaning that the galaxy no longer can enter an EELG phase after several cycles of star formation.

The region of parameter space in our toy model that delivers the greatest time spent as an EELG over 500 Myr has a DC|$\sim 0.15$|⁠. Towards higher DC, the EELG phase is prohibited in later bursts (as a significant older population of stars has built up, diluting the line emission), and at lower DC the galaxy spends such little time in a burst that the total accumulation of time as an EELG is short. Therefore, we are most likely to catch an EELG with a duty cycle of 10–20  per cent (model dependent).

Although the toy-model of repeating bursts presented here is not exactly the same as the SFH adopted for the BEAGLE SED fitting, both approaches indicate that star formation duty cycles of |$\approx 10{\!-\!}20$| per cent are likely to yield galaxies with an EELG phase.

Taken at face value, it seems curious that EELGs comprise 42 per cent of our spectroscopic sample (rising to 62  per cent for the |$z\gt 5.7$| subsample), while from the SFH we infer a short ‘on-phase’ where a young galaxy will spend only |$\sim 10$| per cent of its time as an EELG. This can be reconciled when we consider that our high-redshift sample is selected in the rest-frame UV. In our toy model, we assume top hat bursts with zero star formation in the off-phase. As can be seen from Fig. 17, after a burst of duration |$\sim 10$| Myr the UV flux at 1500 Å  our preferred model will drop by a factor of 5 over the first 10 Myr and by a factor 10 after 20 Myr. With such UV dimming it may not enter our UV-selected sample, as seen in Fig. 8. This implies that galaxies in an EELG phase can be boosted into our sample during an on-duty phase, but would be missed during the off-duty phase. This also explains why the EELGs have lower stellar masses on average than the non-EELGs (which would be detectable in our spectroscopy even if they have much lower sSFRs than the EELGs). Lower mass galaxies are expected to be much more numerous (the stellar mass function may be Schechter in shape, e.g. Weibel et al. 2024), but only those in a transitory phase of intense star formation will be bright enough to appear in our spectroscopic sample.

5 DISCUSSION

Having identified a sample of EELGs and examined how their galaxy properties (line diagnostics, UV-properties, SFHs) compare to non-EELGs, we can now assess what conditions are required for a galaxy to enter an EELG phase. We will discuss what physical characteristics drive a galaxy to have extreme equivalent width line emission in Section 5.1. We will then consider these conditions in the context of the high rate of Ly |$\alpha$| emission line detections in our EELG sample in Section 5.2. Finally, we will discuss in Section 5.3 how the fraction of SFGs in an EELG phase changes with redshift and what evolution in galaxy properties may drive this.

5.1 What causes SFGs to enter an EELG phase?

From our sample of EELGs, we have determined consistent characteristics that all EELGs exhibit, and what observational features are distinct from non-EELGs. We determine that galaxies can enter an EELG phase if they meet the following conditions:

• They contain a young stellar population, with the highest EWs coming from stellar populations caught in the first few Myr. This ensures the presence of massive, hot O-type stars in the stellar population which drive the rest-UV nebular emission through their ionizing UV emission. The blue rest-UV spectral slopes (⁠|$\beta$|⁠) typically seen in EELGs are consistent with this very young stellar population dominating the light, and the SED fitting also finds very low-mass-weighted ages (compared with non-EELGs).

• These young stellar populations must exhibit low levels of chemical enrichment and dust. Many EELGs are consistent with no dust – we note that differential dust reddening between the nebular emission lines (the star-forming regions) and the surrounding continuum from older stars will decrease the observed EW. In terms of metallicity, our EELG sample conforms to the trend that the highest EWs are being produced in galaxies with low oxygen abundance (subsolar, although not extremely metal poor). This reflects that the collisionally excited [O iii]|$\lambda 5007$| emission line luminosity increases with the electron temperature, which increases at lower oxygen abundances due to the restriction in the metal line cooling mechanisms (in extremely metal-poor environments the lack of oxygen atoms will reduce the luminosity). The EELGs also exhibit high O32 ratios (and correspondingly high ionization parameters), typical of low-metallicity galaxies at high redshift.

• These young stellar populations must be part of a rising SFH – either a very recent burst of star formation activity, or a more extended SFH which has recently seen a rapid upturn in SFR. Typically, EELGs in our full sample exhibit a greatly elevated SFR in the recent past compared with the long term average. With EELGs typically having a short-term/long-term SFR ratio |$\gt 3$|⁠.

• The mass created in the recent burst (in the last 10 Myr and responsible for the nebular emission) must make up a significant fraction of the total mass. This is seen in the young mass- and luminosity-weighted-ages among the EELG subsample which are below 10 Myr. These young ages reflect that there is no significant older stellar population that would produce continuum in the rest-optical strong enough to dilute the EW.

From our sample, these are all necessary conditions for an EELG. We note there are galaxies that meet these conditions and yet are not EELGs. Hence, these appear to be necessary but not sufficient conditions for galaxies to be in an EELGs phase. For example, we note that galaxies meeting all of these conditions but also exhibiting very high escape fractions of ionizing photons will show low EWs, as the Lyman continuum photons leak out before generating recombination and nebular lines. Therefore it is plausible that there remains a population of ideal ionizing photon producers that are also ultra-efficient emitters (specifically Lyman-leakers) that will not be captured by an EELG selection.

We note that one EELG in our sample (3 per cent, 1/36, ID: 8083) exhibits a broad H |$\alpha$| emission line and a narrow [O iii]|$\lambda 5007$| emission line, suggesting that the emission we are observing is partially coming from a broad-line region around an AGN (see Maiolino et al. 2023). None of the other EELGs in our sample exhibit a broad line component to any recombination line. Therefore in at least one system the ionizing radiation field is not solely generated by star formation but may have a contribution from other hard ionizing sources (e.g. an AGN).

5.2 Connection to LAEs

As we reported in Section 4.2.1, the fraction of galaxies with detected Ly |$\alpha$| emission is higher in our EELG samples than in our non-EELG samples, with 53 per cent (9/17) of high-redshift EELGs determined to also be Ly |$\alpha$| emitters. This measurement is inline with Chen et al. (2023) who determined that 9/24 (38 per cent) galaxies with [O iii] + H |$\beta$| EW|$\gt 1000$| Å (roughly [O iii]|$\lambda 5007$| EW |$\gt 750$| Å following the Tang et al. 2019 relations) have Ly |$\alpha$| detections, which they use to determine that |$50\pm 11~{{\ \rm per\ cent}}$| of |$z\sim 6$| galaxies exhibit high Ly |$\alpha$| escape fraction (⁠|$f_\mathrm{esc}^{Ly\alpha }\gt 0.2$|⁠). The production of Ly |$\alpha$| photons may be expected to be higher in EELGs considering that within these active star-forming galaxies, the young stellar populations will produce ionizing photons which should yield Ly |$\alpha$| recombination line emission. However, typically Ly |$\alpha$| photons are strongly suppressed within the galaxies through resonant scattering in an optically thick environment. In the case of LAEs where Ly |$\alpha$| emission does escape, our observations are aligned with sight-lines exhibiting a lower column density of neutral hydrogen that exists within these systems and hence allows the Ly |$\alpha$| photons to escape in our direction. The high LAE fraction (53 per cent) in our high-redshift EELG sample, suggests that these EELGs have a low covering fraction of neutral gas, creating channels of lower column density, which may have been generated by a bursty SFH (which may ionise and/or expel gas, e.g. Katz et al. 2023), although we note that there must be some remaining neutral gas in the ISM to produce the line emission that we observe in our EELGs.

Furthermore, for the Ly |$\alpha$| photons that are produced in a galaxy and escape from the ISM and CGM, our observation of any Ly |$\alpha$| emission additionally depends on the transmission of these photons through the local IGM surrounding the galaxy. At high-redshifts these LAEs likely lie within large ionized bubbles, either created by their own stellar populations or from close proximity neighbours. The large ionized bubbles allow the transmission of Ly |$\alpha$| photons such that they can travel unimpeded until they are redshifted to longer wavelengths before they interact and are scattered by any neutral hydrogen in the IGM. Although a visual inspection of the nine high-redshift EELG-LAEs shows that all are compact isolated systems, the JADES survey and FRESCO data (Oesch et al. 2023) have revealed large-scale structures at various redshifts in GOODS-South, and several of our EELGs which exhibit Ly |$\alpha$| emission lie within known overdense regions at |$z=5.8{\!-\!}5.9$| (Stanway et al. 2004, 2007; Helton et al. 2024) and |$z=7.3$| (Helton et al. 2024; Endsley et al. 2024), as can been seen in Fig. 1. As noted by Witstok et al. (2024), these overdensities may lead to a highly ionized IGM and enhanced transmission of Ly |$\alpha$|⁠. However, not all of our EELG-LAEs coincide with these known overdense regions and it is unlikely that these systems created their own ionized bubbles considering how young the stellar populations are. Instead, they may indicate the presence of previously undetected overdense regions populated by galaxies below our observational sensitivity.

In our low-redshift EELG sample the fraction of EELGs with detected Ly |$\alpha$| emission is lower at 33 per cent than the high-redshift EELG subsample. This is consistent with the findings of (Stark et al. 2010; Stark, Ellis & Ouchi 2011) who determine the fraction of LBGs that are Ly |$\alpha$| emitters (⁠|$X_{Ly\alpha }$|⁠) increases with redshift from |$z=4$| to |$z=6$| (the end of the EoR). Since the IGM plays an increasing role as |$z=6$| is approached (and into the EoR), the evolution of the true Ly |$\alpha$| emission, after IGM absorption has been corrected for, may even be more dramatic.

5.3 Evolution of the EELG fraction

As we have discussed, galaxies may spend periods of their SFH in an EELG phase when they meet a set of ideal conditions. We can determine the fraction of galaxies in our survey which are in an EELG phase (hereafter, the EELG fraction) and look at how this fraction may change with redshift.

Our full JWST/NIRSpec sample allows us to measure the EELG fraction both in the EoR and at lower redshifts (⁠|$3\lt z\lt 6$|⁠). In the left panel of Fig. 19, we show the fraction of galaxies in our full sample that exhibit an [O iii]|$\lambda 5007$| EW |$\gt 750$| Å  split into three roughly equal size redshift bins (our high-redshift sample at |$z\gt 5.7$| becomes one bin, and we divide our low-redshift sample equally into two bins at |$3\lt z\lt 4.1$| and |$4.1\lt z\lt 5.7$|⁠). We note that at |$z\gt 5.7$| we have a UV selected sample and we create a similar cut for the two lower redshift bins by removing galaxies fainter than an apparent UV AB magnitude |$= 29$| (see sample selection in Section 2.2 and the M|$_{\rm {UV}}$| distribution of our sample in Fig. 8). To compute the uncertainty on our EELG fractions, we randomly vary the EW of each individual galaxy according to their measurement uncertainty and determine the fraction that have an EW above our EELG threshold, we repeat this |$n=1000$| times and report the standard deviation on the resultant EELG fraction distribution as the fraction uncertainty.

To make a fair comparison over a range of redshifts, we include studies based on a similar UV selection. Prior to JWST, measurements of the EELG fraction in rest-frame UV selected surveys had been made at |$z\sim 2$| using HST/grism spectroscopy (Boyett et al. 2022a, for M|$\rm _{UV}\lt -19$| mag) and at |$z\sim 7$| using Spitzer/IRAC photometric flux excess (Endsley et al. 2021). Comparing these two samples, Boyett et al. (2022a) determined that the EELG fraction increased with redshift by a factor of 10, from |$\sim 4~{{\ \rm per\ cent}}$| at |$z\sim 2$| to |$\sim 40~{{\ \rm per\ cent}}$| at |$z\sim 7$|⁠. We overlay the Boyett et al. (2022a) relation for the redshift evolution of the EELG fraction in Fig. 19. We additionally overlay two recent photometric [O iii] + H|$\beta$| EW distribution studies by Endsley et al. (2023, 2024) at |$6\,\lt\, z\,\lt\, 9$|⁠, and here we determine the EELG fraction adopting the [O iii]-H|$\beta$| EW relation from Tang et al. (2019). We note that the |$z\sim 6$| and |$z\sim 8$| samples from Endsley et al. (2024) have been split into a bright (M|$\rm _{UV}\lt -19.5$|⁠), faint (⁠|$-19.5 \lt $|M|$\rm _{UV} \lt -18$|⁠), and a very faint (M|$\rm _{UV}\gt -18$|⁠) subsample.

We find that the EELG fraction in our spectroscopically confirmed sample increases with redshift, from |$0.23^{+0.04}_{-0.01}$| in our lowest redshift bin (⁠|$3\lt z\lt 4.1$|⁠) to |$0.61\pm 0.04$| in our highest bin (⁠|$z\gt 5.7$|⁠), after we have applied the UV selection (see left panel in Fig. 19). We measure consistent EELG fractions (within |$1\sigma$|⁠) with and without applying our additional UV selection criteria, for both the lowest redshift and intermediate redshift bins. Comparing our results to literature values, the EELG fractions in our full sample exceed the measurements determined in photometric surveys at similar redshifts. We fit a power law |$(1+z)^P$| relation to the spectroscopic EELG fractions from our full sample and those from Boyett et al. (2022a), and determine a best fit |$P=2.6\pm 0.1$| and the expected normalization of the EELG fraction at |$z=0$| to be |$0.30~{{\ \rm per\ cent}}\pm 0.01~{{\ \rm per\ cent}}$|⁠. This relation means that half of all galaxies in a UV-magnitude selection such as ours for which we have spectroscopic redshift are in an EELG phase at |$z\gt 6$|⁠.

The larger EELG fraction in our full sample compared to high-redshift photometric surveys may be down to our initial sample selection requirement of multiple detected emission lines in our NIRSpec spectroscopy. This restriction will remove galaxies from our sample that have low EWs, where the corresponding line emission is below the sensitivity of our observations. To correct for this bias, in the right panel of Fig. 19, we plot the EELG fraction for the same bins where we now consider the complete list of JADES targets for which spectra were obtained. If we make the assumption that the 133 galaxies with expected redshifts |$z\gt 3$| but without multiple NIRSpec line detections (which were removed during our sample selection) have correct photometric redshifts, then the EELG fractions reduce. In the lowest redshift bin (⁠|$3\lt z\lt 4.1$|⁠) the EELG fraction changes from |$0.23^{+0.04}_{-0.01}$| using the sample with robust redshifts from multiple emission lines to |$0.12\pm 0.02$| using the full JADES target list. In the middle redshift bin (⁠|$4.1\lt z\lt 5.7$|⁠) the change in the fraction is from |$0.41\pm 0.03$| to |$0.31\pm 0.03$|⁠. In the highest redshift bin (⁠|$5.7 \lt z \lt 9.5$|⁠), the fraction changes from |$0.61\pm 0.04$| based on the robust redshift sample to |$0.47\pm 0.03$| using the entire JADES targets. While this assumption, that the photometric redshift is correct in the absence of a spectroscopic redshift, may not be robust, and many targeted galaxies may not have line detections for reasons beyond low EW (e.g. incorrect photometric redshifts, placement within the slit), it does demonstrate that by considering non-detections the EELG fraction is lower. In our lowest redshift bin there is better agreement with the Boyett et al. (2022a) best fit relation. However, in our highest redshift bin, we still have a larger EELG fraction than expected from the extrapolation of the Boyett et al. (2022a) relation. We also find good agreement with the |$z\sim 6$| Endsley et al. (2024) ‘bright’ sample (M|$\rm _{UV}\lt -19.5$|⁠). We again determine a best fit power law model, with a measured |$P=2.8\pm 0.1$| and the expected normalization of the EELG fraction at |$z=0$| to be |$0.18~{{\ \rm per\ cent}}\pm 0.01~{{\ \rm per\ cent}}$|⁠. This initial JADES spectroscopy from DR1 used for this paper is largely based on a selection from HST imaging, and was by design and necessity (given the short wavelengths) primarily a rest-frame UV-based selection. Hence, we have compared the EELG fraction in our UV-based selection with comparable selections in the literature, accounting for the impact on completeness introduced by only a subset of our spectroscopically targetted galaxies having multiple emission lines for a robust redshift. JADES spectra taken later in the survey benefit from a selection based on JADES JWST/NIRCam imaging in the GOODS fields (as well as HST photometry), enabling use to select a sample based on longer rest-frame wavelengths which are more sensitive to the underlying stellar mass rather than the recent SFH. In a future paper we will look at the evolution of the EELG fraction with redshift in a stellar-mass limited sample, as well as in a rest-frame UV sample significantly larger than that used in this paper (which will enable consistent UV absolute magnitude cuts to be used in many redshift bins while retaining large numbers of galaxies for robust statistics).

To consider why the observed EELG fraction within our UV-selected sample may change with redshift, we revisit the conditions that must be satisfied for a galaxy to enter an EELG phase. Many of the galaxy properties that we have based criteria on are known to evolve with redshift. At fixed stellar mass, simple galaxy models predict that the sSFR is expected to increase with redshift (e.g. Tacchella et al. 2018), driving higher EWs. The average metallicity and dust content is also lower in the galaxy population at higher redshift (e.g. Curti et al. 2024; Langeroodi et al. 2023), further boosting the [O iii] emission line luminosity. The evolution of these galaxy properties means the average observed [O iii] EW in the galaxy population will increase with redshift and hence a greater proportion of the galaxy population at high redshift will meet the conditions needed to enter an EELG phase.

Although harder to constrain through observations, the star formation duty cycle within the galaxy population is also expected to change with redshift. Evolution in the galaxy merger rate, galaxy dynamical time-scale and rate of gas accretion, all affect the fraction of time a galaxy spends in an active phase of star formation (Kimm et al. 2015; Ceverino et al. 2018; Faucher-Giguère 2018; Furlanetto & Mirocha 2022). These changes may allow galaxies to enter an active star-forming phase more frequently at high redshift, again increasing the fraction of objects observed in an on-phase as EELGs.

Interestingly, Endsley et al. (2024) measure a decrease in the EELG fraction for each of their M|$\rm _{UV}$| bins between |$z\sim 6$| and |$z\sim 8$|⁠. This suggests that the evolution turns over at very high redshift, which is not captured by our best fit power law relation. Such a turn over may be due to a further decrease in the metallicity limiting the [O iii] emission line luminosity or other changes in the galaxy population properties (e.g. |$f_\mathrm{esc}$|⁠). From this initial DR1 data release, we do not have a large enough sample at |$z\sim 8$| to place statistically significant constraints on the EELG fraction to investigate the presence of a turn-over. However, the full JADES spectroscopic survey should have sufficient size to test this, and this will be the subject of a future paper.

Finally, we note that because of the greater fraction of galaxies in an EELG phase at higher redshift, and the evidence for EELGs being productive ionizing photon producers (high |$\xi _\mathrm{ion}^\mathrm{HII}$|⁠, see Section 4.3.1), galaxies in an EELG phase therefore have the potential to play an important role in the reionization of the Universe.

6 CONCLUSIONS

Through a JWST/NIRSpec spectroscopic sample of SFGs over a broad redshift range (⁠|$3\lt z\lt 9.5$|⁠), we examine what galaxy properties cause a galaxy to enter an EELG phase. Out of our sample of 85 SFGs, 42 per cent (36/85) exhibit [O iii]|$\lambda 5007$| rest-frame equivalent widths consistent with being |$\gt 750$| Å. This subsample of EELGs shows high EWs in both the [O iii]|$\lambda 5007$| and H |$\alpha$| emission lines. They are characterized by a hard ionizing radiation field, observed through the ionization state of the ISM (with O32 |$\gtrsim 10$|⁠); a subsolar, although not extremely metal poor, oxygen abundance (⁠|$7\lt \log _{10}(\mathrm{O/H})\lt 8.5$|⁠); a high ionizing photon production efficiency (with mean |$\log _{10}(\xi _\mathrm{ion}^\mathrm{HII}/{\rm erg^{-1}Hz})=25.5\pm 0.2$|⁠); and low levels of dust attenuation (consistent with minimal or no dust, with a mean |$A_{1600}=0.4\pm 0.1$|⁠). These properties boost the [O iii]|$\lambda 5007$| line luminosity, and hence the EW. The high ionization state of the ISM and the high emission line EWs observed in EELGs is driven by the hot O-type stars in young stellar populations, which we also observe in the characteristically blue UV spectral slopes (⁠|$\beta$|⁠) of our EELG sample (with mean observed |$\beta =-2.2\pm 0.1$|⁠, where we quote the standard error on the mean). Through SED modelling of the full prism spectrum using beagle, we constrain the SFHs of our sample, determining that EELGs typically exhibit a young luminosity-weighted ages (⁠|$\lt 5$| Myr) and a recent SFR (either instantaneous, or averaged over a 3 or 5Myr lookback time) that exceeds the long-term average (measured out to a lookback time of 100Myr preceding the recent-measurement time frame) by a factor of 2, with a median excess of 4.3.

We note that these properties that EELGs exhibit are all necessary, but not sufficient conditions for a galaxy to enter an EELG phase, with several lower EW ([O iii]|$\lambda 5007$| EW < 500 Å) systems in our full sample also exhibiting a subset of these characteristics (e.g. blue |$\beta$| slopes or high short-term/long-term SFR ratios). In these cases where low EWs are found in galaxies matching part of our EELG phase criteria, we must consider the roles of dust, metallicity and the escape fraction of ionizing photons in modulating the equivalent width (as it is the UV-ionizing photons that do not leak from the galaxy that drive the nebular line emission).

Stacking the SFHs in bins of [O iii]|$\lambda 5007$| EW we determine that EELGs produce most of their mass within a lookback of 10 Myr, while lower EW stacks favour histories where the peak of star formation occurs at greater lookback times with a drop in the SFR over the last 10 Myr (the time period responsible for the ionizing photon that drive the nebula line emission).

The beagle SED-fitting suggests that SFGs only enter an EELG phase during short-lived periods of rapid star formation. We determine from beagle single-stellar population models²⁰ that the highest probability star-forming duty cycle for an EELG is 10–20 per cent (for a metallicity and ionization parameter: 0.1Z|$_\odot$|⁠, log|$_{10}$|(U) = −2.0), adopting a top hat duty cycle paramatrization. This optimum duty cycle is consistent with the fraction of the last 100 Myr when galaxies in our sample exhibit a short-term/long-term SFR excess of 3 (a threshold which defines an on-phase from our EELG subsample). This suggests that our full spectroscopic sample contains both EELGs currently in an on-phase as well as galaxies with the same duty cycle but captured in an off-phase which are not EELGs but are sufficiently bright in the rest-UV to have entered our spectroscopic selection.

These off-phases make up a greater proportion of the total SFH and are accompanied by a dimming of the UV luminosity. Therefore, SFGs may be boosted into our UV-selected sample during an on-phase as an EELG, confirming why we observe that lower mass galaxies in our sample typically exhibit higher EWs (as a similar mass galaxy in an off-phase may not satisfy our UV selection).

The bursty SFHs in our EELGs may also be responsible for the high detection rate of Ly |$\alpha$| emission, with 53 per cent (9/17) of high-redshift (⁠|$z\gt 5.7$|⁠) EELGs confirmed to be Ly |$\alpha$| emitters, compared to only 18 per cent (2/11) of high-redshift non-EELGs. The rate of LAE detection highlights that the conditions required to enter an EELG phase also promote the escape of Ly |$\alpha$| photons, through the ionization or removal of neutral gas via outflows along our sight line.

Across our full sample (⁠|$3\lt z\lt 9.5$|⁠), 42 per cent of our sample satisfy our [O iii]|$\lambda 5007$| EELG threshold. Splitting our full sample at |$z=5.7$|⁠, we find that the fraction of galaxies in an EELG phase increases with redshift, with |$61\pm 4$| per cent of our full sample at |$z\gt 5.7$| in an EELG phase. At these high redshifts, a greater fraction of galaxies are exhibiting high equivalent width line emission through the combination of the redshift evolution in the metallicity and the duty cycle, along with the relative youth of the Universe meaning that older stellar populations have had less chance to form and contribute to the rest-frame optical continuum.

This combination of the high EELG fraction during the EoR and the high ionization efficiency in EELGs makes these galaxies productive ionizing photon producers. EELGs during the EoR may therefore contribute significantly to the ionizing radiation field responsible for reionizing the Universe.

ACKNOWLEDGEMENTS

This work is based in part on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes (MAST) 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 JWST. KB is supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. AJB, AJC, JC, IEBW, AS, KB, and GCJ 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 no. 789056). ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1). SAr acknowledges support from Grant PID2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/ 10.13039/501100011033). SCa acknowledges support by European Union’s HE ERC Starting grant no. 101040227–WINGS. SAl acknowledges support from the JWST Mid-Infrared Instrument (MIRI) Science Team Lead, grant 80NSSC18K0555, from NASA Goddard Space Flight Center to the University of Arizona. DJE is supported as a Simons Investigator and by JWST/NIRCam contract to the University of Arizona, NAS5-02015. RH acknowledges funding by the Johns Hopkins University, Institute for Data Intensive Engineering and Science (IDIES). KH, BJ, EE, DPS, RE and MR acknowledge the JWST/NIRCam contract to the University of Arizona, NAS5-02015. DPS additionally acknowledges support from the National Science Foundation through the grant AST-2109066. BER acknowledges support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015. 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. RS acknowledges support from an STFC Ernest Rutherford Fellowship (ST/S004831/1). LS, TJL, and JS acknowledges support by the STFC and ERC Advanced Grant 695671 ‘QUENCH’. CS and JW acknowledge support from the Science and Technology Facilities Council (STFC), by the ERC through Advanced Grant 695671 ‘QUENCH’, and by the UKRI Frontier Research grant RISEandFALL. JW gratefully acknowledges support from the Cosmic Dawn Center through the DAWN Fellowship. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. MVM acknowledges support from the National Science Foundation via grant NSF AAG 2205519. and the Wisconsin Alumni Research Foundation via grant MSN251397. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship

DATA AVAILABILITY

The raw observational JWST NIRSpec and NIRCam data, and the corresponding reduced data products, are publicly available on MAST and via the JADES team website https://jades-survey.github.io. The line fluxes measured for our prism spectroscopy are given in Bunker et al. (2023a), when detected above a |$5\sigma$| significance. The individual EELG and beagle galaxy properties (Tables 1 and E1) are made available in a machine readable format.

Footnotes

REFERENCES

APPENDIX A: ALTERNATIVE EW MEASUREMENT ESTIMATES

When computing the EW we use the broad-band photometry to measure the flux density of the continuum. When possible our choice of filter is the one contaminated by the line in question. In this appendix, we check the consistency of our EW with three related methods: using the measured continuum flux density from the spectroscopy, taking the photometry from an adjacent filter – uncontaminated by the emission line in question; and using the line flux from the higher resolution R1000 grating.

A1 Uncontaminated broad-band filter comparison

When we do not have NIRCam imaging of the galaxy in the broad-band filter covering the observed wavelength of the emission line, we determine an estimate of the flux density at the wavelength of the line using the nearest (and where possible contamination free) available filter. Again, we remove any contribution from spectroscopically detected line emission to the selected broad-band filter, accounting for the wavelength-dependent filter transmission. In Fig. A1, we present the comparison of the [O iii]|$\lambda 5007$| EW determined using the continuum flux density as measured in the contaminated filter and the adjacent filter. We find that using the nearest available filter provides good agreement with the measurement made using the line contaminated filter – with a scatter consistent with the variation of |$\beta$| slopes in our sample.

For 10/85 galaxies in our full sample, we do not have NIRCam photometry that covers the observed wavelength of the [O iii]|$\lambda 5007$| emission line. This is typically due to the NIRSpec and NIRCam footprints not being in perfect alignment and due to subtly different footprints in the different filter bands (e.g. F444W covers a larger area than the SW filters, the Medium bands from JEMS cover a subtly different footprint to the JADES broad-band filters).

When measuring the EW using an adjacent filter, the quoted errors present the statistical uncertainty on the line flux measurement and signal to noise of the photometry. The quoted uncertainty does not reflect the systematic uncertainty in the |$\beta$| assumption.

A2 Spectroscopic comparison

We additionally compare the EWs derived from the broad-band photometry with the EWs measured directly from the prism spectrum itself. Here, we present the EWs measured using the continuum in the 1D spectra in Fig. A2, where we show all the spectral measurements where the signal to noise per pixel in the continuum is greater than 1 around the emission line.

When the spectroscopy does not have a significant detection of the continuum, the inferred EW from the spectroscopy is highly uncertain (consistent with infinite EW in the limiting case of zero continuum). We therefore adopt the continuum estimate based on the broad-band photometry, corrected for line contamination using the spectroscopic line fluxes, as this provides better constraints on the EW.

A3 prism and R1000 resolution comparison

We note that the measured emission line fluxes are a factor of |$\sim 1.08$| larger in the R1000 resolution spectroscopy compared to the prism spectroscopy (Bunker et al. 2023b). We see this effect in the measured EW in Fig. A3, with the EW determined using the R1000 spectroscopy larger than those from the prism. We find that the integrated prism spectra best recovers the reported broad-band flux densities, and so we use these measurements as part of our study.

APPENDIX B: M|$\rm _{UV}$| measurement from photometry and spectroscopy

We want to ensure that the slitloss correction applied to the NIRSpec prism spectroscopy, which uses the position of the target in the MSA shutter and assumes a point-like source (see Bunker et al. 2023a for details), performs a sufficient job at recovering the total flux of our targets. This is crucial for our estimate of M|$\rm _{UV}$|⁠. To check the M|$\rm _{UV}$| determined from the NIRSpec spectra in Section 4.3.1, we additionally measure the rest-UV magnitude from direct imaging in the filter closest to rest-frame 1500 Å  (without including Lyman-break). We use the HST/WFC3 F775W filter for galaxies below |$z\lt 4$|⁠, the JWST/NIRCam F090W filter for galaxies |$4\lt z\lt 5.6$|⁠, the JWST/NIRCam F115W filter for galaxies |$5.6\lt z\lt 7.2$| and the JWST/NIRCam F150W filter for any galaxies at |$z\gt 7.2$|⁠. The photometry for each galaxy in our sample in each filter is publicly available as part of the JADES NIRCam data release²¹ (Rieke et al. 2023b).

A variety of aperture photometry measurements are available for each target. These include a range of fixed circular aperture measurements (up to 1 arcsec in diameter), an 80 per cent flux radius aperture, as well as Kron aperture measurements (see Rieke et al. 2023b for more details). We find consistent results for the measured photometry from the Kron and the largest circular apertures (suggesting many of our targets are resolved and are extended). In Fig. B1, we present the comparison of the M|$\rm _{UV}$| measurements from the photometry and spectroscopy. We find good agreement for the majority of bright galaxies (M|$\rm _{UV}$| < −18). Fainter that this, we find a systematic offset with the UV magnitude being underestimated in the spectroscopy by |$\sim 0.5$| magnitude compared to the photometric measurement, although we note that these measurements are still consistent within the uncertainties.

APPENDIX C: UV LUMINOSITY-WEIGHTED SPECTRAL STACK

In Section 4.1, we combined the spectra of multiple galaxies in stacks to study the population properties of EELGs and non-EELGs at high and low redshift. To create the stacks we normalized each spectrum to the rest-UV (1500 Å) and then average-combined the galaxies spectra. This approach provides equal weight to all galaxies, whereas we could alternatively combine our sample by weighting each spectrum by the luminosity, effectively summing the luminosities of the galaxies in the subsamples to obtain an indication of the total spectrum produced by EELGs. This alternative approach means UV luminous galaxies are weighted more strongly in the stacks.

We repeat the non-weighted stacking in Section 4.1, now weighting the 1500Å-normalized spectrum by the UV luminosity. We present the resultant stacks in Fig. C1, which replicates Fig. 3 with the UV-weighted stacks. As expected these stacks show a higher S/N since we prioritize the UV-bright galaxies in the stacks. We find that the measured emission line and UV-spectral properties are in line with those from the non-weighted stacks. In Table C1, we provide the spectral properties for this alternative stacking method. We find consistent stack properties between the two weighting schemes.

APPENDIX D: Dust correction to UV |$\beta$| slopes

We correct our measurements of the observed rest-frame UV slope (⁠|$\beta _{obs}$|⁠, from |$f_\lambda \propto \lambda ^{\beta }$|⁠) for the dust attenuation determined from the measured Balmer decrement, adopting a Calzetti et al. (2000) dust law and R|$_v=4.05$|⁠. To aid us with this correction we first determine a simple relation between dust attenuation and the change in |$\beta$| slope.

We create a grid of stellar²² extinction (⁠|$0\lt A_{1600}\lt 1$|⁠) and observed |$\beta$| slopes (⁠|$-1.5\gt \beta _{obs}\gt -3.5$|⁠), and for each grid point we generate an observed UV slope and correct each data point obeying the Calzetti et al. (2000) law for the given attenuation. We then fit a new |$\beta$| slope (⁠|$f_\lambda \propto \lambda ^{\beta }$|⁠) to the resultant dust corrected data and we compare the observed and dust corrected |$\beta$| values. The dust attenuation and change in |$\beta$| slope is presented in Fig. D1. We fit a relation to |$A_{1600}$| as a function of |$\Delta \beta$| (given as |$\beta _{\rm obs}-\beta _{\rm corr}$|⁠) and |$\beta _{\rm obs}$|

which we employ in our data analysis in Section 4.3.1.

APPENDIX E: SPECTRAL ENERGY DISTRIBUTION RESULTS

In Section 4.4, we modelled the NIRSpec spectroscopy and NIRCam photometry using beagle to constrain the SFH of each galaxy. In Table E1, we present the resultant SFH properties.