© The Authors 2023

1. Introduction

Emission lines can be some of the most prominent features of galaxy spectra. Regions of photoionised gas around young stars (HII regions) routinely exhibit bright emission from hydrogen and helium recombination, as well as collisionally excited emission from numerous metal ions, especially O⁺, O²⁺, Ne²⁺, S⁺, and N⁺ (e.g. Peimbert et al. 2017). Ratios of emission line fluxes are sensitive to abundances of metals and their ionisation states, but also the physical conditions of the emitting nebulae, such as temperature and density, which are in turn linked to the nature of the ionising sources powering the emission (e.g. Osterbrock & Ferland 2006). Emission-line ratios therefore have great diagnostic power for understanding properties of the interstellar medium (ISM) and stellar populations in galaxies (see Kewley et al. 2019 for a review).

However, observed emission lines in galaxies can arise from other phenomena, including active galacitc nuclei (AGN; Kewley et al. 2006; Feltre et al. 2016) or shock-heated gas (Dopita & Sutherland 1996; Allen et al. 2008; Sutherland & Dopita 2017). Meanwhile, stellar populations other than the young O and B stars typically associated with HII regions, such as high-mass X-ray binaries (HMXBs; Senchyna et al. 2020) or post-AGB stars (Byler et al. 2019), can power emission in photoionised regions of gas. Each of these heating sources yields characteristically different emission spectra. An observed global emission spectrum of a star-forming galaxy is typically a superposition of emission from many HII regions with different conditions as well as possible contributions from these additional sources, making the interpretation of emission spectra challenging.

Despite this complexity, galaxies have been observed to follow remarkably tight trends in emission line-ratio space. The so-called BPT and VO87 diagrams, which relate the ratio of [O III] λ5007/Hβ to ratios of either [N II] λ6583/Hα, [S II] λλ6716, 6731/Hα, or [O I] λ6300/Hα (Baldwin et al. 1981; Veilleux & Osterbrock 1987), are perhaps the best known examples of such trends¹. With large surveys like Sloan Digital Sky Survey (SDSS; York et al. 2000), star-forming galaxies at z ∼ 0 have been observed to form a tight sequence in the [N II]–BPT diagram, with a general anti-correlation between the two ratios (e.g. Kauffmann et al. 2003; Kewley et al. 2006). This tight correlation is the result of the fact that the two parameters that most strongly dictate the emission spectra of galaxies, metallicity, and ionisation parameter² (U), are themselves generally correlated in galaxies (Dopita & Evans 1986; Carton et al. 2017). The position within this sequence is thus dominated by these two parameters, although variations in other properties including the N/O abundance ratio, electron density, and radiation hardness also contribute (e.g. Kewley et al. 2013a; Bian et al. 2016; Curti et al. 2022). These correlations are not unique to the BPT and VO87 diagrams, and indeed other combinations of emission lines, such as ([O II] + [O III])/Hβ versus [O III]/[O II] (‘R23–O32’), have been observed to exhibit similar trends (e.g. Veilleux & Osterbrock 1987; Kewley et al. 2006).

Existing studies of earlier cosmic epochs show that, compared to the loci inhabited by local galaxies in line-ratio space, z  ∼  2 galaxies are generally offset to higher excitation and higher ionisation (e.g. Strom et al. 2017; Runco et al. 2021). This is perhaps best exhibited by the widely studied offset in the [N II]–BPT diagram (Kewley et al. 2013b; Steidel et al. 2014, 2016; Runco et al. 2022), where z ∼ 2 galaxies are observed with typically higher [N II]/Hα at a fixed [O III]/Hβ. This is often explained as a combination of higher ionisation parameters, higher densities, and harder ionising spectra at high redshift (Brinchmann et al. 2008; Kewley et al. 2013a; Steidel et al. 2016; Sanders et al. 2016; Hirschmann et al. 2017). Many authors have suggested that super solar α/Fe abundance ratios at these earlier times play an important role by driving harder stellar ionisation fields (e.g. Strom et al. 2017; Topping et al. 2020b,a). However, the N/O abundance ratio is clearly another important parameter to consider, with some studies suggesting this offset can be explained as the result of higher N/O at fixed [O III]/Hβ (Masters et al. 2016).

Despite the wealth of rest-frame optical studies highlighting the different ISM conditions in galaxies at z ∼ 2, emission line studies of galaxies at z ≳ 6 were, until recently, very limited. Broad-band excesses had already suggested an abundance of galaxies at this epoch with very high equivalent width emission from [O III] + Hβ (e.g. Endsley et al. 2021, 2023; Fujimoto et al. 2023). However, constraints on their ISM conditions were limited to studies of far-infrared lines, such as [O III]88 μm and [C II]158 μm from the Atacama Large Millimeter Array (ALMA), which suggested high ionisation parameters and enhanced α/Fe abundance ratios (e.g. Hashimoto et al. 2019; Carniani et al. 2020; Harikane et al. 2020; Fujimoto et al. 2022; Witstok et al. 2022).

The recent commissioning of JWST is transforming our ability to probe ISM physics in the early Universe. Already, emission line measurements from medium-depth observations from the initial months of JWST operations have observed z ≳ 6 galaxies to exhibit generally low metallicity, high ionisation, high excitation, and high temperature (Curti et al. 2023; Katz et al. 2023; Sun et al. 2023; Rhoads et al. 2023; Tacchella et al. 2023; Trump et al. 2023; Langeroodi et al. 2022; Williams et al. 2023; Heintz et al. 2023; Matthee et al. 2023; Fujimoto et al. 2023; Mascia et al. 2023; Sanders et al. 2023; Nakajima et al. 2023; Tang et al. 2023). Furthermore, the observed emission spectra of some of these targets have been found to be consistent with the presence of AGN activity (Brinchmann 2023), while Katz et al. (2023) suggested that the presence of high-mass X-ray binaries could help to explain the observed emission-line ratios.

Although much insight has already been gained from these early observations, they have largely been based on relatively shallow spectroscopy. Here we leverage data taken from the ‘Deep’ spectroscopic tier of the JWST Advanced Deep Extragalactic Survey (JADES), the deepest spectroscopic observations yet taken with NIRSpec, to provide a more detailed look at the emission-line ratios of galaxies at z ∼ 5.5 − 9.5. These multi-object spectroscopy (MOS) observations reached exposure times of up to 28 h in the PRISM/CLEAR (R ∼ 100) and up to 7 h in each of the three medium resolution gratings (R ∼ 1000), providing unprecedented new insights into the ISM properties of galaxies within the first Gyr of the Universe’s history. The structure of this paper is as follows. In Sect. 2 we describe the observations, data reduction and measurement of emission line fluxes. In Sect. 3 we present our emission-line ratios on various diagnostic diagrams and provide comparisons with literature samples. Section 4 then presents a discussion of the implications of these observations. We summarise in Sect. 5. Throughout this paper we adopt the Planck Collaboration VI (2020) cosmology: Ω_Λ = 0.6847, Ω_m = 0.3153, H₀ = 67.36 km s⁻¹ Mpc⁻¹. All magnitues are quoted in the AB magnitude system (Oke & Gunn 1983).

2. Data and analysis

2.1. Observations

The data presented in this paper were obtained via multi-object spectroscopic observations from JWST/NIRSpec using the micro-shutter assembly (MSA; Jakobsen et al. 2022; Ferruit et al. 2022). Observations were carried out as part of JADES in three visits between 21 and 25 Oct. 2022 (Programme ID: 1210; PI: N. Luetzgendorf) in the Great Observatories Origins Deep Survey South (GOODS-S) legacy field in a region overlapping with the Hubble Ultra Deep Field (HUDF; Koekemoer et al. 2013). Each visit consisted of 33 613 s integration in the PRISM/CLEAR low-resolution setting and 8403 s integration in each of G140M/F070LP, G235M/F170LP, G395M/F290LP, and G395H/F290LP filter and grating settings. Across three visits, this totals 28 h of integration in the PRISM, which provides continuous spectral coverage from 0.6 to 5.3 μm at R ∼ 30 − 300, and ∼7 h in each of the medium resolution gratings, which combine to provide R ∼ 700 − 1300 across the full spectral range of NIRSpec, plus 7 h in the high-resolution grating which provides R ∼ 2700 from ∼2.8 to 5.1 μm. In this paper we present only measurements from the PRISM/CLEAR and G395M/F290LP observations. Results obtained from other configurations will be presented in forthcoming work.

Observations within each visit were performed as a 3-shutter nod. The central pointing of each visit was dithered (by < 1 arcsec) such that common targets were observed in different shutters and different detector real-estate. Thus, each visit had a unique MSA configuration, although target allocation (performed with the eMPT³) was optimised for maximising target commonality between all three dither positions. A total of 253 unique targets were observed in the PRISM configuration with the three dithers featuring 145, 155, and 149 targets respectively. 67 targets were triply exposed in the PRISM and reached full depth, while 62 targets featured in only two dithers and a further 124 received only one-third total depth. The hierarchical allocation of shutters means that the double- and triple-exposed targets are biased towards the highest priority (and hence highest redshift) targets. Indeed, of the 27 galaxies presented in this work, 18 received triple coverage, six received double, and three received only single coverage. Details of the target selection scheme are presented in Bunker et al. (2023).

The low dispersion associated with the PRISM mode means that all targets are observed with non-overlapping spectra. However, in the higher resolution modes, individual spectra are dispersed over a much larger extent of detector real-estate. To minimise the possibility of contaminating emission, in our grating observations we isolate our highest priority targets by closing the shutters of low-priority targets on the same row (i.e. targets that could cause overlapping spectra). Thus, for our grating spectra we observe only 198 unique targets (119, 121, and 111 in each dither). However, the high-priority z ∼ 5.5 − 9.5 galaxies presented in this work are almost unaffected by this measure with only two exceptions. One target (JADES-GS+53.11351–27.77284) that was triply exposed in PRISM was only singly exposed in the gratings, while one target (JADES-GS+53.17582–27.77446) was not observed in the medium resolution setting at all. The total integration times obtained on each target are summarised in Table A.1.

2.2. Data reduction

