# Extreme emission line galaxies detected in JADES JWST/NIRSpec I: inferred
galaxy properties

Kit Boyett,1,2 Andrew J. Bunker,3 Emma Curtis-Lake,4 Jacopo Chevallard,3 Alex
J. Cameron,3 Gareth C. Jones,3 Aayush Saxena,3 Stéphane Charlot,5 Mirko
Curti,6,7,8 Imaan E.B. Wallace,3 Santiago Arribas,9 Stefano Carniani,13 Chris
Willott,18 Stacey Alberts,12 Daniel J. Eisenstein,19 Kevin Hainline,12 Ryan
Hausen,20 Benjamin D. Johnson,19 Marcia Rieke,12 Brant Robertson,21 Daniel P.
Stark,12 Sandro Tacchella,7 Christina C. Williams,22 Zuyi Chen,12 Eiichi
Egami,12 Ryan Endsley,23 Isaac Laseter,24Tobias J. Looser,7 Michael V.
Maseda,24 Renske Smit,25 Joris Witstok7,8
1School of Physics, University of Melbourne, Parkville 3010, VIC, Australia
2ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
Australia
3Department of Physics, University of Oxford, Denys Wilkinson Building, Keble
Road, Oxford OX1 3RH, UK
4Centre for Astrophysics Research, Department of Physics, Astronomy and
Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK
5Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98
bis bd Arago, 75014 Paris, France
6European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching,
Germany
7Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge, CB3 0HA, UK.
8Cavendish Laboratory - Astrophysics Group, University of Cambridge, 19 JJ
Thomson Avenue, Cambridge, CB3 0HE, UK.
9Centro de Astrobiología (CAB), CSIC–INTA, Cra. de Ajalvir Km. 4, 28850-
Torrejón de Ardoz, Madrid, Spain
10Cosmic Dawn Center (DAWN), Copenhagen, Denmark
11Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200,
Copenhagen, Denmark
12Steward Observatory University of Arizona 933 N. Cherry Avenue Tucson AZ
85721, USA
13Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
14European Space Agency, European Space Astronomy Centre, Camino Bajo del
Castillo s/n, 28692 Villafranca del Castillo, Madrid, Spain
15ATG Europe for the European Space Agency, ESTEC, Noordwijk, The Netherlands
16Department of Physics and Astronomy, University College London, Gower
Street, London WC1E 6BT, UK
17Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg,
Germany
18NRC Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada
19Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden St., Cambridge
MA 02138 USA
20Department of Physics and Astronomy, The Johns Hopkins University, 3400 N.
Charles St., Baltimore, MD 21218
21Department of Astronomy and Astrophysics University of California, Santa
Cruz, 1156 High Street, Santa Cruz CA 96054, USA
22NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 North
Cherry Avenue, Tucson, AZ 85719, USA
23Department of Astronomy, University of Texas, Austin, TX 78712, USA
24Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter
St., Madison, WI 53706 USA
25Astrophysics Research Institute, Liverpool John Moores University, 146
Brownlow Hill, Liverpool L3 5RF, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Extreme emission line galaxies (EELGs) exhibit large equivalent widths (EW) in
their rest-optical emission lines ([OIII]$\lambda 5007$ or H$\alpha$ rest-
frame EW$>750$Å) which can be tied to a recent upturn in star formation rate,
due to the sensitivity of the nebular line emission and the rest-optical
continuum to young ($<10$Myr) and evolved stellar populations, respectively.
By studying a sample of 85 star forming galaxies (SFGs), spanning the redshift
and magnitude interval $3<z<9.5$ and $-16>$ M$\rm{}_{UV}>-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
mass-weighted ages ($<5$ Myr old) and the highest burst intensity (those with
the highest proportion of their total stellar mass formed in the last 10 Myr).
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\%$ in our $z>5.7$
high-redshift bin, compared to $23^{+4}_{-1}\%$ in our lowest-redshift bin
$3<z<4.1$) due to the combined evolution of metallicity, ionisation parameter
and star formation histories (SFH) with redshift. We report that the EELGs
within our sample exhibit a higher average ionisation 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 ionising photon producers. Additionally, we report that $53\%$
(9/17) of EELGs at $z>5.7$ have observed Lyman-$\alpha$ emission, potentially
lying within large ionised regions. The high detection rate of Lyman-$\alpha$
emitters in our EELG selection suggests that the physical conditions
associated with entering an EELG phase also promote the escape of
Lyman-$\alpha$ photons.

###### keywords:

galaxies: high-redshift – star formation – evolution

††pubyear: 2015††pagerange: Extreme emission line galaxies detected in JADES
JWST/NIRSpec I: inferred galaxy properties–E

##  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 [OIII]$\lambda 5007$ lines in the rest-
frame optical. The nebular emission lines are often driven by the ionising
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., 2023a), 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 et al. 2010; Izotov et al. 2011;
Brunker et al. 2020; Kumari et al. 2023), at moderate redshifts ($z\sim 2$,
van der Wel et al. 2011; Atek 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 reionisation (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 ([OIII]$\lambda 5007$ EW thresholds ranging from $\sim 100-1000$ Å
e.g., Amorín et al. 2015; van der Wel et al. 2011; Tang et al. 2019, 2021a; Du
et al. 2020). 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 James Webb Space Telescope (JWST), studies of EELGs at high redshifts
($z>4$) had been limited to photometric excess measurements, largely from
Spitzer/IRAC, as the rest-optical H$\alpha$ and [OIII]$\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<7$, [OIII]$\lambda 5007$ $z<9.5$). This allows the unambiguous
identification of EELGs and characterisation 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., 2023). The use of NIRCam medium-band imaging to measure the
flux excess between adjacent filters has also enabled EELG identification
(Withers et al., 2023; Endsley et al., 2023a; Tacchella et al., 2023). 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
[OIII]$\lambda 5007$ EW > 750Å, for instance) are observed to exhibit lower
masses, more compact morphologies, younger stellar populations, higher
ionising photon production efficiencies, higher ionisation parameters and
lower gas-phase oxygen abundances (metallicity) than typical star forming
galaxies found at the same epoch (e.g., Tang et al., 2019).

These properties make star-forming galaxies (SFGs) in an EELG phase effective
ionising photon producers. However, at moderate redshifts, their abundance has
been measured to be too low to dominate the ionising output of the star-
forming galaxy population (Boyett et al., 2022a), with the majority of the
ionising output coming from SFGs in more typical modes of star formation. The
bursty star formation histories expected at higher redshifts and lower masses
(Ceverino et al., 2018; Ma et al., 2018; Faucher-Giguère, 2018; Tacchella et
al., 2020), which may make EELGs more common than at lower redshifts, suggests
the potential for high escape fractions of ionising photons111Although we note
that galaxies with near 100$\%$ $f_{\mathrm{esc}}$ would show very low EWs,
since nebular emission lines are powered by the UV-ionising photons which do
not escape into the intergalactic medium (IGM) but are reprocessed locally.
($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 et al. 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 region (HUDF, 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 star-forming galaxies over a redshift range $3<z<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 [OIII]$\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 sub-sample 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 prep.). 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
Multi-object spectroscopy (MOS) mode (Ferruit et al., 2022), with a three
shutters nodding sequence. Three Micro-Shutter Assembly (MSA) configurations
were constructed, each with three nod positions. Targets were assigned to the
configurations following the prioritisation 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 maximises the total number of galaxies
observed whilst retaining the maximum exposure time for the highest priority
sources. The prioritisation 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 prism222The
prism has an average resolving power $\frac{\Delta\lambda}{\lambda}=R\sim
100$, but this varies with wavelength over the range $R=30-300$. (total
integrated exposure time of 100 ksec) and the $R\sim 1000$ medium
resolution333The set of three medium resolution gratings have an average
resolving power $\frac{\Delta\lambda}{\lambda}=R\sim 1000$, but this varies
with wavelength over the range $R=300-1200$. G140M/F070LP, G235M/F170LP and
G395M/F290LP gratings (each with a total integrated exposure time of 25 ksec).
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 prep),
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 ksec (see Rieke et al. 2023b; Eisenstein et al. 2023). 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 ksec and 27.8 ksec. We note that the footprint in the different
filters varies, and a sub-set of the galaxies targeted with NIRSpec have
coverage only in F444W, whilst the majority are covered by all imaging
filters.

###  2.2 Sample selection and selection bias

Figure 1: Redshift distribution of the JADES targeted galaxies (prism in
orange, R1000 grating in green) and our defined full sample (red, $z>3$
galaxies with multiple emission lines, including [OIII]$\lambda 5007$,
detected at $>5\sigma$). Galaxies in our full sample with coverage in the
prism and R1000 gratings are shown in the left and right panels, respectively.
We have an upper redshift bound to our full sample where [OIII]$\lambda 5007$
moves out of the wavelength coverage at $z\sim 9.5$. When spectroscopic
redshifts could not be determined using NIRSpec for the JADES targeted
galaxies, we rely on the photometric redshifts used to selected these
candidates (see Bunker et al. 2023a).

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 $>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>3$, where the
H$\beta$ and [OIII]$\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 $>5\sigma$, including [OIII], in the prism
spectra. These galaxies will become our sample for this study (hereafter "full
sample").

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

This sample was pre-selected for observation through a target prioritisation
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>5.7$ (these include I775, Y105 and J125-band dropouts),
used a rest-frame UV magnitude threshold longward of the Lyman-$\alpha$ break
(rest-frame UV detection band AB magnitude $\lesssim 29$) with prioritisation
based on their photometric redshift. Over a sub-set 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
wavelengths444Those galaxies in Bunker et al. (2023a) with an 8-digit ID
starting with ‘100’ were selected purely from HST imaging, while those with
shorter IDs were selected from HST and JWST/NIRCam images.. We note that, as
expected, all the high-redshift ($z>5.7$) EELGs we identify come from targets
selected to be in the photometric redshift range $5.7<z<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>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 CIII]$\lambda 1909$), we are not
biasing our magnitude limited sample towards strong line emitters, although we
consider the selection effects due to star formation history in Section 4.3.1.

