Open Access
Issue
A&A
Volume 684, April 2024
Article Number A75
Number of page(s) 22
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202346698
Published online 10 April 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The James Webb Space Telescope (JWST) has already begun to revolutionise our view of the high-redshift Universe, opening a new window onto early galaxy formation. In particular, the high-sensitivity Near Infrared Spectrograph (NIRSpec, Jakobsen et al. 2022; Ferruit et al. 2022; Böker et al. 2023) now allows to spectroscopically characterise the interstellar medium (ISM) properties and constrain the nature of the ionising spectra in primordial galaxies. In this context, the study of the gas-phase metallicity1 via rest-frame optical spectroscopy, which has been limited for decades to z ≤ 3.5 due to the intrinsic limitations of ground-based observatories, is now possible out to z ∼ 10. Such observations are providing precious constraints on cosmological simulations and chemical evolution models, enabling scientists to chart the processes that shaped the formation of galaxy structures in the early Universe (Somerville & Davé 2015; Maiolino & Mannucci 2019). For instance, characterising the scaling relation between stellar mass and gas-phase metallicity in galaxies (the so-called mass–metallicity relation (MZR); Lequeux et al. 1979; Tremonti et al. 2004; Lee et al. 2006; Yates et al. 2020; Baker & Maiolino 2023) at high redshift is critical in order to constrain the processes regulating the growth of early galaxies, as this relation is shaped by the interplay between gas accretion, star formation, metal enrichment, and outflows driving the baryon cycle. Furthermore, as a consequence of the long-lasting interplay between these different processes, the metallicity of galaxies has been observed to also correlate with different galactic properties. In particular, an anti-correlation between metallicity and star-formation rate (SFR) at fixed stellar mass has been observed and characterised in detail on both global and local scales in what is usually referred to as the fundamental metallicity relation (FMR; Ellison et al. 2008; Mannucci et al. 2010; Yates et al. 2012; Salim et al. 2015; Telford et al. 2016; Curti et al. 2020b; Baker et al. 2023).

Several studies leveraged the new spectroscopic capabilities of JWST/NIRSpec to start characterising the properties of strong line emitters at z >  3, exploiting the data collected in the framework of Early Release Observations (Pontoppidan et al. 2022, ERO), Early Release Science (ERS), and Cycle 1 General Observer (GO) and Guaranteed Time Observations (GTO) programmes, such as GLASS (Proposal ID: 1324; Treu et al. 2022), Cosmic Evolution Early Release Science (CEERS, Proposal ID: 1345; Finkelstein et al. 2023), and JWST Advanced Deep Extragalactic Survey (JADES, Proposal ID: 1210; Eisenstein et al. 2023; Robertson et al. 2023; Curtis-Lake et al. 2023). Recent works have indeed investigated the ionisation properties of galaxies beyond z = 3, showing that these sources exhibit emission-line ratios consistent with hard ionising spectra and low metallicities (Mascia et al. 2023; Matthee et al. 2023; Sanders et al. 2023; Cameron et al. 2023b). These kinds of analyses have been pushed to some of the highest-redshift galaxies discovered, with Williams et al. (2023) reporting a z = 9.5 galaxy with a very high [O III]/[O II] ratio and relatively low metallicity, while Bunker et al. (2023b) recently reported the detection of multiple emission lines (including several high-ionisation metal lines) in the NIRSpec spectrum of the luminous z ∼ 10.61 galaxy GN-z11 (Oesch et al. 2016). The wealth of emission lines detected in these high-z galaxy spectra can therefore be exploited to assess the chemical enrichment of their ISM.

In the past, several studies have investigated the evolution in the metallicity scaling relations out to z ∼ 3, finding signatures of a clear trend of decreasing metallicity with redshift at fixed stellar mass (e.g. Shapley et al. 2005; Erb et al. 2006; Maiolino et al. 2008; Mannucci et al. 2009; Zahid et al. 2011, 2014; Henry et al. 2013a; Yabe et al. 2014; Wuyts et al. 2014; Guo et al. 2016; Sanders et al. 2021; Topping et al. 2021). However, when reported to the framework of the FMR, no evidence of redshift evolution has been observed, suggesting that, on average, galaxies evolved through smooth secular processes over the past ∼10 Gyr driven by the interplay of gas flows, star formation, and metal enrichment, which is reflected in this latter scaling relation (Mannucci et al. 2010; Belli et al. 2013; Nakajima & Ouchi 2014; Maier et al. 2014; Salim et al. 2015; Hirschauer et al. 2018; Sanders et al. 2021; Hayden-Pawson et al. 2022).

More recent studies attempted to push the investigation of the chemical properties of galaxies up to z ∼ 10, either from the detection of emission lines in individual sources, or from the analysis of composite spectra, finding initial evidence for a relatively mild evolution of the oxygen abundance with redshift at fixed stellar mass (e.g. Curti et al. 2023; Schaerer et al. 2022; Arellano-Córdova et al. 2022; Taylor et al. 2022; Trump et al. 2023; Rhoads et al. 2023; Langeroodi et al. 2023; Matthee et al. 2023; Heintz et al. 2023; Nakajima et al. 2023; Shapley et al. 2023a). Despite these efforts, the majority of studies focused on relatively high-stellar-mass systems, whereas probing the evolution of the metallicity scaling relations at low stellar masses, a regime strongly sensitive to feedback processes due to the shallower gravitational potential, remained challenging (Wuyts et al. 2012; Henry et al. 2013b), that is until the arrival of JWST (Li et al. 2023).

Moreover, we remind the reader that the vast majority of such metallicity determinations are based on strong emission line ratios calibrated on a local sample of galaxies (Pettini & Pagel 2004; Maiolino et al. 2008; Curti et al. 2017) or on local analogues mimicking the physical conditions of high-z sources (Bian et al. 2018; Nakajima et al. 2022), and the applicability of these methods at high redshift still needs to be carefully assessed, as different ionisation properties in the ISM of high-z galaxies could bias the metallicity measurement (Kewley et al. 2013; Steidel et al. 2014; Strom et al. 2017; Sanders et al. 2018, 2023; Cameron et al. 2023b). Nonetheless, the rate of auroral line detections at high z is now rapidly increasing thanks to the JWST, and novel attempts to recalibrate the classical strong line diagnostics for the high-redshift Universe are being provided (Sanders et al. 2024).

In this paper, we aim to further explore the cosmic evolution of the metallicity scaling relations, leveraging deep spectroscopy with JWST/NIRSpec to probe the low-stellar-mass regime M ≈ 106.5–109.5 M from z = 3 out to the highest redshift in which rest-frame optical nebular lines are accessible to NIRSpec (z ∼ 9.5). By combining observations from JADES –which enables us to explore the very low-mass regime down to M ∼106.5– with existing datasets, we assessed the evolution of the MZR slope and normalisation, and its implications on the physical processes in place in early galaxies. Moreover, we tested the FMR framework down to the epoch when the Universe was ≲500 Myr old.

The paper is organised as follows. In Sect. 2, we describe observations, data reduction, and spectral fitting procedures. In Sect. 3, we present the analysis we performed to derive the physical quantities of interest; that is, M, SFR, and metallicity. In Sect. 4, we describe how we assessed the MZR at z >  3, combining our JADES sample with galaxies at higher masses drawn from the CEERS programme, and in Sect. 5 we investigate any evolution in the framework of the fundamental metallicity relation, and discuss possible inferences about the interplay between the gas flows and star formation regulating early galaxy assembly. In Sect. 6, we summarise our conclusions. Throughout this paper, we adopt a Planck Collaboration VI (2020) cosmology, a Chabrier (2003) initial mass function, and a solar metallicity, that is, 12+log(O/H) = 8.69 (Allende Prieto et al. 2001; Asplund et al. 2009).

2. Observations and data processing

2.1. Observations

The data presented in this paper were obtained via multi-object spectroscopy observations conducted with the micro-shutter assembly (MSA) of NIRSpec (Jakobsen et al. 2022; Ferruit et al. 2022; Böker et al. 2023) on the JWST. Observations were performed in three visits carried out between 21 and 25 October 2022 (Program ID: 1210; PI: N. Luetzgendorf) in the Great Observatories Origins Deep Survey South (GOODS-S) legacy field (Giavalisco et al. 2004), as part of one of the ‘deep’ tiers of the JWST-GTO JADES programme (Eisenstein et al. 2023; Bunker et al. 2023a). Each visit consisted of 33 613 s of integration in the PRISM/CLEAR configuration (hereinafter simply ‘PRISM’), and 8403 s integration in each of G140M/F070LP, G235M/F170LP, and G395M/F290LP (hereinafter ‘medium-resolution gratings’). Across three visits, this totals 28 h of integration in the PRISM, providing continuous spectral coverage from 0.6 to 5.3 μm at R ∼ 30 − 300, and ∼7 h in each of the medium-resolution gratings, providing R ∼ 1000 across the full spectral range of NIRSpec2.

Observations within each visit were performed adopting a 3-shutter nodding pattern, with the central pointing of each visit dithered by < 1 arcsec to sample different areas of the detector. A total of 253 unique targets were observed with the PRISM within the three pointings; among these, 67 targets were observed in all three MSA configurations, whereas 62 targets featured in two pointings, and the remaining 124 were observed for only one-third of the total exposure time. Each pointing had a bespoke MSA configuration, and target allocation was performed via the eMPT3 software (Bonaventura et al. 2023) to maximise the number of targets in common between all the three pointings, with special attention for rare objects included in the highest priority classes.

We note that in the medium and high resolution modes, individual spectra are dispersed over a large area on the detector. To minimise the possibility of overlapping emission, in our grating observations we have isolated our highest-priority targets by closing the shutters of low-priority targets on the same row (i.e. targets that could cause overlapping spectra). For this reason, in the grating modes we observe only 198 unique targets.

2.2. Data reduction

Flux-calibrated 2D and 1D spectra have been produced using the pipeline developed by the ESA NIRSpec Science Operations Team (SOT) and the NIRSpec GTO Team. Most of the processing steps in the pipelines adopt the same algorithms included in the official STScI pipeline used to generate the MAST archive products (Ferruit et al. 2022). Initially, we processed the raw data (i.e, level 1a data from the MAST archive) with a ramp-to-slope pipeline that estimates the count rate per pixel by using all unsaturated groups in the ramp, and which is optimised to reject cosmic rays on the basis of the slope of the individual ramps (for more details see Birkmann et al. 2011; Böker et al. 2012; Giardino et al. 2019; Ferruit et al. 2022). All the preprocessed count-rate images were then processed using a data reduction pipeline which includes both ESA NIRSpec SOT codes and bespoke developed NIRSpec GTO algorithms. We briefly outline here the main steps, while for a more detailed description we refer to Bunker et al. (2023a) and to a forthcoming paper of the NIRSpec/GTO collaboration (Carniani et al., in prep.). In brief, the pipeline consists of 11 main steps: (i) identification of non-target galaxies intercepting the open shutters; (ii) pixel-level background subtraction; (iii) extraction of the spectral trace of each target and wavelength and spatial coordinate assignments to each pixel in the 2D maps; (iv) pixel-to-pixel flat-field correction; v) spectrograph optics and disperser correction; (vi) absolute flux calibration; (vii) path-losses correction; (viii) rectification of 2D spectra; (ix) extraction of 1D spectra; (x) combination of 1D spectra generated from each exposure, nod, and pointing; (xi) combination of 2D spectra. Therefore, the data processing workflow returns both a combined 1D and 2D spectrum for each target. We note that the combined 1D spectra are not extracted from the combined 2D maps, but are the result of a weighted average of 1D spectra from all integrations, which allowed us to implement bad pixels and outliers rejection algorithms more efficiently. We adopt an irregular wavelength grid for the 1D and 2D spectra of the PRISM configuration, in order to avoid oversampling of the line spread function at short wavelengths (λ ∼ 1 μm). For the G140M/F070LP grating/filter configuration we extended the calibration of the spectrum up to 1.84 μm, taking into account the transmission filter throughput beyond the nominal wavelength range of this configuration (0.70 μm–1.27 μm). Finally, path-loss corrections are applied by modelling galaxies as point-like sources, taking into account the relative intra-shutter position of each source by leveraging on a full forward-modelling of the telescope and instrument optical paths (Jakobsen et al. 2022; Ferruit et al. 2022).

2.3. Spectral fitting

Both continuum and line emission are modelled simultaneously adopting the penalised pixel fitting algorithm, PPXF (Cappellari 2017, 2023), which 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. We used the high-resolution (R = 10 000) SSP library combining MIST isochrones (Choi et al. 2016) and the C3K theoretical atmospheres (Conroy et al. 2019). The flux blueward of the Lyman break was manually set to zero. These templates are complemented by a fifth-degree multiplicative Legendre polynomial in order to take into account systematic differences between the SSPs and the data (such as dust extinction, mismatch between the SSP models and high-redshift stellar populations, as well as residual flux calibration issues). The emission lines are modelled as pixel-integrated Gaussians, matching the observed spectral resolution. To reduce the number of degrees of freedom, we divide all emission lines in four kinematic groups, constrained to have the same redshift and intrinsic broadening. These are UV lines (blueward of 3000 Å), the Balmer series of hydrogen emission lines, non-hydrogen optical lines (blueward of 9000 Å), and near-infrared (NIR) lines. The stellar continuum has the same kinematics as the Balmer lines. Furthermore, we tie together doublets that have fixed flux ratios as determined by the atomic physics (such as the [O III] λ 5007/4959 ratio, which is set equal to 2.97), and constrain variable-ratio doublets to their physical ranges. We note that, at the resolution of the PRISM/CLEAR mode, the two components of the [O II] λλ 3726, 3729 doublet are completely unresolved and therefore are fit as a single Gaussian component centred at a rest-frame wavelength of 3728.42 Å, whereas for the gratings configurations the two emission lines are fitted separately, and their ratio is constrained to the physical range allowed by the density of the ISM (∈[0.38 − 1.46]); the same applies for the [S II] λλ 6718,6732 doublet (∈[0.44 − 1.45]) (Osterbrock & Ferland 2006).

