We present the evolution of the massmetallicity MZ relation at z  410 derived with 135 galaxies identified in JWSTNIRSpec data taken from the three major public spectroscopy programs of ERO GLASS and CEERS. Because there are many discrepancies between the flux measurements reported by the early ERO studies we first establish our NIRSpec data reduction procedure for reliable emissionline flux measurements and errors successfully explaining Balmer decrements with no statistical tensions thorough comparisons with the early ERO studies. Applying the reduction procedure to the 135 galaxies we obtain emissionline fluxes for physical property measurements. We confirm that 10 out of the 135 galaxies with O III 4363 lines have electron temperatures of 1.12.3  10SUP4SUP K similar to lowerz starforming galaxies which can be explained by heating by young massive stars. We derive the metallicities of the 10 galaxies by a direct method and the rest of the galaxies with strong lines using the metallicity calibrations of Nakajima et al. applicable for these lowmass metalpoor galaxies anchoring the metallicities with the directmethod measurements. We thus obtain the MZ relation and star formation rate SFRMZ relation over z  410. We find that there is a small evolution of the MZ relation from z  23 to z  410 while interestingly the SFRMZ relation shows no evolution up to z  8 but a significant decrease at z gt 8 beyond the errors This SFRMZ relation decrease at z gt 8 may suggest a break of the metallicity equilibrium state via star formation inflow and outflow while further statistical and localbaseline studies are needed for a conclusion.
['Chemical abundances', 'Galaxy chemical evolution', 'Galaxy evolution', 'James Webb Space Telescope', 'High-redshift galaxies', '224', '580', '594', '2291', '734', 'Astrophysics - Astrophysics of Galaxies']
JWST Census for the MassMetallicity Star Formation Relations at z  410 with Selfconsistent Flux Calibration and Proper Metallicity Calibrators
{'JWST Census for the Mass-Metallicity Star-Formation Relations at z = 4 -10 with the Self-Consistent Flux Calibration and the Proper Metallicity Calibrators': 'Kimihiko Nakajima , 1 Masami Ouchi , 1, 2, 3 Yuki Isobe , 2, 4 Yuichi Harikane , 2 Yechi Zhang , 2, 4 Yoshiaki Ono , 2 Hiroya Umeda, 2, 4 and Masamune Oguri 5, 6  \n1 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 2 Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan 3 Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan 4 Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan 5 Center for Frontier Science, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan 6 Department of Physics, Graduate School of Science, Chiba University, 1-33 Yayoi-Cho, Inage-Ku, Chiba 263-8522, Japan  \nSubmitted to ApJ Supplement on Jan 30, 2023 / Accepted on May 8, 2023', 'ABSTRACT': 'We present the evolution of the mass-metallicity (MZ) relations at z = 4 -10 derived with 135 galaxies identified in the JWST/NIRSpec data taken from the three major public spectroscopy programs of ERO, GLASS, and CEERS. Because there are many discrepancies between flux measurements reported by early ERO studies, we first establish our NIRSpec data reduction procedure for reliable emission-line flux measurements and errors successfully explaining Balmer decrements with no statistical tensions via thorough comparisons of the early ERO studies. Applying the reduction procedure to the 135 galaxies, we obtain emission-line fluxes for physical property measurements. We confirm that 10 out of the 135 galaxies with [O iii ] λ 4363-lines have electron temperatures of /similarequal (1 . 1 -2 . 3) × 10 4 K, similar to lowerz star-forming galaxies, that can be explained by heating of young massive stars. We derive metallicities of the 10 galaxies by the direct method and the rest of the galaxies with strong lines by the metallicity calibrations of Nakajima et al. (2022) applicable for these low-mass metal-poor galaxies, anchoring the metallicities with the direct-method measurements. We thus obtain MZ relations and star-formation rate (SFR)-MZ relations over z = 4 -10. We find that there is a small evolution of the MZ relation from z ∼ 2 -3 to z = 4 -10, while interestingly that the SFR-MZ relation shows no evolution up to z ∼ 8 but a significant decrease at z > 8 beyond the error. This SFR-MZ relation decrease at z > 8 may suggest a break of the metallicity equilibrium state via star-formation, inflow, and outflow, while further statistical and local-baseline studies are needed for a conclusion.  \nKeywords: Chemical abundances(224) - Galaxy chemical evolution(580) - Galaxy evolution(594) High-redshift galaxies(734) - James Webb Space Telescope(2291)', '1. INTRODUCTION': "The James Webb Space Telescope (JWST) has dramatically advanced the exploration of the high redshift universe. One important advancement for the high redshift community is provided by the high sensitivity of the Near Infrared Spectrograph (NIRSpec; Jakobsen et al.  \nCorresponding author: Kimihiko Nakajima  \nkimihiko.nakajima@nao.ac.jp  \n2022) at λ /similarequal 2 -5 µ m, allowing us to directly address the key questions relating to the physical conditions of inter-stellar medium (ISM) and the nature of the ionizing spectrum for galaxies in the early universe. The properties of the hot ISM, especially the gasphase metallicity determined by the oxygen abundance (12 + log(O / H)), can be accurately constrained using rest-frame optical emission lines. These emission lines have been well calibrated and understood over decades, as discussed in a recent review by Maiolino & Mannucci (2019). Gas-phase metallicities provide a crucial exper-  \nental tool for studying the early chemical enrichment and feedback processes in galaxy evolution. These measurements can be used in conjunction with cosmological simulations (e.g., Finlator & Dav'e 2008; Ma et al. 2016; Torrey et al. 2019; Langan et al. 2020; Ucci et al. 2021) and chemical evolution models (e.g., Dav'e et al. 2012; Lilly et al. 2013; Dayal & Ferrara 2018) to gain insights into galaxy formation and evolution. Moreover, the efficiency of ionizing photon production, as quantified by ξ ion (the production rate of hydrogen ionizing photons per unit luminosity in the UV-continuum), can be used to study the nature of the ionizing spectrum. The Hydrogen Balmer emission via recombination physics provides the best constraints on ξ ion , while high ionization lines such as He ii offer key spectroscopic tools to investigate the hardness of the ionizing spectrum. These approaches are crucial for understanding the stellar population of early galaxies, diagnosing the presence of accreting black holes in the system, and detecting galaxies hosting the very first generation of stars in the early universe (e.g., Schaerer 2003; Inoue 2011; Kewley et al. 2013; Bouwens et al. 2016; Nakajima & Maiolino 2022; Katz et al. 2022; Trussler et al. 2022).  \nDeterminations of these properties ideally require key optical emission lines such as [O iii ] λλ 5007 , 4959, [O iii ] λ 4363, [O ii ] λλ 3726 , 3729 1 , and the Balmer lines. In particular, [O iii ] λ 4363 is crucial for gas-phase metallicity studies based on electron temperature ( T e ). However, these lines are limited in their applicability to galaxies up to z /lessorsimilar 3 in the pre-JWST era, with the detection of [O iii ] λ 4363 from sources at z /greaterorsimilar 1 still challenging due to faintness (e.g., Christensen et al. 2012; Jones et al. 2015; Sanders et al. 2020, 2021). These lines are beyond the reach of ground-based telescopes at higher redshift. Now JWST/NIRSpec has become online to directly examine these key emission lines in detail for high redshift sources.  \nUsing the Early Release Observations (ERO) of NIRSpec within one month after the data release, several studies report a detection of [O iii ] λ 4363 from three objects at z = 7 . 6 -8 . 5 which are magnified by the strong lensing galaxy cluster SMACS J0723.3-7327 (Schaerer et al. 2022; Curti et al. 2023a; Trump et al. 2022; Rhoads et al. 2023; Arellano-C'ordova et al. 2022; Brinchmann 2022). These objects have provided the first reliable metallicity determinations at high redshift based on the direct T e measurement. The metallicity values from the different teams are overall consistent to each other for each source, resulting in a mass-  \ntallicity relation which is scattered around the extrapolation of the z /similarequal 2 -3 relation towards the lower mass end ( M /star = 10 7 . 5 -10 9 M /circledot ). Curti et al. (2023a) also examine the mass-metallicity-star formation rate (SFR) relation as indicated to be a fundamental relation from z = 0 to z /lessorsimilar 2 -3 (Mannucci et al. 2010; Maiolino & Mannucci 2019; Sanders et al. 2021, and the references therein), suggesting that no clear evolution is seen from z = 0 -3 to z = 7 . 6 -8 . 5, particularly excepting for the z = 8 . 5 object (see below). However, the sample size is obviously too small to conclude the evolution on the mass-metallicity relation and its SFR dependence.  \nMoreover, several pieces of evidence in the previous ERO studies imply that the early release of the NIRSpec products contain some issues. At the early stage of the ERO data release, the reductions including the flux calibration are partly based on the predicted preflight data such as the throughputs of the spectrograph, the optical telescope element, and the NIRSpec fore optics. A sensitivity function, whose shape is strongly dependent on the wavelength, is necessary for appropriate flux measurements as adopted in Curti et al. (2023a) using standard stars taken during the commissioning. Even after the flux-calibration, there remains some tensions for faint emission lines. A notable issue is seen in the Balmer decrements such as H γ /H β and H δ /H β which are (not always but sometimes) unexplainable in a consistent way within the uncertainties according to the Case B recombination. Such tensions, especially seen in the faint emission lines, may be caused by background residuals and/or hot pixels that persist in the released final spectra. Interestingly, the z = 8 . 5 object (ID:04590) is suggested to present a high [O iii ] λλ 4363 / 5007 line ratio, and accordingly a very high T e /greaterorsimilar 25000K, which is hard to be explained by heating of young massive stars alone (Katz et al. 2023). The low metallicity indicates it falls significantly below the mass-metallicity-SFR relation (Curti et al. 2023a) and can be in an early stage of galaxy evolution. Such a high electron temperature and low metallicity need to be confirmed by a more carefully reduced spectrum free from any tensions. The NIRSpec data reduction procedure will be revisited as one of the key parts of this paper for reliable emission-line flux measurements and errors.  \nDespite these possible uncertainties, recent NIRSpec observations are now providing increasing numbers of spectra from high-redshift objects, and the reports of the identifications of strong emission lines such as [O iii ] λλ 5007 , 4959+H β from objects at z > 7 -8 follows. Notable results has come from the Early Release Science (ERS) observations of GLASS (Treu et al.  \n2022; Morishita et al. 2022; Mascia et al. 2023) and CEERS (Finkelstein et al. 2022; Bagley et al. 2022; Sanders et al. 2023; Tang et al. 2023; Fujimoto et al. 2023; Shapley et al. 2023), two Director's Discretionary Time (DDT) programs (Williams et al. 2022; Heintz et al. 2022; Wang et al. 2022; Langeroodi et al. 2022), and the Guaranteed Time Observations of JADES (Robertson et al. 2022; Curtis-Lake et al. 2022; Bunker et al. 2023; Cameron et al. 2023b; Saxena et al. 2023; Curti et al. 2023b). Regarding the ISM properties, Sanders et al. (2023) examine the ionization properties of ∼ 160 galaxies at z = 2 -9 from CEERS by using high-to-low ionization emission line ratios such as [O iii ] λ 5007 / [O ii ] λ 3727, and suggest that galaxies tend to present hard ionizing spectra at high-redshift. The metallicities at z > 6 . 5 are suggested to be sub-solar on average (see also Tang et al. 2023; Shapley et al. 2023). A similarly high ionization state and modest metallicity is also suggested in 29 galaxies at z = 4 . 8 -8 from GLASS (Mascia et al. 2023) and in 26 galaxies at z = 5 . 5 -9 . 5 from JADES (Cameron et al. 2023b, see also Saxena et al. 2023). Fujimoto et al. (2023) explicitly derive metallicities of ∼ 10 z = 8 -9 CEERS galaxies using a metallicity indicator of [O iii ] λ 5007/H β , suggesting that these high-redshift galaxies have a lower metallicity than z = 0 -3 galaxies for a given stellar mass. More extremely, Williams et al. (2022) report a z = 9 . 5 galaxy with [O iii ]+H β , revealing a very high [O iii ] / [O ii ] ratio and a modestly low-metallicity, 12 + log(O / H) = 7 . 48 ± 0 . 08, falling below the z = 0 mass-metallicity relation but still consistent within 2 σ . Including the z = 9 . 5 object, Heintz et al. (2022) analyzed five objects at z > 7 . 8 from the DDT programs and suggest a systematic decrease of metallicity at z > 7 . 8 on the mass-metallicity-SFR relation found in lowerredshift, as seen in the ERO z = 8 . 5 object. A similar conclusion is also derived by Langeroodi et al. (2022). Boyett et al. (2023) report another z = 9 galaxy from GLASS whose [Ne iii ] λ 3869/[O ii ] ratio implies a subsolar metallicity. Finally, Bunker et al. (2023) report the astonishing spectrum of GN-z11 at z = 10 . 6 using the deep JADES observations. The author use several metal lines such as [Ne iii ], [O ii ], as well as the UV lines of N iv ] λ 1486, N iii ] λ 1748, and C iii ] λ 1909 and imply an unusually high nitrogen-to-oxygen abundance ratio for its modest ( ∼ sub-solar) oxygen abundance (see also Cameron et al. 2023a). We note that these metallicities are based on the strong emission line ratios which are empirical indicators and calibrated mostly in the local universe (e.g., Maiolino & Mannucci 2019). The applicability of these methods at high-redshift need to be carefully confirmed (e.g., Curti et al. 2023a), partic-  \nindication of different degrees of ionization in the ISM of galaxies between z = 0 and highredshift for a fixed metallicity (Sanders et al. 2023).  \nIn this paper, we provide a summary of three major public NIRSpec observation programs of ERO, GLASS, and CEERS, aimed at characterizing the massmetallicity relation at z = 4 -10 and examining its evolution over cosmic time. We establish a reliable data reduction procedure for emission-line flux measurements and errors using NIRSpec data, which is applied to construct a large sample of galaxies at z = 4 -10 as detailed in Section 2. In Section 3, we utilize key emission line ratios to derive metallicities for the JWST objects with improved NIRSpec spectra. We then investigate the mass-metallicity relation and its dependence on SFR at high redshift. We begin with the 10 objects in which [O iii ] λ 4363 is identified and metallicity is determined using the direct T e method. The empirical metallicity indicators validated with these 10 objects are then applied to estimate metallicities for the remaining JWST objects. We analyze the mass-metallicity and mass-metallicity-SFR relations at z = 4 -10 using the metallicity measurements obtained for the JWST objects. The results and implications of these metallicity relationships are discussed in Section 4, and we provide a summary of our conclusions. Throughout the paper we adopt a standard ΛCDM cosmology with Ω Λ = 0 . 7, Ω m = 0 . 3, and H 0 = 70kms -1 Mpc -1 .", '2. DATA': 'We assemble the publicly available NIRSpec data of ERO and ERS, re-reduce and analyze the spectra, and identify high-redshift galaxies with the rest-frame optical emission lines to discuss their nebular properties. In the following subsections, we detail our improved reduction procedures, particularly for the ERO data to demonstrate the improvements by comparing with the early results.', '2.1. ERO Data': "2.1.1. NIRSpec observations and Data reduction  \nThe JWST/NIRSpec ERO observations were undertaken on UT 2022 June 30, targeting the SMACS 0723 lensing cluster field (Pontoppidan et al. 2022; Proposal ID: 2736). Multi-slit spectroscopy was taken using the micro-shutter assembly (MSA), with the medium resolution ( R ∼ 1000) gratings/filters of G235M/F170LP and G395M/F290LP sampling the wavelength range of 1 . 7 -3 . 1 and 2 . 9 -5 . 1 µ m, respectively. Two observing blocks were carried out, o007 and o008, each of which consisted of three exposures of 2918 sec in G395M and F235M each, with a 3-shutter slitlet nod pattern, i.e.,  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 1. Observed NIRSpec spectra for the five z = 5 -8 . 5 EROobjects. For each object: (top) 2D spectra for the six individual exposures (three nods each in o007 and o008). (bottom) 1D composite spectrum. We generate the 1D composite spectrum by first median-stacking the individual 2D spectra and then extracting the 1D spectrum with a small aperture (see text). This procedure allows us to minimize the possible impacts of hot pixels and to rescue some individual 2D frames where the spectrum is partly out of the slit (e.g., two nod2 frames of 10612).  \n<!-- image -->  \nslightly nodded positions, nod1, nod2, and nod3, in different shutters within the slitlets. The total integration time was 4 . 86hr in G235M and G395M each.  \nBy visually checking all of the 1D and 2D composite spectra that were publicly available, we identified 5 objects whose [O iii ] λλ 5007 , 4959+H β were clearly detected and that were useful for the investigation of the nebular properties at high-redshift z > 5, with IDs 04590, 05144, 06322, 08140, and 10612. Although there were additional two sources in the ERO sample that were tentatively reported at z = 5 -6 (Brinchmann 2022; Mahler et al. 2022), they were not included in the following analysis because of their lines' poorer confidence levels.  \nFor the five sources, data reduction was re-performed for a more careful analysis and flux calibration. Starting with the Level1 products provided by the observatory, we carried out the Spec2 and Spec3 pipelines per nod and per observing block using the Python library for JWST observations by STScI (ver.1.8.5). We used the latest reference files stored in a pmap file of either 1028 or 1027, where notably the flat files were created from in-flight data for all MOS data taken after the launch, resolving the MOS flux calibration issues reported in the early report of the NIRSpec ERO data (see below for more details about the flux calibration). According to the original reduction procedure for NIRSpec MSA data provided by the CEERS team 2 , during the Spec2 pipeline, a single science image for a particular nodding position was background-subtracted by using the other nods' images. However, we realized some hot pixels including uncleaned cosmic rays in one of the background images fell on the spatial position of the science objects and badly influenced the final 1D spectrum. To avoid such additional noises, we carefully checked the background-subtracted science images and removed one of the other nods' images, if necessary, for the background images and re-ran the Spec2 pipeline in an iterative manner. Furthermore, we subtracted the local residual background at each spectral pixel for each science image by masking the spatial position of the object and taking the median value of the pixels in the range [ -10pix:+10pix] along the spectral direction 3 . We performed the Spec2 pipeline assuming that all of the five sources were point sources, as they all looked compact in the NIRCam images with the half-  \nlight radii smaller than the PSF (e.g., Harikane et al. 2022a; Tacchella et al. 2022; see also Ono et al. 2022 for the size evolution at high-redshift). This assumption affected how the pathloss corrections, which accounted for various types of signal loss in NIRSpec data, were calculated and applied. The correction depended on the position of the source within the configured MSA shutter, the aperture used for the observations, the wavelength, as well as the type of the object. We used the original MSA metadata files to refer to the relative placement of the sources within the MSA shutters assuming the pointing uncertainties during the observations were negligible, but manually changed the values in the column STELLARITY to force the objects to be treated as point sources. Figure 1 summarizes the six 2D outcomes (3 nods for 2 observing blocks) for each of the ERO sources.  \nChecking the six single exposures in both 1D and 2D for each of the sources, we confirmed no signal was detected in the nod2/o007 observing block for 04590, as already pointed out by the earlier studies (e.g., Curti et al. 2023a; Arellano-C'ordova et al. 2022). This exposure obviously needed to be removed from the final composition. Moreover, we noticed the spectra of nod2/o008 of 05144 were badly affected especially around the H γ +[O iii ] λ 4363 wavelength regime. To be conservative, we decided to remove the visit for 05144 to make a composite. In short, the integration time we used for the composite spectra of 04590 and 05144 was 4 . 05hr, and 4 . 86hr for 06355, 08140, and 10612.  \nThe composition was done by median-stacking the available 5-6 2D spectra. The individual 2D spectra were residual background-subtracted, and shifted to have common spatial and spectral ranges for the composition. All the shifts were integer movements. We adopted the 2D median-stacking procedure to further alleviate the possible impacts of hot pixels that persisted in the final spectra, which were produced by the standard procedure by the Spec3 pipeline. Moreover, the 2D median-stacking method allowed us to rescue the frames where the spectrum was partly out of the slit. Examples were found in the two nod2 positions for 10612 whose losses of the fluxes were > 3 σ level if extracted in individual frames and compared with the other nods' 1D. Our method treated the outer regions of the slit as NaN and calculated the median value without NaN at each pixel in 2D. Using the 2D composite, one-dimensional spectrum was produced for each source and grating via the summation of 3 pixels (= 0 . '' 3 /similarequal 2 × FWHM) along the spatial direction centered on the spatial peak position. We adopted the narrow extraction aperture as compared to the previous studies to minimize the effects coming from the noisy regions close to the edge, partic-  \nfor 04590, 08140, and 10612 that were observed with a short MSA slitlet.  \nThe flux calibration was one key procedure to be improved for the early NIRSpec studies, as some of the reference files needed for the successful flux calibration, such as the flat frames for the spectrograph and the fore optics, were based on the predicted pre-flight throughputs. As of the pmap file of 1022, these flat files have been updated and created from in-flight data for all MOS data. The MOS flux calibration accuracy is now estimated at approximately 5 percent or less over the full wavelength range 4 . We therefore used the flux solutions provided by the spec2pipeline together with the latest reference files. Finally, we scaled the spectrum to match the broadband photometry of JWST/NIRCam to combine the spectroscopic and photometric measures based on a common aperture to derive the quantities such as equivalent widths of the optical emission lines and ξ ion . We used the NIRCam broadband that covered the [O iii ]+H β emission for the scaling. The final composite 1D spectrum is presented in the bottom panel for each of the sources in Figure 1.  \nTo evaluate the noise level, we used the 2D noise frames of read-out noise ( VAR\\_RNOISE ) and Poisson noise ( VAR\\_POISSON ) outputted in the Spec2 procedure for each object and nod. In addition, we added the standard deviation of the residual background at each spectral bin in quadrature. We then averaged the noise frames of the available nods and extracted the 1D noise spectra in quadrature by using the same spectral and spatial ranges as adopted for the science spectra. Finally the scaling was performed in the same manner as done for the science spectra.", '2.1.2. Emission line flux measurements': "We measured the key optical emission lines of hydrogen Balmer lines, [O ii ] λ 3727, [Ne iii ] λ 3869, [O iii ] λ 4363, and [O iii ] λλ 5007 , 4959 by performing a Gaussian profile fitting to each line, and used the noise spectrum to weight the fit. In the fitting procedure for the relatively faint objects of 04590, 05144, and 08140, the redshift of the [O iii ] λ 5007 (i.e., the strongest emission line) was adopted as Gaussians for the other lines. We perform all wavelength measurements on the vacuum wavelength scale, although we follow a longstanding convention of specifying the emission lines by their wavelength in air. The error of each flux measurement was calculated by adding the noise levels of the spectral bins in quadrature within the wavelength range of ± FWHM centered on the Gaussian peak. Using the flux measure-  \nFigure 2. Balmer emission line ratios of H γ /H β and H δ /H β for the ERO objects (04590 in red, 05144 in emerald green, 06355 in orange, and 10612 in green) in the main panel. The colored filled circles represent the new emission line flux measurements and errors with our improved reduction, while the other symbols show the measurements based on the early data release, as shown in the legend (Curti et al. 2023a; Rhoads et al. 2023; Arellano-C'ordova et al. 2022). For the measurements of Arellano-C'ordova et al. (2022), we adopt the average values based on the two spectra (o007 and o008) for 06355 and 10612, but refer only to the o008 value for 04590. For those with line ratios outside the ranges of the plot, we place them close to the edge of the plot with an arrow to indicate the direction toward the actually reported line ratio. The lines show expected Balmer decrements as a function of E(B -V) for the three different attenuation curves (Cardelli et al. 1989 with dashed, Calzetti et al. 2000 with dotted, and SMC (Gordon et al. 2003) with dot-dashed) and the two different electron temperatures ( T e = 10000 K in gray and 25000 K in magenta), although the different attenuation curve and T e have a small impact on this diagram. The deviation of the H δ /H β ratio from the theoretically expected value based on the H γ /H β ratio (∆ H δ /H β = ((H δ /H β )obs -(H δ /H β )expected) normalized by (H δ /H β ) expected ) is highlighted in the upper panel, and that of the H γ /H β ratio based on the H δ /H β ratio in the left panel, for each symbol as in the main panel. The corresponding E(B -V) values are printed in these sub-panels, where the Calzetti et al. 2000 attenuation curve and T e = 17500 K are adopted for reference. Our measurements all follow the sequence of the expected Balmer decrements within the uncertainties, while quite a few of the measurements from the earlier studies show the ratios that are not physically and/or consistently explainable.  \n<!-- image -->  \nnd the error evaluations, we confirm that four of the ERO objects, 04590, 05144, 06355, and 10612 present a > 3 σ detection of [O iii ] λ 4363. We also derived equivalent widths of the optical emission lines by combining the fluxes with the continuum level provided from a best-fit SED to the appropriate rest-frame optical broad-band photometry (Section 2.1.3). We summarize the key emission line measurements in Table 1 for the [O iii ] λ 4363-detected sources to discuss the metallicity indicators (Section 3.2). For the full sample, we list the metallicity diagnostic emission line ratios in Table D1 in Appendix D together with the other physical quantities.  \nTo test our new flux measurements with the improved reduction and calibration of the NIRSpec spectra, we show in Figure 2 the Balmer decrements of H γ /H β and H δ /H β for the four ERO objects with a significant detection of the three lines, and compare them with the prediction of Case B recombination. Observed Balmer line ratios should follow the sequence of the expected Balmer decrements as a function of dust attenuation, as shown with the curves adopting different reddening laws (Cardelli et al. 1989; Calzetti et al. 2000; Gordon et al. 2003) and electron temperatures ( T e = 10000 and 25000K). The four ERO objects all follow the sequence and present physically-reasonable Balmer decrements within the uncertainties. The other object (08140) is not detected with H γ and H δ but the non-detections as well as the H α /H β line ratio are consistent with the theoretically expected line ratios at the /lessorsimilar 3 σ level. In contrast, previous studies claim some tensions of the observed Balmer decrements based on the early release products which show inconsistent with the expected line ratios, as referenced with the open and the cross symbols in Figure 2. We speculate such inconsistencies have been improved by implementing the residual background subtractions and the 2-dimensional stacking procedure, where we can ease the potential impacts of hot pixels that could have persisted in the final spectra in the early release. We note that the error-bars of our measurements are slightly larger than those in the previous studies (Curti et al. 2023a; Arellano-C'ordova et al. 2022), mostly due to the addition of the uncertainty of the subtracted residual background into the noise levels.  \nPrior to quantitative analysis, it is necessary to consider corrections for dust reddening. We use the observed Balmer decrement of H γ /H β for the four objects to estimate the dust attenuation assuming the Calzetti et al. (2000) attenuation curve and a Case B recombination with T e = 17500K. Adding a lower signalto-noise ratio emission of H δ does not change the estimate as indicated in Figure 2. For 08140, we use the H α /H β ratio instead. We adopt E(B -V) = 0 . 22 for  \n04590, = 0 . 03 for 05144, = 0 . 08 for 10612, and = 0 . 00 for 08140 and 06355, and do not propagate the uncertainty in E(B -V) in the following analysis.", '2.1.3. NIRCam photometry and SED fit': "All five ERO objects have been observed with JWST/NIRCam in the F090W, F150W, F200W, F277W, F356W, and F444W filters. These photometric data are crucial for deriving stellar properties such as stellar mass, and for evaluating the continuum level to correct for slit-loss in the NIRSpec observations. The derived corrections are important for obtaining key properties, such as the total SFR via recombination physics, ξ ion , and EWs of optical emission lines such as EW(H β ).  \nWe refer to the catalog provided by Harikane et al. (2022a) to extract the necessary photometric data. We adopt the total magnitudes that are estimated from the 0 . '' 3-diameter aperture magnitudes with the aperture corrections. The uncertainties of the magnitudes include the photometric errors from the images as well as the 10 % error floor to account for systematic uncertainties arising from e.g., the zero-point corrections.  \nUsing the total magnitudes and the errors of NIRCam, we perform the Bayesian spectral energy distribution (SED) fitting with Prospector (Johnson et al. 2021) to derive the stellar population properties. The procedure of the SED fitting is the same as that of Harikane et al. (2022a), except for the fixed redshifts based on the NIRSpec emission line measurements. We use the stellar population synthesis package, Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) for stellar SEDs, and include nebular emission from the photoionization models of Cloudy (Byler et al. 2017). We assume the stellar initial mass function (IMF) of Chabrier (2003), the intergalactic medium (IGM) attenuation model of Madau (1995), the Calzetti et al. (2000) dust attenuation law, and a fixed metallicity of 0 . 2Z /circledot . We choose a flexible star formation history as adopted in Harikane et al. (2022a).  \nThe stellar masses we obtain from the SED fitting are the apparent values, i.e., after experiencing the strong lensing magnifications for the ERO objects. To correct for the lensing effects, we refer to the mass model of SMACS J0723 constructed with glafic (Oguri 2010, 2021) updated with the new JWST ERO data, as adopted in Harikane et al. (2022a). The adopted magnification factors as well as the stellar masses after the lensing magnification correction are summarized in Table D1. The errors of the best-fit parameters from the SED fitting correspond to the 1 σ confidence interval (∆ χ 2 < 1) for each parameters.  \n2.2. ERS GLASS Data  \nFigure 3. Redshift distribution of the high-redshift galaxy samples used in this work. The red histogram presents the JWST sample constructed in this study using our analyzed ERO, GLASS, and CEERS spectra. The blue histogram shows the samples compiled from the literature (see Section 3.2.3). The sum of these samples, represented by the gray histogram, is used to investigate the mass-metallicity relationships at z = 4 -10 in this work.  \n<!-- image -->  \nThe ERS GLASS NIRSpec MSA observations were carried out on UT 2022 Nov 10-11 for objects behind the galaxy cluster Abell 2744 (Proposal ID: 1324; Treu et al. 2022). They adopted the R ∼ 2700 high resolution gratings/filters of G140H/F100LP, G235H/F170LP, and G395H/F290LP sampling the wavelength range of 1 . 0 -1 . 9, 1 . 7 -3 . 2, and 2 . 8 -5 . 3 µ m, respectively. Twelve exposures of 1473 sec were conducted in one pointing of NIRSpec MSA for each spectral configuration with the 3-shutter slitlet nod pattern, totaling the on-source integration of 4 . 9hr in each of G140H, G235H, and G395H.  \nThrough the initial screening of all of the public 1D and 2D composite spectra, we identified 15 objects at z > 4 with [O iii ]+H β 5 . The sample contains the members of the protocluster at z = 7 . 89 as presented by Morishita et al. (2022). We re-reduced the spectra of the 15 objects with our custom procedure as described in Section 2.1.1. We carefully checked all of the 12 nods' spectra for each objects to determine the extraction aperture and any nods to be removed from the com-  \nosition. For 10005, the light from another member of the z = 7 . 89 protocluster partly came into the same slit at a spatially-offsetted position. We changed the background frames during the spec2 procedure so that the signal from 10005 was not over-subtracted by the neighboring object in different nods.  \nThe emission line flux measurements were performed in the same manner as explained in Section 2.1.2. Among the 15 objects, we identified 5 with [O iii ] λ 4363. One of the [O iii ] λ 4363-detected source, 40066, lacked the [O ii ] λ 3727 emission which fell in the detector gap. The other 4 objects, whose flux measurements are given in Table 1, are thus used for the metallicity measurements with the direct T e method (Section 3.2). The nebular dust attenuation was evaluated using H α , H β , and/or H γ .  \nThe NIRCam multi-band photometric catalog was generated in the same manner as in Section 2.1.3 based on the images publicly available from the GLASS and UNCOVER (Proposal ID: 2561) programs. We use the data of F115W, F150W, F200W, F277W, F356W, and F444W for the subsequent SED fitting analysis. Note that one of 15 objects (160122) is not covered by NIRCam, and its stellar mass are not determined.  \nAll of the GLASS objects are strongly lensed. To correct for the magnification factors, we use the glafic (Oguri 2010, 2021) mass model of Kawamata et al. (2018) with a minor update based on the Multi Unit Spectroscopic Explorer (MUSE) follow-up spectroscopy (Bergamini et al. 2022). The updated mass model uses 45 multiple image systems, 27 of which have spectroscopic redshifts. This model reproduces the observed multiple image positions with an accuracy of 0 . '' 42. The derived magnification factors, the stellar masses after the lensing magnification correction, as well as the other key quantities for the entire GLASS sample are summarized in Table D1 in Appendix D.", '2.3. ERS CEERS Data': "The ERS CEERS NIRSpec MSA observations were carried out on UT 2022 Dec 20-22, and 24 in the blank field of EGS (Finkelstein et al. 2022; Fujimoto et al. 2023; Proposal ID: 1345). They used both of the Prism ( R ∼ 100) and the medium resolution ( R ∼ 1000) gratings to cover from 1 . 0 to 5 µ min 6 pointings (P4, P5, P7, P8, P9, and P10) with some overlaps of objects. Note that two of the Prism pointings, P9 and P10, were not fully completed. To compensate the two Prism pointings, additional observations were performed on UT 2023 Feb 9-10 (P11 and P12; see also Arrabal Haro et al. 2023). These additional data were also analyzed and presented in this paper to provide a complete spectro-  \n<!-- image -->  \nFigure 4. O32 vs. R23 diagram comparing the JWST sample at z = 4 -9, analyzed in this work, with lower redshift samples at z = 2 -3. The JWST objects are shown as red filled symbols, with circles representing ERO, pentagons representing GLASS, and diamonds representing CEERS. Objects with metallicity measurements obtained using the direct T e method are indicated with red open circles. The gray open symbols represent the averages of lower-redshift samples of continuum-selected galaxies (Nakajima & Ouchi 2014; Troncoso et al. 2014; Sanders et al. 2016; Onodera et al. 2016; Strom et al. 2017) and Ly α emitting galaxies (LAEs; Nakajima & Ouchi 2014; Erb et al. 2016; Nakajima et al. 2020), as indicated in the legend. Arrows indicate 3 σ lower limits. Gray shading illustrates the equivalent distribution for nearby SDSS galaxies.  \n<!-- image -->  \nc sample of CEERS to the community. Three exposures of 1036 sec were conducted for each spectral configuration with the 3-shutter slitlet nod pattern, totaling the on-source integration of 0 . 9hr in each Prism, G140M, G235M, and G395M, excepting for P9 and P11 where another set of three exposures of 1036sec were repeated in the Prism observations.  \nAt redshifts z > 4, we detected a total of 147 Prism spectra and 77 medium grating spectra with [O iii ]+H β emission lines. Among these, three objects (IDs 00618, 00717, and 01027) were observed with two different Prism and medium grating pointings each, resulting in four duplicate observations. Additionally, one object (ID 00603) was observed with two different Prism pointings and one medium grating pointing, resulting in three duplicate observations. Furthermore, 50 objects were observed with both Prism and medium grating, resulting in two duplicate observations with one pointing each. In summary, a total of 163 objects were identified through the ERS CEERS NIRSpec observations, as summarized in Table D1 in Appendix D. For the 163 objects, we performed our re-reduction procedure as described in Section 2.1.1.  \nThe same emission line flux measurements were performed as detailed in Section 2.1.2. We identified two out of the 163 objects with [O iii ] λ 4363 among the  \nCEERS sample. Their flux measurements are added in Table 1. The nebular dust attenuation was evaluated using H α , H β , and/or H γ depending on the detections and the redshift.  \nThe CEERS NIRCam observations consists of 10 pointings, four of which were taken in Jun 2022 and the rest in Dec 2022 (Bagley et al. 2022). The early observations of 4 pointings are presented in Harikane et al. (2022a). We reduce the latter observations' data in the same manner, and generate the NIRCam full photometry catalog, including the seven bands of F115W, F150W, F200W, F277W, F356W, F410M, and F444W (Y. Harikane et al. in prep.). One caveat is that not all of the NIRSpec CEERS objects are covered by NIRCam. We find 91 of the 163 CEERS objects are in the NIRCam coverage. For the 91 objects, we perform the SED fitting as in Section 2.1.3 to derive the stellar masses, and scale the NIRSpec spectra to be consistent with the NIRCam total magnitude photometry (i.e., slit-loss correction). For the others except for a single object (ID 01115), we rely on the Hubble Space Telescope (HST) WFC3 photometry (Stefanon et al. 2017), as detailed in Isobe et al. (2023). Briefly, we derive the UV absolute magnitudes ( M UV ) from the HST/WFC3 photometry, use the empirical relationships between stellar mass and M UV at z = 4 -8 (Song et al. 2016) to estimate the  \nstellar masses, and rescale them to the Chabrier (2003) IMF. We exclude the single object with ID 01115 from the following analysis as its HST photometry is not given. For those not detected in the HST bands, we translate the corresponding 5 σ limiting magnitudes to the 5 σ upper-limits of stellar mass. While acknowledging that the empirical method utilizing M UV for estimating stellar mass may be less reliable compared to the method involving an SED fit to NIRCam photometry, we demonstrate in the subsequent sections (Sections 3.3 and 3.4) that our main findings remain unchanged even when excluding objects without NIRCam photometry. This is partly because the fraction of such objects is small ( ∼ 28%) in deriving the average relationships. Furthermore, we have found reasonable consistency in stellar mass estimates between the two methods by using objects with both NIRCam and HST photometry. Details are described in Appendix A. Note that we cannot correct for the slit-loss of the NIRSpec spectra for the objects without NIRCam photometry. This will not affect the metallicity measurements as we only use the line ratios. However, that does prevent us from deriving the equivalent widths of the optical emission lines such as EW(H β ), and the SFR based on the total H β luminosity.  \nIn summary, we construct a JWST sample of 182 objects at z = 3 . 8 -8 . 9. The redshift distribution of the sample is shown in red in Figure 3.", '3.1. Emission line ratios': "To examine the nebular properties of JWST objects at z = 4 -9 and compare with those of lower-redshift galaxies, we plot them in the line ratio diagnostic diagram of ([O iii ] λλ 5007 , 4959+[O ii ] λ 3727) / H β (R23) and [O iii ] λ 5007 / [O ii ] λ 3727 (O32) in Figure 4 (see also Mascia et al. 2023; Tang et al. 2023; Sanders et al. 2023). This diagram has been widely used to investigate the metallicity and ionization state in the local universe and up to z = 2 -3 (e.g., Kewley & Dopita 2002; Maiolino et al. 2008; Nakajima & Ouchi 2014; Troncoso et al. 2014; Sanders et al. 2016; Onodera et al. 2016; Strom et al. 2017; Nakajima et al. 2020). When compared to the local Sloan Digital Sky Survey (SDSS) galaxies, which are free from AGN contamination as determined by the BPT diagram (e.g., Baldwin et al. 1981; Kauffmann et al. 2003; see below for more details), we observe that the z = 2 -3 galaxies exhibit a higher O32 line ratio, indicating an overall increase in ionization parameter at higher redshifts. Figure 4 confirms the trend, and the z = 4 -9 JWST objects show a fur-  \nFigure 5. The [N ii ] BPT diagram. Red symbols represent the JWST objects, as shown in Figure 4, whose H α +[N ii ] lines are covered in the NIRSpec spectra. Arrows indicate 3 σ limits. Two popular demarcation curves between AGNs and star-forming galaxies, from Kewley et al. (2001) (dotlong dashed; 'Kew01') and Kauffmann et al. (2003) (short dash-long dashed; 'Kau03'), are also shown. None of the objects fall above the demarcation curves significantly beyond the measurement uncertainties.  \n<!-- image -->  \nigh O32 line ratio on average than typically seen in z = 2 -3 continuum-selected galaxies. Remarkably, one of the ERO objects, ERO 04590 at z = 8 . 5 presents a stringent lower-limit of O32 as > 14 . 8 (3 σ ), which corresponds to an ionization parameter of log U > -2 . 21 according to the prescription of Kewley & Dopita (2002) (see also Fujimoto et al. 2022 6 ). The lower-limit is a factor of ∼ 10 higher value than typically seen in the SDSS galaxies.  \nAt z = 2 -3, low-mass Ly α emitting galaxies (LAEs) are suggested to have a remarkably high O32 line ratio, which would be achieved by a hard ionizing radiation field and/or a low metallicity (Nakajima et al. 2016; Trainor et al. 2016; Erb et al. 2016). Its implication regarding the escape of ionizing photons (and that of Ly α photons) are also discussed (e.g., Nakajima & Ouchi 2014; Izotov et al. 2016, 2018; Verhamme et al. 2017; Steidel et al. 2018; Fletcher et al. 2019; Erb et al. 2019; Nakajima et al. 2020; Katz et al. 2020; Flury et al. 2022). The extreme O32 line ratios seen in the z = 4 -9 JWST objects are comparable to those in z = 2 -3 lowmass LAEs and z ∼ 0 Lyman-continuum (LyC) leaking objects, suggesting that high-redshift sources have a  \nnebular and ionization condition of gas that is similar to those typically found in lower-redshift low-mass galaxies having a strong Ly α and/or LyC escape with a hard ionizing spectrum. A similar finding is also discussed by Mascia et al. (2023) and Sanders et al. (2023), and notably by Tang et al. (2023), where the strong Ly α emitting galaxies at z > 7 are directly examined with JWST and suggested to present the largest O32 line ratios.  \nOne uncertainty in interpreting O32 arises as it also depends on the metallicity. In the next section, we present our method to derive the metallicities for the JWST objects using several indicators including R23, particularly by referring to the objects whose [O iii ] λ 4363 is detected and metallicity is precisely determined with electron temperature.  \nBefore moving to the metallicity results, we also examine our objects using another popular emission line diagram, the [N ii ] BPT diagram, which plots the ratios of [O iii ] λ 5007 / H β and [N ii ] λ 6584 / H α (Baldwin et al. 1981). This diagram is commonly used to diagnose the presence of an active galactic nucleus (AGN) and is applicable to our JWST objects up to a redshift of z < 6 . 9, where the H α +[N ii ] emission is covered by NIRSpec spectra. Figure 5 shows the [N ii ] BPT diagram for our JWST objects at z < 6 . 9, along with two curves that discriminate between sources dominated by AGNs and stars (Kewley et al. 2001; Kauffmann et al. 2003). Sources located above the curves are classified as AGNdominated.  \nFigure 5 shows that most JWST sources have an upper limit of [N ii ]. The [N ii ]-detected objects fall either below or on the demarcation curves within the measurement uncertainties. None of the JWST objects are thus classified as obvious AGNs. However, it is worth noting that low-metallicity AGNs with Z /lessorsimilar 0 . 5Z /circledot may contaminate the star-forming galaxy region on the diagram (e.g., Kewley et al. 2013; Nakajima & Maiolino 2022). We cannot fully rule out the presence of such low-metallicity AGNs in our sample based on the [N ii ] BPT result, but we suggest an absence of evolved AGNs as far as we explore ( z < 6 . 9). Accordingly, we do not remove any objects from the sample for the following results based on the [N ii ] BPT diagram. Although there are some objects whose [N ii ] upper-limits are too weak to conclude their BPT diagnostics, we mention in the following sections (Sections 3.3 and 3.4) that our main results are not changed by excluding these unclear objects.  \nAs a complementary approach to our method using the [N ii ] BPT diagram, a companion study by Harikane et al. (2023) conducts a search for faint AGNs  \nFigure 6. Comparison of literature determinations of electron temperature ( T e ) for four ERO galaxies (04590 in red, 06355 in orange, 10612 in green, 05144 in emerald green) where [O iii ] λ 4363 is identified. Different symbols indicate different publications, as shown in the legend. For the measurements of Arellano-C'ordova et al. (2022), we show the average values based on the two spectra (i.e., two observing blocks: o007 and o008) for 06355 and 10612, but refer only to the o008 value for 04590. For each object, the symbols are slightly offset in the horizontal direction for display purposes. Some of the measurements in the literature do not present uncertainties and are plotted without error bars in the vertical direction.  \n<!-- image -->  \nin our full JWST spectroscopic sample by examining the broad component (FWHM > 1000kms -1 ) around the H α emission line. The authors identify 10 objects with signatures of AGNs at redshifts z = 4 . 0 -6 . 9 (see also Kocevski et al. 2023; Ubler et al. 2023). To avoid potential biases in our mass-metallicity relations, we remove these 10 objects from our sample when we examine the mass-metallicity relations. These objects are flagged in Table D1.", '3.2. Gas-phase Metallicity': "3.2.1. Direct T e method  \nTen objects (4 from ERO, 4 from GLASS, and 2 from CEERS) present a significant detection of [O iii ] λ 4363 allowing us to reliably determine the oxygen abundance with the direct T e method at z = 4 . 0 -8 . 5 and use it as a proxy for gas-phase metallicity. We follow the procedure as summarized in Nakajima et al. (2022) to derive T e -based metallicities. Briefly, we first estimate the electron temperature of O 2+ zone ( T e ([O iii ])) using  \nTable 1. Summary of JWST Objects with Direct T e Method  \nNote -( /star ) Observed flux ratios relative to H β ( × 100). Upper-limit values at the 1 σ level. ( † ) The presence of high ionization line of [Ne iv ] is seen (Brinchmann 2022). ( ‡ ) The presence of broad H α is indicated (Harikane et al. 2023)  \nthe reddening-corrected [O iii ] λλ 4363 / 5007 ratio and an assumed electron number density of 100cm -3 with the task getTemDen from PyNeb (Luridiana et al. 2015). Different densities such as 10 -1000cm -3 do not significantly change the results (see Isobe et al. 2023 for the typical electron density in high-redshift galaxies; see also Fujimoto et al. 2022). The temperature of O + zone, T e ([O ii ]), is extrapolated from T e ([O iii ]) employing the prescription of Izotov et al. (2006) (cf., Brinchmann 2022 for a different assumption). We derive the abundance of O 2+ / H + with the fluxes of [O iii ] λλ 4959 , 5007 to H β and T e ([O iii ]), and O + / H + with the [O ii ] to H β reddening-corrected ratio and T e ([O ii ]) using the PyNeb package getIonAbundance . We do not take into account a higher ionization abundance of O 3+ / H + , as we follow the approximation given by Izotov et al. (2006) and confirm no clear presence of the He ii λ 4686 emission line in the spectra. In the spectrum of ERO 04590, the [O ii ] doublet is not detected at the 3 σ level. We therefore count the O 2+ component alone for its oxygen abundance. We have checked that the 3 σ upper-limit of [O ii ] would increase the O/H value by 0 . 05dex, which  \nis small enough with respect to the measurement error ( /similarequal 0 . 15dex). The measured oxygen abundances as well as T e ([O iii ]) are summarized in Table 1.  \nFor the ERO objects, there exist several earlier studies that report the metallicities with the electron temperatures (Schaerer et al. 2022; Curti et al. 2023a; Trump et al. 2022; Rhoads et al. 2023; Arellano-C'ordova et al. 2022; Brinchmann 2022). It is practically useful to compare our measurements with them. Figure 6 compares the T e ([O iii ]) values based on our measurements with the earlier studies for the four [O iii ] λ 4363-detected ERO objects. Different colors correspond to different objects. In some previous studies, the z = 8 . 5 object (ERO 04590 ; red-marks in the plot) is claimed to present an exceptionally high electron temperature at T e > 2 . 5 × 10 4 K, which is not usually observed in the local universe and may require some additional explanations (e.g., Katz et al. 2023; Rhoads et al. 2023). Our re-measurement confirms a high, but not extremely high, electron temperature for ERO 04590, T e = (2 . 08 ± 0 . 26) × 10 4 K, which is explainable by heating of young massive stars without  \nFigure 7. Comparison of metallicity relationships with strong line ratios (R23, R3, R2, O32, and Ne3O2 from top left to bottom right) for ERO (circles), GLASS (pentagons), and CEERS (diamonds) galaxies with metallicity measurements based on the direct T e method. Each data point is color-coded by its EW(H β ) value, if available, as shown in the legend. The black dashed, long-dashed, and dot-dashed curves represent the metallicity relationships from Maiolino et al. (2008), Curti et al. (2017, 2020), and Bian et al. (2018), respectively. The blue and orange curves show the functions derived in Nakajima et al. (2022) for lowmetallicity galaxies, with blue representing high-ionization galaxies (EW(H β ) > 200 ˚ A) and orange representing low-ionization galaxies (EW(H β ) < 100 ˚ A). Each curve is shown only within the metallicity range explored in the original paper.  \n<!-- image -->  \nFigure 8. Empirical metallicity indicator using both R23 and O32 as proposed by Izotov et al. (2019, 2021) in the low-metallicity regime (solid) and its extrapolation (dotted). The red symbols show the JWST objects as in Figure 7.  \n<!-- image -->  \nany special mechanisms. It is thus important to carefully reduce and combine the spectra, as well as to extract the 1D to reliably discuss the nebular properties with the measurements of faint emission lines. For the other objects, their electron temperatures are modest, T e = (1 . 1 -2 . 3) × 10 4 K, fairly consistent with the values in the earlier literature and similar to the electron temperatures seen in lower-redshift star-forming galaxies. We will also discuss comparisons of metallicity determinations later on the mass-metallicity relationship.  \nThe four ERO objects as well as the 6 newly identified GLASS and CEERS objects with [O iii ] λ 4363 provide the opportunity to test the empirical metallicity indicators at such high-redshift of z > 4 (Curti et al. 2023a). In Figure 7, we examine the following five popular metallicity indicators: [O iii ] λ 5007 / H β (R3), [O ii ] λ 3727 / H β (R2), [Ne iii ] λ 3869 / [O ii ] λ 3727 (Ne3O2), as well as R23 and O32 by comparing with the empirical indicators pro-  \ny Maiolino et al. (2008) (see also Nagao et al. 2006), Curti et al. (2017, 2020), Bian et al. (2018), and Nakajima et al. (2022). We draw each empirical relationship only in the calibrated metallicity range without showing any extrapolation. For the relationships of Nakajima et al. (2022), two curves are depicted in each panel, showing the dependence of EW(H β ) on the indicators. The authors use EW(H β ) as a proxy to correct for the degree of ionization state of the gas in each galaxy, because EW(H β ) is sensitive to the current massive-star formation efficiency and known to be well-correlated with the ionization state as probed by e.g., O32 (e.g., Nakajima & Ouchi 2014; Mingozzi et al. 2020; Nakajima et al. 2022). To examine the dependencies of EW(H β ) on the indicators, the JWST objects are color-coded by EW(H β ) except for CEERS 01536, whose EW(H β ) is not derived due to the lack of NIRCam images (Section 2.3).  \nThe JWST objects tend to present high ionization emission lines such as [O iii ] and [Ne iii ] stronger (and low ionization lines such as [O ii ] weaker) than expected from the local empirical relationships as defined by Maiolino et al. (2008) and Curti et al. (2017, 2020). Rather, they present a fairly good agreement with the relationships proposed by Bian et al. (2018) and Nakajima et al. (2022) for the large EW(H β ) objects. Because the relationships of Bian et al. (2018) are based on the highly-ionized objects as typically found at z ∼ 2 -3 on the [N ii ] BPT diagram (e.g., Baldwin et al. 1981; Steidel et al. 2014; Shapley et al. 2015), those consistencies suggest the JWST objects at z = 4 -8 . 5 are highly ionized systems. This trend is consistent with the analysis by Sanders et al. (2023). One caveat is that there is one JWST objects whose EW(H β ) is smaller than 100 ˚ A, GLASS 10021, falling close to the large EW(H β ) relationships of Nakajima et al. (2022) rather than the small EW ones like the other JWST objects despite the relatively small EW(H β ). This suggests we need to test whether the prescriptions work for low ionization sources as well with a larger sample at high-redshift.  \nIn summary, although we cannot fully test the indicators for low ionization sources, the current results suggest that the degree of ionization significantly influences the resulting metallicity value if one relies on an empirical metallicity indicator, and that no strong redshift evolution is seen in the strong line ratios as a function of metallicity at least for highly ionized galaxies. Ionization-corrections, such as those proposed in Nakajima et al. (2022), are thus crucial for metallicity estimations with the strong line methods. In particular, we find R23, R3, and R2 show a good agree-  \nment with the observations and the empirical relationships, and confirm that they are sensitive to a small change in metallicity. The indicators of O32 and Ne3O2 look highly dependent on the ionization state and their plateau-like behaviors against metallicity prevent us from deriving a stable metallicity solution at the metalpoor regime.  \nFigure 8 shows another metallicity indicator proposed by Izotov et al. (2019, 2021) which combines R23 and O32 to correct for the ionization state to improve the accuracy of metallicity in the low-metallicity regime. Unfortunately there is only one object, ERO 04590, whose metallicity is low enough to be fairly compared with the method, and whose non-detection of [O ii ] prevents us from confirming the reliability and accuracy of the method. This is exactly the situation Nakajima et al. (2022) have anticipated, i.e., not all of the lines are spectroscopically available at high redshift, and it is practically useful to use EW(H β ) as a probe of ionization state. Still, many of the JWST objects are located along the simple extrapolation of the relationship, indicating the method can be useful even up to 12+log(O / H) ∼ 7 . 8.  \nOne caveat should be noted for ERO 06355, which has an oxygen abundance of 12+log(O / H) = 8 . 3 and shows an unusually strong [O iii ] line. Curti et al. (2023a) have used the [O iii ] / H β vs. [O ii ] / H β diagram and found that this object is slightly above the limit that can be explained by star-forming galaxies alone. This indicates the possibility of a hidden high-energy ionizing source. In addition, Brinchmann (2022) tentatively suggest the presence of the high ionization line [Ne iv ] λλ 2422 , 2424, indicating a hard ionizing spectrum up to ∼ 60 eV, which cannot be fully explained by a conventional stellar population. These pieces of evidence suggest that the anomalous nature of this object may be explained by the presence of a high-energy ionizing source in the system, which may account for its deviation from the empirical metallicity relationships. Another note is for GLASS 160133 and GLASS 150029 that have a broad component associated with their H α emission line, indicating the presence of faint AGNs in the systems. However, they still present metallicity diagnostic line ratios and EW(H β ) that agree reasonably well with the relationships for star-forming galaxies, implying a minor contribution of AGN to the total emission strengths as well as to the optical continuum for these two objects. In any case, one would need higher ionization lines (e.g., [Ne v ]; Cleri et al. 2023) to further discuss the presence of hard radiation components in these systems.  \nFor the other JWST objects without a T e -based metallicity ( N o = 182 -10 = 172), we adopt the empirical metallicity indicators which have been tested and confirmed to work at high-redshift (Section 3.2.1) to estimate metallicities. Prior to estimating metallicities empirically, we require that the objects have coverage of the [O iii ]+H β emission lines in their spectra, and that the [O iii ] λλ 5007 / 4959 doublet line ratio is consistent with the theoretical value (2 . 98; Storey & Zeippen 2000) at a 2 . 5 σ significance level. We only work with objects that meet these criteria and have appropriate flux measurements. As a result, 11 objects were removed from the subsequent analysis.  \nFor the remaining 161 objects, we use the R23-index as a primary indicator, following the overall consistency of the metallicity indicator at high-redshift as shown in Section 3.2.1. One caveat in using R23 is that two solutions are derived for a given R23 value. We use the prescription of Nakajima et al. (2022) with EW(H β ) for the low-metallicity solutions (12 + log(O / H) /lessorsimilar 8 . 0), and the average relationship (without any ionizationcorrection) for the high-metallicity ones ( /greaterorsimilar 8 . 0). The latter is almost consistent with the indicators presented in Curti et al. (2017, 2020). If the two metallicity solutions are not significantly separated, i.e., the observed R23 spans the peak value that appears around 12 + log(O / H) = 8 . 0 within the uncertainty, we adopt the 1 σ lower-limit of the low-metallicity solution as the lowerlimit, and the upper-limit of the high-metallicity solution as the upper-limit. If there are two distinct metallicity solutions, we rely on the O32 line ratios to distinguish between them. It should be noted that O32 is used only for the purpose of distinguishing between the two branches, and not for calculating the actual metallicity value. Specifically, we choose the R23-based metallicity solution whose expected O32 line ratio according to the O32-metallicity indicator of Nakajima et al. (2022) is more close to the observed O32 than the other. If the [O ii ] is not detected and the resulting lower-limit on O32 is higher than the expected O32 at the highmetallicity branch at the > 3 σ level, we choose the lowmetallicity solution. If the lower-limit on O32 is not high enough to distinguish, we translate EW(H β ) to O32 using the average relationship found in Nakajima et al. (2022): log EW(H β ) = 0 . 64 × log O32+1 . 68, and choose the solution. Finally, the systematic uncertainty of the empirical method is added in quadrature to the metallicity errors. For the other cases, we cannot reliably distinguish between the two solutions and thus leave its metallicity unconstrained based on R23.  \nWe are unable to derive a R23-based metallicity for objects that lack [O ii ] coverage and/or appropriate  \ndust correction due to the absence of multiple Balmer emission lines. For such objects, we estimate metallicity using the R3 indicator, as it does not require a reddening correction. We follow a similar procedure as for R23 to account for the two-branch nature of the R3-index and estimate metallicities. As O32 is not applicable for these objects, we rely on EW(H β ) to correct for the ionization state and attempt to differentiate between the two solutions.  \nIn summary, our sample consists of 86 objects with an R23-based metallicity, 49 objects with an R3-based metallicity, and 26 objects without a metallicity constraint. We note that for the objects in the CEERS sample that have multiple spectra taken with different gratings and/or pointings, we estimate the metallicity for each spectrum and calculate the average value and its uncertainty by combining the individual metallicity measurements. While we understand that some of the spectra (P11 and P12) were obtained with different position angles compared to the others, we assume that the objects are compact and the different observations targeting the same objects at the center of the slit probe the same gas properties. After excluding 10 objects with spectroscopic signatures of AGNs (Section 3.1), comprising of 8 objects with R23-based metallicity and 2 objects with T e -based metallicity, we use 127 objects with an R23- or R3-based metallicity, along with 8 objects with direct T e method, for the subsequent analysis and discussion of mass-metallicity relationships.", '3.2.3. Other high-redshift galaxies from the literature': "In addition to our JWST sample of 135 objects constructed with the ERO, GLASS, and CEERS programs, we have compiled reports of metallicity at highredshift from the literature to discuss with a larger sample. Our compilation includes three objects at z = 8 . 1 -9 . 5 identified with the two DDT programs using the prism grating of NIRSpec (Proposal IDs: DD-2756 and DD-2767; (Williams et al. 2022; Heintz et al. 2022; Langeroodi et al. 2022; Wang et al. 2022)), one stacked point of 117 EIGER objects at z = 5 . 3 -6 . 9 whose spectra are taken with the NIRCam slitless spectroscopy (Proposal ID: 1243; Matthee et al. 2022; Kashino et al. 2022), four objects at z = 6 . 1 -6 . 4 with commissioning data of NIRCam slitless spectroscopy (Sun et al. 2022a,b), and four ALMA objects at z = 7 . 2 -9 . 1 whose metallicity is measured with [O iii ]88 µ m (Jones et al. 2020). Their redshift distributions are illustrated in blue in Figure 3. Note that the EIGER stacked object is treated as a single object in the histogram as well as in the following figures.  \nFor a fair comparison of metallicity in the following analysis, it is important to have metallicities that are estimated in a consistent manner, as done for our sample. All of the measurements compiled above are based on empirical indicators using strong emission lines that are calibrated with the direct T e method. For the NIRSpec DDT objects, five objects at z = 7 . 9 -9 . 5 have reported metallicity measurements in several studies (Williams et al. 2022; Heintz et al. 2022; Langeroodi et al. 2022; Wang et al. 2022). Despite initially adopting the O32-index as the metallicity indicator in the original manuscript, Heintz et al. (2022) have updated their results, now primarily using the R3-index of Nakajima et al. (2022) with correction for the ionization state for large EW(H β ) ( ≥ 200 ˚ A) objects. This approach is supported by Figure 7 in this paper and has been adopted for some of our JWST objects. It should be noted, however, that it would be more appropriate to adopt the EW(H β )-dependence of the indicator following the variation of EW(H β ) as observed in the DDT sample, instead of fixing the relation for large EW(H β ) (Heintz et al. 2022). Large EW(H β ) of ∼ 180 -250 ˚ A are confirmed for two of the five objects (Williams et al. 2022; Wang et al. 2022). Among the five objects reported in Heintz et al. (2022), we use three objects, RXJ-z9500, RXJ-z8152, and RXJ-z8149, in this paper. Their metallicities are estimated to be 12+log(O / H) = 7.56 (+0.16/-0.17), 7.68 (+0.18/-0.19), and 7.29 (+0.22/-0.28), respectively. One of the objects not included in our compilation, Abell-z7878, was originally suggested to have an [O iii ] λλ 5007 / 4959 doublet ratio that is significantly smaller than the theoretical value. We remove this object from our compilation following the same approach used for our JWST sample. The remaining object, Abell-z7885, still relies on the O32-index for estimating its metallicity. We have decided not to include it in our compilation due to several uncertainties associated with using the O32-index for metallicity estimations, such as the strong dependence of O32 on ionization state (Figure 7) and the accuracy of dust reddening correction, as we were careful regarding the results presented in the original manuscript of Heintz et al. (2022).  \nThe metallicity values for RXJ-z9500, RXJ-z8152, and RXJ-z8149 in Heintz et al. (2022) are consistent with those reported in other studies. The metallicity of RXJz9500 is first reported by Williams et al. (2022), where it was estimated to be 12+log(O / H) = 7 . 48 ± 0 . 08 using the R23+O32 method of Izotov et al. (2019, 2021) (Figure 8). An ionization correction is thus taken into account for the metallicity measurement there. Similarly,  \nFigure 9. Comparison between the mass-metallicity relation obtained in this study for ERO objects at z = 6 . 3 -8 . 5 using T e -based metallicities (colored circles) and those reported in the literature. Different symbol colors indicate different objects, as shown in Figure 6, and different symbols represent different publications, as indicated in the legend. For the metallicities from Arellano-C'ordova et al. (2022), we show the average values based on the two spectra (i.e., two observing blocks: o007 and o008) for 06355 and 10612, but refer only to the o008 value for 04590. For the relationships by Schaerer et al. (2022) and Curti et al. (2023a), the stellar masses are adopted as originally estimated in each paper. For the others where stellar masses are referenced from other studies, we show the two stellar masses from Carnall et al. (2023) and Tacchella et al. (2022) for each object (only Carnall et al. (2023) for 05144). The comparison confirms a good agreement between our measurements and earlier results on the MZ relation.  \n<!-- image -->  \nLangeroodi et al. (2022) estimate metallicities for RXJz8152 and RXJ-z8149 (RX2129-ID11002 and RX2129ID11022 in Langeroodi et al.) using the R23+O32 method (Izotov et al. 2019, 2021), and report 12 + log(O / H) = 7 . 65 ± 0 . 07 and < 7 . 51 (1 σ ), respectively. Wang et al. (2022) independently suggest a consistent metallicity for RXJ-z8152, although the ionization correction is not fully taken into account. We use the properties, including metallicities, that are calculated by averaging the values from Williams et al. (2022), Heintz et al. (2022), and Langeroodi et al. (2022) for the three DDT objects in our compilation.  \nFor the other compilation, we adopt the metallicity values as derived in the original papers. The stacked EIGER object has already derived a metallicity fairly consistent with the value based on our method using R23. The four NIRCam objects are all located in the high-metallicity branch based on the [N ii ] λ 6564 detection, and the indicator of Bian et al. (2018) is used. We also note that the stellar masses and SFRs are all corrected for different IMFs to have the same  \nFigure 10 illustrates the relation between stellar mass and metallicity (MZ) with the full sample at z = 4 -10. The MZ relations determined at z = 0, /similarequal 2 . 3, and /similarequal 3 . 3 with the direct T e method are also plotted  \n<!-- image -->  \nFigure 10. The relationship between stellar mass and metallicity. The red points represent galaxies at z = 4-10, including the filled circles, pentagons, and diamonds, which represent the ERO, GLASS, and CEERS objects analyzed in this paper. Red circles denote galaxies whose metallicity is determined using the direct T e method. The large red stars represent the average relationships for the ERO, GLASS, and CEERS objects in three equally separated mass ranges ( M /star =10 7 -10 8 , 10 8 -10 9 , and 10 9 -10 10 M /circledot ), along with the best-fit function shown as a thick long-dashed line (Equation 1). Emerald green open pentagons show the average relations at z = 3 -10 based on the JADES+CEERS sample (Curti et al. 2023b). Other red symbols represent high-redshift objects compiled from the literature, including open right-pointing triangles for NIRSpec DDT objects (Williams et al. 2022; Heintz et al. 2022; Langeroodi et al. 2022), left-pointing triangles for four NIRCam objects (Sun et al. 2022a,b), open hexagram for the stacked EIGER object (Matthee et al. 2022), and crosses for ALMA objects (Jones et al. 2020). Additionally, relationships at lower redshifts are displayed, including SDSS stacked galaxies in the local universe shown as a gray dashed curve (Andrews & Martini 2013), and those at z ∼ 2 . 3 with gray open triangles and ∼ 3 . 3 with gray open squares (best-fits shown as gray dot-dashed curve; Sanders et al. 2021). The curves are displayed in the mass ranges explored in the original papers. These low-redshift metallicities are based on the direct T e method.  \n<!-- image -->  \nChabrier (2003) IMF using the conversion factors shown in Madau & Dickinson (2014).  \nIn summary, our total sample consists of 147 galaxies (135 from our reduction of ERO, GLASS, and CEERS, and 12 from the compilation) with metallicity measurements at z = 4 -10. This sample is used in the following analysis of the metallicity relationships.", '3.3. Mass-Metallicity Relation': 'The main focus of this paper is to discuss the evolution of metallicity in star-forming galaxies. In this section, we specifically present the stellar mass-metallicity relation at high redshift using the JWST observations, as well as data compiled from relevant literature.  \nBefore presenting the full results we firstly focus on the four ERO objects and compare in Figure 9 our mass-metallicity relation with those presented in the earlier studies (Schaerer et al. 2022;', 'Curti et al. 2023a; Trump et al. 2022; Rhoads et al.': "2023; Arellano-C'ordova et al. 2022; Brinchmann 2022). The metallicities are all derived based on the direct T e method as summarized in Figure 6. The stellar masses are corrected for different IMFs to have the same Chabrier (2003) IMF. Schaerer et al. (2022) and Curti et al. (2023a) derive the masses in their own papers, while the others refer to either Carnall et al. (2023) or Tacchella et al. (2022) and hence both values are shown in the plot. We confirm our mass-metallicity relations for the four ERO objects overall show a good agreement with the earlier results. One note is that we confirm the metallicity of ERO 04590 is not extremely low, 12 + log(O / H) = 7 . 26 (+0.15/-0.13) for its stellar mass.  \nTable 2. Average values of stellar mass, metallicity and SFR for our JWST sample  \nNote -The first half of this table displays the average masses and metallicities of the samples used for the MZ relations, while the second half shows the samples used for the SFR-MZ relations. The second half samples exclude objects whose SFRs are not measured or constrained (Section 3.4), resulting in a smaller sample size compared to the first half.  \n(Andrews & Martini 2013; Sanders et al. 2021; see also Nishigaki et al. 2023 that the T e -based MZ relation at z = 0 continues decreasing down to M /star ∼ 10 5 M /circledot ). The red symbols denote the galaxies at z = 4 -10, including our 135 ERO, GLASS, and CEERS galaxies, 3 NIRSpec DDT objects, as well as the massive galaxies provided by NIRCam slitless and ALMA spectroscopy, and the low-mass stacked object as presented in Section 3.2.3. We note again that the 10 objects with spectroscopic signatures of AGNs (Harikane et al. 2023) have been excluded here to conservatively discuss the results free from any AGN biases (Section 3.1). The objects analyzed using the direct T e method are marked with a red open circle, suggesting a positive correlation between mass and metallicity in place in the redshift range z = 5 -8 . 5. Furthermore, we divide our full sample of ERO, GLASS, and CEERS galaxies into three according to the stellar mass: M /star = 10 7 -10 8 , 10 8 -10 9 , and 10 9 -10 10 M /circledot , and obtain the average MZ relations as  \nshown with the large open stars in Figure 10 and as given in Table 2. Note that the compiled objects, as well as the CEERS objects with only an upper-limit on M /star , are not used in deriving the average relations. Following the single power law form of Sanders et al. (2021), the average MZ relation of the z = 4 -10 galaxies can be approximated as:  \n12 + log(O / H) = Z 10 + γ log( M /star / 10 10 M /circledot ) (1)  \nwith the best-fit parameters of Z 10 = 8 . 24 ± 0 . 05 and γ = 0 . 25 ± 0 . 03 in the mass range of M /star ∼ 10 7 . 5 -10 9 . 5 M /circledot , as shown with a gray long-dashed line. The parameter γ corresponds to the slope of the MZ relation, and its best-fit value and uncertainty confirm an increasing trend of metallicity with stellar mass, as tentatively seen with the three ERO objects (e.g., Schaerer et al. 2022; Curti et al. 2023a; Trump et al. 2022), and as widely known at low-redshift.  \nCompared to the z = 0 MZ relation, these highredshift galaxies clearly present a metallicity lower than typical galaxies at z = 0 for a given stellar mass. The decrease is typically ∼ 0 . 5dex around M /star ∼ 10 9 M /circledot , but it becomes smaller at the low-mass end ( ∼ 0 . 3dex). Interestingly, a similar offset of ∼ 0 . 3dex is observed between z = 0 and z ∼ 2 -3, suggesting that the evolution of MZ relation is small from z ∼ 2 -3 to z = 4 -10. Although there may be a decrease of ∼ 0 . 2dex in the typical metallicity at the high-mass end of M /star ∼ 10 9 . 5 M /circledot from z ∼ 2 -3 to z = 4 -10, no strong evolution is found beyond the error. The same conclusion can be drawn from comparisons with the MZ relations at z /lessorsimilar 4 whose metallicities are empirically estimated with the strong line indicators, as presented in Appendix B.  \nFigure 10 also includes a comparison with the latest results from the JADES observations. Curti et al. (2023b) very recently report on the metallicity measurements of galaxies at z = 3 -10 based on deep JADES spectroscopic observations (see also Cameron et al. 2023b). In the figure, their MZ relation is shown using emerald green open pentagons, which represent the combination of the low-mass JADES sample with the high-mass CEERS sample, latter of which is provided in this paper. We find that their MZ relation, despite covering slightly different redshift ranges, agrees well with ours (Eq. 1). It is not surprising to see a good agreement between our results and those reported by Curti et al. (2023b) in the high-mass end ( M /star /greaterorsimilar 10 8 M /circledot ), since the regime is dominated by the CEERS objects presented in this paper. The JADES objects confirm that the MZ relation we obtained in Eq. 1 continues down to M /star /similarequal 10 7 M /circledot on average.  \nFigure 11. The MZ relation in three different redshift bins: (a) z = 4-6, (b) z = 6-8, and (c) z = 8-10. The red symbols and the best-fit function (Eq. 1) are as shown in Figure 10. The large stars represent the average MZ relations re-derived in each redshift bin, by splitting the sample into two groups based on stellar mass to have the equal numbers of galaxies. For the JADES+CEERS relations (Curti et al. 2023b in emerald green), we plot their z = 3 -6 sub-sample's relations in Panel (a), while adopt the z = 6 -10 relations in Panels (b) and (c). In addition, the cosmological simulation results at z = 5, 7, and 9 are displayed in Panel (a), (b), and (c), respectively; FIRE in black (Ma et al. 2016), IllustrisTNG in green (Torrey et al. 2019), FirstLight in blue (Langan et al. 2020 in the low-mass regime with solid curves, and Nakazato et al. 2023 in the high-mass regime with shades), Astraeus in yellow (Ucci et al. 2021), and FLARES in magenta (Wilkins et al. 2022). Some extrapolations are applied, as detailed in the text.  \n<!-- image -->  \nIn Appendix A, we investigate our MZ relation using only galaxies for which the stellar mass is welldetermined through SED fitting to the JWST/NIRCam photometry, after excluding CEERS objects whose masses are estimated empirically from M UV (as discussed in Section 2.3). Because we have already excluded the CEERS objects with only an upper-limit on M /star , all of which are M UV -based, in deriving the average relations, the fraction of objects without NIRCam photometry is small ( ∼ 28%). Indeed, our conclusions remain unchanged when considering only the objects with reliable measures of M /star , as demonstrated in the figures in Appendix A. Moreover, we mention in the [N ii ] BPT diagram in Section 3.1 that some of the JWST objects have [N ii ] upper-limits that are too weak to conclude the ionization nature (stars vs. AGNs). By removing these unclear objects and using only galaxies that are surely diagnosed as star-forming galaxies in the [N ii ] BPT diagram (i.e., those with an (upper-limit on) [N ii ]/H α falling below the demarcation curve), we obtain a fully consistent MZ relation as found in the full sample (Figure 10). This implies a negligible contribution of AGNs in our sample.  \nTo further examine any redshift evolution among our sample from z = 4 to 10, we plot in Figure 11 the MZ relations for the three different redshift bins, z = 4 -6, 6 -8, and 8 -10. In each panel, the sub-sample is further divided into two groups based on their masses, and their average values are shown with large stars as summarized in Table 2. Although there are large error bars, these average points suggest that the slope and normalization of the MZ relation in different redshift bins are consistent with those based on the full sample at z = 4 -10 (Eq. 1). Figure 11 indicates that there is no significant evolution in the MZ relation among our z = 4 -10 sample. We note that there may be a weak trend towards lower metallicity in the highest-redshift bin, albeit with a small sample size. Curti et al. (2023b) also tentatively suggest a similar evolution. That will be further explored in the discussion.  \nThe no/weak evolution is consistent with the predictions of some cosmological simulations, as presented and compared in Figure 11. We compile the hydrodynamic, N -body, and/or semi-numerical simulations showing the MZ relation at high-redshift; FIRE by Ma et al. (2016) in black, IllustrisTNG by Torrey et al. (2019) in green, FirstLight by Langan et al. (2020) and Nakazato et al. (2023) in blue (see also Ceverino et al. 2017), Astraeus by Ucci et al. (2021) in yellow, and FLARES by Wilkins et al. (2022) in magenta (see also Lovell et al. 2021). We plot the theoretical predictions  \nat z = 5, 7, and 9 for the redshift bin of z = 4 -6, 6 -8, and 8 -10, respectively. For the FIRE curves, we extrapolate the result at z = 6 to z = 7 and 9 using their redshift evolution function, although the evolution is tiny (∆ log (O / H) = -0 . 025dex from z = 6 to 7, and -0 . 05dex from z = 6 to 9). Similarly, the IllustrisTNG predictions at z = 5, z = 7, and z = 9 are extrapolated from the results at z = 4, z = 6, and z = 6, respectively, using their redshift evolution function (∆log(O / H) = -0 . 08dex from z = 4 to 5, -0 . 075dex from z = 6 to 7, and -0 . 1dex from z = 6 to 9). For the FirstLight results, we refer to Langan et al. (2020) and Nakazato et al. (2023) to prove the low-mass and the massive regime, respectively. We assume the z = 8 relation of Langan et al. (2020) in panel (c), and the z = 6 relation of Nakazato et al. (2023) in panel(a), assuming no strong evolution at z = 8 -9 and z = 5 -6, respectively. For the FLARES results, we adopt the stellar metallicities of only young ( < 10Myr) star particles and assume the stellar and gas-phase metallicities in the region of massive-star formation are comparable. Likewise, the metallicities for the IllustrisTNG and FirstLight model are SFR-weighted and massweighted of young ( < 100Myr) star particles, respectively, allowing a fair comparison with our metallicity measurements based on the nebular emission lines. On the other hand, we note that the FIRE model adopts the mass-weighted metallicity of all gas particle that belong to the ISM, and the Astraeus model counts the oxygen mass in the halo without any weighting which could result in an inequitable comparison with the observations.  \nComparing the observations with the simulations in Figure 11, we find the observed MZ relation and its weak evolution over z = 4 -10 are in good agreement particularly with IllustrisTNG , FirstLight , and FLARES in the mass range of M /star = 10 8 -10 9 . 5 M /circledot . Notably, the slope of the MZ-relation found in the IllustrisTNG results is likely coincide with the observations. On the other hand, the simulations of FIRE and Astraeus are generally suggested to under-predict metallicities except for some metal-poor galaxies such as ERO 04590 found in the highest-redshift bin at z > 8. That is probably due to their implementation of feedbacks, i.e., galaxies eject too many metals from galaxies and/or accreting gas is too efficient in lowering the metal contents in galaxies (Ma et al. 2016; Ucci et al. 2021). Finally, we note that the average metallicity measured for low-mass galaxies below M /star < 10 8 M /circledot can be slightly higher than the existing predictions by ∼ 0 . 2dex. It is thus necessary to increase the sample size to determine the MZ relation in the low-mass end as well as to extend the  \nFigure 12. The relationship between stellar mass and SFR, with SFR derived based on the total H β luminosity (as described in the text) for use in the SFR-MZ relation. The symbols used are consistent with those in Figure 10. The background lines represent sSFR values ranging from 10 -9 to 10 -6 yr -1 , from bottom to top. The JWST objects are distributed along the sequence of sSFR values around 10 -8 -10 -7 yr -1 .  \n<!-- image -->  \ntheoretical predictions such as IllustrisTNG towards lower-mass to conclude the consistency and to betterunderstand the early chemical enrichment in the lowmass systems at high-redshift.", '3.4. Mass-Metallicity-SFR Relation': "The next important aspect is the SFR-dependence of the MZ relation to discuss the chemical evolution, given the claim of a redshift-invariant fundamental relation between mass, metallicity, and SFR (SFR-MZ relation) out to z ∼ 2 -3 (Mannucci et al. 2010; Sanders et al. 2021). There are several expressions to describe the mutual dependencies between the three quantities (Mannucci et al. 2010; Lara-L'opez et al. 2010; Andrews & Martini 2013; Sanders et al. 2017; Curti et al. 2020; Sanders et al. 2021). In this paper, we primarily use the variable µ α originally proposed by Mannucci et al. (2010): µ α = log(M /star ) -α log(SFR), with α = 0 . 66 being adopted here that minimizes the scatter of the local low-metallicity galaxies with a direct T e -based metallicity in the µ α -metallicity plane (Andrews & Martini 2013) 7 :  \n12 + log(O / H) = 0 . 43 × µ 0 . 66 +4 . 58 . (2)  \nFigure 13. (a: Left:) Star-formation rate dependence of the Mass-Metallicity (SFR-MZ) relation. The symbols and curves are consistent with those in Figure 10, including the average points which are based on the three mass ranges that are equally separated. The coefficient of 0 . 66 for the combination of stellar mass and SFR is as indicated by Andrews & Martini (2013) in the local universe, down to M /star ∼ 10 7 . 4 M /circledot . Stacked galaxies at z = 2 -3 with T e -based metallicities are also suggested to fall on the z ∼ 0 relation in this parameter space (Sanders et al. 2021). Notably, the JWST objects analyzed in this study at z = 4 -10 (red) also fall on the same SFR-MZ relation on average. A consistent view is obtained with the JADES+CEERS average points at z = 3 -10 (emerald green; Curti et al. 2023b). (b: Right:) Metallicity difference from the SFR-MZ relation of Andrews & Martini (2013) (∆ log(O/H) = observed -predicted metallicity) for different redshift sources. The symbols are consistent with Panel (a), except for the average points (large red stars), which are recalculated on this panel to represent the average ∆ log(O/H) values for the objects found in three different redshift bins, z = 4 -6, 6 -8, and 8 -10. For the JADES+CEERS sample, average points of three mass bins are plotted for the two sub-samples at z = 3 -6 and 6 -10, as found in Figure 11. This plot illustrates that no evolution is seen up to z ∼ 8, but a significant decrease in metallicity is observed beyond z > 8, albeit with a small sample size in this highest-redshift bin.  \n<!-- image -->  \nThis is advantageous to be directly compared with a majority of the JWST objects down to M /star ∼ 10 7 . 4 M /circledot . Interestingly, the z /similarequal 2 . 3 and /similarequal 3 . 3 stacked objects of Sanders et al. (2021), whose metallicities are reliably determined with the direct T e method, fall directly on the same relationship. Another popular form of the SFR-MZ relation is provided by Curti et al. (2020). Our choice of the relation of Andrews & Martini (2013) against that of Curti et al. (2020) will be revisited later.  \nTo estimate SFRs for the JWST objects, we adopt the total, reddening-corrected H β luminosity for each object, as it is the best indicator for the on-going ( /lessorsimilar 10Myr) star-formation activity. The SFR relation of Kennicutt (1998) is adopted for the Balmer line with a correction of IMF to Chabrier (2003) using the conversion factor of (Madau & Dickinson 2014). Accordingly, the SFR-MZ relation is examined only for the JWST objects whose reddening correction is successfully applied and whose spectrum is slit-loss corrected. The latter constraint depends on whether the object has the NIRCam coverage (only in the CEERS field; Section 2.3). Among the CEERS objects lacking a slitloss correction, we rescue 42 objects whose UV stellar emission is constrained with HST. For these objects, we  \ntranslate M UV into SFR(H β ) assuming a typical ionizing photon production efficiency as indicated at highredshift (log ξ ion = 25 . 6 ± 0 . 2; e.g., Endsley et al. 2022; Matthee et al. 2022) and the factor of 0 . 44 difference between the reddening E(B -V) for the nebular and stellar emission (Calzetti et al. 2000). The assumption of ξ ion will be revisited elsewhere (K. Nakajima et al. in prep.) using the JWST sample presented in this paper. The objects lacking a proper dust reddening correction or having just upper-limits on both M /star and SFR are excluded. In short, 96 out of the 147 objects are used for the SFR-MZ relation at high-redshift. The second half of Table 2 presents the average SFR values for different masses and redshifts.  \nIn Figure 12, we show the individual and average distributions of the JWST objects on the stellar massSFR plane using the SFR values as derived above (i.e., mainly from H β ). The high-redshift objects are distributed along the sequence of sSFR = 10 -8 -10 -7 yr -1 . This is overall consistent with the starformation main sequence of galaxies at z = 4 -7 where sSFR of 10 -8 . 5 -10 -7 . 5 yr -1 is typically suggested (e.g., Stark et al. 2013; de Barros et al. 2014; Santini et al. 2017; Popesso et al. 2023). We note that some of our  \nFigure 13(a) show the SFR-MZ relation of the z = 4 -10 galaxies and their average points as shown in Figure 10 on the µ 0 . 66 -metallicity plane, together with the z = 0 average relation (Eq. 2) and the data-points of z /similarequal 2 . 3 and 3 . 3 stacked galaxies (Sanders et al. 2021). Adopting the best-fit parameter α = 0 . 66 found in the low-redshift metal-poor star-forming galaxies (Andrews & Martini 2013), the z = 4 -10 objects interestingly fall on the same SFR-MZ relationship. The JADES+CEERS sample, notably including lower-mass galaxies than ours, is also illustrated to show a consistent SFR-MZ relation on average.  \n<!-- image -->  \nFigure 14. (a: Left:) Similar to the left panel of Figure 13, but with the regions of µ 0 . 66 explored and unexplored by the original paper of Andrews & Martini (2013) highlighted. The blue solid curve represents the best-fit 2D projection of the SFR-MZ relation by Andrews & Martini (2013), determined down to µ 0 . 66 ∼ 7 . 5, with extrapolation towards lower mass shown by the dotted curve. Among the JWST objects, 83 % have µ 0 . 66 > 7 . 5 and can be directly compared with Andrews & Martini's relation. (b: Right:) Similar to the left panel, but with a different coefficient of 0 . 55 adopted for the combination of mass and SFR on the abscissa axis, as presented by Curti et al. (2020). The blue solid curve represents the best-fit 2D projection of Curti et al. (2020)'s SFR-MZ relation, determined down to µ 0 . 55 ∼ 8 . 5, and extrapolated towards lower mass with the dotted curve. Only approximately 15 % of the JWST objects have µ 0 . 55 > 8 . 5, and as such, the vast majority of the sample occupies a parameter space that is not explored by Curti et al. (2020).  \n<!-- image -->  \ngalaxies are above sSFR /greaterorsimilar 10 -7 . 5 yr -1 , in particular in the low-mass regime. This can be because the current spectroscopically-confirmed sample is partly biased toward actively star-forming systems, and/or because the sample contains higher-redshift objects at z > 7.  \nWe adopt the SFR-MZ relation of Andrews & Martini (2013) at z = 0 that is compared with the relation at z = 4 -10 because the sample of Andrews & Martini (2013) contains comparably low-mass, actively starforming galaxies. Figure 14(a) clarifies the parameter space that is explored by Andrews & Martini (2013), demonstrating that the majority (83 %) of the JWST objects have µ 0 . 66 > 7 . 5 and can be directly compared with the relation of Andrews & Martini (2013). Note that local extremely metal-poor galaxies with M /star  \n∼ 10 5 -10 7 M /circledot and 12 + log(O / H) = 7 . 0 -7 . 7 are suggested to be reproduced by the chemical evolution models (Lilly et al. 2013) with the same parameters as found for the Andrews & Martini's galaxies (Nishigaki et al. 2023). This evidence supports the simple extrapolation of the relation towards lower masses. On the other hand in Figure 14(b), we show another form of the SFR-MZ relation derived by Curti et al. (2020). The authors find another best-fit α = 0 . 55 as the best 2D projection of the SFR-MZ relation on the µ α -metallicity plane by using the global sample of stacked SDSS galaxies. The relation is determined over µ 0 . 55 ∼ 8 . 5 -11 . 0. Because ∼ 85% of the JWST sample have µ 0 . 55 below 8 . 5, most of the z = 4 -10 galaxies occupies the parameter space that is not explored by Curti et al. (2020). Moreover, Curti et al. (2020) clarify a different degree of SFR-dependence of the MZ relation for the low- and high-sSFR subsamples in the local universe. For the low-sSFR subsample (sSFR < 10 -9 . 5 yr -1 ), a weaker SFR dependence is indicated ( α = 0 . 22). On the other hand, the high-sSFR subsample with sSFR > 10 -9 . 5 yr -1 presents a stronger SFR dependence with α = 0 . 65, in close agreement with the best-fit value found by Andrews & Martini (2013). Given the high sSFRs for the JWST objects (sSFR ∼ 10 -7 -10 -8 yr -1 ; see Figure 12) and also for the z = 2 -3 galaxies (sSFR ∼ 10 -8 . 5 yr -1 , we conclude the strong SFR dependence  \nFigure 13(b) clarifies the redshift evolution on the SFR-MZ relation of Andrews & Martini (2013) by showing the residual metallicity for a given µ 0 . 66 with respect to the z = 0 relation: ∆ log(O / H) = 12 + log(O / H) obs -12+log(O / H)( µ 0 . 66 ; Eq. 2) for each galaxy's redshift. In this panel, the average points are derived for the objects found in the different redshift bins, at z = 4 -6, 6 -8, and 8 -10.  \n<!-- image -->  \n<!-- image -->  \nFigure 15. Metallicity difference from the SFR-MZ relation of Andrews & Martini (2013) as a function of stellar mass in three different redshift bins: (a) z = 4 -6, (b) z = 6 -8, and (c) z = 8 -10. The symbols are consistent with those in Figure 11. No clear dependence of ∆ log(O/H) on stellar mass is evident beyond the statistical errors.  \n<!-- image -->  \nof α = 0 . 66 as found by Andrews & Martini (2013) is  \npreferred and adopted in this paper to be compared with the high-redshift galaxies. We discuss the evolution of the Curti et al.'s SFR-MZ relation later in Section 4.  \nTwo interesting results arise in Figure 13(b). One is that the SFR-MZ relation shows no evolution up to z ∼ 8 within ∆ log O / H /similarequal 0 . 3dex. Another is that a significant decrease of metallicity is found beyond the error at z > 8. The decrease at the highest-redshift bin is not visible in Panel (a) due to the small sample size, but hinted by the MZ relation at z = 8 -10 (Figure 11c). To further examine the evolution of the SFR-MZ relation, we plot in Figure 15 the residual metallicity as a function of stellar mass for each of the three redshift bins. According to the two average points probing different mass ranges, no significant dependency of the residual metallicity on stellar mass is indicated in each redshift bin with the current sample. For the highest redshift bin, the decrease of metallicity is suggested for individual galaxies regardless of their stellar masses. Although the JADES+CEERS sub-sample at z = 6 -10 appears to exhibit a mass dependence in the residual metallicity, we suggest this would be caused by a bias introduced by the dominance of lower-redshift galaxies (i.e., z = 6 -8) in the lower-mass regime. This is consistent with the MZ plane in Figure 11, where the lowest-mass point of the JADES+CEERS sub-sample at z = 6 -10 is more consistent with our z = 6 -8 relation rather than the z = 8 -10 one. These results suggest that chemical properties of star-forming galaxies up to z ∼ 8 are very similar to those of local and z = 2 -3 counterparts, while a break of the metallicity equilibrium state may be in place beyond z ∼ 8. More discussions follow in Section 4.  \nFinally, we have conducted several tests to demonstrate the robustness of the results presented in this section. In Appendix A, we examine the SFR-MZ relation using only the JWST galaxies whose stellar mass is welldetermined with SED fit to the NIRCam photometry and whose SFR is obtained with the slit-loss corrected H β . Our conclusions remain unchanged even with the smaller sample that has good measurements of M /star and SFR, although the scatter of individual data points, as seen in Figure 15 for example, becomes less prominent. In Appendix C, we further investigate the SFR-MZ rela-  \nion by using only the JWST galaxies with µ 0 . 66 > 7 . 5, where the relation is explored by Andrews & Martini (2013), allowing for a more robust comparison with highredshift galaxies to discuss their evolution. We confirm that almost the same SFR-MZ relation and its redshift evolution are obtained based on this subsample, suggesting that the extrapolation of the relation at µ 0 . 66 ∼ 7 -7 . 5 for the low-mass end galaxies in our sample does not introduce biases. Furthermore, we note that a consistent view of the SFR-MZ relations is also supported by using only galaxies that are surely diagnosed as star-formation dominated in the [N ii ] BPT diagram (see also Section 3.3).", '4. DISCUSSION AND SUMMARY': "We have conducted a (re-)analysis of the JWST/NIRSpec ERO data, as well as the ERS data of GLASS and CEERS, to investigate the chemical properties of galaxies at redshifts ranging from z = 4 to 10. Our analysis includes the use of 135 JWST objects based on our improved reduction and calibration of the NIRSpec data, as well as the compilation of 12 objects from the other recent JWST observations and previous literature. We confirm that our new emissionline flux measurements and errors successfully address the issues reported in the literature related to Balmer decrements and electron temperatures. The estimated electron temperatures for the four ERO objects, along with the 6 GLASS + CEERS objects with [O iii ] λ 4363 at z = 4 -8 . 5 range from T e = 1 . 1 × 10 4 to 2 . 3 × 10 4 K. These temperatures are similar to those found in lowerredshift star-forming galaxies and can be fully explained by heating of young massive stars, without the need for additional ionizing mechanisms. This is particularly evident for the z = 8 . 5 object, whose T e has been updated from /greaterorsimilar 25000K to (2 . 08 ± 0 . 26) × 10 4 K.  \nWe have determined the mass-metallicity (MZ) relation for galaxies at z = 4 -10 and find no significant evolution compared to z ∼ 2 -3 when extrapolating to the low-mass regime. This result is consistent with the early JWST study results reported shortly after the ERO data release, except for the z = 8 . 5 object whose deviation from the others on the MZ relation becomes small after re-measurement of T e . Furthermore, our compilation of results does not show any significant evolution among the z = 4 -10 sample. Theoretical simulations also predict a similar trend of small redshift evolution in the MZ relation at z ∼ 5 -9 (∆log(O / H) /lessorsimilar 0 . 1dex; Ma et al. 2016; Langan et al. 2020; Ucci et al. 2021; Nakazato et al. 2023), with some simulations suggesting a potentially observable decrease in metallicity with redshift for a given mass ( ∼ 0 . 26dex  \nfrom z = 5 to 9; Torrey et al. 2019). Although there may be a weak evolution towards lower metallicity at z > 8, the current sample has relatively large uncertainties in metallicity (∆ log(O/H) ∼ 0 . 3dex), which makes it difficult to fully distinguish between different predictions of redshift evolution at high redshifts z > 4, especially for low-mass galaxies with M /star < 10 8 M /circledot . Further detailed observational and theoretical studies are needed to address the apparent discrepancy and understand the early chemical enrichment in low-mass systems at high redshifts.  \nWe have also investigated the SFR dependence of the MZ relation (SFR-MZ relation) at high redshifts, finding (i) no evolution from z = 0 to z = 4 -8 within ∆log(O / H) = 0 . 3dex, and (ii) a significant decrease of metallicity at z > 8. The former finding supports the idea that the SFR-MZ relation, also known as the Fundamental Metallicity Relation (FMR), can indeed describe the properties of galaxies with no redshift evolution. This suggests the existence of a metallicity equilibrium state via mechanisms such as star-formation, metal-poor gas inflow, and outflow (e.g., Lilly et al. 2013) that persists up to z ∼ 8. This sounds inconsistent with the potential deviation from the FMR at z /greaterorsimilar 2 . 5 as originally indicated by Mannucci et al. (2010). The authors have suggested, based on metallicities estimated using local empirical relations, that the same SFR-MZ relation shows no redshift evolution up to z ∼ 2 . 5. They claimed that galaxies at z ∼ 3 fall below the relation by approximately ∼ 0 . 6dex, with a combination of M /star and SFR of α = 0 . 32 in µ α . However, when using T e -based metallicities from Sanders et al. (2021) and adopting a different parameterization of α = 0 . 66 in µ α from Andrews & Martini (2013), no clear evolution from z = 0 up to z /similarequal 3 . 3 is found. Our findings, which utilize the parameterization of α = 0 . 66, support the claim of no average evolution of the SFRMZ relation up to z ∼ 8. This apparent inconsistency is likely attributed to the differences in metallicity estimations. High-redshift galaxies are known to have a higher ionization state of gas, parameterized by the ionization parameter, compared to local galaxies (Figure 4; see also e.g., Nakajima & Ouchi 2014; Sanders et al. 2023). Consequently, local empirical relations that assume a relatively low ionization parameter may result in biased metallicity estimates for high-redshift galaxies, particularly at z /greaterorsimilar 2 . 5 where the H α +[N ii ] λ 6584 lines are not available in the pre-JWST era. In fact, Figure 7 demonstrates that JWST objects with a large EW(H β ) exhibit a systematic offset from the empirical relations of Maiolino et al. (2008), which is used for metallicity estimations in Mannucci et al. (2010). Therefore, accu-  \nFigure 16. Similar to the right panel of Figure 13, but employing the SFR-MZ relation of Curti et al. (2020). As similarly found in Figure 13, our analysis (red) and the JADES+CEERS study by Curti et al. (2023b) (emerald green) suggest a consistent evolution of the SFR-MZ relation at z = 4 -10 on average. If we adopt the formalism of Curti et al. (2020) which is defined at z = 0, our JWST objects begin to deviate from the SFR-MZ relation at z ∼ 6. This finding contrasts with the use of another z = 0 SFR-MZ relation of Andrews & Martini (2013), which suggests deviations starting at z ∼ 8 (Figure 13b). It is worth noting, however, that a simple extrapolation of the relation towards lower mass and higher SFR is assumed for the majority of the JWST objects (see Figure 14(b)).  \n<!-- image -->  \nations using the reliable direct T e method, or accounting for the ionization state evolution as prescribed in Nakajima et al. (2022) and Izotov et al. (2019, 2021), are essential for future metallicity studies of the high-redshift universe with JWST.  \nThe latter finding of the potential evolution of the SFR-MZ relation beyond z ∼ 8 is intriguing. A similar result is supported by Heintz et al. (2022) as independently illustrated with the three DDT objects in our figures (red open right-pointing triangles), and also by the recent study of Curti et al. (2023b) based on the deep JADES observations. However, in contrast to our results using the SFR-MZ relation of Andrews & Martini (2013), Curti et al. (2023b) argue that galaxies begin to deviate at z ∼ 6, if they adopt the SFR-MZ relation of Curti et al. (2020). Because we confirm that our results and the JADES study yield the similar MZ and SFRMZ relations, as evidenced by the consistency between the red and emerald green symbols in Figures 10, 11, and 13, we attribute the redshift differences where deviations occur to the different functional forms of the  \nSFR-MZ relation defined at z = 0. This is also supported by the evolution of the SFR-MZ relation based on the Curti et al. (2020)'s formalism for our JWST objects, as shown in Figure 16, which indicates a consistent evolution with deviations starting at z ∼ 6 as reported by Curti et al. (2023b). While the use of a 3D SFR-MZ relation can account for the SFR-dependence of the MZ relation's slope and capture the metallicity variations more accurately than a simple 2D projection using µ α as assumed in the Andrews & Martini's relation (e.g., Curti et al. 2020), this is probably not an essential issue in this study. This is because the JWST objects analyzed so far present a narrow distribution of sSFR (Figure 12; see also Figure 4 of Curti et al. 2023b). We acknowledge that extended local-baseline studies that cover low-mass and high sSFR galaxies at z = 0 are necessary to discuss the chemical evolution, including the fundamental SFR-MZ relation, and determine when and if high-redshift galaxies begin to deviate from the relation. In this work, however, we regard that the best estimates are the results based on the SFR-MZ relation by Andrews & Martini (2013) as it provides a more appropriate comparison by covering the low-mass regime where most of the JWST objects are found (Section 3.4 and Figure 14).  \nIf we assume that the first star formation began to occur at z = 16 -27 (Harikane et al. 2022b; see also e.g., Abel et al. 2002; Bromm et al. 2002; Dayal & Ferrara 2018), galaxies at z = 8 -10 would be at most ∼ 200 -500Myr old. During such a short timescale after the Big Bang, there could be a higher probability for galaxies at higher redshifts to be in the early stages of formation, where they have not yet reached the metallicity equilibrium state via processes such as star formation, inflow, and outflow. An alternative explanation could be varying degrees of feedback processes, including inflow and outflow, in the early universe. As evident from the comparison of cosmological simulation results in Figure 11, galaxies could exhibit lower metallicities if metals are efficiently ejected from galaxies or diluted by gas accretion from inflow of pristine gas, as observed in simulations such as FIRE and Astraeus . Furthermore, other events such as galaxy mergers and AGN activities may also play a role in modulating star-formation, causing metal redistribution, and potentially influencing the SFR-MZ relationship at highredshift (e.g., Springel et al. 2005; Torrey et al. 2012; Weinberger et al. 2018). Due to the limited sample size and associated metallicity errors, particularly at z > 8, it is challenging to make definitive conclusions from the current study. Further statistical investigations of metallicity are crucial to confirm the suggested lack of  \nevolution in the SFR-MZ relation up to z ∼ 8, followed by a decrease at higher redshifts. These findings need to be rigorously compared with theoretical studies to better understand the underlying physics that govern early galaxy evolution.", 'ACKNOWLEDGEMENTS': 'We are grateful to Hidenobu Yajima, Hajime Fukushima, Chris Lovell, Yurina Nakazato, and the anonymous referee for useful comments and discussions that improved our manuscript. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1324 (ERSGLASS), #1345 (ERS-CEERS), #2561 (UNCOVER), and #2736 (ERO). The authors acknowledge the teams  \nof JWST commissioning, ERO, GLASS, UNCOVER, and CEERS for developing their observing programs with a zero-exclusive-access period. Moreover, this work is based in part on observations taken by the CANDELS Multi-Cycle Treasury Program with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. This paper is supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, as well as the joint research program of the Institute of Cosmic Ray Research (ICRR), the University of Tokyo. This work is supported by KAKENHI (JP19H00697, JP20H00180, and JP21H04467) Grant-in-Aid for Scientific Research (A) through the Japan Society for the Promotion of Science. In addition, KN acknowledges support from JSPS KAKENHI Grant JP20K22373. YI is supported by JSPS KAKENHI Grant JP21J20785, and also acknowledges funding from the Hayakawa Satio Fund awarded by the Astronomical Society of Japan. YH acknowledges support from JSPS KAKENHI Grant JP21K13953.', 'REFERENCES': "Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93, doi: 10.1126/science.295.5552.93  \n```\nAndrews, B. H., & Martini, P. 2013, ApJ, 765, 140, doi: 10.1088/0004-637X/765/2/140 Arellano-C'ordova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, ApJL, 940, L23, doi: 10.3847/2041-8213/ac9ab2 Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, arXiv e-prints, arXiv:2304.05378, doi: 10.48550/arXiv.2304.05378 Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2022, arXiv e-prints, arXiv:2211.02495, doi: 10.48550/arXiv.2211.02495 Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5, doi: 10.1086/130766 Bergamini, P., Acebron, A., Grillo, C., et al. 2022, arXiv e-prints, arXiv:2207.09416, doi: 10.48550/arXiv.2207.09416 Bian, F., Kewley, L. J., & Dopita, M. A. 2018, ApJ, 859, 175, doi: 10.3847/1538-4357/aabd74 Bouwens, R. J., Smit, R., Labb'e, I., et al. 2016, ApJ, 831, 176, doi: 10.3847/0004-637X/831/2/176 Boyett, K., Trenti, M., Leethochawalit, N., et al. 2023, arXiv e-prints, arXiv:2303.00306, doi: 10.48550/arXiv.2303.00306 Brinchmann, J. 2022, arXiv e-prints, arXiv:2208.07467. https://arxiv.org/abs/2208.07467\n```  \n```\nBromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23, doi: 10.1086/323947 Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, arXiv e-prints, arXiv:2302.07256, doi: 10.48550/arXiv.2302.07256 Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66 Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692 Cameron, A. J., Katz, H., Rey, M. P., & Saxena, A. 2023a, arXiv e-prints, arXiv:2302.10142, doi: 10.48550/arXiv.2302.10142 Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023b, arXiv e-prints, arXiv:2302.04298, doi: 10.48550/arXiv.2302.04298 Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900 Carnall, A. C., Begley, R., McLeod, D. J., et al. 2023, MNRAS, 518, L45, doi: 10.1093/mnrasl/slac136 Ceverino, D., Glover, S. C. O., & Klessen, R. S. 2017, MNRAS, 470, 2791, doi: 10.1093/mnras/stx1386 Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392 Christensen, L., Laursen, P., Richard, J., et al. 2012, MNRAS, 427, 1973,\n```  \ndoi: 10.1111/j.1365-2966.2012.22007.x", 'A. MZ AND SFR-MZ RELATIONS BASED ONLY ON GALAXY PROPERTIES DERIVED FROM SED FITTING TO NIRCAM PHOTOMETRY': 'In the main text, we present the mass-metallicity (MZ) and the SFR-MZ relations for the full JWST sample with the NIRSpec spectra. Out of the 135 JWST objects with metallicity measurements, 81 objects have NIRCam imaging data, and their stellar masses M /star are derived based on SED fitting to the NIRCam photometry. The remaining 54 objects, all of which are in the CEERS field, are not covered by NIRCam, and their M /star is empirically estimated using UV luminosities from HST photometry (Section 2.3). Additionally, their SFRs are estimated from M UV with some assumptions for comparison with SFR(H β ) (Section 3.4). We note that 23 out of the 54 objects are not detected in HST and have only upper limits on M /star and SFR. Since objects with only an upper limit on M /star are not considered when deriving the average MZ relations and best-fit, there are practically 31 objects (i.e., = 54 -23), accounting for ∼ 28% of the sample, which may introduce bias due to different methods used for M /star and SFR estimation. In this appendix, we present the main figures with only the JWST objects having NIRCam photometry and properties derived from SED fitting, in order to examine how the main results can be affected by excluding the 31 objects with less certain estimations of M /star and SFR.  \nFigure A1 presents the MZ relation using only the JWST objects with M /star based on NIRCam photometry. Their average points are shown with the inverted purple stars. The red stars and the gray long-dashed line represent the average relations based on the full sample of the JWST objects as shown in Figure 10. We confirm that the difference between the two is negligibly small, and the MZ relation based only on the JWST objects with reliable M /star estimates is fully consistent with that based on the full sample, within the uncertainty.  \nLikewise, Figures A2, A3, and A4 are regenerated plots of Figures 11, 13, and 15 from the main text, respectively, using only the JWST objects with NIRCam photometry. These figures confirm that our results remain unchanged when adopting the subsample, indicating that no clear biases are introduced by using the 31 objects with M /star and SFR estimated in a different manner. We have also checked the consistency of M /star between the two methods by comparing the estimates from the 81 CEERS objects with both measurements, and found that ∼ 90% of them have consistent M /star values at the 2 σ level. The remaining 8 objects show overestimated M UV -based M /star values compared to NIRCam-based M /star at > 2 σ level. Although such a small fraction of outliers may exist among the objects lacking NIRCam photometry and introduce additional scatter in the MZ and MZ-SFR relations (e.g., Fig 11(a) versus Fig A2(a) at the high-mass region), the impact is minor as demonstrated in this appendix.  \nFigure A1. Similar to Figure 10, but excluding the individual JWST objects whose stellar masses are empirically derived, so all the plotted data-points have stellar masses that are derived based on SED fitting to NIRCam photometry. The average of the JWST objects plotted on this figure are represented by purple inverted stars. As a comparison, the average relations based on the full sample are shown with the red stars and the gray long-dashed line, as presented in 10, which are positioned almost behind the purple inverted stars. This confirms that our results remain unchanged when adopting the subsample with certain estimations of M /star .  \n<!-- image -->  \nFigure A2. Similar to Figure 11, but excluding the individual JWST objects without JWST/NIRCam photometry. The average relations of the JWST objects in each panel are represented by purple inverted stars.  \n<!-- image -->  \nFigure A3. Similar to Figure 13, but excluding the individual JWST objects without JWST/NIRCam photometry. The average relations of the JWST objects in each panel are represented by purple inverted stars.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure A4. Similar to Figure 15, but excluding the individual JWST objects without JWST/NIRCam photometry. The average relations of the JWST objects in each panel are represented by purple inverted stars. In the z = 8 -10 panel, only a single average point is displayed due to the limited sample size.  \n<!-- image -->', 'B. COMPARING WITH EMPIRICALLY ESTIMATED METALLICITY RELATIONS AT LOW-REDSHIFTS USING STRONG LINE INDICATORS': 'In the main text we refer to Andrews & Martini (2013) and Sanders et al. (2021) to compare our results with the MZ relations at lower redshifts. This is because their metallicities are reliably determined with the direct T e method. Moreover, the z = 0 relation of Andrews & Martini (2013) is explored down to M /star ∼ 10 7 . 4 M /circledot and can be directly compared with most of the JWST objects. In this appendix, we also present comparisons with the other MZ relations at low-redshifts whose metallicities are empirically estimated with the strong line indicators.  \nFigure B5 shows the strong line-based MZ relations at z /lessorsimilar 3 in addition to our results (the average points and the best-fit) and the T e -based MZ relations at z = 0 -3 (Andrews & Martini 2013; Sanders et al. 2021) as in Figure 10. The strong line-based MZ relations include those for z ∼ 0 . 7 galaxies (Savaglio et al. 2005), z ∼ 1 -2 (Papovich et al. 2022), z ∼ 2 -3 (Erb et al. 2006; Steidel et al. 2014), and z ∼ 3 -4 (Troncoso et al. 2014 (see also Maiolino et al. 2008; Mannucci et al. 2009), Onodera et al. 2016). The results by Savaglio et al. (2005) and Erb et al. (2006) are revisited by Maiolino et al. (2008) and their metallicities are re-measured. We also apply corrections to the mass scale to have the same Chabrier (2003) IMF. For the z = 2 -3 galaxies the [N ii ] / H α index is mainly used for the metallicity estimations, while for the other redshift studies combinations of [O ii ], [O iii ] and H β (i.e., R23-index and O32-index in practice) are utilized, in addition to [Ne iii ] when available. Although there are some differences in the MZ relations between different studies and methods for a similar redshift, several factors account for the scatter, such as different SFR activities for the different samples and variations of the metallicity indicators. Despite the slight differences, these overall confirm a redshift evolution of MZ relation, i.e., a decreasing trend of metallicity in the mass range of 10 9 -10 11 M /circledot . Comparing with these strong line-based MZ relations, our main conclusion does not change that the JWST objects at z = 4 -10 are distributed along the simple extrapolation of the MZ relations at z ∼ 2 -3 towards lower mass.  \nFigure C6 shows the redshift evolution of the SFR-MZ relation by only using the JWST objects having µ 0 . 66 > 7 . 5. Their average are represented by the inverted blue stars. As a comparison, the average relations based on the full sample are shown with the red stars, as presented in Figure 13, which are positioned almost behind the blue stars. This confirms that the same results arise from the subsample as discussed in the main text, i.e., (i) the same Andrews & Martini  \n<!-- image -->  \nFigure B5. Similar to Figure 10, but excluding the individual JWST objects, and instead focusing on comparisons with low-redshift MZ relations based on empirically estimated metallicities using strong line indicators, as indicated in the legend. Each MZ relation is depicted in the mass range explored in the original paper.  \n<!-- image -->', 'C. EVOLUTION OF THE SFR-MZ RELATION WITHOUT THE LOW-MASS END GALAXIES': 'Figure 14(a) highlights that the majority of the JWST sample (83 %) has µ 0 . 66 > 7 . 5 and is directly compared with the lower-redshift galaxies in the SFR-MZ relation of Andrews & Martini (2013). Still, we rely on the simple extrapolation of the relation towards lower µ 0 . 66 for the remaining galaxies and discuss the evolution of the SFR-MZ relation together with the direct comparisons of galaxies at µ 0 . 66 > 7 . 5. In this appendix, we present how the low-mass end galaxies at µ 0 . 66 < 7 . 5 can impact the view of the evolution of the SFR-MZ relation.  \n(2013) relation between mass, metallicity, and SFR can explain the properties of galaxies up to z ∼ 8, and (ii) a deficit in metallicity is indicated at z > 8. We therefore conclude that there is no critical impact of the low-mass end galaxies and the extrapolation of the relation at least in the range µ 0 . 66 ∼ 7 -7 . 5 for the discussion of the evolution of the SFR-MZ relation.  \nFigure C6. Similar to the right panel of Figure 13, but excluding the individual JWST objects at µ 0 . 66 < 7 . 5. The average relations of the subsample are denoted with the blue inverted stars. As a comparison, the average relations based on the full sample are shown with the red stars, as presented in Figure 13, which are positioned almost behind the blue inverted stars. This confirms that our results on the redshift evolution of the SFR-MZ relation remain unchanged when considering only the objects at µ 0 . 66 > 7 . 5.  \n<!-- image -->', 'D. SUMMARY OF PROPERTIES FOR JWST OBJECTS IN A TABULAR FORMAT': 'In Table D1, we have compiled the essential properties of the complete sample of JWST objects presented in this paper, based on our improved reduction and calibration of the NIRSpec data from ERO, ERS-GLASS, and ERSCEERS.  \n34  \nNakajima et al.  \nd  \na  \no  \nr  \nB  \n,  \nd  \na  \no  \nr  \nB  \n,  \n1 . 8 6 1 . 3 8 1 . 4 1 1 . 3 4 1 . 3 3 1 . 7 2 2 . 0 6 1 . 7 6 1 . 4 0 1 . 1 9 1 . 7 0 7 . 9 1 1 . 4 8 1 . 4 0 1 . 5 4 2 . 3 9  \n=  \n0  \n7  \n.  \n8  \n=  \n8  \n1  \n.  \n3  \n=  \n8  \n7  \n.  \n1  \n=  \n3  \n2  \n.  \n3  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n=  \n2  \n=  \no  \n2  \n=  \no  \n2  \n=  \no  \n2  \n=  \no  \nd  \na  \no  \nr  \nB  \n,  \n2  \n=  \no  \n2  \n=  \no  \n1  \n0  \n9  \n5  \n4  \n0  \nO  \nR  \nE  \n1  \n4  \n4  \n1  \n5  \n0  \nO  \nR  \nE  \n1  \n5  \n5  \n3  \n6  \n0  \nO  \nR  \nE  \n1  \n0  \n4  \n1  \n8  \n0  \nO  \nR  \nE  \n1  \n2  \n1  \n6  \n0  \n1  \nO  \nR  \nE  \n2  \n3  \n0  \n0  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n7  \n0  \n0  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n4  \n4  \n0  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n7  \n6  \n0  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n4  \n1  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n3  \n2  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n5  \n5  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n6  \n5  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n2  \n6  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n1  \n8  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n6  \n8  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n7  \n9  \n3  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n3  \n0  \n4  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n7  \n0  \n4  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n9  \n3  \n4  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \ne  \nt  \no  \nN  \nD  \nI  \n)  \n¶  \n(  \n0  \n0  \n0  \n0  \n1  \nS  \nS  \nA  \nL  \nG  \n1  \n0  \n0  \n0  \n0  \n1  \nS  \nS  \nA  \nL  \nG  \n3  \n0  \n0  \n0  \n0  \n1  \nS  \nS  \nA  \nL  \nG  \n5  \n0  \n0  \n0  \n0  \n1  \nS  \nS  \nA  \nL  \nG  \n1  \n2  \n0  \n0  \n1  \nS  \nS  \nA  \nL  \nG  \n8  \n0  \n0  \n0  \n5  \n1  \nS  \nS  \nA  \nL  \nG  \n8  \n0  \n0  \n2  \n2  \n2  \n0  \n9  \n2  \n0 0 3 1 0 0 6 0 1 0 0 0 0 0 7 5 0 0  \n5  \nS  \nS  \nA  \nL  \nG  \n1  \nS  \nS  \nA  \nL  \nG  \n1  \nS  \nS  \nA  \nL  \nG  \n5  \nS  \nS  \nA  \nL  \nG  \n8  \nS  \nS  \nA  \nL  \nG  \n1  \nS  \nS  \nA  \nL  \nG  \n6  \n6  \n0  \n0  \n4  \nS  \nS  \nA  \nL  \nG  \n3  \n3  \n1  \n0  \n6  \n1  \nS  \nS  \nA  \nL  \nG  \n9  \n2  \n0  \n0  \n8  \nS  \nS  \nA  \nL  \nG  \nD  \nI  \nMass-Metallicity Star-Formation Relations at  \nz  \n= 4  \n-  \n10  \n35  \na d a d  \no  \nr  \no  \nr  \n8  \n9  \n4  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n1  \n5  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n3  \n5  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n4  \n5  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n4  \n5  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n7  \n5  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n3  \n0  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n1  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n6  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n0  \n7  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n7  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n8  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n8  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n9  \n6  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n0  \n6  \n1  \n7  \n1  \n6  \n4  \n2  \n9  \n9  \n2  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \n3  \n3  \n9  \n0  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n1  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n3  \n2  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n2  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n2  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n2  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n3  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n6  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n6  \n0  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n0  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n1  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n4  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n3  \n4  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n4  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n0  \n6  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n3  \n6  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n3  \n7  \n1  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n0  \n2  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n36  \nNakajima et al.  \nD  \nI  \nd  \na  \no  \nr  \nd  \na  \no  \nr  \nd  \na  \no  \nr  \n7  \n1  \n2  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n3  \n2  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n4  \n2  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n6  \n2  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n8  \n2  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n9  \n2  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n2  \n3  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n3  \n3  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n5  \n3  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n6  \n3  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n7  \n3  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n8  \n3  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n9  \n3  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n0  \n0  \n4  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n1  \n0  \n0  \n1  \n0  \n2  \n3  \n3  \n9  \n4  \n2  \n5  \n0 1 4 0 1 4 0 1 4 0 1 4 0 1 4 0 1 4  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \n5  \n6  \n4  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n6  \n4  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n7  \n4  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n1  \n5  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n3  \n5  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n3  \n5  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n3  \n5  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n4  \n5  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n1  \n6  \n5  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n6  \n5  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n0  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n1  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n0  \n2  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n2  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n3  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n1  \n5  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n5  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n6  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \nD  \nI  \nMass-Metallicity Star-Formation Relations at  \nz  \n= 4  \n-  \n10  \n37  \nd  \na  \no  \nr  \n2  \n7  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n7  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n1  \n9  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n9  \n6  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n0  \n7  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n3  \n7  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n4  \n7  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n5  \n7  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n9  \n5  \n7  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n7  \n6  \n7  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n3  \n8  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n1  \n9  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n3  \n5  \n9  \n1  \n0  \nS  \nR  \nE  \nE  \nC  \n0  \n0  \n0  \n2  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n3  \n9  \n8  \n6  \n1  \n3  \n2  \n0  \n4  \n8  \n6  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \nS  \nR  \nE  \nE  \nC  \n4  \n7  \n1  \n2  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n6  \n3  \n2  \n0  \nS  \nR  \nE  \nE  \nC  \n2  \n8  \n7  \n2  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n8  \n5  \n3  \n0  \nS  \nR  \nE  \nE  \nC  \n4  \n4  \n1  \n4  \n0  \nS  \nR  \nE  \nE  \nC  \n6  \n9  \n1  \n4  \n0  \nS  \nR  \nE  \nE  \nC  \n0  \n1  \n2  \n4  \n0  \nS  \nR  \nE  \nE  \nC  \n8  \n8  \n2  \n8  \n0  \nS  \nR  \nE  \nE  \nC  \n5  \n9  \n9  \n0  \n1  \nS  \nR  \nE  \nE  \nC  \n7  \n1  \n1  \n1  \n1  \nS  \nR  \nE  \nE  \nC  \n3  \n8  \n3  \n1  \n1  \nS  \nR  \nE  \nE  \nC  \n4  \n2  \n6  \n1  \n1  \nS  \nR  \nE  \nE  \nC  \n6  \n7  \n6  \n1  \n1  \nS  \nR  \nE  \nE  \nC  \n8  \n2  \n7  \n1  \n1  \nS  \nR  \nE  \nE  \nC  \n1  \n2  \n2  \n2  \n1  \nS  \nR  \nE  \nE  \nC  \n6  \n9  \n4  \n2  \n1  \nS  \nR  \nE  \nE  \nC  \n7  \n7  \n7  \n4  \n1  \nS  \nR  \nE  \nE  \nC  \n9  \n2  \n3  \n1  \n3  \nS  \nR  \nE  \nE  \nC  \n38  \nNakajima et al.  \n0  \n9  \n5  \n4  \n1  \n9  \n8  \n.  \n4  \n1  \n2  \n7  \n9  \n6  \n7  \n3  \nS  \nR  \nE  \nE  \nC  \n0  \n0  \n5  \n8  \n0  \n9  \n8  \n.  \n4  \n1  \n2  \n2  \n7  \n0  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n6  \n7  \n2  \n1  \n6  \n9  \n.  \n4  \n1  \n2  \n3  \n8  \n0  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n4  \n5  \n0  \n6  \n9  \n8  \n.  \n4  \n1  \n2  \n9  \n3  \n2  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n8  \n9  \n7  \n7  \n2  \n9  \n.  \n4  \n1  \n2  \n2  \n7  \n3  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n4  \n7  \n0  \n8  \n9  \n8  \n.  \n4  \n1  \n2  \n4  \n7  \n3  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n6  \n5  \n0  \n2  \n1  \n8  \n.  \n4  \n1  \n2  \n2  \n3  \n4  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n5  \n1  \n1  \n3  \n4  \n8  \n.  \n4  \n1  \n2  \n5  \n4  \n4  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n4  \n2  \n9  \n3  \n7  \n7  \n.  \n4  \n1  \n2  \n3  \n7  \n5  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n7  \n2  \n6  \n9  \n6  \n7  \n.  \n4  \n1  \n2  \n6  \n7  \n5  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n5  \n6  \n8  \n1  \n7  \n7  \n.  \n4  \n1  \n2  \n6  \n9  \n5  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n5  \n8  \n9  \n4  \n8  \n8  \n.  \n4  \n1  \n2  \n0  \n1  \n7  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n0  \n3  \n6  \n1  \n9  \n8  \n.  \n4  \n1  \n2  \n6  \n1  \n9  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n3  \n8  \n1  \n8  \n0  \n6  \n0  \n9  \n0  \n1  \n0  \n5  \n0  \n2  \n0  \n0  \n1  \n9  \n.  \n4  \n1  \n2  \n4  \n5  \n9  \n0  \n8  \nS  \nR  \nE  \nE  \nC  \n5 4 7 0 4 4 9 8 4 2 6 2 9 8 2 9 1 1  \n0  \n8  \n.  \n4  \n1  \n2  \n8  \n1  \n0  \n1  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n8  \n.  \n4  \n1  \n2  \n2  \n2  \n0  \n1  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n8  \n.  \n4  \n1  \n2  \n6  \n2  \n0  \n1  \n8  \nS  \nR  \nE  \nE  \nC  \n1  \n8  \n.  \n4  \n1  \n2  \n2  \n3  \n0  \n1  \n8  \nS  \nR  \nE  \nE  \nC  \n8  \n7  \n.  \n4  \n1  \n2  \n9  \n4  \n0  \n1  \n8  \nS  \nR  \nE  \nE  \nC  \n9  \n7  \n.  \n4  \n1  \n2  \n3  \n6  \n0  \n1  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n7  \n0  \n5  \n0  \n2  \n8  \n.  \n4  \n1  \n2  \n8  \n6  \n0  \n1  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n6  \n8  \n9  \n9  \n1  \n7  \n.  \n4  \n1  \n2  \n3  \n4  \n0  \n2  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n1  \n6  \n5  \n6  \n6  \n7  \n.  \n4  \n1  \n2  \n2  \n5  \n0  \n2  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n7  \n5  \n5  \n0  \n0  \n9  \n.  \n4  \n1  \n2  \n0  \n0  \n3  \n2  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n3  \n1  \n2  \n3  \n5  \n8  \n.  \n4  \n1  \n2  \n7  \n0  \n5  \n2  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n9  \n0  \n4  \n0  \n0  \n9  \n.  \n4  \n1  \n2  \n8  \n9  \n3  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n0  \n7  \n0  \n8  \n3  \n9  \n.  \n4  \n1  \n2  \n9  \n3  \n4  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n7  \n4  \n8  \n5  \n0  \n9  \n.  \n4  \n1  \n2  \n2  \n0  \n5  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n3  \n9  \n6  \n6  \n5  \n9  \n.  \n4  \n1  \n2  \n2  \n9  \n5  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n7  \n2  \n7  \n5  \n8  \n7  \n.  \n4  \n1  \n2  \n2  \n7  \n7  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n7  \n1  \n4  \n1  \n2  \n8  \n.  \n4  \n1  \n2  \n9  \n7  \n7  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n2  \n2  \n0  \n3  \n0  \n8  \n.  \n4  \n1  \n2  \n6  \n5  \n8  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n8  \n7  \n3  \n8  \n2  \n8  \n.  \n4  \n1  \n2  \n0  \n6  \n8  \n3  \n8  \nS  \nR  \nE  \nE  \nC  \n0  \n4  \n0  \n2  \n1  \n3  \n8  \n.  \n4  \n1  \n2  \n6  \n3  \n8  \n5  \n8  \nS  \nR  \nE  \nE  \nC  \n.  \nA  \n.  \nR  \nD  \nI  \n)  \ng  \ne  \nd  \n(  \n)  \nd  \ne  \nu  \nn  \ni  \nt  \nn  \no  \nc  \n(  \n1  \nD  \ne  \nl  \nb  \na  \nT  \n3  \n2  \nR  \n3  \nR  \n)  \nβ  \nH  \n(  \nW  \nE  \nR  \nF  \nS  \ng  \no  \nl  \n/star  \nM  \ng  \no  \nl  \nV  \nU  \nM  \nt  \nf  \ni  \nh  \ns  \nd  \ne  \nR  \n.  \nl  \nc  \ne  \nD  \n.  \nA  \n.  \nR  \nD  \nI  \n)  \n/star  \n(  \n)  \n/star  \n(  \n)  \n˚A  \n(  \n)  \n1  \n-  \nr  \ny  \n/circledot  \nM  \n(  \n)  \n/circledot  \nM  \n(  \n)  \ng  \na  \nm  \n(  \n)  \ng  \ne  \nd  \n(  \n)  \ng  \ne  \nd  \n(  \nr  \no  \nf  \n)  \nd  \ne  \nt  \na  \nm  \ni  \nt  \ns  \ne  \nt  \no  \nn  \ne  \nr  \na  \ns  \ne  \ni  \nt  \ni  \nc  \ni  \nl  \nl  \na  \nt  \ne  \nm  \ne  \nh  \nt  \ne  \nc  \nn  \ne  \nh  \nd  \nn  \na  \n(  \nn  \nw  \no  \nh  \ns  \nt  \no  \nn  \ne  \nr  \na  \ns  \ne  \nu  \nl  \na  \nv  \ne  \nh  \nT  \n.  \nl  \ne  \nv  \ne  \nl  \nσ  \n3  \ne  \nh  \nt  \nt  \na  \ns  \ne  \nu  \nl  \na  \nv  \nt  \ni  \nm  \ni  \nl  \n-  \nr  \ne  \nw  \no  \nL  \n)  \n/star  \n(  \n-  \ne  \nt  \no  \nN  \na  \ns  \ne  \ns  \ns  \na  \nm  \nr  \na  \nl  \nl  \ne  \nt  \ns  \ne  \nh  \nt  \ns  \nu  \nh  \nt  \nd  \nn  \na  \n,  \ny  \nr  \nt  \ne  \nm  \no  \nt  \no  \nh  \np  \nm  \na  \nC  \nR  \nI  \nN  \ne  \nh  \nt  \ne  \nv  \na  \nh  \n1  \n=  \ng  \na  \nfl  \nh  \nt  \ni  \nw  \ns  \nt  \nc  \ne  \nj  \nb  \nO  \n.  \ng  \na  \nfl  \ny  \nr  \nt  \ne  \nm  \no  \nt  \no  \nh  \nP  \n)  \n†  \n(  \n.  \nd  \ne  \nt  \na  \nu  \nl  \na  \nv  \ne  \ny  \nl  \nr  \ne  \np  \no  \nr  \np  \ne  \nm  \ne  \nT  \nt  \nc  \ne  \nr  \ni  \nd  \ne  \nh  \nt  \n;  \nd  \ne  \nt  \np  \no  \nd  \na  \ns  \ni  \nd  \no  \nh  \nt  \ne  \nm  \nh  \nc  \ni  \nh  \nw  \ng  \nn  \ni  \ny  \nf  \ni  \nc  \ne  \np  \ns  \n,  \ng  \na  \nfl  \ny  \ni  \nt  \ni  \nc  \nl  \nl  \na  \nt  \ne  \nM  \n)  \n‡  \n(  \n.  \nT  \nS  \nH  \nf  \no  \ns  \ne  \ni  \nt  \ni  \ns  \no  \nn  \ni  \nm  \nu  \nl  \nV  \nU  \ne  \nh  \nt  \nh  \nt  \ni  \nw  \nd  \ne  \nt  \na  \nm  \ni  \nt  \ns  \ne  \nMass-Metallicity Star-Formation Relations at  \nz  \n= 4  \n-  \n10  \n39  \n)  \nO 3 2 l o g ( O / H ) F l a g s ( /star ) + 1 2 P ( † ) M ( ‡ t h o s e w h o s e m e a s u r e d f t h e d u s t -r e d d e n i n g i s r e r e l i a b l y d e t e r m i n e d m a s s e s a n d S F R s a r e t h o d o r t h e e m p i r i c a l t s , t h e v a l u e s o f M U V , p o i n t i n g s / g r a t i n g s , t h e e c t s w i t h s p e c t r o s c o p i c r d e r i v i n g t h e M Z a n d t i m a t i o n b a s e d o n t h e  \ni  \nn  \ne  \nv  \ni  \ng  \ne  \nr  \na  \n2  \n3  \nO  \nd  \nn  \na  \n3  \n2  \nR  \n.  \nl  \ne  \nv  \ne  \nl  \nσ  \n5  \n.  \n2  \ne  \nh  \nt  \nt  \na  \n8  \n9  \n.  \n2  \nf  \no  \ne  \nu  \nl  \na  \nv  \nl  \na  \nc  \ni  \nt  \ne  \nr  \no  \ne  \nh  \nt  \ne  \nh  \nt  \nh  \nt  \ni  \nw  \nt  \nn  \ne  \nt  \ns  \ni  \ns  \nn  \no  \nc  \nt  \no  \nn  \ns  \ni  \no  \ni  \nt  \na  \nr  \ne  \nn  \ni  \nl  \n9  \n5  \n9  \n4  \n/  \n7  \n0  \n0  \n5  \nλ  \nλ  \n]  \ni  \ni  \ni  \nO  \n[  \nr  \ni  \ne  \nh  \nt  \n,  \n0  \n=  \ng  \na  \nfl  \nh  \nt  \ni  \nw  \ns  \nt  \nc  \ne  \nj  \nb  \no  \ne  \nh  \nt  \nr  \no  \nF  \n.  \ns  \ne  \ni  \nt  \ni  \ns  \no  \nn  \ni  \nm  \nu  \nl  \nβ  \nH  \nd  \ne  \nt  \nc  \ne  \nr  \nr  \no  \nc  \ns  \ns  \no  \nl  \n-  \nt  \ni  \nl  \ns  \ne  \nh  \nt  \nh  \nt  \ni  \nw  \nd  \ne  \nv  \ni  \nr  \ne  \nd  \ne  \nr  \na  \ns  \nR  \nF  \nS  \nd  \nn  \na  \nt  \nfi  \nD  \nE  \nS  \ne  \nh  \nt  \nh  \nt  \ni  \nw  \nj  \nb  \no  \ne  \nh  \nT  \n.  \nn  \ne  \nv  \ni  \ng  \ne  \nr  \na  \ny  \nt  \ni  \nc  \ni  \nl  \nl  \na  \nt  \ne  \nm  \nd  \nn  \na  \n,  \ns  \no  \ni  \nt  \na  \nr  \ne  \nn  \ni  \nl  \n,  \n)  \nβ  \nH  \n(  \nW  \nE  \nf  \no  \ns  \ne  \nu  \nl  \na  \nv  \nd  \ne  \ng  \na  \nr  \ne  \nv  \na  \nr  \ni  \ne  \nh  \nt  \nd  \nn  \na  \n,  \nd  \ne  \nt  \no  \nn  \ns  \ni  \ns  \ng  \nn  \ni  \nt  \na  \nr  \ng  \n/  \ns  \ng  \nn  \ni  \nt  \nn  \ni  \no  \np  \nf  \no  \nr  \ne  \nb  \nm  \nu  \nn  \no  \nf  \nd  \ne  \ns  \nu  \nt  \no  \nn  \ne  \nr  \na  \ns  \nt  \nc  \ne  \nj  \nb  \no  \nd  \na  \no  \nr  \nB  \ne  \ns  \ne  \nh  \nT  \n.  \nd  \na  \no  \nr  \nB  \nh  \nt  \ni  \nw  \nd  \ne  \nk  \nr  \na  \nm  \ne  \nr  \na  \n)  \n3  \n2  \n0  \n2  \n.  \nl  \na  \nt  \ne  \ne  \nn  \na  \nk  \ni  \nr  \na  \nH  \n(  \ny  \nb  \nd  \ne  \nt  \na  \nc  \ni  \nd  \nn  \ni  \ns  \na  \ns  \nN  \nG  \nA  \nf  \no  \ns  \ne  \nr  \nu  \nt  \na  \nn  \ng  \ni  \ns  \nc  \ne  \nj  \nb  \no  \nd  \ne  \ns  \nn  \ne  \nl  \ne  \nh  \nt  \nr  \no  \nF  \n.  \nd  \ne  \nt  \no  \nn  \ns  \ni  \nµ  \nr  \no  \nt  \nc  \na  \nf  \nn  \no  \ni  \nt  \na  \nc  \nfi  \ni  \nn  \ng  \na  \nm  \ne  \nh  \nt  \n,  \ns  \nt  \nc  \ne  \nj  \nb  \no  \nd  \ne  \ns  \nn  \ne  \nl  \ny  \nl  \nl  \na  \nn  \no  \ni  \nt  \na  \nt  \ni  \nv  \na  \nr  \ng  \ne  \nh  \nt  \nr  \no  \nF  \n)  \n¶  \n(  \n.  \n3  \nR  \nr  \no  \n3  \n2  \nR  \nf  \no  \nr  \no  \nt  \na  \nc  \ni  \nd  \nn  \ni  \nl  \na  \nr  \ne  \nv  \ne  \ns  \nn  \ni  \nd  \ne  \nv  \nr  \ne  \ns  \nb  \no  \ny  \nl  \np  \ni  \nt  \nl  \nu  \nm  \ne  \nr  \na  \nt  \na  \nh  \nt  \ns  \nt  \nc  \ne  \nj  \nb  \no  \nS  \nR  \nE  \nE  \nC  \ne  \nh  \nt  \nr  \no  \nF  \n.  \nn  \no  \ni  \nt  \na  \nc  \nfi  \ni  \nn  \ng  \na  \nm  \ne  \nh  \nt  \nr  \no  \nf  \nd  \ne  \nt  \nc  \ne  \nr  \nr  \no  \nc  \ny  \nd  \na  \ne  \nr  \nl  \na  \ne  \nr  \na  \nR  \nF  \nS  \nd  \nn  \na  \ns  \ns  \na  \nm  \ns  \ne  \ns  \ns  \na  \nm  \nr  \na  \nl  \nl  \ne  \nt  \ns  \ne  \nh  \nt  \nn  \ni  \ny  \nt  \nn  \ni  \na  \nt  \ne  \nc  \nn  \nu  \n,  \n.  \ne  \n.  \ni  \n(  \nl  \ne  \nv  \ne  \nl  \nm  \nu  \nu  \nn  \ni  \nt  \nn  \no  \nc  \ne  \nh  \nt  \no  \nt  \ns  \nN  \nG  \nA  \nf  \no  \nn  \no  \ni  \nt  \nu  \nb  \ni  \nr  \nt  \nn  \no  \nc  \nn  \nw  \no  \nn  \nk  \nn  \nu  \ne  \nh  \nt  \no  \nt  \ne  \nu  \nd  \ns  \nn  \no  \ni  \nt  \na  \nl  \ne  \nr  \nZ  \nM  \n-  \nR  \nF  \nS  \n.  \n)  \ng  \nn  \ni  \nt  \nt  \nfi  \nD  \nE  \nS'}