JWST Census for the Mass–Metallicity Star Formation Relations at z = 4–10 with Self-consistent Flux Calibration and Proper Metallicity Calibrators

2.1.1. NIRSpec Observations and Data Reduction

The 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). Multislit spectroscopy was taken using the microshutter assembly (MSA), with the medium-resolution (R ∼ 1000) gratings/filters of G235M/F170LP and G395M/F290LP sampling the wavelength ranges 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 s in G395M and F235M each, with a three-shutter slitlet nod pattern, i.e., slightly nodded positions, nod1, nod2, and nod3, in different shutters within the slitlets. The total integration time was 4.86 hr in G235M and G395M each.

By visually checking all of the 1D and 2D composite spectra that were publicly available, we identified five 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 2023; Mahler et al. 2023), they were not included in the following analysis because of their lines' poorer confidence levels.

For the five sources, data reduction was reperformed for a more careful analysis and flux calibration. Starting with the Level 1 products provided by the observatory, we used the Spec2 and Spec3 pipelines per nod and per observing block using the Python library for JWST observations by STScI (v1.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, ⁸ when using the Spec2 pipeline, a single science image for a particular nodding position was background subtracted by using the other nod 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 noise, we carefully checked the background-subtracted science images and removed one of the other nod images, if necessary, for the background images and reran 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 [−10 pixel:+10 pixel] along the spectral direction. ⁹ We applied 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-light radii smaller than the point-spread function (PSF; e.g., Harikane et al. 2023a; Tacchella et al. 2023; see also Ono et al. 2023 for the size evolution at high redshift). This assumption affected how the path-loss corrections, which accounted for various types of signal loss in the 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 (three nods for two observing blocks) for each of the ERO sources.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Checking 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órdova 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 when making the composite. In short, the integration time we used for the composite spectra of 04590 and 05144 was 4.05 hr, and 4.86 hr for 06355, 08140, and 10612.

The composite was made by median stacking the available five to six 2D spectra. The individual 2D spectra were residual background subtracted, and shifted to have common spatial and spectral ranges for the composite. All the shifts were integer movements. We adopted a 2D median-stacking procedure to alleviate further 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 flux losses were >3σ if extracted in individual frames and compared with the other nod 1D spectra. 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, a 1D spectrum was produced for each source and grating via the summation of 3 pixels (= 03 ≃ 2 × FWHM) along the spatial direction centered on the spatial peak position. We adopted a narrow extraction aperture as compared to the previous studies to minimize the effects coming from noisy regions close to the edge, particularly for 04590, 08140, and 10612 that were observed with a short MSA slitlet.

The 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 preflight throughputs. As for 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 quite accurate, with the uncertainties at approximately 5% or less over the entire wavelength range. ¹⁰ We therefore used the flux solutions provided by the Spec2 pipeline 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 quantities such as the equivalent widths (EWs) of the optical emission lines and ξ_ion. We used the NIRCam broadband data 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.

To evaluate the noise level, we used the 2D noise frames of read-out noise (VAR\inferior RNOISE) and Poisson noise (VAR\inferior POISSON) outputted in the Spec2 procedure for each object and nod. In addition, we added the standard deviation of the residual background in 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 [O iii] λ5007 (i.e., the strongest emission line) was adopted for the Gaussians of the other lines. We perform all wavelength measurements on the vacuum wavelength scale, although we follow a long-standing 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 measurements and 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 the EWs of the optical emission lines by combining the fluxes with the continuum level provided from a best-fit spectral energy distribution (SED) to the appropriate rest-frame optical broadband 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.

Table 1. Summary of JWST Objects Studied with the Direct T_e Method

ID	[O ii] a	[Ne iii] a	Hδa	Hγa	[O iii] a	Hβa	[O iii] a	[O iii] a	Redshift	EW(Hβ)	Te (O iii)	12+log (O/H)
	λ3727	λ3869			λ4363		λ4959	λ5007		(Å)	(104 K)
ERO_04590	<5.7	<24.8	21.3 ± 2.3	42.6 ± 2.8	13.8 ± 2.7	100.0 ± 6.4	110.5 ± 6.0	337.8 ± 8.3	8.496	217.8 ± 150.5	2.08 ± 0.26	${7.26}_{-0.13}^{+0.15}$
ERO_05144	49.6 ± 6.2	63.7 ± 8.9	28.5 ± 3.2	46.8 ± 3.1	14.9 ± 2.9	100.0 ± 3.0	238.3 ± 3.9	719.0 ± 5.5	6.378	150.9 ± 51.2	1.64 ± 0.16	${7.82}_{-0.09}^{+0.12}$
ERO_06355 b	99.6 ± 3.7	48.2 ± 2.1	24.6 ± 1.7	47.4 ± 2.1	8.7 ± 1.8	100.0 ± 2.5	267.8 ± 7.9	817.1 ± 5.8	7.665	150.3 ± 4.4	1.13 ± 0.08	${8.35}_{-0.08}^{+0.11}$
ERO_10612	43.8 ± 5.1	54.7 ± 4.0	26.7 ± 3.0	45.7 ± 3.1	22.1 ± 2.8	100.0 ± 3.5	238.6 ± 4.9	712.3 ± 6.7	7.660	210.3 ± 16.4	2.04 ± 0.16	${7.59}_{-0.08}^{+0.08}$
GLASS_100003	37.0 ± 11.9	29.2 ± 8.2	<7.2	51.3 ± 8.0	25.7 ± 7.7	100.0 ± 10.1	259.2 ± 12.1	777.8 ± 16.5	7.877	138.7 ± 19.2	1.98 ± 0.37	${7.64}_{-0.19}^{+0.23}$
GLASS_10021	68.3 ± 5.1	48.4 ± 5.1	24.3 ± 3.5	49.4 ± 4.3	19.8 ± 4.1	100.0 ± 4.6	269.6 ± 5.2	828.4 ± 6.8	7.286	67.9 ± 27.0	1.66 ± 0.17	${7.87}_{-0.10}^{+0.12}$
GLASS_150029 c	42.0 ± 3.3	⋯	22.1 ± 3.4	47.1 ± 4.0	15.3 ± 2.3	100.0 ± 4.5	222.1 ± 5.8	631.6 ± 7.9	4.584	137.9 ± 10.3	1.76 ± 0.14	${7.70}_{-0.08}^{+0.09}$
GLASS_160133 c	47.5 ± 2.1	45.0 ± 2.0	23.8 ± 1.2	46.5 ± 1.5	13.4 ± 1.4	100.0 ± 1.6	255.9 ± 2.0	772.5 ± 3.1	4.015	227.2 ± 89.3	1.48 ± 0.07	${7.95}_{-0.05}^{+0.06}$
CEERS_01027	45.4 ± 9.2	59.5 ± 6.1	25.5 ± 5.0	38.8 ± 4.4	15.6 ± 5.0	100.0 ± 6.4	220.9 ± 7.5	738.1 ± 10.9	7.820	174.7 ± 14.4	1.56 ± 0.24	${7.83}_{-0.15}^{+0.20}$
CEERS_01536	59.5 ± 11.6	43.1 ± 9.0	<12.9	63.4 ± 10.0	30.7 ± 8.2	100.0 ± 9.6	296.1 ± 15.6	805.4 ± 16.5	5.034	⋯	2.25 ± 0.43	${7.61}_{-0.15}^{+0.18}$

Notes.

^a Observed flux ratios relative to Hβ (×100). Upper-limit values at the 1σ level. ^b The presence of a high-ionization line of [Ne iv] is seen (Brinchmann 2023). ^c The presence of broad Hα is indicated (Harikane et al. 2023b).

A machine-readable version of the table is available.

Download table as:  DataTypeset image

To 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 = 10,000 and 25,000 K). The four ERO objects all follow the sequence and present physically reasonable Balmer decrements within the uncertainties. The 3σ upper limits on Hγ/Hβ and Hδ/Hβ, along with the observed Hα/Hβ ratio, are consistent with the sequence of the expected Balmer decrements for 08140. In contrast, previous studies claim some tensions of the observed Balmer decrements based on the early release products, which show inconsistencies 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 two-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 (Arellano-Córdova et al. 2022; Curti et al. 2023a), mostly due to the addition of the uncertainty of the subtracted residual background into the noise levels.