We show an example fit for a z = 4.805 galaxy in Fig. 1. The upper panel shows the fit (in red) to the full PRISM (R100) spectrum in the wavelength region covering the main rest-frame optical emission lines, whereas in the three bottom panels we show a zoom-in of the best-fit to the R1000 spectrum for the same galaxy on the region of [O II] and [Ne III], Hβ and [O III], and Hα, [N II] and [S II]emission lines, respectively.

thumbnail Fig. 1.

Example of the spectral fitting procedure for one of the JADES galaxies in the sample, namely JADES-GS+53.16268-27.80237 at z = 4.805. Top panel: PRISM (R100) spectrum (in black), showing the region of the main rest-frame optical emission lines. The best-fit to the spectrum from PPXF is shown in red. Bottom panels: medium-resolution grating (R1000) spectrum (and its best-fit) for the same galaxy. From left to right, the panels show a zoom onto the regions of the [O II] and [Ne III], Hβ and [O III], and Hα [N II] and [S II] emission lines, respectively. In all panels, the error spectrum is marked by the cyan shaded region.

3. Derived physical quantities

3.1. Gas-phase metallicity

For the purposes of deriving the gas-phase metallicity, we use the observations obtained with both PRISM (R ∼ 100) and medium resolution gratings (R ∼ 1000), treating however the two spectral configurations independently. For each spectral configuration, emission line fluxes as measured from PPXF are reddening corrected on the basis of the decrement measured from the available Balmer lines. More specifically, we exploit the Hα/Hβ ratio at redshifts z <  6.75, where Hα is spectrally covered by NIRSpec, and Hγ/Hβ otherwise, adopting a Gordon et al. (2003) attenuation law and assuming the theoretical ratios from Case B recombination at T = 1.5 × 104 K (i.e. Hα/Hβ = 2.86; Hγ/Hβ = 0.47). In case no Balmer decrement could be measured, we adopt the nebular attenuation inferred from the SED fitting performed on PRISM spectra with the BEAGLE code (Chevallard & Charlot 2016, Chevallard et al., in prep.), which incorporates a two-component dust attenuation model that accounts for the differential attenuation of nebular and stellar emission in a self-consistent way.

The gas-phase metallicity is derived exploiting a revisited version of the calibrations presented in Curti et al. (2020b) (which are anchored to the Te abundance scale), refined using a sample of local, metal-poor galaxies to better sample the low-metallicity regime. Such calibration set is similar to that presented in Nakajima et al. (2022). In addition, a novel diagnostics labeled R̂ (see Table 1), and based on a linear combination of R2 and R3 which further reduces the scatter of the calibration sample at fixed metallicity and span a wider dynamical range compared to the ‘standard’ R23 diagnostic, is implemented in place of the latter. For a detailed description of the rationale behind the new diagnostic, and for a more thorough comparison of classical strong-line calibrations against individual auroral line detections within the hereby presented JADES dataset, we refer to Laseter et al. (2024). The coefficients of the adopted calibrations are tabulated in Appendix B. A minimum signal-to-noise (S/N) of 5 is required on Hα and [O III], whereas of 3 on fainter emission lines, in order to consider a given emission line (and its associated diagnostics) in the metallicity derivation.

Table 1.

Definitions of line ratios adopted throughout the paper.

We visually inspected the spectra to remove sources affected by a poor fitting of the emission lines, as well as two galaxies identified as AGNs based on evidence for broad line emission, as detailed in Maiolino et al. (2023). Moreover, additional 14 AGN candidates identified via significant emission detected in high-ionisation transitions like He IIλ4686 and N IVλ1485 (Scholtz et al. 2023) are also removed from the analysis, as AGN-powered ionisation would compromise the standard metallicity calibrations tuned for star-forming galaxies, as well as complicate the derivation of stellar masses and star-formation rates. However, it is relevant to note that most of these galaxies still show line ratios consistent with the star-forming population according to the ‘BPT’ diagram (Baldwin et al. 1981; Kewley et al. 2001; Kauffmann et al. 2003) , whose applicability in discriminating among different ionising sources at high redshift has been in fact recently questioned, especially at low metallicity (Übler et al. 2023; Maiolino et al. 2023). Considering the challenges in quantitatively assessing the relative contribution of AGN ionisation with respect to the host galaxy in the observed spectra of these sources, we apply a conservative choice by removing the entire sample of narrow-line AGN candidates. Nonetheless, all the conclusions of the present paper are robust against this choice, and in Appendix A we report the results obtained by including also this subsample of galaxies.

In this work, all available strong-line diagnostics are included for each galaxy in the metallicity calculation, in order to reduce potential biases associated with the use of an individual, specific line ratio. The procedure explores the log(O/H) parameter space with a Markov Chain Monte Carlo (MCMC) algorithm exploiting the EMCEE package (Foreman-Mackey et al. 2013), with the logarithmic likelihood defined as

log ( L ) i ( R obs , i R cal , i ) 2 ( σ obs , i 2 + σ cal , i 2 ) , $$ \begin{aligned} {\log }({L}) \propto \sum _{i} \frac{({R}_{\text{obs},i} - R_{\text{cal},i} )^2}{(\sigma ^2_{\text{obs},i} + \sigma ^2_{\text{cal},i})} , \end{aligned} $$(1)

where the sum is performed over the set of available diagnostics used, Robs are the observed line ratios, Rcal are the line ratios predicted by each calibration at a given metallicity, σobs are the uncertainties on the observed line ratios, and σcal are the dispersions of each calibrated line ratio at fixed metallicity.

For the majority of the sources under study however, we note that the [N II] and [S II] emission lines are either undetected, or falls outside of the NIRSpec wavelength coverage (see also Cameron et al. 2023b); moreover, [N II] λ6584 and Hα are severely blended in PRISM spectra at the lowest redshifts probed by our sample. Therefore, we ultimately and effectively involve only diagnostics based on ‘alpha’ elements (i.e. R3, R̂, O32; when available, [Ne III]/[O II] is also included). In some cases, [O III] λ5007 and Hβ are the only emission lines significantly (above three-sigma) detected in the spectrum, which means that R3 is the only available diagnostic. In such cases, further information from upper limits on the [N II] λ6584 (from R1000 spectra, where available) and on the [O II] λ 3727,29 emission lines can be placed on the N2 and O32 diagnostics, helping to discriminate between the two solutions provided by the double-branched R3 calibration. We here further note however that adopting instead a single, double-branched diagnostics (like R3 or R23) for the whole sample, although increasing the self-consistency of the metallicity derivation, might artificially introduce scatter in the log(O/H) distribution, as high-z galaxies typically tend to occupy the high-excitation region of the calibration, close or above the calibration plateau. This means that, even if the degeneracy between the two branches is somehow broken, the calibration itself carries less information about the metallicity in that region (being almost flat and therefore not very sensitive to metallicity variations). Therefore, forcing either the low- or high-solution introduces a ‘gap’ in metallicity between populations of galaxies which share very similar line ratios (see e.g. the discussion in Guo et al. 2016), as the offset between the low- and high-metallicity solution is large at fixed line ratio. This can have a non-negligible impact, for instance, on the inferred slope of the metallicity scaling relations.

The procedure outlined above is applied separately for each galaxy to both PRISM and gratings spectra (where the latter are available). As mentioned above, we note that the N2 diagnostic is not included in any of the metallicity calculations based on PRISM spectra, as the R100 spectral resolution is not enough to resolve [N II] λ6584 from Hα; similarly, S2-based diagnostics are not included for PRISM spectra, as they might be affected by [N II] contamination of the Hα line. Regarding the analysis carried out in the present paper, a fiducial metallicity value is then assumed for each source on the basis of the following scheme:

  • We take the O/H inferred from medium-resolution (R1000) gratings if at least two diagnostics are available; otherwise

  • we take the O/H inferred from the PRISM if at least two diagnostics are available in the PRISM spectrum; otherwise

  • we take the O/H inferred from medium-resolution gratings as based on one individual diagnostic (typically R3), exploiting the information from three-sigma upper and lower limits on other line ratios (for instance, N2 or O32) to break the degeneracy of double-valued calibrations; otherwise

  • we take the O/H inferred from the PRISM as based on one diagnostic (typically R3) in combination with the information from the upper limits on [O II] λλ3727,29.

The hereby adopted scheme leverages the higher resolution grating spectra for galaxies with strong line emission, whereas the much deeper PRISM spectra for fainter targets. We opt for not considering ‘hybrid’ gratings-to-prism line ratios because of potential issues associated with small wavelength or flux inter-calibration uncertainties between the two configurations that are not yet fully understood. However, we compare the metallicity derived from gratings and prism spectra for the 35 galaxies where both measurements are available in the right-hand panel of Fig. 2. The two measurements are in agreement within the 1σ individual uncertainties in ∼75% of the cases (and in ∼95% within 2σ), with negligible systematic offset (i.e. 0.01 dex median offset towards higher grating metallicities) and a standard deviation of 0.15 dex.

thumbnail Fig. 2.

Comparison between the metallicity derived from medium-resolution gratings and PRISM spectra. We plot only objects for which the requirements described in Sect. 3.1 are satisfied for both configurations. The two distributions scatter across the equality line (in red), with a median offset of 0.01 dex and a standard deviation of 0.15 dex. In 75% of the cases the two measurements are consistent within their 1σ uncertainties.

In this work, we focus specifically on galaxies at z ≥ 3, in order to study the evolution of the metallicity scaling relations beyond the epochs probed by previous surveys from the ground. Therefore, modulo the selections and caveats discussed above, our analysis includes a total of 62 galaxies from JADES with a metallicity measurement. The oxygen abundances measured in the selected JADES sample span between 12+log(O/H) = 7.21–8.54, with an average of 7.77 (i.e. 0.12 Z). In Fig. 3, we report the redshift distribution of our final selected sample of galaxies, compared to the distribution of the sample from the JADES deep GOODS-S tier observations with a spectroscopically confirmed redshift4, which in turn constitutes ∼70 per cent of the total sample introduced in Sect. 2.1; for a more detailed description of the full JADES deep sample, the selection function, and galaxy spectroscopic confirmation, we refer to Bunker et al. (2023a).

thumbnail Fig. 3.

Redshift distribution of the selected sample of 54 JADES galaxies with metallicity measurements analysed in this paper. Their histogram is compared to the distribution of the full ‘JADES-Deep’ sample presented in Sect. 2.1 and Bunker et al. (2023a) with a spectroscopic redshift determination.

In particular, the final sub-sample selected for the current analysis represents the ∼57 percent of the JADES deep sample in GOODS-S at z >  3 with a spectroscopic redshift, as remaining galaxies either do not have more than one emission line detected in their spectra, show evidence for contamination from AGN, or are lineless sources identified only through prominent spectral breaks (e.g. Curtis-Lake et al. 2023; Looser et al. 2023).

3.2. Stellar masses and SFRs

To measure the stellar masses for our selected sample of galaxies in JADES, we employ full spectral fitting to PRISM spectra performed with the BEAGLE code (Chevallard & Charlot 2016). The PRISM spectra that have been employed only included slit-loss corrections based on the assumption of a point source, which is not appropriate for all objects in the sample. To address this, where objects have NIRCam photometry from the JADES survey, the BEAGLE fits are run with both PRISM spectroscopy and total, Kron-based photometry (Rieke et al. 2023). In these cases, a low-order calibration polynomial is included in the fit to match the shape and normalisation of the spectrum to the photometry, whereby the total stellar mass and star-formation rate estimates are adjusted to the total fluxes measured for the objects. Where NIRCam photometry is not available we resort to the derived quantities from PRISM spectroscopy alone. We assume a delayed-exponential star-formation history (SFH), a Chabrier (2003) IMF with an upper mass limit of 100 M, and adopt the updated Bruzual & Charlot (2003) stellar population models described in Vidal-García et al. (2017). We assume the total mass currently locked into stars as our fiducial estimate for the stellar mass, which accounts for the fraction of mass returned to the ISM, instead of the integrated SFH. We refer to Chevallard et al. (in prep.) for further information of the BEAGLE fitting procedure.

The star-formation rates of our galaxies are estimated in two different ways, namely from the output of the BEAGLE SED fitting run described above, which provides an estimate of the star-formation rate averaged over the past ten Myr, and from the attenuation-corrected5 Hα luminosity, following the recipe6 outlined in Reddy et al. (2022), Shapley et al. (2023b), better suited for low-metallicity galaxies characterised by an increased rate of ionising photons. For galaxies in which Hα is not observed, either because it is shifted out of the spectral coverage of NIRSpec (i.e. at z ≳ 7), or because it falls in one of the detector gaps (in medium resolution gratings), we exploited the dust-corrected Hβ flux rescaled by the theoretical case B recombination factor at T ∼ 1.5 × 104 K (2.86). For consistency, for any given galaxy we use the Hα(Hβ) flux as measured from either the PRISM or medium resolution gratings spectra according to the configuration exploited for the full metallicity calculation (as described in Sect. 3.1).

On average, the Hα-based SFRs are ∼0.14 dex lower than those derived by BEAGLE, with a dispersion in the deviation between the two measurements of ∼0.3 dex. Such discrepancy is possibly driven by (i) the different underlying stellar population synthesis models (i.e. BPASS for Reddy et al. 2022, revised Bruzual & Charlot 2003 for BEAGLE), (ii) different prescriptions for the star-formation history, (iii) residual contamination of the Hα flux from [N II] in PRISM spectra. We adopt the SED-based measurements as our fiducial value of the total SFR of a galaxy for the analysis presented in the current paper, more consistent with the stellar mass derivation. Nonetheless, we verified that assuming either of the two SFR values does not alter any of the conclusions presented in this work, particularly when discussing the evolution in the fundamental metallicity relation (Sect. 5), as shown in Appendix A.

