# The chemical enrichment in the early Universe as probed by JWST via direct metallicity measurements at z 8

Mirko Curti, Francesco d'Eugenio, Stefano Carniani, Roberto Maiolino, Lester Sandles, Joris Witstok, William M Baker, Jake S. Bennett, Monika Joanna Piotrowska, Sandro Tacchella, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mirko Curti, Francesco d'Eugenio, Stefano Carniani, Roberto Maiolino, Lester Sandles, et al. The chemical enrichment in the early Universe as probed by JWST via direct metallicity measurements at z 8. Monthly Notices of the Royal Astronomical Society, 2023, 518, pp.425-438. 10.1093/mn-ras/stac2737. hal-03839214

HAL Id: hal-03839214

https://hal.science/hal-03839214

Submitted on 7 Nov 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# The chemical enrichment in the early Universe as probed by JWST via direct metallicity measurements at $z \sim 8$

Mirko Curti,¹,²* Francesco D'Eugenio,¹,² Stefano Carniani,³ Roberto Maiolino,¹,²,⁴

Lester Sandles,¹,² Joris Witstok,¹,² William M. Baker,¹,² Jake S. Bennett,¹,⁵

Joanna M. Piotrowska,¹,² Sandro Tacchella,¹,² Stephane Charlot,⁶ Kimihiko Nakajima,⁷

Gabriel Maheson,¹,² Filippo Mannucci,⁸ Amirnezam Amiri,⁹,⁸ Santiago Arribas,¹⁰

Francesco Belfiore,⁸ Nina R. Bonaventura,¹¹ Andrew J. Bunker,¹² Jacopo Chevallard,¹²

Giovanni Cresci,⁸ Emma Curtis-Lake,¹³ Connor Hayden-Pawson,¹,², Gareth C. Jones,¹²

Nimisha Kumari,¹⁴ Isaac Laseter,¹⁵ Tobias J. Looser,¹,² Alessandro Marconi,⁹,⁸

Michael V. Maseda,¹⁵ Jan Scholtz,¹,² Renske Smit,¹⁶ Hannah Übler,¹,² and Imaan E. B. Wallace,¹²

Affiliations are listed at the end of the paper

Accepted 2022 September 22. Received 2022 September 22; in original form 2022 July 26

## ABSTRACT

We analyse the chemical properties of three $z \sim 8$ galaxies behind the galaxy cluster SMACS J0723.3–7327, observed as part of the Early Release Observations programme of the James Webb Space Telescope (JWST). Exploiting [O III]λ4363 auroral line detections in NIRSpec spectra, we robustly apply the direct $T_e$ method for the very first time at such high redshift, measuring metallicities ranging from extremely metal poor (12+log(O/H)≈ 7) to about one-third solar. We also discuss the excitation properties of these sources, and compare them with local strong-line metallicity calibrations. We find that none of the considered diagnostics match simultaneously the observed relations between metallicity and strong-line ratios for the three sources, implying that a proper re-assessment of the calibrations may be needed at these redshifts. On the mass-metallicity plane, the two galaxies at $z \sim 7.6$ (log($M_*/M_\odot$) = 8.1, 8.7) have metallicities that are consistent with the extrapolation of the mass-metallicity relation at $z \sim 2–3$, while the least massive galaxy at $z \sim 8.5$ (log($M_*/M_\odot$) = 7.8) shows instead a significantly lower metallicity. The three galaxies show different level of offset relative to the Fundamental Metallicity Relation, with two of them (at $z \sim 7.6$) being marginally consistent, while the $z \sim 8.5$ source deviating significantly, being probably far from the smooth equilibrium between gas flows, star formation and metal enrichment in place at later epochs.

**Key words:** Galaxies: ISM, galaxies: evolution, galaxies: general, galaxies: abundances

## 1 INTRODUCTION

The study of the abundances of heavy elements (the ‘metallicity’) in the interstellar medium (ISM) of galaxies provides precious insights on the physical processes responsible for their formation, and on how the relative importance of such processes has changed across cosmic time (e.g. Ma et al. 2016; Davé et al. 2017; Torrey et al. 2019; Langan et al. 2020, see also the review by Maiolino & Mannucci 2019).

In the local Universe, the metallicity of the gas-phase is observed to correlate tightly with some of the primary galactic properties. Specifically, the metallicity depends on the stellar mass of the galaxy (mass-metallicity relation, MZR, Tremonti et al. 2004) and has a secondary (inverse) dependence on the star formation rate (SFR, Ellison et al. 2008). This metallicity-$M_*$-SFR relation has been called the Fundamental Metallicity Relation (FMR, Mannucci et al. 2010; Mannucci et al. 2011; Curti et al. 2020a), a signature of the smooth,

long-lasting equilibrium between gas flows and secular evolution (e.g. Bouché et al. 2010; Lilly et al. 2013).

In the past fifteen years, major efforts have been dedicated to explore the metallicity of high-redshift galaxies through massive spectroscopic surveys in the near-IR (e.g. Erb et al. 2016; Maiolino et al. 2008; Mannucci et al. 2009; Zahid et al. 2011; Steidel et al. 2014; Zahid et al. 2014; Strom et al. 2017; Curti et al. 2020b; Sanders et al. 2021, among many others). These works have established that the mass-metallicity relation is already in place up to $z \sim 3.3$, and that it evolves with redshift, in the sense that higher redshift galaxies appear metal deficient compared to galaxies of similar stellar mass at later epochs. However, when taking into account the secondary dependence of the MZR on the SFR, as described by the FMR framework, any evolution is apparently canceled out at least up to $z \sim 3.3$ (e.g. Mannucci et al. 2010; Cresci et al. 2019; Sanders et al. 2021). This finding implies that galaxies up to $z \sim 3.3$ follow, on average, the same smooth evolution as local galaxies, and that the evolution of the MZR is a consequence of $z \sim 2–3$ galaxies

<footnote>* E-mail: mc2041@cam.ac.uk</footnote>
having higher SFR (living at the epoch of the peak of the cosmic
star formation rate density, Madau & Dickinson 2014), therefore
populating preferentially the low-metallicity region of the same, non-
evolving FMR. Together with the existence of the main sequence of
star formation (e.g. Noeske et al. 2007; Sandles et al. 2022; Popesso
et al. 2022, and references therein) this is an indication that secular
rather than stochastic processes dominate the evolution of galaxies
at these epochs.

Pushing the exploration of galaxy metallicities to redshifts higher
than ~ 3.3 has proven very challenging so far due to the intrinsic tech-
nical limitations of current astronomical facilities, as the primary (op-
tical) nebular diagnostics required to measure the gas-phase metal-
licity from galaxy spectra are redshifted outside the wavebands ob-
servable from the ground (but see Troncoso et al. 2014a; Shapley
et al. 2017; Witstok et al. 2021).

An additional problem is that metallicity determinations gener-
ally rely on adopting locally calibrated strong-line diagnostics, as
temperature-sensitive auroral lines are usually too faint to be de-
tected, preventing the use of the more robust and ‘direct’ electron-
temperature (T_e) method. However, it is not clear whether the strong-
line calibrations derived in the local universe are valid for distant
galaxies, which have drastically different properties (Maiolino &
Mannucci 2019). To date, auroral lines have been observed only
in a small number of intermediate-redshift (z~1-3) galaxies, and
constitute mostly serendipitous, low-significance detections (e.g.
Christensen et al. 2012; Patrício et al. 2018; Sanders et al. 2020)

This observational landscape is set to be revolutionised by the
advent of the *James Webb Space Telescope (JWST)* and its near-IR
spectrograph NIRSpec (Jakobsen et al. 2022; Ferruit et al. 2022),
which has opened the capability of obtaining multi-object spectro-
scopic observations in the near-IR from space (up to 5.3µm), com-
bining high multiplexing with a sensitivity much higher than any
current and past facility.

The first JWST/NIRSpec spectra have been recently released within the context of the Early Release Observations, and were obtained by targeting galaxies lensed by the cluster SMACS J0723.3-7327 (Ebeling et al. 2001, 2007, 2010; Mann & Ebeling 2012; Ebeling et al. 2013; Repp et al. 2016; Repp & Ebeling 2018), providing clear detections of nebular lines in galaxies out to z =8.5, which can be used to characterise the chemical enrichment of the ISM at the Epoch of Reionisation (EoR), and beyond. Most remarkably, some of these spectra present clear detections of the [OIII]4363 auroral line, which can be used to robustly measure the gas metallicity with the T<sub>e</sub> method for the very first time at such high redshifts.

In this paper, we extract the direct metallicity of such galaxies at
z ~ 8, which have individual detections of auroral lines, and test
whether they follow the same strong-line metallicity calibration as
local galaxies. In addition, we investigate whether these galaxies,
probing the EoR, follow the same metallicity scaling relations as
lower-redshift galaxies or exhibit evidence for evolution.

This paper is organised as follows. In Section 2 we describe obser-
vations, data processing and data analysis. In particular, we discuss
the different steps we have implemented to properly process and
calibrate the data by using publicly available information. We also
describe our emission lines fitting procedure, and how SED fitting
was performed on NIRCAM photometry to derive different physical
properties of these galaxies. In Section 3, we discuss the emission-
line properties of our targets, in terms of Balmer decrements and dust
attenuation, excitation-diagnostic diagrams, and metallicity proper-
ties as derived with the T<sub>e</sub> method. We also compare the observed line
ratios and metallicity with some of the most widely adopted strong-