These observations were processed by adopting algorithms developed by the ESA NIRSpec Science Operations Team (SOT) and the NIRSpec GTO Team. Details of the data-processing workflow will be presented in a forthcoming NIRSpec GTO collaboration paper. As this work uses the observations carried out with only two spectral configurations of the program, we describe the main steps for the PRISM/CLEAR and G395M/F290LP filter and grating settings. Upon retrieving the level-1a data from the Mikulski Archive for Space Telescopes (MAST), we estimated the count rate per pixel by using the unsaturated groups in the ramp and removing jumps due to cosmic rays identified by estimating the slope of the individual ramps. During this first stage, we also performed the master bias and dark subtraction, corrected snowball artefacts, and flagged saturated pixels. We then performed the pixel-by-pixel background subtraction by combining the three nod exposures of each pointing. We note that for some targets we exclude one of the 3-shutter nods in the background subtraction stage due to the serendipitous presence of a source contaminating the open background shutter. We then created two-dimensional (2D) cutouts of each 3-shutter slit and performed the flat-field, spectrograph optics, and disperser corrections. We then ran the absolute calibration and corrected the 2D spectra for path-losses. Path-loss corrections are dependent on the morphology of the target and its relative position within the shutter. In this paper, we apply the path-loss corrections assuming a point-like source at the location of the target centroid. The high-redshift targets presented in this paper are generally smaller, or comparable in size, to the angular resolution of the telescope at the observed wavelength of emission lines used in this study (λ ≳ 3 μm). We estimate that for the sample presented here, even the ratio with the longest wavelength baseline ([O III] λ5007/[O II] λλ3726, 3729) would be affected only at the 1% (0.01 dex) level by adopting an extended morphology based on fitting to broad-band imaging. Since it is not a priori known how closely the spatial distribution of emission lines trace the starlight observed in broad-band imaging, we decide to adopt the point-like assumption, but acknowledge that this effect introduces some systematic uncertainty to our measurements.

We rectified and interpolated the 2D continuum map onto a regular grid for the G395M/F290LP observation, and an irregular grid for the PRISM/CLEAR to avoid an oversampling of the line spread function at short wavelengths. Finally, the 1D spectra were extracted from the 2D map adopting a box-car aperture as large as the shutter size, centred on the relative position of the target in the shutter. For each target, we combined all 1D spectra and removed any bad pixels with sigma-clipping. Cutouts of five example representative 1D spectra are shown in Fig. 1.

Fig. 1. Example spectra from five galaxies included in our sample. Left: low-resolution PRISM/CLEAR spectra shifted to the rest-frame according to the observed redshift. Flux is shown in fλ but has been renormalised relative to the peak flux from the [O III] λ5007 line. Solid horizontal lines show the zero flux of each spectrum. The shaded region around the ‘zero’ line indicates the 1-σ noise spectrum in the same renormalised units. We show only the subset of spectral coverage that includes the key emission lines used in this study. Vertical dotted lines show the expected centroids of the [O II] λ3727, [Ne III] λ3869, [Ne III] λ3967, Hδ, Hγ, [O III] λ4363, Hβ, [O III] λ4959, [O III] λ5007, [O I] λ6300, Hα, [N II] λ6583, [S II] λ6716, and [S II] λ6731 lines respectively (moving left to right). Right: zoom in on the Hα + [N II] complex as observed in the G395M/F290LP grating spectra of these same galaxies. Observed spectra are shown in the same way as for the left panel (except renormalised to Hα, rather than [O III]). Vertical dotted lines show the expected centroids of Hα and [N II] λ6583.

As a verification of our flux calibration, we compared the NIRSpec spectrophotometry to broad-band photometry from NIRCam and found the differences to be on average less than 10% even before considering the effects of slit-losses due to the MSA. Details of this comparison are presented in Appendix C.

2.3. Emission line ratio sample

Although 253 unique targets were observed in these observations, this study focuses only on a rest-frame ultraviolet (UV) selected sample of confirmed z > 5.5 galaxies with robustly detected emission lines. The criteria used to select targets for this programme are described in detail in Bunker et al. (2023). We describe here only the aspects most relevant to this work.

The highest priorities in target assignment were reserved for galaxies at the highest redshifts, with a preference for the UV-bright galaxies. Candidate z > 5.5 galaxies were primarily drawn from a compilation of existing literature catalogues based on Hubble Space Telescope (HST) broad-band imaging (e.g. Bunker et al. 2004; Bouwens et al. 2015, 2021; Finkelstein et al. 2015; Harikane et al. 2016), with a number of the highest redshift galaxies presented here having been identified with HST (e.g. Bouwens et al. 2011, 2016; Ellis et al. 2013; McLure et al. 2013; Oesch et al. 2013). In fact, only five galaxies presented in this work were newly identified from JWST/NIRCam imaging (see Bunker et al. 2023 for details). In total, 47 galaxies were assigned MSA shutters as z > 5.7 candidates out of 299 possible such candidates within the MSA footprint. Of these, 28 passed redshift-dependent rest-UV cuts to place them in our highest priority classes PC = 1 and PC = 4 (PC = 2 and three contained only five targets between them, none of which featured here; see Bunker et al. 2023 for class descriptions). For this analysis, we consider only galaxies with robust spectroscopically confirmed redshifts in which Hβ is detected with S/N > 5. Of the 28 rest-UV-selected targets in PC = 1 and PC = 4, 18 met these criteria. We also include six (one) galaxies selected as z > 5.7 candidates that had failed the rest-UV selection cut, but were included in PC = 6.1 (PC = 6.2) as a result of meeting the magnitude criterion of m_F160W < 29 (m_F160W > 29) and turned out to match our Hβ S/N criterion for this paper. Finally, three galaxies selected in PC = 7.5 (4.5 < z < 5.7, m_F160W < 29) had z_spec > 5.5 and S/N_Hβ > 5, and we included these in our sample as well. One galaxy exhibited a broad-component under Hα in R ∼ 1000 spectroscopy, suggestive of nuclear activity, and was discarded so as to avoid contamination from emission from AGN. Thus, our final sample included 27 star-forming galaxies, drawn from a largely rest-UV-selected sample.

This work focuses predominantly on the low-resolution PRISM/CLEAR spectra (‘PRISM’ hereafter), since these are significantly deeper than our medium-resolution grating spectra. The spectral resolution of the PRISM varies significantly with wavelength from R ∼ 30 around 1.2 μm up to R ∼ 300 at the long wavelength cutoff around 5 μm (e.g. Fig. 5 in Jakobsen et al. 2022). By limiting ourselves to rest-frame optical emission lines at z ≳ 5.5, the measurements presented in this paper are made at observed wavelengths of λ_obs ≳ 2.4 μm where the resolution is R ≳ 100, and improving steeply as a function of wavelength. Thus, in all cases from the PRISM, the [O III]λλ4959, 5007 doublet is completely resolved from Hβ. The doublet itself is always at least partially resolved, and is completely resolved in our highest redshift targets (Fig. 1). [O II] λλ3726, 3729 and [S II] λλ6716, 6731 are never resolved into two components, but the latter can easily be fully resolved from the Hα + [N II] complex over this redshift range. However, Hα and [N II] λ6583 require resolutions of R ∼ 500 to be fully resolved and thus are significantly blended across the full redshift range. Constraints on this ratio are instead derived from our higher resolution (but shallower) G395M/F290LP grating spectra.

Figure 2 shows the redshift histogram of galaxies in this sample. We divide galaxies into two subsamples, based on the observability of Hα. The purple ‘z ∼ 6’ sample spans over 5.507 ≤ z ≤ 6.931 (z_median = 5.943) and includes 21 galaxies for which Hα is still within the wavelength coverage of NIRSpec. There are an additional six galaxies in our ‘z ∼ 8’ sample, spanning over 7.05 ≤ z ≤ 9.45 (z_median = 7.629), for which we study only emission lines from [O II] λλ3726, 3729 to [O III] λ5007. Since the G395M/F290LP spectra are dispersed across a larger area of detector real estate, they do not provide continuous spectral coverage over the nominal spectral range. For two out of 21 galaxies in the z ∼ 6 sample, we do not have [N II]/Hα constraints because these lines fall into the detector gap of NIRSpec. We additionally have one galaxy (JADES-GS+53.17582–27.77446) which was only observed in the PRISM mode (see Sect. 2.1). For one galaxy in the z ∼ 6 sample, the [S II] doublet falls beyond the observed spectral range.

Fig. 2. Redshift and magnitude distribution of the 27 galaxies in our JADES sample. Top: redshift histogram of galaxies observed in this study. We divide the sample into two sub-samples ‘z∼6’ and ‘z ∼ 8’, based simply on cutting at z = 7 since this is the redshift at which Hα is redshifted beyond the spectral coverage of NIRSpec. Thus, for JADES ‘z ∼ 8’ we only have coverage of lines from [O II] λλ3726, 3729 to [O III] λ5007, while for JADES ‘z ∼ 6’ we can in principle observe all the rest-frame optical lines out to [S II] λλ6716, 6731. The orange hatched histogram indicates galaxies for which we can constrain the [S II] λλ6716, 6731 flux. Green hatching is the same but for [N II] λ6583. Bottom: redshift vs. M1500 for the full z ∼ 5.5 − 9.5 sample. We derived M1500 directly from the PRISM spectra as described in Sect. 2.4.1.

2.4. Spectral analysis and line flux measurements

2.4.1. Calculation of M₁₅₀₀

The continuum is generally well detected in the rest-frame ultraviolet (UV) across our sample. In all except our faintest target, the signal-to-noise on the continuum reaches S/N > 3 for individual spectral pixels longwards of the Lyman break, and in six spectra we reach in excess of S/N ≳ 10 per pixel in the rest-frame UV. For our faintest target, JADES-GS+53.16746–27.77201, we do observe faint continuum in the rest-UV, but the signal-to-noise per pixel is only S/N ≈ 1.8. We derive rest-frame UV absolute magnitudes (M₁₅₀₀) for our sample directly from the spectra. We shift the observed spectra to the rest-frame using the spectroscopically confirmed redshift and integrate across a synthetic 200 Å wide top-hat filter centred on 1500 Å. This integrated flux is then converted into an absolute magnitude. Derived M₁₅₀₀ values are reported in Table A.1. Aside from our faintest target (JADES-GS+53.16746–27.77201 with M₁₅₀₀ = −17.01 ± 0.22), all other targets have magnitudes in the range −17.5 ≤ M₁₅₀₀ ≤ −20.4.