3.3. High-redshift galaxy samples from the literature

The specific selection function of the deep tier of the JADES spectroscopic campaign provides preferential coverage of the low-end of the stellar mass distribution (i.e. M/M ≈ 106.5–109), a regime poorly probed before, especially at such early cosmic epochs. In order to extend the parameter space over which the metallicity scaling relations are analysed, here we complement our JADES sample with a sample of 80, z = 4 − 10 galaxies from the CEERS programme as analysed and presented by Nakajima et al. (2023), which are preferentially distributed across the log(M/M)∼8.5 − 10 regime.

To preserve a good level of self-consistency in the analysis, we have recomputed the metallicity for these galaxies starting from the emission line ratios reported in Table B1 of Nakajima et al. (2023), applying the procedure discussed in Sect. 3.1. Star-formation rates reported in Nakajima et al. (2023) are based on dust-corrected Hβ luminosity and the calibration reported by Kennicutt & Evans (2012). We scale down these values by the median offset (i.e. 0.23 dex) between the Hα-based SFR (this time derived assuming the Kennicutt & Evans 2012 formula) and the BEAGLE-SFR as measured for the JADES galaxies, in order to minimise systematic offsets between the SFRs compiled for CEERS galaxies and the fiducial values assumed for JADES objects.

Stellar masses in the CEERS sample are based on fits to NIRCAM photometry performed with the PROSPECTOR code (Johnson et al. 2021). A full assessment of the systematics involved in the stellar mass determination as based on different spectral fitting codes goes beyond the scope of this paper. Nonetheless, we note that the statistical uncertainties on M quoted by Nakajima et al. (2023) are much larger than those associated with BEAGLE measurements for JADES galaxies. To reduce the impact of such differential measurement uncertainties, we do not apply inverse variance weighting based on M uncertainties in any of our fitting analyses in the paper.

In addition, we include in the analysis four targets whose metallicity have been derived with the ‘direct’ Te-method, namely three galaxies observed in the framework of the Early Release Observations programme as described in Curti et al. (2023), and the galaxy GN-z11 observed in the framework of JADES and presented in Bunker et al. (2023b), Tacchella et al. (2023a), though the nature of its ionising source is debated and the presence of an AGN has been proposed (Maiolino et al. 2024). This brings the total number of galaxies analysed in this paper to 146.

For the three ERO obejcts, M and SFR are derived via BEAGLE fitting to NIRCAM photometry as described in Curti et al. (2023), whereas the metallicities have been recomputed following the same data reduction implemented for the JADES galaxies analysed in this paper, as also described in Laseter et al. (2024). For GN-z11, fiducial values for M and SFR are taken from Bunker et al. (2023b) and are based on BEAGLE fitting to PRISM spectrum assuming a Chabrier (2003) IMF with an upper mass cut-off of 100 M, whereas the oxygen abundance is measured exploiting the detection of [O III] λ4363 in the PRISM spectrum and its ratio over [Ne III] λ3869 and Hγ, as detailed in Cameron et al. (2023a), which quote a fiducial value of 12+log(O/H) =7.82. For comparison, we also obtain an indirect estimate of the [O III] λ5007 flux (which is not covered in the spectra of GN-z11) based on the observed [Ne III]/[O II] (from medium resolution gratings) and the Witstok et al. (2021) [Ne III]/[O II] versus [O III]/[O II] calibration, which we then use to derive the R̂ diagnostic and use it conjunction with Ne3O2 and the methodology described in Sect. 3.1: this approach provides 12+log(O/H) = 7.67 + 0.16 0.10 $ 7.67\substack{+0.16 \\ -0.10} $. At the same time, using the same indirect estimate of the [O III] λ5007 flux to apply the Te method (and adopting Pilyugin et al. 2009 to derive the temperature of the low-ionisation zone) delivers a much higher 12+log(O/H) = 8.42. Given the large systematic uncertainties affecting both methods, we here assume the value reported by Cameron et al. (2023a) (7.82) as our fiducial value, with a conservative uncertainty of 0.35 dex.

Figure 4 shows the distribution of both the JADES sample introduced in this work and the compiled sample from the literature discussed above, in the stellar mass versus SFR diagram; the marginalised histograms of the distribution of the two quantities are also shown in the upper and right-hand inset panels. We note the two samples are quite complementary, and especially including galaxies from Nakajima et al. (2023) allows us to extend the mass regime probed up to M ∼ 109.5 M (⟨log(M/M)⟩CEERS = 8.72; ⟨log(M/M)⟩JADES = 7.83). Galaxies selected from CEERS are also distributed towards higher SFRs (⟨log(SFR)⟩CEERS = 1.3; ⟨log(SFR)⟩JADES = 0.45), though the two galaxy samples probe similar specific star-formation rates (sSFR = SFR/M). For reference, we plot the parametrisation of the ‘Star Forming Main Sequence’ (SFMS, Noeske et al. 2007; Speagle et al. 2014; Renzini & Peng 2015) at z ∼ 0 and z ∼ 6 from Popesso et al. (2023), noting that the steepness of such parametrisation causes its low-mass extrapolation at z ∼ 6 (dashed line) to systematically underestimate the observed star-formation rate for the low-mass JWST galaxy sample (by ∼0.5 dex on average at M = 108 M), given our assumed fiducial M and SFR measurements.

thumbnail Fig. 4.

Distribution in the stellar-mass–SFR plane for our combined JWST sample. This includes both JADES galaxies presented in this work and the literature sample compiled from Nakajima et al. (2023; CEERS), Curti et al. (2023; EROs), and Bunker et al. (2023b; GN-z11), respectively. The top and right-hand inset panels show the histograms of the distribution in both parameters. The SFR for JADES galaxies is derived from BEAGLE fitting to both PRISM spectra and NIRCAM photometry (where available, and to PRISM spectra only otherwise). For CEERS galaxies, the SFRs are compiled from Nakajima et al. (2023) as inferred from the Hβ luminosity and the Kennicutt & Evans (2012) calibration, but have been here scaled down by 0.23 dex to account for the mean offset between the Hα(Hβ)-based and BEAGLE-based SFRs measured for the JADES sample (we refer to Sect. 3.2 for more details). For both ERO galaxies and GN-z11 the M and SFR have been derived via BEAGLE fitting, consistently with what done for JADES sources. The parametrisation of the main sequence of star-formation (SFMS) at z ∼ 0 and z ∼ 6 from Popesso et al. (2023; and their extrapolation at low stellar mass, dashed line) is also shown for reference.

4. Evolution of the mass–metallicity relation beyond z = 3

4.1. The observed MZR at 3 < z < 10

In Fig. 5, we show our combined JWST galaxy sample on the mass-metallicity plane. Each point is marked with a different symbol according to its parent sample (as reported in the legend), and it is colour-coded on the basis of the following redshift binning scheme: yellow at z ∈ [3 − 6], and blue at z ∈ [6 − 10]. The median redshift of the whole sample is ⟨z⟩=5.10, whereas the median redshift of the z ∈ 3 − 6 sub-sample is ⟨z3 − 6 = 4.76 and that of the z ∈ 6 − 10 sub-sample is ⟨z6 − 10 = 6.73. In particular, filled circle points in the plot report galaxies from the JADES survey, while filled crosses symbols mark galaxies from CEERS (Nakajima et al. 2023), whose metallicity has been remeasured as described in Sect. 3.1. For both the total sample and each sub-sample in redshift we show the median (large diamond markers), error on the median (solid errorbar), and standard deviation (dashed errorbar) of the metallicity within three different bins in stellar mass. The size of the stellar mass bins is not uniform to maintain a reasonable (i.e. at least 10) number of galaxies in each bin. The average properties of the binned samples are reported in Table 2.

thumbnail Fig. 5.

Mass-metallicity relation (MZR) for our full JWST sample. Filled circles represent individual galaxies presented in this work from JADES, whereas filled crosses are galaxies observed in the framework of the CEERS programme and compiled from Nakajima et al. (2023), with metallicity derived as detailed in Sect. 3.1. The star symbol reports the JADES/NIRSpec observations of GN-z11 (Bunker et al. 2023b, at z = 10.603), whereas ‘X’ symbols mark galaxies from the EROs as compiled from Curti et al. (2023) and Laseter et al. (2024). Large squared and diamond symbols mark the median values computed in bins of M (full sample, in purple), and (M z), as described in Table 2. An orthogonal linear regression fit to the median values in bins of M for the different redshift sub-samples is shown by the purple (full sample), yellow (z = 3 − 6) and blue (z = 6 − 10) lines, respectively. We include a comparison with previous determinations of the MZR at lower redshifts from Curti et al. (2020b; SDSS at z ∼ 0.07), and Sanders et al. (2021; MOSDEF at z ∼ 2 − 3), as well as the best-fit of the low-mass end of the MZR at z ∼ 3 provided by Li et al. (2023) and based on JWST/NIRISS slitless spectroscopy. The MZR curves at z ∼ 2 − 3 have been scaled down by ∼0.1 dex to account for the systematics differences between the metallicity calibrations used in this work and the Bian et al. (2018) calibrations adopted in the original papers.

Table 2.

Median values, error on the median (and standard deviation) in stellar mass, SFR, and metallicity for our sample of galaxies in each of the (M–redshift) bins considered in this work.

Our galaxies present a median offset of ∼ − 0.5 dex (−0.48 dex and −0.64 dex at z = 3–6 and z = 6–10, respectively) compared to the low-mass end (and its extrapolation) of the MZR in the local Universe (based on individual SDSS galaxies, Curti et al. 2020b). We also compare our results with previous realisations of the MZR at z ∼ 2 − 3. In particular, we first compare our observations with the mass-metallicity relation derived from stacked spectra of galaxies from the MOSDEF Survey as reported by Sanders et al. (2021). To minimise the impact of systematics introduced by the use of different metallicity calibrations (the calibrations from Bian et al. 2018 are adopted in Sanders et al. 2021), we have recomputed the metallicity for the MOSDEF stacked spectra with the same methodology outlined in Sect. 3.1, finding an average offset of ∼0.088 dex towards lower metallicities. Therefore, in Fig. 5 the mass-metallicity relations at z ∼ 2.2 and z ∼ 3.3 from Sanders et al. (2021) are lowered by 0.088 dex compared to their original parametrisation. Overall, the evolution in normalisation probed by our full galaxy sample at z = 3 − 10 appears relatively mild if compared to the extrapolation at low stellar mass of the z ∼ 3.3 MZR from Sanders et al. (2021), with a mean offset for the full sample of 0.045 dex (0.13 dex if we assume the fiducial MZR from Sanders et al. 2021 based on the Bian et al. 2018 calibrations). However, such deviation is mass-dependent and is more prominent at higher M, being almost zero for log(M/M)≲8, while ∼0.15 dex at log(M/M)∼9. This behaviour suggests an evolution in the slope of the mass-metallicity relation in the mass and redshift regimes probed by this work.

In the attempt to characterise the evolution in the slope of the MZR at these redshifts, we perform an orthogonal linear regression fit to both the full JWST sample of individual galaxies, and within the two redshift bins separately, in the functional form

12+log(O/H) = β z log ( M / 10 8 M ) + Z m 8 , $$ \begin{aligned} \text{12+log(O/H)}=\beta _{z} \ {\log }({M}_{\star }/10^8\,{M}_{\odot }) + {Z}_{m8} \, ,\end{aligned} $$(2)

where βz is the slope at a given redshift interval, and Zm8 is the normalisation at log(M/M) = 8. The results are displayed in Fig. 5, as shown by the solid purple (full sample), yellow (z = 3 − 6), and blue (z = 6 − 10) lines, whereas the shaded regions represent the 1-σ confidence interval of each fit, derived from bootstrapping each sample (with replacement) and repeating the fitting procedure 300 times. The best-fit parameters are reported in Table 3. We note that the best-fit MZR (based on individual data points) agrees well with the median values computed in bins of M for each redshift sub-sample.

Table 3.

Best-fit parameters of the mass–metallicity relation from Eq. (2) for the full sample, and the two redshift bins at z = 3 − 6 and z = 6 − 10, for which we report both the slope βz and the normalisation at M = 108 M, Zm8.

We find indications of a flattening of the slope of the low-mass end of the MZR with redshift, with the best-fit slope for the full sample equal to β3 − 10 = 0.17 ± 0.03. If we focus more specifically on the two redshift bins described above, we find β3 − 6 = 0.18 ± 0.03, and β6 − 10 = 0.11 ± 0.05. The low-mass end slope of the MZR as probed by our sample is flatter (at 3.6σ significance) than both observed for the z ∼ 2–3 MZR at higher stellar masses (β2 − 3 = 0.29), and also in the local Universe (0.28, as measured in Curti et al. 2020b). Such apparent invariance between z ∼ 0 and z ∼ 3 suggests that the same physical processes, and in particular how the metal removal efficiency of stellar winds scales with M, governs the MZR slope over the last ∼12 Gyr of cosmic time (Sanders et al. 2021). We note however that most previous assessements of the MZR at high-z are based on near-infrared spectroscopic surveys from the ground, and were therefore limited in the stellar mass range they could probe, struggling to access the regime below log(M)< 109 M and thus to provide a sample matched in stellar mass to the JWST sample discussed in the present work. Further insights on the low-mass end slope of the MZR at z ∼ 2 − 3 come from recent JWST/NIRISS observations within the GLASS-ERS programme, have been presented in Li et al. (2023), and report similar evidence for a shallower slope (β2 − 3 = 0.17 ± 0.03, as shown in magenta in Fig. 5) compared to the simple extrapolation of the z ∼ 3.3 relation inferred by Sanders et al. (2021). In this work, we push the investigation of the MZR in the dwarf regime to much higher redshift, finding a slope at z3 − 6 which is consistent with that reported by Li et al. (2023). Compared to that parametrisation of the MZR, we observe a more prominent evolution in the normalisation of ∼0.25 dex, after once again correcting for the systematics between the Bian et al. (2018) calibrations and the set of calibrations adopted in this work.