At lower redshifts ($z<5.7$), priority was given to obtaining a stellar mass-
limited sample, by selecting on the longest wavelength broadband filter with
good imaging data. In the case of HST selection this was the F160W HST/WFC3
(AB $<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 $<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<z<5.7$,
$3.5<z<4.5$ and $2.5<z<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$,
up-weighted 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, 3/45, see Section 3.1.

We note that the selection on the rest-frame optical means that over specific
redshift ranges, the broadband flux may be boosted by the contribution from
strong line emission (H$\alpha$, [OIII]+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 sub-sample to
address this selection effect.

We also note that the rest-optical selection (at $z<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<6$, seen in Figure 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 $>5\sigma$ significance are reported in the publicly available
catalogue555https://archive.stsci.edu/hlsp/jades, 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 apodised 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 apodised 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 broadband 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 broadband 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 sub-set 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 broadband photometry to provide tighter
constraints on the EWs.

When computing the continuum flux density from the broadband 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 [OIII]$\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 [OIII]$\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\%$ 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 broadband 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>3$ where the H$\beta$ and [OIII] 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 can not 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., [OIII]$\lambda 5007$ or H$\alpha$) or complexes (e.g.,
[OIII]+H$\beta$) depending on the spectroscopic resolution of observations or
if relying on broadband flux-excess. Threshold values for [OIII]$\lambda 5007$
EW in literature have ranged from $\sim 100-1000$ Å (e.g., Amorín et al.,
2015; van der Wel et al., 2011; Tang et al., 2019, 2021a; Du et al., 2020).

In this paper, we will adopt [OIII]$\lambda 5007$ EW $>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 ionisation parameter and ionisation
efficiency as well as potentially high escape fractions of ionising 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 [OIII]$\lambda 5007$ in
our sample, with one additional galaxy consistent within $1\sigma$ of the
[OIII]$\lambda 5007$ threshold (ID 10001892) which also has a H$\alpha$ EW
above $750$Å and hence is an EELG. This makes up 42$\%$ (36/85) of our sample
of emission line galaxies at $z>3$. This sub-sample 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.

ID | zspec | [OIII]$\lambda 5007$ EW | H$\alpha$ EW | Ly$\alpha$ EW | MUV | $A_{1600}$ | $\beta_{\rm{obs}}$ | O32 | R23 |  $\log_{10}(\xi_{\mathrm{ion}}^{\mathrm{HII}}$/ | $\log_{10}\mathrm{(O/H)}$
---|---|---|---|---|---|---|---|---|---|---|---
|  | Å | Å | Å |  |  |  |  |  | $\rm{erg^{-1}Hz})$ |
10058975 | 9.44 | $856\pm 75$ | – | <6 (P) | $-20.4\pm 0.1$ | $0.0\pm 0.1$ | $-2.59\pm 0.01$ | $41.3\pm 8.0$ | $6.0\pm 0.2$ | – | $7.4\pm 0.1$
21842 | 7.98 | $1246\pm 37$ | – |  $24\pm 7$ (M) | $-18.7\pm 0.1$ | – | $-2.23\pm 0.01$ | – | – | – | $7.6\pm 0.1$
10013682∗ | 7.28 | $2767\pm 763$ | – |  $148\pm 20$ (M) | $-15.7\pm 1.4$ | – | $-1.85\pm 0.04$ | – | – | – | $7.7\pm 0.1$
10013905 | 7.21 | $844\pm 77$ | – | <32 (P) | $-18.8\pm 0.1$ | $0.0\pm 0.8$ | $-2.46\pm 0.01$ | $20.0\pm 3.3$ | $6.8\pm 0.4$ | – | –
4297 | 6.72 | $1971\pm 39$ | $1721\pm 88$ | <63 (P) | $-18.5\pm 0.1$ | $0.0\pm 0.3$ | $-2.29\pm 0.02$ | $12.1\pm 2.0$ | $8.9\pm 0.6$ | $25.6\pm 0.1$ | $8.0\pm 0.2$
3334 | 6.71 | $1620\pm 72$ | $1558\pm 307$ | <31 (P) | $-18.1\pm 0.1$ | $1.5\pm 0.6$ | $-2.08\pm 0.02$ | $9.2\pm 4.4$ | $14.0\pm 5.4$ | $25.2^{+0.1}_{-0.2}$ | $7.8\pm 0.1$
16625 | 6.63 | $1444\pm 29$ | $1978\pm 97$ |  $39\pm 17$ (M) | $-18.8\pm 0.1$ | $0.4\pm 0.2$ | $-2.29\pm 0.01$ | $>25.9$ | $<5.1$ | $25.6\pm 0.1$ | $7.3\pm 0.1$
18846 | 6.34 | $1068\pm 21$ | $1185\pm 24$ |  $39\pm 4$ (M) | $-20.2\pm 0.1$ | $0.0\pm 0.1$ | $-2.53\pm 0.01$ | $38.8\pm 7.3$ | $6.4\pm 0.2$ | $25.3\pm 0.1$ | $7.5\pm 0.1$
18976 | 6.33 | $839\pm 30$ | $620\pm 43$ | <20 (P) | $-18.7\pm 0.1$ | $0.0\pm 0.3$ | $-2.30\pm 0.01$ | $>16.2$ | $<6.1$ | $25.4\pm 0.1$ | $7.5\pm 0.1$
19342 | 5.98 | $2090\pm 72$ | $1072\pm 76$ |  $51\pm 22$ (M) | $-18.7\pm 0.1$ | $0.0\pm 0.2$ | $-2.51\pm 0.01$ | $>22.1$ | $<6.0$ | $25.4\pm 0.1$ | $7.5^{+0.2}_{-0.1}$
9422‡ | 5.94 | $3253\pm 65$ | $2214\pm 44$ |  $130\pm 23$ (M) | $-19.7\pm 0.1$ | $0.0\pm 0.1$ | $-2.31\pm 0.01$ | $38.7\pm 4.3$ | $7.7\pm 0.2$ | $25.6\pm 0.1$ | –
6002 | 5.94 | $870\pm 20$ | $619\pm 26$ |  $44\pm 20$ (M) | $-18.8\pm 0.1$ | $0.3\pm 0.2$ | $-2.52\pm 0.01$ | $8.9\pm 1.7$ | $8.6\pm 1.1$ | $25.2\pm 0.1$ | $7.7\pm 0.1$
19606 | 5.89 | $962\pm 21$ | $941\pm 36$ |  $84\pm 10$ (M) | $-18.6^{+0.1}_{-0.2}$ | $0.3\pm 0.3$ | $-2.60\pm 0.02$ | $15.9\pm 6.0$ | $8.1\pm 1.7$ | $25.5\pm 0.1$ | $7.7\pm 0.1$
10056849 | 5.82 | $1022\pm 209$ | $1249\pm 256$ |  $77\pm 20$ (M) | $-18.0\pm 0.1$ | $0.0\pm 0.1$ | $-2.52\pm 0.01$ | $>20.5$ | $<4.4$ | $25.6\pm 0.1$ | $7.4\pm 0.1$
10005113 | 5.82 | $2107\pm 103$ | $1565\pm 109$ | <59 (P) | $-18.0\pm 0.1$ | $0.0\pm 0.2$ | $-2.23\pm 0.02$ | $>12.8$ | $<4.9$ | $25.6\pm 0.1$ | $7.8^{+0.3}_{-0.2}$
22251 | 5.80 | $1318\pm 26$ | $857\pm 18$ | <16 (P) | $-19.0\pm 0.1$ | $0.7\pm 0.2$ | $-2.02\pm 0.01$ | $12.3\pm 1.6$ | $9.8\pm 1.1$ | $25.6\pm 0.1$ | $7.9\pm 0.1$
4404 | 5.78 | $1035\pm 21$ | $793\pm 26$ | <16 (P) | $-19.3\pm 0.1$ | $0.0\pm 0.1$ | $-2.31\pm 0.01$ | $11.7\pm 1.5$ | $6.4\pm 0.3$ | $25.4\pm 0.1$ | –
10016374 | 5.51 | $1758\pm 186$ | $1337\pm 143$ | <19 (P) | $-18.7\pm 0.1$ | $0.4\pm 0.2$ | $-2.21\pm 0.01$ | $10.1\pm 1.6$ | $10.1\pm 1.3$ | $25.5\pm 0.1$ | $7.8\pm 0.1$
9743 | 5.45 | $897\pm 35$ | $996\pm 40$ | <89(P) | $-17.4^{+0.7}_{-4.1}$ | $0.7\pm 0.5$ | $-1.30\pm 0.05$ | $>6.7$ | $<7.3$ | $26.0\pm 0.5$ | –
10015338 | 5.07 | $1639\pm 187$ | $1712\pm 202$ | <18.8 (P) | $-19.5\pm 0.1$ | $0.7\pm 0.4$ | $-2.42\pm 0.01$ | $>14.2$ | $<8.1$ | $25.3\pm 0.1$ | $7.6\pm 0.1$
8113 | 4.90 | $943\pm 19$ | $739\pm 17$ | <49 (P) | $-18.1\pm 0.2$ | $0.5\pm 0.2$ | $-1.69\pm 0.01$ | $6.8\pm 1.0$ | $10.2\pm 1.4$ | $25.6\pm 0.1$ | –
10005217 | 4.89 | $1086\pm 24$ | $1546\pm 37$ | <34 (P) | $-18.1\pm 0.1$ | $0.0\pm 0.1$ | $-2.32\pm 0.01$ | $17.8\pm 2.0$ | $4.6\pm 0.3$ | $25.9\pm 0.1$ | –
7938⋆ | 4.82 | $891\pm 18$ | $620\pm 12$ | <8 (P⋆) | $-19.4\pm 0.1$ | $0.5\pm 0.2$ | $-2.40\pm 0.01$ | $11.6\pm 1.9$ | $11.2\pm 1.3$ | $25.2\pm 0.1$ | $7.8\pm 0.1$
18090 | 4.79 | $976\pm 20$ | $768\pm 15$ | <22 (P) | $-18.9\pm 0.1$ | $0.3\pm 0.2$ | $-1.92\pm 0.01$ | $9.2\pm 1.4$ | $11.6\pm 1.5$ | $25.6\pm 0.1$ | $7.8\pm 0.1$
10001892 | 4.77 | $663\pm 176$ | $1113\pm 308$ | <68 (P) | $-16.0^{+0.5}_{-0.8}$ | – | $-1.69\pm 0.04$ | – | – | – | –
8083† | 4.67 | $1808\pm 36$ | $1362\pm 27$ |  $20\pm 4$ (P) | $-18.7\pm 0.1$ | $1.2\pm 0.1$ | $-1.57\pm 0.01$ | $18.5\pm 2.7$ | $10.2\pm 1.4$ | $25.9\pm 0.1$ | –
10000626 | 4.47 | $763\pm 42$ | $923\pm 45$ |  $185\pm 71$ (P) | $-17.0^{+0.4}_{-0.6}$ | $0.0\pm 0.6$ | $-2.13\pm 0.04$ | $>5.3$ | $<7.2$ | $25.5^{+0.2}_{-0.3}$ | $7.3\pm 0.1$
10001916 | 4.28 | $823\pm 87$ | $1322\pm 202$ | <373 (P) | $-16.7^{+0.5}_{-0.8}$ | – | $-1.80\pm 0.08$ | – | – | – | $7.2\pm 0.2$
7892 | 4.25 | $820\pm 16$ | $677\pm 14$ | <14 (P) | $-18.5\pm 0.1$ | $0.0\pm 0.2$ | $-2.15\pm 0.01$ | $12.4\pm 0.7$ | $7.5\pm 0.3$ | $25.5\pm 0.1$ | $7.8\pm 0.1$
6519 | 4.24 | $1461\pm 31$ | $849\pm 18$ |  $46\pm 7$ (P) | $-19.2\pm 0.1$ | $0.5\pm 0.4$ | $-2.25\pm 0.01$ | $20.0\pm 5.8$ | $12.0\pm 2.8$ | $25.3\pm 0.1$ | $7.7\pm 0.1$
7507 | 4.04 | $765\pm 20$ | $738\pm 15$ |  $44\pm 8$ (P) | $-18.8\pm 0.1$ | $0.0\pm 0.2$ | $-2.44\pm 0.01$ | $>11.1$ | $<6.0$ | $25.3\pm 0.1$ | $7.5\pm 0.1$
4270 | 4.02 | $1030\pm 21$ | $814\pm 16$ | <18 (P) | $-18.9\pm 0.1$ | $1.7\pm 0.2$ | $-1.68\pm 0.01$ | $6.1\pm 1.2$ | $10.3\pm 2.0$ | $25.8\pm 0.1$ | $8.0\pm 0.1$
10035295 | 3.58 | $1255\pm 43$ | $870\pm 25$ | – | $-18.1\pm 0.1$ | $0.7\pm 0.4$ | $-1.92\pm 0.01$ | $32.7\pm 10.1$ | $13.4\pm 3.4$ | $25.7\pm 0.1$ | $7.6^{+0.3}_{-0.1}$
2651 | 3.58 | $840\pm 113$ | $529\pm 36$ | – | $-16.3\pm 1.0$ | – | $-1.42\pm 0.05$ | – | – | – | –
18322 | 3.16 | $1363\pm 41$ | $1432\pm 29$ | – | $-16.9^{+0.4}_{-0.6}$ | $0.6\pm 0.2$ | $-2.19\pm 0.05$ | $9.6\pm 1.9$ | $9.4\pm 1.4$ | $25.9^{+0.2}_{-0.3}$ | –
21150 | 3.09 | $816\pm 16$ | $511\pm 11$ | – | $-19.2^{+0.5}_{-1.0}$ | $0.5\pm 0.3$ | $-2.36\pm 0.01$ | $5.4\pm 1.0$ | $10.3\pm 1.7$ | $25.3^{+0.2}_{-0.5}$ | $8.1\pm 0.1$
Table 1:  Extreme emission line galaxies identified within our sample, in
descending redshift order, selected on a rest-frame [OIII]$\lambda 5007$ EW
$>750$Å. We additionally include one galaxy (ID:10001892) consistent with this
threshold within the measured uncertainties. Lyman-$\alpha$ EWs for our high
redshift sample ($z>5.7$) are taken from Jones et al. (2023) and are marked
(P, M) depending if the line flux was determined from the prism or R1000
medium resolution grating. Lyman-$\alpha$ EWs for galaxies at lower redshifts
are determined for this work using the same method as Jones et al. (2023).
$\ast$ Galaxy ID: 10013682 is an extremely high EW Ly$\alpha$ emitter which
was first discussed in Saxena et al. (2023b). $\ddagger$ ID 9422 may be a
nebular dominated galaxy (see Cameron et al. 2023a). $\bigstar$ We note that
while we can only place an upper limit on the Lyman-$\alpha$ EW of galaxy ID:
7938 from NIRSpec spectroscopy, we confirm Lyman-$\alpha$ detection through
legacy MUSE observations. ${\dagger}$ flagged as a potential AGN (Maiolino et
al., 2023). The galaxy coordinates, photometry and emission line measurements
for these galaxies are publicly available and have been presented in Bunker et
al. (2023b) and Rieke et al. (2023b).

###  3.2 Distribution of equivalent widths

Figure 2: Distribution of equivalents widths. Left panel presents the
[OIII]$\lambda 5007$ EW distribution for the whole sample (black) and split
into a high redshift ($z>5.7$, blue) and low redshift ($z<5.7$, orange) sub-
samples. Right panel shows the same distributions for the H$\alpha$ EW. In
each panel the vertical dashed line represents the EELG threshold of an EW
$>750$Å.

In Figure 2, we present the EW distribution of the [OIII]$\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>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 [OIII]$\lambda 5007$ show similar
profiles (see Figure 2). The high-redshift sample exhibits a flatter EW
distribution which extends to high EW (EW$>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 characterising 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 sub-section 4.1, then examine the galaxy properties in
sub-sections 4.2-4.4.

###  4.1 Stacked spectra of EELGs

Figure 3: Stacks of the galaxy sub-samples split by redshift (at $z=5.7$) and
equivalent width (at [OIII]$\lambda 5007$ EW $=750$Å). The $z>5.7$ EELG sample
is shown in red and the $z<5.7$ EELG sample in orange. This redshift split
naturally separates spectra where the [OIII] doublet is resolved or not. For
the lower redshift sample, the spectral resolving power of the prism around
[OIII] is lower. In blue and green we show the respective non-EELG samples at
high and low redshifts. Top panel: The full spectral range. Middle panel: The
rest-UV. Bottom panel: The rest-optical.

Stack | high-z EELG | low-z EELG | high-z non-EELG | low-z non-EELG
---|---|---|---|---
N galaxies | 17 | 19 | 11 | 37
$[$OIII$]\lambda 5007$ EW | 1524$\pm$201 | 1482$\pm$123 | 572$\pm$59 | 221$\pm$7
H$\alpha$ EW | 1129$\pm$158 | 956$\pm$85 | 400$\pm$45 | 187$\pm$3
Ly$\alpha$ EW | 39$\pm$3 | 21$\pm$5 | $<17$ | $<29$
H$\alpha$/H$\beta$ | 2.98$\pm$0.11 | 3.20$\pm$0.13 | 3.00$\pm$0.32 | 3.40$\pm$0.22
A1600 | 0.15$\pm$0.14 | 0.42$\pm$0.15 | 0.18$\pm$0.41 | 0.65$\pm$0.24
O32∗ | 25.8$\pm$4.7 | 19.3$\pm$4.1 | 6.1$\pm$1.6 | 6.01$\pm$1.11
R23∗ | 7.6$\pm$0.7 | 8.3$\pm$0.9 | 9.8$\pm$2.4 | 9.47$\pm$1.51
Ne3O2∗ | 2.9$\pm$0.6 | 1.0$\pm$0.3 | $<1.4$ | $<0.4$
12+$\log_{10}\mathrm{(O/H)}$ | $7.53^{+0.07}_{-0.10}$ | $7.66\pm 0.14$ | $7.95^{+0.32}_{-0.17}$ | $7.96^{+0.25}_{-0.19}$
$\log_{10}(\xi_{\mathrm{ion}}^{\mathrm{HII}})$ | 25.49$\pm 0.03$ | 25.66$\pm 0.04$ | 25.50${}^{+0.08}_{-0.10}$ | 25.37${}^{+0.05}_{-0.06}$
$\beta_{obs}$ | -2.36$\pm$0.03 | -1.97$\pm$0.04 | -2.18$\pm$0.05 | -1.87$\pm$0.02
$\beta_{int}^{*}$ | -2.45$\pm$0.11 | -2.20$\pm$0.12 | -2.28$\pm$0.28 | -2.22$\pm$0.15

Table 2: Galaxy properties derived for our stacked spectra. Properties marked
with (∗) have been corrected for dust. Upper limits are given at $3\sigma$.

In this Section, we stack the prism spectra for four different sub-samples of
our dataset, split by redshift ($z=5.7$) and whether a galaxy meets the
threshold to be considered an EELG ([OIII]$\lambda 5007$ or H$\alpha$ EW
$>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 $>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 normalise 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 broadband photometry.
At higher redshifts ($z>4$), if the spectroscopic continuum is only weakly
detected in the rest-UV (integrated S/N$<3$ in the range 1400-1600Å) we again
adopt the MUV inferred from the broadband imaging, as described in the
Appendix E, this only applied to 2 of our galaxies (ID 10009693 and 10009320).
We then average these UV normalised spectra; we note that because of the
different redshifts of each galaxy, we do not have uniform coverage of all
wavelengths (as shown in Figure 4) and this is taken into account in the
averaging. The top panel of Figure 3 presents the four sub-sets over the full
wavelength range, and we show individual zooms around the Lyman-$\alpha$ and
[OIII]+H$\beta$ wavelengths in the lower panels. In Appendix D 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 sub-sections, along with a
discussion of the distribution of individual galaxy measurements.

Figure 4: The number of galaxies contributing to each of our four stacks as a
function of wavelength. The colours match those in Figure 3. Each individual
galaxy has a spectrum that covers a unique observed wavelength range, due to
their redshift, the position of the dispersed light on the detector and any
bad pixels/contamination. Therefore the number of galaxies contributing to
each rest-frame wavelength in the stack varies.

###  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 sub-samples.
In sub-section 4.2.1 we examine the detected hydrogen emission lines and we
consider how the dust properties, [OIII]$\lambda 5007$/H$\alpha$ flux ratio,
H$\alpha$ and H$\beta$ EWs, and the detection of Lyman-$\alpha$ emission vary
across our sample. In sub-section 4.2.2 we examine how the ionisation
parameter, metallicity and ionising photon production efficiency change with
[OIII]$\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 $S/N>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}K$ and $n_{e}=10^{2}cm^{-3}$,
respectively). At the highest redshifts ($z>7.1$) we lose H$\alpha$ from our
wavelength coverage and instead we consider the H$\beta/$H$\delta$
decrement666The H$\gamma$ Balmer line is intrinsically stronger than H$\delta$
and hence could provide a more sensitive measurement of the dust reddening.
However, we do not consider the H$\gamma$ emission line because it is blended
in our spectroscopy with the [OIII]$\lambda 4363$ emission line. (with an
intrinsic value of 3.88, Osterbrock & Ferland 2006). This applies to 6
galaxies in our full sample of which 4 are EELGs. For 2 of these galaxies
(both EELGs ID: 10013905 & 10058975) we determine a H$\beta/$H$\delta$
decrement consistent with no dust, but for the other 4 no other Balmer lines
are detected leaving the dust attenuation in these galaxies unconstrained
(hence we will exclude these 4 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
sub-sample 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
[OIII]$\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 sub-set 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 [OIII]$\lambda 5007$/H$\alpha$ line fluxes in Figure 5.
As can be seen, in our sample [OIII]$\lambda 5007$ is nearly always stronger
than H$\alpha$, with an average flux ratio for the full sample of
$f_{\rm{[OIII]}\lambda 5007}/f_{H\alpha}=2.01\pm 0.07$ after correction for
dust reddening. Splitting our full sample into EELGs (coloured blue in Figure
5) and non-EELGs (coloured orange), we see only a marginally higher ratio in
the EELG sub-sample with a $f_{\rm{[OIII]}\lambda 5007}/f_{H\alpha}$ flux
ratio of $2.18\pm 0.15$, compared to $1.93\pm 0.08$ for our non-EELG sample.

Figure 5: Dust corrected flux ratio of the [OIII]$\lambda 5007$/H$\alpha$ line
fluxes against the H$\alpha$ luminosity. We colour code the EELG galaxies with
[OIII]$\lambda 5007$ EW $>750$Å blue and non-EELGs in orange. We overlay the
lower redshift ($0.75<z<1.5$) Domínguez et al. (2013) stacks in grey.

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

We can additionally consider how the EW of these strong Balmer rest-optical
emission lines varies with [OIII]$\lambda 5007$ EW. The hydrogen and [OIII]
nebular emission lines are powered by the same sources, typically the hot
ionising 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 Figure 6 we present the H$\alpha$ and
H$\beta$ EWs against the [OIII]$\lambda 5007$ EWs. The EW of both Balmer lines
show a clear positive trend with the [OIII]$\lambda 5007$ EW, with best fit
relations

$\begin{split}\rm{H}\alpha\,\text{EW/\AA}=(0.71\pm 0.01)\rm{[OIII]}\lambda
5007\text{EW/\AA}+(42\pm 1)\end{split}$

and

$\begin{split}\rm{H}\beta\,\text{EW/\AA}=(0.15\pm 0.01)\rm{[OIII]}\lambda
5007\text{EW/\AA}+(1.1\pm 0.9)\end{split}$

We additionally find that the EELG stacks (defined as having [OIII]$\lambda
5007$EW$>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 Figure 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 [OIII]$\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 10Myrs for an instantaneous
burst, e.g., Leitherer et al. 1999).

Figure 6: Top panels: Comparison of the measured [OIII]$\lambda 5007$ EW with
the H$\alpha$ EW and H$\beta$ EW. Bottom panels: Comparison of the measured
$\log_{10}$[OIII]$\lambda 5007$ EW with the $\log_{10}$ O32 (Those with a
black boarder are corrected for dust) and $\log_{10}$Ne3O2 diagnostics when
all lines were detected at a $>3\sigma$ significance. Both O32 and Ne3O2 are
related to the ionisation parameter $\log_{10}(U)$ and we provide an
additional y-axis in each lower panel to present this, using the Witstok et
al. (2021) relations. In each panel, lines undetected at $3\sigma$ are
replaced with a $3\sigma$ limit. In each panel, we include a shaded region to
mark the EELG EW threshold of [OIII]$\lambda 5007>750$Å, and the measured
quantities for our four stacks are shown in large circles as in the legend in
the top-left panel (coloured to match Figure 3). In the bottom-right panel we
change the stack symbol to triangles when Ne3O2 is given as a 3$\sigma$ upper
limit.

Finally, we turn to the Lyman-$\alpha$ rest-UV emission line, which our
NIRSpec prism spectroscopy has sensitive coverage of down to $z>4$. Although
Lyman-$\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 ionisation of the ISM, and starburst-driven outflows) might
be able to clear escape channels for Lyman-$\alpha$ photons (Ostlin et al.,
2005; Herenz et al., 2017). We now look for any evidence of Lyman-$\alpha$
emitters (LAEs) in our stacks and in individual EELGs, and we report strong
evidence that Lyman-$\alpha$ does indeed escape from EELGs at high redshift.

In the middle panel of Figure 3 we zoom in on the rest-frame UV in the stacked
data to study the emergence or otherwise of Lyman-$\alpha$. We first consider
the two EELG stacks. The high-redshift EELG stack has the strongest
Lyman-$\alpha$ emission (Ly$\alpha$ EW = $39\pm 3$Å), with the low-redshift
EELG stack also showing Lyman-$\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 ($<1\mu m$) means that in our
lower redshift stack we are less sensitive to weak Lyman-$\alpha$ emission
close to the spectral break (see Jones et al., 2023).

Turning to the non-EELG stacks, there is no Lyman-$\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<5.7$). Softening of the spectral break may instead be due to either
damped Lyman-$\alpha$ from HI in the surrounding environment (ISM and CGM) of
the galaxies (e.g., Heintz et al., 2023) or be an artefact, resulting from the
softening of the Lyman-break due to convolution with the Line Spread Function
(LSF) of NIRSpec (see Jones et al. 2023).

Individually, we find that out of 17 EELGs at $z>5.7$, $53\%$ (9/17) have
detected Lyman-$\alpha$ emission, as determined by Jones et al. (2023); Saxena
et al. (2023a); Witstok et al. (2023). This would suggest that selections
based on high [OIII]$\lambda 5007$ EW galaxies are also preferentially
selecting LAEs. This LAE fraction is higher than among the high-redshift non-
EELGs, where only $18\%$ (2/11) have detected Lyman-$\alpha$ emission777Galaxy
IDs 3968 and 18179. These galaxies have measured [OIII]$\lambda 5007$ EWs of
$295\pm 19$Å and $241\pm 6$Å, respectively, and Ly$\alpha$ EWs of $57\pm 17$Å
and $41\pm 9$Å, respectively.. These two non-EELG Lyman-$\alpha$ emitters show
moderate [OIII]$\lambda 5007$ EWs ($\sim 300$Å), in both cases the relative
strength of the nebular emission lines (being driven by recent $<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 4 EELGs in our sample, Lyman-$\alpha$ does not
lie within the wavelength window where we have sensitive spectral coverage
with NIRSpec. Between $4.0<z<5.7$ we have a sample of 15 EELGS with 4 showing
Lyman-$\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 Lyman-$\alpha$ emission which is not visible in our low dispersion
prism spectrum. This leaves us with an LAE fraction of 5/15 ($33\%$) of the
EELG galaxies in the redshift range $4.0<z<5.7$, nominally lower than the
fraction at higher redshifts ($z>5.7$) of $53\%$ (9/17), although we note that
the combination of sensitivity and spectral resolution of the prism means that
NIRSpec’s abilty to detect Lyman-$\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 Lyman-$\alpha$ emitters, whereas for the non-EELGs at these
high redshifts only $18\%$ of our sample have Lyman-$\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\%$) non-EELGs at $4<z<5.7$ show Lyman-$\alpha$
emission, the lowest LAE fraction in any of our sub-samples. Like the high-z
non-EELG LAEs these two888Galaxy IDs 5759 and 7304. These galaxies have
measured [OIII]$\lambda 5007$ EWs of $523\pm 13$Å and $233\pm 21$Å,
respectively, and Ly$\alpha$ EWs of $63\pm 10$Å and $34\pm 9$Å, respectively.
also show rest-optical continuum emission from an evolved stellar population.

####  4.2.2 Ionisation parameter, Excitation state and Metallicity

To investigate how the ionisation state of the ISM changes with [OIII]$\lambda
5007$ EW within our full sample we employ the O32999line flux ratio between
[OIII]$\lambda 4959,5007$ and [OII]$\lambda 3727,3729$ and Ne3O2101010line
flux ratio between [NeIII]$\lambda 3869$ and [OII]$\lambda 3727,3729$
diagnostics (e.g., see Tang et al. 2019), which are sensitive to the
ionisation state of the ISM which is shaped by the ionising flux and hardness
of the ionising spectrum. O32 provides a flux ratio between two different
ionisation states of the same elements, while Ne3O2 provides a flux ratio
between two alpha-elements in different ionisation states. The Ne3O2
diagnostic utilises 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 alpha elements and this
abundance ratio is expected to vary only minimally (e.g. Kobayashi et al.,
2020). Together they provide a complimentary estimate of the ionisation
parameter.

In the bottom-left panel of Figure 6 we plot the O32 diagnostic against the
[OIII]$\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
[OIII]$\lambda 5007$ EW of the full sample in the lower-right hand panel. Both
O32 and Ne3O2 are related to the ionisation parameter, and in both these
Figure panels we additionally provide the ionisation parameter derived
following the Witstok et al. (2021) relations.

We find a clear positive correlation in both O32 and Ne3O2 with [OIII]$\lambda
5007$ EW, obtaining best fit relations:

$\begin{split}\log_{10}(\mathrm{O32})=&(0.69\pm
0.03)\times\log_{10}([OIII]\lambda 5007\,\text{EW/\AA})\\\ &-(0.96\pm
0.10)\end{split}$

and

$\begin{split}\log_{10}(\mathrm{Ne3O2})=&(0.78\pm
0.12)\times\log_{10}([OIII]\lambda 5007\,\text{EW/\AA})\\\ &-(2.4\pm
0.4)\end{split}$

This correlation suggests that galaxies with greater [OIII]$\lambda 5007$ EW
have more highly ionised 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 ionising
spectra, and have higher ionising photon fluxes, making the ISM highly
ionised. 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 sub-samples in Figure 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 Figure 6 compared to the non-EELG stacks.

Figure 7: Top panel: The ionisation and excitation conditions of the full
sample as inferred from the dust corrected O32 and R23 strong line diagnostics
(when either [OII] or H$\beta$ is non-detected we provide a $2\sigma$ limit on
the diagnostics). Each data point is coloured by the [OIII]$\lambda 5007$ EW.
The local SDSS galaxy population is presented in the blue shaded region.
Bottom panel: The Oxygen abundance relation with [OIII]$\lambda 5007$ EW, with
metallicities taken from Curti et al. (2023). We exclude 10 galaxies where
only [OIII] and H$\alpha$ were detected as metallicity constraints could not
be placed using only these lines. We include a shaded region to mark the EELG
EW threshold of [OIII]$\lambda 5007>750$Å, and in both panels the measured
quantities for our four stacks are shown in large circles (coloured to match
Figure 3)

To further examine the state of the ISM we employ the dust-corrected O32 and
R23111111Line flux ratio between [OIII]$\lambda 4959,5007$ \+
[OII]$\lambda\lambda 3727,3729$ and H$\beta$, a proxy for the Oxygen abundance
(metallicity). diagnostic diagram (Figure 7), where galaxies of different
ionisation 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 Figure 7. We find that our full sample ($3<z<9.5$) lies above the
local ($z<0.1$) SDSS galaxy locus121212We use archival data from the Sloan
Digital Sky Survey (SDSS; York et al. 2000) Data-Release 7 (Abazajian et al.,
2009). We set a $\geq 5\sigma$ detection requirement on each emission line for
inclusion in Figure 7 and we reject AGN following the method set out in Runco
et al. (2021), 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; Yang et al., 2017b; Cameron et
al., 2023b). We find a trend that galaxies with larger [OIII]$\lambda 5007$
EWs tend to sit further from the local SDSS galaxy distribution, and have much
higher O32.