line abundance diagnostics calibrated on local galaxies. In Section 4,
we discuss the observed metallicity properties of these targets in the
framework of the most relevant metallicity scaling relations, i.e. the
mass-metallicity and the fundamental metallicity relation. Finally, in
Section 5 we report our conclusions.

Throughout this work, we assume a standard lambda cold dark
matter cosmology with H₀=70 km s⁻<fnref n="1" />, Ω_m=0.3, and Ω_Λ=0.7.

## 2 OBSERVATIONS, DATA PROCESSING AND DATA ANALYSIS

2.1 NIRSpec spectroscopy

2.1.1 NIRSpec observations

We use publicly available data from the Early Release Observations of
the cluster SMACS J0723.3-7327 (Program ID: 2736, Pontoppidan
et al. 2022). The NIRSpec observations consist of two pointings
with the same Multi Shutter Array (MSA) configuration but differ-
ent acquisition stars (and different filters for the acquisition); in the
following, we refer to these two observations as Obs 7 and Obs 8.

NIRSpec observations were carried out by using the disperser-filter
combinations G235M/F170LP and G395M/F290LP, which cover the
wavelength range between 1.66µm and 5.16µm and provide spectra
with a spectral resolution of *R* ~ 1000. For each observation, three
nodding positions of 20 groups and two integrations each were per-
formed for each grating setup. The total exposure time of the two
individual pointings is the same and corresponds to 8,840 seconds
for each grating. While shutters were opened on 35 targets, in this
paper we only focus on the three targets with the highest redshift:
ID 4590 (z=8.4953), ID 6355 (z=7.6643) and ID 10612 (z=7.6592).
These redshifts are based on emission-line velocities measured on
the JWST spectra (see Sect. 2.1.3).

2.1.2 NIRSpec data reduction

We have retrieved the level 2 data (i.e. count rate maps) from the
MAST archive, but we have then reprocessed the data through the
GTO pipeline (NIRSpec/GTO collaboration, in prep.). As most of
the processing steps use the same algorithms that the pipeline used
to generate the MAST archive products (see Fig. 11 and section 4.3
of Ferruit et al. 2022), the resulting 2D spectra are not very different
from those provided by the standard MAST pipeline; however, we
perform our own optimised extraction aperture and bad/cosmic pixel
flagging and masking. Most importantly, to mitigate the fact that
many of the calibration reference files used by the pipeline may still
correspond to ground data or be model based<fnref n="1" />, we have generated a
correction to the response function by analysing the observation of the
spectrophotometric calibration star 2MASS J18083474+6927286,
observed during commissioning, and publicly available through pro-
gramme JWST 1128. We note that this star has been subsequently
removed from the list of primary calibrators for JWST, on account of
its variability (σ = 0.41 per cent, Gordon et al. 2022); however, for
our purposes, this adds a negligible amount of uncertainty. Specific-
ally, the resulting spectrum of the star processed in the same way as
the scientific targets was compared with its intrinsic spectrum from
STScI (Bohlin et al. 2014, 2020, astronomical catalogue CALSPEC)
to derive a more accurate flux calibration.

Given that our sources are marginally resolved, we consider both

<footnote n="1">¹ https://jwst-docs.stsci.edu/jwst-calibration-pipeline-caveats</footnote>
the assumption of point-like source and extended source to correct for path-losses. However, as the sources are fairly compact (especially compared to the point spread function, PSF, at long wavelengths), we use the point-source assumption for our fiducial analysis. We note that at long wavelengths (F290LP filter, with PSF comparable or larger than the shutter size), which are of greatest interest in this paper, the relative flux correction does not strongly depend on which of the two cases is adopted, while providing an offset of ~2 on the absolute calibration; yet, we verified that adopting the extended-source corrections does not significantly affect the main results of the paper. Background subtraction was performed through the standard technique of subtracting the average of two nodding positions from the other position.

We also inspect the exposures of the individual nods in order to identify artefacts that might have escaped automated flagging and other potential issues. We note that the auroral lines are detected in the individual exposures, hence confirming their detection at high confidence level, with the exception of galaxy ID 4590 (z=8.5), for which no emission lines are detected in one of the nods of Obs 7. By inspecting the 2D images before background subtraction, it was found that, although not identified as faulty, shutter [3, 27, 167], on which source ID 4590 was nodded, did not open in Obs 7. As a consequence, for source ID 4590, we use Obs 8 and only the two noddings of Obs 7 for which the source is in open shutters.

Stacking was performed with the GTO pipeline, taking into account the variances and quality arrays of both observations. The resulting spectra are shown in Fig.1. Note that these spectra are in $F_{\lambda}$ (while the spectra provided by MAST archive are in $F_{V}$). These spectra can be retrieved from a publicly accessible repository<fnref n="2" />

### 2.1.3 Spectral fitting

In our analysis, we use only data from the G395M/F290LP disperser/filter combination, because it covers the entire rest-frame wavelength range relevant to this work. We measure line fluxes using *ppxf*, the penalised pixel fitting algorithm of Cappellari (2017). We rebin the spectra to a regular grid in linear velocity space, using the original velocity sampling as uniform pixel size (102 km s⁻¹)<fnref n="3" />. To model the continuum, we use a library of simple stellar-population (SSP) spectra coupled with a 10th-order multiplicative Legendre polynomial. The SSP spectra are the high-resolution (R=10,000) synthetic spectra from the C3K library (Conroy et al. 2019) with MIST isochrones (Choi et al. 2016) and solar abundances. Note that repeating our fit with a polynomial background does not change our results. The emission lines and continuum are tied together to have the same velocity offset and velocity dispersion (while *ppxf* also takes into account the variable spectral resolution). For [OII]λλ3727, 3729, we constrain the doublet ratio to its physical range; for [NeIII]λλ3870, 3969, we fixed the ratio to 0.31. As a safety check, for [OIII]λλ4959, 5007, we leave the ratio a free parameter and check that the recovered value is consistent within the errors with the theoretical expectation of 0.34 (e.g. Osterbrock & Ferland 2006). We perform two fits: after the initial fit, we reject spectral pixels beyond four standard deviations from the best-fit model, then repeat the procedure with outliers masked to obtain the final measurements.

The resulting fitted lines are shown in Fig. 1 (with zoomed insets for lines most relevant for this work) and the resulting line fluxes

relative to $H\beta$ are given in Table 1. The uncertainties are taken from *ppxf* (up-scaled by $\sqrt{\chi^2}$), but we find them consistent with the values obtained from running one hundred Monte Carlo realisations of the best-fit spectrum, with random Gaussian noise taken from the noise vector.

## 2.2 NIRCam imaging

### 2.2.1 NIRCam observations

We use the deep NIRCam imaging data of SMACS J0723.3–7327 from the Early Release Observations (Programme ID 2736) in the F090W, F150W, F200W, F277W, F356W, and F444W filters, which cover an observed wavelength range of $\lambda_{\text{obs}} = 0.8–5\,\mu\text{m}$. We reduce the raw level-1 data products with the public JWST pipeline (v1.6.1),<fnref n="4" /> using the latest available calibration files (CRDS_CTX=jwst_0927.pmap). An additional background subtraction is performed on the final mosaiced images by using PHOTUTILS to mask pixels identified with sources and then measuring the background with PHOTUTILS.BACKGROUND2D.

We conduct aperture photometry on the final mosaics with a range of aperture sizes, given the extended morphology of the three galaxies in this study and their close-by neighbours. Specifically, we use circular apertures with radii from 0.2 – 0.4 arcsec. We perform an aperture correction with the help of WebbPSF (Perrin et al. 2015).<fnref n="5" /> We find that the final SFRs and stellar masses only weakly depend on the aperture size. For simplicity, we therefore assume a fiducial size of 0.3 arcsec. We estimated the uncertainties of these fluxes from the error maps of the mosaic images.

### 2.2.2 Spectral energy distribution fitting

We perform SED modeling of the NIRCam photometry with the Bayesian code BEAGLE (Chevallard & Charlot 2016) with the main aim to derive stellar masses and SFRs for our three galaxies. For consistency with the local MZR and FMR considered here (discussed in Section 4), we use a Chabrier (2003) stellar initial mass function (IMF).<fnref n="6" /> For the star-formation history (SFH), we assume a delayed-exponential form. We find for all three galaxies stellar masses and SFRs of $10^{8–9} M_{\odot}$ and $16–65 M_{\odot}$ yr⁻¹ (Table 2). Furthermore, the inferred stellar (and gas-phase<fnref n="7" />) metallicities are significantly sub-solar (12+log(O/H)=7.25 ± 0.21, 7.53 ± 0.08, and 7.50 ± 0.12) and an overall low attenuation ($A_{V}$=0.37 ± 0.04, 0.50 ± 0.03 and 0.21 ± 0.03) for ID4590, ID6355 and ID10612, respectively. It is interesting that we obtain broad-band metallicities fairly consistent with the $T_{e}$ ones but for ID6355 (which is underestimated by a factor of ~5), confirming in particular the very low metallicity of ID4590 (see section 3.3).

Within the BEAGLE framework, we also explore a constant SFH, which has however little effect on the inferred stellar masses and SFRs. Finally, we explore two other SED fitting codes, namely BAGPIPES (Carnall et al. 2018) and PROSPECTOR (Johnson et al. 2021). We run PROSPECTOR with the same setup as in Tacchella et al. (2022a).

<footnote n="4">⁴ Available at https://github.com/spacetelescope/jwst.</footnote>

<footnote n="5">⁵ Available at https://github.com/spacetelescope/webbpsf.</footnote>