It is also interesting to observe that the slope inferred at z = 6 − 10 (in blue in Fig. 5) is even flatter (and consistent with zero within ∼2σ), suggesting that the mass-metallicity relation could be in an initial build-up phase. However, the high-redshift bins are only sparsely populated and are probably subject to strong selection biases, and it is therefore difficult to draw any strong conclusion at this stage relative to any possible evolution in the MZR slope between z = 3 − 6 and z >  6.

Finally, we note that the distribution of our data in the MZR plane is characterised by a large amount of scatter, at any given mass and redshift. Although the sample considered in the present work is not as complete and unbiased as required to perform a robust assessment of the scatter of the scaling relation, we can tentatively measure the amount of intrinsic scatter in our galaxy sample as σ MZR = σ obs 2 σ meas 2 $ \sigma_{\text{MZR}} = \sqrt{\sigma^{2}_{\text{obs}} - \sigma^{2}_{\text{meas}}} $, where σobs is the standard deviation of the full JWST sample around the best-fit relation and σmeas is the average measurement uncertainty associated with the metallicity determination7. We find σMZR ∼ 0.073 dex for our full sample, consistent with the scatter measured in the local MZR at higher masses (log(M/M)> 9.5) (σMZR, z = 0 = 0.075; Curti et al. 2020b). We note that the scatter at high masses in local galaxies is likely dominated by the flattening of the MZR at high masses (log(M/M)≥10), interpreted in classical gas-regulator models (Finlator & Davé 2008; Lilly et al. 2013) as a consequence of the metallicity approaching the stellar yield, whereas no clear ‘turnover’ nor high-mass flattening is probed at higher redshifts.

4.2. Local analogues of high-redshift galaxies

The best-fit relations (and median values in bins of M and redshift) reported in this work show remarkable agreement with the MZR probed by samples of local star-forming galaxies with extreme line emission properties, such as the ‘Green Pea’ and ‘Blueberry’ galaxies compiled from Yang et al. (2017a,b), whose metallicity has been derived with the ‘direct’ Te method; these objects are marked by green and cyan symbols, respectively, in the left panel of Fig. 6. In particular, the MZR slope inferred from the combined sample of ‘Green Pea’ and ‘Blueberry’ galaxies is 0.18 ± 0.02, fully consistent with the slope measured from our sample, whereas JWST galaxies have, on average, 0.12 dex lower metallicity. These low-redshift sources are well matched in stellar mass with our JWST sample (log(M/M)≈6.5–9.5), are very compact, characterised by very high equivalent width of the emission lines, low metallicities, and high ionisation parameters, and therefore have been long identified as potential analogues of very high-redshift sources. The good agreement between the emission line and metallicity properties in our z ∼ 3 − 10 JWST sample and those observed in this sample of local, extreme line emitters corroborates such interpretation (see also Cameron et al. 2023b).

thumbnail Fig. 6.

Comparison with local analogues and cosmological simulations. Left-hand panel: the best-fit MZR presented in this work (for both the full JWST sample in purple, and in the z = 3 − 6 and z = 6 − 10 redshift binsin yellow and blue, respectively) is compared to the Te-based mass-metallicity relation from local, metal-poor ‘Blueberry’ and ‘Green Pea’ galaxies (Yang et al. 2017a,b), which share similar excitation and emission line properties to our high-redshift galaxy sample (Cameron et al. 2023b). Their distribution agrees well in slope and normalisation (with ∼0.11 dex offset in log(O/H)) with the average relation inferred in this work for our low-mass JWST sample. Right-hand panel: our best-fit MZRs are compared with the predictions from different suites of cosmological simulations at z ∼ 6–8, namely FIRE (Ma et al. 2016), IllustrisTNG (Torrey et al. 2019), FirstLight (Langan et al. 2020), Astraeus (Ucci et al. 2023), SERRA (Pallottini et al. 2022), and with the chemodynamical simulations from Kobayashi & Taylor (2023). In addition, the typical MZR slopes predicted by (semi-) analytical chemical evolution models (e.g. Davé et al. 2012) implementing feedback via either ‘energy-driven’ or ‘momentum-driven’ winds are also shown (the normalisation is assumed arbitrary to aid visualisation). The latter agrees well with the slope observed at low stellar mass for galaxies at z ≥ 3.

4.3. Comparison with models and simulations

Here, we compare the observed MZR in our z = 3 − 10 galaxy sample with the predictions of different cosmological simulations, as depicted in the right-hand panel of Fig. 6. In particular, we report the mass-metallicity relationship at z ∼ 6 (close to the median redshift of our full sample; i.e 5.86) from FIRE (Ma et al. 2016), IllustrisTNG (TNG100) (Torrey et al. 2019), FirstLight (Langan et al. 2020), and from the chemo-dynamical simulations by Kobayashi & Taylor (2023), as well as predictions at higher redshift (z ∼ 8) from Astraeus (Ucci et al. 2023), and Serra (Pallottini et al. 2022). For IllustrisTNG simulations we followed the approach described in Torrey et al. (2019) and already followed in Curti et al. (2023), considering central galaxies and assuming that oxygen comprises 35 per cent of the SFR-weighted metal mass fraction within twice the stellar half-mass radius. Within FIRE simulations the gas-phase metallicity is defined as the mass-weighted metallicity of all gas particles that belong to the ISM, assuming solar abundance ratios (Asplund et al. 2009). We note that the FIRE results at z = 6 are renormalised to match the normalisation of the MZR in the local Universe. Results from the FirstLight simulations at z = 6 are compiled from Langan et al. (2020, and Langan, priv. comm.), and based on the assumption that the unresolved nebular region around each star particle shares the same mass ratio of metals produced in Type II SN explosions as in the star particle. The galaxy metallicity is defined as the mass-weighted average nebular metallicity including all star particles younger than 100 Myr. The predictions from chemo-dynamical simulations by Kobayashi & Taylor (2023) are instead based on the SFR-weighted gas-phase metallicity. Finally, the Astraeus simulations (Ucci et al. 2023) accounts for the oxygen mass in the halo without any weighting, whereas in SERRA gas and stellar metallicity are coupled and the metallicity is tracked as the sum of all heavy elements (Pallottini et al. 2022).

Overall, most simulations appear to reasonably match the normalisation of the observed MZR at M ≈ 108 − 9, though under-predicting the abundances observed at log(M/M)< 8, with the theoretical MZRs characterised by steeper slopes than inferred from our sample. Only the extrapolation of IllustrisTNG predictions to lower M seems to match reasonably well (in normalisation) the average distribution of galaxies in the mass-metallicity plane below M = 108M, but diverges more dramatically at higher masses. Interestingly, the lack of a clear correlation (and the large scatter) at z ∼ 8 predicted by SERRA seems in line with our observations in the highest redshift bin (z >  6).

We also compare our inferred MZR slopes with those predicted by (semi-)analytical chemical evolution models under the assumption of different modes of stellar feedback. In particular, in the framework of ‘equilibrium’ models (e.g. Davé et al. 2012; Lilly et al. 2013), an ‘energy-driven winds’ scenario predict a slope for the MZR of βen ∼ 0.33, whereas for ‘momentum-driven winds’ the predicted slope is shallower, i.e. βmom ∼ 0.17 (Guo et al. 2016); these two cases are shown in the right-hand panel of Fig. 6 with thick black lines (here the normalisation is assumed arbitrary, and the two curves are artificially shifted along the x-axis to aid the visualisation and comparison of the various trends). We observe the former to be in good agreement with the MZR slope probed by galaxies at intermediate and high masses (log(M/M)≥9) and up to z ∼ 3, whereas the latter is more consistent with the observed MZR slope as probed by our JWST sample at the low-mass end and at higher redshift. Although preliminary, these results are consistent with a scenario in which the dominant feedback mechanism in galaxies changes in different mass regimes, producing a shallower MZR slope at low M. At the same time, observing a similar MZR slope at M ∼ 107 − 109 between z ∼ 3 (from Li et al. 2023) and z ∼ 6 (from this work) suggests that the scaling of the mass-loading factor of outflows with stellar mass is invariant in this mass regime during this epoch.

4.4. Can the evolution in the MZR normalisation be explained by evolving the gas fraction?

Another critical parameter which contributes to setting the scaling relation between M and metallicity is the gas fraction. The average gas fraction in galaxies is observed to increase with redshift at fixed M, as inferred from large existing datasets up to z ∼ 2 − 3 (e.g. Saintonge et al. 2016; Scoville et al. 2017; Tacconi et al. 2018, 2020). If the evolution proceeds at a similar rate even at earlier epochs, the impact of metal dilution is expected to become increasingly more important with redshift, overcoming the impact of metal loss due to outflows at fixed stellar mass (Sanders et al. 2021). However, this seems apparently hard to reconcile with the relatively mild evolution in the MZR normalisation observed in our JWST sample. In this section, we briefly test the impact of evolving the gas fraction on the predicted MZR evolution at low M, under a set of simple assumptions.

In the framework of analytical chemical evolution models that allow the gas fraction to vary (see e.g. Peeples & Shankar 2011; Lilly et al. 2013), and under the assumption of (almost) pristine gas inflows8, the metallicity of the ISM is the result of the balance between the nucleosynthetic stellar yield, the gas dilution, and the amount of metals lost due to outflows. Here we adopt the formalism of the Peeples & Shankar (2011) model (their Eqs. (9), (10), and (11)) to predict the evolution in the metallicity at fixed M in the low-mass regime (specifically, at M = 108M) as the consequence of the sole evolution in the gas fraction with redshift, keeping other parameters such as oxygen yields, return mass fraction, and mass-loading factor of outflows, fixed.

As for the oxygen yield yO and return fraction R9, we adopt yO = 0.032 and R = 0.463 as fiducial values, following (Vincenzo et al. 2016) and assuming a Chabrier (2003) IMF with upper-mass cutoff of 100 M, stellar yields from Nomoto et al. (2013), and metallicities typical of our JWST galaxy sample. We assume then that the redshift evolution of the gas fraction μgas (here defined as μgas = Mgas/M) follows with the same scaling up to z ∼ 6 as inferred up to z ∼ 3 (i.e. log(μgas) ∝ 2.49log(1 + z) as based on the linear parametrisation of Eq. (6) and Table 3 in Tacconi et al. 2018 10). Furthermore, in modelling the observed median MZR for the full JWST sample, we assume that galaxies lie on the SFMS (i.e. δ(MS) = 0) and on the M-Re relation from van der Wel et al. (2014) evaluated at the median redshift of the sample at log(M/M)=8 (z ∼ 5.2), the only dependence of the gas fraction being therefore on stellar mass and redshift.

We performed our calculations for two different scenarios. In the first case, we extrapolate the MZR at z ∼ 3.3 from Sanders et al. (2021; with a slope of 0.30) to log(M/M)=8 to fix the ’observed’ metallicity, and solve Eq. (10) from Peeples & Shankar (2011) for the mass-loading factor of outflows ζw, finding log(ζw)=1.46; this value is higher than what predicted by Sanders et al. (2021) if extrapolating their Eq. (16) to lower M, mainly reflecting our different choice of the fiducial yO. Then, assuming ζw invariant with redshift between z = 3 and z = 5.2 (the median redshift of the full JWST sample), we derive the expected metallicity at M = 108M by evolving the gas fraction μgas according to the scaling relations of Tacconi et al. (2018). We predict an evolution in metallicity of 0.17 dex, larger than the (almost negligible) evolution of 0.02 dex probed by the extrapolation of the Sanders et al. (2021) MZR at z ∼ 3.3.

In the second case, we instead solve for ζw in order to match slope (0.17) and normalisation (12+log(O/H) = 7.98 at M = 108M) of the MZR at z = 3 in the dwarf regime from Li et al. (2023), finding log(ζw) = 1.16. Fixing again ζw with redshift, we predict a decrease in log(O/H) as driven purely by an evolution in the gas fraction of 0.22 dex, which is in good agreement with that observed in such mass regime between the z = 3 MZR from Li et al. (2023) and our median MZR at z = 5.2 (i.e. 0.26 dex).

Within the described framework, an increase in the gas reservoir of low-mass galaxies between redshifts ∼3 and ∼5 is almost sufficient to entirely explain the (relatively mild) observed evolution of the MZR normalisation, without any adjustment to the contribution of the metal-loading factor of outflows with redshift at fixed M (whereas the scaling of ζw with M might be varying, as reflected in the apparent flattening of the MZR slope). Indeed, solving directly for ζw at z = 5.2 by imposing the measured slope and normalisation of the MZR at log(M/M)=8 as derived in this work delivers log(ζw) = 1.22, only slightly higher than the observational constraints inferred from the MZR at z ∼ 3 from Li et al. (2023) (log(ζw)=1.16, as derived above). On the other hand, an extrapolation of the Sanders et al. (2021) MZR in the low-mass regime (i.e. fixing the MZR slope) would be consistent with almost no evolution in the normalisation, requiring either a much shallower evolution in the μgas versus redshift scaling relation at z >  3, an increase in the oxygen yields, and/or a non-negligible metallicity of the infalling gas compared to that of the ISM in order to reconcile the model-predicted MZR with the observed median MZR at log(M/M)=8.