Zoom In Zoom Out Reset image size

Prior 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 Case B recombination with T_e = 17,500 K. Adding the lower signal-to-noise ratio emission of Hδ does not change the estimates, as indicated in Figure 2. For 08140, we use the Hα/Hβ ratio instead. We adopt E(B−V) = 0.22 for 04590, 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

The ERS GLASS NIRSpec MSA observations were carried out on UT 2022 November 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 ranges of 1.0–1.9, 1.7–3.2, and 2.8–5.3 μm, respectively. Twelve exposures of 1473 s were conducted in one pointing of NIRSpec MSA for each spectral configuration with the three-shutter slitlet nod pattern, totaling an on-source integration of 4.9 hr in each of G140H, G235H, and G395H.

Through the initial screening of all of the public 1D and 2D composite spectra, we identified 15 objects at z > 4 with [O iii]+Hβ. ¹¹ The sample contains the members of the protocluster at z = 7.89 as presented by Morishita et al. (2023). We rereduced the spectra of the 15 objects with our custom procedure as described in Section 2.1.1. We carefully checked all of the 12 nod spectra for each object to determine the extraction aperture and any nods to be removed from the composition. For 10005, the light from another member of the z = 7.89 protocluster partly came into the same slit at a spatially offset position. We changed the background frames during the spec2 procedure so that the signal from 10005 was not oversubtracted by the neighboring object in the different nods.