2.4.2. PRISM/CLEAR emission line measurements

Emission line measurements and continuum modelling are performed simultaneously using the penalised pixel fitting algorithm, PPXF (Cappellari 2017, 2023). PPXF models the continuum as a linear superposition of simple stellar-population (SSP) spectra, using non-negative weights and matching the spectral resolution of the observed spectrum. As input, we used the high-resolution (R = 10 000) SSP library combining isochrones from MESA Isochrones and Stellar Tracks (MIST; Choi et al. 2016) and the C3K theoretical atmospheres (Conroy et al. 2019). The flux bluewards of the Lyman break was manually set to 0. These templates are complemented by a 5th-degree multiplicative Legendre polynomial, to take into account systematic differences between the SSPs and the data (e.g. dust, mismatch between the SSP models and high-redshift stellar populations, and residual flux calibration problems). The emission lines are modelled as pixel-integrated Gaussians, again matching the observed spectral resolution. To reduce the number of degrees of freedom, we divide the emission lines into four kinematic groups, constrained to have the same redshift and intrinsic broadening. These groups are: UV lines (bluewards of 3000 Å in the rest-frame), the Balmer series of Hydrogen, non-Hydrogen optical lines (bluewards of 9000 Å in the rest-frame), and near-infrared (NIR) lines. The stellar continuum has the same kinematics as the Balmer lines. Furthermore, we tie together doublets that have fixed ratios, and constrain variable-ratio doublets to their physical ranges. In particular, this study focuses on the following lines of interest: [O II] λλ3726, 3729, [Ne III] λ3869, Hγ, Hβ, [O III] λ4959, [O III] λ5007, [O I] λ6300, Hα, [N II] λ6583, [S II] λλ6716, 6731. We note that, at the resolution of the PRISM/CLEAR mode, the two components of the [O II] λλ3726, 3729 doublet are significantly blended and so are fit as a single Gaussian component. We visually inspect these fits to the 1D spectra to ensure they are reliable. For the faintest lines, we also individually inspect the 2D spectrum for each target and confirm that the emission line can be visually identified along the same trace to confirm the measured flux is not affected by uncorrected artefacts. Where a line is not detected, we derive an upper limit directly from the noise spectrum. The 1-σ upper limit is obtained by integrating the variance spectrum over an interval spanning 1.4 × σ_line either side of the expected centroid of the line (where σ_line is the Gaussian line width), and multiplying this by the spectral extent of a pixel.

We note that, at the spectral resolution of our PRISM observations, [Ne III] λ3869 is partially blended with the He Iλ3889 line. From grating spectra, we expect the [Ne III] line to be brighter than He I, typically by at least a factor of two, suggesting that performing a two-component fit can recover the [Ne III] flux (e.g. Appendix C.3 in Cameron et al. 2021). The fitting procedure outlined above simultaneously fits for both lines in this complex, and we confirm from visual inspection that these fits model the complex well in cases with a high signal-to-noise ratio. In this paper, the only ratio for which we consider the [Ne III] line is the [Ne III] λ3869/[O II] λλ3726, 3729 ratio. This ratio has a very short wavelength baseline and hence the uncertainties introduced by partial blending are somewhat offset by having almost no dependence on wavelength dependent corrections or calibrations. Nonetheless, we note that that the [Ne III]-based measurements presented in this work should be treated with caution, and that any inference requiring high precision [Ne III]-based measurements should be performed with medium- to high-resolution spectra.

2.4.3. Constraints on [N II]/Hα from G395M/F290LP

As discussed in Sect. 2.3, the only line ratio for which we consider the data from our higher spectral resolution, but much shallower, G395M/F290LP spectra is [N II] λ6583/Hα. These lines are sufficiently close in wavelength that they remain at least partially blended in PRISM/CLEAR spectra across the entire redshift range considered here. However, in the R ∼ 1000 G395M/F290LP spectra, this complex is always completely resolved.

We performed a visual inspection simultaneously for the PRISM and G395M/F290LP 1D and 2D spectra of each target to identify the possible presence of [N II]. In fact, we find no convincing evidence for the detection of [N II] λ6583 in any of our 18 galaxies, either from the resolved R ∼ 1000 spectra, or by way of the appearance of a ‘red wing’ in the Hα profile in the PRISM spectra. We derive upper-limits on the [N II]/Hα ratio in the same way as described above, fitting the Hα line with pPXF, and integrating the variance spectrum across the expected spectral extent of [N II] λ6583, based on the redshift and line-width measured for Hα.

2.4.4. Emission line ratio measurements

In this paper we focus on the following key diagnostic line ratios:

(1)

Several of these line ratios, such as N2, S2 and R3, are calculated across short wavelength baselines and are thus largely insensitive to dust attenuation or wavelength-dependent calibration issues. However, to accommodate the analysis of the longer baseline ratios, we correct our emission line fluxes and upper limits for dust extinction according to the Small Magellanic Cloud (SMC) dust law from Gordon et al. (2003). We calculate the E(B − V) from the decrements of Balmer lines. For the z ∼ 6 sample where Hα is available, we use the Hα/Hβ decrement, assuming an intrinsic ratio of Hα/Hβ = 2.86 (Osterbrock & Ferland 2006). For the z ∼ 8 sample, where possible, we use Hγ/Hβ assuming an intrinsic ratio of Hγ/Hβ = 0.468. Above z > 7, Hγ is sufficiently deblended from [O III] λ4363 so as to avoid a biased measurement. In two cases, the S/N on Hγ is too low. For these galaxies, we adopt the attenuation obtained from spectral energy distribution (SED) fitting with BEAGLE (Chevallard et al., in prep.), but find that only very small corrections (E(B − V) < 0.02) are required. Overall, the mean E(B − V) from this sample is 0.1. Our final dust-corrected emission-line ratios are presented in Table 1.

2.5. Stacked spectra

Although this paper focuses on presenting line ratio measurements made for individual objects, we also perform a stack of all the spectra in each of our two redshift bins to aid in assessing median quantities of the sample, and push the limits of emission line non-detections. In our z ∼ 6 bin, we perform this stacking for both the medium-resolution G395M/F290LP spectra (for the N2 ratio) and the low-resolution PRISM spectra (for all other line ratios). In the z ∼ 8 bin we stack only the PRISM spectra.

Each spectrum is de-redshifted to the rest-frame and then re-sampled onto a uniform wavelength grid with spectral pixels 3 Å wide in the rest-frame. For the z ∼ 6 stack, we renormalise each spectrum by the Hα flux, while for the z ∼ 8 stack we renormalise by Hβ. We then take the median across all spectra in each spectral bin to form our final composite spectrum.

In order to extract emission line fluxes, we first subtract off the continuum by masking the regions of bright emission lines and fitting a spline. We then fit each emission line in the continuum-subtracted spectra with single-component Gaussian profiles. We correct the measured emission-line ratios for the effects of dust-reddening in the same way as for the individual measurements, described in Sect. 2.4.4, using the measured Hα/Hβ ratio to derive E(B − V) for the z ∼ 6 sample and the Hγ/Hβ ratio for the z ∼ 8 sample.

We repeat the stacking process and line measurements for each sample where each spectrum is instead normalised by the UV continuum flux (taken as the median flux between 1400 and 1600 Å in the rest-frame spectrum). Finally, we repeat these measurements removing the target with the highest average noise level. Although the ratios measurement from these alternative stacks do not vary so much as to change any of the key findings in this paper, the change is larger than the statistical uncertainties on the measurements. For the final dust-corrected emission-line ratios reported in Table B.2, we adopt the values measured from the Balmer-emission-normalised stacks. However, rather than adopting the statistical uncertainty, we report the systematic uncertainty, which we take as the maximum difference between the values derived from any two of the stacking regimes described above.

3. Diagnostic diagrams

3.1. Comparison samples

In this section, we compare the emission-line ratios measured for our JADES z ∼ 6 and z ∼ 8 samples to various literature samples in several diagnostic diagrams. To form our baseline z ∼ 0 comparison sample, we use 0.03 < z < 0.1 galaxies from SDSS MPA-JHU catalogues⁴ (Aihara et al. 2011) to generate 2D probability density functions (PDFs). We consider only galaxies with spectra flagged as ‘reliable’ in the public release catalogue, and select for only galaxies with prominent emission lines by making a signal-to-noise cut of S/N > 40 on Hα. For diagnostic diagrams that include the [Ne III] λ3869 line, we additionally apply a cut of S/N_λ3869 > 5. In the N2–BPT, S2–VO87, and O1–VO87 diagrams we plot all sources that pass these criteria (including any AGN), however, for all other plots, we additionally select for only star-forming galaxies using the Kauffmann et al. (2003) criterion.

We consider z ∼ 2 galaxies from the MOSFIRE Deep Evolution Field (MOSDEF) survey public release emission line catalogue⁵ (Kriek et al. 2015; Reddy et al. 2015). For this sample we impose signal-to-noise cuts of S/N > 5 on each emission line and plot only detections (i.e. we do not show any limits). We additionally select only galaxies with a sky line flag with < 0.2 for all lines considered, as recommended by the catalogue release. For lines that are flagged as being near the edge of the spectrum, we also select only those flagged as having reliable flux. For diagrams involving N2, S2 or O1, the MOSDEF sample spans from 1.2 < z < 2.6 with median z = 2.1, while for diagrams involving only the bluer lines, it probes slightly higher redshifts: typically 2.0 ≲ z ≲ 3.6 with median of z ≈ 2.3.