There are nonetheless a number of strong assumptions in place underlying this simple modelling. First, in analytical chemical evolution models the predicted normalisation of the MZR is highly sensitive to the assumed oxygen yield and the shape of the upper mass cutoff of the IMF (Vincenzo et al. 2016). For instance, assuming different values for the oxygen yield in the range 0.015–0.04 (which are typical values spanned in the literature for different IMFs) propagates into a shift of + 0.04 0.14 $ \substack{+0.04 \\ -0.14} $ dex in the predicted log(O/H) at M = 108M. Furthermore, we have assumed that the gas mass budget is dominated by hydrogen in the molecular phase and took Mgas = MH2 (Tacconi et al. 2018). Although different from what is observed in relatively higher mass star-forming galaxies in the local Universe (Saintonge et al. 2017; Catinella et al. 2018), there is evidence for the HI mass to be the sub-dominant component already at z ∼ 2.5 (contributing no more than ∼20% to the total gas mass), as based on dynamical mass arguments (Wuyts et al. 2016; Price et al. 2016). In addition, we have ignored the fraction of metals expelled from outflows that might be reaccreted onto galaxies via galactic fountains from the enriched CGM (Fraternali & Binney 2008; Anglés-Alcázar et al. 2017; Péroux et al. 2020), although accreted metals are generally modelled to be relevant for setting the MZR shape only at much lower redshift than explored here (Davé et al. 2011). Finally, systematics in the adopted metallicity calibrations can also shift the normalisation of the observed MZR by non-negligible (∼0.1–0.15 dex) amounts. We defer a more detailed modelling of the MZR evolution in the framework of analytical chemical evolution models, as well as a proper assessment of the involved systematics, to a forthcoming paper of this collaboration, which will exploit the improved statistics provided by the full JADES dataset.

5. Possible evidence of evolution in the fundamental metallicity relation

It is known that the observed distribution of galaxies on the mass-metallicity plane reflects selection effects associated with the average star-formation rate probed by the samples under study (Yates et al. 2012; Salim et al. 2015; Telford et al. 2016; Cresci et al. 2019). In this sense, a least biased picture of the chemical evolution stage of these early galaxies might be given by considering the full relationship between M, O/H, and SFR, in what is usually referred to as the fundamental metallicity relation (FMR; Mannucci et al. 2010), which takes into account the impact of the current SFR on the expected level of chemical enrichment of the ISM, as driven by the interplay between gas accretion, star formation, and outflows.

In Fig. 7 we report the deviation of the measured metallicity in our sample, as a function of redshift, from the predictions of the locally defined FMR as parametrised in Curti et al. (2020b). Individual points in Fig. 7 are colour-coded by their specific star-formation-rate (sSFR=SFR/M). The error bars on the y-axis includes not only the uncertainty on the individual metallicity measurements, but also an additional term (co-added in quadrature) associated with the FMR predictions as derived from the standard deviation of the distribution of log(O/H) obtained by varying one hundred times the input M and SFR of each galaxy within their uncertainties. The weighted median values in three bins of redshift (z ∈ [3 − 5]; z ∈ [5 − 7]; z ∈ [7 − 10]) for our combined JWST sample are shown by larger diamond markers to highlight possible trends, with solid and dashed errorbars representing the error on the median (as derived from bootstrapping by varying each point 300 times within its uncertainty) and the dispersion of the sample in each bin, respectively. The average offset from the FMR as measured by galaxy samples at lower redshifts (z <  3) as compiled from the literature (SDSS at z ∼ 0, Curti et al. 2020b; ‘Blueberry’ and ‘Green Pea’ galaxies at z ∼ 0.2, Yang et al. 2017a,b; KLEVER at z ∼ 1.5 − 2.5, Curti et al. 2020a; Hayden-Pawson et al. 2022; MOSDEF at z ∼ 2 − 3, Sanders et al. 2021) are also shown with different symbols as outlined in the legend, and colour-coded by the average sSFR value of the relative sample.

thumbnail Fig. 7.

Deviations of the JWST sample from the predictions of the local FMR plotted as a function of redshift. Marker symbols are as in Fig. 5, with each galaxy colour-coded by its sSFR. The weighted average of the deviation in three redshift bins (colour-coded by the mean sSFR of the binned sample) for both the JWST sample and for galaxy samples at lower redshift (SDSS at z ∼ 0, KLEVER at z ∼ 1.5 − 2.5, MOSDEF at z ∼ 3) are marked by larger symbols. Dashed error bars report the dispersion in Δlog(O/H) of each binned sample, whereas thick, solid bars report the standard error on the mean (which is often smaller than the marker’s size). Overall, z >  3 galaxies observed by JWST are scattered across the local FMR predictions, but show evidence of being preferentially offset towards lower metallicity values than expected with increasing redshift, with an average offset of −0.4 dex at z >  6.

Individual galaxies are mostly offset below the FMR at any redshift z >  3, with 83 over 146 sources deviating from the FMR predictions at more than 1σ (48 percent of galaxies at ⟨z⟩∼4, 54 percent at ⟨z⟩∼6, and 87 percent at ⟨z⟩∼8), and 25 at more than 3σ (19 percent of galaxies at ⟨z⟩∼4, 14 percent at ⟨z⟩∼6, and 19 percent at ⟨z⟩∼8). We perform a two-sample Kolmogorov-Smirnov (KS) test on the two distributions, namely the JWST galaxy sample and the full SDSS sample used to parametrise the FMR from Curti et al. (2020b; which is well approximated by a normal distribution with mean = 0 and σ = 0.054 dex), which are found to be significantly different (p-value ≪ 0.01), both for the global sample and within all individual redshift bins.

Looking more closely at the trends in the redshift bins, the median Δlog(O/H) at ⟨z⟩∼4 is = − 0.23 ± 0.04 dex (and 1 − σ dispersion of 0.29 dex), the median deviation at ⟨z⟩∼4 is = − 0.34 ± 0.05 (with a dispersion of 0.25 dex), whereas the median deviation at ⟨z⟩∼8 is =−0.50 ± 0.10 (with a dispersion of 0.38 dex). A Z-test confirms that these are all significantly offset at more than 5σ from the FMR predictions. We note however that we are here neglecting any contribution from additional systematic uncertainties on the measured quantities (M, SFR, O/H), as well as further uncertainties associated with the extrapolation of the FMR parametrisation outside the parameter space of its calibration sample. However, including a further 0.25 dex systematic uncertainty on the individual metallicity predictions evaluated at the average M and SFR of the JWST sample (following Fig. 11 in Curti et al. 2020b), the significance of the deviation of the median values from the predictions of the local FMR only slightly decreases to ∼4σ in each redshift bin.

Furthermore, we identify a tentative trend in which galaxies are seen to sit preferentially below the FMR predictions with increasing redshift, with in particular sources at z >  6 that are significantly less enriched than their M and SFR would predict were they local galaxies. A Spearman correlation test between Δlog(O/H) and redshift z on the full JWST sample provides evidence for a weak, though significant negative correlation, with ρ = −0.18, p-value = 0.03. Moreover, high redshift galaxies are also characterised, on average, by higher specific star-formation rates (sSFR) than lower redshift systems, forming stars at a higher pace while rapidly building-up their mass. A Spearman correlation test between Δlog(O/H) and sSFR indeed reveals a moderate (ρ = 0.29), significant correlation (p-value < 0.01).

Finally, in the left-hand panel of Fig. 8 we plot the same metallicity deviation from the FMR as a function of stellar mass, reporting both individual galaxies and median values in the same M-z bins described in Table 2, here colour-coded by their (median) redshift. Although no clear trend can be identified from individual points, with galaxies showing different levels of offset from the FMR regardless of their stellar mass, from the median values we find weak evidence for higher redshift galaxies to be preferentially offset from the FMR at fixed stellar mass, especially in the two highest M bins, where the offset between median points at different redshift is significant at 2.3σ and 1.2σ, respectively. No clear trend is instead identified with M, at fixed redshift.

thumbnail Fig. 8.

FMR offset as a function of stellar mass and the projection of minumum scatter. Left-hand panel: deviations from the FMR plotted as a function of stellar mass. Symbols are as in Figs. 5 and 7, with individual points (and M–redshift bins) colour-coded according to redshift. At fixed stellar mass, JWST galaxies are preferentially offset towards lower O/H with incresing redshift, and a similar trend is also seen with increasing M in the highest redshift bin. Right-hand panel: two-dimensional projection of the M-SFR-O/H space in the O/H vs μ = log(M)-αlog(SFR) plane, assuming the parametrisation for high sSFR (log(sSFR/yr−1) > −9.5) galaxies inferred from the analysis of the SDSS sample (α = 0.65, Curti et al. 2020b). Most of the individual galaxies, as well as the median values in the M-redshift bins, are offset below the relation, further suggesting that the local FMR parametrisation does not hold at the very high-z and low-M probed in this work.

Although a more statistically representative sample of early galaxies assembled by JWST is required to either robustly confirm or revisit such trends, these observations already suggest a few possible interpretations. On the one hand, they might reveal prominent accretion of pristine gas at high-redshift, happening on timescales shorter than the gas depletion, star-formation, ISM enrichment and mixing timescales compared to lower redshift galaxies, increasing the gas reservoir and diluting the galaxy metal content at fixed M and SFR (Dekel et al. 2009; Lilly et al. 2013; Somerville & Davé 2015; Davé et al. 2017). This could also be reflected into an increased level of stochasticity in the star-formation history during such early phases of galaxy formation, especially at low stellar masses. If the assembly of early galaxies is primarily driven by stochastic gas accretion from the cosmic web, as often modelled with a scatter term in the baryonic accretion rate parameter (e.g. Forbes et al. 2014), the relatively longer timescales over which the star-formation and enrichment from supernovae can balance the variations in the gas reservoir would make a system more likely to be observed out of the equilibrium which is at the baseline of the scaling relations observed at lower redshift (Davé et al. 2011; Zahid et al. 2012). At the same time, the efficiency of metal removals could also be enhanced, especially in galaxies with high specific star-formation rate, such that the metal loading factor of outflows is higher, at fixed SFR, than that of galaxies following the observed FMR at lower redshifts (Baker et al. 2023).

A possible way to differentiate between these different scenarios is to study the morphology of galaxies with the position relative to the FMR. For instance, Tacchella et al. (2023b) analysed the three galaxies at z ∼ 8 from the EROs, finding that the galaxy with the lowest gas-phase metallicity is very compact and consistent with rapid gas accretion. The two other objects with relatively higher gas-phase metallicity show more complex multi-component morphologies on kpc scales, indicating that their recent increase in star-formation rate is driven by mergers or internal gravitational instabilities. More recently, Langeroodi & Hjorth (2023) investigated the correlation between the deviation from the FMR at high redshift and the ‘compactness’ of galaxies, finding more compact sources to be preferentially offset, and interpreting this result as further evidence for enhanced stochasticity in the accretion and star-formation episodes, happening on timescales much shorter than those associated with the ISM enrichment and therefore diluting the metal content of such compact galaxies more than the average scaling relation would predict.

As seen from a different perspective, the apparent inconsistency between the measured oxygen abundance and the FMR predictions reveals the inadequacy of the local formalism when extrapolated to such early epochs and low stellar masses. Within the framework originally described in Curti et al. (2020b) in fact, the SFR dependence is only embedded in the ‘turnover mass’ term (i.e. the value in M at which the slope flattens), which is described to scale linearly with the star-formation rate. Therefore, no direct dependence of the MZR slope on the SFR is included in such FMR parametrisation, nor it is, by definition, in any of the two-dimensional projections explored in the literature in the form μα = log(M) − αlog(SFR) versus log(O/H) (e.g. Mannucci et al. 2010; Andrews & Martini 2013), which implicitly assumes no SFR-dependence of the MZR slope. In fact, a rotation of the M and SFR axis in the 3D space removes the apparent secondary dependence of the MZR on the star-formation rate only if the slopes of the different mass-metallicity relations are invariant for samples with different average star-formation rates. Then, if galaxies follow the same scaling relations between M, SFR, and O/H over cosmic time (as in the original FMR framework), such a projection is capable to cancel out any apparent evolution of the MZR with redshift. On the contrary, if such an assumption does not hold anymore, and the low-mass end slope of the MZR is SFR-dependent, then a single 2D projection is not capable of capturing the full interplay between gas flows, star-formation an chemical enrichment which manifests in the existence of a relationship between the observational quantities M, SFR, and O/H.

In Curti et al. (2020b), such a possibility was briefly discussed, and two different regimes were identified and explored, providing two different parametrisations of the FMR based on the sSFR of galaxies. In this scenario, the slope of the MZR is still invariant with SFR, but the relative strength of the SFR-dependence is significantly different in the two regimes, being much more prominent for high sSFR galaxies, which translates into two different μα projections, with α = 0.22 at low sSFR, and α = 0.65 at high sSFR11. However we stress that, being based on SDSS galaxies, these parametrisations are still statistically dominated by galaxies at higher M compared to the regime explored in the present work. For reference, we plot our combined JWST sample presented in this work in the O/H versus μα projection, for α = 0.65, in the right-hand panel of Fig. 8: such relation, though calibrated on the highest sSFR galaxies in the SDSS, is not yet capable to account for the full redshift evolution in the MZR as seen at z = 3 − 10, and most of the objects (as well as their median binned values) clearly fall below the best-fit line (in black).

Adopting different realisations of the FMR can impact the level of offset of high-redshift galaxies (and its significance), but does not cancel out nor invert the observed trend. In this sense, a comparison between the observed metallicities for the JWST sample and the FMR parametrisation from Andrews & Martini (2013) is presented and discussed in Appedix A.

6. Summary and conclusions

