We analysed the gasphase metallicity properties of a sample of lowstellarmass log MSUBSUBMSUBSUB  9 galaxies at 3 lt z lt 10 observed with JWSTNIRSpec as part of the JADES programme in its deep GOODSS tier. By combining this sample with more massive galaxies at similar redshifts from other programmes we study the scaling relations between stellar mass MSUBSUB oxygen abundance OH and starformation rate SFR for 146 galaxies spanning three orders of magnitude in stellar mass and out to the epoch of early galaxy assembly. We find evidence for a shallower slope at the lowmass end of the massmetallicity relation MZR with 12logOH  7.72  0.02  0.17  0.03 logMSUBSUB10SUP8SUP MSUBSUB in good agreement with the MZR probed by local analogues of highredshift systems such as the Green Pea and Blueberry galaxies. The inferred slope is well matched by models including momentumdriven supernova SN winds suggesting that feedback mechanisms in dwarf galaxies and at high z might be different from those in place at higher masses. The evolution in the normalisation is observed to be relatively mild compared to previous determinations of the MZR at z  3 0.1  0.2 dex across the explored mass regime. We observe a deviation from the local fundamental metallicity relation FMR for our sample at high redshift especially at z gt 6 with galaxies significantly less enriched than predicted given their MSUBSUB and SFR with a median offset in logOH of 0.5 dex significant at 5. These observations are consistent with an enhanced stochasticity in the gas accretion and starformation history of highredshift systems prompting us to reconsider the nature of the relationship between MSUBSUB OH and SFR in the early Universe. P The highlevel data products exploited in this paper for the deep tier of the JADESGTO Programme in GOODSSouth namely the redshifts stellar masses starformation rates and metallicities are reported in Table C.1. Fully reduced spectra and emission line fluxes are available through the MAST database at A hrefhttpsarchive.stsci.eduhlspjadeshttpsarchive.stsci.eduhlspjadesA and are described in ltxref reftypebibr ridR16gtBunker et al. 2023altxrefgt.
['ISM: abundances', 'Galaxy: evolution', 'Galaxy: fundamental parameters', 'galaxies: abundances', 'galaxies: high-redshift', 'galaxies: ISM', 'Astrophysics - Astrophysics of Galaxies']
JADES Insights into the lowmass end of the massmetallicitySFR relation at 3 lt z lt 10 from deep JWSTNIRSpec spectroscopy
{'JADES: Insights on the low-mass end of the mass - metallicity star-formation rate relation at 3 < z < 10 from deep JWST/NIRSpec spectroscopy': "Mirko Curti 1,2,3 ⋆ , Roberto Maiolino 2,3 , Emma Curtis-Lake 4 , Jacopo Chevallard 5 , Stefano Carniani 6 , Francesco D'Eugenio 2,3 , Tobias J. Looser 2,3 , Jan Scholtz 2,3 , Stephane Charlot 7 , Alex Cameron 5 , Hannah Übler 2,3 , Joris Witstok 2,3 , Kristian Boyett 8,9 , Isaac Laseter 10 , Lester Sandles 2,3 , Santiago Arribas 11 , Andrew Bunker 5 , Giovanna Giardino 12 , Michael V. Maseda 10 , Tim Rawle 13 , Bruno Rodríguez Del Pino 11 , Renske Smit 14 , Chris J. Willott 15 , Daniel J. Eisenstein 16 , Ryan Hausen 17 , Benjamin Johnson 16 , Marcia Rieke 18 , Brant Robertson 19 , Sandro Tacchella 2,3 , Christina C. Williams 20 , Christopher Willmer 18 , William M. Baker 2,3 , Rachana Bhatawdekar 13,21 , Eiichi Egami 18 , Jakob M. Helton 18 , Zhiyuan Ji 18 , Nimisha Kumari 22 , Michele Perna 11 , Irene Shivaei 18 , and Fengwu Sun 18  \n(A ffi liations can be found after the references)", 'ABSTRACT': "We analyse the gas-phase metallicity properties of a sample of low stellar mass (log M ⋆ / M ⊙ ≲ 9) galaxies at 3 < z < 10, observed with JWST / NIRSpec as part of the JADES programme in its deep GOODS-S tier. By combining this sample with more massive galaxies at similar redshifts from other programmes, we study the scaling relations between stellar mass (M ⋆ ), oxygen abundance (O / H), and star-formation rate (SFR) for 146 galaxies, spanning across three orders of magnitude in stellar mass and out to the epoch of early galaxy assembly. We find evidence for a shallower slope at the low-mass-end of the mass-metallicity relation (MZR), with 12 + log(O / H) = (7.72 ± 0.02) + (0.17 ± 0.03) log(M ⋆ / 10 8 M ⊙ ), in good agreement with the MZR probed by local analogues of high-redshift systems like 'Green Pea' and 'Blueberry' galaxies. The inferred slope is well matched by models including 'momentum-driven' SNe winds, suggesting that feedback mechanisms in dwarf galaxies (and at high-z) might be di ff erent from those in place at higher masses. The evolution in the normalisation is observed to be relatively mild compared to previous determinations of the MZR at z ∼ 3 ( ∼ 0 . 1 -0 . 2 dex across the explored mass regime). We observe a deviation from the local fundamental metallicity relation (FMR) for our sample at high redshift, especially at z > 6, with galaxies significantly less enriched (with a median o ff set in log(O / H) of ∼ 0 . 5 dex, significant at ∼ 5 σ ) than predicted given their M ⋆ and SFR. These observations are consistent with an enhanced stochasticity in the star-formation history, and / or with an increased e ffi ciency in metal removals by outflows, prompting us to reconsider the nature of the relationship between M ⋆ , O / H, and SFR in the early Universe.  \nKey words. Galaxies: abundances - Galaxies: ISM - ISM: abundances - Galaxies: evolution - Galaxies: high redshift", '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 the spectroscopic characterisation of the properties of the interstellar medium (ISM) and the nature of ionizing spectra of primordial galaxies. In this context, the study of the gas-phase metallicity 1 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. 2019;  \nBaker & Maiolino 2023) at high redshift is critical 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 di ff erent processes, the metallicity of galaxies has been observed to correlate also with di ff erent galactic properties. In particular, an anti-correlation between metallicity and star-formation rate, 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. 2023b).  \nSeveral 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. 2022, Finkelstein et al., in preparation), and JWST Advanced Deep Extragalactic Survey (JADES, Pro-  \nposal ID: 1210; Robertson et al. 2022; Curtis-Lake et al. 2022, Eisenstein et al., in preparation). Recent works have indeed investigated the ionisation properties of galaxies beyond z = 3, showing that these sources exhibit emission line ratios consistent with hard ionizing spectra and low metallicities (Mascia et al. 2023; Matthee et al. 2022; Sanders et al. 2023b; Cameron et al. 2023b). These kind of analyses have been pushed to some of the highest redshift galaxies discovered, with Williams et al. (2022) 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 reporting 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 hence be exploited to assess the chemical enrichment of their ISM.  \nIn the past, several studies have investigated the evolution in the metallicity scaling relations out to z ∼ 3, finding signatures of clear trend of decresing 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; Henry et al. 2013a; Yabe et al. 2014; Zahid 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 as driven by the interplay of gas flows, starformation, and metal enrichment, which is reflected on such 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).  \nMore 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. 2022; Rhoads et al. 2023; Langeroodi et al. 2022; Matthee et al. 2022; Heintz et al. 2022; Nakajima et al. 2023; Shapley et al. 2023a). Despite the e ff orts, the majority of previous 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) until the arrival of JWST (Li et al. 2022).  \nMoreover, it is worth to recall that the vast majority of such metallicity determinations are based on strong emission line ratios that have been calibrated on 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 di ff erent ionization 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; Sanders et al. 2023b; Cameron et al. 2023b). Nonetheless, the rate of auroral lines detections at high-z is now rapidly increasing thanks to the JWST , and novel attempts to re-calibrate the classical strong line diagnostics for the high-redshift Universe are being provided (Sanders et al. 2023a).  \nIn this paper, we aim to explore further the cosmic evolution of the metallicity scaling relations, leveraging deep spectroscopy  \nTable 1. Definitions of line ratios adopted throughout the paper.  \nwith JWST / NIRSpec to probe the low stellar mass regime M ⋆ ≈ 10 6 . 5 -10 9 . 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 ⋆ ∼ 10 6 . 5 , with existing datasets, we assess the evolution of the MZR slope and normalisation, and its implications on the physical processes in place in early galaxies. Moreover, we test the FMR framework down to the epoch when the Universe was ≲ 500 Myr old.  \nThe paper is organised as follows. In Section 2, we describe observations, data reduction, and spectral fitting procedures. Section 3 we present the analysis performed to derive the physical quantities of interest, i.e., M ⋆ , SFR, and metallicity. In Section 4, we assess the mass-metallicity relation (MZR) at z > 3 combining our JADES sample with galaxies at higher masses drawn from the CEERS programme, whereas in Section 5 we investigate any evolution in the framework of the fundamental metallicity relation, and discuss possible inferences about the interplay between gas flows and star-formation regulating early galaxy assembly. In Section 6, we summarise our conclusions.  \nThroughout this paper, we adopt a Planck Collaboration et al. (2020) cosmology, a Chabrier (2003) initial mass function, and a solar metallicity 12 + log(O / H) = 8.69 (Allende Prieto et al. 2001; Asplund et al. 2009).", '2.1. Observations': "The data presented in this paper have been obtained via multiobject spectroscopy observations conducted with the microshutter 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 the 21st-25th of October 2022 (Program ID: 1210; PI: N. Luetzgendorf) in the Great Observatories Origins Deep Survey South (GOODSS) 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 / CLEARconfiguration (hereinafter just 'PRISM'), and 8,403 s integration in each of G140M / F070LP, G235M / F170LP, G395M / F290LP (hereinafter 'medium resolution gratings'). Across three visits, this totals 28 hours of integration in the PRISM, providing continuous spectral coverage from 0 . 6 -5 . 3 µ m at R ∼ 30 -300, and ∼ 7 hours in each of the medium resolution gratings, providing R ∼ 1000 across the full spectral range of NIRSpec 2 .  \nFig. 1. Top panel: Example fit to the PRISM (R100) spectrum (in black) of the z = 4 . 805 galaxy JADES-GS + 53.16268-27.80237, zoomed-in on the region of the main rest-frame optical emission lines. The best-fit 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-in on the regions of the [O ii ] and [Ne iii ], H β and [O iii ], H α [N ii ] and [S ii ] emission lines, respectively. In all panels, the error spectrum is marked by the cyan shaded region.  \n<!-- image -->  \nObservations within each visit were performed adopting a 3-shutter nodding pattern, with the central pointing of each visit dithered by < 1 arcsec to sample di ff erent 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 onethird of the total exposure time. Each pointing had a bespoke MSA configuration, and target allocation was performed via the eMPT 3 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.  \nWe 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 o ffi cial 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 pre-processed count-rate images were then processed using a data reduction  \npipeline 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 preparation). 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 e ffi ciently. 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 intrashutter position of each source by leveraging on a full forwardmodelling 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, 2022), which models the continuum as a linear super-  \nposition of simple stellar-population (SSP) spectra, using nonnegative 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 C3Ktheoretical atmospheres (Conroy et al. 2019). The flux blueward of the Lyman break was manually set to zero. These templates are complemented by a 5 th -degree multiplicative Legendre polynomial, to take into account systematic di ff erences between the SSPs and the data (e.g., 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 (e.g., 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).  \nWe show an example fit for a z = 4 . 805 galaxy in Figure 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.', '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 × 10 4 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 preparation), which incorporates a two-component dust attenuation model that accounts for the differential attenuation of nebular and stellar emission in a selfconsistent way.  \nThe gas-phase metallicity is derived exploiting a revisited version of the calibrations presented in Curti et al. (2020b) (which are anchored to the T e abundance scale), refined using  \nFig. 2. Comparison between the metallicity derived from medium resolution gratings and PRISM spectra, for objects in which both configurations satisfy the requirements described in Section 3.1. The two distributions scatter across the equality line (in red), with a median o ff set 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.  \n<!-- image -->  \nFig. 3. The redshift distribution of the selected sample of 62 JADES galaxies with metallicity measurements analysed in this paper, compared to the distribution of the total 'JADES-Deep' sample presented in Section 2.1 and Bunker et al. 2023a with spectroscopic redshift determination.  \n<!-- image -->  \na sample of local, metal-poor galaxies to better sample the lowmetallicity 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. (2023). 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.  \nWe visually inspected the spectra to remove sources a ff ected 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. (2023a). 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., in preparation) 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; Kau ff mann et al. 2003) , whose applicability in discriminating among di ff erent ionising sources at high redshift has been in fact recently questioned, especially at low metallicity (Übler et al. 2023; Maiolino et al. 2023a). 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 sub-sample of galaxies.  \nIn 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  \nlog(L) ∝ X i (Robs , i -R cal , i ) 2 ( σ 2 obs , i + σ 2 cal , i ) , (1)  \nwhere 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.  \nFor 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 e ff ectively 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  \n(being almost flat and hence 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 o ff set between the low- and highmetallicity 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.  \nThe 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 a ff ected by [N ii ] contamination of the H α line. As for 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:  \n- -we take the O / H inferred from medium resolution (R1000) gratings if at least two diagnostics were available, else\n- -we take the O / H inferred from the PRISM if at least two diagnostics were available in the PRISM spectrum, else\n- -we take the O / H inferred from medium resolution gratings as based on one individual diagnostics (typically R3), exploiting the information from 3-sigma upper / lower limits on other line ratios (e.g., N2, O32) to break the degeneracy of double-valued calibrations, else\n- -we take the O / H inferred from the PRISM as based on one diagnostics (typically R3) in combination with the information from the upper limits on [O ii ] λλ 3727,29.  \nThe 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 intercalibration uncertainties between the two configurations that are not yet fully understood. However, we compare the metallicity derived from gratings and prism spectra for galaxies where both measurements are available in the right-hand panel of Figure 2, i.e., 35 galaxies. The two measurements are in agreement within the 1 σ individual uncertainties in ∼ 75% of the cases (and in ∼ 95% within 2 σ ), with negligible systematic o ff set (i.e., 0 . 01 dex median o ff set towards higher grating metallicities) and a standard deviation of 0 . 15 dex.  \nIn 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.218.54, with an average of 7.77 (i.e., 0 . 12 Z ⊙ ). In Figure 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 redshift 4 , which in turn constitutes ∼ 70 per cent of the total sample introduced in Section 2.1; for a more detailed description of the full JADES deep galaxy sample, the selection function, and redshifts spectroscopic confirmation, we refer to Bunker et al. (2023a). In particular, the final sub-sample selected for the current analysis represents the ∼ 57 percent of  \nthe 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. 2022; Looser et al. 2023).", '3.2. Stellar masses and star-formation rates': '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 & the JADES Collaboration 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 100M ⊙ , 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 preparation) for further information of the beagle fitting procedure.  \nThe star-formation rates of our galaxies are estimated in two di ff erent ways, i.e., from the output of the beagle SED fitting run described above, which provides an estimate of the starformation rate averaged over the past ten Myr, and from the attenuation-corrected 5 H α luminosity, following the recipe 6 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 re-scaled by the theoretical case B recombination factor at T ∼ 1 . 5 × 10 4 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 Section 3.1).  \nOn 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 di ff erent underlying stellar population synthesis models (i.e., BPASS for Reddy et al. 2022, revised Bruzual & Charlot 2003 for beagle ), ii) di ff erent 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  \nfor 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 (Section 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 ⊙ ≈ 10 6 . 5 -10 9 ), 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.  \nTo preserve a good level of self-consistency in the analysis, we have re-computed 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 Section 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 o ff set (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 o ff sets between the SFRs compiled for CEERS galaxies and the fiducial values assumed for JADES objects.  \nStellar masses in the CEERS sample are based on fits to NIRCAMphotometry performed with the prospector code (Johnson et al. 2021). A full assessment of the systematics involved in the stellar mass determination as based on di ff erent 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 di ff erential measurement uncertainties, we do not apply inverse variance weighting based on M ⋆ uncertainties in any of our fitting analyses in the paper.  \nIn addition, we include in the analysis four targets whose metallicity have been derived with the 'direct' T e-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. (2023), though the nature of its ionising source is debated and the presence of an AGN has been proposed (Maiolino et al. 2023b). This brings the total number of galaxies analysed in this paper to 146.  \nFor the three ERO obejcts, M ⋆ and SFR are derived via bea -gle fitting to NIRCAM photometry as described in Curti et al. (2023), whereas the metallicities have been re-computed following the same data reduction implemented for the JADES galaxies analysed in this paper, as also described in Laseter et al. (2023). 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-o ff of 100M ⊙ , 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  \nFig. 4. The distribution in the stellar mass versus star-formation rate plane for our combined JWST sample, including 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 o ff set between the H α (H β )-based and beagle -based SFRs measured for the JADES sample (we refer to Section 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 (and its extrapolation at low-mass)  \n<!-- image -->  \n⋆  \n¯  \nof the main sequence of star-formation (SFMS) at z ∼ 0 and z ∼ 6 from Popesso et al. 2023 is also shown for reference.  \net 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 GNz11) 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 Section 3.1: this approach provides 12 + log(O / H) = 7 . 67 + 0 . 16 -0 . 10 . At the same time, using the same indirect estimate of the [O iii ] λ 5007 flux to apply the T e 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 a ff ecting 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.  \nIn Figure 4 we show 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 starformation rate 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 ⋆ ∼ 10 9 . 5 M ⊙ ( ⟨ log(M ⋆ / M ⊙ ) ⟩ CEERS = 8.72; ⟨ log(M ⋆ / M ⊙ ) ⟩ JADES  \n= 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 comparison, we plot the parametrisation of the 'Main Sequence of Star Formation' (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 vast majority of our combined, low-mass JWST galaxy sample is o ff set above the low-mass extrapolation of the relation at z ∼ 6, given our fiducial M ⋆ and SFR measurements.", '4. The evolution of the mass-metallicity relation beyond z=3': "In Figure 5, we show our combined JWST galaxy sample on the mass-metallicity plane. Each point is marked with a di ff erent 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 ⟨ z ⟩ 3 -6 = 4 . 76 and that of the z ∈ 6 -10 sub-sample is ⟨ z ⟩ 6 -10 = 6 . 73. In particular, filled circle points in the plot report galaxies from the JADES survey, while filled crosses symbols mark galaxies from  \nFig. 5. The 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 Section 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. 2023. 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 di ff erent 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. 2022 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 di ff erences between the metallicity calibrations used in this work and the Bian et al. 2018 calibrations adopted in the original papers.  \n<!-- image -->  \n¯  \nTable 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. The number of galaxies in each bin and their median redshift is also reported.  \nCEERS (Nakajima et al. 2023), whose metallicity has been remeasured as described in Section 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 di ff erent bins in stellar mass. The size of the stellar mass bins is not uniform to maintain a reasonable (i.e., at least 10) num-  \nber of galaxies in each bin. The average properties of the binned samples are reported in Table 2.  \nOur galaxies present a median o ff set 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  \nFig. 6. 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 T e-based mass-metallicity relation from local, metal-poor 'Blueberry' and 'Green Pea' galaxies (Yang et al. 2017b,a), 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 o ff set 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 di ff erent 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 'energydriven' 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.  \n<!-- image -->  \n¯  \n<!-- image -->  \n¯  \nTable 3. Best-fit parameters of the mass-metallicity relation from equation 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 ⋆ = 10 8 M ⊙ , Zm8.  \nrealisations 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 di ff erent metallicity calibrations (the calibrations from Bian et al. (2018) are adopted in Sanders et al. 2021), we have re-computed the metallicity for the MOSDEF stacked spectra with the same methodology outlined in Section 3.1, finding an average o ff set of ∼ 0 . 088 dex towards lower metallicities. Therefore, in Figure 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 o ff set 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.  \nIn 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  \nwithin the two redshift bins separately, in the functional form  \n12 + log(O / H) = β z log(M ⋆/ 10 8 M ⊙ ) + Zm8 (2)  \nwhere β z is the slope at a given redshift interval, and Zm8 is the normalisation at log(M ⋆ / M ⊙ ) = 8. The results are displayed in Figure 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.  \nWe find indications of a flattening of the slope of the lowmass 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 e ffi ciency of stellar winds scales with M ⋆ , governs the MZR slope over the last ∼ 12Gyr 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 hence limited in the stellar mass range they could probe, struggling to access the regime below log(M ⋆ ) < 10 9 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. (2022), and report similar evidence for a shallower slope ( β 2-3 = 0 . 17 ± 0 . 03, as shown in magenta in Figure 5) compared to the simple extrapolation  \nof 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. (2022). 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.  \nFinally, it is interesting to note that the slope inferred at z = 6 -10 (in blue in Figure 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 probably subject to strong selection biases, hence it is di ffi cult to draw any strong conclusion at this stage relative to any possible evolution in the MZR slope between z = 3 -6 and z > 6.", '4.1. Comparison with 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, e.g., the 'Green Pea' and 'Blueberry' galaxies compiled from Yang et al. (2017a,b), whose metallicity has been derived with the 'direct' T e 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).", '4.2. Comparison with models and simulations': "We here compare the observed MZR in our z = 3 -10 galaxy sample with the predictions of di ff erent cosmological simulations, as depicted in the right-hand panel of Figure 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 SFRweighted 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  \nUniverse. Results from the FirstLight simulations at z = 6 are compiled from (Langan et al. 2020, , and Langan, private communication), 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).  \nOverall, most simulations appear to reasonably match the normalisation of the observed MZR at M ⋆ ≈ 10 8 -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 massmetallicity plane below M ⋆ = 10 8 M ⊙ , 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).  \nWe also compare our inferred MZR slopes with those predicted by (semi-)analytical chemical evolution models under the assumption of di ff erent 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 'momentumdriven 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 Figure 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 di ff erent mass regimes, producing a shallower MZR slope at low M ⋆ . At the same time, observing a similar MZR slope at M ⋆ ∼ 10 7 -10 9 between z ∼ 3 (from Li et al. 2022) and z ∼ 6 (from this work) suggests that the scaling of the massloading factor of outflows with stellar mass is invariant in this mass regime during this epoch.  \nAnother 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 evolve 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; Tacconi et al. 2020). If the evolution proceeds at a similar rate even at earlier epochs, the impact of metal dilution from larger gas reservoirs is expected to become more important as redshift increases. However, this seems to be in tension with the relatively mild evolution in the MZR normalisation observed in our JWST sample.  \nIn the framework of analytical chemical evolution models that allow the gas fraction to vary (see e.g., Peeples & Shankar 2011; Lilly et al. 2013; Sanders et al. 2021), and under the assumption of (nearly-)pristine gas inflows, the metallicity of the  \nISM is the result of the balance between the nucleosynthetic stellar yield, the gas dilution, and the amount of metals lost due to outflows. If we assume that the same scaling of the gas fraction µ gas (here defined as µ gas = Mgas / M ⋆ ) with redshift as inferred up to z ∼ 3 persists up to z ∼ 6 (e.g., log( µ gas) ∝ 2 . 49 log(1 + z) as based on the linear parametrisation of Tacconi et al. 2018), and fixing all the other parameters such as stellar yield, return fraction, and mass-loading factor of outflows, adopting the formalism of the Peeples & Shankar (2011) model (their equations (9) and (10)) one can predict the evolution in the metallicity at fixed M ⋆ in the low-mass regime (e.g., at M ⋆ ∼ 10 8 M ⊙ ) as solely driven by the evolution in the gas fraction with redshift.  \nWe perform these calculations for two di ff erent scenarios. First, we assume the low-mass extrapolation of the MZR at z ∼ 3 . 3 from Sanders et al. (2021) (with a slope of 0.30) to solve for the mass-loading factor of outflows ζ w at M ⋆ = 10 8 M ⊙ ; then, assuming a constant ζ w with redshift between z = 3 and z = 6, we derive the expected metallicity at M ⋆ = 10 8 M ⊙ by evolving the gas fraction µ gas with redshift according to the scaling relations of Tacconi et al. (2018). From this, we predict an evolution in metallicity of 0 . 26 dex, larger than the almost negligible evolution probed by the extrapolation of the Sanders et al. (2021) MZRat 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 ⋆ = 10 8 M ⊙ ) of the MZR at z = 3 in the dwarf regime from Li et al. (2022); fixing ζ w with redshift, we predict a decrease in log(O / H) as driven by the evolution in the gas fraction of 0 . 35 dex, which is in slightly better agreement but still larger than observed in this mass regime (i.e, 0 . 25 dex) between z = 3 and z = 6. In both cases, we adopt the oxygen yield and return fraction provided by Vincenzo et al. (2016) for a Chabrier (2003) IMF with an upper cuto ff of 100M ⊙ , i.e. y = 0 . 030 and R = 0 . 30, respectively. Nonetheless, we note that in such analytical chemical evolution models the predicted normalisation of the MZRishighly sensitive to the large uncertainties on the assumed of oxygen stellar yield and the shape of the upper mass cuto ff of the IMF (Vincenzo et al. 2016), the relative contribution of atomic and molecular gas to the total gas mass at high redshift, which is generally estimated from relatively high-mass galaxy samples (Saintonge et al. 2017; Catinella et al. 2018), the relative importance of star-formation e ffi ciency (Baker et al. 2023a), and the fraction of metals expelled from outflows, which may then be re-accreted onto galaxies via galactic fountains from the enriched CGM (Fraternali & Binney 2008; Péroux et al. 2020), whereas on the other hand systematics in the adopted metallicity calibrations can shift the normalisation of the 'observed' MZR by non-negligible amounts.  \nFinally, 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. We note the sample size considered in the ected by selection biases) to perform a robust assessment of the scatter of the scaling relation, however we can tentatively measure the amount of intrinsic scatter in our galaxy sample as  \nobs is the standard deviation meas is the average measurement uncertainty associated with the metal073 dex for our full sample, consistent with the scatter measured in the local MZR at 0.075 Curti et al. 2020b). We note that the scatter at high masses in local galaxies  \npresent work not being large enough (as well as being a ff σ MZR = q σ 2 obs -σ 2 meas , where σ of the full JWST sample around the best-fit relation and σ licity determination 7 . We find σ MZR ∼ 0 . higher masses (log(M ⋆ / M ⊙ ) > 9.5) ( σ MZR , z = 0 =  \nis 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.", '5. Evidence for evolution in the fundamental metallicity relation ?': "It is known that the observed distribution of galaxies on the mass-metallicity plane reflects selection e ff ects 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. 2018). 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 it is usually referred to as the Fundamental Metallicity Relation (FMR, Mannucci et al. 2010), and which takes observationally into account the impact of the current star formation rate on the expected level of chemical enrichment of the ISM, as driven by the interplay between gas accretion, star formation, and outflows.  \nIn Figure 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 Figure 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 o ff set 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. 2017b,a; 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 di ff erent symbols as outlined in the legend, and colour-coded by the average sSFR value of the relative sample.  \nIndividual galaxies are mostly o ff set 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 KolmogorovSmirnov (KS) test on the two distributions, i.e. 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 di ff erent (p-value << 0 . 01), both for the global sample and within all individual redshift bins.  \nLooking 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  \nFig. 7. Deviations of the JWST sample from the predictions of the local Fundamental Metallicity Relation (FMR), here plotted as a function of redshift. Marker symbols are as in Figure 5, with each galaxy colour-coded by its specific star-formation rate (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 o ff set towards lower metallicity values than expected with increasing redshift, with an average o ff set of -0 . 4 dex at z > 6.  \n<!-- image -->  \nof 0 . 38 dex). A Z-test confirms that these are all significantly o ff set 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 Figure 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.  \nFurthermore, 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).  \nFinally, in the left-hand panel of Figure 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 di ff erent levels of o ff set from the FMR regardless of their stellar mass, from the median values we find weak evidence for higher redshift galaxies to be preferentially o ff set from the FMR at fixed stellar mass, especially in the two highest M ⋆ bins, where the o ff set between median points at di ff erent redshift is significant at 2 . 3 σ and 1 . 2 σ , respectively. No clear trend is instead identified with M ⋆ , at fixed redshift.", '5.1. Interpreting the apparent FMR evolution': "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  \nFig. 8. Left-hand panel: Deviations from the FMR plotted as a function of stellar mass. Symbols are as in Figure 5 and 7, with individual points (and M ⋆ -redshift bins) colour-coded by redshift. At fixed stellar mass, JWST galaxies are preferentially o ff set 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 o ff set 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.  \n<!-- image -->  \n¯  \n<!-- image -->  \n¡  \ncould 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 sometimes modelled with a scatter term in the baryonic accretion rate parameter (e.g., Forbes et al. 2014), the longer timescales over which the enrichment from supernovae can balance accretion and star-formation would make a system more likely to be observed out of the equilibrium which is at the basis of the scaling relations observed at lower redshift. (Davé et al. 2011; Zahid et al. 2012). On the other hand, the e ffi ciency of metal removals could also be enhanced, such that the metal loading factor of outflows is higher at fixed SFR than the one required to reproduce the observed FMR at lower redshifts (Baker et al. 2023b).  \nA possible way to di ff erentiate between these di ff erent scenarios is to study the morphology of galaxies with the position relative to the FMR. For instance, Tacchella et al. (2022) 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 o ff -set, 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.  \nAs seen from a di ff erent 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,  \nthe 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 di ff erent mass-metallicity relations are invariant for samples with di ff erent 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.  \n⋆  \nIn Curti et al. (2020b), such a possibility was briefly discussed, and two di ff erent regimes were identified and explored, providing two di ff erent 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 SFRdependence is significantly di ff erent in the two regimes, being much more prominent for high sSFR galaxies, which translates into two di ff erent µα projections, with α = 0 . 22 at low sSFR, and α = 0 . 65 at high sSFR 8 . However we stress that, being based on SDSSgalaxies, 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 Figure 8: such relation, though calibrated on the highest sSFR galaxies in the SDSS, is  \nnot 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).", '6. Summary and Conclusions': 'We have analysed the metallicity properties of a sample of lowmass (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 gratings 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 epoch of early galaxy assembly. Our main findings are summarised as follows:  \n- -We find evidence for a mild evolution in the normalisation of the mass-metallicity relation (MZR) at z > 3, as previously suggested by similar studies (Langeroodi et al. 2022; Heintz et al. 2022; Matthee et al. 2022; 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 show an average o ff set 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).\n- -An orthogonal distance regression fit performed on the full JWST sample suggests a flattening of the slope of low-mass end of the MZR at z = 3 -10 (Figure 5). This is in agreement with previous findings at z ∼ 3 in the same stellar mass regime as reported by Li et al. (2022), compared to which wemeasure however 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 Figure 6).\n- -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 Figure 6). On the other hand, analytical chemical evolution models (e.g. Peeples & Shankar 2011) in which the gas fraction is left free to vary according to the same redshift scaling as inferred up to z ∼ 3 (Tacconi et al. 2018) would predict (once fixing all the other parameters) a stronger evolution in the MZR normalisation at low M than observed.\n- ⋆ -Notably, our galaxy sample is distributed over a similar region (and with a fully consistent inferred slope) of the massmetallicity plane as low-redshift, low M ⋆ , metal poor systems with extreme ionisation conditions like \'Blueberry\' and \'Green Pea\' galaxies (Yang et al. 2017b,a) (left panel of Figure 6). These objects have been identified as local analogues of high-redshift galaxies, and we observe that they  \nshare similar metallicity properties to those pertaining to the high-z sources analysed in this paper (see also Cameron et al. 2023b).  \n- -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 redshift increases (Figure 7). This suggests that in these objects metals are either more e ffi ciently 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.\n- -This behaviour also highlights potential inconsistencies of 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 starformation (right panel of Figure 8).  \nThe observations presented in this paper represent a step to 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 to further reduce the large systematic uncertainties and provide a more comprehensive picture of galaxy formation across di ff erent cosmic epochs, a goal fully within the reach of JWST capabilities.  \nAcknowledgements. 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 / CSAJames 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.', 'Data Availability': 'The high-level data products exploited in this paper for the deep tier of the JADES-GTO Programme in GOODS-South, i.e., the redshifts, stellar masses, star-formation rates, and metallicities, are reported in Table B.1. Fully reduced spectra and emission line fluxes are available through the MAST database at https://archive.stsci.edu/hlsp/jades , and are described in Bunker et al. (2023a).', 'References': "Allende Prieto, C., Lambert, D. L., & Asplund, M. 2001, ApJ, 556, L63 Andrews, B. H. & Martini, P. 2013, ApJ, 765, 140  \n- Arellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, ApJ, 940, L23 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 Baker, W. M. & Maiolino, R. 2023, MNRAS[ arXiv:2303.08145 ] Baker, W. M., Maiolino, R., Belfiore, F., et al. 2023a, MNRAS, 518, 4767 Baker, W. M., Maiolino, R., Belfiore, F., et al. 2023b, MNRAS, 519, 1149 Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5 Belli, S., Jones, T., Ellis, R. S., & Richard, J. 2013, ApJ, 772, 141 Bian, F., Kewley, L. J., & Dopita, M. A. 2018, ApJ, 859, 175 Birkmann, S. M., Böker, T., Ferruit, P., et al. 2011, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8150, Cryogenic Optical Systems and Instruments XIII, ed. J. B. Heaney & E. T. Kvamme, 81500B Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, arXiv e-prints, arXiv:2301.13766 Böker, T., Birkmann, S., de Marchi, G., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, ed. M. C. Clampin, G. G. Fazio, H. A. MacEwen, & J. Oschmann, Jacobus M., 84423F\n- Bonaventura, N., Jakobsen, P., Ferruit, P., Arribas, S., & Giardino, G. 2023, arXiv e-prints, arXiv:2302.10957\n- Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000\n- Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023a, arXiv e-prints, arXiv:2306.02467\n- Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023b, arXiv e-prints, arXiv:2302.07256\n- Cameron, A. J., Katz, H., Rey, M. P., & Saxena, A. 2023a, MNRAS, 523, 3516 Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023b, arXiv e-prints, arXiv:2302.04298\n- Cappellari, M. 2017, MNRAS, 466, 798\n- Cappellari, M. 2022, MNRAS submitted [ 2208.14974 ]\n- Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875\n- Chabrier, G. 2003, PASP, 115, 763\n- Chevallard, J. & Charlot, S. 2016, MNRAS, 462, 1415\n- Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102\n- Conroy, C., Naidu, R. P., Zaritsky, D., et al. 2019, ApJ, 887, 237\n- Cresci, G., Mannucci, F., & Curti, M. 2018, arXiv e-prints\n- Curti, M., Cresci, G., Mannucci, F., et al. 2017, MNRAS, 465, 1384\n- Curti, M., D'Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425\n- Curti, M., Maiolino, R., Cirasuolo, M., et al. 2020a, Monthly Notices of the Royal Astronomical Society, 492, 821, arXiv: 1910.13451\n- Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020b, Monthly Notices of the Royal Astronomical Society, 491, 944, publisher: Oxford Academic\n- Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2022, arXiv e-prints, arXiv:2212.04568\n- Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98\n- Davé, R., Finlator, K., & Oppenheimer, B. D. 2011, MNRAS, 416, 1354\n- Davé, R., Rafieferantsoa, M. H., Thompson, R. J., & Hopkins, P. F. 2017, MNRAS, 467, 115\n- Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451\n- Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465\n- Ellison, S. L., Patton, D. R., Simard, L., & McConnachie, A. W. 2008, ApJ, 672, L107\n- Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813\n- Ferruit, P., Jakobsen, P., Giardino, G., et al. 2022, A&A, 661, A81\n- Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2022, arXiv e-prints, arXiv:2211.05792\n- Finlator, K. & Davé, R. 2008, MNRAS, 385, 2181\n- Forbes, J. C., Krumholz, M. R., Burkert, A., & Dekel, A. 2014, MNRAS, 443, 168\n- Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306\n- Fraternali, F. & Binney, J. J. 2008, MNRAS, 386, 935\n- Giardino, G., Birkmann, S., Robberto, M., et al. 2019, PASP, 131, 094503  \nGiavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93  \n- Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wol ff , M. J. 2003, ApJ, 594, 279\n- Guo, Y., Koo, D. C., Lu, Y., et al. 2016, ApJ, 822, 103\n- Hayden-Pawson, C., Curti, M., Maiolino, R., et al. 2022, MNRAS, 512, 2867 Heintz, K. E., Brammer, G. B., Giménez-Arteaga, C., et al. 2022, arXiv e-prints, arXiv:2212.02890\n- Henry, A., Martin, C. L., Finlator, K., & Dressler, A. 2013a, ApJ, 769, 148 Henry, A., Scarlata, C., Domínguez, A., et al. 2013b, ApJ, 776, L27\n- Hirschauer, A. S., Salzer, J. J., Janowiecki, S., & Wegner, G. A. 2018, AJ, 155, 82\n- Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80 Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, The Astrophysical Journal Supplement Series, 254, 22\n- Kau ff mann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54 Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531\n- Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121\n- Kewley, L. J., Maier, C., Yabe, K., et al. 2013, ApJ, 774, L10\n- Kobayashi, C. & Taylor, P. 2023, arXiv e-prints, arXiv:2302.07255\n- Langan, I., Ceverino, D., & Finlator, K. 2020, MNRAS, 494, 1988\n- Langeroodi, D. & Hjorth, J. 2023, arXiv e-prints, arXiv:2307.06336\n- Langeroodi, D., Hjorth, J., Chen, W., et al. 2022, arXiv e-prints, arXiv:2212.02491\n- Laseter, I. H., Maseda, M. V., Curti, M., et al. 2023, arXiv e-prints, arXiv:2306.03120\n- Lee, H., Skillman, E. D., Cannon, J. M., et al. 2006, ApJ, 647, 970\n- Lequeux, J., Peimbert, M., Rayo, J. F., Serrano, A., & Torres-Peimbert, S. 1979, A&A, 80, 155\n- Li, M., Cai, Z., Bian, F., et al. 2022, arXiv e-prints, arXiv:2211.01382\n- Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119\n- Looser, T. J., D'Eugenio, F., Maiolino, R., et al. 2023, arXiv e-prints, arXiv:2302.14155\n- Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140 Maier, C., Lilly, S. J., Ziegler, B. L., et al. 2014, ApJ, 792, 3\n- Maiolino, R. & Mannucci, F. 2019, A&A Rev., 27, 3\n- Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463\n- Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023a, arXiv e-prints, arXiv:2308.01230\n- Maiolino, R., Scholtz, J., Witstok, J., et al. 2023b, arXiv e-prints, arXiv:2305.12492\n- Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115\n- Mannucci, F., Cresci, G., Maiolino, R., et al. 2009, MNRAS, 398, 1915  \nMascia,  \nS.,  \nPentericci,  \nL.,  \nCalabro',  \nA.,  \net  \nal.  \n2023,  \narXiv  \ne-prints,  \narXiv:2301.02816  \n- Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2022, arXiv e-prints, arXiv:2211.08255\n- Nakajima, K. & Ouchi, M. 2014, MNRAS, 442, 900\n- Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, arXiv e-prints, arXiv:2301.12825\n- Nakajima, K., Ouchi, M., Xu, Y., et al. 2022, arXiv e-prints, arXiv:2206.02824\n- Noeske, K. G., Faber, S. M., Weiner, B. J., et al. 2007, ApJ, 660, L47\n- Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129 Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and  \nActive Galactic Nuclei, 2nd edn. (University Science Books)  \nPallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621  \nPeeples, M. S. & Shankar, F. 2011, MNRAS, 417, 2962  \n- Péroux, C., Nelson, D., van de Voort, F., et al. 2020, MNRAS, 499, 2462 Pettini, M. & Pagel, B. E. J. 2004, MNRAS, 348, L59\n- Pilyugin, L. S., Mattsson, L., Vílchez, J. M., & Cedrés, B. 2009, MNRAS, 398, 485\n- Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6 Pontoppidan, K., Blome, C., Braun, H., et al. 2022, arXiv e-prints, arXiv:2207.13067\n- Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526  \nReddy, N. A., Topping, M. W., Shapley, A. E., et al. 2022, ApJ, 926, 31  \n- Renzini, A. & Peng, Y.-j. 2015, ApJ, 801, L29\n- Rhoads, J. E., Wold, I. G. B., Harish, S., et al. 2023, ApJ, 942, L14  \nRieke, M. & the JADES Collaboration. 2023, arXiv e-prints, arXiv:2306.02466  \n- Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2022, arXiv e-prints, arXiv:2212.04480  \nSaintonge, A., Catinella, B., Cortese, L., et al. 2016, MNRAS, 462, 1749 Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22 Salim, S., Lee, J. C., Davé, R., & Dickinson, M. 2015, ApJ, 808, 25 Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19 Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2018, ApJ, 858, 99  \n- Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023a, arXiv e-prints, arXiv:2303.08149\n- Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023b, arXiv e-prints, arXiv:2301.06696 Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4 Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150 Shapley, A. E., Coil, A. L., Ma, C.-P., & Bundy, K. 2005, ApJ, 635, 1006 Shapley, A. E., Reddy, N. A., Sanders, R. L., Topping, M. W., & Brammer, G. B. 2023a, arXiv e-prints, arXiv:2303.00410 Shapley, A. E., Sanders, R. L., Reddy, N. A., Topping, M. W., & Brammer, G. B. 2023b, arXiv e-prints, arXiv:2301.03241 Somerville, R. S. & Davé, R. 2015, ARA&A, 53, 51 Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165 Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164 Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023, arXiv e-prints, arXiv:2302.07234 Tacchella, S., Johnson, B. D., Robertson, B. E., et al. 2022, arXiv e-prints, arXiv:2208.03281 Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157 Taylor, A. J., Barger, A. J., & Cowie, L. L. 2022, ApJ, 939, L3 Telford, O. G., Dalcanton, J. J., Skillman, E. D., & Conroy, C. 2016, ApJ, 827, 35 Topping, M. W., Shapley, A. E., Sanders, R. L., et al. 2021, MNRAS, 506, 1237 Torrey, P., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 484, 5587 Tremonti, C. A., Heckman, T. M., Kau ff mann, G., et al. 2004, ApJ, 613, 898 Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110 Trump, J. R., Arrabal Haro, P., Simons, R. C., et al. 2022, arXiv e-prints, arXiv:2207.12388 Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2302.06647 Ucci, G., Dayal, P., Hutter, A., et al. 2023, MNRAS, 518, 3557 Vidal-García, A., Charlot, S., Bruzual, G., & Hubeny, I. 2017, MNRAS, 470, 3532 Vincenzo, F., Matteucci, F., Belfiore, F., & Maiolino, R. 2016, MNRAS, 455, 4183 Williams, H., Kelly, P. L., Chen, W., et al. 2022, arXiv e-prints, arXiv:2210.15699 Witstok, J., Smit, R., Maiolino, R., et al. 2021, MNRAS, 508, 1686 Wuyts, E., Kurk, J., Förster Schreiber, N. M., et al. 2014, ApJ, 789, L40 Wuyts, E., Rigby, J. R., Sharon, K., & Gladders, M. D. 2012, ApJ, 755, 73 Yabe, K., Ohta, K., Iwamuro, F., et al. 2014, MNRAS, 437, 3647 Yang, H., Malhotra, S., Gronke, M., et al. 2017a, ApJ, 844, 171 Yang, H., Malhotra, S., Rhoads, J. E., & Wang, J. 2017b, ApJ, 847, 38 Yates, R. M., Kau ff mann, G., & Guo, Q. 2012, MNRAS, 422, 215 Yates, R. M., Schady, P., Chen, T.-W., Schweyer, T., & Wiseman, P. 2019, arXiv e-prints Zahid, H. J., Bresolin, F., Kewley, L. J., Coil, A. L., & Davé, R. 2012, ApJ, 750, 120 Zahid, H. J., Dima, G. I., Kudritzki, R.-P., et al. 2014, ApJ, 791, 130 Zahid, H. J., Kewley, L. J., & Bresolin, F. 2011, ApJ, 730, 137\n- 1 European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany\n- 2 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 3 Cavendish Laboratory - Astrophysics Group, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK\n- 4 Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK\n- 5 Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK\n- 6 Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy\n- 7 Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France  \n- 8 School of Physics, University of Melbourne, Parkville 3010, VIC, Australia\n- 9 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia\n- 10 Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706 USA\n- 11 Centro de Astrobiología (CAB), CSIC-INTA, Cra. de Ajalvir Km. 4, 28850- Torrejón de Ardoz, Madrid, Spain\n- 12 ATG Europe for the European Space Agency, ESTEC, Noordwijk, The Netherlands\n- 13 European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s / n, 28692 Villanueva de la Cañada, Madrid, Spain\n- 14 Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK\n- 15 NRC Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada\n- 16 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge MA 02138 USA\n- 17 Department of Physics and Astronomy, The Johns Hopkins University, 3400 N. Charles St. Baltimore, MD 21218\n- 18 Steward Observatory University of Arizona 933 N. Cherry Avenue Tucson AZ 85721, USA\n- 19 Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA\n- 20 NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA\n- 21 European Space Agency, ESA / ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, NL\n- 22 AURA for European Space Agency, Space Telescope Science Institute, 3700 San Martin Drive. Baltimore, MD, 21210", 'Appendix A: Metallicity scaling relations under different assumptions on sample selection and derivation of galaxy properties': "In the main body of the paper, we discuss the mass-metallicity (MZR) and the fundamental metallicity (FMR) relations 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 di ff erent choices than our set of fiducial assumptions, particularly regarding sample selection, SFR measurement, and FMR parametrisation.  \nIn 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. 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' starforming galaxies can be applied. Stellar masses and SFRs are also derived assuming negligible AGN contribution to the continuum level. No significant change to the observed evolution in the metallicity scaling relations is found by including these objects.  \nIn Figure A.2 instead we report, similarly to Figure 7, the deviations from the predictions of the local FMR, but this time assuming di ff erent 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 (hence ∆ 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.  \nFinally, in the three panels of Figure A.3 we compare (from left to right) the observed metallicity for our JWST sample to the predictions of the FMR adopting the parametrisation of Andrews &Martini (2013) (based on stacked spectra of SDSS galaxies in bins of M ⋆ and SFR), and in particular considering the projection on the log(O / H) versus µα plane, under di ff erent estimates for the SFR of galaxies, i.e., fiducial, H α -based with low-metallicity calibration, H α -based with solar metallicity calibration. In all cases, the relation from Andrews & Martini (2013) fails to match the bulk of individual galaxies and the median values in bins of M ⋆ and redshift, despite embedding a stronger secondary dependence on SFR ( α = 0.66) into µα , and a steeper slope, compared to di ff erent realisations of the FMR (Mannucci et al. 2010; Curti et al. 2020b, but see the discussion in Section 5). However, we note that the case with H α -based SFR and the calibration from Kennicutt & Evans (2012), which significantly reduces the number of ionizing photons per unit SFR compared to both our fiducial and the H α -based, low-metallicity calibration scenario, provides the best agreement with the FMR predictions of Andrews & Martini (2013), as driven by the higher estimates of the SFR which reduce µα and move the points towards the left on the xaxis of the diagram.", 'Appendix B: Summary of properties for JADES galaxies': "<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nNotes: a deri v ed via bea gle fitting to PRISM spectra and NIRCAM photometry (where a v ailable, see 'NC' column). b deri v ed from H α (H β ) flux adopting the calibration for lo w-metallicity g alaxies from Reddy et al. ( 2022 ). ⋆ Spectral configuration (P = prism; M = medium resolution gratings) considered for deri ving the metallicity . See Section 3.1 for more details. † NIRCAM photometry: for g alaxies with flag = 1 the photometry has been included in the M ⋆ and SFR deri v ation with bea gle ; objects with flag = 0 ha v e M ⋆ and SFR measured only from fitting to the PRISM spectra. ‡ Galaxies identified as candidate (narro w-lines) A GN based on the criteria described in Scholtz et al. (in preparation). These objects are remo v ed from the fiducial analysis presented in the main body of the paper due to the unkno wn A GN contrib ution and the potential biases introduced in the M ⋆ , SFR, and metallicity deri v ation, nonetheless their main properties (estimated self-consistently as for the other sources) are included here (and sho wn in Section A ).  \nFig. A.3. Same as Figure 8, but adopting the formalism from Andrews & Martini 2013, based on SDSS stacked spectra in bins of M ⋆ and SFR. The three panels show, from left to right, the measured versus FMR-predicted metallicities based on three di ff erent SFR estimates, i.e., the fiducial value adopted throughout the paper (and based on beagle fitting to both R100 spectra and photometry), the SFR from the H α flux and the calibration for low-metallicity galaxies from Reddy et al. 2022; Shapley et al. 2023b, 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 calibration for metal poor galaxies, decreasing µα and improving the agreement with the predictions of the Andrews & Martini 2013 parametrisation of the FMR.  \n<!-- image -->  \n¯  \nFig. 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-deep (Scholtz et al., in preparation), here marked as dashed hexagons. No clear change in the median trends of the metallicity scaling relations is found, and the main results reported in the paper are una ff ected.  \n<!-- image -->  \n<!-- image -->  \nFig. 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.  \n<!-- image -->  \n<!-- image -->  \n¡  \n<!-- image -->  \n¡  \n<!-- image -->  \n¡"}