A number of line ratio measurements have already been reported from NIRSpec at z ≳ 6. The JWST Early Release Observations (ERO; Pontoppidan et al. 2022) provided emission line measurements for three galaxies at z > 7.5, one of which was reported to have an extremely low metallicity of 12 + log(O/H) = 6.99 ± 0.11 (Curti et al. 2023). We use these three galaxies as comparison points in oxygen-based diagrams, adopting the recently re-measured ratio values from Nakajima et al. (2023). We compare our individual galaxy measurements with individual 4.5 ≤ z ≤ 8.0 galaxies targeted in the GLASS survey, which leverages gravitational lensing (Mascia et al. 2023). Additionally, we compare to composite spectra presented in Sanders et al. (2023) and Tang et al. (2023) from the Cosmic Evolution Early Release Science (CEERS) survey. The z ∼ 5.6 and z ∼ 7.5 redshift bins in Sanders et al. (2023) are well matched in redshift to the sample presented here. Those authors report a median stellar mass of log for their z ∼ 5.6 composite, while they do not estimate masses in their z ∼ 7.5 bin. Tang et al. (2023) present one composite spectrum of 16 z ≳ 7 galaxies from CEERS with median redshift z = 7.7. Those authors report a median UV magnitude of M_UV = −20.6. This is considerably brighter than our median magnitude of M₁₅₀₀ = −18.71, indicating that the shallower CEERS spectroscopy is probing generally brighter galaxies that our JADES data.

We also consider a number of possible ‘high-redshift analogues’ identified from low-redshift samples by various authors. ‘Green pea’ galaxies are a population of galaxies at z ∼ 0 with unusually high equivalent-width emission of [O III] λ5007 (∼1000 Å; Cardamone et al. 2009; Yang et al. 2017a). These galaxies are characterised by masses of M_* ∼ 10^{8.5 − 9.5} M_⊙ and high star formation rates (∼10 M_⊙ yr⁻¹), analogous to galaxies at higher redshifts. We consider also ‘blueberry’ galaxies, which are a population of z ∼ 0 dwarfs with even more extreme properties than green peas, having masses of M_* ∼ 10^{6.5 − 7.5} M_⊙, very high ionisation ([O III]/[O II] ≳10) and very low metallicities (12 + log(O/H)≈7.1 − 7.8), (Yang et al. 2017b). We also consider samples of the most metal-poor galaxies known in the low-redshift Universe (Izotov et al. 2018, 2019).

The z ∼ 0 comparison samples are corrected for dust assuming a Cardelli et al. (1989) law with R_V = 3.1 where the A_V is derived from the Hα/Hβ decrement assuming an intrinsic ratio of 2.86. Meanwhile, the z ∼ 2 MOSDEF sample is corrected in an identical manner to that used for our z ∼ 6 sample, described in Sect. 2.4.

3.2. Excitation properties of z  ∼ 6 galaxies: BPT and VO87 diagrams

The left panel of Fig. 3 shows our z ∼ 6 JADES sample plotted onto the so-called N2–BPT diagram (Baldwin et al. 1981). As described in Sect. 2.4.3, we find no convincing detections of [N II] λ6583 in any of our JADES z ∼ 6 sample (purple points). Figure 3 also shows comparison samples of z ∼ 0 galaxies from SDSS (background 2D PDF), extremely metal-poor galaxies at z ∼ 0 (Izotov et al. 2018, 2019, yellow pentagons), z ∼ 2 galaxies from MOSDEF (green circles), and the z ∼ 5.6 composite CEERS spectrum from Sanders et al. (2023).

Fig. 3. BPT and VO87 diagrams showing the comparison between JADES z ∼ 6 galaxies and samples across various redshifts. Individual measurements are shown as solid purple points, while values derived from stacked JADES spectra are shown as the open magenta diamonds. The background grey 2D PDF shows z ∼ 0 galaxies from SDSS. Extremely metal-poor z ∼ 0 galaxies from Izotov et al. (2018, 2019) are shown as yellow pentagons. MOSDEF z ∼ 2 galaxies are shown as green circles (Kriek et al. 2015). The z ∼ 5.6 composite CEERS spectra presented by Sanders et al. (2023) are shown as the blue open circles. The solid navy lines in each panel show the theoretical maximum starbust lines from Kewley et al. (2001), while the dotted line in the left panel shows the empirical demarcation derived by Kauffmann et al. (2003). Left: classical BPT or ‘N2–BPT’ – Purple triangles provide 3-σ upper limits on the locations of our JADES galaxies since we are unable to detect [N II] λ6583 in any of the 18 galaxies for which we have G395M/F290LP spectral coverage of the Hα complex, although the R3 ratio is very well constrained in all cases from PRISM data. Even after stacking the grating spectra of all galaxies, we do not recover a robust detection of [N II] λ6583. Centre: in the S2–VO87 diagram (often referred to as ‘S2–BPT’), three out of 20 of our JADES galaxies show detections of [S II] in individual spectra. The tight 3-σ upper limits we place on our non-detections highlight that there must be more than an order of magnitude in scatter in the S2 ratio within the sample. Right: we find one tentative detection of [O I] λ6300 in our O1–VO87 diagram.

Based on the R3 ratios alone, we can already see that these JADES z ∼ 6 galaxies span the upper extreme of what is observed for galaxies in SDSS and MOSDEF, with the 18 galaxies shown here having a median R3 of 0.74, and all lying above R3 > 0.53. Considering the N2 ratio, SDSS galaxies lying along the star-forming BPT sequence with 0.72 < R3 < 0.76, have a median value of N2_median = −1.4. None of our JADES galaxies have upper limits on the N2 ratio which are inconsistent with this value at the 3-σ level, while the majority of our sample have individual upper limits of only ≲ − 0.9 (see also Table 1). Beyond our individual non-detections, we find that we are unable to robustly detect [N II] λ6583 even in a stack of all 18 spectra, placing a 3-σ upper limit of N2 < −1.44. Previous studies have widely observed that z ∼ 2 galaxies are generally offset from the z ∼ 0 N2–BPT sequence, towards higher values of N2 at fixed R3 (Kewley et al. 2013b; Steidel et al. 2014; Strom et al. 2017; Runco et al. 2022; green MOSDEF points in Fig. 3). Given the anti-correlation between N2 and R3 in this sequence, and the high R3 values measured in our sample, it is not possible to characterise from individual measurements what fraction of the JADES sample would lie significantly ‘above’ the z ∼ 0 sequence. However, the upper limit imposed from our stacked spectrum suggests the median value at this redshift might be more aligned with the z ∼ 0 sequence than measurements from z ∼ 2 samples.

Sun et al. (2023) made emission line measurements at z ∼ 6 with NIRCam slitless spectroscopy, with two modest detections of [N II] λ6583 showing much higher N2 values of N2  =   − 0.77 and N2 = −0.69 (see also Sun et al. 2022). However, those two galaxies were more than an order of magnitude brighter (m_F444W < 25) than the JADES galaxies probed here (m_{F444W; median} = 28.1, for galaxies with NIRCam coverage), suggesting Sun et al. (2023) are probing more massive systems. Sanders et al. (2023) present six individual galaxies at z > 5 with detections of [N II] λ6583 from CEERS spectroscopy, also showing much higher N2 values (N2 ≳ −1.0). Taking a composite spectrum of all 38 galaxies in the CEERS sample from 5.0 < z < 6.5 they find , quoting a relatively high median stellar mass of log . Our inability to detect N2 in any of our sample may indicate none of our galaxies have reached the same levels of metal-enrichment as some of the brighter individual galaxies observed in Sun et al. (2022) or Sanders et al. (2023). However, comparison with low-redshift galaxies with extremely low metallicity (Izotov et al. 2018, 2019, yellow pentagons) highlights that the 7 h deep G395M/F290LP spectra presented here are well short of being sufficient to probe [N II] emission from the most metal-poor systems in the early Universe. We see in the left panel of Fig. 3 that, at z ∼ 0, these most metal poor systems often show N2 ratios well below N2 < −2, and have generally milder values of R3 (0.0 < R3 < 0.5). At z ∼ 0, the sequence of R3 in galaxies is known to be double-valued with metallicity, increasing from the lowest metallicities up to a turnover around 12 + log(O/H)≈8, before decreasing towards higher metallicity (e.g. Curti et al. 2020). The fact that our JADES R3 ratios are notably higher than those from the Izotov et al. (2018, 2019) galaxies suggests we are not probing galaxies below this turnover and are more likely sampling galaxies only down to around ∼0.1 × Z_⊙, where the R3 sequence is relatively flat (Curti et al. 2020). Understanding N/O abundance ratios in metal-poor galaxies in the Universe will be critical to developing chemical evolution models of galaxies (Maiolino & Mannucci 2019; Hayden-Pawson et al. 2022). This study demonstrates that probing N/O abundances in the most metal poor galaxies at z ≳ 5.5 with emission lines will likely be very challenging even with very deep spectra.

The centre panel of Fig. 3 shows our JADES sample plotted onto the S2–VO87 diagram (Veilleux & Osterbrock 1987)⁶. Because the [S II] λλ6716, 6731 doublet is more widely separated in wavelength from Hα, we can measure this ratio directly from the much deeper PRISM spectrum⁷. Consequently, we observe three [S II] detections at greater than 3-σ. These detections are broadly consistent with lying along the upper extension of the z ∼ 0 sequence and z ∼ 2 sequence. The deeper PRISM spectra also allow us to place much tighter upper limits on the S2 ratio for the 17 out of 20 non-detections in our sample. In some cases the 3-σ upper limits are more than an order of magnitude lower in S2 than our three detections, suggesting there is significant scatter in S2 ratios within the z ∼ 6 population. The composite z ∼ 5.6 CEERS spectrum from Sanders et al. (2023) yields , which is broadly in line with the three JADES galaxies with measurable [S II]. This suggests that stacks based on the shallower, but wider area, CEERS spectroscopy are biased towards galaxies that are more in line with the tip of the z ∼ 0 sequence, which are likely more evolved systems.

In the right panel of Fig. 3 we next plot O1 = [O I] λ6300/Hα against the R3 ratio. As with the S2 ratio, we can again derive [O I]/Hα constraints directly from the PRISM spectra. We find one tentative [O I] detection in our sample (JADES-GS+53.15613−27.77584), indicating the presence of neutral gas in this galaxy. The measured O1 ratio for this galaxy places it above the theoretical maximum starburst excitation limit from Kewley et al. (2001), which could suggest the presence of shock-heated gas in the ISM (Sutherland & Dopita 2017). Interestingly, this galaxy is not one of the three galaxies with a [S II] detection from Fig. 3. However, it does have the second lowest O32 ratio of galaxies in this sample (Table 1). Overall, the individual O1 constraints presented here are not strongly constraining relative to the Kewley et al. (2001) demarcation line. The measurement from our stacked spectrum suggests the median O1 in this sample could be consistent with the upper extension of the z ∼ 0 star-forming sequence. However, we showed above that there is significant diversity in S2 ratios within our sample. If this diversity also holds for O1, it is clear that larger samples probing even deeper would be required to characterise the sample properties of O1 at z ∼ 6.