We analysed the metallicity properties of a sample of low-mass (log(M/M)≲8.5) galaxies at 3 <  z <  10, exploiting deep NIRSpec spectroscopic observations from the JADES programme. Detection of multiple emission lines in both PRISM and medium-resolution grating spectra allowed us to derive the gas-phase metallicity in these sources, which we complemented with a sample of sources at log(M/M)∼8.5 − 10 from other programmes and compiled from the literature (CEERS, Nakajima et al. 2023; ERO, Curti et al. 2023; and GN-z11, Bunker et al. 2023b) to provide insights into the cosmic evolution of the metallicity scaling relations in the low-mass regime up to the epoch of early galaxy assembly. Our main findings can be summarised as follows:

  • We find evidence for a mild evolution in the normalisation of the MZR at z >  3, as previously suggested by similar studies (Langeroodi et al. 2023; Heintz et al. 2023; Matthee et al. 2023; Nakajima et al. 2023; Shapley et al. 2023a). Our sample of galaxies spans a wide range in 12+log(O/H) (corresponding to Z ∼ 0.03 − 0.6 Z), with a median metallicity of 12+log(O/H) = 7.77 (Z = 0.12 Z). Compared to a simple extrapolation in the low-M regime of the MZR z ∼ 3.3 from (Sanders et al. 2021), the JWST sample shows an average offset in log(O/H) of only 0.05 dex towards lower metallicity (negligible at log(M/M)≲8, more prominent at log(M/M)∼9).

  • An orthogonal distance regression fit performed on the full JWST sample suggests a flattening of the slope at the low-mass end of the MZR at z = 3 − 10 (Fig. 5). This is in agreement with previous findings at z ∼ 3 in the same stellar-mass regime as reported by Li et al. (2023), compared to which we nevertheless measure a more prominent evolution in the MZR normalisation of ∼0.25 dex. The flattening of the MZR slope at low masses could indicate a change in the dominant mechanisms regulating supernovae-driven outflows compared to higher mass galaxies. The slope measured for our combined sample of galaxies at z = 3 − 10, β3 − 10 = 0.17 ± 0.03, is in good agreement with predictions from chemical evolution models implementing a ‘momentum-driven’ feedback mode for SNe winds (right panel of Fig. 6).

  • Theoretical simulations of the high-redshift Universe generally predict steeper MZR slopes than observed, broadly matching the normalisation of the relation at log(M/M)∼8 − 9 but struggling to reproduce the metallicity at the lowest masses probed by the present work (right panel of Fig. 6). In the framework of analytical chemical evolution models (e.g. Peeples & Shankar 2011), evolving the gas fraction with redshift according to the same scaling relations inferred up to z ∼ 3 (Tacconi et al. 2018) would almost entirely account (fixing all other parameters) for the evolution in the MZR normalisation at log(M/M)∼8 as observed between z ∼ 3 and z ∼ 5.2, once the flattening of the MZR slope in this mass and redshift regime is considered.

  • Notably, our galaxy sample is distributed over a similar region (and with a fully consistent inferred slope) of the mass–metallicity plane as low-redshift, low-M, metal-poor systems with extreme ionisation conditions, such as the ‘Blueberry’ and ‘Green Pea’ galaxies (Yang et al. 2017a,b) (left panel of Fig. 6). These objects have been identified as local analogues of high-redshift galaxies, and we observe that they share similar metallicity properties with those pertaining to the high-z sources analysed in this paper (see also Cameron et al. 2023b).

  • We find evidence for a deviation of our JWST sample from the framework of the fundamental metallicity relation (FMR) as described for local galaxies. In particular, galaxies appear significantly metal deficient compared to local galaxies matched in M and SFR, and this trend is observed more prominently as increased redshift (Fig. 7). This suggests that in these objects metals are either more efficiently removed by SNe-driven outflows, or that the gas is strongly diluted by stochastic accretion of (nearly) pristine gas, happening on timescales shorter than those associated with the enrichment from subsequent star formation.

  • This behaviour also highlights potential inconsistencies in the assumed FMR parametrisation at very high z, with the strength of the metallicity–SFR dependence at fixed M being no longer redshift invariant. In fact, if the low-mass end slope of the MZR depends on the average SFR of galaxies, a simple two-dimensional projection of the FMR in the log(O/H) versus μα plane is not capable of capturing the metallicity variations as caused by the interplay of dilution, enrichment, and metal-loaded outflows driven by star formation (right panel of Fig. 8).

The observations presented in this paper represent a step towards understanding the chemical enrichment processes in place during early galaxy formation. Deep JWST observations have already proven to be capable of extending the mass and luminosity regimes probed by previous facilities by orders of magnitudes, extending our investigations up to the highest redshift galaxies ever discovered. Improving the statistical robustness of these samples in the near future will be crucial for forthcoming studies in order to further reduce the large systematic uncertainties and to provide a more comprehensive picture of galaxy formation across different cosmic epochs, a goal fully within the reach of JWST capabilities.


1

As traced by the oxygen abundance, and usually quoted as 12+log(O/H).

2

Further 7 h were spent in the high-resolution G395H/F290LP grating, providing R ∼ 2700 between ∼2.8 − 5.1 μm; however, we are not exploiting such observations in this specific work.

4

i.e. with a redshift confirmed from either R1000 or R100 spectra.

5

As described in Sect. 3.1, we exploited the measured Balmer decrements with a Gordon et al. (2003) law where available, whereas the output Av from BEAGLE otherwise.

6

Assuming a conversion factor of 10−41.67 (M yr−1)/(erg s−1).

7

The dispersion of the adopted strong-line calibrations are included directly at the metallicity derivation step.

8

Or that, at least, the metallicity of the infalling gas is much lower than that of the outflowing gas.

9

Defined as the total mass fraction returned into the ISM by a stellar generation.

10

We note that assuming the quadratic form of the scaling relation between μgas and redshift instead would produce a turnover when extrapolating at z >  3, that is the gas fraction would be decreasing with redshift.

11

The best-fit α value for the global sample is α = 0.55.

Acknowledgments

We thank the referee for the insightful comments that helped improving this paper. MC acknowledges support from the ESO Fellowship Programme. MC, RM, FDE, TJL, JW, JS, LS, and WB acknowledge support by the Science and Technology Facilities Council (STFC), ERC Advanced Grant 695671 “QUENCH”. JW also acknowledges support from the Fondation MERAC. SC acknowledges support by European Union’s HE ERC Starting Grant No. 101040227 – WINGS. ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1). JC, AJB, AJC, and GCJ acknowledge funding from the “FirstGalaxies” Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 789056). HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. SA and BRDP acknowledge support from the research project PID2021-127718NB- I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI). RS acknowledges support from a STFC Ernest Rutherford Fellowship (ST/S004831/1). DJE is supported as a Simons Investigator and by JWST/NIRCam contract to the University of Arizona, NAS5- 02015. BER, BDJ, EE, MR, and CNAW acknowledge support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015. The work of RH is supported by the Johns Hopkins University, Institute for Data Intensive Engineering and Science (IDIES). 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. KB acknowledges support from the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. MC is grateful to Kimihiko Nakajima for making his measurements of line ratios from the CEERS sample public. MC also thanks Ivanna Langan and Daniel Ceverino for sharing the data from FIRSTLIGHT simulations. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1180 and #1210.