<footnote n="6">⁶ Arguably, this IMF may not be appropriate for these young, low-metallicity galaxies. However, it is not established what the correct IMF at such early times should be and an IMF conversion only introduces a constant scaling factor.</footnote>

<footnote n="7">⁷ In SED modeling, we assume that the stellar and gas-phase metallicity are the same.</footnote>

<footnote n="2">² Available at https://doi.org/10.5281/zenodo.6940561.</footnote>

<footnote n="3">³ We verified that e.g. the $H\delta/H\beta$ line ratios are on average unchanged when fitting the data without rebinning.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** JWST/NIRSpec spectra of the three sources analysed in this work. Several Hydrogen, Helium and metals emission lines are detected; the best-fit models are highlighted in blue. The data and best-fit continuum are traced by the solid grey and red lines, respectively; the dots are the residuals. The bottom panels show a zoom-in on the spectral region of three groups of lines; from left to right they are: the [O II]λ3727, 29 doublet, Hδ, Hγ and [O III]4363, and Hβ and [O III]λ5007. Grey regions have been masked due to artefacts (either in the spectrum or in the noise) or sigma clipping; the wide grey region at 4450-4750 Å in the top panel of ID 4590 falls in the detector gap.</figcaption></figure>

With BAGPIPES, we explored three alternative SFHs: either constant, delayed exponential or delayed exponential with burst. Reassuringly, we find that the derived stellar masses and SFRs are consistent within 1σ uncertainties between the different codes, as well as with those obtained by Tacchella et al. (2022b) via a PROSPECTOR run including emission lines and assuming a bursty prior on the SFH. Moreover, the

stellar masses (converted to the same IMF) agree within 1σ uncer-
tainties with the results presented by Carnall et al. (2022). However, a
more thorough comparison between these different SED fitting codes
goes beyond the scope of the present paper.
<figure><img src="image_4.png" /><figcaption>Figure 1. (continued)</figcaption></figure>

<table><caption>**Table 1.** Observed line intensities normalised to Hβ = 1 and errors on the ratio. The Hβ emission line fluxes (in physical units of 10⁻¹⁸ erg s⁻¹ cm⁻²) measured from the spectra are as follows: ID 4590 =1.54 ± 0.06; ID 6355 =2.11 ± 0.05; ID 10612 =1.19 ± 0.04.</caption><thead><tr><th>Galaxy ID</th><th>[O III]λλ3727, 29</th><th>[Ne III]λ3869</th><th>Hζ</th><th>Hε</th><th>Hδ</th><th>Hγ</th><th>[O III]λ4363</th><th>[O III]λ4959</th><th>[O III]λ5007</th></tr></thead><tbody><tr><td>4590</td><td>0.15 ± 0.07</td><td>0.17 ± 0.03</td><td>0.13 ± 0.02</td><td>0.13 ± 0.02</td><td>0.16 ± 0.02</td><td>0.41 ± 0.03</td><td>0.14 ± 0.02</td><td>0.96 ± 0.05</td><td>3.08 ± 0.13</td></tr><tr><td>6355</td><td>0.90 ± 0.10</td><td>0.45 ± 0.03</td><td>0.18 ± 0.02</td><td>0.18 ± 0.02</td><td>0.26 ± 0.02</td><td>0.46 ± 0.02</td><td>0.10 ± 0.02</td><td>2.66 ± 0.07</td><td>8.29 ± 0.21</td></tr><tr><td>10612</td><td>0.26 ± 0.08</td><td>0.50 ± 0.04</td><td>0.19 ± 0.03</td><td>0.14 ± 0.03</td><td>0.17 ± 0.03</td><td>0.44 ± 0.03</td><td>0.18 ± 0.03</td><td>2.34 ± 0.09</td><td>7.11 ± 0.24</td></tr></tbody></table>

### 2.2.3 Magnification factors

The galaxies analysed in the present paper are background galaxies of the SMACS J0723.3–7327 lensing cluster with publicly available lens models. In this work, we exploit the models recently provided by Mahler et al. (2022), which combine ancillary HST with novel JWST/NIRCAM data to better constrain the cluster mass distribution. Magnification maps are derived for each target redshift, and the average value from within a 1''-wide box around the central coordinates of each galaxy is assumed as the fiducial magnification factor. The associated uncertainty is computed from the standard deviation of one hundred realisations of the magnification maps obtained through Monte Carlo simulations by perturbing the input model parameters. These results in the magnification factors listed in Table 2.

Additional models for this cluster have been recently published by Pascale et al. (2022) and Caminha et al. (2022). A full assessment of the systematics introduced by the adoption of different lensing models is far beyond the scope of this work. However we remark that, for the highest magnification galaxy ID4590, the different values obtained by the available models can impact the stellar mass and SFR determination for this source up to a factor of 30 – 50%.

## 3 LINE RATIOS AND ABUNDANCES AT Z ~ 8

### 3.1 Balmer decrements and dust attenuation

Fig. 2 shows a comparison between the observed and theoretical fluxes relative to Hβ for different Balmer emission lines, namely Hζ, Hε, Hδ, and Hγ. The theoretical values of the corresponding Balmer decrements computed using PyNeb (Luridiana et al. 2012, 2015), assuming Case B recombination, an electronic temperature $T_e = 1.5 \times 10^4$ K and an electron density $N_e = 300$ cm⁻³ (i.e., Hζ/Hβ ≈ 0.105, Hε/Hβ ≈ 0.159, Hδ/Hβ ≈ 0.259, Hγ/Hβ ≈ 0.468), are marked by the horizontal blue lines. Under realistic assumptions, these ratios are fairly insensitive to the precise physical conditions of the ISM. For example, varying the temperature between $T_e = 0.5 \times 10^4$ K and $T_e = 3 \times 10^4$ K and the density between $N_e = 10$ cm⁻³ and $N_e = 10,000$ cm⁻³ would make these ratios change by at most ~ 3.5 per cent.

We note that while a poorly corrected stellar absorption could potentially bias the Balmer ratios, the continuum is nearly undetected in the spectra of these three galaxies. Moreover these are extremely young galaxies that are not expected to have significant Balmer absorption based on the level of continuum detection. While Balmer-line absorption is expected to be weaker for Hβ than for Hγ and Hδ at least, even assuming an equivalent width of 4 Å for Hβ absorption and 0 Å for the other Balmer lines would lead to a correction of less than 5 per cent.
The Balmer ratios shown in Fig. 2 are consistent with the theo-
retical value, or slightly lower. The exceptions are Hζ for galaxies
ID10612 and ID6355, however, Hζ is very faint and possibly con-
taminated by HeI. The other exception is Hδ, which is significantly
lower than the theoretical value for galaxies ID4590 and ID10612,
while is well in agreement with the theoretical value for ID6355;
however, the former two galaxies are those with the lowest S/N and,
in addition, some background residual problems seem to be present
at these wavelengths for these galaxies. We note that, as already re-
ported by Schäerer et al. (2022), the use of the level 3 products from
the MAST archive would result in Balmer ratios much higher than the
theoretical values, confirming the presence of residual issues in the
data processing and/or flux calibration suggested both in the JWST
documentation<fnref n="8" /> and by the significant levels of spectrophotometric
correction we had to introduce (see Sect. 2.1.2).

We simultaneously fit the ratios of Hγ/Hβ and Hε/Hβ to infer the dust attenuation by assuming the SMC curve from Gordon et al. (2003) with R<sub>V</sub> = 2.505; as mentioned, Hδ/Hβ appears strongly underestimated in both ID4590 and ID10612 compared to the other two Balmer decrements at fixed attenuation, possibly due to poor continuum subtraction, hence we consider this ratio only for ID6355. The inferred values of nebular A<sub>V</sub> for the three sources are reported in Table 2. The use of the locally-derived Gordon et al. (2003) SMC curve may be arguable for such distant galaxies, however it as been shown to be adequate for z ~ 2 low-metallicity galaxies (Shivaei et al. 2020). Moreover, since the [O III]λ4363 emission line is close in wavelength to the Hγ line and [O III]λ5007 is close to Hβ, and considering that the attenuation is mostly driven by Hγ/Hβ, the choice of the attenuation curve has a minor effect on the dust-corrected value of the [O III]4363/[O III]5007 line ratio (and hence on the inferred temperature and metallicity).

We compare the nebular attenuation to the attenuation of the con-
tinuum derived from the BEAGLE SED fitting, once scaling stellar
to nebular attenuation by a factor of 2.11⁹, following Shivaei et al.
(2020). We find agreement within the uncertainties for ID4590 and
ID10612, whereas a significantly higher attenuation from the SED
fitting in ID6355.

3.2 Excitation diagnostics

At the redshifts of our targets Hα and adjacent nebular lines are
shifted outside the NIRSpec band, so they cannot be used to con-
struct the classical BPT diagnostic diagrams (Baldwin et al. 1981).
However, interesting information about the excitation properties
of these sources can be inferred from the [O III]λ5007/Hβ versus
[O III]λ3727, 29/Hβ diagram. This diagram is shown in the left-hand
panel of Fig.3, where the shaded region marks the distribution of
local galaxies in the MPA-JHU catalog from the Sloan Digital Sky
Survey (Tremonti et al. 2004; Brinchmann et al. 2004) and other
surveys at z ~ 1 – 3 (Troncoso et al. 2014a; Onodera et al. 2016;
Sanders et al. 2021; Hayden-Pawson et al. 2022). The dashed line
indicates the local dividing line between star-forming galaxies and
active galactic nuclei (AGN) identified by Lamareille (2010).