3.3. Excitation properties from z∼5.5−9.5: R2–R3 diagram

In Fig. 4 we instead consider the R2 = [O II] λλ3726, 3729/Hβ ratio. The longer wavelength baseline of this ratio has the drawback of being much more sensitive to the reddening correction applied compared to the N2–BPT, S2–VO87, and O1–VO87 diagrams discussed in Sect. 3.2. However, the shorter rest-frame wavelength of these emission lines means that the R2–R3 diagram can be used to explore the ISM conditions of galaxies up to z ∼ 9.5 with JWST/NIRSpec. [O II] is typically intrinsically much brighter than [N II], [S II], or [O I]. Moreover, like [S II] and [O I], ratio measurements can be taken directly from the PRISM spectrum. As a result, not only does Fig. 4 include the additional six galaxies from the z ∼ 8 sample, we also have a much higher detection rate, with 24 out of 27 galaxies in our combined sample having detections of [O II].

Fig. 4. R2–R3 diagram showing our two JADES sub-samples at z ∼ 6 (purple points) and z ∼ 8 (maroon points). Comparison samples from SDSS and MOSDEF are shown as in Fig. 3. We additionally show comparison with two populations of z ∼ 0 extreme starbursts: green peas (Yang et al. 2017a; orange plusses) and blueberries (Yang et al. 2017b; blue crosses). The three well-studied z > 7.5 galaxies from the SMACS 0723 ERO observations are shown as lime green pentagons, adopting the ratio values presented by Nakajima et al. (2023). Unlike the BPT and VO87 diagrams, the majority of our sample have > 3σ detections of all lines in individual spectra, revealing a large diversity within the sample. Ratios measured from stacked JADES spectra are shown as the open magenta and red diamonds (see Sect. 2.5).

We see immediately in Fig. 4 that there is more than an order of magnitude in scatter in the R2 ratio, suggesting there is significant diversity in the ISM conditions of galaxies within this sample. We find a median value of R2 = −0.23 with standard deviation of 0.38 among our 24 detections, significantly offset from the median value from MOSDEF of R2 = 0.38. Instead, we find that the range of R2 values observed for this sample are much more aligned with those observed in extreme z ∼ 0 dwarf starbursts such as green peas (Cardamone et al. 2009; Yang et al. 2017a; orange plusses) and blueberries (Yang et al. 2017b; blue crosses).

Both R2 and R3 are known to form a double-valued sequence in metallicity (e.g. Curti et al. 2020). In z ∼ 0 galaxies, R3 does not change strongly with metallicity between ∼0.1 − 0.5 × Z_⊙. On the other hand, R2 turns over at higher metallicities between ∼0.5 − 0.8 × Z_⊙, at a maximum value of R2 ≈ 0.5 (Curti et al. 2020). The combination of highly consistent R3 ratios but wildly varying R2 ratios indcates that our sample likely falls into the metallicity range of ∼0.1 − 0.5 × Z_⊙. We have five galaxies in our JADES sample with R2 > 0.0, overlapping with the region inhabited by typical z ∼ 2 MOSDEF galaxies. These may represent systems which are more evolved and more metal-rich. Below its metallicity turnover, R2 values in the z ∼ 0 sample from Curti et al. (2020) drop steeply with decreasing metallicity. The lowest R2 values measured in our sample lie below R2 < −0.75, below the range probed by the Curti et al. (2020) calibration. However, we note that the low R2 portion of our sample exhibits a fair resemblence to blueberry galaxies from Yang et al. (2017b) which have been show to have very low-metallicties ≲0.1 × Z_⊙. However, given that R2 is strongly affected by the ionisation parameter, inferring metallicity values from this ratio is likely not robust. Nonetheless, there is clearly a large range of metallicity being probed in this sample.

We see no evidence for evolution of the R2 ratio between our z ∼ 8 and z ∼ 6 samples, albeit with only six z ∼ 8 galaxies. We note that all of our R2 > 0 galaxies are from the z ∼ 6 sample; however, the measurements from stacked spectra in the two samples actually yield a higher value for the z ∼ 8 sample in Fig. 4 (open diamonds), likely affected by uncertainty due to small sample size. This suggests that there is no rapid time evolution in this diagram at this epoch and that the scatter is instead driven by sample diversity.

3.4. Ionisation-excitation diagrams: R23–O32 and R23–Ne3O2

Figure 5 shows the R23–O32 diagram with our JADES sample plotted against the same literature comparison samples shown in Fig. 4. R23 is often taken as an indication of total excitation, as it encompasses emission lines from both singly and doubly ionised oxygen, while the O32 ratio is sensitive to ionisation parameter (e.g. Kewley et al. 2019).

Fig. 5. R23–O32 diagram, probing the ionisation and excitation of the ISM. The JADES sub-samples at z ∼ 6 (purple points) and z ∼ 8 (maroon points) lie at generally high excitation an high ionisation. Comparison samples, including SDSS, MOSDEF, green peas, blueberries, and SMACS 0723, are shown as in Fig. 4. z ∼ 4.5 − 8 galaxies from GLASS are shown as grey diamond (Mascia et al. 2023). Blue and orange open circles show the z ∼ 5.6 and z ∼ 7.5 composite CEERS spectra from Sanders et al. (2023), while the z ∼ 7.7 composite CEERS spectrum from Tang et al. (2023) is shown as the green open hexagon. Stacked JADES spectra are shown as the open magenta and red diamonds.

It is in this diagram that the offset of z ≳ 6 from typical z ∼ 0 and z ∼ 2 galaxies is perhaps most striking: O32 ratios are clearly much higher in JADES z ∼ 5.5 − 9.5 galaxies compared to MOSDEF galaxies at z ∼ 2. Across our combined sample we find a median O32 value of 0.98 with a standard deviation of 0.36, with 18 of our galaxies have robust detections showing O32 > 0.75. This is much higher than MOSDEF galaxies plotted in Fig. 5 which have a median O32 value 0.18 and are already tracing the upper end of the z ∼ 0 relation. This indicates that galaxies across our sample exhibit very high ionisation parameters, much larger than what are seen in typical galaxies at z ∼ 0 − 3. However, we can see in Fig. 5 that our JADES sample largely overlaps with the region spanned by blueberries (Yang et al. 2017b) and green peas (Yang et al. 2017a), which are examples of z ∼ 0 galaxies with extreme star formation activity. Comparable O32 ratios have also been observed in z ∼ 2 − 4 Lyman-α emitters (Nakajima et al. 2016; Tang et al. 2021) and z ∼ 1.3 − 2.4 extreme [O III] emitters (Tang et al. 2019). It is therefore likely that the high O32 ratios observed in our sample are indicative of intense star formation activity.

We find that our z ∼ 6 sample is significantly offset from the z ∼ 5.6 stack presented in Sanders et al. (2023), with our composite O32 value found to be more than a factor of four larger. Indeed where individual detections from CEERS are reported in Sanders et al. (2023) they typically have O32 < 0.75. Mirroring what was observed in our R2 sequence from Fig. 4, we find six galaxies with much lower O32, more in line with the composite spectrum from Sanders et al. (2023). Given the resemblance between the upper end of our O32 sample and dwarf starbursts such as blueberries, it is likely that this portion of our sample is probing higher star formation rate and lower-metallicity galaxies (analogous to blueberries), while the handful of lower O32 galaxies presented here in alignment with the Sanders et al. (2023)z ∼ 5.6 composite could be probing more evolved systems.

On the other hand, the composite CEERS values at z ∼ 7.5 and z ∼ 7.7, from Sanders et al. (2023) and Tang et al. (2023) respectively, are in closer agreement with our composite spectra, albeit still offset towards slightly higher R23. This could suggest that the fraction of more evolved, metal-enriched systems is increasing rapidly from z ∼ 8 to z ∼ 6, despite the bulk of our UV-selected z ∼ 6 sample showing no significant evolution from our z ∼ 8 sample. Addressing this question in more detail requires larger samples of deep spectroscopy where properties of individual galaxies can be isolated.

By leveraging gravitational lensing, the measurements of 4.5 ≲ z ≲ 8 galaxies from Mascia et al. (2023) can more readily probe fainter and lower mass systems. These are similarly offset from the Sanders et al. (2023)z ∼ 5.6 composite and largely overlap with our JADES sample, spanning a large range in R23. Two galaxies are reported by Mascia et al. (2023) as having R23 ≲ 0.6 despite having high O32 comparable with our JADES sample. This indicates that gravitational lensing may be important for study low metallicity galaxies (below the R23 turnover) and not simply ultra-deep observations like those presented here.

We note that, as with Fig. 4, this diagram relies on ratios with long baselines in wavelength and that dust reddening can have a significant effect. The dust-correction (described in Sect. 2.4.4) is most significant on the O32 ratio here, with larger corrections serving to reduce O32 ratios. From the corrections we apply, all but two galaxies have corrections that shift O32 by less than 0.15. The offset between the Sanders et al. (2023)z ∼ 5.6 composite and our z ∼ 6 stack cannot be explained by the dust correction; ‘correcting’ our value down to O32 ≈ 0.5 would require E(B − V)≳1, significantly larger than our measured value of E(B − V) = 0.06.