References

  1. Allende Prieto, C., Lambert, D. L., & Asplund, M. 2001, ApJ, 556, L63 [Google Scholar]
  2. Andrews, B. H., & Martini, P. 2013, ApJ, 765, 140 [NASA ADS] [CrossRef] [Google Scholar]
  3. Anglés-Alcázar, D., Faucher-Giguère, C.-A., Kereš, D., et al. 2017, MNRAS, 470, 4698 [CrossRef] [Google Scholar]
  4. Arellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 940, L23 [CrossRef] [Google Scholar]
  5. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  6. Baker, W. M., & Maiolino, R. 2023, MNRAS, 521, 4173 [CrossRef] [Google Scholar]
  7. Baker, W. M., Maiolino, R., Belfiore, F., et al. 2023, MNRAS, 519, 1149 [Google Scholar]
  8. Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5 [Google Scholar]
  9. Belli, S., Jones, T., Ellis, R. S., & Richard, J. 2013, ApJ, 772, 141 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bian, F., Kewley, L. J., & Dopita, M. A. 2018, ApJ, 859, 175 [CrossRef] [Google Scholar]
  11. Birkmann, S. M., Böker, T., Ferruit, P., et al. 2011, SPIE Conf. Ser., 8150, 81500B [NASA ADS] [Google Scholar]
  12. Böker, T., Birkmann, S., de Marchi, G., et al. 2012, SPIE Conf. Ser., 8442, 84423F [Google Scholar]
  13. Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, PASP, 135, 038001 [CrossRef] [Google Scholar]
  14. Bonaventura, N., Jakobsen, P., Ferruit, P., Arribas, S., & Giardino, G. 2023, A&A, 672, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023a, A&A, submitted [arXiv:2306.02467] [Google Scholar]
  17. Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023b, A&A, 677, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Cameron, A. J., Katz, H., Rey, M. P., & Saxena, A. 2023a, MNRAS, 523, 3516 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023b, A&A, 677, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Cappellari, M. 2017, MNRAS, 466, 798 [Google Scholar]
  21. Cappellari, M. 2023, MNRAS, 526, 3273 [NASA ADS] [CrossRef] [Google Scholar]
  22. Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875 [NASA ADS] [CrossRef] [Google Scholar]
  23. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  24. Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415 [NASA ADS] [CrossRef] [Google Scholar]
  25. Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
  26. Conroy, C., Naidu, R. P., Zaritsky, D., et al. 2019, ApJ, 887, 237 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cresci, G., Mannucci, F., & Curti, M. 2019, A&A, 627, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Curti, M., Cresci, G., Mannucci, F., et al. 2017, MNRAS, 465, 1384 [Google Scholar]
  29. Curti, M., Maiolino, R., Cirasuolo, M., et al. 2020a, MNRAS, 492, 821 [CrossRef] [Google Scholar]
  30. Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020b, MNRAS, 491, 944 [Google Scholar]
  31. Curti, M., D’Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425 [Google Scholar]
  32. Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nat. Astron., 7, 622 [NASA ADS] [CrossRef] [Google Scholar]
  33. Davé, R., Finlator, K., & Oppenheimer, B. D. 2011, MNRAS, 416, 1354 [CrossRef] [Google Scholar]
  34. Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98 [Google Scholar]
  35. Davé, R., Rafieferantsoa, M. H., Thompson, R. J., & Hopkins, P. F. 2017, MNRAS, 467, 115 [NASA ADS] [Google Scholar]
  36. Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451 [Google Scholar]
  37. Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, ApJS, submitted, [arXiv:2306.02465] [Google Scholar]
  38. Ellison, S. L., Patton, D. R., Simard, L., & McConnachie, A. W. 2008, ApJ, 672, L107 [NASA ADS] [CrossRef] [Google Scholar]
  39. Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813 [NASA ADS] [CrossRef] [Google Scholar]
  40. Ferruit, P., Jakobsen, P., Giardino, G., et al. 2022, A&A, 661, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13 [NASA ADS] [CrossRef] [Google Scholar]
  42. Finlator, K., & Davé, R. 2008, MNRAS, 385, 2181 [NASA ADS] [CrossRef] [Google Scholar]
  43. Forbes, J. C., Krumholz, M. R., Burkert, A., & Dekel, A. 2014, MNRAS, 443, 168 [NASA ADS] [CrossRef] [Google Scholar]
  44. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  45. Fraternali, F., & Binney, J. J. 2008, MNRAS, 386, 935 [CrossRef] [Google Scholar]
  46. Giardino, G., Birkmann, S., Robberto, M., et al. 2019, PASP, 131 [Google Scholar]
  47. Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93 [NASA ADS] [CrossRef] [Google Scholar]
  48. Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279 [NASA ADS] [CrossRef] [Google Scholar]
  49. Guo, Y., Koo, D. C., Lu, Y., et al. 2016, ApJ, 822, 103 [NASA ADS] [CrossRef] [Google Scholar]
  50. Hayden-Pawson, C., Curti, M., Maiolino, R., et al. 2022, MNRAS, 512, 2867 [NASA ADS] [CrossRef] [Google Scholar]
  51. Heintz, K. E., Brammer, G. B., Giménez-Arteaga, C., et al. 2023, Nat. Astron., 7, 1517 [NASA ADS] [CrossRef] [Google Scholar]
  52. Henry, A., Martin, C. L., Finlator, K., & Dressler, A. 2013a, ApJ, 769, 148 [NASA ADS] [CrossRef] [Google Scholar]
  53. Henry, A., Scarlata, C., Domínguez, A., et al. 2013b, ApJ, 776, L27 [NASA ADS] [CrossRef] [Google Scholar]
  54. Hirschauer, A. S., Salzer, J. J., Janowiecki, S., & Wegner, G. A. 2018, AJ, 155, 82 [NASA ADS] [CrossRef] [Google Scholar]
  55. Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22 [NASA ADS] [CrossRef] [Google Scholar]
  57. Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54 [Google Scholar]
  58. Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [NASA ADS] [CrossRef] [Google Scholar]
  59. Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121 [Google Scholar]
  60. Kewley, L. J., Maier, C., Yabe, K., et al. 2013, ApJ, 774, L10 [Google Scholar]
  61. Kobayashi, C., & Taylor, P. 2023, arXiv e-prints [arXiv:2302.07255] [Google Scholar]
  62. Langan, I., Ceverino, D., & Finlator, K. 2020, MNRAS, 494, 1988 [NASA ADS] [CrossRef] [Google Scholar]
  63. Langeroodi, D., & Hjorth, J. 2023, arXiv e-prints [arXiv:2307.06336] [Google Scholar]
  64. Langeroodi, D., Hjorth, J., Chen, W., et al. 2023, ApJ, 957, 39 [NASA ADS] [CrossRef] [Google Scholar]
  65. Laseter, I. H., Maseda, M. V., Curti, M., et al. 2024, A&A, 681, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Lee, H., Skillman, E. D., Cannon, J. M., et al. 2006, ApJ, 647, 970 [NASA ADS] [CrossRef] [Google Scholar]
  67. Lequeux, J., Peimbert, M., Rayo, J. F., Serrano, A., & Torres-Peimbert, S. 1979, A&A, 80, 155 [Google Scholar]
  68. Li, M., Cai, Z., Bian, F., et al. 2023, ApJ, 955, L18 [CrossRef] [Google Scholar]
  69. Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119 [NASA ADS] [CrossRef] [Google Scholar]
  70. Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2023, Nature, submitted [arXiv:2302.14155] [Google Scholar]
  71. Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140 [NASA ADS] [CrossRef] [Google Scholar]
  72. Maier, C., Lilly, S. J., Ziegler, B. L., et al. 2014, ApJ, 792, 3 [NASA ADS] [CrossRef] [Google Scholar]
  73. Maiolino, R., & Mannucci, F. 2019, A&A Rev., 27, 3 [NASA ADS] [CrossRef] [Google Scholar]
  74. Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023, A&A, submitted, [arXiv:2308.01230] [Google Scholar]
  76. Maiolino, R., Scholtz, J., Witstok, J., et al. 2024, Nature, 627, 59 [NASA ADS] [CrossRef] [Google Scholar]
  77. Mannucci, F., Cresci, G., Maiolino, R., et al. 2009, MNRAS, 398, 1915 [Google Scholar]
  78. Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115 [NASA ADS] [CrossRef] [Google Scholar]
  79. Mascia, S., Pentericci, L., Calabrò, A., et al. 2023, A&A, 672, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2023, ApJ, 950, 67 [NASA ADS] [CrossRef] [Google Scholar]
  81. Nakajima, K., & Ouchi, M. 2014, MNRAS, 442, 900 [Google Scholar]
  82. Nakajima, K., Ouchi, M., Xu, Y., et al. 2022, ApJS, 262, 3 [CrossRef] [Google Scholar]
  83. Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33 [NASA ADS] [CrossRef] [Google Scholar]
  84. Noeske, K. G., Faber, S. M., Weiner, B. J., et al. 2007, ApJ, 660, L47 [NASA ADS] [CrossRef] [Google Scholar]
  85. Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457 [CrossRef] [Google Scholar]
  86. Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129 [NASA ADS] [CrossRef] [Google Scholar]
  87. Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, 2nd edn. (University Science Books) [Google Scholar]
  88. Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621 [NASA ADS] [Google Scholar]
  89. Peeples, M. S., & Shankar, F. 2011, MNRAS, 417, 2962 [NASA ADS] [CrossRef] [Google Scholar]
  90. Péroux, C., Nelson, D., van de Voort, F., et al. 2020, MNRAS, 499, 2462 [CrossRef] [Google Scholar]
  91. Pettini, M., & Pagel, B. E. J. 2004, MNRAS, 348, L59 [NASA ADS] [CrossRef] [Google Scholar]
  92. Pilyugin, L. S., Mattsson, L., Vílchez, J. M., & Cedrés, B. 2009, MNRAS, 398, 485 [NASA ADS] [CrossRef] [Google Scholar]
  93. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14 [NASA ADS] [CrossRef] [Google Scholar]
  95. Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526 [Google Scholar]
  96. Price, S. H., Kriek, M., Shapley, A. E., et al. 2016, ApJ, 819, 80 [Google Scholar]
  97. Reddy, N. A., Topping, M. W., Shapley, A. E., et al. 2022, ApJ, 926, 31 [NASA ADS] [CrossRef] [Google Scholar]
  98. Renzini, A., & Peng, Y.-J. 2015, ApJ, 801, L29 [Google Scholar]
  99. Rhoads, J. E., Wold, I. G. B., Harish, S., et al. 2023, ApJ, 942, L14 [NASA ADS] [CrossRef] [Google Scholar]
  100. Rieke, M. J., Robertson, B., Tacchella, S., et al. 2023, ApJS, 269, 16 [NASA ADS] [CrossRef] [Google Scholar]
  101. Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, Nat. Astron., 7, 611 [NASA ADS] [CrossRef] [Google Scholar]
  102. Saintonge, A., Catinella, B., Cortese, L., et al. 2016, MNRAS, 462, 1749 [Google Scholar]
  103. Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22 [Google Scholar]
  104. Salim, S., Lee, J. C., Dav’e, R., & Dickinson, M. 2015, ApJ, 808, 25 [NASA ADS] [CrossRef] [Google Scholar]
  105. Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2018, ApJ, 858, 99 [NASA ADS] [CrossRef] [Google Scholar]
  106. Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19 [CrossRef] [Google Scholar]
  107. Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023, ApJ, 955, 54 [NASA ADS] [CrossRef] [Google Scholar]
  108. Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2024, ApJ, 962, 24 [NASA ADS] [CrossRef] [Google Scholar]
  109. Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Scholtz, J., Maiolino, R., D’Eugenio, F., et al. 2023, A&A, submitted, [arXiv:2311.18731] [Google Scholar]
  111. Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150 [Google Scholar]
  112. Shapley, A. E., Coil, A. L., Ma, C.-P., & Bundy, K. 2005, ApJ, 635, 1006 [CrossRef] [Google Scholar]
  113. Shapley, A. E., Reddy, N. A., Sanders, R. L., Topping, M. W., & Brammer, G. B. 2023a, ApJ, 950, L1 [NASA ADS] [CrossRef] [Google Scholar]
  114. Shapley, A. E., Sanders, R. L., Reddy, N. A., Topping, M. W., & Brammer, G. B. 2023b, ApJ, 954, 157 [NASA ADS] [CrossRef] [Google Scholar]
  115. Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51 [Google Scholar]
  116. Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 [Google Scholar]
  117. Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165 [Google Scholar]
  118. Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164 [NASA ADS] [CrossRef] [Google Scholar]
  119. Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023a, ApJ, 952, 74 [NASA ADS] [CrossRef] [Google Scholar]
  120. Tacchella, S., Johnson, B. D., Robertson, B. E., et al. 2023b, MNRAS, 522, 6236 [NASA ADS] [CrossRef] [Google Scholar]
  121. Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 [Google Scholar]
  122. Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157 [NASA ADS] [CrossRef] [Google Scholar]
  123. Taylor, A. J., Barger, A. J., & Cowie, L. L. 2022, ApJ, 939, L3 [NASA ADS] [CrossRef] [Google Scholar]
  124. Telford, O. G., Dalcanton, J. J., Skillman, E. D., & Conroy, C. 2016, ApJ, 827, 35 [NASA ADS] [CrossRef] [Google Scholar]
  125. Topping, M. W., Shapley, A. E., Sanders, R. L., et al. 2021, MNRAS, 506, 1237 [CrossRef] [Google Scholar]
  126. Torrey, P., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 484, 5587 [NASA ADS] [Google Scholar]
  127. Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [Google Scholar]
  128. Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110 [NASA ADS] [CrossRef] [Google Scholar]
  129. Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2023, ApJ, 945, 35 [NASA ADS] [CrossRef] [Google Scholar]
  130. Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023, A&A, 677, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  131. Ucci, G., Dayal, P., Hutter, A., et al. 2023, MNRAS, 518, 3557 [Google Scholar]
  132. van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28 [Google Scholar]
  133. Vidal-García, A., Charlot, S., Bruzual, G., & Hubeny, I. 2017, MNRAS, 470, 3532 [Google Scholar]
  134. Vincenzo, F., Matteucci, F., Belfiore, F., & Maiolino, R. 2016, MNRAS, 455, 4183 [Google Scholar]
  135. Williams, H., Kelly, P. L., Chen, W., et al. 2023, Science, 380, 416 [NASA ADS] [CrossRef] [Google Scholar]
  136. Witstok, J., Smit, R., Maiolino, R., et al. 2021, MNRAS, 508, 1686 [NASA ADS] [CrossRef] [Google Scholar]
  137. Wuyts, E., Rigby, J. R., Sharon, K., & Gladders, M. D. 2012, ApJ, 755, 73 [NASA ADS] [CrossRef] [Google Scholar]
  138. Wuyts, E., Kurk, J., Förster Schreiber, N. M., et al. 2014, ApJ, 789, L40 [NASA ADS] [CrossRef] [Google Scholar]
  139. Wuyts, S., Förster Schreiber, N. M., Wisnioski, E., et al. 2016, ApJ, 831, 149 [Google Scholar]
  140. Yabe, K., Ohta, K., Iwamuro, F., et al. 2014, MNRAS, 437, 3647 [NASA ADS] [CrossRef] [Google Scholar]
  141. Yang, H., Malhotra, S., Gronke, M., et al. 2017a, ApJ, 844, 171 [NASA ADS] [CrossRef] [Google Scholar]
  142. Yang, H., Malhotra, S., Rhoads, J. E., & Wang, J. 2017b, ApJ, 847, 38 [NASA ADS] [CrossRef] [Google Scholar]
  143. Yates, R. M., Kauffmann, G., & Guo, Q. 2012, MNRAS, 422, 215 [NASA ADS] [CrossRef] [Google Scholar]
  144. Yates, R. M., Schady, P., Chen, T. W., Schweyer, T., & Wiseman, P. 2020, A&A, 634, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Zahid, H. J., Kewley, L. J., & Bresolin, F. 2011, ApJ, 730, 137 [NASA ADS] [CrossRef] [Google Scholar]
  146. Zahid, H. J., Bresolin, F., Kewley, L. J., Coil, A. L., & Davé, R. 2012, ApJ, 750, 120 [CrossRef] [Google Scholar]
  147. Zahid, H. J., Dima, G. I., Kudritzki, R.-P., et al. 2014, ApJ, 791, 130 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Metallicity scaling relations under different assumptions as to sample selection and derivation of galaxy properties

In the main body of the paper, we discuss the mass–metallicity relation (MZR) and the fundamental metallicity relation (FMR) for a combined sample of galaxies observed with the JWST from different observational programmes. Here, we investigate whether any of the main results presented in the paper change if we adopt different choices from our set of fiducial assumptions, particularly regarding sample selection, SFR measurement, and FMR parametrisation.

In Figure A.1 we show the MZR (left-hand panel) and deviations from the FMR (right-hand panel) for our sample of galaxies which now includes 14 objects from JADES identified as AGN candidates ((Scholtz et al. 2023), highlighted in cyan). For these galaxies, we here assume a mild contribution of AGN ionisation to the emission line spectrum, so that the standard metallicity calibrations adopted for ‘normal’ star-forming galaxies can be applied. Stellar masses and SFRs are also derived assuming negligible AGN contribution to the continuum level; for what concern the SFR however, any contribution to Hα emission from the narrow-line region would imply the inferred SFR have to be considered as upper limits. No significant change to the observed evolution in the metallicity scaling relations is found when including these objects in the sample, with a Spearman correlation ρ = −0.18, p-value=0.017.

In Figure A.2 instead we report, similarly to Figure 7, the deviations from the predictions of the local FMR, but this time assuming different estimates for the SFR of our entire JWST sample. In particular, we adopt the calibration for metal poor galaxies described in Reddy et al. (2022), Shapley et al. (2023b) in the left-hand panel, whereas the calibration from Kennicutt & Evans (2012), and suited for solar metallicities, is assumed in the right-hand panel. As already discussed in Section 3.2 of the main text, once estimated on JADES galaxies the former provides on average lower (s)SFRs by 0.14 dex, whereas the latter delivers higher (s)SFR by 0.23 dex, compared to our fiducial values based on BEAGLE fitting to both PRISM spectra and photometry. This means that the predicted metallicities from the FMR (i.e. Δlog(O/H)) are, on average, higher (lower) than our fiducial case when the low-metallicity (solar metallicity) calibration to convert the Hα flux to SFR is adopted. Nonetheless, the median trends and the observed evolution in the FMR at high-z are robust against the choice of the SFR calibration (Spearman ρ = −0.17, p-value=0.04 in both cases.)

Furthermore, in Figure A.3 we compare the metallicity for our JWST sample with the predictions based on the FMR described as in Andrews & Martini (2013), which is parametrisation in the form O/H versus μα, with μα=0.66. Compared to Curti et al. (2020b), the Andrews & Martini (2013) FMR is based on ‘direct’ (Te-based) metallicities derived from stacked spectra in bins of M and SFR, rather than on individual galaxies with strong-line-based (though Te-calibrated) metallicities. On the one hand, this allows to extend the parametrisation towards lower masses, while on the other, sample selection and weighting effects can possibly bias low the metallicity of stacked spectra dominated by few bright [O III] λ4363 emitters. As a consequence, the inferred SFR-dependence of O/H at (fixed) low M is stronger for Andrews & Martini (2013; with a slope of the linear O/H vs. μα relation of 0.43), which manifests as a lower predicted O/H at low μα, and therefore a smaller offset, for high-z galaxies (see also Nakajima et al. 2023).

Nonetheless, z > 3 galaxies still appear to lie preferentially below the FMR predictions, as shown in the left-hand panel of Figure A.3, with median values in redshift bins offset by -0.07 dex at z∼4 (though now marginally consistent at 2σ with the FMR given the error on the median), -0.21 dex at z∼6 (6σ significance), and -0.39 dex at z∼8 (5σ significance).

In the middle panel of Figure A.3 we show instead the median values in M and redshift bins for the JWST sample on the O/H vs μα projection. Median points at z∼3 − 6 are still offset below the FMR by 0.135 dex, 0.073 dex, and 0.178 dex at log(M)=7.4, 8.1, and 8.9, respectively (at 2.1σ, 1.5σ, and 3.6σ significance, respectively), whereas higher redshift bins (z∼6 − 10) show larger offsets (0.22 dex and 0.25 dex, at 3σ significance) at log(M)>7.5, then only exception being the high-z bin at log(M)=7.35 which is in agreement with the Andrews & Martini (2013) FMR.