We observe no clear turn-over in the R23 diagnostic for the highest EWs,
suggesting that while these galaxies are all metal-poor (sub-solar), 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 turn-over 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. (2023). These measurements rely on multiple strong-
line calibrations utilising 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 [OIII] and H$\alpha$ were
detected in these cases). The Curti et al. (2023) metallicities for our EELGs
lie in the range $7.3<12+\log_{10}(O/H)<8$ (where solar metallicity is $8.69$
in these units), suggesting that while these systems are all sub-solar (all
$<0.2Z_{\odot}$) none show extremely low metallicity (below $<0.04Z_{\odot}$).

We find a best-fit relation between the metallicity and the [OIII]$\lambda
5007$ EW of:

$\begin{split}12+\log_{10}(\rm{O/H})=&(-1.4\pm 0.6)\times
10^{-4}(\rm{[OIII]}\lambda 5007\,\text{EW/\AA})\\\ &+(7.86\pm
0.06)\end{split}$

We find a trend that the higher the [OIII]$\lambda 5007$ EW of a galaxy the
lower the metallicity, albeit with some degree of scatter. This trend is
additionally seen in the stacked sub-samples in the lower panel of Figure 7,
with the metallicity measured to be lower in the stacks with the highest
[OIII]$\lambda 5007$ EW. This is expected, since [OIII]$\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 sub-solar oxygen
abundance (up to a point when there are too few oxygen atoms). We note that
the highest EW EELGs in our sample have metallicities of $\sim 0.1Z_{\odot}$,
with highest EWs found in galaxies with measured metallicities above and below
this value not reaching as high. This turn-over is consistent with
photoionisation models, such as in Tang et al. (2021b), that predict only
marginal increase in the [OIII]+H$\beta$ EW between 0.3-0.1$Z_{\odot}$, with a
turn-over around $0.1Z_{\odot}$, after which the EW begins to decrease with
further decrease in metallicity.