Nonetheless, in Fig. 6 we show our JADES z ∼ 6 and z ∼ 8 samples with Ne3O2 replacing O32 to provide the ionisation axis, which provides the bonus that Ne3O2 is essentially unaffected by dust extinction and the associated correction uncertainties⁸. Despite the lower detection rate, and generally much larger measurement uncertainties associated with switching from [O III] λ5007 to the intrinsically much fainter [Ne III] λ3869, we see the same story reflected here with the JADES sample exhibiting large scatter (median Ne3O2 = −0.11 with standard deviation 0.29) and generally offset from SDSS and MOSDEF galaxies, instead being more aligned with blueberry dwarf starbursts. This again supports the conclusion that our JADES sample is probing a diverse sample of galaxies, many of which exhibit very high ionisation parameters analogous to those observed in local compact starbursts. Tang et al. (2023) presented Ne3O2 measurements from their z ∼ 7.7 composite CEERS spectrum. These are marginally lower than our Ne3O2 composite values. As in Fig. 5, the larger offset is again the 0.1 dex difference in R23. This offset could reflect an increase in metallicity from our sample to the Tang et al. (2023) galaxies which are reported has having typically brighter UV magnitudes (median M_UV = −20.6 from Tang et al. 2023).

Fig. 6. R23–Ne3O2 diagram, probing the ionisation and excitation of the ISM. R23–Ne3O2 closely resembles the R23–O32 diagram with the additional advantage that the [Ne III] λ3869/[O II] λλ3726,329 ratio is less sensitive to the applied dust-correction, or any wavelength-dependent flux calibrations. JADES galaxies are plotted as in Fig. 5. Fewer literature samples present [Ne III] ratios so we show only SDSS, blueberries, MOSDEF, and the z ∼ 7.7 composite CEERS spectrum from Tang et al. (2023). These are each plotted in the same way as in previous figures.

Finally, we note that [Ne III] λ3869 emission may become important for confirming spectroscopic redshifts at z ≳ 10 where [O III] λλ4959, 5007 is no longer observable within the wavelength coverage of NIRSpec. We find the [Ne III] λ3869 emission to be brighter than the [O II] doublet in at least 8 out of 27 galaxies in our sample. Given that these cases should represent the more extreme examples (highest ionisation parameters and lowest metallicities), [Ne III] λ3869 may be one of the most accessible emission lines in searches for the highest redshift galaxies. However, existing deep spectroscopic observations of such galaxies have already demonstrated that observing any emission lines at those redshifts can be challenging (Curtis-Lake et al. 2023; Robertson et al. 2023).

4. Discussion

4.1. ISM conditions in the epoch of reionisation

Section 3 showed that our z ∼ 5.5 − 9.5 sample from JADES exhibits significant diversity in ISM conditions. The measured ratios highlight a clear difference between this high-redshift sample and typical galaxies at z ∼ 0 − 3, although parallels can be drawn instead with extreme starburst populations of low-redshift galaxies, such as blueberries and green peas. In this section we use photoionisation models to explore the parameter space spanned by our sample.

Metallicity and ionisation parameter are the dominant parameters driving the emission-line ratios of H II regions (Kewley et al. 2019). In Fig. 7, we compare our measured ratios to photoionisation model grids in four of the diagnostic plots highlighted in Sect. 3. These models are calculated with CLOUDY v17 (Ferland et al. 2017) adopting input Binary Population and Spectral Synthesis (BPASS) v2.0 ionising spectra (Eldridge & Stanway 2016; Stanway et al. 2016), assuming an instantaneous burst with an age of 3 Myr. The models are plane-parallel and assume a 1 kpc thick shell with n_H = 100 cm⁻³. The inner radius is set to 10³⁰ cm and ionisation parameter is varied between −4.0 < log U < −0.5 in steps of 0.5. Models are calculated at 11 values of metallicity: Z = [1, 2, 3, 4, 6, 8, 10, 14, 20, 30, 40] × 10⁻³ (or Z/Z_⊙ = 0.07, 0.15, 0.22, 0.30, 0.45, 0.60, 0.75, 1.0, 1.5, 2.2, 3.0, assuming Z_⊙ = 0.0134; Asplund et al. 2009). Metal abundances ratios are assumed to follow solar abundance patterns (Asplund et al. 2009), and the depletion of metals onto dust is accounted for.

Fig. 7. JADES z ∼ 5.5 − 9.5 galaxies compared to photoionisation models. Measured ratios are shown as purple points, as in Figs. 3–6 except that all galaxies are plotted in the same colour, regardless of redshift. SDSS galaxies are shown in the same manner as Figs. 3–6. Each line shows a set of models with constant metallicity with values of Z = [1, 2, 3, 4, 6, 8, 10, 14, 20, 30, 40] × 10−3 increasing from black to yellow. Assuming Z⊙ = 0.0134 (Asplund et al. 2009), these correspond to Z/Z⊙ = 0.07, 0.15, 0.22, 0.30, 0.45, 0.60, 0.75, 1.0, 1.5, 2.2, 3.0. Ionisation parameter varies along the grid line, increasing with marker size, with values of log U = −4.0, −3.5, −3.0, −2.5, −2.0, −1.5, −1.0, −0.5. Details of the photoionisation models are described in Sect. 4.1.

Figure 7 highlights that photoionisation models driven by SEDs of young stars can reproduce the parameter space populated by our JADES galaxies with only fairly standard assumptions. This indicates we do not need to invoke alternative heating sources (such as AGN or shocks) to explain our findings. We see from this simple comparison that the large scatter observed in our emission line ratio sample likely reflects a diversity in metallicity and ionisation parameter among these galaxies. Our most extreme galaxies (i.e. highest O32) are only reproduced by the models with very low-metallicity (0.07 − 0.15 × Z_⊙) and very high ionisation parameters (−2.0 ≲ log U ≲ −1.0). These ionisation parameters are significantly higher than those typically observed in z ∼ 0 H II regions (Dopita et al. 2000) or z ∼ 2 − 3 galaxies (Sanders et al. 2016). We note that we have not explored the effect of harder radiation fields, which would increase O32 at a fixed ionisation parameter (e.g. Sanders et al. 2016). Indeed, many studies at z ∼ 2 − 3 have argued that harder radiation fields, driven by super-solar α/Fe abundance ratios, are a dominant effect in driving the z ∼ 2 N2–BPT offset (Steidel et al. 2016; Strom et al. 2017; Topping et al. 2020b,a; Cullen et al. 2021). Given the long timescale of iron enrichment via Type Ia supernovae, one would of course expect that these z ≳ 5.5 galaxies would be α-enhanced, and models with non-solar abundance ratios would be more appropriate for parameter estimation.

We can see from these comparisons that our sample likely spans quite a range in metallicity. At lower O32 values, we see a turnover in the models where metallicity and ionisation parameter become highly degenerate. However, model metallicities of at least 0.45 × Z_⊙ are required to approach our highest R2 and R23 values. The R2, R3, and R23 ratios are all known to exhibit double-valued relations (e.g. Curti et al. 2020) in metallicity. We tentatively observe that our highest O32 galaxies tend to have the lowest R23, which reverses the trend seen in SDSS galaxies (grey 2D PDF), suggesting we are sampling the turnover of this relation. Some of our ‘mildest’ galaxies (i.e. highest R2, lowest O32) have R2 ≈ 0.4 and R23 ≈ 1.0 which are broadly consistent with the turnover observed for these ratios in z ∼ 0 galaxies from Curti et al. (2020). This suggests these galaxies are the most chemically evolved systems in our sample and may have metallicities as high as 0.3 − 0.5 × Z_⊙. Indeed, these are also the points that are in best agreement with the Sanders et al. (2023)z ∼ 5.6 composite spectrum. This suggests the shallower CEERS spectroscopy is probing on average more chemically evolved systems than those targeted in JADES.

Overall, it is clear from this simple model comparison that a range of ISM conditions is required to explain the diversity of emission-line ratios observed in this sample. As noted in Sect. 3, our sample largely spans the same regions of emission-line-ratio-space inhabited by low-redshift extreme starbursts such as blueberries and green peas. Thus, the emerging picture is that the galaxies in our sample are likely low-metallicity (∼0.1 − 0.3 × Z_⊙) systems undergoing intense periods of star formation driving strong radiation fields with high ionisation parameters (−2.0 ≲ log U ≲ −1.0) in the ISM.

4.2. High O32: implications for escape fraction and reionisation

The high O32 ratios presented in Fig. 5 could have implications for interpreting the ionising photon escape fraction of these z ∼ 5.5 − 9.5 galaxies. When and how the Universe was reionised remains one of the most important open questions in high-redshift astrophysics. The efficiency with which ionising photons are produced within galaxies and the fraction of these which escape (f_esc) are two important quantities required to understand the contribution of galaxies to this process. Unfortunately, f_esc can never be directly measured within the epoch of reionisation itself because any light escaping the galaxy bluewards of the Lyman-limit will interact with intervening neutral hydrogen. Thus indirect probes of f_esc are highly sought after.

High O32 has been proposed as being an indicator of higher f_esc (e.g. Nakajima & Ouchi 2014; Faisst 2016). The reasoning is that the high O32 ratio selects for highly ionised systems which have a higher likelihood of having density-bounded channels through which ionising photons can escape. However, some studies have shown that O32 does not necessarily correlate well with f_esc; instead finding that results are very dependent on the observed line of sight, and that high O32 can be produced in systems without Lyman continuum (LyC) escape (Paalvast et al. 2018; Bassett et al. 2019; Katz et al. 2020). Nonetheless, samples of Lyman continuum (LyC) leaking galaxies at low-redshift generally show that the fraction of galaxies with high f_esc increases towards higher O32, even if the correlation is not tight (e.g. Izotov et al. 2016; Flury et al. 2022). Indeed, the limited number of known LyC leakers at higher redshift typically have high O32 ratios (Vanzella et al. 2016, 2020, 2022; Fletcher et al. 2019; Saxena et al. 2022). Extreme [O III] emitters at z ∼ 2 − 3 have also been shown to exhibit high O32 ratios ([O III]/[O II] ∼ 10), while also preferentially selecting for high equivalent-width Lyman-α emission (Tang et al. 2019, 2021). This evidence therefore seems to suggest that high O32 is a necessary but not sufficient condition for high f_esc.

We have shown in Sect. 3.4 that our JADES z ∼ 5.5 − 9.5 galaxies have characteristically high O32 ratios, with 21 out of 27 of the sample having [O III]/[O II] > 5 and a median value of 9.4. Owing to the large scatter in O32–f_esc relations, we do not attempt to derive constraints on the f_esc of these galaxies. However, the high fraction of our sample that meets the criterion of high O32 suggests that the fraction of potential LyC leakers might be increasing to high redshift.