Intermediate-redshift star-forming galaxies are characterized by a
large scatter, but they also tend to be shifted towards higher values of
[O III]/Hβ. This is an effect similar to that seen in the classical BPT
diagram (see the discussion in Maiolino & Mannucci 2019), which
has been attributed to a combination of higher ionization parameter

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Flux ratios between different Balmer lines and the Hβ as observed in all galaxies' spectra. The theoretical values expected from atomic physics (assuming the case B recombination at T= 1.5 × 10⁴K and N_e = 300 cm⁻3) are marked by the blue lines. All observed ratios are consistent with little dust attenuation. Grey shaded areas mark regions with unreliable ratios: non-physical Hζ/Hβ ratio might be driven by Hζ flux being contaminated by He I emission, whereas residuals in the background subtraction may affect the Hδ emission in ID4590 and ID10612.</figcaption></figure>

and harder radiation field associated with α-enhanced stellar popula-
tions (e.g. Strom et al. 2018; Topping et al. 2020, but see also Curti
et al. 2022). Therefore, finding that the two galaxies at z ~ 7.6 are loc-
ated above Lamareille (2010)'s line does not necessarily imply that
these are AGN, but that they follow the same trend as intermediate-
redshift galaxies.

Galaxy ID4590 at *z* ~ 8.5 is located in a region of the diagram poorly populated by both local and intermediate-redshift galaxies, with log([O III]/Hβ) ~ 0.5 and log([O II]/Hβ) ~ -0.75. According to the photoionisation models presented in Nakajima & Maiolino (2022), and shown in Fig. 4, this region is populated by young, metal poor galaxies. However, we also note that the presence of an AGN with low metallicity (Z ~ 10⁻³) can not be completely ruled out by these models.

The right-hand panel of Fig. 3 shows instead the excitation diagram
O32=[O III]λ5007/[O II]λ3727, 29 versus R23 = ([O III]λ4959, 5007 +
[O II]λ3727, 29)/Hβ, where the same local and intermediate-redshift
galaxies are shown as in the left-hand panel. This diagram can be
considered as a proxy of ionization parameter (primarily traced by
O32, Díaz et al. (2000)) versus metallicity (primary traced by R23,
e.g. Nagao et al. (2006)). However, high O32 values are sometimes
also considered as indicative of density-bounded clouds, possibly
associated with high escape fractions (f<sub>esc</sub>) of Lyman-continuum
(LyC) photons (Nakajima & Ouchi 2014; Nakajima et al. 2020; Bar-
row et al. 2020). All three galaxies at z>7 considered here have very
high values of O32. Especially, ID10612 at z=7.6 and ID4590 at
z=8.5 have extremely high values of O32 that are rarely seen even in
galaxies at z~2-3, if not in extreme [O III]λ5007 emitters with very
high equivalent widths (Tang et al. 2019). The fact that these galax-
ies have such high values of O32 can potentially indicate very high
ionisation parameters and/or large LyC escape fractions (Chisholm
et al. 2022). Based on Lyα emitters (LAEs) probed by different low-
and intermediate-redshift galaxy surveys with similar line ratios and
metallicities to our JWST sample we can provide (following e.g., Na-
kajima et al. 2020; Izotov et al. 2018b) a rough estimate of the f<sub>esc</sub> in
our galaxies. In particular, the O32 value for ID6355 (log O32 = 0.96)
is consistent with LyC leakers with f<sub>esc</sub> ≲ 0.1, whereas ID4590 and

<footnote n="8">⁸ see https://jwst-docs.stsci.edu/jwst-calibration-pipeline-caveats</footnote>

<footnote n="9">⁹ we note the exact determination of this factor is quite uncertain</footnote>
ID10612 show $O_{32}$ values (log $O_{32}$ = 1.18; 1.36, respectively) comparable with $f_{esc}$ from ~0.1 up to ~0.5. We stress however that the relationship between $O_{32}$ and $f_{esc}$ is quite scattered, and that a large $O_{32}$ may be a necessary but not sufficient condition for the presence of relevant LyC emission (Izotov et al. 2017).

### 3.3 Oxygen abundance determination

We exploit the highly significant detection of the [O III]λ4363 auroral line in JWST/NIRSpec spectra of the three galaxies studied here to measure the oxygen abundance (a proxy of the total gas-phase metallicity) by means of the electron-temperature ($T_e$) method. More specifically, [O III]λ4363 is detected at $8.6\sigma$, $7.3\sigma$, and $6.4\sigma$ in the spectra of ID4590, ID10612, and ID6355, respectively.

For the purposes of temperature and abundance computation, we model each galaxy as a single HII region constituted by two ionisation zones: the high-ionization one, traced by the $O^{++}$ ion; and the low-ionization one, traced by the $O^{+}$ ion. We have direct access to the temperature of the high-ionisation zone (t3) from the ratio between the [O III]λ4363 and [O III]λ45007 emission lines.

We implement the iterative procedure described in Nicholls et al. (2013) to infer t3 from this auroral-to-nebular diagnostic line ratio, adopting the atomic data from Palay et al. (2012), and assuming an electron density $N_e = 300 \text{ cm}^{-3}$ (a value typically measured in the ISM of $z \sim 2-3$ galaxies, e.g. Sanders et al. 2016); different values of $N_e$ would not impact significantly our results, as the [O III]λ4363/[O III]λ45007 ratio depends only mildly on gas density. The derived t3 temperatures for our galaxies are reported in Table 2. The temperatures derived for ID10612 and ID6355 (~17600 K and ~12,200 K) are consistent with the typical values of low- and intermediate-metallicity HII regions, whereas the ID4590 source at $z \sim 8.5$ is characterised by a much higher temperature (~28,600 K), possibly driven by low metallicity and high ionisation parameter, and similar to what observed in local extremely metal-poor galaxies (e.g. Izotov et al. 2018a, 2021a).

Since we do not have access to temperature diagnostics for the low-ionisation species (i.e., no auroral line emission is detected from $N^+$ or $S^+$ species in the spectra), here we adopt the temperature-temperature relation from Pilyugin et al. (2009) to infer the temperature of the low-ionisation zone (t2): $t_2 = 2640 + 0.835 \times t_3$. While the applicability of such relation at $z \sim 8$ is questionable, adopting a different locally calibrated t2-t3 relation does not significantly impact the total oxygen abundance determination since in these low-metallicity environments the $O^{++}$ abundance is largely dominant over $O^+$ (Andrews & Martini 2013; Curti et al. 2017). For comparison, adopting the t2-t3 relation from Garnett (1990) impacts the final metallicity calculation by less than 1 per cent.

Finally, we compute the ionic abundances of $O^+$ and $O^{++}$ using the *getIonAbundance* routine of PYNEB (v 1.1.10, Luridiana et al. 2012, 2015, adopting the atomic data from Palay et al. 2012), which provides the abundance of ionic species given the electron temperature (of the corresponding ionisation zone), the gas density, and the flux of a strong emission line (e.g., [O III]λ5007 for $O^{++}$, [O III]λ3727, 29 for $O^+$) relative to Hβ, and assuming the same atomic data from Palay et al. (2012). Here, we assume that the total oxygen abundance is the sum of the abundances of the singly- and doubly-ionised oxygen species, i.e., $O/H = O^+/H + O^{++}/H$, neglecting contributions from higher ionisation states.

The derived electron temperatures and oxygen abundances for our galaxy sample are reported in Table 2; uncertainties on both quantities are computed through Monte Carlo simulations by randomly perturbing (assuming a Gaussian noise distribution) one hundred

times all measured line fluxes by their measurement errors<fnref n="10" /> (including the uncertainty from reddening correction), randomly varying the density between 100 and $1000 \text{ cm}^{-3}$, and repeating the full temperature and abundance calculation for each iteration, taking the standard deviation of the distribution of temperatures and abundances obtained as the representative uncertainties on the measured values. Nonetheless, we note that these values probably underestimate the true uncertainties, because the systemic contributions are not considered. In particular, temperature stratification in HII regions can significantly halter the derived metallicity (Stasińska 2002; Kewley et al. 2019).

As briefly mentioned above, the high t3 measured in ID4590 places this object at the limit of what has ever been observed in the local Universe for star-forming galaxies. Such high temperature might also reveal the presence of a much harder ionising continuum associated with AGN activity (Dors et al. 2021; Riffel et al. 2021), which is indeed not ruled out also by photoionisation models, as already shown by comparing the observed line ratios with the predictions of the Nakajima & Maiolino (2022) models in Fig. 4. Although we favour the scenario of star-forming galaxy with low metallicity and high ionisation parameter, for comparison we verified that the analytic approach presented in Dors et al. (2021, equations 9 and 10), tuned for deriving $T_e$-based metallicities in AGNs, provides a comparable estimate of the total oxygen abundance, finding temperatures and metallicities consistent within the uncertainties with our fiducial $T_e$-method. We simply note that, in case of the presence of an AGN, the metallicity of ID4590 might be underestimated by not accounting for the abundance of the $O^{+++}$ species, and that the stellar mass and SFR determination would require further care in decomposing the stellar from the AGN contribution to the photometry. However, we also note that obscuration from the dusty tori in type II AGNs could weaken the [O III]λ4363 emission (e.g Nagao et al. 2001).

It is also interesting to note that the galaxy with the lowest metallicity, ID4590, is also the one with the highest dust attenuation. While one should not over-interpret this result, as the uncertainty on attenuation is large, this is something that is observed in some local dwarf, low metallicity galaxies (e.g. SBS 0335-052), which can have large dust masses and high dust attenuation. At a given gas mass, compactness and geometrical effect can strongly contribute to the higher dust attenuation; moreover, high gas density can foster the rapid growth of dust in the ISM (Schneider et al. 2016).