###  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. To examine the non-
ionising 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 MUV

Figure 8: UV absolute magnitude determined for our full sample. Left panel:
MUV plotted against spectroscopic redshift. We overlay as a black line the
apparent m$\rm{}_{UV}=29$ magnitude used in the selection of the high-redshift
targets, plotted as dashed ($z>5.7$) and dotted ($z<5.7$) lines to indicate
the UV selection within our full sample. We use galaxies below this line in
our redshift comparison to obtain a consistent selection across our sample.
Right panel: MUV plotted against the [OIII]$\lambda 5007$ EW. We additionally
overlay dashed lines at MUV = -18.0 and -19.5, which correspond to the UV
absolute magnitude thresholds used in Endsley et al. (2023a). In both panels,
we identify our EELG sample in blue and non-EELG sample in orange. Galaxies
are presented as diamonds when MUV was measured directly from the spectroscopy
and as squares when MUV was measured from photometry.  Figure 9: The UV
spectral slope ($\beta$, where $f_{\lambda}\propto\lambda^{\beta}$) plotted
against the [OIII]$\lambda 5007$ EW. Left panel: The observed $\beta$ slope.
Centre panel: $\beta$ slope corrected for dust attenuation. Right panel:
$\beta$ slope plotted against [OIII]$\lambda 5007$ EW, both corrected for dust
attenuation (taking into account the differential extinction of the nebular
line relative to the continuum as in Calzetti et al. 2000). In each panel we
mark galaxies where the $\beta$ slope was measured directly from the prism
spectroscopy with diamonds, and plot as squares those $z<4$ where we rely on
the spectrum at a slight longer rest-frame wavelength and photometry for our
$\beta$ measurement. In all panels, we include a shaded region to mark the
EELG EW threshold of [OIII]$\lambda 5007>750$Å, and the measured quantities
for our four stacks are shown in large circles (coloured to match Figure 3)