5. Conclusions

We have presented rest-frame optical emission line ratio measurements in 27 individual galaxies at 5.5 < z < 9.5 spanning over −17.0 < M₁₅₀₀ < −20.4, observed as part of the ‘Deep’ tier of the JADES program, which represent the deepest integrations yet taken with JWST/NIRSpec. We find that galaxies in this sample occupy regions of line-ratio space that are offset from those inhabited by ‘typical’ galaxies at z ∼ 0 or z ∼ 2, although generally aligned with more extreme low-redshift populations such as ‘blueberry’ and ‘green pea’ dwarf starbursts (Yang et al. 2017a,b).

The high signal-to-noise achieved in individual spectra from our 28 h deep PRISM/CLEAR integrations reveal significant intrinsic scatter in the observed line ratios, with our JADES z ∼ 5.5 − 9.5 galaxies spanning more than an order of magnitude in [O II]/Hβ and [O III]/[O II] ratios. This highlights the diversity in ISM conditions of galaxies at this epoch.

Despite the depth of our observations, we find that low ionisation lines such as [N II] λ6583, [S II] λλ6716, 6731 and [O I] λ6300 are very challenging to detect in these systems. We find no convincing detections of [N II] λ6583 in 7 h deep G395M/F290LP spectra for any of the 18 galaxies in which we have adequate spectral coverage. Stacking these 18 spectra is still unable to deliver a nitrogen detection. However, the 3-σ upper limit suggests our sample median is not offset from the z ∼ 0 N2–BPT sequence contrary to what has been observed in z ∼ 2 samples (e.g. Steidel et al. 2014). In our PRISM spectra we find three modest detections of [S II] λλ6716, 6731, while the tight upper limits we derive on other targets suggests more than an order of magnitude in intrinsic scatter in the S2 ratio within the sample. The fainter [O I] λ6300 line yields only one tentative detection.

The bluer [O II] λλ3726, 3729 and [Ne III] λ3869 are much more readily detected throughout the sample. Measured R2, O32 and Ne3O2 ratios particularly highlight the diversity in this sample; we find median values and standard deviations of R2 = −0.23  ±  0.38, O32 = 0.98  ±  0.36, and Ne3O2 = −0.03  ±  0.27. We find that these ratios are generally consistent with very high ionisation parameters up to log U ∼ −1.0. The sample appears to span a significant range in metallicity of ∼0.1 − 0.3 × Z_⊙, depending on metallicity calibrator assumptions. These ratios are found to be generally consistent with extreme low-redshift populations, particularly green pea and blueberry dwarf starburst samples. This is consistent with the suggestion that these z ∼ 5.5 − 9.5 galaxies are low-mass, low-metallicity galaxies undergoing periods of rapid star formation, driving strong radiation fields.

We compare our sample to composite spectra at z ∼ 5.6 and z ∼ 7.5 from CEERS (Sanders et al. 2023; Tang et al. 2023). We find significant offset between our z ∼ 6 sample median and the z ∼ 5.6 CEERS composite. We instead observe good agreement between this z ∼ 5.6 CEERS composite measurement and the low O32 end of our sample. Assuming these galaxies represent the more evolved, metal-enriched extent of our sample, it is clear that the shallower CEERS spectra are preferentially picking up on the more evolved portion of the z ∼ 6 population. We find no evidence for significant evolution between our z ∼ 6 and z ∼ 8 samples (albeit with a very small sample). Both our z ∼ 6 and z ∼ 8 sample medians are found to be in somewhat closer agreement with the CEERS z ∼ 7.5 composite spectra, which showed large offset from the CEERS z ∼ 5.6 composite spectrum. This could indicate that the fraction of moderately metal-enriched systems increases rapidly from z ∼ 8 to z ∼ 5.5, even if our sample median does not show evidence of significant redshift evolution.

Within its first six months of observations, JWST/NIRSpec has already shed considerable light on the ISM conditions of z ≳ 5.5 galaxies. The diversity of ISM conditions presented in this work highlights the need for observations that span the full range of galaxy properties. This will require a combination of deep observations reaching down to even fainter brightness limits than those presented here, as well as wider area surveys probing rarer, more evolved systems at early epochs. Characterising the ISM conditions in larger samples of galaxies across this broader parameter space will help chart galaxy assembly across the first ∼1 Gyr of cosmic time.

Acknowledgments

A.J.C., A.S., A.J.B., J.C., and G.C.J. 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). F.D.E., R.M., M.C., T.J.L., J.W., L.S., and J.S. acknowledge support by the Science and Technology Facilities Council (STFC), ERC Advanced Grant 695671 “QUENCH”. R.M. also acknowledges funding from a research professorship from the Royal Society. J.W. also acknowledges support from the Fondation MERAC. S.C. acknowledges support by European Union’s HE ERC Starting Grant No. 101040227 – WINGS. E.C.L. acknowledges support of an STFC Webb Fellowship (ST/W001438/1). The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant no.140. S.A. and B.R.D.P. acknowledge support from the research project PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI). R.S. acknowledges support from a STFC Ernest Rutherford Fellowship (ST/S004831/1). H.Ü. gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. D.J.E. is supported as a Simons Investigator and by JWST/NIRCam contract to the University of Arizona, NAS5-02015. E.E., M.R., B.J.D., B.E.R., and F.S. acknowledge support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015. The work 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. R.B. acknowledges support from an STFC Ernest Rutherford Fellowship [grant number ST/T003596/1]. This research is supported in part by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.

References

Appendix A: Table of selected targets

Table A.1 gives an overview of the targets selected in our JADES z ∼ 5.5 − 9.5 sample.

Appendix B: Emission line ratio measurements from stacked spectra

In Section 2.5 we described our procedure for measuring emission-line ratios from stacked spectra. Here we briefly summarise the values obtained from this analysis. Table B.1 outlines the number of galaxies included in each stack and median properties of these samples. Table B.2 provides the dust-corrected emission-line ratios obtained from the resultant stacked spectra. These values are shown as open diamonds in Figures 3–6.

Appendix C: Verification of flux calibration

As a verification of our spectral flux calibration, we compared the fluxes obtained from NIRCam photometry with those obtained from the PRISM/CLEAR spectra. Using imaging data from the JADES programme (Eisenstein et al. 2023; Rieke et al. 2023), we extract apodised fluxes for each filter by extracting the flux within the projected shutter position of each target. Assuming the MSA pointing uncertainty is negligible, this should represent the total amount of flux in each spectral region incident on the shutter area. However, the point spread functions (PSFs) of NIRCam and NIRSpec are different, and in reality the final flux recorded by NIRSpec is additionally reduced by losses due to diffraction caused by the shutter itself. Our data reduction pipeline includes a step which corrects for these slitlosses, meaning that the output spectrum should be scaled to the same flux density as the total magnitude of the source (correcting for position within the slit and flux falling outside the slit). As a further check on the relative flux calibration between NIRCam and NIRSpec, we also consider the flux from the NIRCam imaging which would fall within the NIRSpec microshutter aperture. To make a fair comparison, we first remove the slitloss correction applied by the pipeline to the NIRSpec spectrum (to consider just the uncorrected flux within the shutter) and convolve the resultant PRISM/CLEAR spectrum with the NIRCam filter transmission curves.

Figure C.1 shows the comparison between the NIRCam and NIRSpec fluxes derived in this way, with grey points representing each individual pointing⁹ of each target galaxy. Typical statistical uncertainties for individual measurements are shown by the horizontal grey marks, which include the contributions from both the NIRCam and the NIRSpec uncertainties. For a handful of galaxies, the NIRSpec continuum flux was low enough that the final measured flux did not have S/N > 3, and we removed these this comparison. From the large red points and errorbars, showing the median and standard deviation across all observations in a particular filter, we see that for all filters the median of the comparison points is well within one standard deviation of unity. However, on average many of these indicate that the NIRCam fluxes are up to ∼10% lower than the NIRSpec fluxes. In particular, we see evidence of a subtle trend with the NIRCam-to-NIRSpec flux ratio increasing slightly to long wavelengths. This arises as a result of the NIRSpec PSF degrading more significantly with increasing wavelength beyond 2.5 μm.

Fig. C.1. Ratio of apodised broad-band fluxes extracted from the projected MSA shutter regions on the NIRCam imaging compared to fluxes obtained from convolving PRISM/CLEAR spectra with the NIRCam filer transmission curves. Refer to the text for more details. Grey points show values derived for each independent pointing of each galaxy targeted in this work. Horizontal grey marks indicate the typical 1-σ measurement uncertainty on the individual measurements in each filter (which encompasses uncertainties from both NIRCam and NIRSpec). Red diamonds and errorbars show the median and standard deviation for each filter.

Fig. C.2. Comparison of Hα flux measurements made from low-resolution PRISM/CLEAR and R ∼ 1000 G395M/F290LP measurements for galaxies in this sample. Fluxes reported in the R ∼ 1000 data are systematically higher by a median factor of 1.16.

Our NIRSpec data reduction pipeline does take this effect into account, in the form of the slitloss correction (which was undone for Figure C.2). These slitlosses will depend on the morphology of the galaxy and its location within the shutter, but the ∼10% drift seen with wavelength in the comparison in Figure C.1 is entirely consistent with the differing PSFs of NIRCam and NIRSpec as a function of wavelength. Slitlosses presented in this work were calculated assuming each galaxy is a point source at the location of the centroid derived from imaging. Table C.1 gives the half-light radii from the JADES public catalogue (Rieke et al. 2023)¹⁰ for the subset of our sample for which NIRCam imaging was available. In this study, we focused only on rest-frame optical emission lines at z > 5.5, meaning that the observed-frame wavelengths were all λ ≳ 2.5 μm. At these wavelengths, the NIRCam PSF is typically worse than ( in F277W, in F356W¹¹). Thus, at the wavelengths considered in this work, the assumption of a point source profile should be reasonable for all but a few of our largest targets.