The emission-line flux measurements were performed in the same manner as explained in Section 2.1.2. Among the 15 objects, we identified five with [O iii] λ4363. One of the [O iii] λ4363-detected source, 40066, lacked [O ii] λ3727 emission which fell in the detector gap. The other four 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γ.

The NIRCam multiband 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 is not determined.

All 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 Multi Unit Spectroscopic Explorer (MUSE) follow-up spectroscopy (Bergamini et al. 2023). 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 042. 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.

The ERS CEERS NIRSpec MSA observations were carried out on UT 2022 December 20–22 and 24 in the blank field of EGS (Finkelstein et al. 2023; Fujimoto et al. 2023; Proposal ID: 1345). They used both prism (R ∼ 100) and medium-resolution (R ∼ 1000) gratings to cover from 1.0 to 5 μm in six 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 February 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 spectroscopic sample of CEERS to the community. Three exposures of 1036 s were conducted for each spectral configuration with the three-shutter slitlet nod pattern, totaling an on-source integration of 0.9 hr in each prism, G140M, G235M, and G395M observation, except for P9 and P11 where another set of three exposures of 1036 s were repeated in the prism observations.

At 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 the 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 rereduction procedure as described in Section 2.1.1.

The 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 CEERS 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.

The CEERS NIRCam observations consists of 10 pointings, four of which were taken in 2022 June and the rest in 2022 December (Bagley et al. 2023). The early observations of four pointings are presented in Harikane et al. (2023a). 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. 2023, in preparation). One caveat is that not all of the NIRSpec CEERS objects are covered by NIRCam. We find 91 of the 163 CEERS objects have NIRCam coverage. For the 91 objects, we perform 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 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 stellar 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 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 the 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 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 EWs of the optical emission lines such as EW(Hβ), and the SFR based on the total Hβ luminosity.

In 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.

Zoom In Zoom Out Reset image size

To examine the nebular properties of JWST objects at z = 4–9 and compare with those of lower-redshift galaxies, we plot them in a 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 kind of 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 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 redshift. Figure 4 confirms the trend, and the z = 4–9 JWST objects show a further high 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 $\mathrm{log}U\gt -2.21$ according to the prescription of Kewley & Dopita (2002; see also Fujimoto et al. 2022). ¹² The lower limit is a factor of ∼10 higher value than typically seen in SDSS galaxies.