Finally, we note that, in a recent paper, Schaerer et al. (2022) analyses the same three galaxies, obtaining different results in terms of gas temperature and metallicity. They notice non-physical Balmer decrements in their spectra, and apply a power-law, wavelength dependent correction. They also infer un-physically high temperature for ID4590. As discussed in the previous sections, we obtain reasonable Balmer decrements, without introducing any *a posteriori* adjustment, and we can reliably measure the temperature in ID4590, which results consistent with the ISM properties of extremely metal-poor galaxies. We believe that these differences can be ascribed to our additional processing of the data, as discussed in detail in Sect. 2.1.2.

### 3.4 Comparison with strong-line metallicity calibrations

Metallicity diagnostics based on strong nebular emission-line ratios are widely adopted by the astronomical community to infer the chemical abundances of large samples of galaxies, for which auroral

<footnote n="10">¹⁰ we verified that errors from the *ppxf* fitting procedure are consistent with those inferred from bootstrapping</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** [O III]/Hβ versus [O II]/Hβ diagram (left-hand panel) and [O III]/[O II] versus R<sub>23</sub> diagram (right-hand panel) for local SDSS galaxies (grey contours marking the 70% and 90% of the distribution), galaxies at z ~ 1–3 from the literature, and the JWST sample at z ~ 8.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** [O III]/Hβ vs [O II]/Hβ diagram illustrating the photoionization models for low metallicity PopII, PopIII and AGNs from Nakajima & Maiolino 2022. Different symbols mark different population of objects (as explained in the legend), while different colours belong to different values of the gas-phase metallicity. The line ratios observed in our z ~ 8 galaxies are marked with red symbols.</figcaption></figure>

<table><caption>**Table 2.** Derived galaxy properties.</caption><thead><tr><th>Galaxy ID</th><th>4590</th><th>6355</th><th>10612</th></tr></thead><tbody><tr><td>Redshift</td><td>8.496</td><td>7.665</td><td>7.658</td></tr><tr><td>μ<sup>a</sup></td><td>3.74 ± 0.07</td><td>1.231 ± 0.002</td><td>1.339 ± 0.003</td></tr><tr><td>log(M<sub>⋆</sub>/M<sub>☉</sub>)<sup>b</sup></td><td>7.75 ± 0.07</td><td>8.72 ± 0.04</td><td>8.08 ± 0.04</td></tr><tr><td>log(SFR [M<sub>☉</sub>yr<sup>-1</sup>])<sup>b</sup></td><td>0.35 ± 0.07</td><td>1.47 ± 0.04</td><td>0.90 ± 0.04</td></tr><tr><td>A<sub>V</sub> [nebular]<sup>c</sup></td><td>0.68<sup>+0.34</sup><sub>-0.25</sub></td><td>0.0<sup>+0.1</sup><sub>0.0</sub></td><td>0.40<sup>+0.46</sup><sub>-0.27</sub></td></tr><tr><td>A<sub>V</sub> [stellar]<sup>d</sup></td><td>0.37 ± 0.04</td><td>0.50 ± 0.03</td><td>0.21 ± 0.03</td></tr><tr><td>T<sub>e</sub> ([O III][10<sup>-4</sup> K])</td><td>2.77 ± 0.42</td><td>1.20 ± 0.07</td><td>1.75 ± 0.16</td></tr><tr><td>12+log(O/H)</td><td>6.99 ± 0.11</td><td>8.24 ± 0.07</td><td>7.73 ± 0.12</td></tr></tbody></table>

<sup>a</sup> derived from the lens models presented in Mahler et al. 2022

<sup>b</sup> values are corrected for magnification (errors on μ are propagated)

<sup>c</sup> inferred from ratios of Balmer lines

<sup>d</sup> inferred from SED fitting

lines are often not detected. Several authors have calibrated such diagnostics by exploiting samples of local star-forming galaxies and HII regions with T<sub>e</sub>-based metallicity determinations (e.g. Pettini & Pagel 2004; Maiolino et al. 2008; Marino et al. 2013; Pilyugin & Grebel 2016; Curti et al. 2017; Bian et al. 2018; Sanders et al. 2021). However, an evolution in the ISM conditions of high-redshift galaxies compared to the local Universe might impact the intrinsic dependence of strong-line ratios on gas-phase metallicity, potentially hampering their use as abundance diagnostics at high redshift and thus biasing the assessment and interpretation of the chemical evolution history of galaxies.

Here, we exploit the T<sub>e</sub>-based abundances delivered by JWST/NIRSpec spectra to test the rest-frame optical, strong-line metallicity calibrations for the first time at z ~ 8. In Fig. 5, the T<sub>e</sub>-based log(O/H) measurements of the three galaxies in our sample are plotted against some of the most widely adopted strong-line metallicity diagnostics, specifically: R<sub>3</sub>=[O III]λ5007/Hβ, R<sub>2</sub>=[O III]λ3727, 29/Hβ, R<sub>23</sub>=[O III]λ4959, 5007+[O III]λ3727, 29/Hβ, O32=[O III]λ5007/[O II]λ3727, 29, and Ne3O2=[Ne III]λ3869/[O II]λ3727, 29.

Different calibration curves (in the form of polynomial relations between metallicity and strong-line ratios) are shown for these diagnostics in each panel, namely from Maiolino et al. (2008) Curti et al. (2020a), Sanders et al. (2021), Bian et al. (2018) and Nakajima et al. (2022). In brief, Curti et al. (2020a) complemented the Curti et al. (2017) calibrations based on T<sub>e</sub> metallicity measurements performed on stacked spectra of SDSS galaxies in the [O II]λ3727, 29/Hβ versus [O III]λ5007/Hβ diagram; here, we also include a new calibration of the Ne3O2 diagnostic based on the same methodology and combined sample of stacks and individual SDSS galaxies with [O III]λ4363 detections described in Curti et al. (2017), and which was not published before.<fnref n="11" /> Sanders et al. (2021) combined the measurements performed by Andrews & Martini (2013) in stellar mass (M<sub>⋆</sub>)-SFR stacks of SDSS galaxies with a sample of low-metallicity dwarfs from Berg et al. (2012), whereas Maiolino et al. (2008) combined a

<footnote n="11">¹¹ The calibration is presented in the usual form [log(R)=∑<sub>i=0</sub><sup>n</sup> a<sub>i</sub> x<sup>i</sup>], where x=12+log(O/H)-8.69 (solar metallicity from Allende Prieto et al. 2001), and the best-fit coefficients are: c<sub>0</sub> = -1.632; c<sub>1</sub> = -2.061; c<sub>2</sub> = -0.461.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** The relationships between $T_e$ metallicity and strong-line ratios for the JWST galaxy sample are compared with some widely adopted abundance calibrations, namely from Maiolino et al. 2008, Bian et al. 2018, Curti et al. 2017, 2020a, Sanders et al. 2021, and Nakajima et al. 2022. All these calibration relations are built from samples of local star-forming galaxies with $T_e$-based measurements, but are based on slightly different assumptions (see the text for more details). Solid lines represent the relationships in the metallicity calibration range as provided by the authors, whereas dashed lines mark extrapolations of the polynomial fit outside that range. Shaded areas encompass the region between the calibrations based on the low- and high-EW($H\beta$) sub-samples of metal-poor galaxies as presented in Nakajima et al. 2022.</figcaption></figure>

sample of galaxies with auroral line detections at low metallicity from Nagao et al. (2006) with predictions from photoionisation models in the high-metallicity regime. Nakajima et al. (2022) instead extended the calibrations based on the Curti et al. (2017) SDSS stacks to the low-metallicity regime leveraging $T_e$-measurements in extremely metal-poor galaxies from the EMPRESS survey (Kojima et al. 2020); in particular, in Fig. 5 we show the calibrations best-fit obtained by Nakajima et al. (2022) from the full galaxy sample with the solid line, whereas the shaded areas encompass the region between the calibrations based on the sub-samples of metal-poor galaxies characterised by high- and low-EW($H\beta$) (i.e., $EW(H\beta) > 200$ Å and $< 100$ Å, respectively). Finally, Bian et al. (2018) built a set of calibrations on a sample of local analogs of high-redshift galaxies (for which auroral lines are detected in stacked spectra), selected to match the location of observed $z \sim 2$ sources in the [N II]-BPTdiagram (Shapley et al. 2015; Strom et al. 2017); such calibrations, though still based on local galaxies, are thus specifically built for being applied to high-redshift galaxy spectra.

The two $z \sim 7.6$ sources, ID10612 and ID6533, present a higher level of nebular excitation than lower redshift calibration samples, placing them above the upper plateau of R23, R3 and O32 curves, and below the expected R2 value at given metallicity, in the Curti et al. (2020a), Maiolino et al. (2008), Nakajima et al. (2022), and Sanders et al. (2021) calibrations, whereas the Bian et al. (2018) curves better cover that region of the line-ratios parameter space. The extremely metal-poor $z \sim 8.5$ ID4590 source, instead, is located in a region outside the validity range of most calibrations under study; nonetheless, it exhibits line-ratios properties resembling those

of extremely metal-poor galaxies in the local Universe, and especially Maiolino et al. (2008) and Nakajima et al. (2022) calibrations, as well as some of the extrapolations of Curti et al. 2020a and Sanders et al. 2021 curves (shown as dashed lines<fnref n="12" />), provide a reasonable match to most of the line ratios observed in this galaxy.