The majority of the emission-line ratios studied in this work have relatively short wavelength baselines, for which the precise slitloss corrections will make a negligible difference. The O32 and R23 ratios are the notable exceptions to this. However, even before slitloss corrections are applied, Figure C.2 demonstrates that, on average, there is only a ∼10% effect across the full wavelength range. Given that [O II] λλ3726,3729 and [O III] λ5007 are only separated by ∼1.4 μm in the observed frame, the effect would be even smaller than this quoted 10%. Furthermore, the point-source slitloss corrections we do apply should be adequate for most sources. Thus, the final flux calibration uncertainty on the O32 ratio is likely well below 10% which is smaller than the quoted measurements uncertainties throughout this paper.

C.1. Comparison of grating flux calibration

This paper almost exclusively makes use of the PRISM/CLEAR mode. This is because, at redshifts beyond z > 5, the rest-frame optical emission lines fall beyond λ_obs ≳ 2.5 μm where the PRISM/CLEAR spectral resolution is R ∼ 200 − 300, meaning most of these lines are well resolved. The notable exception to this is the [N II] λ6583/Hα ratio, for which we require the R ∼ 1000 G395M/F290LP grating. Although reported ratios are only ever measured within the same spectral observation, we briefly consider the relative flux calibration of the PRISM/CLEAR and G395M/F290LP. Figure C.1 shows the ratio of Hα fluxes measured in the G395M/F290LP to those measured in the PRISM/CLEAR mode for galaxies in this sample which have spectral coverage in both settings. We note that the fluxes measured here are not those reported in the main text, where the continuum is modelled with stellar templates. Rather the continuum here is modelled with a simple spline fit. This is done to minimise any model-dependent effects in this comparison. The dashed magenta line shows the median value of 1.16, indicating that fluxes measured from the grating are systematically larger. This effect has been noted in Bunker et al. (2023). We note that in the PRISM/CLEAR, the reported Hα flux would be blended with [N II], whereas these lines are resolved in the G395M/F290LP measurements. However, as pointed out in Section 3, [N II] λ6583 is not clearly detected in any of these galaxies, and this would contribute at the < 10% level. Furthermore, the effect of this blending would be to increase the PRISM/CLEAR flux, but it is the G395M/F290LP flux that is determined to be larger. As outlined above, the flux calibration in the PRISM/CLEAR mode is in good agreement with broad-band photometry. Thus, we consider that the PRISM/CLEAR flux calibration is likely reliable, while the grating flux may be overestimated. However, given we only make use of the G395M/F290LP data to measure the ratio of two lines very close in wavelength within the same observation, the systematically high fluxes recorded should have no effect on any of the findings reported here.

All Tables

All Figures

Fig. 1. Example spectra from five galaxies included in our sample. Left: low-resolution PRISM/CLEAR spectra shifted to the rest-frame according to the observed redshift. Flux is shown in fλ but has been renormalised relative to the peak flux from the [O III] λ5007 line. Solid horizontal lines show the zero flux of each spectrum. The shaded region around the ‘zero’ line indicates the 1-σ noise spectrum in the same renormalised units. We show only the subset of spectral coverage that includes the key emission lines used in this study. Vertical dotted lines show the expected centroids of the [O II] λ3727, [Ne III] λ3869, [Ne III] λ3967, Hδ, Hγ, [O III] λ4363, Hβ, [O III] λ4959, [O III] λ5007, [O I] λ6300, Hα, [N II] λ6583, [S II] λ6716, and [S II] λ6731 lines respectively (moving left to right). Right: zoom in on the Hα + [N II] complex as observed in the G395M/F290LP grating spectra of these same galaxies. Observed spectra are shown in the same way as for the left panel (except renormalised to Hα, rather than [O III]). Vertical dotted lines show the expected centroids of Hα and [N II] λ6583.

In the text

Fig. 2. Redshift and magnitude distribution of the 27 galaxies in our JADES sample. Top: redshift histogram of galaxies observed in this study. We divide the sample into two sub-samples ‘z∼6’ and ‘z ∼ 8’, based simply on cutting at z = 7 since this is the redshift at which Hα is redshifted beyond the spectral coverage of NIRSpec. Thus, for JADES ‘z ∼ 8’ we only have coverage of lines from [O II] λλ3726, 3729 to [O III] λ5007, while for JADES ‘z ∼ 6’ we can in principle observe all the rest-frame optical lines out to [S II] λλ6716, 6731. The orange hatched histogram indicates galaxies for which we can constrain the [S II] λλ6716, 6731 flux. Green hatching is the same but for [N II] λ6583. Bottom: redshift vs. M1500 for the full z ∼ 5.5 − 9.5 sample. We derived M1500 directly from the PRISM spectra as described in Sect. 2.4.1.

In the text

Fig. 3. BPT and VO87 diagrams showing the comparison between JADES z ∼ 6 galaxies and samples across various redshifts. Individual measurements are shown as solid purple points, while values derived from stacked JADES spectra are shown as the open magenta diamonds. The background grey 2D PDF shows z ∼ 0 galaxies from SDSS. Extremely metal-poor z ∼ 0 galaxies from Izotov et al. (2018, 2019) are shown as yellow pentagons. MOSDEF z ∼ 2 galaxies are shown as green circles (Kriek et al. 2015). The z ∼ 5.6 composite CEERS spectra presented by Sanders et al. (2023) are shown as the blue open circles. The solid navy lines in each panel show the theoretical maximum starbust lines from Kewley et al. (2001), while the dotted line in the left panel shows the empirical demarcation derived by Kauffmann et al. (2003). Left: classical BPT or ‘N2–BPT’ – Purple triangles provide 3-σ upper limits on the locations of our JADES galaxies since we are unable to detect [N II] λ6583 in any of the 18 galaxies for which we have G395M/F290LP spectral coverage of the Hα complex, although the R3 ratio is very well constrained in all cases from PRISM data. Even after stacking the grating spectra of all galaxies, we do not recover a robust detection of [N II] λ6583. Centre: in the S2–VO87 diagram (often referred to as ‘S2–BPT’), three out of 20 of our JADES galaxies show detections of [S II] in individual spectra. The tight 3-σ upper limits we place on our non-detections highlight that there must be more than an order of magnitude in scatter in the S2 ratio within the sample. Right: we find one tentative detection of [O I] λ6300 in our O1–VO87 diagram.

In the text

Fig. 4. R2–R3 diagram showing our two JADES sub-samples at z ∼ 6 (purple points) and z ∼ 8 (maroon points). Comparison samples from SDSS and MOSDEF are shown as in Fig. 3. We additionally show comparison with two populations of z ∼ 0 extreme starbursts: green peas (Yang et al. 2017a; orange plusses) and blueberries (Yang et al. 2017b; blue crosses). The three well-studied z > 7.5 galaxies from the SMACS 0723 ERO observations are shown as lime green pentagons, adopting the ratio values presented by Nakajima et al. (2023). Unlike the BPT and VO87 diagrams, the majority of our sample have > 3σ detections of all lines in individual spectra, revealing a large diversity within the sample. Ratios measured from stacked JADES spectra are shown as the open magenta and red diamonds (see Sect. 2.5).

In the text

Fig. 5. R23–O32 diagram, probing the ionisation and excitation of the ISM. The JADES sub-samples at z ∼ 6 (purple points) and z ∼ 8 (maroon points) lie at generally high excitation an high ionisation. Comparison samples, including SDSS, MOSDEF, green peas, blueberries, and SMACS 0723, are shown as in Fig. 4. z ∼ 4.5 − 8 galaxies from GLASS are shown as grey diamond (Mascia et al. 2023). Blue and orange open circles show the z ∼ 5.6 and z ∼ 7.5 composite CEERS spectra from Sanders et al. (2023), while the z ∼ 7.7 composite CEERS spectrum from Tang et al. (2023) is shown as the green open hexagon. Stacked JADES spectra are shown as the open magenta and red diamonds.

In the text

Fig. 6. R23–Ne3O2 diagram, probing the ionisation and excitation of the ISM. R23–Ne3O2 closely resembles the R23–O32 diagram with the additional advantage that the [Ne III] λ3869/[O II] λλ3726,329 ratio is less sensitive to the applied dust-correction, or any wavelength-dependent flux calibrations. JADES galaxies are plotted as in Fig. 5. Fewer literature samples present [Ne III] ratios so we show only SDSS, blueberries, MOSDEF, and the z ∼ 7.7 composite CEERS spectrum from Tang et al. (2023). These are each plotted in the same way as in previous figures.

In the text

Fig. 7. JADES z ∼ 5.5 − 9.5 galaxies compared to photoionisation models. Measured ratios are shown as purple points, as in Figs. 3–6 except that all galaxies are plotted in the same colour, regardless of redshift. SDSS galaxies are shown in the same manner as Figs. 3–6. Each line shows a set of models with constant metallicity with values of Z = [1, 2, 3, 4, 6, 8, 10, 14, 20, 30, 40] × 10−3 increasing from black to yellow. Assuming Z⊙ = 0.0134 (Asplund et al. 2009), these correspond to Z/Z⊙ = 0.07, 0.15, 0.22, 0.30, 0.45, 0.60, 0.75, 1.0, 1.5, 2.2, 3.0. Ionisation parameter varies along the grid line, increasing with marker size, with values of log U = −4.0, −3.5, −3.0, −2.5, −2.0, −1.5, −1.0, −0.5. Details of the photoionisation models are described in Sect. 4.1.

In the text

Fig. C.1. Ratio of apodised broad-band fluxes extracted from the projected MSA shutter regions on the NIRCam imaging compared to fluxes obtained from convolving PRISM/CLEAR spectra with the NIRCam filer transmission curves. Refer to the text for more details. Grey points show values derived for each independent pointing of each galaxy targeted in this work. Horizontal grey marks indicate the typical 1-σ measurement uncertainty on the individual measurements in each filter (which encompasses uncertainties from both NIRCam and NIRSpec). Red diamonds and errorbars show the median and standard deviation for each filter.

In the text

	Fig. C.2. Comparison of Hα flux measurements made from low-resolution PRISM/CLEAR and R ∼ 1000 G395M/F290LP measurements for galaxies in this sample. Fluxes reported in the R ∼ 1000 data are systematically higher by a median factor of 1.16.
In the text