Zoom In Zoom Out Reset image size

At z = 2–3, low-mass LAEs are suggested to have a remarkably high O32 line ratio, which would be achieved by a hard ionizing radiation field and/or 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 low-mass LAEs and z ∼ 0 Lyman-continuum (LyC) leaking objects, suggesting that high-redshift sources have a nebular and ionization condition of gas that is similar to those typically found in lower-redshift low-mass galaxies having strong Lyα and/or LyC escape with a hard ionizing spectrum. Similar findings are also discussed by Mascia et al. (2023) and Sanders et al. (2023), and notably by Tang et al. (2023), where the strong LAE galaxies at z > 7 are directly examined with JWST and suggested to present the largest O32 line ratios.

One 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.

Before 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 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 the 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 AGN dominated.

Zoom In Zoom Out Reset image size

Figure 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 ≲ 0.5 Z_⊙ 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.

As a complementary approach to our method using the [N ii] BPT diagram, a companion study by Harikane et al. (2023b) conducts a search for faint AGNs in our full JWST spectroscopic sample by examining the broad component (FWHM >1000 km s⁻¹) 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; Übler et al. 2023). To avoid potential biases in our MZ relation, we remove these 10 objects from our sample when we examine the MZ relation. These objects are flagged in Table D1.

3.2.1. Direct T_e Method

Ten objects (four from ERO, four from GLASS, and two from CEERS) present a significant detection of [O iii] λ4363 allowing us to determine the oxygen abundance reliably with the direct T_e method at z = 4.0–8.5 and use it as a proxy for the 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 the O²⁺ zone (T_e ([O iii])) using the reddening-corrected [O iii] λ λ4363/5007 ratio and an assumed electron number density of 100 cm⁻³ with the task getTemDen from PyNeb (Luridiana et al. 2015). Different densities such as 10–1000 cm⁻³ 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 the O⁺ zone, T_e ([O ii]), is extrapolated from T_e ([O iii]) employing the prescription of Izotov et al. (2006; see Brinchmann 2023 for a different assumption). We derive the abundance of O²⁺/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³⁺/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²⁺ 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.05 dex, which is small enough with respect to the measurement error (≃0.15 dex). The measured oxygen abundances as well as T_e ([O iii]) are summarized in Table 1.

For the ERO objects, several earlier studies exist that report metallicities with electron temperatures (Schaerer et al. 2022; Curti et al. 2023a; Trump et al. 2023; Rhoads et al. 2023; Arellano-Córdova et al. 2022; Brinchmann 2023). It is practically useful to compare our measurements with theirs. 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⁴ K, which is not usually observed in the local Universe and may require some additional explanations (e.g., Katz et al. 2023b; Rhoads et al. 2023). Our remeasurement confirms a high, but not extremely high, electron temperature for ERO_04590, T_e = (2.08 ± 0.26) × 10⁴ K, which is explainable by heating by young massive stars without any special mechanisms. It is thus important to reduce and combine the spectra carefully, as well as to extract the 1D spectra to discuss the nebular properties reliably with the measurements of faint emission lines. For the other objects, their electron temperatures are modest, T_e = (1.1–2.3) × 10⁴ 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 MZ relationship.

Zoom In Zoom Out Reset image size

The four ERO objects as well as the six 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 proposed by 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).

Zoom In Zoom Out Reset image size

The 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 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 object whose EW(Hβ) is smaller than 100 Å, 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.

In 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 good agreement 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 metal-poor regime.

Figure 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 nondetection 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 the 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+\mathrm{log}({\rm{O}}/{\rm{H}})$ ∼ 7.8.

Zoom In Zoom Out Reset image size