Some of the observed deviations are actually still within the scatter of the calibration relations (typically ranging between 0.1-0.2 dex in a given line ratio at fixed metallicity, depending on the diagnostic and the calibration sample). In Table 3, we compare the observed and predicted (given the measured metallicity) line ratios for each considered strong-line diagnostics (i.e., the vertical offset from the calibration relation), and report the significance of the deviation (in units of $\sigma$) taking into account both measurements uncertainties and the intrinsic dispersion of each individual calibration<fnref n="13" />. In many cases, the measurements clearly fail to simultaneously match all the relations, prompting for a self-consistent revision of the calibrations at such early epochs which could be applied to galaxy spectra at different redshifts and/or observed in different filters, hence probing different emission lines.

<footnote n="12">¹² It should be stressed that sometimes high-order polynomials express non-physical behaviours outside their fitting range</footnote>

<footnote n="13">¹³ Since Bian et al. 2018 do not provide an estimate of the dispersion of the calibrations, we assume a typical dispersion of 0.15 dex for all diagnostics. We further note that in some cases, e.g., the $Ne_3O_2$ calibration from Nakajima et al. 2022, a very large intrinsic dispersion of the calibration may artificially lower the significance of the computed deviation, which is therefore not a good representation of the calibration accuracy.</footnote>
<table><caption>**Table 3.** Significance of the deviation (in units of $\sigma$) of the observed line ratios in the JWST z~8 galaxy sample from the predictions of each of the strong-line calibrations presented in Fig. 5. Both the measurements uncertainties on the line ratios and the intrinsic dispersion of the calibrations ($\sigma_{cal}$, reported here as provided by the various authors) are taken into account.</caption><thead><tr><th>Galaxy ID</th><th>R2</th><th>R3</th><th>O32</th><th>R23</th><th>Ne3O2</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center; font-weight: bold;">Maiolino et al. 2008</td></tr><tr><td>σ<sub>cal</sub></td><td>0.10</td><td>0.10</td><td>0.20</td><td>0.06</td><td>0.06</td></tr><tr><td>4590</td><td>1.68</td><td>3.03</td><td>0.08</td><td>1.77</td><td>0.39</td></tr><tr><td>6355</td><td>2.31</td><td>1.88</td><td>4.39</td><td>0.96</td><td>3.16</td></tr><tr><td>10612</td><td>2.20</td><td>0.86</td><td>2.97</td><td>0.41</td><td>3.18</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Curti et al. 2017, 2020a</td></tr><tr><td>σ<sub>cal</sub></td><td>0.11</td><td>0.09</td><td>0.15</td><td>0.06</td><td>0.12</td></tr><tr><td>4590</td><td>14.94</td><td>3.22</td><td>2.54</td><td>0.59</td><td>1.97</td></tr><tr><td>6355</td><td>2.99</td><td>2.62</td><td>3.80</td><td>1.80</td><td>3.82</td></tr><tr><td>10612</td><td>1.89</td><td>1.80</td><td>2.11</td><td>0.64</td><td>2.04</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Bian et al. 2018</td></tr><tr><td>σ<sub>cal</sub><sup>1</sup></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>4590</td><td>-</td><td>0.35</td><td>5.23</td><td>1.34</td><td>4.7</td></tr><tr><td>6355</td><td>-</td><td>1.50</td><td>3.71</td><td>0.95</td><td>2.52</td></tr><tr><td>10612</td><td>-</td><td>0.01</td><td>0.57</td><td>0.004</td><td>0.91</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Sanders et al. 2021</td></tr><tr><td>σ<sub>cal</sub></td><td>0.13</td><td>0.10</td><td>0.19</td><td>0.08</td><td>0.20</td></tr><tr><td>4590</td><td>3.10</td><td>2.31</td><td>1.98</td><td>1.37</td><td>0.92</td></tr><tr><td>6355</td><td>3.65</td><td>2.86</td><td>4.05</td><td>1.77</td><td>2.64</td></tr><tr><td>10612</td><td>3.29</td><td>3.86</td><td>4.32</td><td>1.49</td><td>4.28</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Nakajima et al. 2022</td></tr><tr><td>σ<sub>cal</sub></td><td>0.27</td><td>0.16</td><td>0.39</td><td>0.10</td><td>0.42</td></tr><tr><td>4590</td><td>0.95</td><td>2.54</td><td>1.62</td><td>1.18</td><td>0.95</td></tr><tr><td>6355</td><td>1.26</td><td>1.75</td><td>1.72</td><td>1.38</td><td>0.18</td></tr><tr><td>10612</td><td>1.29</td><td>0.74</td><td>1.32</td><td>0.53</td><td>1.58</td></tr></tbody></table>

# 4 COSMIC EVOLUTION OF THE METALLICITY SCALING RELATIONS

## 4.1 The Mass-metallicity relation

In this section we investigate the evolution in the metallicity scaling relations of galaxies as probed by the $T_e$-based O/H measurements enabled by JWST/NIRSpec at $z \sim 8$. In Fig. 6, we plot the three galaxies analysed in the current paper on the mass-metallicity plane, along with the local relation inferred for the SDSS sample by Curti et al. 2020a from $T_e$-based calibrations (contours in grey, best-fit MZR in black). Compared to the MZR in the local Universe probed by SDSS galaxies, galaxies at $z \sim 8$ appear metal deficient at fixed stellar mass, in agreement with the expected redshift evolution of the MZR (see Maiolino & Mannucci 2019, and references therein). More specifically, ID6355 is offset by 0.18 dex from the $z \sim 0$ MZR, whereas ID10612 is offset by 0.52 dex. We also compare these observations at $z \sim 8$ with the MZR at $z \sim 2$ and $z \sim 3.3$, here parametrised by a linear regression fit performed on stacked spectra from the MOSDEF survey as presented in Sanders et al. (2021). These

authors adopt the Bian et al. (2018) $T_e$-empirical calibrations based on local analogues of high-$z$ galaxies in order to reduce the potential biases affecting the use of local strong-line metallicity diagnostics at high redshift. Extrapolating the Sanders et al. (2021) relation to the low-mass range, we find that the metallicity of ID6355 is broadly consistent with the $z \sim 2.2$ MZR, whereas ID10612 is offset by 0.12 dex from the $z \sim 3.3$ MZR extrapolation. Naively, this could be interpreted as a signature of little or no evolution of the low-mass end of the MZR (either in slope or normalisation) between $z \sim 3.3$ and $z \sim 7.6$, possibly implying a mild cosmic evolution in both the gas fraction (and consequent metal dilution) and the metal-loading factor of star-formation-driven outflows over $\sim 1.2$ Gyr (Sanders et al. 2021).

In contrast, we find that galaxy ID4560 is more than 1 dex offset from the $z \sim 0$ MZR, and deviates strongly ($\sim 0.8$ dex) also from the extrapolation of the MZR at $z \sim 3.3$. Rather than very rapid redshift evolution of the MZR, this could indicate that ID4560 is in a very early evolutionary stage, in the process of rapidly building up its metals and approaching the MZR within a timeframe of just a few tens Myr (e.g. Fig. 11 in Maiolino et al. 2008). Alternatively, this finding could indicate that the slope of the MZR becomes steeper at very low masses (not probed by the $z \sim 2.2-3.3$ surveys).

## 4.2 The fundamental metallicity relation

In this section we explore whether these galaxies at $z \sim 8$ are consistent with the Fundamental Metallicity Relation (FMR), which (as discussed in the introduction) describes the correlation between metallicity, stellar mass, and SFR, observational consequence of the interplay between gas accretion, star formation, and outflows which govern the secular evolution of galaxies, and which was found not to evolve (or marginally evolve), out to $z \sim 3$. Fig. 7 shows the deviation of the measured galaxy metallicities from the predictions of the FMR, as parameterized by equation 5 of Curti et al. (2020a). Different samples with available rest-frame optical spectroscopy compiled from the literature at various redshifts are included in the attempt to trace the cosmic evolution of this scaling relation across almost the entire history of the Universe. Specifically, these galaxies are compiled from Cresci et al. (2012, zCOSMOS at $z \sim 0.3-0.6$), Hayden-Pawson et al. (2022, KLEVER at $z \sim 2.2$), Sanders et al. (2021, MOSDEF at $z \sim 2.3-3.3$), Onodera et al. (2016, $z \sim 3.3$), and Troncoso et al. (2014b, AMAZE at $z \sim 3.3$). The metallicities of galaxies at $z > 1$ have been consistently re-computed adopting the Bian et al. (2018) $T_e$-based strong-line calibrations (to account for potential evolution in the diagnostics), whereas we adopt the Curti et al. (2020a) calibrations for SDSS and zCOSMOS galaxies. For comparison, adopting instead the Curti et al. (2020a) calibrations for all galaxy samples produces, on average, $a \sim 0.1$ dex lower metallicity in $z \sim 2-3$ galaxies. The various symbols in Fig. 7 mark, for each considered literature dataset, the average offset in O/H (and standard deviation as errorbars) from the FMR predictions for individual galaxies within the sample. Red symbols mark instead individual $T_e$-based measurements for the JWST galaxy sample at $z \sim 8$; the errorbars on such datapoints are obtained by co-adding in quadrature the uncertainty on the log(O/H) measurement from the $T_e$ method and the uncertainty on the metallicity predicted by the FMR, evaluated from the standard deviation of the distribution of metallicities obtained by varying one hundred times the input $M_\star$ and SFR of galaxies within their errors.

The two galaxies at $z \sim 7.6$ show offsets of about $\sim 0.2$ dex (significant at $\sim 2\sigma$) from the FMR, though with different signs, suggesting different evolutionary stages. However, we also note that considering additional sources of uncertainty, like the errors associated to the