We determine the absolute magnitude in the rest-frame UV (MUV 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 MUV. 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 CIV$\lambda 1550$, HeII$\lambda 1640$, [OIII]$\lambda
1660$ and CIII]$\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
shallower 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. We note that for 24 galaxies in our full sample at $z<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 MUV at 1500Å. For these objects, we additionally measure MUV from
the JADES aperture photometry131313We adopt the ‘KRON’ Diameter Rieke et al.
(2023b) circular apertures 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 2
galaxies at $z>4$ (ID: 10009693 and 10009320) with a S/N$<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 MUV.

The MUV and $\beta$ slope measurements for our full sample are shown against
the [OIII]$\lambda 5007$ EW in Figures 8 and 9. Our full sample covers a broad
MUV range, between $-16>$M$\rm{}_{UV}>-20$. As noted in Section 2.2, at
$z>5.7$ the HST-Deep spectroscopy is essentially a rest-UV-based selection (a
proxy for star formation rate) with a magnitude cut of $AB\sim 29$ mag in the
filter covering wavelengths just redward the Lyman-$\alpha$ break, but
uncontaminated by the break. We plot in Figure 8 the corresponding MUV
selection at $z>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 MUV cut, although in Figure 8 we
show where such a cut would have fallen adopting the same rest-UV apparent
magnitude limit as in the $z>5.7$ selection (dotted line). About 12% of the
$z<5.7$ sample (7/57) are UV-faint ($m_{\rm{UV}}>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 Figure 9. The full sample
covers a broad range of $\beta$ slopes ($-0.5>\beta>-3$). The average observed
$\beta$ slope of EELGs of $-2.2\pm 0.5$ is steeper than non-EELGs at $-1.9\pm
0.7$ albeit consistent within the uncertainties.

The UV spectral slope can depend on the level of dust attenuation and the star
formation history. Nebular continuum emission (e.g. Topping et al. 2022;
Cameron et al. 2023a) 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 B 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 [OIII]$\lambda 5007$ EW in the central panel of
Figure 9. As noted in Section 4.2.1, we are unable to constrain the dust
attenuation for 4 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<-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<-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 ionising
photons, however they show low [OIII]$\lambda 5007$ EWs ($<500$Å) and also low
EWs of the Balmer lines $H\alpha$ and H$\beta$. In the right hand panel of
Figure 9 we additionally correct our [OIII]$\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 [OIII]$\lambda 5007$ EWs$>750$Å, there
still remains a sub-sample 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. 2023b), where the lack of
high EWs has been attributed to either extremely-low oxygen abundance or a
high escape fraction of ionising radiation. Extremely low metallicities and
high $f_{\mathrm{esc}}$ can both reduce the luminosity of the [OIII]$\lambda
5007$ line. However, we do not see strong H$\alpha$ emission in these galaxies
with low [OIII]$\lambda 5007$ EWs and blue $\beta$ slopes, which dis-favours
the extremely low metallicity scenario. Hence these particular galaxies may be
“Lyman Continuum leakers" with non-negligible escape fractions of ionising
photons.

To further examine the effect of dust and star formation history on the
observed $\beta$ slope, we compare the with the Balmer decrement for both the
EELG and non-EELG samples in Figure 10. We determine best fit relations

$H\alpha/H\beta=(-1.3\pm 0.3)*\beta+(5.8\pm 0.8)$

for the EELG sample and

$H\alpha/H\beta=(-0.5\pm 0.3)*\beta+(4.1\pm 0.6)$

for the non-EELG sample. For the EELG sample (shown in blue in Figure 10),
there seems to be a tight correlation 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 Figure 10). However, the non-EELG population displays a
higher scatter and typically redder observed $\beta$ slopes that do no exhibit
such a tight correlation with the Balmer decrement. 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 Figure 10, where the high-redshift EELG
stack (shown in blue) 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 Figure 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.

Figure 10: The observed Balmer decrement plotted against the observed UV-slope
($\beta$). Our full sample is split into EELGs (blue) and non-EELGs (orange),
and are modelled by best fit relations shown in the matching colours. The best
fit relation to the full sample is shown in pink. Both sub-samples favour an
intrinsic UV-slope of $\beta\sim-2.3$. When the theoretical curve is plotted
(black) we find good agreement with the EELG sub-sample and poorer agreement
with the non-EELG sub-sample.

####  4.3.2 Ionisation efficiency

We determine the ionisation efficiency ($\xi_{\mathrm{ion}}^{\mathrm{HII}}$,
the ionising photon production rate per unit UV luminosity) of our full sample
using the dust corrected H$\alpha$ flux and MUV for each galaxy. We follow the
method of Tang et al. (2019) to determine $\xi_{\mathrm{ion}}^{\mathrm{HII}}$.
We split our sample by redshift ($z=5.7$) and plot ionisation efficiency as a
function of [OIII]$\lambda 5007$ EW, along with the best fit relations, on
Figure 11. We over plot the $z\sim 2$ best fit from Tang et al. (2019) and we
find a consistent relation with larger $\xi_{\mathrm{ion}}^{\mathrm{HII}}$
being found at higher [OIII]$\lambda 5007$ EWs. This correlation should not be
surprising as the calculation of $\xi_{\mathrm{ion}}^{\mathrm{HII}}$ and EW
are closely related quantities: $\xi_{\mathrm{ion}}^{\mathrm{HII}}$ comes from
the ratio of H$\alpha$ flux to the rest-UV continuum and EW is related to
ratio of the [OIII]$\lambda 5007$ flux to the rest-optical continuum. We note
that the best fit relation for our lower redshift sample lies at higher
$\xi_{\mathrm{ion}}^{\mathrm{HII}}$ than our high-redshift sample, but these
are both consistent within the uncertainties. We determine a best fit relation

$\begin{split}\log_{10}(\xi_{\mathrm{ion}}^{\mathrm{HII}}/{\rm
erg^{-1}Hz})=&(0.59\pm 0.03)\log_{10}([OIII]\lambda 5007EW/{\text{\AA}})\\\
&+(23.6\pm 0.1)\end{split}$

for the high-redshift sample and

$\begin{split}\log_{10}(\xi_{\mathrm{ion}}^{\mathrm{HII}}/{\rm
erg^{-1}Hz})=&(0.73\pm 0.03)\log_{10}([OIII]\lambda 5007EW/{\text{\AA}})\\\
&+(23.4\pm 0.1)\end{split}$

for the low-redshift sample.