One caveat should be noted for ERO_06355, which has an oxygen abundance of $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ = 8.3 and shows an unusually strong [O iii] line. Curti et al. (2023a) have used the [O iii]/Hβ versus [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 (2023) tentatively suggest the presence of the high-ionization lines [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β) values that agree reasonably well with the relationships for star-forming galaxies, implying a minor contribution from the 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 discuss further the presence of hard radiation components in these systems.

3.2.2. Empirical Method

3.2.3. Other High-redshift Galaxies from the Literature

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 MZ relation at high redshift using the JWST observations, as well as data compiled from relevant literature.

Before presenting the full results we first focus on the four ERO objects and compare in Figure 9 our MZ relation with those presented in the earlier studies (Arellano-Córdova et al. 2022; Schaerer et al. 2022; Curti et al. 2023a; Trump et al. 2023; Rhoads et al. 2023; Brinchmann 2023). 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. (2023) and hence both values are shown in the plot. We confirm that our MZ relation for the four ERO objects overall shows good agreement with the earlier results. One note is that we confirm the metallicity of ERO_04590 is not extremely low, $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ = 7.26 (+0.15/−0.13) for its stellar mass.

Zoom In Zoom Out Reset image size

Figure 10 illustrates the relation between stellar mass and metallicity (MZ) with the full sample at z = 4–10. The MZ relation determined at z = 0, ≃2.3, and ≃3.3 with the direct T_e method is also plotted (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_⋆ ∼ 10⁵ M_⊙). The red symbols denote the galaxies at z = 4–10, including our 135 ERO, GLASS, and CEERS galaxies, the three 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. 2023b) have been excluded here to discuss conservatively the results free from any AGN biases (Section 3.1). The objects analyzed using the direct T_e method are marked with red open circles, suggesting a positive correlation between mass and metallicity in the redshift range z = 5–8.5. Furthermore, we divide our full sample of ERO, GLASS, and CEERS galaxies into three subsamples according to stellar mass: M_⋆ = 10⁷–10⁸, 10⁸–10⁹, and 10⁹–10¹⁰ M_⊙, and obtain the average MZ relation as shown 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_⋆, are not used in deriving the average relation. 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

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O}}/{\rm{H}})={Z}_{10}+\gamma \mathrm{log}({M}_{\star }/{10}^{10}\,{M}_{\odot }),\end{eqnarray} \tag{ 1 }$$

with the best-fit parameters of Z₁₀ = 8.24 ± 0.05 and γ = 0.25 ± 0.03 in the mass range of M_⋆ ∼10^7.5–10^9.5 M_⊙, 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. 2023), and as widely known at low redshift.

Zoom In Zoom Out Reset image size

Table 2. Average Values of the Stellar Mass, Metallicity, and SFR for our JWST Sample

Sample	$\mathrm{log}{M}_{\star }$	$12+\mathrm{log}({\rm{O}}/{\rm{H}})$	log SFR
	(M⊙)		(M⊙ yr−1)
For the MZ relation (including objects without an SFR measurement)
z = 4−10	7.60 ± 0.21	7.61 ± 0.30	⋯
	8.56 ± 0.27	7.95 ± 0.20	⋯
	9.27 ± 0.21	8.03 ± 0.26	⋯
z = 4−6	8.30 ± 0.47	7.87 ± 0.28	⋯
	9.15 ± 0.26	8.04 ± 0.23	⋯
z = 6−8	7.85 ± 0.29	7.77 ± 0.27	⋯
	8.76 ± 0.38	7.90 ± 0.22	⋯
z = 8−10	8.29 ± 0.49	7.63 ± 0.26	⋯
	9.29 ± 0.59	7.96 ± 0.23	⋯
For the SFR–MZ relation
z = 4−10	7.56 ± 0.21	7.50 ± 0.28	0.55 ± 0.38
	8.55 ± 0.27	7.91 ± 0.21	1.02 ± 0.31
	9.27 ± 0.22	8.05 ± 0.26	1.45 ± 0.44
z = 4−6	8.32 ± 0.52	7.81 ± 0.31	0.84 ± 0.40
	9.21 ± 0.25	8.06 ± 0.25	1.45 ± 0.42