<footnote n="1">¹ For Bian et al. (2018) we assume a dispersion of 0.15 dex on each calibration.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Mass-Metallicity relation (MZR) for the JWST sample. The MZR for SDSS galaxies from Curti et al. 2020a and its best-fit are shown by the grey contours and black curve, respectively. The best-fit to stacked spectra of MOSDEF galaxies from Sanders et al. 2021 (based on the Bian et al. 2018 calibrations tuned for high-z) are shown as representative for the MZR at z ~ 2.3 and z ~ 3.3. Moreover, the MZR at z ~ 8 predicted by TNG (Nelson et al. 2019), EAGLE (Schaye et al. 2015), FIRE (Ma et al. 2016) (extrapolated to z=8), and SERRA (Pallottini et al. 2022) simulations are shown as dashed purple, orange, brown (light brown for re-scaled relation to match both local and z=2–3 observations) and blue lines, respectively, with the shaded areas marking the region between the 16th and 84th percentiles of the predicted relation.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Deviations in the observed log(O/H) from the predictions of the local Fundamental Metallicity Relation (FMR, from Curti et al. 2020a, equation 5), plotted as a function of redshift (look back time is reported on the top axis). High-redshift galaxies with $T_e$-based measurements from JWST/NIRSpec are shown as red symbols. The average (and standard deviation) offset computed for different datasets complied from the literature at different redshifts are also shown for comparison. All metallicities from the literature samples are self-consistently re-computed adopting the $T_e$-based Curti et al. 2020a calibrations at $z< 1$ and the Bian et al. 2018 at $z> 1$. At $z \sim 8$, galaxies appear offset from the predictions of the local FMR (though with different levels of significance), suggesting that they are far from the smooth equilibrium between chemical enrichment and gas flows that sets the average scaling relations in local galaxies (and up to $z \sim 3$).</figcaption></figure>

• We robustly detect the [O III]λ4363 auroral line in all three galaxies, with ratios relative to the strong [O III]λ4959,5007 doublet that are large relative to local galaxies (0.01–0.04).

• The inferred gas temperatures range from 1.2 10^4 K to 2.8 10^4 K, with the $T_e$-based metallicities ranging from extremely metal poor (12+log(O/H)~7) to about one third solar.

• None of the local strong-line metallicity calibrations seem to provide a good prediction of the observed metallicities at $z \sim 8$ simultaneously across all diagnostics and for all galaxies (Fig. 5). Different calibrations seem to perform better for some galaxies or for some diagnostic ratios, suggesting that a systemic re-calibration of the strong-line diagnostics is needed for these early epochs.

• The two most “massive” galaxies of the sample at $z \sim 7.6$ (log($M_*/M_\odot$)=8.1-8.7) are scattered around the mass-metallicity relation at $z \sim 2-3$ (Fig. 6), potentially suggestive of little evolution of the MZR from $z \sim 2-3$ to $z \sim 8$. On the contrary, the least massive and most metal-poor galaxy (log($M_*/M_\odot$)=7.7, 12+log(O/H)=6.99) deviates significantly from the relation at lower redshifts, suggesting that this galaxy is in a phase of rapid evolution or that the MZR has a much steeper slope at such low masses. Different cosmological simulations can reproduce some of the galaxies on the MZR, but none of them seem capable of simultaneously matching all of them.

• The three galaxies are widely scattered relative to the Fundamental Metallicity Relation, with two of them marginally consistent and one (ID4590 at $z \sim 8.5$) strongly deviating (Fig. 7). This finding suggests that galaxies at such early epochs are undergoing rapidly evolving mechanisms, and not yet settled on the local scaling relations that are instead associated with smooth, secular evolutionary processes. In any case, much larger samples are required in order to draw strong conclusions on the evolution of such scaling relations.

In conclusion, these first datasets already highlight the giant leap in the redshift parameter space brought to the field by the advent of rest-frame optical spectroscopy enabled by JWST, which is about to open a window on new modes of galaxy formation and evolution, helping to discover objects with properties radically different than all the other known galaxies. Several forthcoming observational programmes will target the ‘redshift desert’ between $z \sim 3–8$ (and beyond), allowing us to obtain a complete and coherent picture of the evolution of the chemical properties of galaxies across the entire cosmic history.

Nonetheless, we also stress that assessing the applicability of the strong-line metallicity diagnostics at high redshift remains an open problem, which requires both more a carefully selected sample of galaxies and deeper observations to detect auroral lines over a wider range of metallicities. While the observations discussed in the present paper provide a preview of the possibilities opened by JWST in this sense, a variety of programmes scheduled in Cycle 1 primarily designed to address this science topic (e.g. Curti et al. 2021; Strom et al. 2021; Shapley et al. 2021) will soon provide a more robust re-calibration of the strong-line metallicity diagnostics for the high-redshift Universe.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for his/her comments which contributed to improving this paper.