We find that the average $\xi_{\mathrm{ion}}^{\mathrm{HII}}$ in the EELG
sample (high and low redshift combined) is
$\log_{10}(\xi_{\mathrm{ion}}^{\mathrm{HII}}/{\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 ionising spectra
(bluer UV spectral slopes) of the EELGs (in line with results from Simmonds et
al. 2024 using the JADES photometric dataset). We find this pattern is
repeated in Figure 11 in the stacks of EELGs (both high and low redshift)
which are either consistent or have higher ionisation efficiencies than the
non-EELG stacks.

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

Figure 11: $\log_{10}(\xi_{\mathrm{ion}}^{\mathrm{HII}}/{\rm erg^{-1}Hz})$
ionisation efficiency (corrected for dust attenuation) plotted against the
[OIII]$\lambda 5007$ EW, split between the high-redshift (blue, $z>5.7$) and
low-redshift (orange, $z<5.7$) samples. We overlaid Tang et al. (2019) $z\sim
2$ best fit relation. The shaded region marks the EELG EW threshold of
[OIII]$\lambda 5007>750$Å and the measured quantities for our four stacks are
shown in large circles (coloured to match Figure 3).

###  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 star forming history. 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 Spectrum modelling

We model the full prism spectrum ($0.6\lesssim\lambda_{obs}/\mu m\lesssim
5.2$) for each galaxy using BEAGLE (version 0.27.10, Chevallard & Charlot,
2016) to constrain the galaxies’ physical properties and star formation
histories. Briefly, for each galaxy the modelled redshift is set to the
spectroscopic value, and the formation redshift is set to $z=30$
(corresponding to an age of the Universe of 100 Myr for the adopted
cosmology). We utilise the photoionisation models from Gutkin et al. (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 prep). 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 component141414$SFR\propto t\times
exp(-t/\tau)$, a function of time ($t$) and the e-folding time ($\tau$), where
we fit for the $e$-folding time, $\tau$, the maximum stellar age, and the
normalisation of the star formation rate and a constant 10Myr burst (starting
at a lookback time of 10Myr), where the two components are allowed to coincide
simultaneously in time. The choice of a 10Myr burst 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
10Myr, and the disruption timescale of stellar birth clouds (Murray et al.,
2010; Murray, 2011). For the single stellar models used here, more than 99% of
the ionising flux is emitted by stars younger than 10Myr.

For 64 galaxies in our full sample, deep NIRCam photometry from JADES provides
complementary constraints on the spectral energy distribution and are included
as part of the SED fitting. For the remaining 21 galaxies, which lie outside
the footprint of the JADES data release NIRCam imaging, we constrain the SED
using the prism spectroscopy alone. The full BEAGLE data products used in this
paper will be presented in Chevallard in prep. 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.

Parameter | Range | Description
---|---|---
$\log_{10}(\tau/\rm{yr})$ | [6, 12] | Star formation timescale
$\log_{10}(Z_{\rm{ISM}}/Z_{\odot})$ | [-2.2, 0.4] | Metallicity, $Z_{\odot}=0.0152$
$\log_{10}(M_{*}/M_{\odot})$ | [6, 12] | Mass (integral of SFH)
$\log_{10}(t_{\rm{max}}/\rm{yr})$ | [6, 10.8] | Maximum stellar age
$\log_{10}(\rm{SFR/M_{\odot}\rm{yr}^{-1}})$ | [-4, 4] | Star formation rate
z | [z${}_{\rm{spec}},\sigma=0.01$]⋆ | Redshift
$\rm{z}_{form}$ | $\rm{z}_{form}=30$ | formation redshift
$\log_{10}(U)$ | [-4, -1] | Galaxy-wide ionization parameter
$\log_{10}(\xi_{\mathrm{d}})$ | [0.1, 0.5] | Galaxy-wide dust-to-metal mass ratio
$\tau_{\rm{v,eff}}$ | [0, 6] | V-band attenuation optical depth
$\mu$ | $\mu=0.4$ |  | Relative V-band attenuation from
---
diffuse ISM to total

Table 3: BEAGLE free parameter variable ranges set for spectral model fitting.
Uniform priors are set over each range, except for the redshift(⋆) which
follows a Gaussian prior with mean set to $z_{\rm{spec}}$ and $\sigma=0.01$.

####  4.4.2 Stellar Mass and SFH burst intensity

Figure 12: Total mass (Left panel) and Mass-weighted-age (Right panel) against
the [OIII]$\lambda 5007$ EW of our full sample. Galaxies where BEAGLE utilised
both the spectroscopy and the photometry are shown as diamonds and those using
just the spectroscopy as squares. Two potential AGN are shown in green. The
shaded region in each panel marks the EELG EW threshold of [OIII]$\lambda
5007>750$Å. Here, EELGs show typically lower masses and mass-weighted-ages.
Figure 13: Burst fractions against $\log_{10}($[OIII]$\lambda 5007EW)$ for our
full sample. Left panel: the ratio of the average SFR measured over the last
10Myr to the average SFR measured over the last 100Myr. Right panel: fraction
of total stellar mass (stars and stellar remnants) formed in the last 10Myr.
In both panels the burst intensity typically needs to be high for an object to
be an EELG (i.e., strong recent star formation). The two plots plateau at a
burst fraction of 1.0 and 10.0, respectively, when the maximum light-weighted-
age of the galaxy is less than 10Myrs old. We bin the data into 4 EW bins of
equal number and over plot the best fit to the data. Marker shapes follow
those in Figure 12. The shaded region marks the EELG EW threshold of
[OIII]$\lambda 5007>750$Å.

We determine galaxy properties for 84/85 galaxies in our full sample, with
only one galaxy151515NIRSpec ID: 3184 unable to be fit by BEAGLE due to low
signal-to-noise in the continuum. We present the full range of star formation
history property constraints in Tables 4 and 5 in Appendix C. In Figure 12 we
present the BEAGLE best fit total stellar mass (including stars and stellar
remnants) and mass-weighted-age of the stellar population as a function of the
[OIII]$\lambda 5007$ EW.

We highlight the region identifying EELGs (shaded region in Figure 12) and
find that the EELG sample on average exhibits lower masses (typically
$\lesssim 10^{8}M_{\odot}$) and lower mass-weighted-ages ($\lesssim 10$Myr)
than the non-EELG sample. While both EELG and non-EELG sub-samples cover a
broad range of masses, the dynamic range of mass-weighted-ages in the EELG
sub-sample is considerably smaller than the non-EELG sample. Here, the high
sSFR needed in EELGs for the nebular lines to dominate over the rest-optical
stellar continuum (to produce a large EW), drives the mass-weighted-age to low
values. On the other hand, the larger dynamic range among non-EELGs reflects
the greater variety of SFHs in this sub-sample.

To further examine the evolution of these properties with time we can use the
star forming histories of each galaxy to assess their duty cycle (what
fraction of time a galaxy spends in an active star forming phase). The
intensity of star formation over the last 10Myr (SFR10) is responsible for
nebular emission. To quantify if the star formation rate has recently spiked,
we compare SFR10 to the average star formation rate over a longer timescale of
100Myr (SFR100) as shown in the left panel of Figure 13. The majority of EELGs
(bar two outliers161616These are seen in the right hand panel of Figure 13; ID
21150 (EW=816$\pm 16$Å, sSFR=$11\pm 1$Gyr-1) and ID 19342 (EW=2090$\pm 72$Å,
sSFR=$52\pm 13$Gyr-1), which have uncharacteristically low sSFR for their EW.
Inspection of the BEAGLE output shows that 21150 has a SED that fits well the
photometry but has a poor fit to the spectrum, potentially linked to modelling
of the line spread function as this is the lowest redshift object in our
sample. The SED of 19342 is fit well by both the photometry and spectroscopy,
suggesting here that the chosen SFH may be too simplistic. ) have recent star
formation rates (within the past 10Myr) which are $>8$ times higher than the
longer-term average (over the past 100Myr), implying that a recent sharp
increase in the star formation rate is required to produce an EELG (indeed
most EELGs are close to the upper bound ratio of 10, where all the recent star
formation occurred within the past 10Myr). We also show in the right panel of
Figure 13 the fraction of the total mass formed within the past 10Myr, which
is proportional to the sSFR (closely related to the mass doubling time).

We find that both the SFR10/SFR100 ratio and the sSFR show a positive
correlation with $\log_{10}$ ([OIII]$\lambda 5007$ EW), with the galaxies with
the highest EWs showing burst intensities close or on the upper limit
(indicating very young stellar ages $\leq 10$Myr with most or all of the
stellar mass formed within this time). This positive trend is seen clearly
when we bin our galaxies into four [OIII]$\lambda 5007$ EW quartiles of equal
number, and overlay these in Figure 13. The full sample is best fit by the
relations

$\begin{split}\mathrm{SFR_{10}/SFR_{100}}\,=\,&(5.6\pm
0.6)\,\log_{10}([\mathrm{OIII}]\lambda 5007\,\mathrm{EW}/{\text{\AA}})\\\
&-(7.5\pm 1.7)\end{split}$

and

$\begin{split}\mathrm{sSFR(Gyr^{-1})}\,=\,&(62\pm
7)\,\log_{10}([\mathrm{OIII}]\lambda 5007\,\mathrm{EW}/{\text{\AA}})\\\
&-(96\pm 19)\end{split}$

where the fraction of the total mass produced in the last 10Myr (using SFR10)
is equal to sSFR (Gyr-1)/100.

We note considerable scatter around these relations (shown in Figure 13),
which may mostly comes from constraining the SFH, but may also indicate that
it might not be star formation history alone which governs the EW of the rest
optical lines, but other factors such as dust attenuation and metallicity may
play a role.

This correlation between the [OIII]$\lambda 5007$ EW and the burst strength
(or sSFR) suggests that galaxies going through a period of active star
formation may only enter an EELG phase if the stellar mass formed in the burst
makes up a significant proportion of the total stellar mass of the galaxy; the
BEAGLE SED-fitting of our sample indicates that EELGs have $\gtrsim 80\%$ of
their total light-weighted-mass formed in last 10Myr. This threshold for an
EELG phase is reflected in the low mass-weighted ages ($<10$Myr) found from
the BEAGLE fitting of our EELGs compared to the broad range of mass-weighted
ages for the non-EELG sample (see Figure 12).

While these burst mass fraction (sSFR >80 Gyr-1) and burst SFR
(SFR10/SFR${}_{100}>8$) thresholds seem to be required for a galaxy to be an
EELG, there are some non-EELGs which also meet these criteria (i.e., this
‘burstiness’ appears to be a necessary but not sufficient condition for a
galaxy to be an EELG). In each panel of Figure 13 there are galaxies that have
high star formation burst intensity values yet have relatively low
[OIII]$\lambda 5007$ EWs of $<500$Å. These low EW outliers at high burst
intensity may be generated by a number of mechanisms: as already discussed in
Section 4.3.1, high escape fractions of ionising photons will reduce the
luminosity (and EW) of the nebular lines; alternatively, low EWs may arise if
a galaxy is observed very soon after an intense burst of star formation has
ended (within the last 3Myrs), when it will still have a high time-averaged
SFR10 (and hence a large sSFR and SFR${}_{10}/$SFR100 ratio) but the most
massive O-stars which produce many of the ionsing photons will have died,
leading to a rapid decline in the EW.

The high SFR10/SFR${}_{100}>8$ seen in our EELG sample suggests that we are
catching these galaxies in an “on phase" of active star formation which is
short ($<10\%$) 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).

Taken at face value, it seems curious that EELGs comprise 42% of our
spectroscopic sample (rising to 62% for the $z>5.7$ sub-sample), while from
the SFH we infer a very short “on phase” where a young galaxy ($<100$Myr) will
spend only $\sim 10$% 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. If a
high-redshift EELG were to enter an off-duty phase with a factor of $\sim 8$
lower SFR, the rest-frame UV magnitude would be more than two mags fainter
(assuming the UV luminosity is proportional to the SFR, e.g. Kennicutt 1998).
As can be seen from Figure 8, with such UV dimming it may not enter our UV-
selected sample. 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, 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 Lyman-$\alpha$ emission line
detections in our EELG sample in Section 5.2. Finally, we will discuss in
Section 5.3 how the fraction of star-forming galaxies 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:

  1. 1.

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

  2. 2.

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 (sub-
solar, although not extremely metal poor). This reflects that the
collisionally excited [OIII]$\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 ionisation
parameters), typical of low-metallicity galaxies at high redshift.

  3. 3.

These young stellar populations must be part of a rising star forming history
\- either a very recent burst of star formation activity, or a more extended
star formation history which has recently seen a rapid upturn in star
formation rate. Typically EELGs in our full sample exhibit much higher star
formation rates averaged over the most recent 10Myr (SFR10) than the average
over a longer 100Myr interval (SFR100), with EELGs typically having
SFR10/SFR${}_{100}>8$.

  4. 4.

The mass created in the recent burst (in the last 10Myr and responsible for
the nebular emission) must make up a significant fraction of the total mass
formed in the last 100Myr. This is seen in the trend that high EW galaxies
have the highest sSFRs ($>80$ Gyr-1). Typically EELGs have SFHs where $>80\%$
of their total mass has been formed in the last 10Myr. This ensures 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 ionising 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 ionising 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$\%$, 1/36, ID:8083) exhibits a broad
H$\alpha$ emission line and a narrow [OIII]$\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 ionising radiation field is not solely
generated by star formation but may have a contribution from other hard
ionising sources (e.g., an AGN).

###  5.2 Connection to LAEs