z = 6−8	7.94 ± 0.29	7.72 ± 0.26	0.77 ± 0.29
	8.79 ± 0.34	7.90 ± 0.21	1.19 ± 0.32
z = 8−10	8.13 ± 0.53	7.54 ± 0.28	1.03 ± 0.37
	9.53 ± 0.59	7.84 ± 0.20	1.37 ± 0.60

Note. The first half of this table displays the average masses and metallicities of the samples used for the the MZ relation, while the second half shows the samples used for the SFR–MZ relation. The second group of samples excludes objects whose SFRs are not measured or constrained (Section 3.4), resulting in a smaller sample size compared to the first half.

Download table as:  ASCIITypeset image

Compared to the z = 0 MZ relation, these high-redshift galaxies clearly present a metallicity lower than typical galaxies at z = 0 for a given stellar mass. The decrease is typically ∼0.5 dex around M_⋆ ∼10⁹ M_⊙, but it becomes smaller at the low-mass end (∼0.3 dex). Interestingly, a similar offset of ∼0.3 dex is observed between z = 0 and z ∼ 2–3, suggesting that the evolution of the MZ relation is small from z ∼ 2–3 to z = 4–10. Although there may be a decrease of ∼0.2 dex in the typical metallicity at the high-mass end of M_⋆ ∼10^9.5 M_⊙ from z ∼ 2–3 to z = 4–10, no strong evolution is found beyond the errors. The same conclusion can be drawn from comparisons with the MZ relation at z ≲ 4, where the metallicities are empirically estimated with the strong line indicators, as presented in Appendix B.

Figure 10 also includes a comparison with the latest results from the JADES observations. Curti et al. (2023b) very recently reported 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, the latter of which is provided in this paper. We find that their MZ relation, despite covering slightly different redshift ranges, agrees well with ours (Equation (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_⋆ ≳ 10⁸ M_⊙), since this regime is dominated by the CEERS objects presented in this paper. The JADES objects confirm that the MZ relation we obtained in Equation (1) continues down to M_⋆ ≃ 10⁷ M_⊙ on average.

In Appendix A, we investigate our MZ relation using only galaxies for which the stellar mass is well determined 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 CEERS objects with only an upper limit on M_⋆, all of which are M_UV based, in deriving the average relation, the fraction of objects without NIRCam photometry is small (∼28%). Indeed, our conclusions remain unchanged when considering only objects with reliable measures of M_⋆, 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 versus 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.

To examine further 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 subsample 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 the different redshift bins are consistent with those based on the full sample at z = 4–10 (Equation (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 toward 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.

Zoom In Zoom Out Reset image size

The 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 seminumerical 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. (2023) in yellow, and FLARES by Wilkins et al. (2023) in magenta (see also Lovell et al. 2021). We plot the theoretical predictions at z = 5, 7, and 9 for the redshift bins 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 (${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})=-0.025$ dex from z = 6 to 7, and −0.05 dex 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 (${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})=-0.08$ dex from z = 4 to 5, −0.075 dex from z = 6 to 7, and −0.1 dex 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 regimes, 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 (<10 Myr) 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 models are SFR weighted and mass weighted by young (<100 Myr) 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.

Comparing 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_⋆ = 10⁸–10^9.5 M_⊙. Notably, the slope of the MZ relation found in the IllustrisTNG results likely coincides with the observations. On the other hand, the simulations of FIRE and Astraeus are generally suggested to underpredict 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 feedback, i.e., their galaxies eject too much metal and/or accreting gas is too efficient in lowering the metal content in the galaxies (Ma et al. 2016; Ucci et al. 2023). Finally, we note that the average metallicity measured for low-mass galaxies below M_⋆ < 10⁸ M_⊙ can be slightly higher than the existing predictions by ∼0.2 dex. 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 theoretical predictions such as IllustrisTNG toward lower mass to conclude the consistency and to understand better the early chemical enrichment in the low-mass systems at high redshift.

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ópez 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): ${\mu }_{\alpha }=\mathrm{log}({M}_{\star })-\alpha \mathrm{log}(\mathrm{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): ¹³

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O}}/{\rm{H}})=0.43\times {\mu }_{0.66}+4.58.\end{eqnarray} \tag{ 2 }$$