M.C., R.M., F.D'E., L.S., J.S., W.B., J.W., T.J.L., C.H-P. and J.S.B. acknowledge support by the Science and Technology Facilities Council (STFC) and ERC Advanced Grant 695671 "QUENCH". RM also acknowledges funding from a research professorship from the Royal Society. S.C. acknowledges support from the ERC Advanced Grant
INTERSTELLAR H2020/740120. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. SA acknowledges funding from grant PID2021-127718NB-I00 by the Spanish Ministry of Science and Innovation/State Agency of Research (MCIN/AEI). J.M.P and J.W. gratefully acknowledge support from the MERAC Foundation. AJB acknowledges 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). RS acknowledges support from an STFC Ernest Rutherford Fellowship (ST/S004831/1). ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1). We are grateful to Pierre Ferruit, Peter Jakobsen, and Nora Lützgendorf for sharing their expertise on NIRSpec and the processing of its unique data.

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. Some of the data presented in this paper were obtained from the Multimission Archive at the Space Telescope Science Institute (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts.

Finally, M.C. is heartily grateful to Vittoria for showing him a completely new and wonderful, yet unknown and mysterious, part of the Universe.

## DATA AVAILABILITY

The raw (level 2) data underlying this article are publicly available from the MAST archive. The processed and calibrated 1D spectra analysed in this work are publicly available at https://doi.org/10.5281/zenodo.6940561. Other data are available upon reasonable request to the corresponding author.

## References

* Allende Prieto C., Lambert D. L., Asplund M., 2001, ApJ, 556, L63

* Andrews B. H., Martini P., 2013, ApJ, 765, 140

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

* Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5

* Barrow K. S. S., Robertson B. E., Ellis R. S., Nakajima K., Saxena A., Stark D. P., Tang M., 2020, ApJ, 902, L39

* Berg D. A., et al., 2012, ApJ, 754, 98

* Bian F., Kewley L. J., Dopita M. A., 2018, ApJ, 859, 175

* Bohlin R. C., Gordon K. D., Tremblay P. E., 2014, PASP, 126, 711

* Bohlin R. C., Hubeny I., Rauch T., 2020, AJ, 160, 21

* Bouché N., et al., 2010, ApJ, 718, 1001

* Brinchmann J., Charlot S., White S. D. M., Tremonti C., Kauffmann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 1151

* Caminha G. B., Suyu S. H., Mercurio A., Brammer G., Bergamini P., Vanzella E., Acebron A., 2022, arXiv e-prints, p. arXiv:2207.07567

* Cappellari M., 2017, MNRAS, 466, 798

* Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379

* Carnall A. C., et al., 2022, arXiv e-prints, p. arXiv:2207.08778

* Chabrier G., 2003, PASP, 115, 763

* Chevallard J., Charlot S., 2016, MNRAS, 462, 1415

* Chisholm J., et al., 2022, arXiv e-prints, p. arXiv:2207.05771

* Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102

* Christensen L., et al., 2012, MNRAS, 427, 1973

* Conroy C., Naidu R. P., Zaritsky D., Bonaca A., Cargile P., Johnson B. D., Caldwell N., 2019, ApJ, 887, 237

* Crain R. A., et al., 2015, MNRAS, 450, 1937

* Cresci G., Mannucci F., Sommariva V., Maiolino R., Marconi A., Brusa M., 2012, MNRAS, 421, 262

* Cresci G., Mannucci F., Curti M., 2019, A&A, 627, A42

* Curti M., Cresci G., Mannucci F., Marconi A., Maiolino R., Esposito S., 2017, MNRAS, 465, 1384

* Curti M., Mannucci F., Cresci G., Maiolino R., 2020a, Monthly Notices of the Royal Astronomical Society, 491, 944

* Curti M., et al., 2020b, Monthly Notices of the Royal Astronomical Society, 492, 821

* Curti M., et al., 2021, Opening the era of direct metallicity measurements in high redshift galaxies, JWST Proposal. Cycle 1, ID. #1879

* Curti M., et al., 2022, MNRAS, 512, 4136

* Davé R., Rafieferantsoa M. H., Thompson R. J., Hopkins P. F., 2017, MNRAS, 467, 115

* Díaz A. I., Castellanos M., Terlevich E., Luisa García-Vargas M., 2000, MNRAS, 318, 462

* Dors O. L., Contini M., Riffel R. A., Pérez-Montero E., Krabbe A. C., Cardaci M. V., Hägele G. F., 2021, MNRAS, 501, 1370

* Ebeling H., Edge A. C., Henry J. P., 2001, ApJ, 553, 668

* Ebeling H., Barrett E., Donovan D., Ma C. J., Edge A. C., van Speybroeck L., 2007, ApJ, 661, L33

* Ebeling H., Edge A. C., Mantz A., Barrett E., Henry J. P., Ma C. J., van Speybroeck L., 2010, MNRAS, 407, 83

* Ebeling H., et al., 2013, MNRAS, 432, 62

* Ellison S. L., Patton D. R., Simard L., McConnachie A. W., 2008, ApJ, 672, L107

* Erb D. K., Pettini M., Steidel C. C., Strom A. L., Rudie G. C., Trainor R. F., Shapley A. E., Reddy N. A., 2016, ApJ, 830, 52

* Ferruit P., et al., 2022, A&A, 661, A81

* Garnett D. R., 1990, ApJ, 363, 142

* Gordon K. D., Clayton G. C., Misselt K. A., Landolt A. U., Wolff M. J., 2003, ApJ, 594, 279

* Gordon K. D., et al., 2022, AJ, 163, 267

* Hayden-Pawson C., et al., 2022, MNRAS, 512, 2867

* Izotov Y. I., Thuan T. X., Guseva N. G., 2017, MNRAS, 471, 548

* Izotov Y. I., Thuan T. X., Guseva N. G., Liss S. E., 2018a, MNRAS, 473, 1956

* Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Thuan T. X., Fricke Verhamme A., Orlitová I., 2018b, MNRAS, 478, 4851

* Izotov Y. I., Guseva N. G., Fricke K. J., Henkel C., 2019, arXiv e-prints

* Izotov Y. I., Thuan T. X., Guseva N. G., 2021a, MNRAS, 504, 3996

* Izotov Y. I., Thuan T. X., Guseva N. G., 2021b, MNRAS, 504, 3996

* Jakobsen P., et al., 2022, A&A, 661, A80

* Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021, ApJS, 254, 22

* Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, Annual Review of Astronomy and Astrophysics, 57, 511

* Kojima T., et al., 2020, ApJ, 898, 142

* Lamareille F., 2010, A&A, 509, A53

* Langan I., Ceverino D., Finlator K., 2020, MNRAS, 494, 1988

* Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013, ApJ, 772, 119

* Luridiana V., Morisset C., Shaw R. A., 2012, in IAU Symposium. pp 422–423, doi:10.1017/S1743921312011738

* Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42

* Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016, MNRAS, 456, 2140

* Madau P., Dickinson M., 2014, ARA&A, 52, 415

* Mahler G., et al., 2022, arXiv e-prints, p. arXiv:2207.07101

* Maiolino R., Mannucci F., 2019, A&ARv, 27, 3

* Maiolino R., et al., 2008, A&A, 488, 463

* Mann A. W., Ebeling H., 2012, MNRAS, 420, 2120

* Mannucci F., et al., 2009, MNRAS, 398, 1915

* Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010, MNRAS, 408, 2115

* Mannucci F., Salvaterra R., Campisi M. A., 2011, MNRAS, 414, 1263
Marinacci F., et al., 2018, MNRAS, 480, 5113

Marino R. A., et al., 2013, A&A, 559, A114

McAlpine S., et al., 2016, Astronomy and Computing, 15, 72

Nagao T., Murayama T., Taniguchi Y., 2001, ApJ, 549, 155

Nagao T., Maiolino R., Marconi A., 2006, A&A, 459, 85

Naiman J. P., et al., 2018, MNRAS, 477, 1206

Nakajima K., Maiolino R., 2022, MNRAS, 513, 5134

Nakajima K., Ouchi M., 2014, MNRAS, 442, 900

Nakajima K., Ellis R. S., Robertson B. E., Tang M., Stark D. P., 2020, ApJ, 889, 161

Nakajima K., et al., 2022, arXiv e-prints, p. arXiv:2206.02824

Nelson D., et al., 2018, MNRAS, 475, 624

Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2

Nicholls D. C., Dopita M. A., Sutherland R. S., Kewley L. J., Palay E., 2013, ApJS, 207, 21

Noeske K. G., et al., 2007, ApJ, 660, L47

Onodera M., et al., 2016, ApJ, 822, 42

Osterbrock D. E., Ferland G. J., 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, 2nd edn. University Science Books

Palay E., Nahar S. N., Pradhan A. K., Eissner W., 2012, MNRAS, 423, L35

Pallottini A., et al., 2022, MNRAS, 513, 5621

Pascale M., et al., 2022, arXiv e-prints, p. arXiv:2207.07102

Patrício V., Christensen L., Rhodin H., Cañamas R., Lara-López M. A., 2018, MNRAS, 481, 3520

Perrin M. D., Long J., Sivaramakrishnan A., Lajoie C.-P., Elliot E., Pueyo L., Albert L., 2015, WebbPSF: James Webb Space Telescope PSF Simulation Tool, Astrophysics Source Code Library, record ascl:1504.007 (ascl:1504.007)

Pettini M., Pagel B. E. J., 2004, MNRAS, 348, L59

Pillepich A., et al., 2018, MNRAS, 475, 648

Pilyugin L. S., Grebel E. K., 2016, MNRAS, 457, 3678

Pilyugin L. S., Mattsson L., Vílchez J. M., Cedrés B., 2009, MNRAS, 398, 485

Pontoppidan K., et al., 2022, arXiv e-prints, p. arXiv:2207.13067

Popesso P., et al., 2022, arXiv e-prints, p. arXiv:2203.10487

Repp A., Ebeling H., 2018, MNRAS, 479, 844

Repp A., Ebeling H., Richard J., 2016, MNRAS, 457, 1399

Riffel R. A., et al., 2021, MNRAS, 501, L54

Sanders R. L., et al., 2016, ApJ, 816, 23

Sanders R. L., et al., 2020, MNRAS, 491, 1427

Sanders R. L., et al., 2021, ApJ, 914, 19

Sandles L., Curtis-Lake E., Charlot S., Chevallard J., Maiolino R., 2022, MNRAS, 515, 2951

Schaerer D., Marques-Chaves R., Oesch P., Naidu R., Barrufet L., Izotov Y. I., Guseva N. G., Brammer G., 2022, arXiv e-prints, p. arXiv:2207.10034

Schaye J., et al., 2015, MNRAS, 446, 521

Schneider R., Hunt L., Valiante R., 2016, MNRAS, 457, 1842

Shapley A. E., et al., 2015, ApJ, 801, 88

Shapley A. E., et al., 2017, ApJ, 846, L30

Shapley A. E., et al., 2021, The AURORA Survey: First Direct Metallicity Calibrations at High Redshift, JWST Proposal. Cycle 1, ID. #1914

Shivaei I., et al., 2020, The Astrophysical Journal, 899, 117

Springel V., et al., 2018, MNRAS, 475, 676

Stasińska G., 2002, Cosmochemistry. The melting pot of the elements, pp 115–170

Steidel C. C., et al., 2014, ApJ, 795, 165

Strom A. L., Steidel C. C., Rudie G. C., Trainor R. F., Pettini M., Reddy N. A., 2017, ApJ, 836, 164

Strom A. L., Steidel C. C., Rudie G. C., Trainor R. F., Pettini M., 2018, ApJ, 868, 117

Strom A. L., Rudie G. C., Chen Y., Law D. R., Maseda M., Steidel C. C., Trainor R. F., 2021, CECILIA: A direct-method metallicity calibration for Cosmic Noon through the Epoch of Reionization, JWST Proposal. Cycle 1, ID. #2593

Tacchella S., et al., 2022a, ApJ, 927, 170

Tacchella S., et al., 2022b, ApJ, 927, 170

Tang M., Stark D. P., Chevallard J., Charlot S., 2019, MNRAS, 489, 2572

Topping M. W., Shapley A. E., Reddy N. A., Sanders R. L., Coil A. L., Kriek M., Mobasher B., Siana B., 2020, MNRAS, 495, 4430

Torrey P., et al., 2019, MNRAS, 484, 5587

Tremonti C. A., et al., 2004, ApJ, 613, 898

Troncoso P., et al., 2014a, A&A, 563, A58

Troncoso P., et al., 2014b, A&A, 563, A58

Witstok J., Smit R., Maiolino R., Curti M., Laporte N., Massey R., Richard J., Swinbank M., 2021, MNRAS, 508, 1686

Zahid H. J., Kewley L. J., Bresolin F., 2011, ApJ, 730, 137

Zahid H. J., Dima G. I., Kudritzki R.-P., Kewley L. J., Geller M. J., Hwang H. S., Silverman J. D., Kashino D., 2014, ApJ, 791, 130

## AFFILIATIONS

¹Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK

²Cavendish Laboratory - Astrophysics Group, University of Cambridge, 19 JJ Thompson Avenue, Cambridge, CB3 0HE, UK

³Scuola Normale Superiore, Università di Pisa, Piazza dei Cavalieri 7, I-56126 Pisa, Italy

⁴Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

⁵Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK

⁶Sorbonne Université, UPMC-CNRS, UMR7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

⁷National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

⁸INAF — Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125, Florence, Italy

⁹Dipartimento di Fisica e Astronomia, Università di Firenze, Via G. Sansone 1, I-50125, Sesto Fiorentino (Florence), Italy

¹⁰Centro de Astrobiología (CAB), CSIC-INTA, Cra. de Ajalvir Km. 4, 28850- Torrejón de Ardoz, Madrid, Spain

¹¹Cosmic Dawn Center, Niels Bohr Institute, University of Copenhagen, Rådmandsgade 62, 2200 Copenhagen N, DK

¹²Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

¹³Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield, AL10 9AB, UK

¹⁴AURA for the European Space Agency, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁵Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706 USA

¹⁶Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK

This paper has been typeset from a TEX/LATEX file prepared by the author.