As we reported in Section 4.2.1, the fraction of galaxies with detected
Lyman-$\alpha$ emission is higher in our EELG samples than in our non-EELG
samples, with $53\%$ (9/17) of high-redshift EELGs determined to also be
Lyman-$\alpha$ emitters. This measurement is inline with Chen et al. (2023)
who determined that 9/24 (38$\%$) galaxies with [OIII]+H$\beta$ EW$>1000$Å
(roughly [OIII]$\lambda 5007$ EW $>750$Å following the Tang et al. 2019
relations) have Lyman-$\alpha$ detections, which they use to determine that
$50\pm 11\%$ of $z\sim 6$ galaxies exhibit high Lyman-$\alpha$ escape fraction
($f_{\mathrm{esc}}^{Ly\alpha}>0.2$). The production of Lyman-$\alpha$ photons
may be expected to be higher in EELGs considering that within these active
star-forming galaxies, the young stellar populations will produce ionising
photons which should yield Lyman-$\alpha$ recombination line emission.
However, typically Lyman-$\alpha$ photons are strongly suppressed within the
galaxies through resonant scattering in an optically thick environment. In the
case of LAEs where Lyman-$\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 Lyman-$\alpha$ photons
to escape in our direction. The high LAE fraction (53$\%$) 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 Lyman-$\alpha$ photons that are produced in a galaxy and
escape from the interstellar and circum-galactic medium (ISM and CGM), our
observation of any Lyman-$\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 ionised bubbles, either
created by their own stellar populations or from close proximity neighbours.
The large ionised bubbles allow the transmission of Lyman-$\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 9 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 Lyman-$\alpha$ emission
lie within known overdense regions at $z=5.8-5.9$ (Helton et al., 2023;
Stanway et al., 2004, 2007) and $z=7.3$ (Helton et al., 2023; Endsley et al.,
2023a), as can been seen in Figure 1. As noted by Witstok et al. (2023), these
overdensities may lead to a highly-ionised IGM and enhanced transmission of
Lyman-$\alpha$. However, not all of our EELG-LAEs coincide with these known
overdense regions and it is unlikely that these systems created their own
ionised 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
Lyman-$\alpha$ emission is lower at 33$\%$ than the high-redshift EELG sub-
sample. This is consistent with the findings of (Stark et al., 2010; Stark et
al., 2011) who determine the fraction of LBGs that are Lyman-$\alpha$ emitters
($X_{Ly\alpha}$) increases with redshfit 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 Lyman-$\alpha$ emission, after IGM
absorption has been corrected for, may even be more dramatic.

###  5.3 Evolution of the EELG fraction

Figure 14: The fraction of galaxies in an EELG phase ([OIII]$\lambda 5007$
EW$>750$Å) as a function of redshift. In the left panel we plot the EELG
fraction within our full sample (those with spectroscopically confirmed
redshfit), and in the right hand panel we plot the EELG fraction additionally
considering JADES targets which had no emission line detections and were
excluded from our full sample (under the assumption that the photometric
redshift is correct and the lack of any detected line emission means these are
not EELGs). In both panels, our full sample is split into three redshift bins
of roughly equal size (orange diamonds, $\blacklozenge$). The two lower
redshift bins ($z<5.7$) are shown before and after an apparent UV magnitude
cut (mag$=29$, see Figure 8) is applied to match our high redshift bin (shown
without and with a black border, respectively). We overlay the
spectroscopically determined $z\sim 2$ EELG fraction from Boyett et al.
(2022a) (red diamonds, $\blacklozenge$) and combine these with the data points
in our work to determine the best fit power law relation and $1\sigma$
uncertainty (blue line and shaded region). We additionally overlay in grey
EELG fractions determined from photometric surveys. The Endsley et al. (2021)
$z\sim 7$ IRAC/Spitzer EELG fraction is shown as a square ($\blacksquare$),
with the Boyett et al. (2022a) best fit relation using this point shown as the
black dashed line. Crosses, stars and triangles (✖$,\bigstar,\blacktriangle$)
present the Endsley et al. (2023a) JADES/NIRCam EELG fractions for their
bright, faint and very faint MUV samples. The plus (✚) shows the Endsley et
al. (2023b) CEERS/NIRCam EELG fraction.

As we have discussed, galaxies may spend periods of their star forming history
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<z<6$). In the left hand panel of Figure 14
we show the fraction of galaxies in our full sample that exhibit an
[OIII]$\lambda 5007$ EW $>750$Å split into 3 roughly equal size redshift bins
(our high-redshift sample at $z>5.7$ becomes one bin, and we divide our low-
redshift sample equally into two bins at $3<z<4.1$ and $4.1<z<5.7$). We note
that at $z>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$. We consider the uncertainty on each EW measurement when
computing an uncertainty on our EELG fractions.

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) 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\%$ at $z\sim 2$ to $\sim 40\%$ at
$z\sim 7$. We overlay the Boyett et al. (2022a) relation for the redshift
evolution of the EELG fraction in Figure 14. We additionally overlay two
recent [OIII]+H$\beta$ EW distribution studies by Endsley et al. (2023b, a) at
$6<z<9$, and here we determine the EELG fraction adopting the [OIII]-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. (2023a) have been split into a bright
(M$\rm{}_{UV}<-19.5$), faint ($-19.5<$M$\rm{}_{UV}<-18$, ) and a very faint
(M$\rm{}_{UV}>-18$) sub-sample.

We find that the EELG fraction in our spectroscopically confirmed sample
increases with redshift, from $0.23^{+0.04}_{-0.01}$ at $z=3.6\pm 0.1$ in our
lowest redshift bin ($3<z<4.1$) to $0.61\pm 0.04$ at $z=6.5\pm 0.2$ in our
highest ($z>5.7$), after we have applied the UV-selection (see left panel in
Figure 14). We note that in each of our two lower redshift bins, we measure
consistent EELG fractions (within $1\sigma$) with and without applying our
additional UV selection criteria. 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 normalisation of the EELG fraction at $z=0$ to be $0.30\%\pm 0.01\%$.
This relation means that half of all galaxies in a UV-magnitude selection such
as ours are in an EELG phase at $z>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. In
the right hand panel of Figure 14, we plot the EELG fraction for the same bins
where we now consider the complete list of JADES targets. If we make the
assumption that the 133 galaxies with expected redshifts $z>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<z<4.1$) the EELG fraction changes from
$0.23^{+0.04}_{-0.01}$ using the spectroscopic sample to $0.12\pm 0.02$ using
the full JADES target list. In the middle redshift bin ($4.1<z<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<z<9.5$) the fraction changes from $0.61\pm 0.04$
based on the full spectroscopic 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. (2023a) ‘bright’sample (M$\rm{}_{UV}<-19.5$). We again determine a best
fit power law model, with a measured $P=2.8\pm 0.1$ and the expected
normalisation of the EELG fraction at $z=0$ to be $0.18\%\pm 0.01\%$.

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., Langeroodi et al., 2023; Curti et al., 2023), further boosting
the [OIII] emission line luminosity. The evolution of these galaxy properties
means the average observed [OIII] EW in the galaxy population will increase
with redshfit 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 timescale 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. (2023a) measure a decrease in the EELG fraction
for each of their MUV 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 [OIII] emission line luminosity or other
changes in the galaxy population properties (e.g., $f_{\mathrm{esc}}$). 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.

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 ionising
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 reionisation of the Universe.

##  6 Conclusions

Through a JWST/NIRSpec spectroscopic sample of star forming galaxies over a
broad redshift range ($3<z<9.5$), we examine what galaxy properties cause a
galaxy to enter an EELG phase. Out of our sample of 85 SFGs, 42% (36/85)
exhibit [OIII]$\lambda 5007$ rest-frame equivalent widths consistent with
being $>750$ Å. This sub-sample of EELGs shows high EWs in both the
[OIII]$\lambda 5007$ and H$\alpha$ emission lines. They are characterised by a
hard ionising radiation field, observed through the ionisation state of the
ISM (with O32 $\gtrsim 10$); a sub-solar, although not extremely metal poor,
oxygen abundance ($7<\log_{10}(\mathrm{O/H})<8.5$); a high ionising 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 [OIII]$\lambda 5007$
line luminosity, and hence the EW. The high ionisation 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.5$). Through further modelling of the full prism spectrum, we
constrain the star forming histories of our full sample, determining that
EELGs typically exhibit large sSFRs ($>80$Gyr-1), low mass-weighted-ages
($\lesssim 10$Myr) and a high ratio of the star formation rate averaged over
the last 10Myr compared to the longer timescale 100Myr average
(SFR${}_{10}/$SFR${}_{100}>8$).

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 low EW
([OIII]$\lambda 5007$ EW < 500Å) systems in our full sample also exhibiting a
sub-set of these characteristics (e.g., blue $\beta$ slopes or high
SFR${}_{10}/$SFR100). 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 ionising photons in modulating the
equivalent width (as it is the UV-ionising photons that do not leak from the
galaxy that drive the nebular line emission).

The high star forming burst intensity (SFR${}_{10}/$SFR100) in our BEAGLE SED-
fitting of our EELGs suggests that SFGs only enter an EELG phase during short-
lived periods of rapid star formation. These on-phases of the duty cycle are
mirrored by longer term off-phase periods and a dimming of their UV
luminosity. SFGs therefore 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 star forming histories in our EELGs may also be responsible for the
high detection rate of Lyman-$\alpha$ emission, with $53\%$ (9/17) of high
redshift ($z>5.7$) EELGs confirmed to be Lyman-$\alpha$ emitters, compared to
only 18% (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 Lyman-$\alpha$ photons, through the ionisation or removal of
neutral gas via outflows along our sight line.

Across our full sample ($3<z<9.5$), 42% of our sample satisfy our
[OIII]$\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$% of our full sample at $z>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
ionisation efficiency in EELGs makes these galaxies productive ionising photon
producers. EELGs during the EoR may therefore contribute significantly to the
ionising radiation field responsible for reionising the Universe.

## Acknowledgements

This work is based in part on observations made with the NASA/ESA/CSA James
Webb Space Telescope. The data were obtained from the Mikulski Archive for
Space Telescopes at the Space Telescope Science Institute, which is operated
by the Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-03127 for 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 $\&$ 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). SC acknowledges support by European
Union’s HEERC Starting Grant No. 101040227-WINGS. MC acknowledge support by
the Science and Technology Facilities Council (STFC),ERC Advanced Grant 695671
"QUENCH". S.A. 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). S.C acknowledges support by European
Union’s HE ERC Starting Grant No. 101040227 - WINGS. S. Alberts 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 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. JW acknowledges
support from the Science and Technology Facilities Council (STFC), by the ERC
through Advanced Grant 695671 “QUENCH”, by the UKRI Frontier Research grant
RISEandFALL. RS acknowledges support from a STFC Ernest Rutherford Fellowship
(ST/S004831/1). T.J.L. acknowledges support by the Science and Technology
Facilities Council (STFC) and ERC Advanced Grant 695671 "QUENCH". M.V.M.
acknowledges support from the National Science Foundation via grant NSF AAG
2205519. and the Wisconsin Alumni Research Foundation via grant MSN251397.
L.S. acknowledges support by the Science and Technology Facilities Council
(STFC) and ERC Advanced Grant 695671 "QUENCH".

## 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 stacked sample galaxy
properties (Tables 1 and 2) are made available in a machine readable format
upon publication.