Finally, it is important to remark once again that a different choice of the SFR calibration can also affect the predicted FMR-metallicity: this is further shown in the right-hand panel of Fig. A.3, where the same O/H vs μα projection for the Andrews & Martini (2013) FMR is shown when adopting the Hα-based SFR with the solar metallicity calibration from Kennicutt & Evans (2012). We note that adopting the Hα-based SFR and the calibration from Kennicutt & Evans (2012) for all our galaxies, which significantly reduces the number of ionising photons per unit SFR compared to both our fiducial and the Hα-based, low-metallicity calibration scenario, provides a better agreement with the FMR predictions of Andrews & Martini (2013). This is in fact driven by the increase in the inferred SFR, which reduces μα and consequently moves the points towards the left along the x-axis of the diagram, and it observed in particular for median-binned points at 3 <  z <  6 (which are now only mildy offset from the FMR predictions, and consistent with it within ∼2 − −3σ), whereas higher redshift bins (z> 6) at log(M)> 7.5 keep showing an offset below the FMR at > 3σ significance. The choice of the SFR-calibration should therefore be evaluated as carefully as that of the metallicity calibrations, as it might introduce biases and produce tensions among different assessments of the evolution of these important scaling relations (Nakajima et al. 2023).

thumbnail Fig. A.1.

Similar to Figure 5 (left-hand panel) and Figure 7 (right-hand panel), but including the sample of 14 (narrow-line) AGN candidates in JADES-deepScholtz et al. 2023, here highlighted with cyan contours. No clear change in the median trends of the metallicity scaling relations is found, and the main results reported in the paper are unaffected.

thumbnail Fig. A.2.

Same as Figure 7, but with SFR inferred from the (attenuation corrected) Hα flux (or Hβ, at z ≳ 7). The calibration for low metallicity galaxies described in Reddy et al. 2022; Shapley et al. 2023b is adopted for the entire JWST sample in the left-hand panel, while in the right-hand panel the SFR is derived following the calibration from Kennicutt & Evans 2012 for solar metallicity.

thumbnail Fig. A.3.

Comparison between observed metallicity for the JWST sample and the predictions of the FMR, here adopting the formalism from Andrews & Martini 2013. The three panels show, from left to right, the offset from the FMR as a function of redshift, and the measured versus FMR-predicted metallicities based on two different SFR estimates, namely the fiducial value adopted throughout the paper (and based on BEAGLE fitting to both R100 spectra and photometry), and the Hα-based SFR based on the calibration from Kennicutt & Evans 2012, respectively. In the latter case, the SFRs are higher on average by ∼0.23 dex than our fiducial values, and by ∼0.37 dex compared to the SFR-calibration for metal poor galaxies, decreasing μα and improving the agreement with the predictions of the Andrews & Martini 2013 parametrisation of the FMR, especially at z∼4; nonetheless, higher redshift (z > 6) galaxies still sit preferentially below the FMR, appearing metal deficient compared to local galaxies at fixed M and SFR.

Appendix B: Adopted metallicity calibrations

As discussed in Section 3.1, in this work we have adopted the set of strong-line metallicity calibrations presented in Curti et al. (2017) and Curti et al. (2020b), whose polynomial functions have been however slightly revised to better probe the low-metallicity regime. In addition, we have exploited the novel diagnostics R̂ introduced by Laseter et al. (2024). In table B.1 we report the best-fit coefficients for each of the polynomial calibrations adopted in the paper. Each calibrations is parametrised in the form

log(R)= i = 0 n c i x i , $$ \begin{aligned} \text{ log(R)=}\sum _{i=0}^n c_ix^i \ , \end{aligned} $$(B.1)

where x=12+log(O/H)-8.69, ci are the coefficients, and R is the diagnostic line ratio (see Table 1).

Table B.1.

Best-fit polynomial coefficients for the strong-line metallicity calibrations adopted in this paper.

Appendix C: Summary of properties for JADES galaxies

We here report the derived redshift, stellar mass, star-formation rate, and metallicity for the full JADES galaxy sample analysed in this work. In addition, we indicate whether the metallicity has been derived from PRISM or grating spectra, which diagnostics have been used, if NIRCAM photometry has been included in the M and SFR derivation with BEAGLE, and whether a galaxy has been identified as a type-2 AGN candidate (Scholtz et al. 2023).

Table C.1.

Summary of derived properties for JADES galaxies analysed in this paper.

All Tables

Table 1.

Definitions of line ratios adopted throughout the paper.

Table 2.

Median values, error on the median (and standard deviation) in stellar mass, SFR, and metallicity for our sample of galaxies in each of the (M–redshift) bins considered in this work.

Table 3.

Best-fit parameters of the mass–metallicity relation from Eq. (2) for the full sample, and the two redshift bins at z = 3 − 6 and z = 6 − 10, for which we report both the slope βz and the normalisation at M = 108 M, Zm8.

Table B.1.

Best-fit polynomial coefficients for the strong-line metallicity calibrations adopted in this paper.

Table C.1.

Summary of derived properties for JADES galaxies analysed in this paper.

All Figures

thumbnail Fig. 1.

Example of the spectral fitting procedure for one of the JADES galaxies in the sample, namely JADES-GS+53.16268-27.80237 at z = 4.805. Top panel: PRISM (R100) spectrum (in black), showing the region of the main rest-frame optical emission lines. The best-fit to the spectrum from PPXF is shown in red. Bottom panels: medium-resolution grating (R1000) spectrum (and its best-fit) for the same galaxy. From left to right, the panels show a zoom onto the regions of the [O II] and [Ne III], Hβ and [O III], and Hα [N II] and [S II] emission lines, respectively. In all panels, the error spectrum is marked by the cyan shaded region.

In the text
thumbnail Fig. 2.

Comparison between the metallicity derived from medium-resolution gratings and PRISM spectra. We plot only objects for which the requirements described in Sect. 3.1 are satisfied for both configurations. The two distributions scatter across the equality line (in red), with a median offset of 0.01 dex and a standard deviation of 0.15 dex. In 75% of the cases the two measurements are consistent within their 1σ uncertainties.

In the text
thumbnail Fig. 3.

Redshift distribution of the selected sample of 54 JADES galaxies with metallicity measurements analysed in this paper. Their histogram is compared to the distribution of the full ‘JADES-Deep’ sample presented in Sect. 2.1 and Bunker et al. (2023a) with a spectroscopic redshift determination.

In the text
thumbnail Fig. 4.

Distribution in the stellar-mass–SFR plane for our combined JWST sample. This includes both JADES galaxies presented in this work and the literature sample compiled from Nakajima et al. (2023; CEERS), Curti et al. (2023; EROs), and Bunker et al. (2023b; GN-z11), respectively. The top and right-hand inset panels show the histograms of the distribution in both parameters. The SFR for JADES galaxies is derived from BEAGLE fitting to both PRISM spectra and NIRCAM photometry (where available, and to PRISM spectra only otherwise). For CEERS galaxies, the SFRs are compiled from Nakajima et al. (2023) as inferred from the Hβ luminosity and the Kennicutt & Evans (2012) calibration, but have been here scaled down by 0.23 dex to account for the mean offset between the Hα(Hβ)-based and BEAGLE-based SFRs measured for the JADES sample (we refer to Sect. 3.2 for more details). For both ERO galaxies and GN-z11 the M and SFR have been derived via BEAGLE fitting, consistently with what done for JADES sources. The parametrisation of the main sequence of star-formation (SFMS) at z ∼ 0 and z ∼ 6 from Popesso et al. (2023; and their extrapolation at low stellar mass, dashed line) is also shown for reference.

In the text
thumbnail Fig. 5.

Mass-metallicity relation (MZR) for our full JWST sample. Filled circles represent individual galaxies presented in this work from JADES, whereas filled crosses are galaxies observed in the framework of the CEERS programme and compiled from Nakajima et al. (2023), with metallicity derived as detailed in Sect. 3.1. The star symbol reports the JADES/NIRSpec observations of GN-z11 (Bunker et al. 2023b, at z = 10.603), whereas ‘X’ symbols mark galaxies from the EROs as compiled from Curti et al. (2023) and Laseter et al. (2024). Large squared and diamond symbols mark the median values computed in bins of M (full sample, in purple), and (M z), as described in Table 2. An orthogonal linear regression fit to the median values in bins of M for the different redshift sub-samples is shown by the purple (full sample), yellow (z = 3 − 6) and blue (z = 6 − 10) lines, respectively. We include a comparison with previous determinations of the MZR at lower redshifts from Curti et al. (2020b; SDSS at z ∼ 0.07), and Sanders et al. (2021; MOSDEF at z ∼ 2 − 3), as well as the best-fit of the low-mass end of the MZR at z ∼ 3 provided by Li et al. (2023) and based on JWST/NIRISS slitless spectroscopy. The MZR curves at z ∼ 2 − 3 have been scaled down by ∼0.1 dex to account for the systematics differences between the metallicity calibrations used in this work and the Bian et al. (2018) calibrations adopted in the original papers.

In the text
thumbnail Fig. 6.

Comparison with local analogues and cosmological simulations. Left-hand panel: the best-fit MZR presented in this work (for both the full JWST sample in purple, and in the z = 3 − 6 and z = 6 − 10 redshift binsin yellow and blue, respectively) is compared to the Te-based mass-metallicity relation from local, metal-poor ‘Blueberry’ and ‘Green Pea’ galaxies (Yang et al. 2017a,b), which share similar excitation and emission line properties to our high-redshift galaxy sample (Cameron et al. 2023b). Their distribution agrees well in slope and normalisation (with ∼0.11 dex offset in log(O/H)) with the average relation inferred in this work for our low-mass JWST sample. Right-hand panel: our best-fit MZRs are compared with the predictions from different suites of cosmological simulations at z ∼ 6–8, namely FIRE (Ma et al. 2016), IllustrisTNG (Torrey et al. 2019), FirstLight (Langan et al. 2020), Astraeus (Ucci et al. 2023), SERRA (Pallottini et al. 2022), and with the chemodynamical simulations from Kobayashi & Taylor (2023). In addition, the typical MZR slopes predicted by (semi-) analytical chemical evolution models (e.g. Davé et al. 2012) implementing feedback via either ‘energy-driven’ or ‘momentum-driven’ winds are also shown (the normalisation is assumed arbitrary to aid visualisation). The latter agrees well with the slope observed at low stellar mass for galaxies at z ≥ 3.

In the text
thumbnail Fig. 7.

Deviations of the JWST sample from the predictions of the local FMR plotted as a function of redshift. Marker symbols are as in Fig. 5, with each galaxy colour-coded by its sSFR. The weighted average of the deviation in three redshift bins (colour-coded by the mean sSFR of the binned sample) for both the JWST sample and for galaxy samples at lower redshift (SDSS at z ∼ 0, KLEVER at z ∼ 1.5 − 2.5, MOSDEF at z ∼ 3) are marked by larger symbols. Dashed error bars report the dispersion in Δlog(O/H) of each binned sample, whereas thick, solid bars report the standard error on the mean (which is often smaller than the marker’s size). Overall, z >  3 galaxies observed by JWST are scattered across the local FMR predictions, but show evidence of being preferentially offset towards lower metallicity values than expected with increasing redshift, with an average offset of −0.4 dex at z >  6.

In the text
thumbnail Fig. 8.

FMR offset as a function of stellar mass and the projection of minumum scatter. Left-hand panel: deviations from the FMR plotted as a function of stellar mass. Symbols are as in Figs. 5 and 7, with individual points (and M–redshift bins) colour-coded according to redshift. At fixed stellar mass, JWST galaxies are preferentially offset towards lower O/H with incresing redshift, and a similar trend is also seen with increasing M in the highest redshift bin. Right-hand panel: two-dimensional projection of the M-SFR-O/H space in the O/H vs μ = log(M)-αlog(SFR) plane, assuming the parametrisation for high sSFR (log(sSFR/yr−1) > −9.5) galaxies inferred from the analysis of the SDSS sample (α = 0.65, Curti et al. 2020b). Most of the individual galaxies, as well as the median values in the M-redshift bins, are offset below the relation, further suggesting that the local FMR parametrisation does not hold at the very high-z and low-M probed in this work.

In the text
thumbnail Fig. A.1.

Similar to Figure 5 (left-hand panel) and Figure 7 (right-hand panel), but including the sample of 14 (narrow-line) AGN candidates in JADES-deepScholtz et al. 2023, here highlighted with cyan contours. No clear change in the median trends of the metallicity scaling relations is found, and the main results reported in the paper are unaffected.

In the text
thumbnail Fig. A.2.

Same as Figure 7, but with SFR inferred from the (attenuation corrected) Hα flux (or Hβ, at z ≳ 7). The calibration for low metallicity galaxies described in Reddy et al. 2022; Shapley et al. 2023b is adopted for the entire JWST sample in the left-hand panel, while in the right-hand panel the SFR is derived following the calibration from Kennicutt & Evans 2012 for solar metallicity.

In the text
thumbnail Fig. A.3.

Comparison between observed metallicity for the JWST sample and the predictions of the FMR, here adopting the formalism from Andrews & Martini 2013. The three panels show, from left to right, the offset from the FMR as a function of redshift, and the measured versus FMR-predicted metallicities based on two different SFR estimates, namely the fiducial value adopted throughout the paper (and based on BEAGLE fitting to both R100 spectra and photometry), and the Hα-based SFR based on the calibration from Kennicutt & Evans 2012, respectively. In the latter case, the SFRs are higher on average by ∼0.23 dex than our fiducial values, and by ∼0.37 dex compared to the SFR-calibration for metal poor galaxies, decreasing μα and improving the agreement with the predictions of the Andrews & Martini 2013 parametrisation of the FMR, especially at z∼4; nonetheless, higher redshift (z > 6) galaxies still sit preferentially below the FMR, appearing metal deficient compared to local galaxies at fixed M and SFR.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.