This is advantageous to be directly compared with a majority of the JWST objects down to M_⋆ ∼ 10^7.4 M_⊙. Interestingly, the z ≃ 2.3 and ≃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.

To estimate the SFRs for the JWST objects, we adopt the total, reddening-corrected Hβ luminosity for each object, as it is the best indicator for ongoing (≲10 Myr) star formation activity. The SFR relation of Kennicutt (1998) is adopted for the Balmer lines with a correction of the 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 corrections are successfully applied and whose spectra are slit-loss corrected. The latter constraint depends on whether the object has NIRCam coverage (only in the CEERS field; Section 2.3). Among the CEERS objects lacking a slit-loss correction, we rescue 42 objects whose UV stellar emission is constrained with HST. For these objects, we translate M_UV into SFR(Hβ) assuming a typical ionizing photon production efficiency as indicated at high redshift ($\mathrm{log}{\xi }_{\mathrm{ion}}=25.6\pm 0.2;$ e.g., Endsley et al. 2023; Matthee et al. 2023) and a 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. 2023, in preparation) using the JWST sample presented in this paper. The objects lacking a proper dust reddening correction or having just upper limits on both M_⋆ 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.

In Figure 12, we show the individual and average distributions of the JWST objects on the stellar mass–SFR plane using the SFR values as derived above (i.e., mainly from Hβ). The high-redshift objects are distributed along the sequence of specific SFR (sSFR) = 10⁻⁸–10⁻⁷ yr⁻¹. This is overall consistent with the star formation main sequence of galaxies at z = 4–7, where an sSFR of 10^−8.5–10^−7.5 yr⁻¹ 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 galaxies are above sSFR ≳10^−7.5 yr⁻¹, 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.

Zoom In Zoom Out Reset image size

Figure 13(a) shows 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 (Equation (2)) and the data points of z ≃ 2.3 and 3.3 stacked galaxies (Sanders et al. 2021). Adopting the best-fit parameter α = 0.66 found in 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.

Zoom In Zoom Out Reset image size

We 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 star-forming 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_⋆ ∼ 10⁵–10⁷ M_⊙ and $12+\mathrm{log}({\rm{O}}/{\rm{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 (2013) galaxies (Nishigaki et al. 2023). This evidence supports the simple extrapolation of the relation toward 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 occupy 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⁻¹), a weaker SFR dependence is indicated (α = 0.22). On the other hand, the high-sSFR subsample with sSFR > 10^−9.5 yr⁻¹ 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⁻⁷–10⁻⁸ yr⁻¹; see Figure 12) and also for the z = 2–3 galaxies (sSFR ∼ 10^−8.5 yr⁻¹, we conclude the strong SFR dependence of α = 0.66 as found by Andrews & Martini (2013) is preferred and adopted in this paper to be compared with the high-redshift galaxies. We discuss the evolution of the Curti et al. (2020) SFR–MZ relation later in Section 4.

Zoom In Zoom Out Reset image size

Figure 13(b) clarifies the redshift evolution of 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: ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ = $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ _obs − $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ _μ0.66 (Equation (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.

Two 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 ≃ 0.3 dex. Another is that a significant decrease of metallicity is found beyond the errors 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 11(c)). To examine further 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 the 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 subsample 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 subsample 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.

Zoom In Zoom Out Reset image size

Finally, 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 masses are well determined with SED fits to the NIRCam photometry and whose SFRs are obtained with the slit loss–corrected Hβ. Our conclusions remain unchanged even with the smaller sample that has good measurements of M_⋆ 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 relation 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 high-redshift 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 of galaxies in our sample does not introduce biases. Furthermore, we note that a consistent view of the SFR–MZ relation 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).