## References

  * Abazajian et al. (2009) Abazajian K. N., et al., 2009, ApJS, 182, 543
  * Amorín et al. (2010) Amorín R. O., Pérez-Montero E., Vílchez J. M., 2010, The Astrophysical Journal Letters, 715, L128
  * Amorín et al. (2015) Amorín R., et al., 2015, Astronomy and Astrophysics, 578, A105
  * Atek et al. (2011) Atek H., et al., 2011, The Astrophysical Journal, 743, 121
  * Bacon et al. (2023) Bacon R., et al., 2023, A&A, 670, A4
  * Beckwith et al. (2006) Beckwith S. V. W., et al., 2006, The Astronomical Journal, 132, 1729
  * Boyett et al. (2022a) Boyett K. N. K., Stark D. P., Bunker A. J., Tang M., Maseda M. V., 2022a, MNRAS, 513, 4451
  * Boyett et al. (2022b) Boyett K., et al., 2022b, ApJ, 940, L52
  * Brunker et al. (2020) Brunker S. W., Salzer J. J., Janowiecki S., Finn R. A., Helou G., 2020, The Astrophysical Journal, 898, 68
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Bunker et al. (2020) Bunker A. J., NIRSPEC Instrument Science Team JAESs Collaboration 2020, in da Cunha E., Hodge J., Afonso J., Pentericci L., Sobral D., eds, Proceedings of the International Astronomical Union Vol. 352, Uncovering Early Galaxy Evolution in the ALMA and JWST Era. pp 342–346, doi:10.1017/S1743921319009463
  * Bunker et al. (2023a) Bunker A. J., et al., 2023a, arXiv e-prints, p. arXiv:2306.02467
  * Bunker et al. (2023b) Bunker A. J., et al., 2023b, A&A, 677, A88
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, The Astrophysical Journal, 533, 682
  * Cameron et al. (2023a) Cameron A. J., Katz H., Witten C., Saxena A., Laporte N., Bunker A. J., 2023a, arXiv e-prints, p. arXiv:2311.02051
  * Cameron et al. (2023b) Cameron A. J., et al., 2023b, A&A, 677, A115
  * Cardamone et al. (2009) Cardamone C., et al., 2009, MNRAS, 399, 1191
  * Ceverino et al. (2018) Ceverino D., Klessen R. S., Glover S. C. O., 2018, MNRAS, 480, 4842
  * Charlot & Fall (2000) Charlot S., Fall S. M., 2000, The Astrophysical Journal, 539, 718
  * Chen et al. (2023) Chen Z., Stark D. P., Mason C., Topping M. W., Whitler L., Tang M., Endsley R., Charlot S., 2023, arXiv e-prints, p. arXiv:2311.13683
  * Chevallard & Charlot (2016) Chevallard J., Charlot S., 2016, MNRAS, 462, 1415
  * Cullen et al. (2016) Cullen F., Cirasuolo M., Kewley L. J., McLure R. J., Dunlop J. S., Bowler R. A. A., 2016, MNRAS, 460, 3002
  * Curti et al. (2020) Curti M., Mannucci F., Cresci G., Maiolino R., 2020, MNRAS, 491, 944
  * Curti et al. (2023) Curti M., et al., 2023, arXiv e-prints, p. arXiv:2304.08516
  * De Barros et al. (2019) De Barros S., Oesch P. A., Labbé I., Stefanon M., González V., Smit R., Bouwens R. J., Illingworth G. D., 2019, MNRAS, 489, 2355
  * Domínguez et al. (2013) Domínguez A., et al., 2013, The Astrophysical Journal, 763, 145
  * Du et al. (2020) Du X., Shapley A. E., Tang M., Stark D. P., Martin C. L., Mobasher B., Topping M. W., Chevallard J., 2020, The Astrophysical Journal, 890, 65
  * Eisenstein et al. (2023) Eisenstein D. J., et al., 2023, arXiv e-prints, p. arXiv:2306.02465
  * Eldridge & Stanway (2022) Eldridge J. J., Stanway E. R., 2022, ARA&A, 60, 455
  * Endsley et al. (2021) Endsley R., Stark D. P., Chevallard J., Charlot S., 2021, MNRAS, 500, 5229
  * Endsley et al. (2023a) Endsley R., et al., 2023a, arXiv e-prints, p. arXiv:2306.05295
  * Endsley et al. (2023b) Endsley R., Stark D. P., Whitler L., Topping M. W., Chen Z., Plat A., Chisholm J., Charlot S., 2023b, MNRAS, 524, 2312
  * Faucher-Giguère (2018) Faucher-Giguère C.-A., 2018, MNRAS, 473, 3717
  * Ferland et al. (2017) Ferland G. J., et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Ferruit et al. (2022) Ferruit P., et al., 2022, A&A, 661, A81
  * Furlanetto & Mirocha (2022) Furlanetto S. R., Mirocha J., 2022, MNRAS, 511, 3895
  * Giavalisco et al. (2004) Giavalisco M., et al., 2004, The Astrophysical Journal Letters, 600, L93
  * Gutkin et al. (2016) Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757
  * Heintz et al. (2023) Heintz K. E., et al., 2023, arXiv e-prints, p. arXiv:2306.00647
  * Helton et al. (2023) Helton J. M., et al., 2023, arXiv e-prints, p. arXiv:2302.10217
  * Herenz et al. (2017) Herenz E. C., Hayes M., Papaderos P., Cannon J. M., Bik A., Melinder J., Östlin G., 2017, A&A, 606, L11
  * Izotov et al. (2011) Izotov Y. I., Guseva N. G., Thuan T. X., 2011, The Astrophysical Journal, 728, 161
  * Jakobsen et al. (2022) Jakobsen P., et al., 2022, A&A, 661, A80
  * Jones et al. (2023) Jones G. C., et al., 2023, arXiv e-prints, p. arXiv:2306.02471
  * Kashino et al. (2023) Kashino D., Lilly S. J., Matthee J., Eilers A.-C., Mackenzie R., Bordoloi R., Simcoe R. A., 2023, ApJ, 950, 66
  * Katz et al. (2023) Katz H., et al., 2023, MNRAS, 518, 270
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, Annual Review of Astronomy and Astrophysics, 36, 189
  * Kewley & Dopita (2002) Kewley L. J., Dopita M. A., 2002, The Astrophysical Journal Supplement Series, 142, 35
  * Kimm et al. (2015) Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900
  * Kobayashi et al. (2020) Kobayashi C., Karakas A. I., Lugaro M., 2020, ApJ, 900, 179
  * Kumari et al. (2023) Kumari N., Smit R., Leitherer C., Witstok J., Irwin M. J., Sirianni M., Aloisi A., 2023, arXiv e-prints, p. arXiv:2307.00059
  * Langeroodi et al. (2023) Langeroodi D., et al., 2023, ApJ, 957, 39
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Ma et al. (2018) Ma X., et al., 2018, MNRAS, 478, 1694
  * Maiolino et al. (2023) Maiolino R., et al., 2023, arXiv e-prints, p. arXiv:2308.01230
  * Mármol-Queraltó et al. (2016) Mármol-Queraltó E., McLure R. J., Cullen F., Dunlop J. S., Fontana A., McLeod D. J., 2016, MNRAS, 460, 3587
  * Maseda et al. (2018) Maseda M. V., et al., 2018, The Astrophysical Journal, 854, 29
  * Matthee et al. (2023) Matthee J., Mackenzie R., Simcoe R. A., Kashino D., Lilly S. J., Bordoloi R., Eilers A.-C., 2023, ApJ, 950, 67
  * Miralda-Escudé (1998) Miralda-Escudé J., 1998, ApJ, 501, 15
  * Murray (2011) Murray N., 2011, ApJ, 729, 133
  * Murray et al. (2010) Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191
  * Oesch et al. (2023) Oesch P. A., et al., 2023, MNRAS, 525, 2864
  * Oke & Gunn (1983) Oke J. B., Gunn J. E., 1983, The Astrophysical Journal, 266, 713
  * Onodera et al. (2020) Onodera M., et al., 2020, ApJ, 904, 180
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Ostlin et al. (2005) Ostlin G., Hayes M., Mas-Hesse J. M., Kunth D., Leitherer C., Petrosian A., 2005, arXiv e-prints, pp astro–ph/0507373
  * Rieke (2020) Rieke M., 2020, in da Cunha E., Hodge J., Afonso J., Pentericci L., Sobral D., eds, Proceedings of the International Astronomical Union Vol. 352, Uncovering Early Galaxy Evolution in the ALMA and JWST Era. pp 337–341, doi:10.1017/S1743921319008950
  * Rieke et al. (2023a) Rieke M. J., et al., 2023a, PASP, 135, 028001
  * Rieke et al. (2023b) Rieke M. J., et al., 2023b, ApJS, 269, 16
  * Runco et al. (2021) Runco J. N., et al., 2021, MNRAS, 502, 2600
  * Saxena et al. (2023a) Saxena A., et al., 2023a, arXiv e-prints, p. arXiv:2306.04536
  * Saxena et al. (2023b) Saxena A., et al., 2023b, A&A, 678, A68
  * Simmonds et al. (2023) Simmonds C., et al., 2023, MNRAS, 523, 5468
  * Simmonds et al. (2024) Simmonds C., et al., 2024, MNRAS, 527, 6139
  * Smit et al. (2015) Smit R., et al., 2015, The Astrophysical Journal, 801, 122
  * Stanway et al. (2004) Stanway E. R., et al., 2004, ApJ, 604, L13
  * Stanway et al. (2007) Stanway E. R., et al., 2007, MNRAS, 376, 727
  * Stark et al. (2010) Stark D. P., Ellis R. S., Chiu K., Ouchi M., Bunker A., 2010, MNRAS, 408, 1628
  * Stark et al. (2011) Stark D. P., Ellis R. S., Ouchi M., 2011, ApJ, 728, L2
  * Sun et al. (2023) Sun F., et al., 2023, ApJ, 953, 53
  * Tacchella et al. (2018) Tacchella S., Bose S., Conroy C., Eisenstein D. J., Johnson B. D., 2018, ApJ, 868, 92
  * Tacchella et al. (2020) Tacchella S., Forbes J. C., Caplar N., 2020, MNRAS, 497, 698
  * Tacchella et al. (2023) Tacchella S., et al., 2023, MNRAS, 522, 6236
  * Tang et al. (2019) Tang M., Stark D. P., Chevallard J., Charlot S., 2019, MNRAS, 489, 2572
  * Tang et al. (2021a) Tang M., Stark D. P., Chevallard J., Charlot S., Endsley R., Congiu E., 2021a, MNRAS, 501, 3238
  * Tang et al. (2021b) Tang M., Stark D. P., Chevallard J., Charlot S., Endsley R., Congiu E., 2021b, MNRAS, 503, 4105
  * Topping et al. (2022) Topping M. W., Stark D. P., Endsley R., Plat A., Whitler L., Chen Z., Charlot S., 2022, ApJ, 941, 153
  * Wilkins et al. (2011) Wilkins S. M., Bunker A. J., Stanway E., Lorenzoni S., Caruana J., 2011, MNRAS, 417, 717
  * Williams et al. (2023) Williams C. C., et al., 2023, ApJS, 268, 64
  * Willott et al. (2022) Willott C. J., et al., 2022, PASP, 134, 025002
  * Withers et al. (2023) Withers S., et al., 2023, ApJ, 958, L14